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.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 IntegerTridiagonalMatrix class provides an object for encapsulating integer tridiagonal matrices.
0018: * @version 2.2
0019: * @author Mark Hale
0020: */
0021: public class IntegerTridiagonalMatrix extends
0022: AbstractIntegerSquareMatrix implements TridiagonalMatrix {
0023: /**
0024: * Tridiagonal data.
0025: */
0026: protected final int ldiag[];
0027: protected final int diag[];
0028: protected final int udiag[];
0029:
0030: /**
0031: * Constructs an empty matrix.
0032: * @param size the number of rows/columns
0033: */
0034: public IntegerTridiagonalMatrix(final int size) {
0035: super (size);
0036: ldiag = new int[size];
0037: diag = new int[size];
0038: udiag = new int[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 IntegerTridiagonalMatrix(final int 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 int matrix
0066: */
0067: public boolean equals(AbstractIntegerMatrix m, double tol) {
0068: if (m instanceof TridiagonalMatrix) {
0069: if (numRows != m.rows() || numCols != m.columns())
0070: return false;
0071: int sumSqr = 0;
0072: int ldelta;
0073: int delta = diag[0] - m.getElement(0, 0);
0074: int 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 a double matrix.
0110: * @return a double matrix
0111: */
0112: public AbstractDoubleMatrix toDoubleMatrix() {
0113: final DoubleTridiagonalMatrix m = new DoubleTridiagonalMatrix(
0114: numRows);
0115: m.diag[0] = diag[0];
0116: m.udiag[0] = udiag[0];
0117: int i = 1;
0118: for (; i < numRows - 1; i++) {
0119: m.ldiag[i] = ldiag[i];
0120: m.diag[i] = diag[i];
0121: m.udiag[i] = udiag[i];
0122: }
0123: m.ldiag[i] = ldiag[i];
0124: m.diag[i] = 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 int 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 int 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 int trace() {
0210: int 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 int infNorm() {
0221: int result = Math.abs(diag[0]) + Math.abs(udiag[0]);
0222: int 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 int matrix
0272: * @exception MatrixDimensionException If the matrices are different sizes.
0273: */
0274: public AbstractIntegerSquareMatrix add(
0275: final AbstractIntegerSquareMatrix m) {
0276: if (m instanceof IntegerTridiagonalMatrix)
0277: return add((IntegerTridiagonalMatrix) m);
0278: if (m instanceof TridiagonalMatrix)
0279: return addTridiagonal(m);
0280: if (m instanceof IntegerSquareMatrix)
0281: return add((IntegerSquareMatrix) m);
0282:
0283: if (numRows == m.rows() && numCols == m.columns()) {
0284: final int array[][] = new int[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 IntegerSquareMatrix(array);
0291: } else {
0292: throw new MatrixDimensionException(
0293: "Matrices are different sizes.");
0294: }
0295: }
0296:
0297: public IntegerSquareMatrix add(final IntegerSquareMatrix m) {
0298: if (numRows == m.numRows && numCols == m.numCols) {
0299: final int array[][] = new int[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 IntegerSquareMatrix(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 int tridiagonal matrix
0321: * @exception MatrixDimensionException If the matrices are different sizes.
0322: */
0323: public IntegerTridiagonalMatrix add(final IntegerTridiagonalMatrix m) {
0324: int mRow = numRows;
0325: if (mRow == m.numRows) {
0326: final IntegerTridiagonalMatrix ans = new IntegerTridiagonalMatrix(
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 IntegerTridiagonalMatrix addTridiagonal(
0345: final AbstractIntegerSquareMatrix m) {
0346: int mRow = numRows;
0347: if (mRow == m.rows()) {
0348: final IntegerTridiagonalMatrix ans = new IntegerTridiagonalMatrix(
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 int matrix
0373: * @exception MatrixDimensionException If the matrices are different sizes.
0374: */
0375: public AbstractIntegerSquareMatrix subtract(
0376: final AbstractIntegerSquareMatrix m) {
0377: if (m instanceof IntegerTridiagonalMatrix)
0378: return subtract((IntegerTridiagonalMatrix) m);
0379: if (m instanceof TridiagonalMatrix)
0380: return subtractTridiagonal(m);
0381: if (m instanceof IntegerSquareMatrix)
0382: return subtract((IntegerSquareMatrix) m);
0383:
0384: if (numRows == m.rows() && numCols == m.columns()) {
0385: final int array[][] = new int[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 IntegerSquareMatrix(array);
0392: } else {
0393: throw new MatrixDimensionException(
0394: "Matrices are different sizes.");
0395: }
0396: }
0397:
0398: public IntegerSquareMatrix subtract(final IntegerSquareMatrix m) {
0399: if (numRows == m.numRows && numCols == m.numCols) {
0400: final int array[][] = new int[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 IntegerSquareMatrix(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 int tridiagonal matrix
0415: * @exception MatrixDimensionException If the matrices are different sizes.
0416: */
0417: public IntegerTridiagonalMatrix subtract(
0418: final IntegerTridiagonalMatrix m) {
0419: int mRow = numRows;
0420: if (mRow == m.numRows) {
0421: final IntegerTridiagonalMatrix ans = new IntegerTridiagonalMatrix(
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 IntegerTridiagonalMatrix subtractTridiagonal(
0440: final AbstractIntegerSquareMatrix m) {
0441: int mRow = numRows;
0442: if (mRow == m.rows()) {
0443: final IntegerTridiagonalMatrix ans = new IntegerTridiagonalMatrix(
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 int.
0468: * @return a int square matrix.
0469: */
0470: public AbstractIntegerMatrix scalarMultiply(final int x) {
0471: final IntegerTridiagonalMatrix ans = new IntegerTridiagonalMatrix(
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: // SCALAR PRODUCT
0489:
0490: /**
0491: * Returns the scalar product of this matrix and another.
0492: * @param m a int matrix.
0493: * @exception MatrixDimensionException If the matrices are different sizes.
0494: */
0495: public int scalarProduct(final AbstractIntegerSquareMatrix m) {
0496: if (m instanceof IntegerTridiagonalMatrix)
0497: return scalarProduct((IntegerTridiagonalMatrix) m);
0498: if (m instanceof IntegerSquareMatrix)
0499: return scalarProduct((IntegerSquareMatrix) m);
0500:
0501: if (numRows == m.rows() && numCols == m.columns()) {
0502: int ans = diag[0] * m.getElement(0, 0) + udiag[0]
0503: * m.getElement(0, 1);
0504: final int lastRow = numRows - 1;
0505: for (int i = 1; i < lastRow; i++)
0506: ans += ldiag[i] * m.getElement(i, i - 1) + diag[i]
0507: * m.getElement(i, i) + udiag[i]
0508: * m.getElement(i, i + 1);
0509: ans += ldiag[lastRow] * m.getElement(lastRow, lastRow - 1)
0510: + diag[lastRow] * m.getElement(lastRow, lastRow);
0511: return ans;
0512: } else {
0513: throw new MatrixDimensionException(
0514: "Matrices are different sizes.");
0515: }
0516: }
0517:
0518: public int scalarProduct(final IntegerSquareMatrix m) {
0519: if (numRows == m.numRows && numCols == m.numCols) {
0520: int ans = diag[0] * m.matrix[0][0] + udiag[0]
0521: * m.matrix[0][1];
0522: final int lastRow = numRows - 1;
0523: for (int i = 1; i < lastRow; i++)
0524: ans += ldiag[i] * m.matrix[i][i - 1] + diag[i]
0525: * m.matrix[i][i] + udiag[i]
0526: * m.matrix[i][i + 1];
0527: ans += ldiag[lastRow] * m.matrix[lastRow][lastRow - 1]
0528: + diag[lastRow] * m.matrix[lastRow][lastRow];
0529: return ans;
0530: } else
0531: throw new MatrixDimensionException(
0532: "Matrices are different sizes.");
0533: }
0534:
0535: public int scalarProduct(final IntegerTridiagonalMatrix m) {
0536: if (numRows == m.numRows) {
0537: int ans = diag[0] * m.diag[0] + udiag[0] * m.udiag[0];
0538: final int lastRow = numRows - 1;
0539: for (int i = 1; i < lastRow; i++)
0540: ans += ldiag[i] * m.ldiag[i] + diag[i] * m.diag[i]
0541: + udiag[i] * m.udiag[i];
0542: ans += ldiag[lastRow] * m.ldiag[lastRow] + diag[lastRow]
0543: * m.diag[lastRow];
0544: return ans;
0545: } else
0546: throw new MatrixDimensionException(
0547: "Matrices are different sizes.");
0548: }
0549:
0550: // MATRIX MULTIPLICATION
0551:
0552: /**
0553: * Returns the multiplication of a vector by this matrix.
0554: * @param v a int vector.
0555: * @exception DimensionException If the matrix and vector are incompatible.
0556: */
0557: public AbstractIntegerVector multiply(final AbstractIntegerVector v) {
0558: if (numCols == v.dimension()) {
0559: final int array[] = new int[numRows];
0560: final int lastRow = numRows - 1;
0561: array[0] = diag[0] * v.getComponent(0) + udiag[0]
0562: * v.getComponent(1);
0563: for (int i = 1; i < lastRow; i++)
0564: array[i] = ldiag[i] * v.getComponent(i - 1) + diag[i]
0565: * v.getComponent(i) + udiag[i]
0566: * v.getComponent(i + 1);
0567: array[lastRow] = ldiag[lastRow]
0568: * v.getComponent(lastRow - 1) + diag[lastRow]
0569: * v.getComponent(lastRow);
0570: return new IntegerVector(array);
0571: } else {
0572: throw new DimensionException(
0573: "Matrix and vector are incompatible.");
0574: }
0575: }
0576:
0577: /**
0578: * Returns the multiplication of this matrix and another.
0579: * @param m a int matrix
0580: * @return a AbstractIntegerMatrix or a AbstractIntegerSquareMatrix as appropriate
0581: * @exception MatrixDimensionException If the matrices are incompatible.
0582: */
0583: public AbstractIntegerSquareMatrix multiply(
0584: final AbstractIntegerSquareMatrix m) {
0585: if (m instanceof IntegerTridiagonalMatrix)
0586: return multiply((IntegerTridiagonalMatrix) m);
0587: if (m instanceof TridiagonalMatrix)
0588: return multiplyTridiagonal(m);
0589: if (m instanceof IntegerSquareMatrix)
0590: return multiply((IntegerSquareMatrix) m);
0591:
0592: if (numCols == m.rows()) {
0593: final int mColumns = m.columns();
0594: final int array[][] = new int[numRows][mColumns];
0595: final int lastRow = numRows - 1;
0596: for (int j = 0; j < numRows; j++) {
0597: array[0][j] = diag[0] * m.getElement(0, j) + udiag[0]
0598: * m.getElement(1, j);
0599: for (int i = 1; i < lastRow; i++)
0600: array[i][j] = ldiag[i] * m.getElement(i - 1, j)
0601: + diag[i] * m.getElement(i, j) + udiag[i]
0602: * m.getElement(i + 1, j);
0603: array[lastRow][j] = ldiag[lastRow]
0604: * m.getElement(lastRow - 1, j) + diag[lastRow]
0605: * m.getElement(lastRow, j);
0606: }
0607: return new IntegerSquareMatrix(array);
0608: } else {
0609: throw new MatrixDimensionException("Incompatible matrices.");
0610: }
0611: }
0612:
0613: public IntegerSquareMatrix multiply(final IntegerSquareMatrix m) {
0614: if (numCols == m.numRows) {
0615: final int array[][] = new int[numRows][m.numCols];
0616: final int lastRow = numRows - 1;
0617: for (int j = 0; j < numRows; j++) {
0618: array[0][j] = diag[0] * m.matrix[0][j] + udiag[0]
0619: * m.matrix[1][j];
0620: for (int i = 1; i < lastRow; i++)
0621: array[i][j] = ldiag[i] * m.matrix[i - 1][j]
0622: + diag[i] * m.matrix[i][j] + udiag[i]
0623: * m.matrix[i + 1][j];
0624: array[lastRow][j] = ldiag[lastRow]
0625: * m.matrix[lastRow - 1][j] + diag[lastRow]
0626: * m.matrix[lastRow][j];
0627: }
0628: return new IntegerSquareMatrix(array);
0629: } else
0630: throw new MatrixDimensionException("Incompatible matrices.");
0631: }
0632:
0633: public IntegerSquareMatrix multiply(final IntegerTridiagonalMatrix m) {
0634: int mRow = numRows;
0635: if (numCols == m.numRows) {
0636: final int array[][] = new int[mRow][mRow];
0637: array[0][0] = diag[0] * m.diag[0] + udiag[0] * m.ldiag[1];
0638: array[0][1] = diag[0] * m.udiag[0] + udiag[0] * m.diag[1];
0639: array[0][2] = udiag[0] * m.udiag[1];
0640: if (mRow > 3) {
0641: array[1][0] = ldiag[1] * m.diag[0] + diag[1]
0642: * m.ldiag[1];
0643: array[1][1] = ldiag[1] * m.udiag[0] + diag[1]
0644: * m.diag[1] + udiag[1] * m.ldiag[2];
0645: array[1][2] = diag[1] * m.udiag[1] + udiag[1]
0646: * m.diag[2];
0647: array[1][3] = udiag[1] * m.udiag[2];
0648: }
0649: if (mRow == 3) {
0650: array[1][0] = ldiag[1] * m.diag[0] + diag[1]
0651: * m.ldiag[1];
0652: array[1][1] = ldiag[1] * m.udiag[0] + diag[1]
0653: * m.diag[1] + udiag[1] * m.ldiag[2];
0654: array[1][2] = diag[1] * m.udiag[1] + udiag[1]
0655: * m.diag[2];
0656: } else if (mRow > 4) {
0657: for (int i = 2; i < mRow - 2; i++) {
0658: array[i][i - 2] = ldiag[i] * m.ldiag[i - 1];
0659: array[i][i - 1] = ldiag[i] * m.diag[i - 1]
0660: + diag[i] * m.ldiag[i];
0661: array[i][i] = ldiag[i] * m.udiag[i - 1] + diag[i]
0662: * m.diag[i] + udiag[i] * m.ldiag[i + 1];
0663: array[i][i + 1] = diag[i] * m.udiag[i] + udiag[i]
0664: * m.diag[i + 1];
0665: array[i][i + 2] = udiag[i] * m.udiag[i + 1];
0666: }
0667: }
0668: if (mRow > 3) {
0669: array[mRow - 2][mRow - 4] = ldiag[mRow - 2]
0670: * m.ldiag[mRow - 3];
0671: array[mRow - 2][mRow - 3] = ldiag[mRow - 2]
0672: * m.diag[mRow - 3] + diag[mRow - 2]
0673: * m.ldiag[mRow - 2];
0674: array[mRow - 2][mRow - 2] = ldiag[mRow - 2]
0675: * m.udiag[mRow - 3] + diag[mRow - 2]
0676: * m.diag[mRow - 2] + udiag[mRow - 2]
0677: * m.ldiag[mRow - 1];
0678: array[mRow - 2][mRow - 1] = diag[mRow - 2]
0679: * m.udiag[mRow - 2] + udiag[mRow - 2]
0680: * m.diag[mRow - 1];
0681: }
0682: mRow--;
0683: array[mRow][mRow - 2] = ldiag[mRow] * m.ldiag[mRow - 1];
0684: array[mRow][mRow - 1] = ldiag[mRow] * m.diag[mRow - 1]
0685: + diag[mRow] * m.ldiag[mRow];
0686: array[mRow][mRow] = ldiag[mRow] * m.udiag[mRow - 1]
0687: + diag[mRow] * m.diag[mRow];
0688: return new IntegerSquareMatrix(array);
0689: } else
0690: throw new MatrixDimensionException("Incompatible matrices.");
0691: }
0692:
0693: private IntegerSquareMatrix multiplyTridiagonal(
0694: final AbstractIntegerSquareMatrix m) {
0695: int mRow = numRows;
0696: if (numCols == m.rows()) {
0697: final int array[][] = new int[mRow][mRow];
0698: array[0][0] = diag[0] * m.getElement(0, 0) + udiag[0]
0699: * m.getElement(1, 0);
0700: array[0][1] = diag[0] * m.getElement(0, 1) + udiag[0]
0701: * m.getElement(1, 1);
0702: array[0][2] = udiag[0] * m.getElement(1, 2);
0703: if (mRow > 3) {
0704: array[1][0] = ldiag[1] * m.getElement(0, 0) + diag[1]
0705: * m.getElement(1, 0);
0706: array[1][1] = ldiag[1] * m.getElement(0, 1) + diag[1]
0707: * m.getElement(1, 1) + udiag[1]
0708: * m.getElement(2, 1);
0709: array[1][2] = diag[1] * m.getElement(1, 2) + udiag[1]
0710: * m.getElement(2, 2);
0711: array[1][3] = udiag[1] * m.getElement(2, 3);
0712: }
0713: if (mRow == 3) {
0714: array[1][0] = ldiag[1] * m.getElement(0, 0) + diag[1]
0715: * m.getElement(1, 0);
0716: array[1][1] = ldiag[1] * m.getElement(0, 1) + diag[1]
0717: * m.getElement(1, 1) + udiag[1]
0718: * m.getElement(2, 1);
0719: array[1][2] = diag[1] * m.getElement(1, 2) + udiag[1]
0720: * m.getElement(2, 2);
0721: } else if (mRow > 4) {
0722: for (int i = 2; i < mRow - 2; i++) {
0723: array[i][i - 2] = ldiag[i]
0724: * m.getElement(i - 1, i - 2);
0725: array[i][i - 1] = ldiag[i]
0726: * m.getElement(i - 1, i - 1) + diag[i]
0727: * m.getElement(i, i - 1);
0728: array[i][i] = ldiag[i] * m.getElement(i - 1, i)
0729: + diag[i] * m.getElement(i, i) + udiag[i]
0730: * m.getElement(i + 1, i);
0731: array[i][i + 1] = diag[i] * m.getElement(i, i + 1)
0732: + udiag[i] * m.getElement(i + 1, i + 1);
0733: array[i][i + 2] = udiag[i]
0734: * m.getElement(i + 1, i + 2);
0735: }
0736: }
0737: if (mRow > 3) {
0738: array[mRow - 2][mRow - 4] = ldiag[mRow - 2]
0739: * m.getElement(mRow - 3, mRow - 4);
0740: array[mRow - 2][mRow - 3] = ldiag[mRow - 2]
0741: * m.getElement(mRow - 3, mRow - 3)
0742: + diag[mRow - 2]
0743: * m.getElement(mRow - 2, mRow - 3);
0744: array[mRow - 2][mRow - 2] = ldiag[mRow - 2]
0745: * m.getElement(mRow - 3, mRow - 2)
0746: + diag[mRow - 2]
0747: * m.getElement(mRow - 2, mRow - 2)
0748: + udiag[mRow - 2]
0749: * m.getElement(mRow - 1, mRow - 2);
0750: array[mRow - 2][mRow - 1] = diag[mRow - 2]
0751: * m.getElement(mRow - 2, mRow - 1)
0752: + udiag[mRow - 2]
0753: * m.getElement(mRow - 1, mRow - 1);
0754: }
0755: mRow--;
0756: array[mRow][mRow - 2] = ldiag[mRow]
0757: * m.getElement(mRow - 1, mRow - 2);
0758: array[mRow][mRow - 1] = ldiag[mRow]
0759: * m.getElement(mRow - 1, mRow - 1) + diag[mRow]
0760: * m.getElement(mRow, mRow - 1);
0761: array[mRow][mRow] = ldiag[mRow]
0762: * m.getElement(mRow - 1, mRow) + diag[mRow]
0763: * m.getElement(mRow, mRow);
0764: return new IntegerSquareMatrix(array);
0765: } else {
0766: throw new MatrixDimensionException("Incompatible matrices.");
0767: }
0768: }
0769:
0770: // TRANSPOSE
0771:
0772: /**
0773: * Returns the transpose of this matrix.
0774: * @return a int matrix
0775: */
0776: public Matrix transpose() {
0777: final IntegerTridiagonalMatrix ans = new IntegerTridiagonalMatrix(
0778: numRows);
0779: System.arraycopy(ldiag, 1, ans.udiag, 0, ldiag.length - 1);
0780: System.arraycopy(diag, 0, ans.diag, 0, diag.length);
0781: System.arraycopy(udiag, 0, ans.ldiag, 1, udiag.length - 1);
0782: return ans;
0783: }
0784:
0785: // CHOLESKY DECOMPOSITION
0786:
0787: /**
0788: * Returns the Cholesky decomposition of this matrix.
0789: * Matrix must be symmetric and positive definite.
0790: * @return an array with [0] containing the L-matrix and [1] containing the U-matrix.
0791: */
0792: public AbstractDoubleSquareMatrix[] choleskyDecompose() {
0793: final int N = numRows;
0794: final double arrayL[][] = new double[N][N];
0795: final double arrayU[][] = new double[N][N];
0796: double tmp = Math.sqrt(diag[0]);
0797: arrayL[0][0] = arrayU[0][0] = tmp;
0798: arrayL[1][0] = arrayU[0][1] = ldiag[1] / tmp;
0799: for (int j = 1; j < N; j++) {
0800: tmp = diag[j];
0801: for (int i = 0; i < j; i++)
0802: tmp -= arrayL[j][i] * arrayL[j][i];
0803: arrayL[j][j] = arrayU[j][j] = Math.sqrt(tmp);
0804: if (j + 1 < N) {
0805: tmp = ldiag[j + 1];
0806: for (int k = 0; k < j + 1; k++)
0807: tmp -= arrayL[j][k] * arrayU[k][j + 1];
0808: arrayL[j + 1][j] = arrayU[j][j + 1] = tmp
0809: / arrayU[j][j];
0810: for (int i = j + 2; i < N; i++) {
0811: tmp = 0.0;
0812: for (int k = 0; k < i; k++)
0813: tmp -= arrayL[j][k] * arrayU[k][i];
0814: arrayL[i][j] = arrayU[j][i] = tmp / arrayU[j][j];
0815: }
0816: }
0817: }
0818: final AbstractDoubleSquareMatrix lu[] = new AbstractDoubleSquareMatrix[2];
0819: lu[0] = new DoubleSquareMatrix(arrayL);
0820: lu[1] = new DoubleSquareMatrix(arrayU);
0821: return lu;
0822: }
0823:
0824: // QR DECOMPOSITION
0825:
0826: /**
0827: * Returns the QR decomposition of this matrix.
0828: * Based on the code from <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> (public domain).
0829: * @return an array with [0] containing the Q-matrix and [1] containing the R-matrix.
0830: * @jsci.planetmath QRDecomposition
0831: */
0832: public AbstractDoubleSquareMatrix[] qrDecompose() {
0833: final int N = numRows;
0834: final double array[][] = new double[N][N];
0835: final double arrayQ[][] = new double[N][N];
0836: final double arrayR[][] = new double[N][N];
0837: // copy matrix
0838: array[0][0] = diag[0];
0839: array[0][1] = udiag[0];
0840: for (int i = 1; i < N - 1; i++) {
0841: array[i][i - 1] = ldiag[i];
0842: array[i][i] = diag[i];
0843: array[i][i + 1] = udiag[i];
0844: }
0845: array[N - 1][N - 2] = ldiag[N - 1];
0846: array[N - 1][N - 1] = diag[N - 1];
0847: for (int k = 0; k < N; k++) {
0848: // compute l2-norm of kth column
0849: double norm = array[k][k];
0850: for (int i = k + 1; i < N; i++)
0851: norm = ExtraMath.hypot(norm, array[i][k]);
0852: if (norm != 0.0) {
0853: // form kth Householder vector
0854: if (array[k][k] < 0.0)
0855: norm = -norm;
0856: for (int i = k; i < N; i++)
0857: array[i][k] /= norm;
0858: array[k][k] += 1.0;
0859: // apply transformation to remaining columns
0860: for (int j = k + 1; j < N; j++) {
0861: double s = array[k][k] * array[k][j];
0862: for (int i = k + 1; i < N; i++)
0863: s += array[i][k] * array[i][j];
0864: s /= array[k][k];
0865: for (int i = k; i < N; i++)
0866: array[i][j] -= s * array[i][k];
0867: }
0868: }
0869: arrayR[k][k] = -norm;
0870: }
0871: for (int k = N - 1; k >= 0; k--) {
0872: arrayQ[k][k] = 1.0;
0873: for (int j = k; j < N; j++) {
0874: if (array[k][k] != 0.0) {
0875: double s = array[k][k] * arrayQ[k][j];
0876: for (int i = k + 1; i < N; i++)
0877: s += array[i][k] * arrayQ[i][j];
0878: s /= array[k][k];
0879: for (int i = k; i < N; i++)
0880: arrayQ[i][j] -= s * array[i][k];
0881: }
0882: }
0883: }
0884: for (int i = 0; i < N; i++) {
0885: for (int j = i + 1; j < N; j++)
0886: arrayR[i][j] = array[i][j];
0887: }
0888: final AbstractDoubleSquareMatrix qr[] = new AbstractDoubleSquareMatrix[2];
0889: qr[0] = new DoubleSquareMatrix(arrayQ);
0890: qr[1] = new DoubleSquareMatrix(arrayR);
0891: return qr;
0892: }
0893:
0894: // SINGULAR VALUE DECOMPOSITION
0895:
0896: /**
0897: * Returns the singular value decomposition of this matrix.
0898: * Based on the code from <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> (public domain).
0899: * @return an array with [0] containing the U-matrix, [1] containing the S-matrix and [2] containing the V-matrix.
0900: */
0901: public AbstractDoubleSquareMatrix[] singularValueDecompose() {
0902: final int N = numRows;
0903: final int Nm1 = N - 1;
0904: final double array[][] = new double[N][N];
0905: final double arrayU[][] = new double[N][N];
0906: final double arrayS[] = new double[N];
0907: final double arrayV[][] = new double[N][N];
0908: final double e[] = new double[N];
0909: final double work[] = new double[N];
0910: // copy matrix
0911: array[0][0] = diag[0];
0912: array[0][1] = udiag[0];
0913: for (int i = 1; i < Nm1; i++) {
0914: array[i][i - 1] = ldiag[i];
0915: array[i][i] = diag[i];
0916: array[i][i + 1] = udiag[i];
0917: }
0918: array[Nm1][Nm1 - 1] = ldiag[Nm1];
0919: array[Nm1][Nm1] = diag[Nm1];
0920: // reduce matrix to bidiagonal form
0921: for (int k = 0; k < Nm1; k++) {
0922: // compute the transformation for the kth column
0923: // compute l2-norm of kth column
0924: arrayS[k] = array[k][k];
0925: for (int i = k + 1; i < N; i++)
0926: arrayS[k] = ExtraMath.hypot(arrayS[k], array[i][k]);
0927: if (arrayS[k] != 0.0) {
0928: if (array[k][k] < 0.0)
0929: arrayS[k] = -arrayS[k];
0930: for (int i = k; i < N; i++)
0931: array[i][k] /= arrayS[k];
0932: array[k][k] += 1.0;
0933: }
0934: arrayS[k] = -arrayS[k];
0935: for (int j = k + 1; j < N; j++) {
0936: if (arrayS[k] != 0.0) {
0937: // apply the transformation
0938: double t = array[k][k] * array[k][j];
0939: for (int i = k + 1; i < N; i++)
0940: t += array[i][k] * array[i][j];
0941: t /= array[k][k];
0942: for (int i = k; i < N; i++)
0943: array[i][j] -= t * array[i][k];
0944: }
0945: e[j] = array[k][j];
0946: }
0947: for (int i = k; i < N; i++)
0948: arrayU[i][k] = array[i][k];
0949: if (k < N - 2) {
0950: // compute the kth row transformation
0951: // compute l2-norm of kth column
0952: e[k] = e[k + 1];
0953: for (int i = k + 2; i < N; i++)
0954: e[k] = ExtraMath.hypot(e[k], e[i]);
0955: if (e[k] != 0.0) {
0956: if (e[k + 1] < 0.0)
0957: e[k] = -e[k];
0958: for (int i = k + 1; i < N; i++)
0959: e[i] /= e[k];
0960: e[k + 1] += 1.0;
0961: }
0962: e[k] = -e[k];
0963: if (e[k] != 0.0) {
0964: // apply the transformation
0965: for (int i = k + 1; i < N; i++) {
0966: work[i] = 0.0;
0967: for (int j = k + 1; j < N; j++)
0968: work[i] += e[j] * array[i][j];
0969: }
0970: for (int j = k + 1; j < N; j++) {
0971: double t = e[j] / e[k + 1];
0972: for (int i = k + 1; i < N; i++)
0973: array[i][j] -= t * work[i];
0974: }
0975: }
0976: for (int i = k + 1; i < N; i++)
0977: arrayV[i][k] = e[i];
0978: }
0979: }
0980: // setup the final bidiagonal matrix of order p
0981: int p = N;
0982: arrayS[Nm1] = array[Nm1][Nm1];
0983: e[N - 2] = array[N - 2][Nm1];
0984: e[Nm1] = 0.0;
0985: arrayU[Nm1][Nm1] = 1.0;
0986: for (int k = N - 2; k >= 0; k--) {
0987: if (arrayS[k] != 0.0) {
0988: for (int j = k + 1; j < N; j++) {
0989: double t = arrayU[k][k] * arrayU[k][j];
0990: for (int i = k + 1; i < N; i++)
0991: t += arrayU[i][k] * arrayU[i][j];
0992: t /= arrayU[k][k];
0993: for (int i = k; i < N; i++)
0994: arrayU[i][j] -= t * arrayU[i][k];
0995: }
0996: for (int i = k; i < N; i++)
0997: arrayU[i][k] = -arrayU[i][k];
0998: arrayU[k][k] += 1.0;
0999: for (int i = 0; i < k - 1; i++)
1000: arrayU[i][k] = 0.0;
1001: } else {
1002: for (int i = 0; i < N; i++)
1003: arrayU[i][k] = 0.0;
1004: arrayU[k][k] = 1.0;
1005: }
1006: }
1007: for (int k = Nm1; k >= 0; k--) {
1008: if (k < N - 2 && e[k] != 0.0) {
1009: for (int j = k + 1; j < N; j++) {
1010: double t = arrayV[k + 1][k] * arrayV[k + 1][j];
1011: for (int i = k + 2; i < N; i++)
1012: t += arrayV[i][k] * arrayV[i][j];
1013: t /= arrayV[k + 1][k];
1014: for (int i = k + 1; i < N; i++)
1015: arrayV[i][j] -= t * arrayV[i][k];
1016: }
1017: }
1018: for (int i = 0; i < N; i++)
1019: arrayV[i][k] = 0.0;
1020: arrayV[k][k] = 1.0;
1021: }
1022: final double eps = Math.pow(2.0, -52.0);
1023: int iter = 0;
1024: while (p > 0) {
1025: int k, action;
1026: // action = 1 if arrayS[p] and e[k-1] are negligible and k<p
1027: // action = 2 if arrayS[k] is negligible and k<p
1028: // action = 3 if e[k-1] is negligible, k<p, and arrayS[k], ..., arrayS[p] are not negligible (QR step)
1029: // action = 4 if e[p-1] is negligible (convergence)
1030: for (k = p - 2; k >= -1; k--) {
1031: if (k == -1)
1032: break;
1033: if (Math.abs(e[k]) <= eps
1034: * (Math.abs(arrayS[k]) + Math
1035: .abs(arrayS[k + 1]))) {
1036: e[k] = 0.0;
1037: break;
1038: }
1039: }
1040: if (k == p - 2) {
1041: action = 4;
1042: } else {
1043: int ks;
1044: for (ks = p - 1; ks >= k; ks--) {
1045: if (ks == k)
1046: break;
1047: double t = (ks != p ? Math.abs(e[ks]) : 0.0)
1048: + (ks != k + 1 ? Math.abs(e[ks - 1]) : 0.0);
1049: if (Math.abs(arrayS[ks]) <= eps * t) {
1050: arrayS[ks] = 0.0;
1051: break;
1052: }
1053: }
1054: if (ks == k) {
1055: action = 3;
1056: } else if (ks == p - 1) {
1057: action = 1;
1058: } else {
1059: action = 2;
1060: k = ks;
1061: }
1062: }
1063: k++;
1064: switch (action) {
1065: // deflate negligible arrayS[p]
1066: case 1: {
1067: double f = e[p - 2];
1068: e[p - 2] = 0.0;
1069: for (int j = p - 2; j >= k; j--) {
1070: double t = ExtraMath.hypot(arrayS[j], f);
1071: final double cs = arrayS[j] / t;
1072: final double sn = f / t;
1073: arrayS[j] = t;
1074: if (j != k) {
1075: f = -sn * e[j - 1];
1076: e[j - 1] *= cs;
1077: }
1078: for (int i = 0; i < N; i++) {
1079: t = cs * arrayV[i][j] + sn * arrayV[i][p - 1];
1080: arrayV[i][p - 1] = -sn * arrayV[i][j] + cs
1081: * arrayV[i][p - 1];
1082: arrayV[i][j] = t;
1083: }
1084: }
1085: }
1086: break;
1087: // split at negligible arrayS[k]
1088: case 2: {
1089: double f = e[k - 1];
1090: e[k - 1] = 0.0;
1091: for (int j = k; j < p; j++) {
1092: double t = ExtraMath.hypot(arrayS[j], f);
1093: final double cs = arrayS[j] / t;
1094: final double sn = f / t;
1095: arrayS[j] = t;
1096: f = -sn * e[j];
1097: e[j] *= cs;
1098: for (int i = 0; i < N; i++) {
1099: t = cs * arrayU[i][j] + sn * arrayU[i][k - 1];
1100: arrayU[i][k - 1] = -sn * arrayU[i][j] + cs
1101: * arrayU[i][k - 1];
1102: arrayU[i][j] = t;
1103: }
1104: }
1105: }
1106: break;
1107: // perform one QR step
1108: case 3: {
1109: // calculate the shift
1110: final double scale = Math.max(Math.max(Math.max(Math
1111: .max(Math.abs(arrayS[p - 1]), Math
1112: .abs(arrayS[p - 2])), Math
1113: .abs(e[p - 2])), Math.abs(arrayS[k])), Math
1114: .abs(e[k]));
1115: double sp = arrayS[p - 1] / scale;
1116: double spm1 = arrayS[p - 2] / scale;
1117: double epm1 = e[p - 2] / scale;
1118: double sk = arrayS[k] / scale;
1119: double ek = e[k] / scale;
1120: double b = ((spm1 + sp) * (spm1 - sp) + epm1 * epm1) / 2.0;
1121: double c = (sp * epm1) * (sp * epm1);
1122: double shift = 0.0;
1123: if (b != 0.0 || c != 0.0) {
1124: shift = Math.sqrt(b * b + c);
1125: if (b < 0.0)
1126: shift = -shift;
1127: shift = c / (b + shift);
1128: }
1129: double f = (sk + sp) * (sk - sp) + shift;
1130: double g = sk * ek;
1131: // chase zeros
1132: for (int j = k; j < p - 1; j++) {
1133: double t = ExtraMath.hypot(f, g);
1134: double cs = f / t;
1135: double sn = g / t;
1136: if (j != k)
1137: e[j - 1] = t;
1138: f = cs * arrayS[j] + sn * e[j];
1139: e[j] = cs * e[j] - sn * arrayS[j];
1140: g = sn * arrayS[j + 1];
1141: arrayS[j + 1] *= cs;
1142: for (int i = 0; i < N; i++) {
1143: t = cs * arrayV[i][j] + sn * arrayV[i][j + 1];
1144: arrayV[i][j + 1] = -sn * arrayV[i][j] + cs
1145: * arrayV[i][j + 1];
1146: arrayV[i][j] = t;
1147: }
1148: t = ExtraMath.hypot(f, g);
1149: cs = f / t;
1150: sn = g / t;
1151: arrayS[j] = t;
1152: f = cs * e[j] + sn * arrayS[j + 1];
1153: arrayS[j + 1] = -sn * e[j] + cs * arrayS[j + 1];
1154: g = sn * e[j + 1];
1155: e[j + 1] *= cs;
1156: if (j < Nm1) {
1157: for (int i = 0; i < N; i++) {
1158: t = cs * arrayU[i][j] + sn
1159: * arrayU[i][j + 1];
1160: arrayU[i][j + 1] = -sn * arrayU[i][j] + cs
1161: * arrayU[i][j + 1];
1162: arrayU[i][j] = t;
1163: }
1164: }
1165: }
1166: e[p - 2] = f;
1167: iter++;
1168: }
1169: break;
1170: // convergence
1171: case 4: {
1172: // make the singular values positive
1173: if (arrayS[k] <= 0.0) {
1174: arrayS[k] = -arrayS[k];
1175: for (int i = 0; i < p; i++)
1176: arrayV[i][k] = -arrayV[i][k];
1177: }
1178: // order the singular values
1179: while (k < p - 1) {
1180: if (arrayS[k] >= arrayS[k + 1])
1181: break;
1182: double tmp = arrayS[k];
1183: arrayS[k] = arrayS[k + 1];
1184: arrayS[k + 1] = tmp;
1185: if (k < Nm1) {
1186: for (int i = 0; i < N; i++) {
1187: tmp = arrayU[i][k + 1];
1188: arrayU[i][k + 1] = arrayU[i][k];
1189: arrayU[i][k] = tmp;
1190: tmp = arrayV[i][k + 1];
1191: arrayV[i][k + 1] = arrayV[i][k];
1192: arrayV[i][k] = tmp;
1193: }
1194: }
1195: k++;
1196: }
1197: iter = 0;
1198: p--;
1199: }
1200: break;
1201: }
1202: }
1203: final AbstractDoubleSquareMatrix svd[] = new AbstractDoubleSquareMatrix[3];
1204: svd[0] = new DoubleSquareMatrix(arrayU);
1205: svd[1] = new DoubleDiagonalMatrix(arrayS);
1206: svd[2] = new DoubleSquareMatrix(arrayV);
1207: return svd;
1208: }
1209:
1210: }
|