Math Is Fun Forum

  Discussion about math, puzzles, games and fun.   Useful symbols: ÷ × ½ √ ∞ ≠ ≤ ≥ ≈ ⇒ ± ∈ Δ θ ∴ ∑ ∫ • π ƒ -¹ ² ³ °

You are not logged in.

#1 2012-12-25 05:49:02

Giewont
Member
Registered: 2012-12-25
Posts: 3

Gram-Schmidt matrix factorization in C++

Hello !

I'm very beginner in numerical methods... I have the pseudo-code
(look at image), but I don't know how intrepret this and how
I can implement this in C++...

Can You help me ?

Offline

#2 2012-12-25 05:50:16

bobbym
bumpkin
From: Bumpkinland
Registered: 2009-04-12
Posts: 109,606

Re: Gram-Schmidt matrix factorization in C++

Hi;

I need a little bit more than that.

#include <iostream>
#include <glm/glm.hpp>

glm::vec3 sum_over_e(glm::vec3* e, glm::vec3* e_prime, int& i)
{
    int k = 0;
    glm::vec3 result;

    while (k < i-1)
    {
        float e_prime_k_squared = glm::dot(e_prime[k], e_prime[k]);
        result += ((glm::dot(e[i], e_prime[k]) / e_prime_k_squared) * e_prime[k]);
        k++;
    }

    return result;
}

int main(int argc, char** argv)
{
    int n = 3;  // number of vectors we're working with
    glm::vec3 e[] = {
        glm::vec3(sqrt(2)/2, sqrt(2)/2, 0),
        glm::vec3(-1, 1, -1),
        glm::vec3(0, -2, -2)
    };

    glm::vec3 e_prime[n];
    e_prime[0] = e[0];  // step A

    int i = 0;  // step B

    do  // step C
    {
        e_prime[i] = e[i] - sum_over_e(e, e_prime, i);

        i++;    // step D
    } while (i < n);

    for (int loop_count = 0; loop_count < n; loop_count++)
    {
        std::cout << "Vector e_prime_" << loop_count+1 << ": < " 
                  << e_prime[loop_count].x << ", " 
                  << e_prime[loop_count].y << ", " 
                  << e_prime[loop_count].z << " >" << std::endl;
    }

    return 0;

That is supposed to orthogonalize those three vectors using Gram Schmidt. I have not tried it but it is supposed to work.

I must point out that the above method can be numerically unstable and the modified method will produce better results.


In mathematics, you don't understand things. You just get used to them.
If it ain't broke, fix it until it is.
Always satisfy the Prime Directive of getting the right answer above all else.

Offline

#3 2013-01-06 23:21:28

Giewont
Member
Registered: 2012-12-25
Posts: 3

Re: Gram-Schmidt matrix factorization in C++

I found a solution:

    int k, i, j;
    for (k=0; k<3; k++)
   {
      r[k][k]=0; // equivalent to sum = 0
      for (i=0; i<3; i++)
        r[k][k] = r[k][k] + a[i][k] * a[i][k]; // rkk = sqr(a0k) + sqr(a1k) + sqr(a2k) 
      r[k][k] = sqrt(r[k][k]);  // ||a||
      
      for (i=0; i<3; i++)
          q[i][k] = a[i][k]/r[k][k];
     
      for(j=k+1; j<3; j++) {
        r[k][j]=0;
        for(i=0; i<3; i++) r[k][j] += q[i][k] * a[i][j];
        for (i=0; i<3; i++) a[i][j] = a[i][j] - r[k][j]*q[i][k];
        
      }

I have tested this and it's working good.
But now I must parallelize this code using OpenMP (it's the project for my college) and
I have a new problem. I try parallelize the main "for":

int k, i, j;
#pragma omp parallel for private (k, i, j) shared (a, q, r)
// ........

but it doesn't work correctly. I noticed that the problematic fragment is:

for (i=0; i<3; i++) 
          q[i][k] = a[i][k]/r[k][k];

but I don't know why it makes a problem...

Anybody help ?

Offline

#4 2013-01-07 06:53:27

bobbym
bumpkin
From: Bumpkinland
Registered: 2009-04-12
Posts: 109,606

Re: Gram-Schmidt matrix factorization in C++

Hi;

What was the error message the compiler gave?


In mathematics, you don't understand things. You just get used to them.
If it ain't broke, fix it until it is.
Always satisfy the Prime Directive of getting the right answer above all else.

Offline

#5 2013-01-09 04:32:12

Giewont
Member
Registered: 2012-12-25
Posts: 3

Re: Gram-Schmidt matrix factorization in C++

bobbym wrote:

Hi;

What was the error message the compiler gave?

There aren't compiler's errors, but the application produces bad, meaningless results.
(completely different than non-parallel version) sad

Offline

#6 2013-01-09 04:39:54

bobbym
bumpkin
From: Bumpkinland
Registered: 2009-04-12
Posts: 109,606

Re: Gram-Schmidt matrix factorization in C++

Hi;

That is weird. Did you try to use the code I provided?


In mathematics, you don't understand things. You just get used to them.
If it ain't broke, fix it until it is.
Always satisfy the Prime Directive of getting the right answer above all else.

Offline

Board footer

Powered by FluxBB