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. */
10 crludcmp(mat
, n
, indx
, mode
, d
)
16 int i
, imax
, j
, k
, singular
= FALSE
;
20 CMatrix cmat
= (CMatrix
) mat
;
21 RMatrix rmat
= (RMatrix
) mat
;
27 /* set up the pivot permutation vector */
28 for (i
= 0; i
< n
; i
++) indx
[i
] = i
;
30 /* get scaling information for implicit pivoting */
31 for (i
= 0; i
< n
; i
++) {
33 for (j
= 0; j
< n
; j
++) {
34 temp
= (mode
== RE
) ? fabs(rmat
[i
][j
]) : modulus(cmat
[i
][j
]);
35 if (temp
> big
) big
= temp
;
38 vv
[i
] = 1.0; /* no scaling for zero rows */
41 else vv
[i
] = 1.0 / big
;
44 /* loop over columns for Crout's method */
45 for (j
= 0; j
< n
; j
++) {
46 for (i
= 0; i
< j
; i
++) {
47 if (mode
== RE
) rsum
= rmat
[i
][j
];
48 else csum
= cmat
[i
][j
];
50 for (k
= 0; k
< i
; k
++)
51 if (mode
== RE
) rsum
-= rmat
[i
][k
] * rmat
[k
][j
];
52 else csum
= csub(csum
, cmul(cmat
[i
][k
], cmat
[k
][j
]));
54 if (mode
== RE
) rmat
[i
][j
] = rsum
;
55 else cmat
[i
][j
] = csum
;
58 for (i
= j
; i
< n
; i
++) {
59 if (mode
== RE
) rsum
= rmat
[i
][j
];
60 else csum
= cmat
[i
][j
];
62 for (k
= 0; k
< j
; k
++)
63 if (mode
== RE
) rsum
-= rmat
[i
][k
] * rmat
[k
][j
];
64 else csum
= csub(csum
, cmul(cmat
[i
][k
], cmat
[k
][j
]));
66 if (mode
== RE
) rmat
[i
][j
] = rsum
;
67 else cmat
[i
][j
] = csum
;
69 temp
= vv
[i
] * ((mode
== RE
) ? fabs(rsum
) : modulus(csum
));
70 if (temp
>= big
) { big
= temp
; imax
= i
; }
73 /* interchange rows if needed */
75 for (k
= 0; k
< n
; k
++) {
78 rmat
[imax
][k
] = rmat
[j
][k
];
83 cmat
[imax
][k
] = cmat
[j
][k
];
92 /* divide by the pivot element */
93 temp
= (mode
== RE
) ? fabs(rmat
[j
][j
]) : modulus(cmat
[j
][j
]);
94 if (temp
== 0.0) singular
= TRUE
;
97 rdum
= 1.0 / rmat
[j
][j
];
98 for (i
= j
+ 1; i
< n
; i
++) rmat
[i
][j
] *= rdum
;
101 cdum
= cdiv(cart2complex(1.0, 0.0), cmat
[j
][j
]);
102 for (i
= j
+ 1; i
< n
; i
++) cmat
[i
][j
] = cmul(cmat
[i
][j
], cdum
);
110 crlubksb(a
, n
, indx
, b
, mode
)
116 int i
, ii
, ip
, j
, singular
= FALSE
;
117 CMatrix ca
= (CMatrix
) a
;
118 CVector cb
= (CVector
) b
;
119 RMatrix ra
= (RMatrix
) a
;
120 RVector rb
= (RVector
) b
;
124 /* forward substitute using L part */
125 for (i
= 0, ii
= -1; i
< n
; i
++) {
136 for (j
= ii
; j
<= i
- 1; j
++)
137 if (mode
== RE
) rsum
-= ra
[i
][j
] * rb
[j
];
138 else csum
= csub(csum
, cmul(ca
[i
][j
], cb
[j
]));
141 if (rsum
!= 0.0) ii
= i
;
143 else if (csum
.real
!= 0.0 || csum
.imag
!= 0.0) ii
= i
;
145 if (mode
== RE
) rb
[i
] = rsum
;
149 /* back substitute using the U part */
150 for (i
= n
- 1; i
>= 0; i
--) {
153 for (j
= i
+ 1; j
< n
; j
++) rsum
-= ra
[i
][j
] * rb
[j
];
154 if (ra
[i
][i
] == 0.0) {
158 else rb
[i
] = rsum
/ ra
[i
][i
];
162 for (j
= i
+ 1; j
< n
; j
++) csum
= csub(csum
, cmul(ca
[i
][j
], cb
[j
]));
163 if (modulus(ca
[i
][i
]) == 0.0) {
167 else cb
[i
] = cdiv(csum
, ca
[i
][i
]);