1    /*
2    Copyright 2000 by Ralph Hartley
3    This software is licenced under the terms of the
4    Gnu Public Licence
5    */
6    import java.lang.Math;
7    
8    class Linset {
9      public double A[][];
10     public double b[];
11     public int n;
12     int indx[];
13   
14     static final double TINY = 1.0e-20;
15     static boolean DEBUG = false;
16   
17     public Linset(int N) {
18           n=N;
19           A = new double[N][N];
20           b = new double[N];
21           indx = new int[N];
22     }    
23     public void lubksb() {
24           int i,j,ip;
25           int ii= -1;
26           double sum;
27   
28           for (i=0;i<n;i++) {
29             ip = indx[i];
30             sum = b[ip];
31             b[ip] = b[i];
32             if (ii != -1)
33           for (j=ii;j<=i-1;j++) sum -= A[i][j]*b[j];
34             else if (sum!=0) ii=i;
35             b[i] = sum;
36           }
37           for (i=n-1;i>=0;i--) {
38             sum = b[i];
39             for (j=i+1;j<n;j++) sum -= A[i][j]*b[j];
40             b[i] = sum/A[i][i];
41           }
42     }        
43     public void ludcmp() {
44           int i,j,k;
45           int imax=0;
46           double big,dum,sum;
47           double[] vv = new double[n];
48           double[] temp;
49   
50           if (DEBUG) {
51             System.out.println("N = "+n);
52             for (i=0;i<n;i++) {
53               String str = "";
54               for (j=0;j<n;j++) str = str + " " + A[i][j];
55               System.out.println(str);
56             }
57           }
58   
59           for (i=0;i<n;i++) {
60             big = 0;
61             for (j=0;j<n;j++)
62           if (Math.abs(A[i][j]) > big) big = Math.abs(A[i][j]);
63             if (big==0.0) ;//@@@throw singular
64             vv[i] = 1.0/big;
65           }
66   
67           for (j=0;j<n;j++) {
68             for (i=0;i<j;i++) {
69           sum = A[i][j];
70           for (k=0;k<i;k++) sum -= A[i][k]*A[k][j];
71           A[i][j] = sum;
72             }
73             big = 0.0;
74             for (i=j;i<n;i++) {
75           sum = A[i][j];
76           for (k=0;k<j;k++) sum -= A[i][k]*A[k][j];
77           A[i][j] = sum;
78           if (vv[i]*Math.abs(sum)>=big) {
79             big = vv[i]*Math.abs(sum);
80             imax = i;
81           }
82             }
83             if (j!=imax) {
84           temp = A[imax];
85           A[imax] = A[j];
86           A[j] = temp;
87           vv[imax] = vv[j];
88             }
89             indx[j] = imax;
90             if (A[j][j] == 0) A[j][j] = TINY;
91             if (j!=n-1) {
92           dum = 1.0/(A[j][j]);
93           for (i=j+1;i<n;i++) A[i][j] *= dum;
94             }
95           }
96     }    
97   }
98