You are not logged in.

- Topics: Active | Unanswered

Pages: **1**

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

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

**bobbym****Administrator**- From: Bumpkinland
- Registered: 2009-04-12
- Posts: 107,190

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

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

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

**bobbym****Administrator**- From: Bumpkinland
- Registered: 2009-04-12
- Posts: 107,190

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

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

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)

Offline

**bobbym****Administrator**- From: Bumpkinland
- Registered: 2009-04-12
- Posts: 107,190

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

Pages: **1**