0001: /* AUTO-GENERATED */
0002: package JSci.maths.matrices;
0003:
0004: import JSci.maths.ArrayMath;
0005: import JSci.maths.ExtraMath;
0006: import JSci.maths.LinearMath;
0007: import JSci.maths.Mapping;
0008: import JSci.maths.DimensionException;
0009: import JSci.maths.MaximumIterationsExceededException;
0010: import JSci.maths.vectors.AbstractIntegerVector;
0011: import JSci.maths.vectors.IntegerVector;
0012: import JSci.maths.groups.AbelianGroup;
0013: import JSci.maths.algebras.*;
0014: import JSci.maths.fields.*;
0015:
0016: /**
0017: * The IntegerSquareMatrix class provides an object for encapsulating integer square matrices.
0018: * @version 2.2
0019: * @author Mark Hale
0020: */
0021: public class IntegerSquareMatrix extends AbstractIntegerSquareMatrix {
0022: /**
0023: * Array containing the elements of the matrix.
0024: */
0025: protected final int matrix[][];
0026:
0027: /**
0028: * Constructs a matrix by wrapping an array.
0029: * @param array an assigned value
0030: */
0031: public IntegerSquareMatrix(final int array[][]) {
0032: super (array.length);
0033: if (!ArrayMath.isSquare(array))
0034: throw new MatrixDimensionException("Array is not square.");
0035: matrix = array;
0036: }
0037:
0038: /**
0039: * Constructs an empty matrix.
0040: */
0041: public IntegerSquareMatrix(final int size) {
0042: this (new int[size][size]);
0043: }
0044:
0045: /**
0046: * Constructs a matrix from an array of vectors (columns).
0047: * @param array an assigned value
0048: */
0049: public IntegerSquareMatrix(final AbstractIntegerVector array[]) {
0050: this (array.length);
0051: for (int i = 0; i < numRows; i++) {
0052: for (int j = 0; j < numCols; j++)
0053: matrix[i][j] = array[j].getComponent(i);
0054: }
0055: }
0056:
0057: /**
0058: * Compares two ${nativeTyp} matrices for equality.
0059: * @param m a int matrix
0060: */
0061: public boolean equals(AbstractIntegerMatrix m, double tol) {
0062: if (m != null && numRows == m.rows() && numCols == m.columns()) {
0063: int sumSqr = 0;
0064: for (int i = 0; i < numRows; i++) {
0065: for (int j = 0; j < numCols; j++) {
0066: int delta = matrix[i][j] - m.getElement(i, j);
0067: sumSqr += delta * delta;
0068: }
0069: }
0070: return (sumSqr <= tol * tol);
0071: } else {
0072: return false;
0073: }
0074: }
0075:
0076: /**
0077: * Returns a string representing this matrix.
0078: */
0079: public String toString() {
0080: final StringBuffer buf = new StringBuffer(5 * numRows * numCols);
0081: for (int i = 0; i < numRows; i++) {
0082: for (int j = 0; j < numCols; j++) {
0083: buf.append(matrix[i][j]);
0084: buf.append(' ');
0085: }
0086: buf.append('\n');
0087: }
0088: return buf.toString();
0089: }
0090:
0091: /**
0092: * Converts this matrix to a double matrix.
0093: * @return a double matrix
0094: */
0095: public AbstractDoubleMatrix toDoubleMatrix() {
0096: final double ans[][] = new double[numRows][numCols];
0097: for (int i = 0; i < numRows; i++) {
0098: for (int j = 0; j < numCols; j++)
0099: ans[i][j] = matrix[i][j];
0100: }
0101: return new DoubleSquareMatrix(ans);
0102: }
0103:
0104: /**
0105: * Converts this matrix to a complex matrix.
0106: * @return a complex matrix
0107: */
0108: public AbstractComplexMatrix toComplexMatrix() {
0109: ComplexSquareMatrix cm = new ComplexSquareMatrix(numRows);
0110: for (int i = 0; i < numRows; i++) {
0111: for (int j = 0; j < numCols; j++)
0112: cm.setElement(i, j, matrix[i][j], 0.0);
0113: }
0114: return cm;
0115: }
0116:
0117: /**
0118: * Returns an element of the matrix.
0119: * @param i row index of the element
0120: * @param j column index of the element
0121: * @exception MatrixDimensionException If attempting to access an invalid element.
0122: */
0123: public int getElement(int i, int j) {
0124: if (i >= 0 && i < numRows && j >= 0 && j < numCols)
0125: return matrix[i][j];
0126: else
0127: throw new MatrixDimensionException(getInvalidElementMsg(i,
0128: j));
0129: }
0130:
0131: /**
0132: * Sets the value of an element of the matrix.
0133: * Should only be used to initialise this matrix.
0134: * @param i row index of the element
0135: * @param j column index of the element
0136: * @param x a number
0137: * @exception MatrixDimensionException If attempting to access an invalid element.
0138: */
0139: public void setElement(int i, int j, int x) {
0140: if (i >= 0 && i < numRows && j >= 0 && j < numCols)
0141: matrix[i][j] = x;
0142: else
0143: throw new MatrixDimensionException(getInvalidElementMsg(i,
0144: j));
0145: }
0146:
0147: /**
0148: * Returns the l<sup><img border=0 alt="infinity" src="doc-files/infinity.gif"></sup>-norm.
0149: * @author Taber Smith
0150: */
0151: public int infNorm() {
0152: int result = 0;
0153: for (int i = 0; i < numRows; i++) {
0154: int tmpResult = 0;
0155: for (int j = 0; j < numCols; j++)
0156: tmpResult += Math.abs(matrix[i][j]);
0157: if (tmpResult > result)
0158: result = tmpResult;
0159: }
0160: return result;
0161: }
0162:
0163: /**
0164: * Returns the Frobenius or Hilbert-Schmidt (l<sup>2</sup>) norm.
0165: * @jsci.planetmath FrobeniusMatrixNorm
0166: */
0167: public double frobeniusNorm() {
0168: double result = 0.0;
0169: for (int j, i = 0; i < numRows; i++) {
0170: for (j = 0; j < numCols; j++)
0171: result = ExtraMath.hypot(result, matrix[i][j]);
0172: }
0173: return result;
0174: }
0175:
0176: /**
0177: * Returns the determinant.
0178: */
0179: public int det() {
0180: if (numRows == 2) {
0181: return matrix[0][0] * matrix[1][1] - matrix[0][1]
0182: * matrix[1][0];
0183: } else {
0184: final DoubleSquareMatrix lu[] = (DoubleSquareMatrix[]) this
0185: .luDecompose(null);
0186: double det = lu[1].matrix[0][0];
0187: for (int i = 1; i < numRows; i++)
0188: det *= lu[1].matrix[i][i];
0189: return Math.round((float) det) * LUpivot[numRows];
0190: }
0191: }
0192:
0193: /**
0194: * Returns the trace.
0195: */
0196: public int trace() {
0197: int result = matrix[0][0];
0198: for (int i = 1; i < numRows; i++)
0199: result += matrix[i][i];
0200: return result;
0201: }
0202:
0203: //============
0204: // OPERATIONS
0205: //============
0206:
0207: /**
0208: * Returns the negative of this matrix.
0209: */
0210: public AbelianGroup.Member negate() {
0211: final int array[][] = new int[numRows][numCols];
0212: for (int i = 0; i < numRows; i++) {
0213: array[i][0] = -matrix[i][0];
0214: for (int j = 1; j < numCols; j++)
0215: array[i][j] = -matrix[i][j];
0216: }
0217: return new IntegerSquareMatrix(array);
0218: }
0219:
0220: // ADDITION
0221:
0222: /**
0223: * Returns the addition of this matrix and another.
0224: * @param m a int matrix
0225: * @exception MatrixDimensionException If the matrices are different sizes.
0226: */
0227: public AbstractIntegerSquareMatrix add(
0228: final AbstractIntegerSquareMatrix m) {
0229: if (m instanceof IntegerSquareMatrix)
0230: return add((IntegerSquareMatrix) m);
0231:
0232: if (numRows == m.rows() && numCols == m.columns()) {
0233: final int array[][] = new int[numRows][numCols];
0234: for (int i = 0; i < numRows; i++) {
0235: array[i][0] = matrix[i][0] + m.getElement(i, 0);
0236: for (int j = 1; j < numCols; j++)
0237: array[i][j] = matrix[i][j] + m.getElement(i, j);
0238: }
0239: return new IntegerSquareMatrix(array);
0240: } else {
0241: throw new MatrixDimensionException(
0242: "Matrices are different sizes.");
0243: }
0244: }
0245:
0246: public IntegerSquareMatrix add(final IntegerSquareMatrix m) {
0247: if (numRows == m.numRows && numCols == m.numCols) {
0248: final int array[][] = new int[numRows][numCols];
0249: for (int j, i = 0; i < numRows; i++) {
0250: array[i][0] = matrix[i][0] + m.matrix[i][0];
0251: for (j = 1; j < numCols; j++)
0252: array[i][j] = matrix[i][j] + m.matrix[i][j];
0253: }
0254: return new IntegerSquareMatrix(array);
0255: } else
0256: throw new MatrixDimensionException(
0257: "Matrices are different sizes.");
0258: }
0259:
0260: // SUBTRACTION
0261:
0262: /**
0263: * Returns the subtraction of this matrix by another.
0264: * @param m a int matrix
0265: * @exception MatrixDimensionException If the matrices are different sizes.
0266: */
0267: public AbstractIntegerSquareMatrix subtract(
0268: final AbstractIntegerSquareMatrix m) {
0269: if (m instanceof IntegerSquareMatrix)
0270: return subtract((IntegerSquareMatrix) m);
0271:
0272: if (numRows == m.rows() && numCols == m.columns()) {
0273: final int array[][] = new int[numRows][numCols];
0274: for (int i = 0; i < numRows; i++) {
0275: array[i][0] = matrix[i][0] - m.getElement(i, 0);
0276: for (int j = 1; j < numCols; j++)
0277: array[i][j] = matrix[i][j] - m.getElement(i, j);
0278: }
0279: return new IntegerSquareMatrix(array);
0280: } else {
0281: throw new MatrixDimensionException(
0282: "Matrices are different sizes.");
0283: }
0284: }
0285:
0286: public IntegerSquareMatrix subtract(final IntegerSquareMatrix m) {
0287: if (numRows == m.numRows && numCols == m.numCols) {
0288: final int array[][] = new int[numRows][numCols];
0289: for (int i = 0; i < numRows; i++) {
0290: array[i][0] = matrix[i][0] - m.matrix[i][0];
0291: for (int j = 1; j < numCols; j++)
0292: array[i][j] = matrix[i][j] - m.matrix[i][j];
0293: }
0294: return new IntegerSquareMatrix(array);
0295: } else
0296: throw new MatrixDimensionException(
0297: "Matrices are different sizes.");
0298: }
0299:
0300: // SCALAR MULTIPLICATION
0301:
0302: /**
0303: * Returns the multiplication of this matrix by a scalar.
0304: * @param x a int.
0305: * @return a int square matrix.
0306: */
0307: public AbstractIntegerMatrix scalarMultiply(final int x) {
0308: final int array[][] = new int[numRows][numCols];
0309: for (int i = 0; i < numRows; i++) {
0310: array[i][0] = x * matrix[i][0];
0311: for (int j = 1; j < numCols; j++)
0312: array[i][j] = x * matrix[i][j];
0313: }
0314: return new IntegerSquareMatrix(array);
0315: }
0316:
0317: // SCALAR DIVISON
0318:
0319: // SCALAR PRODUCT
0320:
0321: /**
0322: * Returns the scalar product of this matrix and another.
0323: * @param m a int matrix.
0324: * @exception MatrixDimensionException If the matrices are different sizes.
0325: */
0326: public int scalarProduct(final AbstractIntegerSquareMatrix m) {
0327: if (m instanceof IntegerSquareMatrix)
0328: return scalarProduct((IntegerSquareMatrix) m);
0329:
0330: if (numRows == m.rows() && numCols == m.columns()) {
0331: int ans = 0;
0332: for (int i = 0; i < numRows; i++) {
0333: ans += matrix[i][0] * m.getElement(i, 0);
0334: for (int j = 1; j < numCols; j++)
0335: ans += matrix[i][j] * m.getElement(i, j);
0336: }
0337: return ans;
0338: } else {
0339: throw new MatrixDimensionException(
0340: "Matrices are different sizes.");
0341: }
0342: }
0343:
0344: public int scalarProduct(final IntegerSquareMatrix m) {
0345: if (numRows == m.numRows && numCols == m.numCols) {
0346: int ans = 0;
0347: for (int i = 0; i < numRows; i++) {
0348: ans += matrix[i][0] * m.matrix[i][0];
0349: for (int j = 1; j < numCols; j++)
0350: ans += matrix[i][j] * m.matrix[i][j];
0351: }
0352: return ans;
0353: } else
0354: throw new MatrixDimensionException(
0355: "Matrices are different sizes.");
0356: }
0357:
0358: // MATRIX MULTIPLICATION
0359:
0360: /**
0361: * Returns the multiplication of a vector by this matrix.
0362: * @param v a int vector.
0363: * @exception DimensionException If the matrix and vector are incompatible.
0364: */
0365: public AbstractIntegerVector multiply(final AbstractIntegerVector v) {
0366: if (numCols == v.dimension()) {
0367: final int array[] = new int[numRows];
0368: for (int i = 0; i < numRows; i++) {
0369: array[i] = matrix[i][0] * v.getComponent(0);
0370: for (int j = 1; j < numCols; j++)
0371: array[i] += matrix[i][j] * v.getComponent(j);
0372: }
0373: return new IntegerVector(array);
0374: } else {
0375: throw new DimensionException(
0376: "Matrix and vector are incompatible.");
0377: }
0378: }
0379:
0380: /**
0381: * Returns the multiplication of this matrix and another.
0382: * @param m a int matrix
0383: * @return a AbstractIntegerMatrix or a AbstractIntegerSquareMatrix as appropriate
0384: * @exception MatrixDimensionException If the matrices are incompatible.
0385: */
0386: public AbstractIntegerSquareMatrix multiply(
0387: final AbstractIntegerSquareMatrix m) {
0388: if (m instanceof IntegerSquareMatrix)
0389: return multiply((IntegerSquareMatrix) m);
0390:
0391: if (numCols == m.rows()) {
0392: final int mColumns = m.columns();
0393: final int array[][] = new int[numRows][mColumns];
0394: for (int j = 0; j < numRows; j++) {
0395: for (int k = 0; k < mColumns; k++) {
0396: array[j][k] = matrix[j][0] * m.getElement(0, k);
0397: for (int n = 1; n < numCols; n++)
0398: array[j][k] += matrix[j][n]
0399: * m.getElement(n, k);
0400: }
0401: }
0402: return new IntegerSquareMatrix(array);
0403: } else {
0404: throw new MatrixDimensionException("Incompatible matrices.");
0405: }
0406: }
0407:
0408: public IntegerSquareMatrix multiply(final IntegerSquareMatrix m) {
0409: if (numCols == m.numRows) {
0410: final int array[][] = new int[numRows][m.numCols];
0411: for (int j = 0; j < numRows; j++) {
0412: for (int k = 0; k < m.numCols; k++) {
0413: array[j][k] = matrix[j][0] * m.matrix[0][k];
0414: for (int n = 1; n < numCols; n++)
0415: array[j][k] += matrix[j][n] * m.matrix[n][k];
0416: }
0417: }
0418: return new IntegerSquareMatrix(array);
0419: } else
0420: throw new MatrixDimensionException("Incompatible matrices.");
0421: }
0422:
0423: // DIRECT SUM
0424:
0425: /**
0426: * Returns the direct sum of this matrix and another.
0427: */
0428: public AbstractIntegerSquareMatrix directSum(
0429: final AbstractIntegerSquareMatrix m) {
0430: final int array[][] = new int[numRows + m.numRows][numCols
0431: + m.numCols];
0432: for (int i = 0; i < numRows; i++) {
0433: for (int j = 0; j < numCols; j++)
0434: array[i][j] = matrix[i][j];
0435: }
0436: for (int i = 0; i < m.numRows; i++) {
0437: for (int j = 0; j < m.numCols; j++)
0438: array[i + numRows][j + numCols] = m.getElement(i, j);
0439: }
0440: return new IntegerSquareMatrix(array);
0441: }
0442:
0443: // TENSOR PRODUCT
0444:
0445: /**
0446: * Returns the tensor product of this matrix and another.
0447: */
0448: public AbstractIntegerSquareMatrix tensor(
0449: final AbstractIntegerSquareMatrix m) {
0450: final int array[][] = new int[numRows * m.numRows][numCols
0451: * m.numCols];
0452: for (int i = 0; i < numRows; i++) {
0453: for (int j = 0; j < numCols; j++) {
0454: for (int k = 0; k < m.numRows; j++) {
0455: for (int l = 0; l < m.numCols; l++)
0456: array[i * m.numRows + k][j * m.numCols + l] = matrix[i][j]
0457: * m.getElement(k, l);
0458: }
0459: }
0460: }
0461: return new IntegerSquareMatrix(array);
0462: }
0463:
0464: // TRANSPOSE
0465:
0466: /**
0467: * Returns the transpose of this matrix.
0468: * @return a int matrix
0469: */
0470: public Matrix transpose() {
0471: final int array[][] = new int[numCols][numRows];
0472: for (int i = 0; i < numRows; i++) {
0473: array[0][i] = matrix[i][0];
0474: for (int j = 1; j < numCols; j++)
0475: array[j][i] = matrix[i][j];
0476: }
0477: return new IntegerSquareMatrix(array);
0478: }
0479:
0480: // INVERSE
0481:
0482: /**
0483: * Returns the inverse of this matrix.
0484: * @return a double square matrix.
0485: */
0486: public AbstractDoubleSquareMatrix inverse() {
0487: final int N = numRows;
0488: final double arrayL[][] = new double[N][N];
0489: final double arrayU[][] = new double[N][N];
0490: final DoubleSquareMatrix lu[] = (DoubleSquareMatrix[]) this
0491: .luDecompose(null);
0492: arrayL[0][0] = 1.0 / lu[0].matrix[0][0];
0493: arrayU[0][0] = 1.0 / lu[1].matrix[0][0];
0494: for (int i = 1; i < N; i++) {
0495: arrayL[i][i] = 1.0 / lu[0].matrix[i][i];
0496: arrayU[i][i] = 1.0 / lu[1].matrix[i][i];
0497: }
0498: for (int i = 0; i < N - 1; i++) {
0499: for (int j = i + 1; j < N; j++) {
0500: double tmpL = 0.0, tmpU = 0.0;
0501: for (int k = i; k < j; k++) {
0502: tmpL -= lu[0].matrix[j][k] * arrayL[k][i];
0503: tmpU -= arrayU[i][k] * lu[1].matrix[k][j];
0504: }
0505: arrayL[j][i] = tmpL / lu[0].matrix[j][j];
0506: arrayU[i][j] = tmpU / lu[1].matrix[j][j];
0507: }
0508: }
0509: // matrix multiply arrayU x arrayL
0510: final double inv[][] = new double[N][N];
0511: for (int i = 0; i < N; i++) {
0512: for (int j = 0; j < i; j++) {
0513: for (int k = i; k < N; k++)
0514: inv[i][LUpivot[j]] += arrayU[i][k] * arrayL[k][j];
0515: }
0516: for (int j = i; j < N; j++) {
0517: for (int k = j; k < N; k++)
0518: inv[i][LUpivot[j]] += arrayU[i][k] * arrayL[k][j];
0519: }
0520: }
0521: return new DoubleSquareMatrix(inv);
0522: }
0523:
0524: // LU DECOMPOSITION
0525:
0526: /**
0527: * Returns the LU decomposition of this matrix.
0528: * @param pivot an empty array of length <code>rows()+1</code>
0529: * to hold the pivot information (null if not interested).
0530: * The last array element will contain the parity.
0531: * @return an array with [0] containing the L-matrix
0532: * and [1] containing the U-matrix.
0533: * @jsci.planetmath LUDecomposition
0534: */
0535: public final AbstractDoubleSquareMatrix[] luDecompose(int pivot[]) {
0536: if (LU != null) {
0537: if (pivot != null)
0538: System.arraycopy(LUpivot, 0, pivot, 0, pivot.length);
0539: return LU;
0540: }
0541: int pivotrow;
0542: final int N = numRows;
0543: final double arrayL[][] = new double[N][N];
0544: final double arrayU[][] = new double[N][N];
0545: double tmp;
0546: if (pivot == null)
0547: pivot = new int[N + 1];
0548: for (int i = 0; i < N; i++)
0549: pivot[i] = i;
0550: pivot[N] = 1;
0551: // LU decomposition to arrayU
0552: for (int j = 0; j < N; j++) {
0553: for (int i = 0; i < j; i++) {
0554: tmp = matrix[pivot[i]][j];
0555: for (int k = 0; k < i; k++)
0556: tmp -= arrayU[i][k] * arrayU[k][j];
0557: arrayU[i][j] = tmp;
0558: }
0559: double max = 0.0;
0560: pivotrow = j;
0561: for (int i = j; i < N; i++) {
0562: tmp = matrix[pivot[i]][j];
0563: for (int k = 0; k < j; k++)
0564: tmp -= arrayU[i][k] * arrayU[k][j];
0565: arrayU[i][j] = tmp;
0566: // while we're here search for a pivot for arrayU[j][j]
0567: tmp = Math.abs(tmp);
0568: if (tmp > max) {
0569: max = tmp;
0570: pivotrow = i;
0571: }
0572: }
0573: // swap row j with pivotrow
0574: if (pivotrow != j) {
0575: double[] tmprow = arrayU[j];
0576: arrayU[j] = arrayU[pivotrow];
0577: arrayU[pivotrow] = tmprow;
0578: int k = pivot[j];
0579: pivot[j] = pivot[pivotrow];
0580: pivot[pivotrow] = k;
0581: // update parity
0582: pivot[N] = -pivot[N];
0583: }
0584: // divide by pivot
0585: tmp = arrayU[j][j];
0586: for (int i = j + 1; i < N; i++)
0587: arrayU[i][j] /= tmp;
0588: }
0589: // move lower triangular part to arrayL
0590: for (int j = 0; j < N; j++) {
0591: arrayL[j][j] = 1.0;
0592: for (int i = j + 1; i < N; i++) {
0593: arrayL[i][j] = arrayU[i][j];
0594: arrayU[i][j] = 0.0;
0595: }
0596: }
0597: LU = new DoubleSquareMatrix[2];
0598: LU[0] = new DoubleSquareMatrix(arrayL);
0599: LU[1] = new DoubleSquareMatrix(arrayU);
0600: LUpivot = new int[pivot.length];
0601: System.arraycopy(pivot, 0, LUpivot, 0, pivot.length);
0602: return LU;
0603: }
0604:
0605: /**
0606: * Returns the LU decomposition of this matrix.
0607: * Warning: no pivoting.
0608: * @return an array with [0] containing the L-matrix
0609: * and [1] containing the U-matrix.
0610: * @jsci.planetmath LUDecomposition
0611: */
0612: public final AbstractDoubleSquareMatrix[] luDecompose() {
0613: final int N = numRows;
0614: final double arrayL[][] = new double[N][N];
0615: final double arrayU[][] = new double[N][N];
0616: double tmp;
0617: // LU decomposition to arrayU
0618: for (int j = 0; j < N; j++) {
0619: for (int i = 0; i < j; i++) {
0620: tmp = matrix[i][j];
0621: for (int k = 0; k < i; k++)
0622: tmp -= arrayU[i][k] * arrayU[k][j];
0623: arrayU[i][j] = tmp;
0624: }
0625: for (int i = j; i < N; i++) {
0626: tmp = matrix[i][j];
0627: for (int k = 0; k < j; k++)
0628: tmp -= arrayU[i][k] * arrayU[k][j];
0629: arrayU[i][j] = tmp;
0630: }
0631: // divide
0632: tmp = arrayU[j][j];
0633: for (int i = j + 1; i < N; i++)
0634: arrayU[i][j] /= tmp;
0635: }
0636: // move lower triangular part to arrayL
0637: for (int j = 0; j < N; j++) {
0638: arrayL[j][j] = 1.0;
0639: for (int i = j + 1; i < N; i++) {
0640: arrayL[i][j] = arrayU[i][j];
0641: arrayU[i][j] = 0.0;
0642: }
0643: }
0644: DoubleSquareMatrix[] lu = new DoubleSquareMatrix[2];
0645: lu[0] = new DoubleSquareMatrix(arrayL);
0646: lu[1] = new DoubleSquareMatrix(arrayU);
0647: return lu;
0648: }
0649:
0650: // CHOLESKY DECOMPOSITION
0651:
0652: /**
0653: * Returns the Cholesky decomposition of this matrix.
0654: * Matrix must be symmetric and positive definite.
0655: * @return an array with [0] containing the L-matrix and [1] containing the U-matrix.
0656: * @jsci.planetmath CholeskyDecomposition
0657: */
0658: public AbstractDoubleSquareMatrix[] choleskyDecompose() {
0659: final int N = numRows;
0660: final double arrayL[][] = new double[N][N];
0661: final double arrayU[][] = new double[N][N];
0662: double tmp = Math.sqrt(matrix[0][0]);
0663: arrayL[0][0] = arrayU[0][0] = tmp;
0664: for (int i = 1; i < N; i++)
0665: arrayL[i][0] = arrayU[0][i] = matrix[i][0] / tmp;
0666: for (int j = 1; j < N; j++) {
0667: tmp = matrix[j][j];
0668: for (int i = 0; i < j; i++)
0669: tmp -= arrayL[j][i] * arrayL[j][i];
0670: arrayL[j][j] = arrayU[j][j] = Math.sqrt(tmp);
0671: for (int i = j + 1; i < N; i++) {
0672: tmp = matrix[i][j];
0673: for (int k = 0; k < i; k++)
0674: tmp -= arrayL[j][k] * arrayU[k][i];
0675: arrayL[i][j] = arrayU[j][i] = tmp / arrayU[j][j];
0676: }
0677: }
0678: final AbstractDoubleSquareMatrix lu[] = new AbstractDoubleSquareMatrix[2];
0679: lu[0] = new DoubleSquareMatrix(arrayL);
0680: lu[1] = new DoubleSquareMatrix(arrayU);
0681: return lu;
0682: }
0683:
0684: // QR DECOMPOSITION
0685:
0686: /**
0687: * Returns the QR decomposition of this matrix.
0688: * Based on the code from <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> (public domain).
0689: * @return an array with [0] containing the Q-matrix and [1] containing the R-matrix.
0690: * @jsci.planetmath QRDecomposition
0691: */
0692: public AbstractDoubleSquareMatrix[] qrDecompose() {
0693: final int N = numRows;
0694: final double array[][] = new double[N][N];
0695: final double arrayQ[][] = new double[N][N];
0696: final double arrayR[][] = new double[N][N];
0697: // copy matrix
0698: for (int i = 0; i < N; i++) {
0699: array[i][0] = matrix[i][0];
0700: for (int j = 1; j < N; j++)
0701: array[i][j] = matrix[i][j];
0702: }
0703: for (int k = 0; k < N; k++) {
0704: // compute l2-norm of kth column
0705: double norm = array[k][k];
0706: for (int i = k + 1; i < N; i++)
0707: norm = ExtraMath.hypot(norm, array[i][k]);
0708: if (norm != 0.0) {
0709: // form kth Householder vector
0710: if (array[k][k] < 0.0)
0711: norm = -norm;
0712: for (int i = k; i < N; i++)
0713: array[i][k] /= norm;
0714: array[k][k] += 1.0;
0715: // apply transformation to remaining columns
0716: for (int j = k + 1; j < N; j++) {
0717: double s = array[k][k] * array[k][j];
0718: for (int i = k + 1; i < N; i++)
0719: s += array[i][k] * array[i][j];
0720: s /= array[k][k];
0721: for (int i = k; i < N; i++)
0722: array[i][j] -= s * array[i][k];
0723: }
0724: }
0725: arrayR[k][k] = -norm;
0726: }
0727: for (int k = N - 1; k >= 0; k--) {
0728: arrayQ[k][k] = 1.0;
0729: for (int j = k; j < N; j++) {
0730: if (array[k][k] != 0.0) {
0731: double s = array[k][k] * arrayQ[k][j];
0732: for (int i = k + 1; i < N; i++)
0733: s += array[i][k] * arrayQ[i][j];
0734: s /= array[k][k];
0735: for (int i = k; i < N; i++)
0736: arrayQ[i][j] -= s * array[i][k];
0737: }
0738: }
0739: }
0740: for (int i = 0; i < N; i++) {
0741: for (int j = i + 1; j < N; j++)
0742: arrayR[i][j] = array[i][j];
0743: }
0744: final AbstractDoubleSquareMatrix qr[] = new AbstractDoubleSquareMatrix[2];
0745: qr[0] = new DoubleSquareMatrix(arrayQ);
0746: qr[1] = new DoubleSquareMatrix(arrayR);
0747: return qr;
0748: }
0749:
0750: // SINGULAR VALUE DECOMPOSITION
0751:
0752: /**
0753: * Returns the singular value decomposition of this matrix.
0754: * Based on the code from <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> (public domain).
0755: * @return an array with [0] containing the U-matrix, [1] containing the S-matrix and [2] containing the V-matrix.
0756: * @jsci.planetmath SingularValueDecomposition
0757: */
0758: public AbstractDoubleSquareMatrix[] singularValueDecompose() {
0759: final int N = numRows;
0760: final int Nm1 = N - 1;
0761: final double array[][] = new double[N][N];
0762: final double arrayU[][] = new double[N][N];
0763: final double arrayS[] = new double[N];
0764: final double arrayV[][] = new double[N][N];
0765: final double e[] = new double[N];
0766: final double work[] = new double[N];
0767: // copy matrix
0768: for (int i = 0; i < N; i++) {
0769: array[i][0] = matrix[i][0];
0770: for (int j = 1; j < N; j++)
0771: array[i][j] = matrix[i][j];
0772: }
0773: // reduce matrix to bidiagonal form
0774: for (int k = 0; k < Nm1; k++) {
0775: // compute the transformation for the kth column
0776: // compute l2-norm of kth column
0777: arrayS[k] = array[k][k];
0778: for (int i = k + 1; i < N; i++)
0779: arrayS[k] = ExtraMath.hypot(arrayS[k], array[i][k]);
0780: if (arrayS[k] != 0.0) {
0781: if (array[k][k] < 0.0)
0782: arrayS[k] = -arrayS[k];
0783: for (int i = k; i < N; i++)
0784: array[i][k] /= arrayS[k];
0785: array[k][k] += 1.0;
0786: }
0787: arrayS[k] = -arrayS[k];
0788: for (int j = k + 1; j < N; j++) {
0789: if (arrayS[k] != 0.0) {
0790: // apply the transformation
0791: double t = array[k][k] * array[k][j];
0792: for (int i = k + 1; i < N; i++)
0793: t += array[i][k] * array[i][j];
0794: t /= array[k][k];
0795: for (int i = k; i < N; i++)
0796: array[i][j] -= t * array[i][k];
0797: }
0798: e[j] = array[k][j];
0799: }
0800: for (int i = k; i < N; i++)
0801: arrayU[i][k] = array[i][k];
0802: if (k < N - 2) {
0803: // compute the kth row transformation
0804: // compute l2-norm of kth column
0805: e[k] = e[k + 1];
0806: for (int i = k + 2; i < N; i++)
0807: e[k] = ExtraMath.hypot(e[k], e[i]);
0808: if (e[k] != 0.0) {
0809: if (e[k + 1] < 0.0)
0810: e[k] = -e[k];
0811: for (int i = k + 1; i < N; i++)
0812: e[i] /= e[k];
0813: e[k + 1] += 1.0;
0814: }
0815: e[k] = -e[k];
0816: if (e[k] != 0.0) {
0817: // apply the transformation
0818: for (int i = k + 1; i < N; i++) {
0819: work[i] = 0.0;
0820: for (int j = k + 1; j < N; j++)
0821: work[i] += e[j] * array[i][j];
0822: }
0823: for (int j = k + 1; j < N; j++) {
0824: double t = e[j] / e[k + 1];
0825: for (int i = k + 1; i < N; i++)
0826: array[i][j] -= t * work[i];
0827: }
0828: }
0829: for (int i = k + 1; i < N; i++)
0830: arrayV[i][k] = e[i];
0831: }
0832: }
0833: // setup the final bidiagonal matrix of order p
0834: int p = N;
0835: arrayS[Nm1] = array[Nm1][Nm1];
0836: e[N - 2] = array[N - 2][Nm1];
0837: e[Nm1] = 0.0;
0838: arrayU[Nm1][Nm1] = 1.0;
0839: for (int k = N - 2; k >= 0; k--) {
0840: if (arrayS[k] != 0.0) {
0841: for (int j = k + 1; j < N; j++) {
0842: double t = arrayU[k][k] * arrayU[k][j];
0843: for (int i = k + 1; i < N; i++)
0844: t += arrayU[i][k] * arrayU[i][j];
0845: t /= arrayU[k][k];
0846: for (int i = k; i < N; i++)
0847: arrayU[i][j] -= t * arrayU[i][k];
0848: }
0849: for (int i = k; i < N; i++)
0850: arrayU[i][k] = -arrayU[i][k];
0851: arrayU[k][k] += 1.0;
0852: for (int i = 0; i < k - 1; i++)
0853: arrayU[i][k] = 0.0;
0854: } else {
0855: for (int i = 0; i < N; i++)
0856: arrayU[i][k] = 0.0;
0857: arrayU[k][k] = 1.0;
0858: }
0859: }
0860: for (int k = Nm1; k >= 0; k--) {
0861: if (k < N - 2 && e[k] != 0.0) {
0862: for (int j = k + 1; j < N; j++) {
0863: double t = arrayV[k + 1][k] * arrayV[k + 1][j];
0864: for (int i = k + 2; i < N; i++)
0865: t += arrayV[i][k] * arrayV[i][j];
0866: t /= arrayV[k + 1][k];
0867: for (int i = k + 1; i < N; i++)
0868: arrayV[i][j] -= t * arrayV[i][k];
0869: }
0870: }
0871: for (int i = 0; i < N; i++)
0872: arrayV[i][k] = 0.0;
0873: arrayV[k][k] = 1.0;
0874: }
0875: final double eps = Math.pow(2.0, -52.0);
0876: int iter = 0;
0877: while (p > 0) {
0878: int k, action;
0879: // action = 1 if arrayS[p] and e[k-1] are negligible and k<p
0880: // action = 2 if arrayS[k] is negligible and k<p
0881: // action = 3 if e[k-1] is negligible, k<p, and arrayS[k], ..., arrayS[p] are not negligible (QR step)
0882: // action = 4 if e[p-1] is negligible (convergence)
0883: for (k = p - 2; k >= -1; k--) {
0884: if (k == -1)
0885: break;
0886: if (Math.abs(e[k]) <= eps
0887: * (Math.abs(arrayS[k]) + Math
0888: .abs(arrayS[k + 1]))) {
0889: e[k] = 0.0;
0890: break;
0891: }
0892: }
0893: if (k == p - 2) {
0894: action = 4;
0895: } else {
0896: int ks;
0897: for (ks = p - 1; ks >= k; ks--) {
0898: if (ks == k)
0899: break;
0900: double t = (ks != p ? Math.abs(e[ks]) : 0.0)
0901: + (ks != k + 1 ? Math.abs(e[ks - 1]) : 0.0);
0902: if (Math.abs(arrayS[ks]) <= eps * t) {
0903: arrayS[ks] = 0.0;
0904: break;
0905: }
0906: }
0907: if (ks == k) {
0908: action = 3;
0909: } else if (ks == p - 1) {
0910: action = 1;
0911: } else {
0912: action = 2;
0913: k = ks;
0914: }
0915: }
0916: k++;
0917: switch (action) {
0918: // deflate negligible arrayS[p]
0919: case 1: {
0920: double f = e[p - 2];
0921: e[p - 2] = 0.0;
0922: for (int j = p - 2; j >= k; j--) {
0923: double t = ExtraMath.hypot(arrayS[j], f);
0924: final double cs = arrayS[j] / t;
0925: final double sn = f / t;
0926: arrayS[j] = t;
0927: if (j != k) {
0928: f = -sn * e[j - 1];
0929: e[j - 1] *= cs;
0930: }
0931: for (int i = 0; i < N; i++) {
0932: t = cs * arrayV[i][j] + sn * arrayV[i][p - 1];
0933: arrayV[i][p - 1] = -sn * arrayV[i][j] + cs
0934: * arrayV[i][p - 1];
0935: arrayV[i][j] = t;
0936: }
0937: }
0938: }
0939: break;
0940: // split at negligible arrayS[k]
0941: case 2: {
0942: double f = e[k - 1];
0943: e[k - 1] = 0.0;
0944: for (int j = k; j < p; j++) {
0945: double t = ExtraMath.hypot(arrayS[j], f);
0946: final double cs = arrayS[j] / t;
0947: final double sn = f / t;
0948: arrayS[j] = t;
0949: f = -sn * e[j];
0950: e[j] *= cs;
0951: for (int i = 0; i < N; i++) {
0952: t = cs * arrayU[i][j] + sn * arrayU[i][k - 1];
0953: arrayU[i][k - 1] = -sn * arrayU[i][j] + cs
0954: * arrayU[i][k - 1];
0955: arrayU[i][j] = t;
0956: }
0957: }
0958: }
0959: break;
0960: // perform one QR step
0961: case 3: {
0962: // calculate the shift
0963: final double scale = Math.max(Math.max(Math.max(Math
0964: .max(Math.abs(arrayS[p - 1]), Math
0965: .abs(arrayS[p - 2])), Math
0966: .abs(e[p - 2])), Math.abs(arrayS[k])), Math
0967: .abs(e[k]));
0968: double sp = arrayS[p - 1] / scale;
0969: double spm1 = arrayS[p - 2] / scale;
0970: double epm1 = e[p - 2] / scale;
0971: double sk = arrayS[k] / scale;
0972: double ek = e[k] / scale;
0973: double b = ((spm1 + sp) * (spm1 - sp) + epm1 * epm1) / 2.0;
0974: double c = (sp * epm1) * (sp * epm1);
0975: double shift = 0.0;
0976: if (b != 0.0 || c != 0.0) {
0977: shift = Math.sqrt(b * b + c);
0978: if (b < 0.0)
0979: shift = -shift;
0980: shift = c / (b + shift);
0981: }
0982: double f = (sk + sp) * (sk - sp) + shift;
0983: double g = sk * ek;
0984: // chase zeros
0985: for (int j = k; j < p - 1; j++) {
0986: double t = ExtraMath.hypot(f, g);
0987: double cs = f / t;
0988: double sn = g / t;
0989: if (j != k)
0990: e[j - 1] = t;
0991: f = cs * arrayS[j] + sn * e[j];
0992: e[j] = cs * e[j] - sn * arrayS[j];
0993: g = sn * arrayS[j + 1];
0994: arrayS[j + 1] *= cs;
0995: for (int i = 0; i < N; i++) {
0996: t = cs * arrayV[i][j] + sn * arrayV[i][j + 1];
0997: arrayV[i][j + 1] = -sn * arrayV[i][j] + cs
0998: * arrayV[i][j + 1];
0999: arrayV[i][j] = t;
1000: }
1001: t = ExtraMath.hypot(f, g);
1002: cs = f / t;
1003: sn = g / t;
1004: arrayS[j] = t;
1005: f = cs * e[j] + sn * arrayS[j + 1];
1006: arrayS[j + 1] = -sn * e[j] + cs * arrayS[j + 1];
1007: g = sn * e[j + 1];
1008: e[j + 1] *= cs;
1009: if (j < Nm1) {
1010: for (int i = 0; i < N; i++) {
1011: t = cs * arrayU[i][j] + sn
1012: * arrayU[i][j + 1];
1013: arrayU[i][j + 1] = -sn * arrayU[i][j] + cs
1014: * arrayU[i][j + 1];
1015: arrayU[i][j] = t;
1016: }
1017: }
1018: }
1019: e[p - 2] = f;
1020: iter++;
1021: }
1022: break;
1023: // convergence
1024: case 4: {
1025: // make the singular values positive
1026: if (arrayS[k] <= 0.0) {
1027: arrayS[k] = -arrayS[k];
1028: for (int i = 0; i < p; i++)
1029: arrayV[i][k] = -arrayV[i][k];
1030: }
1031: // order the singular values
1032: while (k < p - 1) {
1033: if (arrayS[k] >= arrayS[k + 1])
1034: break;
1035: double tmp = arrayS[k];
1036: arrayS[k] = arrayS[k + 1];
1037: arrayS[k + 1] = tmp;
1038: if (k < Nm1) {
1039: for (int i = 0; i < N; i++) {
1040: tmp = arrayU[i][k + 1];
1041: arrayU[i][k + 1] = arrayU[i][k];
1042: arrayU[i][k] = tmp;
1043: tmp = arrayV[i][k + 1];
1044: arrayV[i][k + 1] = arrayV[i][k];
1045: arrayV[i][k] = tmp;
1046: }
1047: }
1048: k++;
1049: }
1050: iter = 0;
1051: p--;
1052: }
1053: break;
1054: }
1055: }
1056: final AbstractDoubleSquareMatrix svd[] = new AbstractDoubleSquareMatrix[3];
1057: svd[0] = new DoubleSquareMatrix(arrayU);
1058: svd[1] = new DoubleDiagonalMatrix(arrayS);
1059: svd[2] = new DoubleSquareMatrix(arrayV);
1060: return svd;
1061: }
1062:
1063: }
|