title: LU Decomposition from Numerical Recipes not working; what am I doing wrong?
body: <p>I've literally copied and pasted from the supplied source code for Numerical Recipes for C for in-place LU Matrix Decomposition, problem is its not working. </p>
  
  <p>I'm sure I'm doing something stupid but would appreciate anyone being able to point me in the right direction on this; I've been working on its all day and can't see what I'm doing wrong.</p>
  
  <p>POST-ANSWER UPDATE: The project is <a href="http://www.andrewbolster.info/blog/2011/04/lu-decomposition-in-c-and-under-cuda/" rel="nofollow">finished and working</a>. Thanks to everyone for their guidance.</p>
  
  <pre><code>#include &lt;stdlib.h&gt;
  #include &lt;stdio.h&gt;
  #include &lt;math.h&gt;
  #define MAT1 3
  #define TINY 1e-20
  
  int h_NR_LU_decomp(float *a, int *indx){
    //Taken from Numerical Recipies for C
    int i,imax,j,k;
    float big,dum,sum,temp;
    int n=MAT1;
  
    float vv[MAT1];
    int d=1.0;
    //Loop over rows to get implicit scaling info
    for (i=0;i&lt;n;i++) {
      big=0.0;
      for (j=0;j&lt;n;j++)
        if ((temp=fabs(a[i*MAT1+j])) &gt; big) 
          big=temp;
      if (big == 0.0) return -1; //Singular Matrix
      vv[i]=1.0/big;
    }
    //Outer kij loop
    for (j=0;j&lt;n;j++) {
      for (i=0;i&lt;j;i++) {
        sum=a[i*MAT1+j];
        for (k=0;k&lt;i;k++) 
          sum -= a[i*MAT1+k]*a[k*MAT1+j];
        a[i*MAT1+j]=sum;
      }
      big=0.0;
      //search for largest pivot
      for (i=j;i&lt;n;i++) {
        sum=a[i*MAT1+j];
        for (k=0;k&lt;j;k++) sum -= a[i*MAT1+k]*a[k*MAT1+j];
        a[i*MAT1+j]=sum;
        if ((dum=vv[i]*fabs(sum)) &gt;= big) {
          big=dum;
          imax=i;
        }
      }
      //Do we need to swap any rows?
      if (j != imax) {
        for (k=0;k&lt;n;k++) {
          dum=a[imax*MAT1+k];
          a[imax*MAT1+k]=a[j*MAT1+k];
          a[j*MAT1+k]=dum;
        }
        d = -d;
        vv[imax]=vv[j];
      }
      indx[j]=imax;
      if (a[j*MAT1+j] == 0.0) a[j*MAT1+j]=TINY;
      for (k=j+1;k&lt;n;k++) {
        dum=1.0/(a[j*MAT1+j]);
        for (i=j+1;i&lt;n;i++) a[i*MAT1+j] *= dum;
      }
    }
    return 0;
  }
  
  void main(){
      //3x3 Matrix
      float exampleA[]={1,3,-2,3,5,6,2,4,3};
      //pivot array (not used currently)
      int* h_pivot = (int *)malloc(sizeof(int)*MAT1);
      int retval = h_NR_LU_decomp(&amp;exampleA[0],h_pivot);
      for (unsigned int i=0; i&lt;3; i++){
        printf("\n%d:",h_pivot[i]);
        for (unsigned int j=0;j&lt;3; j++){
          printf("%.1lf,",exampleA[i*3+j]);
        }
      }
  }
  </code></pre>
  
  <p><a href="http://www.wolframalpha.com/input/?i=LU%20decomposition%20of%20%7B%7B1,3,-2%7D,%7B3,5,6%7D,%7B2,4,3%7D%7D" rel="nofollow">WolframAlpha</a> says the answer should be</p>
  
  <pre><code>1,3,-2
  2,-2,7
  3,2,-2
  </code></pre>
  
  <p>I'm getting:</p>
  
  <pre><code>2,4,3
  0.2,2,-2.8
  0.8,1,6.5
  </code></pre>
  
  <p>And so far I have found at least 3 different versions of the 'same' algorithm, so I'm completely confused.</p>
  
  <p>PS yes I know there are at least a dozen different libraries to do this, but I'm more interested in understanding what I'm doing wrong than the right answer.</p>
  
  <p>PPS since in LU Decomposition the lower resultant matrix is unity, and using Crouts algorithm as (i think) implemented, array index access is still safe, both L and U can  be superimposed on each other in-place; hence the single resultant matrix for this.  </p>