Gram Schmidt implementation works partially, outputs a bad number

Question

I'm trying to implement Gram Schmidt in code for a school project and I have a problem.It outputs a weird number and I don't know why.The rest of them are good but this one is bad.Sorry if the post is bad but it's my first post here and I really need help.Thanks

#include <iostream>
#define dim 100

using namespace std;



void inputMatrix(int *n, int *m, double x[dim][dim]){
    int i,j;
    for(i = 0; i < *n; i++){
        cout << "V" << i << ":";
        for(j = 0; j < *m; j++)
            cin >> x[i][j];
    }

}

void outputMatrix(int *n, int *m, double x[dim][dim]){
    int i,j;
    for(i=0;i<*n;i++){
        for(j=0;j<*m;j++)
            cout<<x[i][j]<<"  ";
        cout<<endl;
    }
}
void initialize(int *m,double v[dim]){
    int i;
    for(i=0;i<*m;i++){
         v[i]=0;
    }    
}
int main(){
    double v[dim][dim], f[dim][dim], e[dim][dim],p1,p2,a[dim];
    int n,m,i,j,z;
    cout << "Introduceti numarul de vectori: ";cin >> n;
    cout << "Introduceti numarul de elemente: ";cin >> m;
    inputMatrix(&n,&m,v);
    double div = 0;
    for(i = 0; i < m; i++){
        f[0][i] = v[0][i];
    }
    outputMatrix(&n,&m,v);
    cout << endl;
    outputMatrix(&n,&m,f);
    for(i = 1;i < n; i++){
        z = 0;
        initialize(&m,a);
        mk1:
        p1 = 0;
        p2 = 0;
        for(j = 0; j < m; j++){
            p1 = f[z][j] * v[i][j] + p1;
            p2 = f[z][j] * f[z][j] + p2;
        }
        div = p1 / p2;
        for(j = 0; j < m; j++){
            a[j] = f[z][j] * div + a[j];
        }
        z++;
        if( z < i){
            goto mk1;
        }
        else{
            for(j = 0; j < m;j++){
                f[i][j] = v[i][j] - a[j];
            }
        }
    cout << endl;
    outputMatrix(&n,&m,f);
    return 0;
}

Output:

Introduceti numarul de vectori: 3
Introduceti numarul de elemente: 4
V0:1
2
3
0
V1:1
2
0
0
V2:1
0
0
1
1  2  3  0  
1  2  0  0  
1  0  0  1  

1  2  3  0  
0  0  0  0  
0  0  0  0  

1  2  3  0  
0.642857  1.28571  -1.07143  0  
0.8  -0.4  5.55112e-17  1

I don't understand why it outputs this "5.55112e-17" Thanks for help!

Thanks for the answer, the correct output should be 0 there, so that's an artifact of calculating, but it can be resolved ? — Viorel Andrei Chis, Mar 20 '18 at 12:20
Instead of `#define dim 100`, use a const, i.e., `const int dim_matrix = 100`. — Xam, Mar 20 '18 at 17:17

score 2 · Accepted Answer · answered Mar 20 '18 at 12:15

2

I don't understand why it outputs this "5.55112e-17" Thanks for help!

This is the result of rounding errors, this notation means 5.55112 times 10 to the power of -17 and can be considered an artifact of calculating with computer representations of (real) numbers. For some more elaboration, consider this answer.

answered Mar 20 '18 at 12:15

Jodocus

7,493
1
29
45

I understand the problem here, the thing is that it should output 0 instead of that, it is possible somehow to fix it ? – Viorel Andrei Chis Mar 20 '18 at 12:22
1

@ViorelAndreiChis Of course you could just purge every number below a certain threshold to be zero: `if(abs(n) < 1.e-14) n = 0.;` This might produce exact results here, but may be problematic for problem where e.g. the vectors are near linear dependent. – Jodocus Mar 20 '18 at 12:38
Thanks for help, now it works just fine, thank you again! – Viorel Andrei Chis Mar 20 '18 at 12:49

score 2 · Answer 2 · answered Mar 20 '18 at 13:03

5.55112e-17 is very close to zero. It is not exactly zero, due to the numerical resolution of your used binary representation of numbers (e.g. double).

You can use a higher resolution for the calculation than for the printing of the results. Alternatively, you can round at the end. Nevertheless, you can not get rid of the numerics.

Here a few other remarks on your code:

Please use proper containers, here std::vector. This also enables you to use dynamic memory allocation (Your code will crash if the matrix becomes bigger than 100).
Please dont use goto. Below you find an alternative implementation.
It is not meaningful to hand pointers to integers to a function. You don't gain anything from that, but the readability is reduced.

Code without goto

int main() {
    double v[dim][dim], f[dim][dim], e[dim][dim], p1, p2, a[dim];
    int n, m, i, j, z;
    cout << "Introduceti numarul de vectori: ";
    cin >> n;
    cout << "Introduceti numarul de elemente: ";
    cin >> m;
    inputMatrix(&n, &m, v);
    double div = 0;
    for (i = 0; i < m; i++) {
        f[0][i] = v[0][i];
    }
    outputMatrix(&n, &m, v);
    cout << endl;
    outputMatrix(&n, &m, f);
    for (i = 1; i < n; i++) {
        z = 0;
        initialize(&m, a);
        do {
            p1 = 0;
            p2 = 0;
            for (j = 0; j < m; j++) {
                p1 = f[z][j] * v[i][j] + p1;
                p2 = f[z][j] * f[z][j] + p2;
            }
            div = p1 / p2;
            for (j = 0; j < m; j++) {
                a[j] = f[z][j] * div + a[j];
            }
            z++;
        } while (z < i);
        for (j = 0; j < m; j++) {
            f[i][j] = v[i][j] - a[j];
        }
    }
    cout << endl;
    outputMatrix(&n, &m, f);
    return 0;
}

For matrices I would not make use `std::vector` at all: matrices are not resizeable. I agree on the dynamic memory part, though. Using a matrix with a memory size close to the number of elements will benefit the processor's cache performance of the program. — Jorge Bellon, Mar 20 '18 at 13:11

Gram Schmidt implementation works partially, outputs a bad number

2 Answers2