Added unit tests

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

@@ -1,20 +1,39 @@
 package com.wildbitsfoundry.etk4j.math.linearalgebra;
 
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
 import com.wildbitsfoundry.etk4j.constants.ConstantsETK;
 import com.wildbitsfoundry.etk4j.math.MathETK;
-import com.wildbitsfoundry.etk4j.math.complex.Complex;
 
 /**
- * Schur implements the Schur decomposition of a matrix.  Specifically,
- * given a square matrix A, there is a unitary matrix U such that
- * <pre>
- *      T = U^H AU
- * </pre>
- * is upper triangular.  Schur represents T as a Zutmat and U as a Zmat.
- *
- * @author G. W. Stewart
- * @version Pre-alpha
- */
+ * Class transforming a general real matrix to Schur form.
+ * <p>A m × m matrix A can be written as the product of three matrices: A = P
+ * × T × P<sup>T with P an orthogonal matrix and T an quasi-triangular
+ * matrix. Both P and T are m × m matrices.</p>
+ * <p>Transformation to Schur form is often not a goal by itself, but it is an
+ * intermediate step in more general decomposition algorithms like
+ * {@link EigenvalueDecompositionDense eigen decomposition}. This class is therefore
+ * intended for internal use by the library and is not public. As a consequence
+ * of this explicitly limited scope, many methods directly returns references to
+ * internal arrays, not copies.</p>
+ * <p>This class is based on the method hqr2 in class EigenvalueDecomposition
+ * from the <a href="http://math.nist.gov/javanumerics/jama/">JAMA library.
+ * */
 public class SchurDecompositionDense {
 
     /**
@@ -25,11 +44,11 @@ public class SchurDecompositionDense {
     /**
      * P matrix.
      */
-    private final double matrixP[][];
+    private final double[][] matrixP;
     /**
      * T matrix.
      */
-    private final double matrixT[][];
+    private final double[][] matrixT;
     /**
      * Cached value of P.
      */
@@ -449,21 +468,4 @@ private static class ShiftInfo {
     }
 
     // TODO align return matrices name with Octave, Matlab
-    public static void main(String[] args) {
-        double[][] matrix = {
-                {65, 35, 40, 69},
-                {99, 64, 37, 2},
-                {39, 48, 35, 90},
-                {30, 93, 87, 17}
-        };
-
-        MatrixDense A = MatrixDense.from2DArray(matrix);
-        SchurDecompositionDense schur = new SchurDecompositionDense(A);
-
-        System.out.println(schur.getT());
-        System.out.println();
-        System.out.println(schur.getP());
-        System.out.println();
-        System.out.println(schur.getP().multiply(schur.getT()).multiply(schur.getPT()));
-    }
 }

src/test/java/com/wildbitsfoundry/etk4j/math/linearalgebra/MatrixDenseTest.java

@@ -68,9 +68,9 @@ public void testExpm() {
     @Test
     public void testSetRow() {
         MatrixDense m = MatrixDense.Factory.magic(3);
-		double[] row = {-1, -1, -1};
+        double[] row = {-1, -1, -1};
         m.setRow(1, row);
-		assertArrayEquals(row, m.getRow(1), 1e-12);
+        assertArrayEquals(row, m.getRow(1), 1e-12);
     }
 
     @Test
@@ -102,6 +102,46 @@ public void testHessembergDecomposition() {
         assertArrayEquals(expectedU, hess.getU().getArray(), 1e-12);
     }
 
+    @Test
+    public void testSchurDecomposition() {
+        double[][] matrix = {
+                {65, 35, 40, 69},
+                {99, 64, 37, 2},
+                {39, 48, 35, 90},
+                {30, 93, 87, 17}
+        };
+
+        MatrixDense A = MatrixDense.from2DArray(matrix);
+        SchurDecompositionDense schur = new SchurDecompositionDense(A);
+
+        double[] expectedT = {211.53821949564204, 27.963104879902804, 11.529192222479173, -11.83821931259974,
+                -2.8398992587956425E-29, 27.374680489564316, 64.83795403092125, 3.6006930114899935, 0.0,
+                -31.22064867637267, -13.419034896325854, 61.10465539906722, 0.0, 0.0, 2.2836324137572356E-24,
+                -44.49386508888031};
+        assertArrayEquals(expectedT, schur.getT().getArray(), 1e-12);
+        double[] expectedP = {0.49782038512887483, 0.15442410294358475, 0.7718547804323013, 0.3640992426578396,
+                0.4680103341616455, 0.7666501603471897, -0.30669802271073127, -0.3148812182758139, 0.505725691055896,
+                -0.5430518840085535, 0.0954550851000274, -0.6634941623023844, 0.5266713554713984, -0.3057701413553275,
+                -0.5486937647884095, 0.5727626528185856};
+        assertArrayEquals(expectedP, schur.getP().getArray(), 1e-12);
+    }
+
+    @Test
+    public void testVerifySchurDecomposition() {
+        double[][] matrix = {
+                {65, 35, 40, 69},
+                {99, 64, 37, 2},
+                {39, 48, 35, 90},
+                {30, 93, 87, 17}
+        };
+
+        MatrixDense A = MatrixDense.from2DArray(matrix);
+        SchurDecompositionDense schur = new SchurDecompositionDense(A);
+        MatrixDense actual = schur.getP().multiply(schur.getT()).multiply(schur.getPT());
+
+        assertArrayEquals(A.getArray(), actual.getArray(), 1e-12);
+    }
+
     @Test
     public void allTests() {
         MatrixDense A, B, C, Z, O, I, R, S, X, SUB, M, T, SQ, DEF, SOL;