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.AbstractDoubleVector;
0011: import JSci.maths.vectors.DoubleVector;
0012: import JSci.maths.groups.AbelianGroup;
0013: import JSci.maths.algebras.*;
0014: import JSci.maths.fields.*;
0015:
0016: /**
0017: * The DoubleSquareMatrix class provides an object for encapsulating double square matrices.
0018: * @version 2.2
0019: * @author Mark Hale
0020: */
0021: public class DoubleSquareMatrix extends AbstractDoubleSquareMatrix {
0022: /**
0023: * Array containing the elements of the matrix.
0024: */
0025: protected final double matrix[][];
0026:
0027: /**
0028: * Constructs a matrix by wrapping an array.
0029: * @param array an assigned value
0030: */
0031: public DoubleSquareMatrix(final double 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 DoubleSquareMatrix(final int size) {
0042: this (new double[size][size]);
0043: }
0044:
0045: /**
0046: * Constructs a matrix from an array of vectors (columns).
0047: * @param array an assigned value
0048: */
0049: public DoubleSquareMatrix(final AbstractDoubleVector 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 double matrix
0060: */
0061: public boolean equals(AbstractDoubleMatrix m, double tol) {
0062: if (m != null && numRows == m.rows() && numCols == m.columns()) {
0063: double sumSqr = 0;
0064: for (int i = 0; i < numRows; i++) {
0065: for (int j = 0; j < numCols; j++) {
0066: double 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 an integer matrix.
0093: * @return an integer matrix
0094: */
0095: public AbstractIntegerMatrix toIntegerMatrix() {
0096: final int ans[][] = new int[numRows][numCols];
0097: for (int i = 0; i < numRows; i++) {
0098: for (int j = 0; j < numCols; j++)
0099: ans[i][j] = Math.round((float) matrix[i][j]);
0100: }
0101: return new IntegerSquareMatrix(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 double 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, double 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 double infNorm() {
0152: double result = 0;
0153: for (int i = 0; i < numRows; i++) {
0154: double 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 double 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 det * LUpivot[numRows];
0190: }
0191: }
0192:
0193: /**
0194: * Returns the trace.
0195: */
0196: public double trace() {
0197: double 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 double array[][] = new double[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 DoubleSquareMatrix(array);
0218: }
0219:
0220: // ADDITION
0221:
0222: /**
0223: * Returns the addition of this matrix and another.
0224: * @param m a double matrix
0225: * @exception MatrixDimensionException If the matrices are different sizes.
0226: */
0227: public AbstractDoubleSquareMatrix add(
0228: final AbstractDoubleSquareMatrix m) {
0229: if (m instanceof DoubleSquareMatrix)
0230: return add((DoubleSquareMatrix) m);
0231:
0232: if (numRows == m.rows() && numCols == m.columns()) {
0233: final double array[][] = new double[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 DoubleSquareMatrix(array);
0240: } else {
0241: throw new MatrixDimensionException(
0242: "Matrices are different sizes.");
0243: }
0244: }
0245:
0246: public DoubleSquareMatrix add(final DoubleSquareMatrix m) {
0247: if (numRows == m.numRows && numCols == m.numCols) {
0248: final double array[][] = new double[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 DoubleSquareMatrix(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 double matrix
0265: * @exception MatrixDimensionException If the matrices are different sizes.
0266: */
0267: public AbstractDoubleSquareMatrix subtract(
0268: final AbstractDoubleSquareMatrix m) {
0269: if (m instanceof DoubleSquareMatrix)
0270: return subtract((DoubleSquareMatrix) m);
0271:
0272: if (numRows == m.rows() && numCols == m.columns()) {
0273: final double array[][] = new double[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 DoubleSquareMatrix(array);
0280: } else {
0281: throw new MatrixDimensionException(
0282: "Matrices are different sizes.");
0283: }
0284: }
0285:
0286: public DoubleSquareMatrix subtract(final DoubleSquareMatrix m) {
0287: if (numRows == m.numRows && numCols == m.numCols) {
0288: final double array[][] = new double[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 DoubleSquareMatrix(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 double.
0305: * @return a double square matrix.
0306: */
0307: public AbstractDoubleMatrix scalarMultiply(final double x) {
0308: final double array[][] = new double[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 DoubleSquareMatrix(array);
0315: }
0316:
0317: // SCALAR DIVISON
0318:
0319: /**
0320: * Returns the division of this matrix by a scalar.
0321: * @param x a double.
0322: * @return a double square matrix.
0323: */
0324: public AbstractDoubleMatrix scalarDivide(final double x) {
0325: final double array[][] = new double[numRows][numCols];
0326: for (int i = 0; i < numRows; i++) {
0327: array[i][0] = matrix[i][0] / x;
0328: for (int j = 1; j < numCols; j++)
0329: array[i][j] = matrix[i][j] / x;
0330: }
0331: return new DoubleSquareMatrix(array);
0332: }
0333:
0334: // SCALAR PRODUCT
0335:
0336: /**
0337: * Returns the scalar product of this matrix and another.
0338: * @param m a double matrix.
0339: * @exception MatrixDimensionException If the matrices are different sizes.
0340: */
0341: public double scalarProduct(final AbstractDoubleSquareMatrix m) {
0342: if (m instanceof DoubleSquareMatrix)
0343: return scalarProduct((DoubleSquareMatrix) m);
0344:
0345: if (numRows == m.rows() && numCols == m.columns()) {
0346: double ans = 0;
0347: for (int i = 0; i < numRows; i++) {
0348: ans += matrix[i][0] * m.getElement(i, 0);
0349: for (int j = 1; j < numCols; j++)
0350: ans += matrix[i][j] * m.getElement(i, j);
0351: }
0352: return ans;
0353: } else {
0354: throw new MatrixDimensionException(
0355: "Matrices are different sizes.");
0356: }
0357: }
0358:
0359: public double scalarProduct(final DoubleSquareMatrix m) {
0360: if (numRows == m.numRows && numCols == m.numCols) {
0361: double ans = 0;
0362: for (int i = 0; i < numRows; i++) {
0363: ans += matrix[i][0] * m.matrix[i][0];
0364: for (int j = 1; j < numCols; j++)
0365: ans += matrix[i][j] * m.matrix[i][j];
0366: }
0367: return ans;
0368: } else
0369: throw new MatrixDimensionException(
0370: "Matrices are different sizes.");
0371: }
0372:
0373: // MATRIX MULTIPLICATION
0374:
0375: /**
0376: * Returns the multiplication of a vector by this matrix.
0377: * @param v a double vector.
0378: * @exception DimensionException If the matrix and vector are incompatible.
0379: */
0380: public AbstractDoubleVector multiply(final AbstractDoubleVector v) {
0381: if (numCols == v.dimension()) {
0382: final double array[] = new double[numRows];
0383: for (int i = 0; i < numRows; i++) {
0384: array[i] = matrix[i][0] * v.getComponent(0);
0385: for (int j = 1; j < numCols; j++)
0386: array[i] += matrix[i][j] * v.getComponent(j);
0387: }
0388: return new DoubleVector(array);
0389: } else {
0390: throw new DimensionException(
0391: "Matrix and vector are incompatible.");
0392: }
0393: }
0394:
0395: /**
0396: * Returns the multiplication of this matrix and another.
0397: * @param m a double matrix
0398: * @return a AbstractDoubleMatrix or a AbstractDoubleSquareMatrix as appropriate
0399: * @exception MatrixDimensionException If the matrices are incompatible.
0400: */
0401: public AbstractDoubleSquareMatrix multiply(
0402: final AbstractDoubleSquareMatrix m) {
0403: if (m instanceof DoubleSquareMatrix)
0404: return multiply((DoubleSquareMatrix) m);
0405:
0406: if (numCols == m.rows()) {
0407: final int mColumns = m.columns();
0408: final double array[][] = new double[numRows][mColumns];
0409: for (int j = 0; j < numRows; j++) {
0410: for (int k = 0; k < mColumns; k++) {
0411: array[j][k] = matrix[j][0] * m.getElement(0, k);
0412: for (int n = 1; n < numCols; n++)
0413: array[j][k] += matrix[j][n]
0414: * m.getElement(n, k);
0415: }
0416: }
0417: return new DoubleSquareMatrix(array);
0418: } else {
0419: throw new MatrixDimensionException("Incompatible matrices.");
0420: }
0421: }
0422:
0423: public DoubleSquareMatrix multiply(final DoubleSquareMatrix m) {
0424: if (numCols == m.numRows) {
0425: final double array[][] = new double[numRows][m.numCols];
0426: for (int j = 0; j < numRows; j++) {
0427: for (int k = 0; k < m.numCols; k++) {
0428: array[j][k] = matrix[j][0] * m.matrix[0][k];
0429: for (int n = 1; n < numCols; n++)
0430: array[j][k] += matrix[j][n] * m.matrix[n][k];
0431: }
0432: }
0433: return new DoubleSquareMatrix(array);
0434: } else
0435: throw new MatrixDimensionException("Incompatible matrices.");
0436: }
0437:
0438: // DIRECT SUM
0439:
0440: /**
0441: * Returns the direct sum of this matrix and another.
0442: */
0443: public AbstractDoubleSquareMatrix directSum(
0444: final AbstractDoubleSquareMatrix m) {
0445: final double array[][] = new double[numRows + m.numRows][numCols
0446: + m.numCols];
0447: for (int i = 0; i < numRows; i++) {
0448: for (int j = 0; j < numCols; j++)
0449: array[i][j] = matrix[i][j];
0450: }
0451: for (int i = 0; i < m.numRows; i++) {
0452: for (int j = 0; j < m.numCols; j++)
0453: array[i + numRows][j + numCols] = m.getElement(i, j);
0454: }
0455: return new DoubleSquareMatrix(array);
0456: }
0457:
0458: // TENSOR PRODUCT
0459:
0460: /**
0461: * Returns the tensor product of this matrix and another.
0462: */
0463: public AbstractDoubleSquareMatrix tensor(
0464: final AbstractDoubleSquareMatrix m) {
0465: final double array[][] = new double[numRows * m.numRows][numCols
0466: * m.numCols];
0467: for (int i = 0; i < numRows; i++) {
0468: for (int j = 0; j < numCols; j++) {
0469: for (int k = 0; k < m.numRows; j++) {
0470: for (int l = 0; l < m.numCols; l++)
0471: array[i * m.numRows + k][j * m.numCols + l] = matrix[i][j]
0472: * m.getElement(k, l);
0473: }
0474: }
0475: }
0476: return new DoubleSquareMatrix(array);
0477: }
0478:
0479: // TRANSPOSE
0480:
0481: /**
0482: * Returns the transpose of this matrix.
0483: * @return a double matrix
0484: */
0485: public Matrix transpose() {
0486: final double array[][] = new double[numCols][numRows];
0487: for (int i = 0; i < numRows; i++) {
0488: array[0][i] = matrix[i][0];
0489: for (int j = 1; j < numCols; j++)
0490: array[j][i] = matrix[i][j];
0491: }
0492: return new DoubleSquareMatrix(array);
0493: }
0494:
0495: // INVERSE
0496:
0497: /**
0498: * Returns the inverse of this matrix.
0499: * @return a double square matrix.
0500: */
0501: public AbstractDoubleSquareMatrix inverse() {
0502: final int N = numRows;
0503: final double arrayL[][] = new double[N][N];
0504: final double arrayU[][] = new double[N][N];
0505: final DoubleSquareMatrix lu[] = (DoubleSquareMatrix[]) this
0506: .luDecompose(null);
0507: arrayL[0][0] = 1.0 / lu[0].matrix[0][0];
0508: arrayU[0][0] = 1.0 / lu[1].matrix[0][0];
0509: for (int i = 1; i < N; i++) {
0510: arrayL[i][i] = 1.0 / lu[0].matrix[i][i];
0511: arrayU[i][i] = 1.0 / lu[1].matrix[i][i];
0512: }
0513: for (int i = 0; i < N - 1; i++) {
0514: for (int j = i + 1; j < N; j++) {
0515: double tmpL = 0.0, tmpU = 0.0;
0516: for (int k = i; k < j; k++) {
0517: tmpL -= lu[0].matrix[j][k] * arrayL[k][i];
0518: tmpU -= arrayU[i][k] * lu[1].matrix[k][j];
0519: }
0520: arrayL[j][i] = tmpL / lu[0].matrix[j][j];
0521: arrayU[i][j] = tmpU / lu[1].matrix[j][j];
0522: }
0523: }
0524: // matrix multiply arrayU x arrayL
0525: final double inv[][] = new double[N][N];
0526: for (int i = 0; i < N; i++) {
0527: for (int j = 0; j < i; j++) {
0528: for (int k = i; k < N; k++)
0529: inv[i][LUpivot[j]] += arrayU[i][k] * arrayL[k][j];
0530: }
0531: for (int j = i; j < N; j++) {
0532: for (int k = j; k < N; k++)
0533: inv[i][LUpivot[j]] += arrayU[i][k] * arrayL[k][j];
0534: }
0535: }
0536: return new DoubleSquareMatrix(inv);
0537: }
0538:
0539: // LU DECOMPOSITION
0540:
0541: /**
0542: * Returns the LU decomposition of this matrix.
0543: * @param pivot an empty array of length <code>rows()+1</code>
0544: * to hold the pivot information (null if not interested).
0545: * The last array element will contain the parity.
0546: * @return an array with [0] containing the L-matrix
0547: * and [1] containing the U-matrix.
0548: * @jsci.planetmath LUDecomposition
0549: */
0550: public final AbstractDoubleSquareMatrix[] luDecompose(int pivot[]) {
0551: if (LU != null) {
0552: if (pivot != null)
0553: System.arraycopy(LUpivot, 0, pivot, 0, pivot.length);
0554: return LU;
0555: }
0556: int pivotrow;
0557: final int N = numRows;
0558: final double arrayL[][] = new double[N][N];
0559: final double arrayU[][] = new double[N][N];
0560: double tmp;
0561: if (pivot == null)
0562: pivot = new int[N + 1];
0563: for (int i = 0; i < N; i++)
0564: pivot[i] = i;
0565: pivot[N] = 1;
0566: // LU decomposition to arrayU
0567: for (int j = 0; j < N; j++) {
0568: for (int i = 0; i < j; i++) {
0569: tmp = matrix[pivot[i]][j];
0570: for (int k = 0; k < i; k++)
0571: tmp -= arrayU[i][k] * arrayU[k][j];
0572: arrayU[i][j] = tmp;
0573: }
0574: double max = 0.0;
0575: pivotrow = j;
0576: for (int i = j; i < N; i++) {
0577: tmp = matrix[pivot[i]][j];
0578: for (int k = 0; k < j; k++)
0579: tmp -= arrayU[i][k] * arrayU[k][j];
0580: arrayU[i][j] = tmp;
0581: // while we're here search for a pivot for arrayU[j][j]
0582: tmp = Math.abs(tmp);
0583: if (tmp > max) {
0584: max = tmp;
0585: pivotrow = i;
0586: }
0587: }
0588: // swap row j with pivotrow
0589: if (pivotrow != j) {
0590: double[] tmprow = arrayU[j];
0591: arrayU[j] = arrayU[pivotrow];
0592: arrayU[pivotrow] = tmprow;
0593: int k = pivot[j];
0594: pivot[j] = pivot[pivotrow];
0595: pivot[pivotrow] = k;
0596: // update parity
0597: pivot[N] = -pivot[N];
0598: }
0599: // divide by pivot
0600: tmp = arrayU[j][j];
0601: for (int i = j + 1; i < N; i++)
0602: arrayU[i][j] /= tmp;
0603: }
0604: // move lower triangular part to arrayL
0605: for (int j = 0; j < N; j++) {
0606: arrayL[j][j] = 1.0;
0607: for (int i = j + 1; i < N; i++) {
0608: arrayL[i][j] = arrayU[i][j];
0609: arrayU[i][j] = 0.0;
0610: }
0611: }
0612: LU = new DoubleSquareMatrix[2];
0613: LU[0] = new DoubleSquareMatrix(arrayL);
0614: LU[1] = new DoubleSquareMatrix(arrayU);
0615: LUpivot = new int[pivot.length];
0616: System.arraycopy(pivot, 0, LUpivot, 0, pivot.length);
0617: return LU;
0618: }
0619:
0620: /**
0621: * Returns the LU decomposition of this matrix.
0622: * Warning: no pivoting.
0623: * @return an array with [0] containing the L-matrix
0624: * and [1] containing the U-matrix.
0625: * @jsci.planetmath LUDecomposition
0626: */
0627: public final AbstractDoubleSquareMatrix[] luDecompose() {
0628: final int N = numRows;
0629: final double arrayL[][] = new double[N][N];
0630: final double arrayU[][] = new double[N][N];
0631: double tmp;
0632: // LU decomposition to arrayU
0633: for (int j = 0; j < N; j++) {
0634: for (int i = 0; i < j; i++) {
0635: tmp = matrix[i][j];
0636: for (int k = 0; k < i; k++)
0637: tmp -= arrayU[i][k] * arrayU[k][j];
0638: arrayU[i][j] = tmp;
0639: }
0640: for (int i = j; i < N; i++) {
0641: tmp = matrix[i][j];
0642: for (int k = 0; k < j; k++)
0643: tmp -= arrayU[i][k] * arrayU[k][j];
0644: arrayU[i][j] = tmp;
0645: }
0646: // divide
0647: tmp = arrayU[j][j];
0648: for (int i = j + 1; i < N; i++)
0649: arrayU[i][j] /= tmp;
0650: }
0651: // move lower triangular part to arrayL
0652: for (int j = 0; j < N; j++) {
0653: arrayL[j][j] = 1.0;
0654: for (int i = j + 1; i < N; i++) {
0655: arrayL[i][j] = arrayU[i][j];
0656: arrayU[i][j] = 0.0;
0657: }
0658: }
0659: DoubleSquareMatrix[] lu = new DoubleSquareMatrix[2];
0660: lu[0] = new DoubleSquareMatrix(arrayL);
0661: lu[1] = new DoubleSquareMatrix(arrayU);
0662: return lu;
0663: }
0664:
0665: // CHOLESKY DECOMPOSITION
0666:
0667: /**
0668: * Returns the Cholesky decomposition of this matrix.
0669: * Matrix must be symmetric and positive definite.
0670: * @return an array with [0] containing the L-matrix and [1] containing the U-matrix.
0671: * @jsci.planetmath CholeskyDecomposition
0672: */
0673: public AbstractDoubleSquareMatrix[] choleskyDecompose() {
0674: final int N = numRows;
0675: final double arrayL[][] = new double[N][N];
0676: final double arrayU[][] = new double[N][N];
0677: double tmp = Math.sqrt(matrix[0][0]);
0678: arrayL[0][0] = arrayU[0][0] = tmp;
0679: for (int i = 1; i < N; i++)
0680: arrayL[i][0] = arrayU[0][i] = matrix[i][0] / tmp;
0681: for (int j = 1; j < N; j++) {
0682: tmp = matrix[j][j];
0683: for (int i = 0; i < j; i++)
0684: tmp -= arrayL[j][i] * arrayL[j][i];
0685: arrayL[j][j] = arrayU[j][j] = Math.sqrt(tmp);
0686: for (int i = j + 1; i < N; i++) {
0687: tmp = matrix[i][j];
0688: for (int k = 0; k < i; k++)
0689: tmp -= arrayL[j][k] * arrayU[k][i];
0690: arrayL[i][j] = arrayU[j][i] = tmp / arrayU[j][j];
0691: }
0692: }
0693: final AbstractDoubleSquareMatrix lu[] = new AbstractDoubleSquareMatrix[2];
0694: lu[0] = new DoubleSquareMatrix(arrayL);
0695: lu[1] = new DoubleSquareMatrix(arrayU);
0696: return lu;
0697: }
0698:
0699: // QR DECOMPOSITION
0700:
0701: /**
0702: * Returns the QR decomposition of this matrix.
0703: * Based on the code from <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> (public domain).
0704: * @return an array with [0] containing the Q-matrix and [1] containing the R-matrix.
0705: * @jsci.planetmath QRDecomposition
0706: */
0707: public AbstractDoubleSquareMatrix[] qrDecompose() {
0708: final int N = numRows;
0709: final double array[][] = new double[N][N];
0710: final double arrayQ[][] = new double[N][N];
0711: final double arrayR[][] = new double[N][N];
0712: // copy matrix
0713: for (int i = 0; i < N; i++) {
0714: array[i][0] = matrix[i][0];
0715: for (int j = 1; j < N; j++)
0716: array[i][j] = matrix[i][j];
0717: }
0718: for (int k = 0; k < N; k++) {
0719: // compute l2-norm of kth column
0720: double norm = array[k][k];
0721: for (int i = k + 1; i < N; i++)
0722: norm = ExtraMath.hypot(norm, array[i][k]);
0723: if (norm != 0.0) {
0724: // form kth Householder vector
0725: if (array[k][k] < 0.0)
0726: norm = -norm;
0727: for (int i = k; i < N; i++)
0728: array[i][k] /= norm;
0729: array[k][k] += 1.0;
0730: // apply transformation to remaining columns
0731: for (int j = k + 1; j < N; j++) {
0732: double s = array[k][k] * array[k][j];
0733: for (int i = k + 1; i < N; i++)
0734: s += array[i][k] * array[i][j];
0735: s /= array[k][k];
0736: for (int i = k; i < N; i++)
0737: array[i][j] -= s * array[i][k];
0738: }
0739: }
0740: arrayR[k][k] = -norm;
0741: }
0742: for (int k = N - 1; k >= 0; k--) {
0743: arrayQ[k][k] = 1.0;
0744: for (int j = k; j < N; j++) {
0745: if (array[k][k] != 0.0) {
0746: double s = array[k][k] * arrayQ[k][j];
0747: for (int i = k + 1; i < N; i++)
0748: s += array[i][k] * arrayQ[i][j];
0749: s /= array[k][k];
0750: for (int i = k; i < N; i++)
0751: arrayQ[i][j] -= s * array[i][k];
0752: }
0753: }
0754: }
0755: for (int i = 0; i < N; i++) {
0756: for (int j = i + 1; j < N; j++)
0757: arrayR[i][j] = array[i][j];
0758: }
0759: final AbstractDoubleSquareMatrix qr[] = new AbstractDoubleSquareMatrix[2];
0760: qr[0] = new DoubleSquareMatrix(arrayQ);
0761: qr[1] = new DoubleSquareMatrix(arrayR);
0762: return qr;
0763: }
0764:
0765: // SINGULAR VALUE DECOMPOSITION
0766:
0767: /**
0768: * Returns the singular value decomposition of this matrix.
0769: * Based on the code from <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> (public domain).
0770: * @return an array with [0] containing the U-matrix, [1] containing the S-matrix and [2] containing the V-matrix.
0771: * @jsci.planetmath SingularValueDecomposition
0772: */
0773: public AbstractDoubleSquareMatrix[] singularValueDecompose() {
0774: final int N = numRows;
0775: final int Nm1 = N - 1;
0776: final double array[][] = new double[N][N];
0777: final double arrayU[][] = new double[N][N];
0778: final double arrayS[] = new double[N];
0779: final double arrayV[][] = new double[N][N];
0780: final double e[] = new double[N];
0781: final double work[] = new double[N];
0782: // copy matrix
0783: for (int i = 0; i < N; i++) {
0784: array[i][0] = matrix[i][0];
0785: for (int j = 1; j < N; j++)
0786: array[i][j] = matrix[i][j];
0787: }
0788: // reduce matrix to bidiagonal form
0789: for (int k = 0; k < Nm1; k++) {
0790: // compute the transformation for the kth column
0791: // compute l2-norm of kth column
0792: arrayS[k] = array[k][k];
0793: for (int i = k + 1; i < N; i++)
0794: arrayS[k] = ExtraMath.hypot(arrayS[k], array[i][k]);
0795: if (arrayS[k] != 0.0) {
0796: if (array[k][k] < 0.0)
0797: arrayS[k] = -arrayS[k];
0798: for (int i = k; i < N; i++)
0799: array[i][k] /= arrayS[k];
0800: array[k][k] += 1.0;
0801: }
0802: arrayS[k] = -arrayS[k];
0803: for (int j = k + 1; j < N; j++) {
0804: if (arrayS[k] != 0.0) {
0805: // apply the transformation
0806: double t = array[k][k] * array[k][j];
0807: for (int i = k + 1; i < N; i++)
0808: t += array[i][k] * array[i][j];
0809: t /= array[k][k];
0810: for (int i = k; i < N; i++)
0811: array[i][j] -= t * array[i][k];
0812: }
0813: e[j] = array[k][j];
0814: }
0815: for (int i = k; i < N; i++)
0816: arrayU[i][k] = array[i][k];
0817: if (k < N - 2) {
0818: // compute the kth row transformation
0819: // compute l2-norm of kth column
0820: e[k] = e[k + 1];
0821: for (int i = k + 2; i < N; i++)
0822: e[k] = ExtraMath.hypot(e[k], e[i]);
0823: if (e[k] != 0.0) {
0824: if (e[k + 1] < 0.0)
0825: e[k] = -e[k];
0826: for (int i = k + 1; i < N; i++)
0827: e[i] /= e[k];
0828: e[k + 1] += 1.0;
0829: }
0830: e[k] = -e[k];
0831: if (e[k] != 0.0) {
0832: // apply the transformation
0833: for (int i = k + 1; i < N; i++) {
0834: work[i] = 0.0;
0835: for (int j = k + 1; j < N; j++)
0836: work[i] += e[j] * array[i][j];
0837: }
0838: for (int j = k + 1; j < N; j++) {
0839: double t = e[j] / e[k + 1];
0840: for (int i = k + 1; i < N; i++)
0841: array[i][j] -= t * work[i];
0842: }
0843: }
0844: for (int i = k + 1; i < N; i++)
0845: arrayV[i][k] = e[i];
0846: }
0847: }
0848: // setup the final bidiagonal matrix of order p
0849: int p = N;
0850: arrayS[Nm1] = array[Nm1][Nm1];
0851: e[N - 2] = array[N - 2][Nm1];
0852: e[Nm1] = 0.0;
0853: arrayU[Nm1][Nm1] = 1.0;
0854: for (int k = N - 2; k >= 0; k--) {
0855: if (arrayS[k] != 0.0) {
0856: for (int j = k + 1; j < N; j++) {
0857: double t = arrayU[k][k] * arrayU[k][j];
0858: for (int i = k + 1; i < N; i++)
0859: t += arrayU[i][k] * arrayU[i][j];
0860: t /= arrayU[k][k];
0861: for (int i = k; i < N; i++)
0862: arrayU[i][j] -= t * arrayU[i][k];
0863: }
0864: for (int i = k; i < N; i++)
0865: arrayU[i][k] = -arrayU[i][k];
0866: arrayU[k][k] += 1.0;
0867: for (int i = 0; i < k - 1; i++)
0868: arrayU[i][k] = 0.0;
0869: } else {
0870: for (int i = 0; i < N; i++)
0871: arrayU[i][k] = 0.0;
0872: arrayU[k][k] = 1.0;
0873: }
0874: }
0875: for (int k = Nm1; k >= 0; k--) {
0876: if (k < N - 2 && e[k] != 0.0) {
0877: for (int j = k + 1; j < N; j++) {
0878: double t = arrayV[k + 1][k] * arrayV[k + 1][j];
0879: for (int i = k + 2; i < N; i++)
0880: t += arrayV[i][k] * arrayV[i][j];
0881: t /= arrayV[k + 1][k];
0882: for (int i = k + 1; i < N; i++)
0883: arrayV[i][j] -= t * arrayV[i][k];
0884: }
0885: }
0886: for (int i = 0; i < N; i++)
0887: arrayV[i][k] = 0.0;
0888: arrayV[k][k] = 1.0;
0889: }
0890: final double eps = Math.pow(2.0, -52.0);
0891: int iter = 0;
0892: while (p > 0) {
0893: int k, action;
0894: // action = 1 if arrayS[p] and e[k-1] are negligible and k<p
0895: // action = 2 if arrayS[k] is negligible and k<p
0896: // action = 3 if e[k-1] is negligible, k<p, and arrayS[k], ..., arrayS[p] are not negligible (QR step)
0897: // action = 4 if e[p-1] is negligible (convergence)
0898: for (k = p - 2; k >= -1; k--) {
0899: if (k == -1)
0900: break;
0901: if (Math.abs(e[k]) <= eps
0902: * (Math.abs(arrayS[k]) + Math
0903: .abs(arrayS[k + 1]))) {
0904: e[k] = 0.0;
0905: break;
0906: }
0907: }
0908: if (k == p - 2) {
0909: action = 4;
0910: } else {
0911: int ks;
0912: for (ks = p - 1; ks >= k; ks--) {
0913: if (ks == k)
0914: break;
0915: double t = (ks != p ? Math.abs(e[ks]) : 0.0)
0916: + (ks != k + 1 ? Math.abs(e[ks - 1]) : 0.0);
0917: if (Math.abs(arrayS[ks]) <= eps * t) {
0918: arrayS[ks] = 0.0;
0919: break;
0920: }
0921: }
0922: if (ks == k) {
0923: action = 3;
0924: } else if (ks == p - 1) {
0925: action = 1;
0926: } else {
0927: action = 2;
0928: k = ks;
0929: }
0930: }
0931: k++;
0932: switch (action) {
0933: // deflate negligible arrayS[p]
0934: case 1: {
0935: double f = e[p - 2];
0936: e[p - 2] = 0.0;
0937: for (int j = p - 2; j >= k; j--) {
0938: double t = ExtraMath.hypot(arrayS[j], f);
0939: final double cs = arrayS[j] / t;
0940: final double sn = f / t;
0941: arrayS[j] = t;
0942: if (j != k) {
0943: f = -sn * e[j - 1];
0944: e[j - 1] *= cs;
0945: }
0946: for (int i = 0; i < N; i++) {
0947: t = cs * arrayV[i][j] + sn * arrayV[i][p - 1];
0948: arrayV[i][p - 1] = -sn * arrayV[i][j] + cs
0949: * arrayV[i][p - 1];
0950: arrayV[i][j] = t;
0951: }
0952: }
0953: }
0954: break;
0955: // split at negligible arrayS[k]
0956: case 2: {
0957: double f = e[k - 1];
0958: e[k - 1] = 0.0;
0959: for (int j = k; j < p; j++) {
0960: double t = ExtraMath.hypot(arrayS[j], f);
0961: final double cs = arrayS[j] / t;
0962: final double sn = f / t;
0963: arrayS[j] = t;
0964: f = -sn * e[j];
0965: e[j] *= cs;
0966: for (int i = 0; i < N; i++) {
0967: t = cs * arrayU[i][j] + sn * arrayU[i][k - 1];
0968: arrayU[i][k - 1] = -sn * arrayU[i][j] + cs
0969: * arrayU[i][k - 1];
0970: arrayU[i][j] = t;
0971: }
0972: }
0973: }
0974: break;
0975: // perform one QR step
0976: case 3: {
0977: // calculate the shift
0978: final double scale = Math.max(Math.max(Math.max(Math
0979: .max(Math.abs(arrayS[p - 1]), Math
0980: .abs(arrayS[p - 2])), Math
0981: .abs(e[p - 2])), Math.abs(arrayS[k])), Math
0982: .abs(e[k]));
0983: double sp = arrayS[p - 1] / scale;
0984: double spm1 = arrayS[p - 2] / scale;
0985: double epm1 = e[p - 2] / scale;
0986: double sk = arrayS[k] / scale;
0987: double ek = e[k] / scale;
0988: double b = ((spm1 + sp) * (spm1 - sp) + epm1 * epm1) / 2.0;
0989: double c = (sp * epm1) * (sp * epm1);
0990: double shift = 0.0;
0991: if (b != 0.0 || c != 0.0) {
0992: shift = Math.sqrt(b * b + c);
0993: if (b < 0.0)
0994: shift = -shift;
0995: shift = c / (b + shift);
0996: }
0997: double f = (sk + sp) * (sk - sp) + shift;
0998: double g = sk * ek;
0999: // chase zeros
1000: for (int j = k; j < p - 1; j++) {
1001: double t = ExtraMath.hypot(f, g);
1002: double cs = f / t;
1003: double sn = g / t;
1004: if (j != k)
1005: e[j - 1] = t;
1006: f = cs * arrayS[j] + sn * e[j];
1007: e[j] = cs * e[j] - sn * arrayS[j];
1008: g = sn * arrayS[j + 1];
1009: arrayS[j + 1] *= cs;
1010: for (int i = 0; i < N; i++) {
1011: t = cs * arrayV[i][j] + sn * arrayV[i][j + 1];
1012: arrayV[i][j + 1] = -sn * arrayV[i][j] + cs
1013: * arrayV[i][j + 1];
1014: arrayV[i][j] = t;
1015: }
1016: t = ExtraMath.hypot(f, g);
1017: cs = f / t;
1018: sn = g / t;
1019: arrayS[j] = t;
1020: f = cs * e[j] + sn * arrayS[j + 1];
1021: arrayS[j + 1] = -sn * e[j] + cs * arrayS[j + 1];
1022: g = sn * e[j + 1];
1023: e[j + 1] *= cs;
1024: if (j < Nm1) {
1025: for (int i = 0; i < N; i++) {
1026: t = cs * arrayU[i][j] + sn
1027: * arrayU[i][j + 1];
1028: arrayU[i][j + 1] = -sn * arrayU[i][j] + cs
1029: * arrayU[i][j + 1];
1030: arrayU[i][j] = t;
1031: }
1032: }
1033: }
1034: e[p - 2] = f;
1035: iter++;
1036: }
1037: break;
1038: // convergence
1039: case 4: {
1040: // make the singular values positive
1041: if (arrayS[k] <= 0.0) {
1042: arrayS[k] = -arrayS[k];
1043: for (int i = 0; i < p; i++)
1044: arrayV[i][k] = -arrayV[i][k];
1045: }
1046: // order the singular values
1047: while (k < p - 1) {
1048: if (arrayS[k] >= arrayS[k + 1])
1049: break;
1050: double tmp = arrayS[k];
1051: arrayS[k] = arrayS[k + 1];
1052: arrayS[k + 1] = tmp;
1053: if (k < Nm1) {
1054: for (int i = 0; i < N; i++) {
1055: tmp = arrayU[i][k + 1];
1056: arrayU[i][k + 1] = arrayU[i][k];
1057: arrayU[i][k] = tmp;
1058: tmp = arrayV[i][k + 1];
1059: arrayV[i][k + 1] = arrayV[i][k];
1060: arrayV[i][k] = tmp;
1061: }
1062: }
1063: k++;
1064: }
1065: iter = 0;
1066: p--;
1067: }
1068: break;
1069: }
1070: }
1071: final AbstractDoubleSquareMatrix svd[] = new AbstractDoubleSquareMatrix[3];
1072: svd[0] = new DoubleSquareMatrix(arrayU);
1073: svd[1] = new DoubleDiagonalMatrix(arrayS);
1074: svd[2] = new DoubleSquareMatrix(arrayV);
1075: return svd;
1076: }
1077:
1078: // POLAR DECOMPOSITION
1079:
1080: /**
1081: * Returns the polar decomposition of this matrix.
1082: */
1083: public AbstractDoubleSquareMatrix[] polarDecompose() {
1084: final int N = numRows;
1085: final AbstractDoubleVector evec[] = new AbstractDoubleVector[N];
1086: double eval[];
1087: try {
1088: eval = LinearMath.eigenSolveSymmetric(this , evec);
1089: } catch (MaximumIterationsExceededException e) {
1090: return null;
1091: }
1092: final double tmpa[][] = new double[N][N];
1093: final double tmpm[][] = new double[N][N];
1094: double abs;
1095: for (int i = 0; i < N; i++) {
1096: abs = Math.abs(eval[i]);
1097: tmpa[i][0] = eval[i] * evec[i].getComponent(0) / abs;
1098: tmpm[i][0] = abs * evec[i].getComponent(0);
1099: for (int j = 1; j < N; j++) {
1100: tmpa[i][j] = eval[i] * evec[i].getComponent(j) / abs;
1101: tmpm[i][j] = abs * evec[i].getComponent(j);
1102: }
1103: }
1104: final double arg[][] = new double[N][N];
1105: final double mod[][] = new double[N][N];
1106: for (int i = 0; i < N; i++) {
1107: for (int j = 0; j < N; j++) {
1108: arg[i][j] = evec[0].getComponent(i) * tmpa[0][j];
1109: mod[i][j] = evec[0].getComponent(i) * tmpm[0][j];
1110: for (int k = 1; k < N; k++) {
1111: arg[i][j] += evec[k].getComponent(i) * tmpa[k][j];
1112: mod[i][j] += evec[k].getComponent(i) * tmpm[k][j];
1113: }
1114: }
1115: }
1116: final AbstractDoubleSquareMatrix us[] = new AbstractDoubleSquareMatrix[2];
1117: us[0] = new DoubleSquareMatrix(arg);
1118: us[1] = new DoubleSquareMatrix(mod);
1119: return us;
1120: }
1121:
1122: // MAP ELEMENTS
1123:
1124: /**
1125: * Applies a function on all the matrix elements.
1126: * @param f a user-defined function
1127: * @return a double matrix
1128: */
1129: public AbstractDoubleMatrix mapElements(final Mapping f) {
1130: final double array[][] = new double[numRows][numCols];
1131: for (int i = 0; i < numRows; i++) {
1132: array[i][0] = f.map(matrix[i][0]);
1133: for (int j = 1; j < numCols; j++)
1134: array[i][j] = f.map(matrix[i][j]);
1135: }
1136: return new DoubleSquareMatrix(array);
1137: }
1138: }
|