0001: /* AUTO-GENERATED */
0002: package JSci.maths.matrices;
0003:
0004: import JSci.maths.ArrayMath;
0005: import JSci.maths.ExtraMath;
0006: import JSci.maths.Mapping;
0007: import JSci.maths.LinearMath;
0008: import JSci.maths.MaximumIterationsExceededException;
0009: import JSci.maths.DimensionException;
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 DoubleTridiagonalMatrix class provides an object for encapsulating double tridiagonal matrices.
0018: * @version 2.2
0019: * @author Mark Hale
0020: */
0021: public class DoubleTridiagonalMatrix extends AbstractDoubleSquareMatrix
0022: implements TridiagonalMatrix {
0023: /**
0024: * Tridiagonal data.
0025: */
0026: protected final double ldiag[];
0027: protected final double diag[];
0028: protected final double udiag[];
0029:
0030: /**
0031: * Constructs an empty matrix.
0032: * @param size the number of rows/columns
0033: */
0034: public DoubleTridiagonalMatrix(final int size) {
0035: super (size);
0036: ldiag = new double[size];
0037: diag = new double[size];
0038: udiag = new double[size];
0039: }
0040:
0041: /**
0042: * Constructs a matrix from an array.
0043: * Any non-tridiagonal elements in the array are ignored.
0044: * @param array an assigned value
0045: * @exception MatrixDimensionException If the array is not square.
0046: */
0047: public DoubleTridiagonalMatrix(final double array[][]) {
0048: this (array.length);
0049: if (!ArrayMath.isSquare(array))
0050: throw new MatrixDimensionException("Array is not square.");
0051: diag[0] = array[0][0];
0052: udiag[0] = array[0][1];
0053: int i = 1;
0054: for (; i < array.length - 1; i++) {
0055: ldiag[i] = array[i][i - 1];
0056: diag[i] = array[i][i];
0057: udiag[i] = array[i][i + 1];
0058: }
0059: ldiag[i] = array[i][i - 1];
0060: diag[i] = array[i][i];
0061: }
0062:
0063: /**
0064: * Compares two ${nativeTyp} matrices for equality.
0065: * @param m a double matrix
0066: */
0067: public boolean equals(AbstractDoubleMatrix m, double tol) {
0068: if (m instanceof TridiagonalMatrix) {
0069: if (numRows != m.rows() || numCols != m.columns())
0070: return false;
0071: double sumSqr = 0;
0072: double ldelta;
0073: double delta = diag[0] - m.getElement(0, 0);
0074: double udelta = udiag[0] - m.getElement(0, 1);
0075: sumSqr += delta * delta + udelta * udelta;
0076: int i = 1;
0077: for (; i < numRows - 1; i++) {
0078: ldelta = ldiag[i] - m.getElement(i, i - 1);
0079: delta = diag[i] - m.getElement(i, i);
0080: udelta = udiag[i] - m.getElement(i, i + 1);
0081: sumSqr += ldelta * ldelta + delta * delta + udelta
0082: * udelta;
0083: }
0084: ldelta = ldiag[i] - m.getElement(i, i - 1);
0085: delta = diag[i] - m.getElement(i, i);
0086: sumSqr += ldelta * ldelta + delta * delta;
0087: return (sumSqr <= tol * tol);
0088: } else {
0089: return false;
0090: }
0091: }
0092:
0093: /**
0094: * Returns a string representing this matrix.
0095: */
0096: public String toString() {
0097: final StringBuffer buf = new StringBuffer(5 * numRows * numCols);
0098: for (int i = 0; i < numRows; i++) {
0099: for (int j = 0; j < numCols; j++) {
0100: buf.append(getElement(i, j));
0101: buf.append(' ');
0102: }
0103: buf.append('\n');
0104: }
0105: return buf.toString();
0106: }
0107:
0108: /**
0109: * Converts this matrix to an integer matrix.
0110: * @return an integer matrix
0111: */
0112: public AbstractIntegerMatrix toIntegerMatrix() {
0113: final IntegerTridiagonalMatrix m = new IntegerTridiagonalMatrix(
0114: numRows);
0115: m.diag[0] = Math.round((float) diag[0]);
0116: m.udiag[0] = Math.round((float) udiag[0]);
0117: int i = 1;
0118: for (; i < numRows - 1; i++) {
0119: m.ldiag[i] = Math.round((float) ldiag[i]);
0120: m.diag[i] = Math.round((float) diag[i]);
0121: m.udiag[i] = Math.round((float) udiag[i]);
0122: }
0123: m.ldiag[i] = Math.round((float) ldiag[i]);
0124: m.diag[i] = Math.round((float) diag[i]);
0125: return m;
0126: }
0127:
0128: /**
0129: * Converts this matrix to a complex matrix.
0130: * @return a complex matrix
0131: */
0132: public AbstractComplexMatrix toComplexMatrix() {
0133: final ComplexTridiagonalMatrix m = new ComplexTridiagonalMatrix(
0134: numRows);
0135: m.diagRe[0] = diag[0];
0136: m.udiagRe[0] = udiag[0];
0137: int i = 1;
0138: for (; i < numRows - 1; i++) {
0139: m.ldiagRe[i] = ldiag[i];
0140: m.diagRe[i] = diag[i];
0141: m.udiagRe[i] = udiag[i];
0142: }
0143: m.ldiagRe[i] = ldiag[i];
0144: m.diagRe[i] = diag[i];
0145: return m;
0146: }
0147:
0148: /**
0149: * Returns an element of the matrix.
0150: * @param i row index of the element
0151: * @param j column index of the element
0152: * @exception MatrixDimensionException If attempting to access an invalid element.
0153: */
0154: public double getElement(int i, int j) {
0155: if (i >= 0 && i < numRows && j >= 0 && j < numCols) {
0156: if (j == i - 1)
0157: return ldiag[i];
0158: else if (j == i)
0159: return diag[i];
0160: else if (j == i + 1)
0161: return udiag[i];
0162: else
0163: return 0;
0164: } else
0165: throw new MatrixDimensionException(getInvalidElementMsg(i,
0166: j));
0167: }
0168:
0169: /**
0170: * Sets the value of an element of the matrix.
0171: * Should only be used to initialise this matrix.
0172: * @param i row index of the element
0173: * @param j column index of the element
0174: * @param x a number
0175: * @exception MatrixDimensionException If attempting to access an invalid element.
0176: */
0177: public void setElement(int i, int j, final double x) {
0178: if (i >= 0 && i < numRows && j >= 0 && j < numCols) {
0179: if (j == i - 1)
0180: ldiag[i] = x;
0181: else if (j == i)
0182: diag[i] = x;
0183: else if (j == i + 1)
0184: udiag[i] = x;
0185: else
0186: throw new MatrixDimensionException(
0187: getInvalidElementMsg(i, j));
0188: } else
0189: throw new MatrixDimensionException(getInvalidElementMsg(i,
0190: j));
0191: }
0192:
0193: /**
0194: * Returns true if this matrix is symmetric.
0195: */
0196: public boolean isSymmetric() {
0197: if (ldiag[1] != udiag[0])
0198: return false;
0199: for (int i = 1; i < numRows - 1; i++) {
0200: if (ldiag[i + 1] != udiag[i])
0201: return false;
0202: }
0203: return true;
0204: }
0205:
0206: /**
0207: * Returns the trace.
0208: */
0209: public double trace() {
0210: double tr = diag[0];
0211: for (int i = 1; i < numRows; i++)
0212: tr += diag[i];
0213: return tr;
0214: }
0215:
0216: /**
0217: * Returns the l<sup><img border=0 alt="infinity" src="doc-files/infinity.gif"></sup>-norm.
0218: * @author Taber Smith
0219: */
0220: public double infNorm() {
0221: double result = Math.abs(diag[0]) + Math.abs(udiag[0]);
0222: double tmpResult;
0223: int i = 1;
0224: for (; i < numRows - 1; i++) {
0225: tmpResult = Math.abs(ldiag[i]) + Math.abs(diag[i])
0226: + Math.abs(udiag[i]);
0227: if (tmpResult > result)
0228: result = tmpResult;
0229: }
0230: tmpResult = Math.abs(ldiag[i]) + Math.abs(diag[i]);
0231: if (tmpResult > result)
0232: result = tmpResult;
0233: return result;
0234: }
0235:
0236: /**
0237: * Returns the Frobenius (l<sup>2</sup>) norm.
0238: * @author Taber Smith
0239: */
0240: public double frobeniusNorm() {
0241: double result = diag[0] * diag[0] + udiag[0] * udiag[0];
0242: int i = 1;
0243: for (; i < numRows - 1; i++)
0244: result += ldiag[i] * ldiag[i] + diag[i] * diag[i]
0245: + udiag[i] * udiag[i];
0246: result += ldiag[i] * ldiag[i] + diag[i] * diag[i];
0247: return Math.sqrt(result);
0248: }
0249:
0250: /**
0251: * Returns the operator norm.
0252: * @exception MaximumIterationsExceededException If it takes more than 50 iterations to determine an eigenvalue.
0253: */
0254: public double operatorNorm()
0255: throws MaximumIterationsExceededException {
0256: return Math
0257: .sqrt(ArrayMath
0258: .max(LinearMath
0259: .eigenvalueSolveSymmetric((DoubleTridiagonalMatrix) (this
0260: .transpose().multiply(this )))));
0261: }
0262:
0263: //============
0264: // OPERATIONS
0265: //============
0266:
0267: // ADDITION
0268:
0269: /**
0270: * Returns the addition of this matrix and another.
0271: * @param m a double matrix
0272: * @exception MatrixDimensionException If the matrices are different sizes.
0273: */
0274: public AbstractDoubleSquareMatrix add(
0275: final AbstractDoubleSquareMatrix m) {
0276: if (m instanceof DoubleTridiagonalMatrix)
0277: return add((DoubleTridiagonalMatrix) m);
0278: if (m instanceof TridiagonalMatrix)
0279: return addTridiagonal(m);
0280: if (m instanceof DoubleSquareMatrix)
0281: return add((DoubleSquareMatrix) m);
0282:
0283: if (numRows == m.rows() && numCols == m.columns()) {
0284: final double array[][] = new double[numRows][numCols];
0285: for (int i = 0; i < numRows; i++) {
0286: array[i][0] = getElement(i, 0) + m.getElement(i, 0);
0287: for (int j = 1; j < numCols; j++)
0288: array[i][j] = getElement(i, j) + m.getElement(i, j);
0289: }
0290: return new DoubleSquareMatrix(array);
0291: } else {
0292: throw new MatrixDimensionException(
0293: "Matrices are different sizes.");
0294: }
0295: }
0296:
0297: public DoubleSquareMatrix add(final DoubleSquareMatrix m) {
0298: if (numRows == m.numRows && numCols == m.numCols) {
0299: final double array[][] = new double[numRows][numCols];
0300: for (int i = 0; i < numRows; i++)
0301: System.arraycopy(m.matrix[i], 0, array[i], 0, numRows);
0302: array[0][0] += diag[0];
0303: array[0][1] += udiag[0];
0304: int n = numRows - 1;
0305: for (int i = 1; i < n; i++) {
0306: array[i][i - 1] += ldiag[i];
0307: array[i][i] += diag[i];
0308: array[i][i + 1] += udiag[i];
0309: }
0310: array[n][n - 1] += ldiag[n];
0311: array[n][n] += diag[n];
0312: return new DoubleSquareMatrix(array);
0313: } else
0314: throw new MatrixDimensionException(
0315: "Matrices are different sizes.");
0316: }
0317:
0318: /**
0319: * Returns the addition of this matrix and another.
0320: * @param m a double tridiagonal matrix
0321: * @exception MatrixDimensionException If the matrices are different sizes.
0322: */
0323: public DoubleTridiagonalMatrix add(final DoubleTridiagonalMatrix m) {
0324: int mRow = numRows;
0325: if (mRow == m.numRows) {
0326: final DoubleTridiagonalMatrix ans = new DoubleTridiagonalMatrix(
0327: mRow);
0328: ans.diag[0] = diag[0] + m.diag[0];
0329: ans.udiag[0] = udiag[0] + m.udiag[0];
0330: mRow--;
0331: for (int i = 1; i < mRow; i++) {
0332: ans.ldiag[i] = ldiag[i] + m.ldiag[i];
0333: ans.diag[i] = diag[i] + m.diag[i];
0334: ans.udiag[i] = udiag[i] + m.udiag[i];
0335: }
0336: ans.ldiag[mRow] = ldiag[mRow] + m.ldiag[mRow];
0337: ans.diag[mRow] = diag[mRow] + m.diag[mRow];
0338: return ans;
0339: } else
0340: throw new MatrixDimensionException(
0341: "Matrices are different sizes.");
0342: }
0343:
0344: private DoubleTridiagonalMatrix addTridiagonal(
0345: final AbstractDoubleSquareMatrix m) {
0346: int mRow = numRows;
0347: if (mRow == m.rows()) {
0348: final DoubleTridiagonalMatrix ans = new DoubleTridiagonalMatrix(
0349: mRow);
0350: ans.diag[0] = diag[0] + m.getElement(0, 0);
0351: ans.udiag[0] = udiag[0] + m.getElement(0, 1);
0352: mRow--;
0353: for (int i = 1; i < mRow; i++) {
0354: ans.ldiag[i] = ldiag[i] + m.getElement(i, i - 1);
0355: ans.diag[i] = diag[i] + m.getElement(i, i);
0356: ans.udiag[i] = udiag[i] + m.getElement(i, i + 1);
0357: }
0358: ans.ldiag[mRow] = ldiag[mRow]
0359: + m.getElement(mRow, mRow - 1);
0360: ans.diag[mRow] = diag[mRow] + m.getElement(mRow, mRow);
0361: return ans;
0362: } else {
0363: throw new MatrixDimensionException(
0364: "Matrices are different sizes.");
0365: }
0366: }
0367:
0368: // SUBTRACTION
0369:
0370: /**
0371: * Returns the subtraction of this matrix by another.
0372: * @param m a double matrix
0373: * @exception MatrixDimensionException If the matrices are different sizes.
0374: */
0375: public AbstractDoubleSquareMatrix subtract(
0376: final AbstractDoubleSquareMatrix m) {
0377: if (m instanceof DoubleTridiagonalMatrix)
0378: return subtract((DoubleTridiagonalMatrix) m);
0379: if (m instanceof TridiagonalMatrix)
0380: return subtractTridiagonal(m);
0381: if (m instanceof DoubleSquareMatrix)
0382: return subtract((DoubleSquareMatrix) m);
0383:
0384: if (numRows == m.rows() && numCols == m.columns()) {
0385: final double array[][] = new double[numRows][numCols];
0386: for (int i = 0; i < numRows; i++) {
0387: array[i][0] = getElement(i, 0) - m.getElement(i, 0);
0388: for (int j = 1; j < numCols; j++)
0389: array[i][j] = getElement(i, j) - m.getElement(i, j);
0390: }
0391: return new DoubleSquareMatrix(array);
0392: } else {
0393: throw new MatrixDimensionException(
0394: "Matrices are different sizes.");
0395: }
0396: }
0397:
0398: public DoubleSquareMatrix subtract(final DoubleSquareMatrix m) {
0399: if (numRows == m.numRows && numCols == m.numCols) {
0400: final double array[][] = new double[numRows][numCols];
0401: for (int i = 0; i < numRows; i++) {
0402: array[i][0] = getElement(i, 0) - m.matrix[i][0];
0403: for (int j = 1; j < numCols; j++)
0404: array[i][j] = getElement(i, j) - m.matrix[i][j];
0405: }
0406: return new DoubleSquareMatrix(array);
0407: } else
0408: throw new MatrixDimensionException(
0409: "Matrices are different sizes.");
0410: }
0411:
0412: /**
0413: * Returns the subtraction of this matrix and another.
0414: * @param m a double tridiagonal matrix
0415: * @exception MatrixDimensionException If the matrices are different sizes.
0416: */
0417: public DoubleTridiagonalMatrix subtract(
0418: final DoubleTridiagonalMatrix m) {
0419: int mRow = numRows;
0420: if (mRow == m.numRows) {
0421: final DoubleTridiagonalMatrix ans = new DoubleTridiagonalMatrix(
0422: mRow);
0423: ans.diag[0] = diag[0] - m.diag[0];
0424: ans.udiag[0] = udiag[0] - m.udiag[0];
0425: mRow--;
0426: for (int i = 1; i < mRow; i++) {
0427: ans.ldiag[i] = ldiag[i] - m.ldiag[i];
0428: ans.diag[i] = diag[i] - m.diag[i];
0429: ans.udiag[i] = udiag[i] - m.udiag[i];
0430: }
0431: ans.ldiag[mRow] = ldiag[mRow] - m.ldiag[mRow];
0432: ans.diag[mRow] = diag[mRow] - m.diag[mRow];
0433: return ans;
0434: } else
0435: throw new MatrixDimensionException(
0436: "Matrices are different sizes.");
0437: }
0438:
0439: private DoubleTridiagonalMatrix subtractTridiagonal(
0440: final AbstractDoubleSquareMatrix m) {
0441: int mRow = numRows;
0442: if (mRow == m.rows()) {
0443: final DoubleTridiagonalMatrix ans = new DoubleTridiagonalMatrix(
0444: mRow);
0445: ans.diag[0] = diag[0] - m.getElement(0, 0);
0446: ans.udiag[0] = udiag[0] - m.getElement(0, 1);
0447: mRow--;
0448: for (int i = 1; i < mRow; i++) {
0449: ans.ldiag[i] = ldiag[i] - m.getElement(i, i - 1);
0450: ans.diag[i] = diag[i] - m.getElement(i, i);
0451: ans.udiag[i] = udiag[i] - m.getElement(i, i + 1);
0452: }
0453: ans.ldiag[mRow] = ldiag[mRow]
0454: - m.getElement(mRow, mRow - 1);
0455: ans.diag[mRow] = diag[mRow] - m.getElement(mRow, mRow);
0456: return ans;
0457: } else {
0458: throw new MatrixDimensionException(
0459: "Matrices are different sizes.");
0460: }
0461: }
0462:
0463: // SCALAR MULTIPLICATION
0464:
0465: /**
0466: * Returns the multiplication of this matrix by a scalar.
0467: * @param x a double.
0468: * @return a double square matrix.
0469: */
0470: public AbstractDoubleMatrix scalarMultiply(final double x) {
0471: final DoubleTridiagonalMatrix ans = new DoubleTridiagonalMatrix(
0472: numRows);
0473: final int lastRow = numRows - 1;
0474: ans.diag[0] = x * diag[0];
0475: ans.udiag[0] = x * udiag[0];
0476: for (int i = 1; i < lastRow; i++) {
0477: ans.ldiag[i] = x * ldiag[i];
0478: ans.diag[i] = x * diag[i];
0479: ans.udiag[i] = x * udiag[i];
0480: }
0481: ans.ldiag[lastRow] = x * ldiag[lastRow];
0482: ans.diag[lastRow] = x * diag[lastRow];
0483: return ans;
0484: }
0485:
0486: // SCALAR DIVISON
0487:
0488: /**
0489: * Returns the division of this matrix by a scalar.
0490: * @param x a double.
0491: * @return a double square matrix.
0492: */
0493: public AbstractDoubleMatrix scalarDivide(final double x) {
0494: final DoubleTridiagonalMatrix ans = new DoubleTridiagonalMatrix(
0495: numRows);
0496: final int lastRow = numRows - 1;
0497: ans.diag[0] = diag[0] / x;
0498: ans.udiag[0] = udiag[0] / x;
0499: for (int i = 1; i < lastRow; i++) {
0500: ans.ldiag[i] = ldiag[i] / x;
0501: ans.diag[i] = diag[i] / x;
0502: ans.udiag[i] = udiag[i] / x;
0503: }
0504: ans.ldiag[lastRow] = ldiag[lastRow] / x;
0505: ans.diag[lastRow] = diag[lastRow] / x;
0506: return ans;
0507: }
0508:
0509: // SCALAR PRODUCT
0510:
0511: /**
0512: * Returns the scalar product of this matrix and another.
0513: * @param m a double matrix.
0514: * @exception MatrixDimensionException If the matrices are different sizes.
0515: */
0516: public double scalarProduct(final AbstractDoubleSquareMatrix m) {
0517: if (m instanceof DoubleTridiagonalMatrix)
0518: return scalarProduct((DoubleTridiagonalMatrix) m);
0519: if (m instanceof DoubleSquareMatrix)
0520: return scalarProduct((DoubleSquareMatrix) m);
0521:
0522: if (numRows == m.rows() && numCols == m.columns()) {
0523: double ans = diag[0] * m.getElement(0, 0) + udiag[0]
0524: * m.getElement(0, 1);
0525: final int lastRow = numRows - 1;
0526: for (int i = 1; i < lastRow; i++)
0527: ans += ldiag[i] * m.getElement(i, i - 1) + diag[i]
0528: * m.getElement(i, i) + udiag[i]
0529: * m.getElement(i, i + 1);
0530: ans += ldiag[lastRow] * m.getElement(lastRow, lastRow - 1)
0531: + diag[lastRow] * m.getElement(lastRow, lastRow);
0532: return ans;
0533: } else {
0534: throw new MatrixDimensionException(
0535: "Matrices are different sizes.");
0536: }
0537: }
0538:
0539: public double scalarProduct(final DoubleSquareMatrix m) {
0540: if (numRows == m.numRows && numCols == m.numCols) {
0541: double ans = diag[0] * m.matrix[0][0] + udiag[0]
0542: * m.matrix[0][1];
0543: final int lastRow = numRows - 1;
0544: for (int i = 1; i < lastRow; i++)
0545: ans += ldiag[i] * m.matrix[i][i - 1] + diag[i]
0546: * m.matrix[i][i] + udiag[i]
0547: * m.matrix[i][i + 1];
0548: ans += ldiag[lastRow] * m.matrix[lastRow][lastRow - 1]
0549: + diag[lastRow] * m.matrix[lastRow][lastRow];
0550: return ans;
0551: } else
0552: throw new MatrixDimensionException(
0553: "Matrices are different sizes.");
0554: }
0555:
0556: public double scalarProduct(final DoubleTridiagonalMatrix m) {
0557: if (numRows == m.numRows) {
0558: double ans = diag[0] * m.diag[0] + udiag[0] * m.udiag[0];
0559: final int lastRow = numRows - 1;
0560: for (int i = 1; i < lastRow; i++)
0561: ans += ldiag[i] * m.ldiag[i] + diag[i] * m.diag[i]
0562: + udiag[i] * m.udiag[i];
0563: ans += ldiag[lastRow] * m.ldiag[lastRow] + diag[lastRow]
0564: * m.diag[lastRow];
0565: return ans;
0566: } else
0567: throw new MatrixDimensionException(
0568: "Matrices are different sizes.");
0569: }
0570:
0571: // MATRIX MULTIPLICATION
0572:
0573: /**
0574: * Returns the multiplication of a vector by this matrix.
0575: * @param v a double vector.
0576: * @exception DimensionException If the matrix and vector are incompatible.
0577: */
0578: public AbstractDoubleVector multiply(final AbstractDoubleVector v) {
0579: if (numCols == v.dimension()) {
0580: final double array[] = new double[numRows];
0581: final int lastRow = numRows - 1;
0582: array[0] = diag[0] * v.getComponent(0) + udiag[0]
0583: * v.getComponent(1);
0584: for (int i = 1; i < lastRow; i++)
0585: array[i] = ldiag[i] * v.getComponent(i - 1) + diag[i]
0586: * v.getComponent(i) + udiag[i]
0587: * v.getComponent(i + 1);
0588: array[lastRow] = ldiag[lastRow]
0589: * v.getComponent(lastRow - 1) + diag[lastRow]
0590: * v.getComponent(lastRow);
0591: return new DoubleVector(array);
0592: } else {
0593: throw new DimensionException(
0594: "Matrix and vector are incompatible.");
0595: }
0596: }
0597:
0598: /**
0599: * Returns the multiplication of this matrix and another.
0600: * @param m a double matrix
0601: * @return a AbstractDoubleMatrix or a AbstractDoubleSquareMatrix as appropriate
0602: * @exception MatrixDimensionException If the matrices are incompatible.
0603: */
0604: public AbstractDoubleSquareMatrix multiply(
0605: final AbstractDoubleSquareMatrix m) {
0606: if (m instanceof DoubleTridiagonalMatrix)
0607: return multiply((DoubleTridiagonalMatrix) m);
0608: if (m instanceof TridiagonalMatrix)
0609: return multiplyTridiagonal(m);
0610: if (m instanceof DoubleSquareMatrix)
0611: return multiply((DoubleSquareMatrix) m);
0612:
0613: if (numCols == m.rows()) {
0614: final int mColumns = m.columns();
0615: final double array[][] = new double[numRows][mColumns];
0616: final int lastRow = numRows - 1;
0617: for (int j = 0; j < numRows; j++) {
0618: array[0][j] = diag[0] * m.getElement(0, j) + udiag[0]
0619: * m.getElement(1, j);
0620: for (int i = 1; i < lastRow; i++)
0621: array[i][j] = ldiag[i] * m.getElement(i - 1, j)
0622: + diag[i] * m.getElement(i, j) + udiag[i]
0623: * m.getElement(i + 1, j);
0624: array[lastRow][j] = ldiag[lastRow]
0625: * m.getElement(lastRow - 1, j) + diag[lastRow]
0626: * m.getElement(lastRow, j);
0627: }
0628: return new DoubleSquareMatrix(array);
0629: } else {
0630: throw new MatrixDimensionException("Incompatible matrices.");
0631: }
0632: }
0633:
0634: public DoubleSquareMatrix multiply(final DoubleSquareMatrix m) {
0635: if (numCols == m.numRows) {
0636: final double array[][] = new double[numRows][m.numCols];
0637: final int lastRow = numRows - 1;
0638: for (int j = 0; j < numRows; j++) {
0639: array[0][j] = diag[0] * m.matrix[0][j] + udiag[0]
0640: * m.matrix[1][j];
0641: for (int i = 1; i < lastRow; i++)
0642: array[i][j] = ldiag[i] * m.matrix[i - 1][j]
0643: + diag[i] * m.matrix[i][j] + udiag[i]
0644: * m.matrix[i + 1][j];
0645: array[lastRow][j] = ldiag[lastRow]
0646: * m.matrix[lastRow - 1][j] + diag[lastRow]
0647: * m.matrix[lastRow][j];
0648: }
0649: return new DoubleSquareMatrix(array);
0650: } else
0651: throw new MatrixDimensionException("Incompatible matrices.");
0652: }
0653:
0654: public DoubleSquareMatrix multiply(final DoubleTridiagonalMatrix m) {
0655: int mRow = numRows;
0656: if (numCols == m.numRows) {
0657: final double array[][] = new double[mRow][mRow];
0658: array[0][0] = diag[0] * m.diag[0] + udiag[0] * m.ldiag[1];
0659: array[0][1] = diag[0] * m.udiag[0] + udiag[0] * m.diag[1];
0660: array[0][2] = udiag[0] * m.udiag[1];
0661: if (mRow > 3) {
0662: array[1][0] = ldiag[1] * m.diag[0] + diag[1]
0663: * m.ldiag[1];
0664: array[1][1] = ldiag[1] * m.udiag[0] + diag[1]
0665: * m.diag[1] + udiag[1] * m.ldiag[2];
0666: array[1][2] = diag[1] * m.udiag[1] + udiag[1]
0667: * m.diag[2];
0668: array[1][3] = udiag[1] * m.udiag[2];
0669: }
0670: if (mRow == 3) {
0671: array[1][0] = ldiag[1] * m.diag[0] + diag[1]
0672: * m.ldiag[1];
0673: array[1][1] = ldiag[1] * m.udiag[0] + diag[1]
0674: * m.diag[1] + udiag[1] * m.ldiag[2];
0675: array[1][2] = diag[1] * m.udiag[1] + udiag[1]
0676: * m.diag[2];
0677: } else if (mRow > 4) {
0678: for (int i = 2; i < mRow - 2; i++) {
0679: array[i][i - 2] = ldiag[i] * m.ldiag[i - 1];
0680: array[i][i - 1] = ldiag[i] * m.diag[i - 1]
0681: + diag[i] * m.ldiag[i];
0682: array[i][i] = ldiag[i] * m.udiag[i - 1] + diag[i]
0683: * m.diag[i] + udiag[i] * m.ldiag[i + 1];
0684: array[i][i + 1] = diag[i] * m.udiag[i] + udiag[i]
0685: * m.diag[i + 1];
0686: array[i][i + 2] = udiag[i] * m.udiag[i + 1];
0687: }
0688: }
0689: if (mRow > 3) {
0690: array[mRow - 2][mRow - 4] = ldiag[mRow - 2]
0691: * m.ldiag[mRow - 3];
0692: array[mRow - 2][mRow - 3] = ldiag[mRow - 2]
0693: * m.diag[mRow - 3] + diag[mRow - 2]
0694: * m.ldiag[mRow - 2];
0695: array[mRow - 2][mRow - 2] = ldiag[mRow - 2]
0696: * m.udiag[mRow - 3] + diag[mRow - 2]
0697: * m.diag[mRow - 2] + udiag[mRow - 2]
0698: * m.ldiag[mRow - 1];
0699: array[mRow - 2][mRow - 1] = diag[mRow - 2]
0700: * m.udiag[mRow - 2] + udiag[mRow - 2]
0701: * m.diag[mRow - 1];
0702: }
0703: mRow--;
0704: array[mRow][mRow - 2] = ldiag[mRow] * m.ldiag[mRow - 1];
0705: array[mRow][mRow - 1] = ldiag[mRow] * m.diag[mRow - 1]
0706: + diag[mRow] * m.ldiag[mRow];
0707: array[mRow][mRow] = ldiag[mRow] * m.udiag[mRow - 1]
0708: + diag[mRow] * m.diag[mRow];
0709: return new DoubleSquareMatrix(array);
0710: } else
0711: throw new MatrixDimensionException("Incompatible matrices.");
0712: }
0713:
0714: private DoubleSquareMatrix multiplyTridiagonal(
0715: final AbstractDoubleSquareMatrix m) {
0716: int mRow = numRows;
0717: if (numCols == m.rows()) {
0718: final double array[][] = new double[mRow][mRow];
0719: array[0][0] = diag[0] * m.getElement(0, 0) + udiag[0]
0720: * m.getElement(1, 0);
0721: array[0][1] = diag[0] * m.getElement(0, 1) + udiag[0]
0722: * m.getElement(1, 1);
0723: array[0][2] = udiag[0] * m.getElement(1, 2);
0724: if (mRow > 3) {
0725: array[1][0] = ldiag[1] * m.getElement(0, 0) + diag[1]
0726: * m.getElement(1, 0);
0727: array[1][1] = ldiag[1] * m.getElement(0, 1) + diag[1]
0728: * m.getElement(1, 1) + udiag[1]
0729: * m.getElement(2, 1);
0730: array[1][2] = diag[1] * m.getElement(1, 2) + udiag[1]
0731: * m.getElement(2, 2);
0732: array[1][3] = udiag[1] * m.getElement(2, 3);
0733: }
0734: if (mRow == 3) {
0735: array[1][0] = ldiag[1] * m.getElement(0, 0) + diag[1]
0736: * m.getElement(1, 0);
0737: array[1][1] = ldiag[1] * m.getElement(0, 1) + diag[1]
0738: * m.getElement(1, 1) + udiag[1]
0739: * m.getElement(2, 1);
0740: array[1][2] = diag[1] * m.getElement(1, 2) + udiag[1]
0741: * m.getElement(2, 2);
0742: } else if (mRow > 4) {
0743: for (int i = 2; i < mRow - 2; i++) {
0744: array[i][i - 2] = ldiag[i]
0745: * m.getElement(i - 1, i - 2);
0746: array[i][i - 1] = ldiag[i]
0747: * m.getElement(i - 1, i - 1) + diag[i]
0748: * m.getElement(i, i - 1);
0749: array[i][i] = ldiag[i] * m.getElement(i - 1, i)
0750: + diag[i] * m.getElement(i, i) + udiag[i]
0751: * m.getElement(i + 1, i);
0752: array[i][i + 1] = diag[i] * m.getElement(i, i + 1)
0753: + udiag[i] * m.getElement(i + 1, i + 1);
0754: array[i][i + 2] = udiag[i]
0755: * m.getElement(i + 1, i + 2);
0756: }
0757: }
0758: if (mRow > 3) {
0759: array[mRow - 2][mRow - 4] = ldiag[mRow - 2]
0760: * m.getElement(mRow - 3, mRow - 4);
0761: array[mRow - 2][mRow - 3] = ldiag[mRow - 2]
0762: * m.getElement(mRow - 3, mRow - 3)
0763: + diag[mRow - 2]
0764: * m.getElement(mRow - 2, mRow - 3);
0765: array[mRow - 2][mRow - 2] = ldiag[mRow - 2]
0766: * m.getElement(mRow - 3, mRow - 2)
0767: + diag[mRow - 2]
0768: * m.getElement(mRow - 2, mRow - 2)
0769: + udiag[mRow - 2]
0770: * m.getElement(mRow - 1, mRow - 2);
0771: array[mRow - 2][mRow - 1] = diag[mRow - 2]
0772: * m.getElement(mRow - 2, mRow - 1)
0773: + udiag[mRow - 2]
0774: * m.getElement(mRow - 1, mRow - 1);
0775: }
0776: mRow--;
0777: array[mRow][mRow - 2] = ldiag[mRow]
0778: * m.getElement(mRow - 1, mRow - 2);
0779: array[mRow][mRow - 1] = ldiag[mRow]
0780: * m.getElement(mRow - 1, mRow - 1) + diag[mRow]
0781: * m.getElement(mRow, mRow - 1);
0782: array[mRow][mRow] = ldiag[mRow]
0783: * m.getElement(mRow - 1, mRow) + diag[mRow]
0784: * m.getElement(mRow, mRow);
0785: return new DoubleSquareMatrix(array);
0786: } else {
0787: throw new MatrixDimensionException("Incompatible matrices.");
0788: }
0789: }
0790:
0791: // TRANSPOSE
0792:
0793: /**
0794: * Returns the transpose of this matrix.
0795: * @return a double matrix
0796: */
0797: public Matrix transpose() {
0798: final DoubleTridiagonalMatrix ans = new DoubleTridiagonalMatrix(
0799: numRows);
0800: System.arraycopy(ldiag, 1, ans.udiag, 0, ldiag.length - 1);
0801: System.arraycopy(diag, 0, ans.diag, 0, diag.length);
0802: System.arraycopy(udiag, 0, ans.ldiag, 1, udiag.length - 1);
0803: return ans;
0804: }
0805:
0806: // CHOLESKY DECOMPOSITION
0807:
0808: /**
0809: * Returns the Cholesky decomposition of this matrix.
0810: * Matrix must be symmetric and positive definite.
0811: * @return an array with [0] containing the L-matrix and [1] containing the U-matrix.
0812: */
0813: public AbstractDoubleSquareMatrix[] choleskyDecompose() {
0814: final int N = numRows;
0815: final double arrayL[][] = new double[N][N];
0816: final double arrayU[][] = new double[N][N];
0817: double tmp = Math.sqrt(diag[0]);
0818: arrayL[0][0] = arrayU[0][0] = tmp;
0819: arrayL[1][0] = arrayU[0][1] = ldiag[1] / tmp;
0820: for (int j = 1; j < N; j++) {
0821: tmp = diag[j];
0822: for (int i = 0; i < j; i++)
0823: tmp -= arrayL[j][i] * arrayL[j][i];
0824: arrayL[j][j] = arrayU[j][j] = Math.sqrt(tmp);
0825: if (j + 1 < N) {
0826: tmp = ldiag[j + 1];
0827: for (int k = 0; k < j + 1; k++)
0828: tmp -= arrayL[j][k] * arrayU[k][j + 1];
0829: arrayL[j + 1][j] = arrayU[j][j + 1] = tmp
0830: / arrayU[j][j];
0831: for (int i = j + 2; i < N; i++) {
0832: tmp = 0.0;
0833: for (int k = 0; k < i; k++)
0834: tmp -= arrayL[j][k] * arrayU[k][i];
0835: arrayL[i][j] = arrayU[j][i] = tmp / arrayU[j][j];
0836: }
0837: }
0838: }
0839: final AbstractDoubleSquareMatrix lu[] = new AbstractDoubleSquareMatrix[2];
0840: lu[0] = new DoubleSquareMatrix(arrayL);
0841: lu[1] = new DoubleSquareMatrix(arrayU);
0842: return lu;
0843: }
0844:
0845: // QR DECOMPOSITION
0846:
0847: /**
0848: * Returns the QR decomposition of this matrix.
0849: * Based on the code from <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> (public domain).
0850: * @return an array with [0] containing the Q-matrix and [1] containing the R-matrix.
0851: * @jsci.planetmath QRDecomposition
0852: */
0853: public AbstractDoubleSquareMatrix[] qrDecompose() {
0854: final int N = numRows;
0855: final double array[][] = new double[N][N];
0856: final double arrayQ[][] = new double[N][N];
0857: final double arrayR[][] = new double[N][N];
0858: // copy matrix
0859: array[0][0] = diag[0];
0860: array[0][1] = udiag[0];
0861: for (int i = 1; i < N - 1; i++) {
0862: array[i][i - 1] = ldiag[i];
0863: array[i][i] = diag[i];
0864: array[i][i + 1] = udiag[i];
0865: }
0866: array[N - 1][N - 2] = ldiag[N - 1];
0867: array[N - 1][N - 1] = diag[N - 1];
0868: for (int k = 0; k < N; k++) {
0869: // compute l2-norm of kth column
0870: double norm = array[k][k];
0871: for (int i = k + 1; i < N; i++)
0872: norm = ExtraMath.hypot(norm, array[i][k]);
0873: if (norm != 0.0) {
0874: // form kth Householder vector
0875: if (array[k][k] < 0.0)
0876: norm = -norm;
0877: for (int i = k; i < N; i++)
0878: array[i][k] /= norm;
0879: array[k][k] += 1.0;
0880: // apply transformation to remaining columns
0881: for (int j = k + 1; j < N; j++) {
0882: double s = array[k][k] * array[k][j];
0883: for (int i = k + 1; i < N; i++)
0884: s += array[i][k] * array[i][j];
0885: s /= array[k][k];
0886: for (int i = k; i < N; i++)
0887: array[i][j] -= s * array[i][k];
0888: }
0889: }
0890: arrayR[k][k] = -norm;
0891: }
0892: for (int k = N - 1; k >= 0; k--) {
0893: arrayQ[k][k] = 1.0;
0894: for (int j = k; j < N; j++) {
0895: if (array[k][k] != 0.0) {
0896: double s = array[k][k] * arrayQ[k][j];
0897: for (int i = k + 1; i < N; i++)
0898: s += array[i][k] * arrayQ[i][j];
0899: s /= array[k][k];
0900: for (int i = k; i < N; i++)
0901: arrayQ[i][j] -= s * array[i][k];
0902: }
0903: }
0904: }
0905: for (int i = 0; i < N; i++) {
0906: for (int j = i + 1; j < N; j++)
0907: arrayR[i][j] = array[i][j];
0908: }
0909: final AbstractDoubleSquareMatrix qr[] = new AbstractDoubleSquareMatrix[2];
0910: qr[0] = new DoubleSquareMatrix(arrayQ);
0911: qr[1] = new DoubleSquareMatrix(arrayR);
0912: return qr;
0913: }
0914:
0915: // SINGULAR VALUE DECOMPOSITION
0916:
0917: /**
0918: * Returns the singular value decomposition of this matrix.
0919: * Based on the code from <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> (public domain).
0920: * @return an array with [0] containing the U-matrix, [1] containing the S-matrix and [2] containing the V-matrix.
0921: */
0922: public AbstractDoubleSquareMatrix[] singularValueDecompose() {
0923: final int N = numRows;
0924: final int Nm1 = N - 1;
0925: final double array[][] = new double[N][N];
0926: final double arrayU[][] = new double[N][N];
0927: final double arrayS[] = new double[N];
0928: final double arrayV[][] = new double[N][N];
0929: final double e[] = new double[N];
0930: final double work[] = new double[N];
0931: // copy matrix
0932: array[0][0] = diag[0];
0933: array[0][1] = udiag[0];
0934: for (int i = 1; i < Nm1; i++) {
0935: array[i][i - 1] = ldiag[i];
0936: array[i][i] = diag[i];
0937: array[i][i + 1] = udiag[i];
0938: }
0939: array[Nm1][Nm1 - 1] = ldiag[Nm1];
0940: array[Nm1][Nm1] = diag[Nm1];
0941: // reduce matrix to bidiagonal form
0942: for (int k = 0; k < Nm1; k++) {
0943: // compute the transformation for the kth column
0944: // compute l2-norm of kth column
0945: arrayS[k] = array[k][k];
0946: for (int i = k + 1; i < N; i++)
0947: arrayS[k] = ExtraMath.hypot(arrayS[k], array[i][k]);
0948: if (arrayS[k] != 0.0) {
0949: if (array[k][k] < 0.0)
0950: arrayS[k] = -arrayS[k];
0951: for (int i = k; i < N; i++)
0952: array[i][k] /= arrayS[k];
0953: array[k][k] += 1.0;
0954: }
0955: arrayS[k] = -arrayS[k];
0956: for (int j = k + 1; j < N; j++) {
0957: if (arrayS[k] != 0.0) {
0958: // apply the transformation
0959: double t = array[k][k] * array[k][j];
0960: for (int i = k + 1; i < N; i++)
0961: t += array[i][k] * array[i][j];
0962: t /= array[k][k];
0963: for (int i = k; i < N; i++)
0964: array[i][j] -= t * array[i][k];
0965: }
0966: e[j] = array[k][j];
0967: }
0968: for (int i = k; i < N; i++)
0969: arrayU[i][k] = array[i][k];
0970: if (k < N - 2) {
0971: // compute the kth row transformation
0972: // compute l2-norm of kth column
0973: e[k] = e[k + 1];
0974: for (int i = k + 2; i < N; i++)
0975: e[k] = ExtraMath.hypot(e[k], e[i]);
0976: if (e[k] != 0.0) {
0977: if (e[k + 1] < 0.0)
0978: e[k] = -e[k];
0979: for (int i = k + 1; i < N; i++)
0980: e[i] /= e[k];
0981: e[k + 1] += 1.0;
0982: }
0983: e[k] = -e[k];
0984: if (e[k] != 0.0) {
0985: // apply the transformation
0986: for (int i = k + 1; i < N; i++) {
0987: work[i] = 0.0;
0988: for (int j = k + 1; j < N; j++)
0989: work[i] += e[j] * array[i][j];
0990: }
0991: for (int j = k + 1; j < N; j++) {
0992: double t = e[j] / e[k + 1];
0993: for (int i = k + 1; i < N; i++)
0994: array[i][j] -= t * work[i];
0995: }
0996: }
0997: for (int i = k + 1; i < N; i++)
0998: arrayV[i][k] = e[i];
0999: }
1000: }
1001: // setup the final bidiagonal matrix of order p
1002: int p = N;
1003: arrayS[Nm1] = array[Nm1][Nm1];
1004: e[N - 2] = array[N - 2][Nm1];
1005: e[Nm1] = 0.0;
1006: arrayU[Nm1][Nm1] = 1.0;
1007: for (int k = N - 2; k >= 0; k--) {
1008: if (arrayS[k] != 0.0) {
1009: for (int j = k + 1; j < N; j++) {
1010: double t = arrayU[k][k] * arrayU[k][j];
1011: for (int i = k + 1; i < N; i++)
1012: t += arrayU[i][k] * arrayU[i][j];
1013: t /= arrayU[k][k];
1014: for (int i = k; i < N; i++)
1015: arrayU[i][j] -= t * arrayU[i][k];
1016: }
1017: for (int i = k; i < N; i++)
1018: arrayU[i][k] = -arrayU[i][k];
1019: arrayU[k][k] += 1.0;
1020: for (int i = 0; i < k - 1; i++)
1021: arrayU[i][k] = 0.0;
1022: } else {
1023: for (int i = 0; i < N; i++)
1024: arrayU[i][k] = 0.0;
1025: arrayU[k][k] = 1.0;
1026: }
1027: }
1028: for (int k = Nm1; k >= 0; k--) {
1029: if (k < N - 2 && e[k] != 0.0) {
1030: for (int j = k + 1; j < N; j++) {
1031: double t = arrayV[k + 1][k] * arrayV[k + 1][j];
1032: for (int i = k + 2; i < N; i++)
1033: t += arrayV[i][k] * arrayV[i][j];
1034: t /= arrayV[k + 1][k];
1035: for (int i = k + 1; i < N; i++)
1036: arrayV[i][j] -= t * arrayV[i][k];
1037: }
1038: }
1039: for (int i = 0; i < N; i++)
1040: arrayV[i][k] = 0.0;
1041: arrayV[k][k] = 1.0;
1042: }
1043: final double eps = Math.pow(2.0, -52.0);
1044: int iter = 0;
1045: while (p > 0) {
1046: int k, action;
1047: // action = 1 if arrayS[p] and e[k-1] are negligible and k<p
1048: // action = 2 if arrayS[k] is negligible and k<p
1049: // action = 3 if e[k-1] is negligible, k<p, and arrayS[k], ..., arrayS[p] are not negligible (QR step)
1050: // action = 4 if e[p-1] is negligible (convergence)
1051: for (k = p - 2; k >= -1; k--) {
1052: if (k == -1)
1053: break;
1054: if (Math.abs(e[k]) <= eps
1055: * (Math.abs(arrayS[k]) + Math
1056: .abs(arrayS[k + 1]))) {
1057: e[k] = 0.0;
1058: break;
1059: }
1060: }
1061: if (k == p - 2) {
1062: action = 4;
1063: } else {
1064: int ks;
1065: for (ks = p - 1; ks >= k; ks--) {
1066: if (ks == k)
1067: break;
1068: double t = (ks != p ? Math.abs(e[ks]) : 0.0)
1069: + (ks != k + 1 ? Math.abs(e[ks - 1]) : 0.0);
1070: if (Math.abs(arrayS[ks]) <= eps * t) {
1071: arrayS[ks] = 0.0;
1072: break;
1073: }
1074: }
1075: if (ks == k) {
1076: action = 3;
1077: } else if (ks == p - 1) {
1078: action = 1;
1079: } else {
1080: action = 2;
1081: k = ks;
1082: }
1083: }
1084: k++;
1085: switch (action) {
1086: // deflate negligible arrayS[p]
1087: case 1: {
1088: double f = e[p - 2];
1089: e[p - 2] = 0.0;
1090: for (int j = p - 2; j >= k; j--) {
1091: double t = ExtraMath.hypot(arrayS[j], f);
1092: final double cs = arrayS[j] / t;
1093: final double sn = f / t;
1094: arrayS[j] = t;
1095: if (j != k) {
1096: f = -sn * e[j - 1];
1097: e[j - 1] *= cs;
1098: }
1099: for (int i = 0; i < N; i++) {
1100: t = cs * arrayV[i][j] + sn * arrayV[i][p - 1];
1101: arrayV[i][p - 1] = -sn * arrayV[i][j] + cs
1102: * arrayV[i][p - 1];
1103: arrayV[i][j] = t;
1104: }
1105: }
1106: }
1107: break;
1108: // split at negligible arrayS[k]
1109: case 2: {
1110: double f = e[k - 1];
1111: e[k - 1] = 0.0;
1112: for (int j = k; j < p; j++) {
1113: double t = ExtraMath.hypot(arrayS[j], f);
1114: final double cs = arrayS[j] / t;
1115: final double sn = f / t;
1116: arrayS[j] = t;
1117: f = -sn * e[j];
1118: e[j] *= cs;
1119: for (int i = 0; i < N; i++) {
1120: t = cs * arrayU[i][j] + sn * arrayU[i][k - 1];
1121: arrayU[i][k - 1] = -sn * arrayU[i][j] + cs
1122: * arrayU[i][k - 1];
1123: arrayU[i][j] = t;
1124: }
1125: }
1126: }
1127: break;
1128: // perform one QR step
1129: case 3: {
1130: // calculate the shift
1131: final double scale = Math.max(Math.max(Math.max(Math
1132: .max(Math.abs(arrayS[p - 1]), Math
1133: .abs(arrayS[p - 2])), Math
1134: .abs(e[p - 2])), Math.abs(arrayS[k])), Math
1135: .abs(e[k]));
1136: double sp = arrayS[p - 1] / scale;
1137: double spm1 = arrayS[p - 2] / scale;
1138: double epm1 = e[p - 2] / scale;
1139: double sk = arrayS[k] / scale;
1140: double ek = e[k] / scale;
1141: double b = ((spm1 + sp) * (spm1 - sp) + epm1 * epm1) / 2.0;
1142: double c = (sp * epm1) * (sp * epm1);
1143: double shift = 0.0;
1144: if (b != 0.0 || c != 0.0) {
1145: shift = Math.sqrt(b * b + c);
1146: if (b < 0.0)
1147: shift = -shift;
1148: shift = c / (b + shift);
1149: }
1150: double f = (sk + sp) * (sk - sp) + shift;
1151: double g = sk * ek;
1152: // chase zeros
1153: for (int j = k; j < p - 1; j++) {
1154: double t = ExtraMath.hypot(f, g);
1155: double cs = f / t;
1156: double sn = g / t;
1157: if (j != k)
1158: e[j - 1] = t;
1159: f = cs * arrayS[j] + sn * e[j];
1160: e[j] = cs * e[j] - sn * arrayS[j];
1161: g = sn * arrayS[j + 1];
1162: arrayS[j + 1] *= cs;
1163: for (int i = 0; i < N; i++) {
1164: t = cs * arrayV[i][j] + sn * arrayV[i][j + 1];
1165: arrayV[i][j + 1] = -sn * arrayV[i][j] + cs
1166: * arrayV[i][j + 1];
1167: arrayV[i][j] = t;
1168: }
1169: t = ExtraMath.hypot(f, g);
1170: cs = f / t;
1171: sn = g / t;
1172: arrayS[j] = t;
1173: f = cs * e[j] + sn * arrayS[j + 1];
1174: arrayS[j + 1] = -sn * e[j] + cs * arrayS[j + 1];
1175: g = sn * e[j + 1];
1176: e[j + 1] *= cs;
1177: if (j < Nm1) {
1178: for (int i = 0; i < N; i++) {
1179: t = cs * arrayU[i][j] + sn
1180: * arrayU[i][j + 1];
1181: arrayU[i][j + 1] = -sn * arrayU[i][j] + cs
1182: * arrayU[i][j + 1];
1183: arrayU[i][j] = t;
1184: }
1185: }
1186: }
1187: e[p - 2] = f;
1188: iter++;
1189: }
1190: break;
1191: // convergence
1192: case 4: {
1193: // make the singular values positive
1194: if (arrayS[k] <= 0.0) {
1195: arrayS[k] = -arrayS[k];
1196: for (int i = 0; i < p; i++)
1197: arrayV[i][k] = -arrayV[i][k];
1198: }
1199: // order the singular values
1200: while (k < p - 1) {
1201: if (arrayS[k] >= arrayS[k + 1])
1202: break;
1203: double tmp = arrayS[k];
1204: arrayS[k] = arrayS[k + 1];
1205: arrayS[k + 1] = tmp;
1206: if (k < Nm1) {
1207: for (int i = 0; i < N; i++) {
1208: tmp = arrayU[i][k + 1];
1209: arrayU[i][k + 1] = arrayU[i][k];
1210: arrayU[i][k] = tmp;
1211: tmp = arrayV[i][k + 1];
1212: arrayV[i][k + 1] = arrayV[i][k];
1213: arrayV[i][k] = tmp;
1214: }
1215: }
1216: k++;
1217: }
1218: iter = 0;
1219: p--;
1220: }
1221: break;
1222: }
1223: }
1224: final AbstractDoubleSquareMatrix svd[] = new AbstractDoubleSquareMatrix[3];
1225: svd[0] = new DoubleSquareMatrix(arrayU);
1226: svd[1] = new DoubleDiagonalMatrix(arrayS);
1227: svd[2] = new DoubleSquareMatrix(arrayV);
1228: return svd;
1229: }
1230:
1231: // MAP ELEMENTS
1232:
1233: /**
1234: * Applies a function on all the matrix elements.
1235: * @param f a user-defined function
1236: * @return a double matrix
1237: */
1238: public AbstractDoubleMatrix mapElements(final Mapping f) {
1239: double zeroValue = f.map(0.0);
1240: if (Math.abs(zeroValue) <= JSci.GlobalSettings.ZERO_TOL)
1241: return tridiagonalMap(f);
1242: else
1243: return generalMap(f, zeroValue);
1244: }
1245:
1246: private AbstractDoubleMatrix tridiagonalMap(Mapping f) {
1247: int mRow = numRows;
1248: final DoubleTridiagonalMatrix ans = new DoubleTridiagonalMatrix(
1249: mRow);
1250: ans.diag[0] = f.map(diag[0]);
1251: ans.udiag[0] = f.map(udiag[0]);
1252: mRow--;
1253: for (int i = 1; i < mRow; i++) {
1254: ans.ldiag[i] = f.map(ldiag[i]);
1255: ans.diag[i] = f.map(diag[i]);
1256: ans.udiag[i] = f.map(udiag[i]);
1257: }
1258: ans.ldiag[mRow] = f.map(ldiag[mRow]);
1259: ans.diag[mRow] = f.map(diag[mRow]);
1260: return ans;
1261: }
1262:
1263: private AbstractDoubleMatrix generalMap(Mapping f, double zeroValue) {
1264: final double array[][] = new double[numRows][numRows];
1265: for (int i = 0; i < numRows; i++) {
1266: for (int j = 0; j < numRows; j++) {
1267: array[i][j] = zeroValue;
1268: }
1269: }
1270: int mRow = numRows;
1271: array[0][0] = f.map(diag[0]);
1272: array[0][1] = f.map(udiag[0]);
1273: mRow--;
1274: for (int i = 1; i < mRow; i++) {
1275: array[i][i - 1] = f.map(ldiag[i]);
1276: array[i][i] = f.map(diag[i]);
1277: array[i][i + 1] = f.map(udiag[i]);
1278: }
1279: array[mRow][mRow - 1] = f.map(ldiag[mRow]);
1280: array[mRow][mRow] = f.map(diag[mRow]);
1281: return new DoubleSquareMatrix(array);
1282: }
1283: }
|