Math Is Fun Forum

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 ?

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.

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 ?

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.

Giewont

Member

Registered: 2012-12-25

Posts: 3

Re: Gram-Schmidt matrix factorization in C++

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)

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.