diff --git a/src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/SchurDecompositionDense.java b/src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/SchurDecompositionDense.java
index a921b67..10b0641 100644
--- a/src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/SchurDecompositionDense.java
+++ b/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
- *
- * T = U^H AU
- *
- * 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.
+ * A m × m matrix A can be written as the product of three matrices: A = P
+ * × T × PT with P an orthogonal matrix and T an quasi-triangular
+ * matrix. Both P and T are m × m matrices.
+ * 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.
+ * This class is based on the method hqr2 in class EigenvalueDecomposition
+ * from the 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()));
- }
}
diff --git a/src/test/java/com/wildbitsfoundry/etk4j/math/linearalgebra/MatrixDenseTest.java b/src/test/java/com/wildbitsfoundry/etk4j/math/linearalgebra/MatrixDenseTest.java
index c3c134e..f562d7c 100644
--- a/src/test/java/com/wildbitsfoundry/etk4j/math/linearalgebra/MatrixDenseTest.java
+++ b/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;