1 #ifdef HAVE_STDLIB_H
2 #include <stdlib.h>
3 #endif
4 #include <math.h>
5 #define NRANSI
6 #include "nrutil.h"
7 /* #include "yagi.h" */
8 #define TINY 1.0e-20;
9
ludcmp(double ** a,int n,int * indx,double * d)10 void ludcmp(double **a, int n, int *indx, double *d)
11 {
12 int i,imax,j,k;
13 double big,dum,sum,temp;
14 double *vv;
15
16 vv=dvector(1,n);
17 *d=1.0;
18 for (i=1;i<=n;i++) {
19 big=0.0;
20 for (j=1;j<=n;j++)
21 if ((temp=fabs(a[i][j])) > big) big=temp;
22 if (big == 0.0) nrerror("Singular matrix in routine ludcmp");
23 vv[i]=1.0/big;
24 }
25 for (j=1;j<=n;j++) {
26 for (i=1;i<j;i++) {
27 sum=a[i][j];
28 for (k=1;k<i;k++) sum -= a[i][k]*a[k][j];
29 a[i][j]=sum;
30 }
31 big=0.0;
32 for (i=j;i<=n;i++) {
33 sum=a[i][j];
34 for (k=1;k<j;k++)
35 sum -= a[i][k]*a[k][j];
36 a[i][j]=sum;
37 if ( (dum=vv[i]*fabs(sum)) >= big) {
38 big=dum;
39 imax=i;
40 }
41 }
42 if (j != imax) {
43 for (k=1;k<=n;k++) {
44 dum=a[imax][k];
45 a[imax][k]=a[j][k];
46 a[j][k]=dum;
47 }
48 *d = -(*d);
49 vv[imax]=vv[j];
50 }
51 indx[j]=imax;
52 if (a[j][j] == 0.0) a[j][j]=TINY;
53 if (j != n) {
54 dum=1.0/(a[j][j]);
55 for (i=j+1;i<=n;i++) a[i][j] *= dum;
56 }
57 }
58 free_dvector(vv,1,n);
59 }
60 #undef TINY
61 #undef NRANSI
62 /* (C) Copr. 1986-92 Numerical Recipes Software 5#,. */
63