refactoring

graph/new/nr_lautil.c

@@ -0,0 +1,113 @@
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include "nrutil.h"
+#define TINY 1.0e-20
+// #define N 2
+
+
+void ludcmp (float **a, int n, int *indx, float *d)
+{
+  int i, imax, j, k;
+  float big, dum, sum, temp;
+  float *vv;
+
+  vv = vector (1,n);
+  *d = 1.0;
+  for (i=1; i<=n; i++) 
+  {
+    big = 0.0;
+    for (j=1; j<=n; j++)
+      if ((temp = fabs(a[i][j])) > big) big = temp;
+    if (big == 0.0) nrerror("Singular matrix in routine ludcmp");
+    vv[i] = 1.0 / big;
+  }
+  
+  for (j=1; j<=n; j++) 
+  {
+    for (i=1; i<j; i++) 
+    {
+      sum = a[i][j];
+      for (k=1; k<i; k++) sum -= a[i][k] * a[k][j];
+      a[i][j] = sum;
+    }
+    big = 0.0;
+    for (i=j; i<=n; i++) 
+    {
+      sum = a[i][j];
+      for (k=1; k<j; k++)
+        sum -= a[i][k] * a[k][j];
+      a[i][j] = sum;
+      if ((dum = vv[i] * fabs(sum)) >= big) 
+      {
+        big = dum;
+        imax = i;
+      }
+    }
+    if (j != imax) 
+    {
+      for (k=1; k<=n; k++) 
+      {
+        dum = a[imax][k];
+        a[imax][k] = a[j][k];
+        a[j][k] = dum;
+      }
+      *d = -(*d);
+      vv[imax] = vv[j];
+    }
+    indx[j] = imax;
+    if (a[j][j] == 0.0) a[j][j] = TINY;
+    if (j != n) 
+    {
+      dum = 1.0 / a[j][j];
+      for (i=j+1; i<=n; i++) a[i][j] *= dum;
+    }
+  }
+  free_vector (vv, 1, n);
+}
+
+
+void lubksb (float **a, int n, int *indx, float b[])
+{
+  int i, ii=0, ip, j;
+  float sum;
+  
+  for (i=1; i<=n; i++)
+  {
+    ip = indx[i];
+    sum = b[ip];
+    b[ip] = b[i];
+    if (ii)
+      for (j=ii; j<=i-1; j++) sum -= a[i][j] * b[j];
+    else if (sum) ii = i;
+    b[i] = sum;
+  }
+  for (i=n; i>=1; i--) 
+  {
+    sum = b[i];
+    for (j=i+1; j<=n; j++) sum -= a[i][j] * b[j];
+    b[i] = sum / a[i][i];
+  }
+}
+
+
+// int main ()
+// {
+//   float **a, d;
+//   int i, j, *indx;
+//   
+//   float a_array[N][N] = {
+//     { 1.0, 2.0 },
+//     { 3.0, 4.0 }
+//   };
+//   a = matrix (1, N, 1, N);
+//   for (i=1; i<N+1; i++)
+//     for (j=1; j<N+1; j++)
+//       a[i][j] = a_array[i-1][j-1];
+//   indx = ivector (1, N);
+//   ludcmp (a, N, indx, &d);
+//   for(j=1; j<=N; j++) d *= a[j][j];
+//   printf ("%lf\n", d);
+// 
+//   return 0;
+// }