Gram-Schmidt matrix factorization

[Giewont]

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 ?

image of pseudo-code:
http://postimage.org/image/o6sitd3l9/

[hamsterman]

Factorization? Gram-Schmidt is for orthonormalizing a basis. Start with a (numbered) set of vectors. Take the first one, normalize it (divide it by it's length) then fix the remaining ones so that their scalar product with the first one is 0 (using the last formula in the picture). Repeat for the remaining vectors.

If unclear, the algorithm should be well explained in wikipedia.

[FredFlintstone]

I looked for a numerical example of the method.

I wrote the C program and I have verified that the results correspond to those of the example.

#include <cstdlib>
#include <iostream>
#include <math.h>
using namespace std;

// example: http://www.mia.uni-saarland.de/Teaching/NAVC-SS11/sol_c8.pdf
// page 5

double a[3][3] = {
    {1.0, 2.0, 1.0},
    {0.0, 1.0, 2.0},
    {1.0, 2.0, 0.0}
};
// any column of a is a vector

double r[3][3], q[3][3];

int main(int argc, char *argv[]) {
    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||
      cout << endl << "R"<<k<<k<<": " << r[k][k];
      
      for (i=0; i<3; i++) {
          q[i][k] = a[i][k]/r[k][k];
          cout << " q"<<i<<k<<": "<<q[i][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];
        cout << endl << "r"<<k<<j<<": " <<r[k][j] <<endl;
        
        for (i=0; i<3; i++) a[i][j] = a[i][j] - r[k][j]*q[i][k];
        
        for (i=0; i<3; i++) cout << "a"<<j<<": " << a[i][j]<< " ";
      }
}

    system("PAUSE");
    return EXIT_SUCCESS;
}

[Giewont]

I wrote the C program and I have verified that the results correspond to those of the example.

Thank You very much ! It works really good :)

I have one more, additional question: do You know (or somebody...) how can I convert it to parallel version using OpenMP ?

I tried to use:

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...

[hamsterman]

How would you parallelize this? It's a pretty sequential algorithm. I guess if you have a ridiculous number of dimensions, it might be doable...

[Giewont]

It's an exercise for my college. I must parallelize the sequential Gram-Schmidt's algorithm and compare execution time of both version (of course for large matrix).