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