title: Jacobi method works with double, fails with float. What is wrong?
body: <p>I coded a little program to solve a system of n equations with the Jacobi (iterational) method. Below is the code:</p>
  
  <pre><code>#include &lt;math.h&gt;
  #include &lt;stdio.h&gt;
  #include &lt;stdlib.h&gt;
  #include &lt;string.h&gt;
  #include &lt;time.h&gt;
  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&lt;=n-1;i++) {
      for (j=0;j&lt;=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&lt;=n-1;i++) {
      temp=0;
      for (j=0;j&lt;=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&lt;=maxiter-1;k++) {
      for (i=0;i&lt;=n-1;i++) {
          temp=0;
          for (j=0;j&lt;=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&lt;=n-1;i++) {
          temp = temp + fabs(x[i]-xnew[i]);
          x[i]=xnew[i];
      }
      if (temp&lt;0.0001) {
          break;
      }
  }
  
  printf("Iterations = %d\n",k-1);
  
  return 0;
  }
  </code></pre>
  
  <p>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.</p>