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 ?

View Image: DSCN7410 - Kopia.JPG

Offline

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

bobbym
Administrator
From: Bumpkinland
Registered: 2009-04-12
Posts: 83,211

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.
I have the result, but I do not yet know how to get it.
All physicists, and a good many quite respectable mathematicians are contemptuous about proof.

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
Administrator
From: Bumpkinland
Registered: 2009-04-12
Posts: 83,211

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.
I have the result, but I do not yet know how to get it.
All physicists, and a good many quite respectable mathematicians are contemptuous about proof.

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
Administrator
From: Bumpkinland
Registered: 2009-04-12
Posts: 83,211

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.
I have the result, but I do not yet know how to get it.
All physicists, and a good many quite respectable mathematicians are contemptuous about proof.

Offline

Board footer

Powered by FluxBB