001: /* AUTO-GENERATED */
002: package JSci.maths.matrices;
003:
004: import JSci.GlobalSettings;
005: import JSci.maths.ArrayMath;
006: import JSci.maths.ExtraMath;
007: import JSci.maths.LinearMath;
008: import JSci.maths.MaximumIterationsExceededException;
009: import JSci.maths.Mapping;
010: import JSci.maths.vectors.AbstractIntegerVector;
011: import JSci.maths.groups.AbelianGroup;
012: import JSci.maths.algebras.*;
013: import JSci.maths.fields.*;
014:
015: /**
016: * The AbstractIntegerSquareMatrix class provides an object for encapsulating integer square matrices.
017: * @version 2.2
018: * @author Mark Hale
019: */
020: public abstract class AbstractIntegerSquareMatrix extends
021: AbstractIntegerMatrix implements SquareMatrix {
022: protected transient AbstractDoubleSquareMatrix[] LU;
023: protected transient int[] LUpivot;
024:
025: /**
026: * Constructs a matrix.
027: */
028: protected AbstractIntegerSquareMatrix(final int size) {
029: super (size, size);
030: }
031:
032: /**
033: * Converts this matrix to a double matrix.
034: * @return a double square matrix
035: */
036: public AbstractDoubleMatrix toDoubleMatrix() {
037: final double ans[][] = new double[numRows][numCols];
038: for (int i = 0; i < numRows; i++) {
039: for (int j = 0; j < numCols; j++)
040: ans[i][j] = getElement(i, j);
041: }
042: return new DoubleSquareMatrix(ans);
043: }
044:
045: /**
046: * Converts this matrix to a complex matrix.
047: * @return a complex square matrix
048: */
049: public AbstractComplexMatrix toComplexMatrix() {
050: ComplexSquareMatrix cm = new ComplexSquareMatrix(numRows);
051: for (int i = 0; i < numRows; i++) {
052: for (int j = 0; j < numCols; j++)
053: cm.setElement(i, j, getElement(i, j), 0.0);
054: }
055: return cm;
056: }
057:
058: /**
059: * Returns true if this matrix is symmetric.
060: */
061: public boolean isSymmetric() {
062: return this .equals(this .transpose());
063: }
064:
065: /**
066: * Returns true if this matrix is unitary.
067: */
068: public boolean isUnitary() {
069: return this .multiply(this .transpose()).equals(
070: DoubleDiagonalMatrix.identity(numRows));
071: }
072:
073: /**
074: * Returns the determinant.
075: */
076: public int det() {
077: if (numRows == 2) {
078: return getElement(0, 0) * getElement(1, 1)
079: - getElement(0, 1) * getElement(1, 0);
080: } else {
081: final AbstractDoubleSquareMatrix lu[] = this
082: .luDecompose(null);
083: double det = lu[1].getElement(0, 0);
084: for (int i = 1; i < numRows; i++)
085: det *= lu[1].getElement(i, i);
086: return Math.round((float) det) * LUpivot[numRows];
087: }
088: }
089:
090: /**
091: * Returns the trace.
092: */
093: public int trace() {
094: int result = getElement(0, 0);
095: for (int i = 1; i < numRows; i++)
096: result += getElement(i, i);
097: return result;
098: }
099:
100: //============
101: // OPERATIONS
102: //============
103:
104: /**
105: * Returns the negative of this matrix.
106: */
107: public AbelianGroup.Member negate() {
108: final int array[][] = new int[numRows][numCols];
109: for (int i = 0; i < numRows; i++) {
110: array[i][0] = -getElement(i, 0);
111: for (int j = 1; j < numCols; j++)
112: array[i][j] = -getElement(i, j);
113: }
114: return new IntegerSquareMatrix(array);
115: }
116:
117: // ADDITION
118:
119: /**
120: * Returns the addition of this matrix and another.
121: * @param m a int square matrix
122: * @exception MatrixDimensionException If the matrices are not square or different sizes.
123: */
124: public final AbstractIntegerMatrix add(final AbstractIntegerMatrix m) {
125: if (m instanceof AbstractIntegerSquareMatrix)
126: return add((AbstractIntegerSquareMatrix) m);
127: else if (numRows == m.rows() && numCols == m.columns())
128: return add(new SquareMatrixAdaptor(m));
129: else
130: throw new MatrixDimensionException(
131: "Matrices are different sizes.");
132: }
133:
134: /**
135: * Returns the addition of this matrix and another.
136: * @param m a int square matrix
137: * @exception MatrixDimensionException If the matrices are different sizes.
138: */
139: public AbstractIntegerSquareMatrix add(
140: final AbstractIntegerSquareMatrix m) {
141: if (numRows == m.rows() && numCols == m.columns()) {
142: final int array[][] = new int[numRows][numCols];
143: for (int i = 0; i < numRows; i++) {
144: array[i][0] = getElement(i, 0) + m.getElement(i, 0);
145: for (int j = 1; j < numCols; j++)
146: array[i][j] = getElement(i, j) + m.getElement(i, j);
147: }
148: return new IntegerSquareMatrix(array);
149: } else {
150: throw new MatrixDimensionException(
151: "Matrices are different sizes.");
152: }
153: }
154:
155: // SUBTRACTION
156:
157: /**
158: * Returns the subtraction of this matrix and another.
159: * @param m a int square matrix
160: * @exception MatrixDimensionException If the matrices are not square or different sizes.
161: */
162: public final AbstractIntegerMatrix subtract(
163: final AbstractIntegerMatrix m) {
164: if (m instanceof AbstractIntegerSquareMatrix)
165: return subtract((AbstractIntegerSquareMatrix) m);
166: else if (numRows == m.rows() && numCols == m.columns())
167: return subtract(new SquareMatrixAdaptor(m));
168: else
169: throw new MatrixDimensionException(
170: "Matrices are different sizes.");
171: }
172:
173: /**
174: * Returns the subtraction of this matrix by another.
175: * @param m a int square matrix
176: * @exception MatrixDimensionException If the matrices are different sizes.
177: */
178: public AbstractIntegerSquareMatrix subtract(
179: final AbstractIntegerSquareMatrix m) {
180: if (numRows == m.rows() && numCols == m.columns()) {
181: final int array[][] = new int[numRows][numCols];
182: for (int i = 0; i < numRows; i++) {
183: array[i][0] = getElement(i, 0) - m.getElement(i, 0);
184: for (int j = 1; j < numCols; j++)
185: array[i][j] = getElement(i, j) - m.getElement(i, j);
186: }
187: return new IntegerSquareMatrix(array);
188: } else {
189: throw new MatrixDimensionException(
190: "Matrices are different sizes.");
191: }
192: }
193:
194: // SCALAR MULTIPLICATION
195:
196: /**
197: * Returns the multiplication of this matrix by a scalar.
198: * @param x a int.
199: * @return a int square matrix.
200: */
201: public AbstractIntegerMatrix scalarMultiply(final int x) {
202: final int array[][] = new int[numRows][numCols];
203: for (int i = 0; i < numRows; i++) {
204: array[i][0] = x * getElement(i, 0);
205: for (int j = 1; j < numCols; j++)
206: array[i][j] = x * getElement(i, j);
207: }
208: return new IntegerSquareMatrix(array);
209: }
210:
211: // SCALAR DIVISON
212:
213: // SCALAR PRODUCT
214:
215: /**
216: * Returns the scalar product of this matrix and another.
217: * @param m a int square matrix.
218: * @exception MatrixDimensionException If the matrices are not square or different sizes.
219: */
220: public final int scalarProduct(final AbstractIntegerMatrix m) {
221: if (m instanceof AbstractIntegerSquareMatrix)
222: return scalarProduct((AbstractIntegerSquareMatrix) m);
223: else if (numRows == m.rows() && numCols == m.columns())
224: return scalarProduct(new SquareMatrixAdaptor(m));
225: else
226: throw new MatrixDimensionException(
227: "Matrices are different sizes.");
228: }
229:
230: /**
231: * Returns the scalar product of this matrix and another.
232: * @param m a int square matrix.
233: * @exception MatrixDimensionException If the matrices are different sizes.
234: */
235: public int scalarProduct(final AbstractIntegerSquareMatrix m) {
236: if (numRows == m.rows() && numCols == m.columns()) {
237: int ans = 0;
238: for (int i = 0; i < numRows; i++) {
239: ans += getElement(i, 0) * m.getElement(i, 0);
240: for (int j = 1; j < numCols; j++)
241: ans += getElement(i, j) * m.getElement(i, j);
242: }
243: return ans;
244: } else {
245: throw new MatrixDimensionException(
246: "Matrices are different sizes.");
247: }
248: }
249:
250: // MATRIX MULTIPLICATION
251:
252: /**
253: * Returns the multiplication of this matrix and another.
254: * @param m a int square matrix
255: * @exception MatrixDimensionException If the matrices are different sizes.
256: */
257: public AbstractIntegerSquareMatrix multiply(
258: final AbstractIntegerSquareMatrix m) {
259: if (numCols == m.rows()) {
260: final int mColumns = m.columns();
261: final int array[][] = new int[numRows][mColumns];
262: for (int j = 0; j < numRows; j++) {
263: for (int k = 0; k < mColumns; k++) {
264: array[j][k] = getElement(j, 0) * m.getElement(0, k);
265: for (int n = 1; n < numCols; n++)
266: array[j][k] += getElement(j, n)
267: * m.getElement(n, k);
268: }
269: }
270: return new IntegerSquareMatrix(array);
271: } else {
272: throw new MatrixDimensionException("Incompatible matrices.");
273: }
274: }
275:
276: // DIRECT SUM
277:
278: /**
279: * Returns the direct sum of this matrix and another.
280: */
281: public AbstractIntegerSquareMatrix directSum(
282: final AbstractIntegerSquareMatrix m) {
283: final int array[][] = new int[numRows + m.numRows][numCols
284: + m.numCols];
285: for (int i = 0; i < numRows; i++) {
286: for (int j = 0; j < numCols; j++)
287: array[i][j] = getElement(i, j);
288: }
289: for (int i = 0; i < m.numRows; i++) {
290: for (int j = 0; j < m.numCols; j++)
291: array[i + numRows][j + numCols] = m.getElement(i, j);
292: }
293: return new IntegerSquareMatrix(array);
294: }
295:
296: // TENSOR PRODUCT
297:
298: /**
299: * Returns the tensor product of this matrix and another.
300: */
301: public AbstractIntegerSquareMatrix tensor(
302: final AbstractIntegerSquareMatrix m) {
303: final int array[][] = new int[numRows * m.numRows][numCols
304: * m.numCols];
305: for (int i = 0; i < numRows; i++) {
306: for (int j = 0; j < numCols; j++) {
307: for (int k = 0; k < m.numRows; j++) {
308: for (int l = 0; l < m.numCols; l++)
309: array[i * m.numRows + k][j * m.numCols + l] = getElement(
310: i, j)
311: * m.getElement(k, l);
312: }
313: }
314: }
315: return new IntegerSquareMatrix(array);
316: }
317:
318: // TRANSPOSE
319:
320: /**
321: * Returns the transpose of this matrix.
322: * @return a int square matrix
323: */
324: public Matrix transpose() {
325: final int array[][] = new int[numCols][numRows];
326: for (int i = 0; i < numRows; i++) {
327: array[0][i] = getElement(i, 0);
328: for (int j = 1; j < numCols; j++)
329: array[j][i] = getElement(i, j);
330: }
331: return new IntegerSquareMatrix(array);
332: }
333:
334: // INVERSE
335:
336: /**
337: * Returns the inverse of this matrix.
338: * @return a double square matrix.
339: */
340: public AbstractDoubleSquareMatrix inverse() {
341: final int N = numRows;
342: final double arrayL[][] = new double[N][N];
343: final double arrayU[][] = new double[N][N];
344: final AbstractDoubleSquareMatrix lu[] = this .luDecompose(null);
345: arrayL[0][0] = 1.0 / lu[0].getElement(0, 0);
346: arrayU[0][0] = 1.0 / lu[1].getElement(0, 0);
347: for (int i = 1; i < N; i++) {
348: arrayL[i][i] = 1.0 / lu[0].getElement(i, i);
349: arrayU[i][i] = 1.0 / lu[1].getElement(i, i);
350: }
351: for (int i = 0; i < N - 1; i++) {
352: for (int j = i + 1; j < N; j++) {
353: double tmpL = 0.0, tmpU = 0.0;
354: for (int k = i; k < j; k++) {
355: tmpL -= lu[0].getElement(j, k) * arrayL[k][i];
356: tmpU -= arrayU[i][k] * lu[1].getElement(k, j);
357: }
358: arrayL[j][i] = tmpL / lu[0].getElement(j, j);
359: arrayU[i][j] = tmpU / lu[1].getElement(j, j);
360: }
361: }
362: // matrix multiply arrayU x arrayL
363: final double inv[][] = new double[N][N];
364: for (int i = 0; i < N; i++) {
365: for (int j = 0; j < i; j++) {
366: for (int k = i; k < N; k++)
367: inv[i][LUpivot[j]] += arrayU[i][k] * arrayL[k][j];
368: }
369: for (int j = i; j < N; j++) {
370: for (int k = j; k < N; k++)
371: inv[i][LUpivot[j]] += arrayU[i][k] * arrayL[k][j];
372: }
373: }
374: return new DoubleSquareMatrix(inv);
375: }
376:
377: // LU DECOMPOSITION
378:
379: /**
380: * Returns the LU decomposition of this matrix.
381: * @param pivot an empty array of length <code>rows()+1</code>
382: * to hold the pivot information (null if not interested).
383: * The last array element will contain the parity.
384: * @return an array with [0] containing the L-matrix
385: * and [1] containing the U-matrix.
386: * @jsci.planetmath LUDecomposition
387: */
388: public AbstractDoubleSquareMatrix[] luDecompose(int pivot[]) {
389: if (LU != null) {
390: if (pivot != null)
391: System.arraycopy(LUpivot, 0, pivot, 0, pivot.length);
392: return LU;
393: }
394: int pivotrow;
395: final int N = numRows;
396: final double arrayL[][] = new double[N][N];
397: final double arrayU[][] = new double[N][N];
398: if (pivot == null)
399: pivot = new int[N + 1];
400: for (int i = 0; i < N; i++)
401: pivot[i] = i;
402: pivot[N] = 1;
403: // LU decomposition to arrayU
404: for (int j = 0; j < N; j++) {
405: for (int i = 0; i < j; i++) {
406: double tmp = getElement(pivot[i], j);
407: for (int k = 0; k < i; k++)
408: tmp -= arrayU[i][k] * arrayU[k][j];
409: arrayU[i][j] = tmp;
410: }
411: double max = 0.0;
412: pivotrow = j;
413: for (int i = j; i < N; i++) {
414: double tmp = getElement(pivot[i], j);
415: for (int k = 0; k < j; k++)
416: tmp -= arrayU[i][k] * arrayU[k][j];
417: arrayU[i][j] = tmp;
418: // while we're here search for a pivot for arrayU[j][j]
419: tmp = Math.abs(tmp);
420: if (tmp > max) {
421: max = tmp;
422: pivotrow = i;
423: }
424: }
425: // swap row j with pivotrow
426: if (pivotrow != j) {
427: double[] tmprow = arrayU[j];
428: arrayU[j] = arrayU[pivotrow];
429: arrayU[pivotrow] = tmprow;
430: int k = pivot[j];
431: pivot[j] = pivot[pivotrow];
432: pivot[pivotrow] = k;
433: // update parity
434: pivot[N] = -pivot[N];
435: }
436: // divide by pivot
437: double tmp = arrayU[j][j];
438: for (int i = j + 1; i < N; i++)
439: arrayU[i][j] /= tmp;
440: }
441: // move lower triangular part to arrayL
442: for (int j = 0; j < N; j++) {
443: arrayL[j][j] = 1.0;
444: for (int i = j + 1; i < N; i++) {
445: arrayL[i][j] = arrayU[i][j];
446: arrayU[i][j] = 0.0;
447: }
448: }
449: LU = new AbstractDoubleSquareMatrix[2];
450: LU[0] = new DoubleSquareMatrix(arrayL);
451: LU[1] = new DoubleSquareMatrix(arrayU);
452: LUpivot = new int[pivot.length];
453: System.arraycopy(pivot, 0, LUpivot, 0, pivot.length);
454: return LU;
455: }
456:
457: /**
458: * Returns the LU decomposition of this matrix.
459: * Warning: no pivoting.
460: * @return an array with [0] containing the L-matrix
461: * and [1] containing the U-matrix.
462: * @jsci.planetmath LUDecomposition
463: */
464: public AbstractDoubleSquareMatrix[] luDecompose() {
465: final int N = numRows;
466: final double arrayL[][] = new double[N][N];
467: final double arrayU[][] = new double[N][N];
468: // LU decomposition to arrayU
469: for (int j = 0; j < N; j++) {
470: for (int i = 0; i < j; i++) {
471: double tmp = getElement(i, j);
472: for (int k = 0; k < i; k++)
473: tmp -= arrayU[i][k] * arrayU[k][j];
474: arrayU[i][j] = tmp;
475: }
476: for (int i = j; i < N; i++) {
477: double tmp = getElement(i, j);
478: for (int k = 0; k < j; k++)
479: tmp -= arrayU[i][k] * arrayU[k][j];
480: arrayU[i][j] = tmp;
481: }
482: // divide
483: double tmp = arrayU[j][j];
484: for (int i = j + 1; i < N; i++)
485: arrayU[i][j] /= tmp;
486: }
487: // move lower triangular part to arrayL
488: for (int j = 0; j < N; j++) {
489: arrayL[j][j] = 1.0;
490: for (int i = j + 1; i < N; i++) {
491: arrayL[i][j] = arrayU[i][j];
492: arrayU[i][j] = 0.0;
493: }
494: }
495: AbstractDoubleSquareMatrix[] lu = new AbstractDoubleSquareMatrix[2];
496: lu[0] = new DoubleSquareMatrix(arrayL);
497: lu[1] = new DoubleSquareMatrix(arrayU);
498: return lu;
499: }
500:
501: // CHOLESKY DECOMPOSITION
502:
503: /**
504: * Returns the Cholesky decomposition of this matrix.
505: * Matrix must be symmetric and positive definite.
506: * @return an array with [0] containing the L-matrix and [1] containing the U-matrix.
507: * @jsci.planetmath CholeskyDecomposition
508: */
509: public AbstractDoubleSquareMatrix[] choleskyDecompose() {
510: final int N = numRows;
511: final double arrayL[][] = new double[N][N];
512: final double arrayU[][] = new double[N][N];
513: double tmp = Math.sqrt(getElement(0, 0));
514: arrayL[0][0] = arrayU[0][0] = tmp;
515: for (int i = 1; i < N; i++)
516: arrayL[i][0] = arrayU[0][i] = getElement(i, 0) / tmp;
517: for (int j = 1; j < N; j++) {
518: tmp = getElement(j, j);
519: for (int i = 0; i < j; i++)
520: tmp -= arrayL[j][i] * arrayL[j][i];
521: arrayL[j][j] = arrayU[j][j] = Math.sqrt(tmp);
522: for (int i = j + 1; i < N; i++) {
523: tmp = getElement(i, j);
524: for (int k = 0; k < i; k++)
525: tmp -= arrayL[j][k] * arrayU[k][i];
526: arrayL[i][j] = arrayU[j][i] = tmp / arrayU[j][j];
527: }
528: }
529: final AbstractDoubleSquareMatrix lu[] = new AbstractDoubleSquareMatrix[2];
530: lu[0] = new DoubleSquareMatrix(arrayL);
531: lu[1] = new DoubleSquareMatrix(arrayU);
532: return lu;
533: }
534:
535: // QR DECOMPOSITION
536:
537: /**
538: * Returns the QR decomposition of this matrix.
539: * Based on the code from <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> (public domain).
540: * @return an array with [0] containing the Q-matrix and [1] containing the R-matrix.
541: * @jsci.planetmath QRDecomposition
542: */
543: public AbstractDoubleSquareMatrix[] qrDecompose() {
544: final int N = numRows;
545: final double array[][] = new double[N][N];
546: final double arrayQ[][] = new double[N][N];
547: final double arrayR[][] = new double[N][N];
548: // copy matrix
549: for (int i = 0; i < N; i++) {
550: array[i][0] = getElement(i, 0);
551: for (int j = 1; j < N; j++)
552: array[i][j] = getElement(i, j);
553: }
554: for (int k = 0; k < N; k++) {
555: // compute l2-norm of kth column
556: double norm = array[k][k];
557: for (int i = k + 1; i < N; i++)
558: norm = ExtraMath.hypot(norm, array[i][k]);
559: if (norm != 0.0) {
560: // form kth Householder vector
561: if (array[k][k] < 0.0)
562: norm = -norm;
563: for (int i = k; i < N; i++)
564: array[i][k] /= norm;
565: array[k][k] += 1.0;
566: // apply transformation to remaining columns
567: for (int j = k + 1; j < N; j++) {
568: double s = array[k][k] * array[k][j];
569: for (int i = k + 1; i < N; i++)
570: s += array[i][k] * array[i][j];
571: s /= array[k][k];
572: for (int i = k; i < N; i++)
573: array[i][j] -= s * array[i][k];
574: }
575: }
576: arrayR[k][k] = -norm;
577: }
578: for (int k = N - 1; k >= 0; k--) {
579: arrayQ[k][k] = 1.0;
580: for (int j = k; j < N; j++) {
581: if (array[k][k] != 0.0) {
582: double s = array[k][k] * arrayQ[k][j];
583: for (int i = k + 1; i < N; i++)
584: s += array[i][k] * arrayQ[i][j];
585: s /= array[k][k];
586: for (int i = k; i < N; i++)
587: arrayQ[i][j] -= s * array[i][k];
588: }
589: }
590: }
591: for (int i = 0; i < N; i++) {
592: for (int j = i + 1; j < N; j++)
593: arrayR[i][j] = array[i][j];
594: }
595: final AbstractDoubleSquareMatrix qr[] = new AbstractDoubleSquareMatrix[2];
596: qr[0] = new DoubleSquareMatrix(arrayQ);
597: qr[1] = new DoubleSquareMatrix(arrayR);
598: return qr;
599: }
600:
601: // SINGULAR VALUE DECOMPOSITION
602:
603: /**
604: * Returns the singular value decomposition of this matrix.
605: * Based on the code from <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> (public domain).
606: * @return an array with [0] containing the U-matrix, [1] containing the S-matrix and [2] containing the V-matrix.
607: * @jsci.planetmath SingularValueDecomposition
608: */
609: public AbstractDoubleSquareMatrix[] singularValueDecompose() {
610: final int N = numRows;
611: final int Nm1 = N - 1;
612: final double array[][] = new double[N][N];
613: final double arrayU[][] = new double[N][N];
614: final double arrayS[] = new double[N];
615: final double arrayV[][] = new double[N][N];
616: final double e[] = new double[N];
617: final double work[] = new double[N];
618: // copy matrix
619: for (int i = 0; i < N; i++) {
620: array[i][0] = getElement(i, 0);
621: for (int j = 1; j < N; j++)
622: array[i][j] = getElement(i, j);
623: }
624: // reduce matrix to bidiagonal form
625: for (int k = 0; k < Nm1; k++) {
626: // compute the transformation for the kth column
627: // compute l2-norm of kth column
628: arrayS[k] = array[k][k];
629: for (int i = k + 1; i < N; i++)
630: arrayS[k] = ExtraMath.hypot(arrayS[k], array[i][k]);
631: if (arrayS[k] != 0.0) {
632: if (array[k][k] < 0.0)
633: arrayS[k] = -arrayS[k];
634: for (int i = k; i < N; i++)
635: array[i][k] /= arrayS[k];
636: array[k][k] += 1.0;
637: }
638: arrayS[k] = -arrayS[k];
639: for (int j = k + 1; j < N; j++) {
640: if (arrayS[k] != 0.0) {
641: // apply the transformation
642: double t = array[k][k] * array[k][j];
643: for (int i = k + 1; i < N; i++)
644: t += array[i][k] * array[i][j];
645: t /= array[k][k];
646: for (int i = k; i < N; i++)
647: array[i][j] -= t * array[i][k];
648: }
649: e[j] = array[k][j];
650: }
651: for (int i = k; i < N; i++)
652: arrayU[i][k] = array[i][k];
653: if (k < N - 2) {
654: // compute the kth row transformation
655: // compute l2-norm of kth column
656: e[k] = e[k + 1];
657: for (int i = k + 2; i < N; i++)
658: e[k] = ExtraMath.hypot(e[k], e[i]);
659: if (e[k] != 0.0) {
660: if (e[k + 1] < 0.0)
661: e[k] = -e[k];
662: for (int i = k + 1; i < N; i++)
663: e[i] /= e[k];
664: e[k + 1] += 1.0;
665: }
666: e[k] = -e[k];
667: if (e[k] != 0.0) {
668: // apply the transformation
669: for (int i = k + 1; i < N; i++) {
670: work[i] = 0.0;
671: for (int j = k + 1; j < N; j++)
672: work[i] += e[j] * array[i][j];
673: }
674: for (int j = k + 1; j < N; j++) {
675: double t = e[j] / e[k + 1];
676: for (int i = k + 1; i < N; i++)
677: array[i][j] -= t * work[i];
678: }
679: }
680: for (int i = k + 1; i < N; i++)
681: arrayV[i][k] = e[i];
682: }
683: }
684: // setup the final bidiagonal matrix of order p
685: int p = N;
686: arrayS[Nm1] = array[Nm1][Nm1];
687: e[N - 2] = array[N - 2][Nm1];
688: e[Nm1] = 0.0;
689: arrayU[Nm1][Nm1] = 1.0;
690: for (int k = N - 2; k >= 0; k--) {
691: if (arrayS[k] != 0.0) {
692: for (int j = k + 1; j < N; j++) {
693: double t = arrayU[k][k] * arrayU[k][j];
694: for (int i = k + 1; i < N; i++)
695: t += arrayU[i][k] * arrayU[i][j];
696: t /= arrayU[k][k];
697: for (int i = k; i < N; i++)
698: arrayU[i][j] -= t * arrayU[i][k];
699: }
700: for (int i = k; i < N; i++)
701: arrayU[i][k] = -arrayU[i][k];
702: arrayU[k][k] += 1.0;
703: for (int i = 0; i < k - 1; i++)
704: arrayU[i][k] = 0.0;
705: } else {
706: for (int i = 0; i < N; i++)
707: arrayU[i][k] = 0.0;
708: arrayU[k][k] = 1.0;
709: }
710: }
711: for (int k = Nm1; k >= 0; k--) {
712: if (k < N - 2 && e[k] != 0.0) {
713: for (int j = k + 1; j < N; j++) {
714: double t = arrayV[k + 1][k] * arrayV[k + 1][j];
715: for (int i = k + 2; i < N; i++)
716: t += arrayV[i][k] * arrayV[i][j];
717: t /= arrayV[k + 1][k];
718: for (int i = k + 1; i < N; i++)
719: arrayV[i][j] -= t * arrayV[i][k];
720: }
721: }
722: for (int i = 0; i < N; i++)
723: arrayV[i][k] = 0.0;
724: arrayV[k][k] = 1.0;
725: }
726: final double eps = Math.pow(2.0, -52.0);
727: int iter = 0;
728: while (p > 0) {
729: int k, action;
730: // action = 1 if arrayS[p] and e[k-1] are negligible and k<p
731: // action = 2 if arrayS[k] is negligible and k<p
732: // action = 3 if e[k-1] is negligible, k<p, and arrayS[k], ..., arrayS[p] are not negligible (QR step)
733: // action = 4 if e[p-1] is negligible (convergence)
734: for (k = p - 2; k >= -1; k--) {
735: if (k == -1)
736: break;
737: if (Math.abs(e[k]) <= eps
738: * (Math.abs(arrayS[k]) + Math
739: .abs(arrayS[k + 1]))) {
740: e[k] = 0.0;
741: break;
742: }
743: }
744: if (k == p - 2) {
745: action = 4;
746: } else {
747: int ks;
748: for (ks = p - 1; ks >= k; ks--) {
749: if (ks == k)
750: break;
751: double t = (ks != p ? Math.abs(e[ks]) : 0.0)
752: + (ks != k + 1 ? Math.abs(e[ks - 1]) : 0.0);
753: if (Math.abs(arrayS[ks]) <= eps * t) {
754: arrayS[ks] = 0.0;
755: break;
756: }
757: }
758: if (ks == k) {
759: action = 3;
760: } else if (ks == p - 1) {
761: action = 1;
762: } else {
763: action = 2;
764: k = ks;
765: }
766: }
767: k++;
768: switch (action) {
769: // deflate negligible arrayS[p]
770: case 1: {
771: double f = e[p - 2];
772: e[p - 2] = 0.0;
773: for (int j = p - 2; j >= k; j--) {
774: double t = ExtraMath.hypot(arrayS[j], f);
775: final double cs = arrayS[j] / t;
776: final double sn = f / t;
777: arrayS[j] = t;
778: if (j != k) {
779: f = -sn * e[j - 1];
780: e[j - 1] *= cs;
781: }
782: for (int i = 0; i < N; i++) {
783: t = cs * arrayV[i][j] + sn * arrayV[i][p - 1];
784: arrayV[i][p - 1] = -sn * arrayV[i][j] + cs
785: * arrayV[i][p - 1];
786: arrayV[i][j] = t;
787: }
788: }
789: }
790: break;
791: // split at negligible arrayS[k]
792: case 2: {
793: double f = e[k - 1];
794: e[k - 1] = 0.0;
795: for (int j = k; j < p; j++) {
796: double t = ExtraMath.hypot(arrayS[j], f);
797: final double cs = arrayS[j] / t;
798: final double sn = f / t;
799: arrayS[j] = t;
800: f = -sn * e[j];
801: e[j] *= cs;
802: for (int i = 0; i < N; i++) {
803: t = cs * arrayU[i][j] + sn * arrayU[i][k - 1];
804: arrayU[i][k - 1] = -sn * arrayU[i][j] + cs
805: * arrayU[i][k - 1];
806: arrayU[i][j] = t;
807: }
808: }
809: }
810: break;
811: // perform one QR step
812: case 3: {
813: // calculate the shift
814: final double scale = Math.max(Math.max(Math.max(Math
815: .max(Math.abs(arrayS[p - 1]), Math
816: .abs(arrayS[p - 2])), Math
817: .abs(e[p - 2])), Math.abs(arrayS[k])), Math
818: .abs(e[k]));
819: double sp = arrayS[p - 1] / scale;
820: double spm1 = arrayS[p - 2] / scale;
821: double epm1 = e[p - 2] / scale;
822: double sk = arrayS[k] / scale;
823: double ek = e[k] / scale;
824: double b = ((spm1 + sp) * (spm1 - sp) + epm1 * epm1) / 2.0;
825: double c = (sp * epm1) * (sp * epm1);
826: double shift = 0.0;
827: if (b != 0.0 || c != 0.0) {
828: shift = Math.sqrt(b * b + c);
829: if (b < 0.0)
830: shift = -shift;
831: shift = c / (b + shift);
832: }
833: double f = (sk + sp) * (sk - sp) + shift;
834: double g = sk * ek;
835: // chase zeros
836: for (int j = k; j < p - 1; j++) {
837: double t = ExtraMath.hypot(f, g);
838: double cs = f / t;
839: double sn = g / t;
840: if (j != k)
841: e[j - 1] = t;
842: f = cs * arrayS[j] + sn * e[j];
843: e[j] = cs * e[j] - sn * arrayS[j];
844: g = sn * arrayS[j + 1];
845: arrayS[j + 1] *= cs;
846: for (int i = 0; i < N; i++) {
847: t = cs * arrayV[i][j] + sn * arrayV[i][j + 1];
848: arrayV[i][j + 1] = -sn * arrayV[i][j] + cs
849: * arrayV[i][j + 1];
850: arrayV[i][j] = t;
851: }
852: t = ExtraMath.hypot(f, g);
853: cs = f / t;
854: sn = g / t;
855: arrayS[j] = t;
856: f = cs * e[j] + sn * arrayS[j + 1];
857: arrayS[j + 1] = -sn * e[j] + cs * arrayS[j + 1];
858: g = sn * e[j + 1];
859: e[j + 1] *= cs;
860: if (j < Nm1) {
861: for (int i = 0; i < N; i++) {
862: t = cs * arrayU[i][j] + sn
863: * arrayU[i][j + 1];
864: arrayU[i][j + 1] = -sn * arrayU[i][j] + cs
865: * arrayU[i][j + 1];
866: arrayU[i][j] = t;
867: }
868: }
869: }
870: e[p - 2] = f;
871: iter++;
872: }
873: break;
874: // convergence
875: case 4: {
876: // make the singular values positive
877: if (arrayS[k] <= 0.0) {
878: arrayS[k] = -arrayS[k];
879: for (int i = 0; i < p; i++)
880: arrayV[i][k] = -arrayV[i][k];
881: }
882: // order the singular values
883: while (k < p - 1) {
884: if (arrayS[k] >= arrayS[k + 1])
885: break;
886: double tmp = arrayS[k];
887: arrayS[k] = arrayS[k + 1];
888: arrayS[k + 1] = tmp;
889: if (k < Nm1) {
890: for (int i = 0; i < N; i++) {
891: tmp = arrayU[i][k + 1];
892: arrayU[i][k + 1] = arrayU[i][k];
893: arrayU[i][k] = tmp;
894: tmp = arrayV[i][k + 1];
895: arrayV[i][k + 1] = arrayV[i][k];
896: arrayV[i][k] = tmp;
897: }
898: }
899: k++;
900: }
901: iter = 0;
902: p--;
903: }
904: break;
905: }
906: }
907: final AbstractDoubleSquareMatrix svd[] = new AbstractDoubleSquareMatrix[3];
908: svd[0] = new DoubleSquareMatrix(arrayU);
909: svd[1] = new DoubleDiagonalMatrix(arrayS);
910: svd[2] = new DoubleSquareMatrix(arrayV);
911: return svd;
912: }
913:
914: private static class SquareMatrixAdaptor extends
915: AbstractIntegerSquareMatrix {
916: private final AbstractIntegerMatrix matrix;
917:
918: private SquareMatrixAdaptor(AbstractIntegerMatrix m) {
919: super (m.rows());
920: matrix = m;
921: }
922:
923: public int getElement(int i, int j) {
924: return matrix.getElement(i, j);
925: }
926:
927: public void setElement(int i, int j, int x) {
928: matrix.setElement(i, j, x);
929: }
930: }
931: }
|