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
source to share
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;
source to share