user5887077 user5887077 - 17 days ago 5
C Question

Jacobi method works with double, fails with float. What is wrong?

I coded a little program to solve a system of n equations with the Jacobi (iterational) method. Below is the code:

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
int main() {

float *a, *b, *x, *xnew, temp;
int i, j, k, maxiter=10000000, n=4;

a = malloc(n*n*sizeof(*a));
b = malloc(n*sizeof(*b));
x = malloc(n*sizeof(*x));
xnew = malloc(n*sizeof(*xnew));

srand((unsigned) time(NULL));

// Filling the matrix
for (i=0;i<=n-1;i++) {
for (j=0;j<=n-1;j++) {
a[n*i+j] = rand()%60;
}
b[i] = rand();
x[i] = rand();
xorg[i]=x[i];
}

// Establishing diagonal dominance
for (i=0;i<=n-1;i++) {
temp=0;
for (j=0;j<=n-1;j++) {
if (j==i) {continue;}
temp = temp + a[n*i+j];
}
a[n*i+i] = temp+1;
}

// Solve the system. Break when residue is low
for (k=0;k<=maxiter-1;k++) {
for (i=0;i<=n-1;i++) {
temp=0;
for (j=0;j<=n-1;j++) {
if (j==i) {continue;}
temp = temp + a[n*i+j]*x[j];
}
xnew[i] = (b[i]-temp)/a[n*i+i];
}
temp=0;
for (i=0;i<=n-1;i++) {
temp = temp + fabs(x[i]-xnew[i]);
x[i]=xnew[i];
}
if (temp<0.0001) {
break;
}
}

printf("Iterations = %d\n",k-1);

return 0;
}


The break-out-of-the-loop criterion is ridiculously easy. This program should never fail. Yet it apparently does not converge (it uses up all iterations in the loop), UNLESS I change the floats to doubles. Floats have much greater precision than this. What is wrong?
Compiled with CodeBlocks 16.01 under Windows 7, if that even matters.

Answer

It appears that the float datatype does not have the required precision for the above algorithm, given the way is was coded. The algorithm does converge, but the "residue" will never be low enough to exit the loop.

The way I understand this is that due to the way float variables are stored internally, you cannot do computation with extremely small (0.0001) and extremely large (RAND_MAX) numbers and expect reasonable accuracy, as in the above example (temp grows to a huge number in the innermost loop).

Consequently, setting b[i] = rand()%60; and x[i] = rand()%60; will alleviate the problem.

Setting b[i] = rand()%6; x[i] = rand()%6; and a[n*i+j] = rand()%6 will allow an even tighter exit-loop condition to be eventually met.

Interestingly, establishing greater diagonal dominance (change a[n*i+i] = temp+1' to a[n*i+i] = temp+10; will also make the program converge, when it would previously not.

I was not familiar with the ULP conditions described by others, but will invest some time into that

If future readers of this have the time and energy, maybe they should read "What Every Computer Scientist Should Know About Floating-Point Arithmetic", even though I have not.

BTW, xorg was to store the original x vector, for debugging purposes, because I was having a terribly hard time getting CodeBlocks to debug

Thanks to all for contributing.

Comments