001: // -- @(#)LUDecomp.java 1.4 99/03/22
002:
003: /**
004: * A special form of Matrix that holds a copy of the row interchange
005: * information generated during a LU decompose operation. It also holds
006: * the value of a flg that indicates the odd/even number of row interchanges.
007: */package uk.org.ponder.matrix;
008:
009: public class LUDecomp extends Matrix {
010:
011: public final static double TINY = Double.NaN;
012:
013: int indx[]; // Table of row interchanges.
014: double d; // +1.0, -1.0 (row interchanges even, odd)
015: double twod[][]; // 2D version of the matrix
016:
017: /**
018: * LU Decomposition of a Matrix a. Specalised case of a Matrix
019: * operator. This constructor performs the LU decomposition of
020: * the supplied matrix and creates an instance of its decomposition.
021: */
022:
023: public LUDecomp(Matrix m) throws SingularException,
024: NotSquareException {
025:
026: // Make a new matrix of the correct size but initalised
027: // with zeros. Call deepCopy to duplicate the values
028: // from the supplied matrix. These get overwritten by
029: // the decomposition routine.
030: super (m.rows, m.cols);
031:
032: if (m.rows != m.cols) {
033: throw new NotSquareException();
034: }
035: twod = m.asArray();
036: indx = new int[m.rows];
037:
038: LUdecomp();
039:
040: twod = null; // Dont need it any more.
041: }
042:
043: /**
044: * A relatively faithful port of the ludcmp routine from numerical
045: * recipies in 'C'. (See chapter 2). Principal changes are to make
046: * the arrays run over the range 0..n-1 rather than the botched
047: * version in the book. Changes have also been made to allow for
048: * the encapsulation of the mechanisim within an object framework.
049: */
050: private void LUdecomp() throws SingularException {
051:
052: // As in the book this routine overwrites the copy of the
053: // data array with the LU decomposition of the original
054: // data. The main difference is that this is a local copy
055: // of the data that was supplied with the constructor.
056:
057: int i, imax = -1, j, k;
058: int n = rows; // Matrix dimension
059: double big, dum, sum, temp;
060: double[] vv = new double[rows];
061: double[][] a = twod;
062:
063: d = 1.0; // No row interchanges yet.
064:
065: for (i = 0; i < n; i++) { // Loop over rows to get the implicit
066: big = 0.0; // scaling information.
067:
068: for (j = 0; j < n; j++) {
069: if ((temp = Math.abs(a[i][j])) > big)
070: big = temp;
071: }
072:
073: if (big == 0.0) {
074: throw new SingularException();
075: }
076:
077: // No nonzero largest element.
078:
079: // Save the scaling.
080: vv[i] = 1.0 / big;
081:
082: }
083:
084: for (j = 0; j < n; j++) { // Loop over col. in Crout's
085: for (i = 0; i < j; i++) { // method.
086: sum = a[i][j]; // eqn. 2.3.12 except for i=j
087: for (k = 0; k < i; k++) {
088: sum -= a[i][k] * a[k][j];
089: }
090: a[i][j] = sum;
091: }
092: big = 0.0; // Init search for largest
093: for (i = j; i < n; i++) { // pivot element.
094: sum = a[i][j];
095: for (k = 0; k < j; k++) {
096: sum -= a[i][k] * a[k][j];
097: }
098: a[i][j] = sum;
099: if ((dum = vv[i] * Math.abs(sum)) >= big) {
100: // Is the figure of merit for this
101: // pivot better than that found so far?
102: big = dum;
103: imax = i;
104: }
105: }
106: if (j != imax) { // Do we need to interchange rows?
107: for (k = 0; k < n; k++) {
108: dum = a[imax][k];
109: a[imax][k] = a[j][k];
110: a[j][k] = dum;
111: }
112: d = -d; // Change the parity of d.
113: vv[imax] = vv[j]; // Also interchange the scale.
114: }
115:
116: indx[j] = imax;
117:
118: if (a[j][j] == 0.0) {
119: // If the pivot element is zero the matrix is
120: // singular (at least to the precision of the
121: // algorithim) For some applications on
122: // singular matricies, it is desirable to
123: // substitute TINY for zero.
124: a[j][j] = TINY;
125: }
126:
127: if (j != n) {
128:
129: // Now, finally, divide by the pivot element.
130: dum = 1.0 / (a[j][j]);
131: for (i = j + 1; i < n; i++) {
132: a[i][j] *= dum;
133: }
134: }
135: }
136:
137: for (i = 0; i < rows; i++) {
138: for (j = 0; j < cols; j++) {
139: setMval(i, j, a[i][j]);
140: }
141: }
142:
143: }
144:
145: /**
146: * Perform back substitution based on this LU decomposed matrix. This
147: * effectively solves a set of N linear equations A.X = B.
148: * @param data The 1 x N data vector. (RHS vector B)
149: * @returns solution vector X.
150: */
151: public double[] backSubstitute(double[] data)
152: throws SizeMismatchException {
153:
154: int i, ii = -1, ip, j;
155: int n = rows;
156: double sum;
157:
158: if (n != data.length) {
159: throw new SizeMismatchException();
160: }
161:
162: double[] b = new double[n];
163:
164: for (i = 0; i < n; i++) {
165: b[i] = data[i];
166: }
167:
168: for (i = 0; i < n; i++) { // When ii is set to a positive value
169: ip = indx[i]; // it will become the index of the
170: sum = b[ip]; // first nonvanishing element of b.
171: b[ip] = b[i]; // We now do the forward substitution
172: // (eqn. 2.3.6) The only new wrinkle
173: if (ii != -1) { // is to unscramble the permutation
174: // as we go.
175: for (j = ii; j <= i - 1; j++)
176: sum -= getMval(i, j) * b[j];
177: } else if (sum != 0.0) {
178: // A nonzero element was encountered
179: // so from now on we have to do the
180: ii = i; // sums in the loop above.
181: }
182: b[i] = sum;
183: }
184:
185: for (i = n - 1; i >= 0; i--) { // Now we do the back substitution,
186: sum = b[i]; // (eqn. 2.3.7)
187: for (j = i + 1; j < n; j++)
188: sum -= getMval(i, j) * b[j];
189: b[i] = sum / getMval(i, i);// Store a component of the
190: // solution vector X.
191: }
192:
193: return b;
194: }
195:
196: /**
197: * Return the inverse of the original matrix used in the construction
198: * of this LUDecomp object. Simply involves back substitution of the
199: * identity matrix.
200: * @returns The inverse matrix of that used in the creation of this
201: * LUDecomp object.
202: */
203: public Matrix luinvert() throws SizeMismatchException {
204:
205: int i, j, n = rows;
206: Matrix ret = new Matrix(rows, cols); // nb. class col.
207: double[] tmp, col = new double[n];
208:
209: for (j = 0; j < n; j++) {
210: for (i = 0; i < n; i++)
211: col[i] = 0.0;
212: col[j] = 1.0;
213: tmp = backSubstitute(col);
214:
215: for (i = 0; i < n; i++)
216: ret.setMval(i, j, tmp[i]);
217: }
218: return ret;
219: }
220:
221: public double ludeterminant() {
222:
223: int n = rows;
224: double res = d; // d is a class variable
225:
226: for (int i = 0; i < n; i++)
227: res *= getMval(i, i);
228:
229: return res;
230: }
231:
232: /**
233: * Render the contents of this LU decomposition into a string
234: * primarily for debugging use.
235: */
236: public String toString() {
237:
238: String tmp = super .toString();
239:
240: tmp = tmp + "\n" + "[";
241: for (int i = 0; i < indx.length; i++) {
242: if (i != 0)
243: tmp = tmp + ", " + indx[i];
244: else
245: tmp = tmp + indx[i];
246: }
247: tmp = tmp + "] (" + d + ")\n";
248:
249: return (tmp);
250: }
251: }
|