1 /* ludecomp - LU decomposition and backsolving routines. */
2 /* Taken from Numerical Recipies. */
3 /* XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney */
4 /* Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz */
5 /* You may give out copies of this software; for conditions see the */
6 /* file COPYING included with this distribution. */
11 crludcmp(Matrix mat
, int n
, IVector indx
, int mode
, double *d
)
13 int i
, imax
, j
, k
, singular
= FALSE
;
17 CMatrix cmat
= (CMatrix
) mat
;
18 RMatrix rmat
= (RMatrix
) mat
;
24 /* set up the pivot permutation vector */
25 for (i
= 0; i
< n
; i
++) indx
[i
] = i
;
27 /* get scaling information for implicit pivoting */
28 for (i
= 0; i
< n
; i
++) {
30 for (j
= 0; j
< n
; j
++) {
31 temp
= (mode
== RE
) ? fabs(rmat
[i
][j
]) : modulus(cmat
[i
][j
]);
32 if (temp
> big
) big
= temp
;
35 vv
[i
] = 1.0; /* no scaling for zero rows */
38 else vv
[i
] = 1.0 / big
;
41 /* loop over columns for Crout's method */
42 for (j
= 0; j
< n
; j
++) {
43 for (i
= 0; i
< j
; i
++) {
44 if (mode
== RE
) rsum
= rmat
[i
][j
];
45 else csum
= cmat
[i
][j
];
47 for (k
= 0; k
< i
; k
++)
48 if (mode
== RE
) rsum
-= rmat
[i
][k
] * rmat
[k
][j
];
49 else csum
= csub(csum
, cmul(cmat
[i
][k
], cmat
[k
][j
]));
51 if (mode
== RE
) rmat
[i
][j
] = rsum
;
52 else cmat
[i
][j
] = csum
;
55 for (i
= j
; i
< n
; i
++) {
56 if (mode
== RE
) rsum
= rmat
[i
][j
];
57 else csum
= cmat
[i
][j
];
59 for (k
= 0; k
< j
; k
++)
60 if (mode
== RE
) rsum
-= rmat
[i
][k
] * rmat
[k
][j
];
61 else csum
= csub(csum
, cmul(cmat
[i
][k
], cmat
[k
][j
]));
63 if (mode
== RE
) rmat
[i
][j
] = rsum
;
64 else cmat
[i
][j
] = csum
;
66 temp
= vv
[i
] * ((mode
== RE
) ? fabs(rsum
) : modulus(csum
));
67 if (temp
>= big
) { big
= temp
; imax
= i
; }
70 /* interchange rows if needed */
72 for (k
= 0; k
< n
; k
++) {
75 rmat
[imax
][k
] = rmat
[j
][k
];
80 cmat
[imax
][k
] = cmat
[j
][k
];
89 /* divide by the pivot element */
90 temp
= (mode
== RE
) ? fabs(rmat
[j
][j
]) : modulus(cmat
[j
][j
]);
91 if (temp
== 0.0) singular
= TRUE
;
94 rdum
= 1.0 / rmat
[j
][j
];
95 for (i
= j
+ 1; i
< n
; i
++) rmat
[i
][j
] *= rdum
;
98 cdum
= cdiv(cart2complex(1.0, 0.0), cmat
[j
][j
]);
99 for (i
= j
+ 1; i
< n
; i
++) cmat
[i
][j
] = cmul(cmat
[i
][j
], cdum
);
108 crlubksb(double **a
, int n
, int *indx
, double *b
, int mode
)
110 int i
, ii
, ip
, j
, singular
= FALSE
;
111 CMatrix ca
= (CMatrix
) a
;
112 CVector cb
= (CVector
) b
;
113 RMatrix ra
= (RMatrix
) a
;
114 RVector rb
= (RVector
) b
;
118 /* forward substitute using L part */
119 for (i
= 0, ii
= -1; i
< n
; i
++) {
130 for (j
= ii
; j
<= i
- 1; j
++)
131 if (mode
== RE
) rsum
-= ra
[i
][j
] * rb
[j
];
132 else csum
= csub(csum
, cmul(ca
[i
][j
], cb
[j
]));
135 if (rsum
!= 0.0) ii
= i
;
137 else if (csum
.real
!= 0.0 || csum
.imag
!= 0.0) ii
= i
;
139 if (mode
== RE
) rb
[i
] = rsum
;
143 /* back substitute using the U part */
144 for (i
= n
- 1; i
>= 0; i
--) {
147 for (j
= i
+ 1; j
< n
; j
++) rsum
-= ra
[i
][j
] * rb
[j
];
148 if (ra
[i
][i
] == 0.0) {
152 else rb
[i
] = rsum
/ ra
[i
][i
];
156 for (j
= i
+ 1; j
< n
; j
++) csum
= csub(csum
, cmul(ca
[i
][j
], cb
[j
]));
157 if (modulus(ca
[i
][i
]) == 0.0) {
161 else cb
[i
] = cdiv(csum
, ca
[i
][i
]);