fix packaging (into central file) and copyright dates.
[CommonLispStat.git] / lib / ludecomp.c
blobe659d1cb5b4dfb192365a96533bf84d1438cf926
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. */
8 #include "linalg.h"
10 int
11 crludcmp(Matrix mat, size_t n, IVector indx, int mode, double *d)
13 int i, imax, j, k;
14 int singular = FALSE;
15 double big, temp;
16 Complex cdum, csum;
17 double rdum, rsum;
18 CMatrix cmat = (CMatrix) mat;
19 RMatrix rmat = (RMatrix) mat;
20 RVector vv;
22 vv = rvector(n);
23 *d = 1.0;
25 /* set up the pivot permutation vector */
26 for (i = 0; i < n; i++) indx[i] = i;
28 /* get scaling information for implicit pivoting */
29 for (i = 0; i < n; i++) {
30 big = 0.0;
31 for (j = 0; j < n; j++) {
32 temp = (mode == RE) ? fabs(rmat[i][j]) : modulus(cmat[i][j]);
33 if (temp > big) big = temp;
35 if (big == 0.0) {
36 vv[i] = 1.0; /* no scaling for zero rows */
37 singular = TRUE;
39 else vv[i] = 1.0 / big;
42 /* loop over columns for Crout's method */
43 for (j = 0; j < n; j++) {
44 for (i = 0; i < j; i++) {
45 if (mode == RE) rsum = rmat[i][j];
46 else csum = cmat[i][j];
48 for (k = 0; k < i; k++)
50 if (mode == RE) rsum -= rmat[i][k] * rmat[k][j];
51 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;
57 big = 0.0;
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++)
64 if (mode == RE) rsum -= rmat[i][k] * rmat[k][j];
65 else csum = csub(csum, cmul(cmat[i][k], cmat[k][j]));
68 if (mode == RE) rmat[i][j] = rsum;
69 else cmat[i][j] = csum;
71 temp = vv[i] * ((mode == RE) ? fabs(rsum) : modulus(csum));
72 if (temp >= big) { big = temp; imax = i; }
75 /* interchange rows if needed */
76 if (j != imax) {
77 for (k = 0; k < n; k++) {
78 if (mode == RE) {
79 rdum = rmat[imax][k];
80 rmat[imax][k] = rmat[j][k];
81 rmat[j][k] = rdum;
83 else {
84 cdum = cmat[imax][k];
85 cmat[imax][k] = cmat[j][k];
86 cmat[j][k] = cdum;
89 *d = -(*d);
90 vv[imax] = vv[j];
92 indx[j] = imax;
94 /* divide by the pivot element */
95 temp = (mode == RE) ? fabs(rmat[j][j]) : modulus(cmat[j][j]);
96 if (temp == 0.0) singular = TRUE;
97 else if (j < n - 1) {
98 if (mode == RE) {
99 rdum = 1.0 / rmat[j][j];
100 for (i = j + 1; i < n; i++) rmat[i][j] *= rdum;
102 else {
103 cdum = cdiv(cart2complex(1.0, 0.0), cmat[j][j]);
104 for (i = j + 1; i < n; i++)
106 cmat[i][j] = cmul(cmat[i][j], cdum);
112 free_vector(vv);
113 return(singular);
116 int
117 crlubksb(double **a, int n, int *indx, double *b, int mode)
119 int i, ii, ip, j, singular = FALSE;
120 CMatrix ca = (CMatrix) a;
121 CVector cb = (CVector) b;
122 RMatrix ra = (RMatrix) a;
123 RVector rb = (RVector) b;
124 double rsum;
125 Complex csum;
127 /* forward substitute using L part */
128 for (i = 0, ii = -1; i < n; i++) {
129 ip = indx[i];
130 if (mode == RE) {
131 rsum = rb[ip];
132 rb[ip] = rb[i];
134 else {
135 csum = cb[ip];
136 cb[ip] = cb[i];
138 if (ii >= 0)
139 for (j = ii; j <= i - 1; j++)
140 if (mode == RE) rsum -= ra[i][j] * rb[j];
141 else csum = csub(csum, cmul(ca[i][j], cb[j]));
142 else {
143 if (mode == RE) {
144 if (rsum != 0.0) ii = i;
146 else if (csum.real != 0.0 || csum.imag != 0.0) ii = i;
148 if (mode == RE) rb[i] = rsum;
149 else cb[i] = csum;
152 /* back substitute using the U part */
153 for (i = n - 1; i >= 0; i--) {
154 if (mode == RE) {
155 rsum = rb[i];
156 for (j = i + 1; j < n; j++) rsum -= ra[i][j] * rb[j];
157 if (ra[i][i] == 0.0) {
158 singular = TRUE;
159 break;
161 else rb[i] = rsum / ra[i][i];
163 else {
164 csum = cb[i];
165 for (j = i + 1; j < n; j++) csum = csub(csum, cmul(ca[i][j], cb[j]));
166 if (modulus(ca[i][i]) == 0.0) {
167 singular = TRUE;
168 break;
170 else cb[i] = cdiv(csum, ca[i][i]);
174 return(singular);