Possible loss of precision with Gram-Schmidt

I have some code that uses Gram-Schmidt inside a loop. I want to reduce the number of calls to this algorithm as much as possible, but the fact is that despite getting the same result before and after the call, when I print the results of some operations using these values, they are different. For example, in the code below, the result abs(muGS[k][0]) - abs(before2)

should be 0 or very close to 0, since the printable values ​​of this variable (before and after the call) are the same. However, this is not what is happening. muGS

is a double matrix, and its values ​​are usually between 0 and 1.

int k = 1;
double before2;

while(k < end) {

    before2 = muGS[k][0];

    gramSchmidt(b, muGS, cGS, k);

    //prints for debug
    if (abs(muGS[k][0]) - abs(before2) > 0.1) {

        if (abs(muGS[k][0]) - abs(before2) > 0.1)  {
            cout << "1 muGS[k] diff:" << abs(muGS[k][0]) - abs(before2) << endl;
            cout << "1 muGS[k] before:" << muGS[k][0] << endl;
            cout << "1 muGS[k] after:" << muGS[k][0] << endl;
            cout << "1 muGS[k] mult before:" << before2 * before2 << endl;
            cout << "1 muGS[k] mult after:" << muGS[k][0] * muGS[k][0] << endl;
            cout << "1 muGS[k] abs before:" << abs(before2) << endl;
            cout << "1 muGS[k] abs after:" << abs(muGS[k][0]) << endl;
        }
        getchar();
    }

    for (i = k-1; i >= 0; i--) {
        for (j = 0; j < i; j++) {
            muGS[k][j] -= round(muGS[k][i]) * muGS[i][j];
        }
    }

    //some other operations that don't change the value of muGS
    k++;
}

      

Output:

1 muGS[k] diff:0.157396
1 muGS[k] before:0.288172
1 muGS[k] after:0.288172
1 muGS[k] mult before:0.0171023
1 muGS[k] mult after:0.083043
1 muGS[k] abs before:0.130776
1 muGS[k] abs after:0.288172

      

Another thing is that the absolute value is before2

very different from the value before2

. Is it possible that I have some loss of precision or why is this happening?

thank

+3


source to share


1 answer


There are no exact losses. You are simply wrong in your code:

        cout << "1 muGS[k] before:" << muGS[k][0] << endl;
        cout << "1 muGS[k] after:" << muGS[k][0] << endl;

      



You print the same value both before and after. But it should be:

        cout << "1 muGS[k] before:" << before2 << endl;
        cout << "1 muGS[k] after:" << muGS[k][0] << endl;

      

+4


source







All Articles