Merge pull request #36 from StaticBeagle/add-qr-complex-matrix-inversion

Add qr complex matrix inversion

src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/ComplexLUDecompositionDense.java

@@ -175,6 +175,21 @@ public Complex det() {
 	 */
 
 	public ComplexMatrixDense solve(MatrixDense B) {
+		return solve(ComplexMatrixDense.fromRealMatrix(B));
+	}
+
+	/**
+	 * Solve A*X = B
+	 *
+	 * @param B
+	 *            A Matrix with as many rows as A and any number of columns.
+	 * @return X so that L*U*X = B(piv,:)
+	 * @exception IllegalArgumentException
+	 *                Matrix row dimensions must agree.
+	 * @exception RuntimeException
+	 *                Matrix is singular.
+	 */
+	public ComplexMatrixDense solve(ComplexMatrixDense B) {
 		if (B.getRowCount() != _rows) {
 			throw new IllegalArgumentException("Matrix row dimensions must agree.");
 		}
@@ -184,7 +199,7 @@ public ComplexMatrixDense solve(MatrixDense B) {
 
 		// Copy right hand side with pivoting
 		int nx = B.getColumnCount();
-		ComplexMatrixDense Xmat = ComplexMatrixDense.fromRealMatrix(B.subMatrix(_pivot, 0, nx - 1));
+		ComplexMatrixDense Xmat = B.subMatrix(_pivot, 0, nx - 1);
 		Complex[] X = Xmat.getArray();
 
 		final int cols = _cols;

src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/ComplexMatrixDense.java

@@ -296,18 +296,17 @@ public ComplexMatrixDense inv() {
 
     // region solve
     public ComplexMatrixDense solve(MatrixDense B) {
-        return new ComplexLUDecompositionDense(this).solve(B);
-        // Only implemented for square matrices for now.
-//        if (rows == cols) { // Matrix is Squared
-//            return new LUDecomposition(this).solve(B);
-//        } else if (rows > cols) { // Matrix is thin (Overdetermined system)
-//            return new QRDecomposition(this).solve(B);
-//        } else { // Matrix is fat (Under-determined system)
-//            QRDecomposition qr = this.transpose().QR();
-//            Matrix R1 = fwdSubsSolve(qr.getRT(), B);
-//            R1.appendRows(cols - R1.rows);
-//            return qr.QmultiplyX(R1);
-//        }
+        return solve(fromRealMatrix(B));
+    }
+
+    public ComplexMatrixDense solve(ComplexMatrixDense B) {
+        if (rows == cols) { // Matrix is Squared
+            return new ComplexLUDecompositionDense(this).solve(B);
+        } else if (rows > cols) { // Matrix is tall and narrow (Overdetermined system)
+            return new ComplexQRDecompositionDense(this).solve(B);
+        } else { // Matrix is short and wide (Under-determined system)
+            throw new UnsupportedOperationException("Method not implemented yet for short and wide Matrices");
+        }
     }
 
     // endregion
@@ -359,6 +358,35 @@ public ComplexMatrixDense conjugateTranspose() {
         return new ComplexMatrixDense(result, trows, tcols);
     }
 
+    /***
+     * Get sub-matrix
+     *
+     * @param rows The array of row indices.
+     * @param col0 The initial column index.
+     * @param col1 The final column index.
+     * @return {@code A(rows(:), col0 : col1)}.
+     */
+    public ComplexMatrixDense subMatrix(int[] rows, int col0, int col1) {
+        if (col0 < 0 || col1 < 0) {
+            throw new IllegalArgumentException("The column indexes col0 and col1 must be greater than zero.");
+        }
+        if (Arrays.stream(rows).anyMatch(i -> i >= this.rows || i < 0)) {
+            throw new ArrayIndexOutOfBoundsException("The row indexes cannot be greater than the number of rows in the Matrix. and must be greater than zero.");
+        }
+        if (col1 >= cols) {
+            throw new ArrayIndexOutOfBoundsException("The final column index cannot be greater than the number of columns in the Matrix.");
+        }
+        int rowDim = rows.length;
+        int colDim = col1 - col0 + 1;
+        Complex[] data = new Complex[rowDim * colDim];
+        for (int i = 0; i < rowDim; ++i) {
+            for (int j = 0; j < colDim; ++j) {
+                data[i * colDim + j] = this.data[rows[i] * cols + (j + col0)];
+            }
+        }
+        return new ComplexMatrixDense(data, rowDim, colDim);
+    }
+
     public static final class Factory {
         private Factory() {}
 

src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/ComplexQRDecompositionDense.java

@@ -303,6 +303,19 @@ public ComplexMatrixDense solve(ComplexMatrixDense B) {
         return backSubstitutionSolve(R, Y);
     }
 
+    /**
+     * Least squares solution of A*X = B
+     *
+     * @param B A Matrix with as many rows as A and any number of columns.
+     * @return X that minimizes the two norm of Q*R*X-B.
+     * @throws IllegalArgumentException Matrix row dimensions must agree.
+     * @throws RuntimeException         Matrix is rank deficient.
+     */
+
+    public ComplexMatrixDense solve(MatrixDense B) {
+        return solve(ComplexMatrixDense.fromRealMatrix(B));
+    }
+
     @Override
     public String toString() {
         StringBuilder sb = new StringBuilder();

src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/MatrixDense.java

@@ -931,7 +931,7 @@ public ComplexMatrixDense multiply(ComplexMatrixDense matrix) {
     public MatrixDense solve(MatrixDense b) {
         if (rows == cols) { // Matrix is Squared
             return new LUDecompositionDense(this).solve(b);
-        } else if (rows > cols) { // Matrix is thin (Overdetermined system)
+        } else if (rows > cols) { // Matrix is tall and narrow (Overdetermined system)
             return new QRDecompositionDense(this).solve(b);
         } else { // Matrix is short and wide (Under-determined system)
             return this.pinv().multiply(b);

src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/QRDecomposition.java

@@ -3,7 +3,7 @@
 public abstract class QRDecomposition<T extends Matrix> {
     protected final int rows;
     protected final int cols;
-
+    // TODO decide on what to call the V vectors matrix and also compact Q and compact R getter interface
     public QRDecomposition(T matrix) {
         this.rows = matrix.getRowCount();
         this.cols = matrix.getColumnCount();