Repository: Hipparchus-Math/hipparchus
Pull request #17 (2016-11-21 10:07:21)

FieldQRdecomposition.
Added the FieldQRDecomposition with its tests. Based completely on the QRDecomposition,
just did the translation.
Added because I had the need on Orekit, I thought it could be useful also for Hipparchus.


Patch:
diff --git a/hipparchus-core/src/main/java/org/hipparchus/linear/FieldQRDecomposition.java b/hipparchus-core/src/main/java/org/hipparchus/linear/FieldQRDecomposition.java
new file mode 100644
index 000000000..e9280c2db
--- /dev/null
+++ b/hipparchus-core/src/main/java/org/hipparchus/linear/FieldQRDecomposition.java
@@ -0,0 +1,495 @@
+/*
+ * 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.
+ */
+
+package org.hipparchus.linear;
+    
+import java.util.Arrays;
+
+import org.hipparchus.RealFieldElement;
+import org.hipparchus.exception.LocalizedCoreFormats;
+import org.hipparchus.exception.MathIllegalArgumentException;
+import org.hipparchus.util.FastMath;
+import org.hipparchus.util.MathArrays;
+
+
+/**
+ * Calculates the QR-decomposition of a field matrix.
+ * <p>The QR-decomposition of a matrix A consists of two matrices Q and R
+ * that satisfy: A = QR, Q is orthogonal (Q<sup>T</sup>Q = I), and R is
+ * upper triangular. If A is m&times;n, Q is m&times;m and R m&times;n.</p>
+ * <p>This class compute the decomposition using Householder reflectors.</p>
+ * <p>For efficiency purposes, the decomposition in packed form is transposed.
+ * This allows inner loop to iterate inside rows, which is much more cache-efficient
+ * in Java.</p>
+ * <p>This class is based on the class {@link QRDecomposition}.</p>
+ *
+ * @see <a href="http://mathworld.wolfram.com/QRDecomposition.html">MathWorld</a>
+ * @see <a href="http://en.wikipedia.org/wiki/QR_decomposition">Wikipedia</a>
+ *
+ */
+public class FieldQRDecomposition<T extends RealFieldElement<T>> {
+    /**
+     * A packed TRANSPOSED representation of the QR decomposition.
+     * <p>The elements BELOW the diagonal are the elements of the UPPER triangular
+     * matrix R, and the rows ABOVE the diagonal are the Householder reflector vectors
+     * from which an explicit form of Q can be recomputed if desired.</p>
+     */
+    private T[][] qrt;
+    /** The diagonal elements of R. */
+    private T[] rDiag;
+    /** Cached value of Q. */
+    private FieldMatrix<T> cachedQ;
+    /** Cached value of QT. */
+    private FieldMatrix<T> cachedQT;
+    /** Cached value of R. */
+    private FieldMatrix<T> cachedR;
+    /** Cached value of H. */
+    private FieldMatrix<T> cachedH;
+    /** Singularity threshold. */
+    private final T threshold;
+
+    /**
+     * Calculates the QR-decomposition of the given matrix.
+     * The singularity threshold defaults to zero.
+     *
+     * @param matrix The matrix to decompose.
+     *
+     * @see #QRDecomposition(FieldMatrix<T>,T)
+     */
+    public FieldQRDecomposition(FieldMatrix<T> matrix) {
+        this(matrix, matrix.getField().getZero());
+    }
+
+    /**
+     * Calculates the QR-decomposition of the given matrix.
+     *
+     * @param matrix The matrix to decompose.
+     * @param threshold Singularity threshold.
+     */
+    public FieldQRDecomposition(FieldMatrix<T> matrix,
+                           T threshold) {
+        this.threshold = threshold;
+
+        final int m = matrix.getRowDimension();
+        final int n = matrix.getColumnDimension();
+        qrt = matrix.transpose().getData();
+        rDiag = MathArrays.buildArray(threshold.getField(),FastMath.min(m, n));
+        cachedQ  = null;
+        cachedQT = null;
+        cachedR  = null;
+        cachedH  = null;
+
+        decompose(qrt);
+
+    }
+
+    /** Decompose matrix.
+     * @param matrix transposed matrix
+     */
+    protected void decompose(T[][] matrix) {
+        for (int minor = 0; minor < FastMath.min(matrix.length, matrix[0].length); minor++) {
+            performHouseholderReflection(minor, matrix);
+        }
+    }
+
+    /** Perform Householder reflection for a minor A(minor, minor) of A.
+     * @param minor minor index
+     * @param matrix transposed matrix
+     */
+    protected void performHouseholderReflection(int minor, T[][] matrix) {
+
+        final T[] qrtMinor = matrix[minor];
+        final T zero = threshold.getField().getZero();
+        /*
+         * Let x be the first column of the minor, and a^2 = |x|^2.
+         * x will be in the positions qr[minor][minor] through qr[m][minor].
+         * The first column of the transformed minor will be (a,0,0,..)'
+         * The sign of a is chosen to be opposite to the sign of the first
+         * component of x. Let's find a:
+         */
+        T xNormSqr = zero;
+        for (int row = minor; row < qrtMinor.length; row++) {
+            final T c = qrtMinor[row];
+            xNormSqr = xNormSqr.add(c.multiply(c));
+        }
+        final T a = (qrtMinor[minor].getReal() > 0) ? xNormSqr.sqrt().negate() : xNormSqr.sqrt();
+        rDiag[minor] = a;
+
+        if (a.getReal() != 0.0) {
+
+            /*
+             * Calculate the normalized reflection vector v and transform
+             * the first column. We know the norm of v beforehand: v = x-ae
+             * so |v|^2 = <x-ae,x-ae> = <x,x>-2a<x,e>+a^2<e,e> =
+             * a^2+a^2-2a<x,e> = 2a*(a - <x,e>).
+             * Here <x, e> is now qr[minor][minor].
+             * v = x-ae is stored in the column at qr:
+             */
+            qrtMinor[minor] = qrtMinor[minor].subtract(a); // now |v|^2 = -2a*(qr[minor][minor])
+
+            /*
+             * Transform the rest of the columns of the minor:
+             * They will be transformed by the matrix H = I-2vv'/|v|^2.
+             * If x is a column vector of the minor, then
+             * Hx = (I-2vv'/|v|^2)x = x-2vv'x/|v|^2 = x - 2<x,v>/|v|^2 v.
+             * Therefore the transformation is easily calculated by
+             * subtracting the column vector (2<x,v>/|v|^2)v from x.
+             *
+             * Let 2<x,v>/|v|^2 = alpha. From above we have
+             * |v|^2 = -2a*(qr[minor][minor]), so
+             * alpha = -<x,v>/(a*qr[minor][minor])
+             */
+            for (int col = minor+1; col < matrix.length; col++) {
+                final T[] qrtCol = matrix[col];
+                T alpha = zero;
+                for (int row = minor; row < qrtCol.length; row++) {
+                    alpha = alpha.subtract(qrtCol[row].multiply(qrtMinor[row]));
+                }
+                alpha = alpha.divide(a.multiply(qrtMinor[minor]));
+
+                // Subtract the column vector alpha*v from x.
+                for (int row = minor; row < qrtCol.length; row++) {
+                    qrtCol[row] = qrtCol[row].subtract(alpha.multiply(qrtMinor[row]));
+                }
+            }
+        }
+    }
+
+
+    /**
+     * Returns the matrix R of the decomposition.
+     * <p>R is an upper-triangular matrix</p>
+     * @return the R matrix
+     */
+    public FieldMatrix<T> getR() {
+
+        if (cachedR == null) {
+
+            // R is supposed to be m x n
+            final int n = qrt.length;
+            final int m = qrt[0].length;
+            T[][] ra = MathArrays.buildArray(threshold.getField(), m, n);
+            // copy the diagonal from rDiag and the upper triangle of qr
+            for (int row = FastMath.min(m, n) - 1; row >= 0; row--) {
+                ra[row][row] = rDiag[row];
+                for (int col = row + 1; col < n; col++) {
+                    ra[row][col] = qrt[col][row];
+                }
+            }
+            cachedR = MatrixUtils.createFieldMatrix(ra);
+        }
+
+        // return the cached matrix
+        return cachedR;
+    }
+
+    /**
+     * Returns the matrix Q of the decomposition.
+     * <p>Q is an orthogonal matrix</p>
+     * @return the Q matrix
+     */
+    public FieldMatrix<T> getQ() {
+        if (cachedQ == null) {
+            cachedQ = getQT().transpose();
+        }
+        return cachedQ;
+    }
+
+    /**
+     * Returns the transpose of the matrix Q of the decomposition.
+     * <p>Q is an orthogonal matrix</p>
+     * @return the transpose of the Q matrix, Q<sup>T</sup>
+     */
+    public FieldMatrix<T> getQT() {
+        if (cachedQT == null) {
+
+            // QT is supposed to be m x m
+            final int n = qrt.length;
+            final int m = qrt[0].length;
+            T[][] qta = MathArrays.buildArray(threshold.getField(), m, m);
+
+            /*
+             * Q = Q1 Q2 ... Q_m, so Q is formed by first constructing Q_m and then
+             * applying the Householder transformations Q_(m-1),Q_(m-2),...,Q1 in
+             * succession to the result
+             */
+            for (int minor = m - 1; minor >= FastMath.min(m, n); minor--) {
+                qta[minor][minor] = threshold.getField().getOne();
+            }
+
+            for (int minor = FastMath.min(m, n)-1; minor >= 0; minor--){
+                final T[] qrtMinor = qrt[minor];
+                qta[minor][minor] = threshold.getField().getOne();
+                if (qrtMinor[minor].getReal() != 0.0) {
+                    for (int col = minor; col < m; col++) {
+                        T alpha = threshold.getField().getZero();
+                        for (int row = minor; row < m; row++) {
+                            alpha = alpha.subtract(qta[col][row].multiply(qrtMinor[row]));
+                        }
+                        alpha = alpha.divide(rDiag[minor].multiply(qrtMinor[minor]));
+
+                        for (int row = minor; row < m; row++) {
+                            qta[col][row] = qta[col][row].add(alpha.negate().multiply(qrtMinor[row]));
+                        }
+                    }
+                }
+            }
+            cachedQT = MatrixUtils.createFieldMatrix(qta);
+        }
+
+        // return the cached matrix
+        return cachedQT;
+    }
+
+    /**
+     * Returns the Householder reflector vectors.
+     * <p>H is a lower trapezoidal matrix whose columns represent
+     * each successive Householder reflector vector. This matrix is used
+     * to compute Q.</p>
+     * @return a matrix containing the Householder reflector vectors
+     */
+    public FieldMatrix<T> getH() {
+        if (cachedH == null) {
+
+            final int n = qrt.length;
+            final int m = qrt[0].length;
+            T[][] ha = MathArrays.buildArray(threshold.getField(), m, n);
+            for (int i = 0; i < m; ++i) {
+                for (int j = 0; j < FastMath.min(i + 1, n); ++j) {
+                    ha[i][j] = qrt[j][i].divide(rDiag[j].negate());
+                }
+            }
+            cachedH = MatrixUtils.createFieldMatrix(ha);
+        }
+
+        // return the cached matrix
+        return cachedH;
+    }
+
+    /**
+     * Get a solver for finding the A &times; X = B solution in least square sense.
+     * <p>
+     * Least Square sense means a solver can be computed for an overdetermined system,
+     * (i.e. a system with more equations than unknowns, which corresponds to a tall A
+     * matrix with more rows than columns). In any case, if the matrix is singular
+     * within the tolerance set at {@link FieldQRDecomposition#QRDecomposition(FieldMatrix<T>,
+     * T) construction}, an error will be triggered when
+     * the {@link DecompositionSolver#solve(RealVector) solve} method will be called.
+     * </p>
+     * @return a solver
+     */
+    public FieldDecompositionSolver<T> getSolver() {
+        return new FieldSolver<T>(qrt, rDiag, threshold);
+    }
+
+    /** Specialized solver. */
+    private static class FieldSolver<T extends RealFieldElement<T>> implements FieldDecompositionSolver<T>{
+        /**
+         * A packed TRANSPOSED representation of the QR decomposition.
+         * <p>The elements BELOW the diagonal are the elements of the UPPER triangular
+         * matrix R, and the rows ABOVE the diagonal are the Householder reflector vectors
+         * from which an explicit form of Q can be recomputed if desired.</p>
+         */
+        private final T[][] qrt;
+        /** The diagonal elements of R. */
+        private final T[] rDiag;
+        /** Singularity threshold. */
+        private final T threshold;
+
+        /**
+         * Build a solver from decomposed matrix.
+         *
+         * @param qrt Packed TRANSPOSED representation of the QR decomposition.
+         * @param rDiag Diagonal elements of R.
+         * @param threshold Singularity threshold.
+         */
+        private FieldSolver(final T[][] qrt,
+                       final T[] rDiag,
+                       final T threshold) {
+            this.qrt   = qrt;
+            this.rDiag = rDiag;
+            this.threshold = threshold;
+        }
+
+        /** {@inheritDoc} */
+        @Override
+        public boolean isNonSingular() {
+            return !checkSingular(rDiag, threshold, false);
+        }
+
+        /** {@inheritDoc} */
+        @Override
+        public FieldVector<T> solve(FieldVector<T> b) {
+            final int n = qrt.length;
+            final int m = qrt[0].length;
+            if (b.getDimension() != m) {
+                throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
+                                                       b.getDimension(), m);
+            }
+            checkSingular(rDiag, threshold, true);
+
+            final T[] x =MathArrays.buildArray(threshold.getField(),n);
+            final T[] y = b.toArray();
+
+            // apply Householder transforms to solve Q.y = b
+            for (int minor = 0; minor < FastMath.min(m, n); minor++) {
+
+                final T[] qrtMinor = qrt[minor];
+                T dotProduct = threshold.getField().getZero();
+                for (int row = minor; row < m; row++) {
+                    dotProduct = dotProduct.add(y[row].multiply(qrtMinor[row]));
+                }
+                dotProduct =  dotProduct.divide(rDiag[minor].multiply(qrtMinor[minor]));
+
+                for (int row = minor; row < m; row++) {
+                    y[row] = y[row].add(dotProduct.multiply(qrtMinor[row]));
+                }
+            }
+
+            // solve triangular system R.x = y
+            for (int row = rDiag.length - 1; row >= 0; --row) {
+                y[row] = y[row].divide(rDiag[row]);
+                final T yRow = y[row];
+                final T[] qrtRow = qrt[row];
+                x[row] = yRow;
+                for (int i = 0; i < row; i++) {
+                    y[i] = y[i].subtract(yRow.multiply(qrtRow[i]));
+                }
+            }
+
+            return new ArrayFieldVector<T>(x, false);
+        }
+
+        /** {@inheritDoc} */
+        @Override
+        public FieldMatrix<T> solve(FieldMatrix<T> b) {
+            final int n = qrt.length;
+            final int m = qrt[0].length;
+            if (b.getRowDimension() != m) {
+                throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
+                                                       b.getRowDimension(), m);
+            }
+            checkSingular(rDiag, threshold, true);
+
+            final int columns        = b.getColumnDimension();
+            final int blockSize      = BlockFieldMatrix.BLOCK_SIZE;
+            final int cBlocks        = (columns + blockSize - 1) / blockSize;
+            final T[][] xBlocks = BlockFieldMatrix.createBlocksLayout(threshold.getField(),n, columns);
+            final T[][] y       = MathArrays.buildArray(threshold.getField(), b.getRowDimension(), blockSize);
+            final T[]   alpha   = MathArrays.buildArray(threshold.getField(), blockSize);
+
+            for (int kBlock = 0; kBlock < cBlocks; ++kBlock) {
+                final int kStart = kBlock * blockSize;
+                final int kEnd   = FastMath.min(kStart + blockSize, columns);
+                final int kWidth = kEnd - kStart;
+
+                // get the right hand side vector
+                b.copySubMatrix(0, m - 1, kStart, kEnd - 1, y);
+
+                // apply Householder transforms to solve Q.y = b
+                for (int minor = 0; minor < FastMath.min(m, n); minor++) {
+                    final T[] qrtMinor = qrt[minor];
+                    final T factor     = rDiag[minor].multiply(qrtMinor[minor]).reciprocal();
+
+                    Arrays.fill(alpha, 0, kWidth, threshold.getField().getZero());
+                    for (int row = minor; row < m; ++row) {
+                        final T   d    = qrtMinor[row];
+                        final T[] yRow = y[row];
+                        for (int k = 0; k < kWidth; ++k) {
+                            alpha[k] = alpha[k].add(d.multiply(yRow[k]));
+                        }
+                    }
+                    
+                    for (int k = 0; k < kWidth; ++k) {
+                        alpha[k] = alpha[k].multiply(factor);
+                    }
+
+                    for (int row = minor; row < m; ++row) {
+                        final T   d    = qrtMinor[row];
+                        final T[] yRow = y[row];
+                        for (int k = 0; k < kWidth; ++k) {
+                            yRow[k] = yRow[k].add(alpha[k].multiply(d));
+                        }
+                    }
+                }
+
+                // solve triangular system R.x = y
+                for (int j = rDiag.length - 1; j >= 0; --j) {
+                    final int      jBlock = j / blockSize;
+                    final int      jStart = jBlock * blockSize;
+                    final T   factor = rDiag[j].reciprocal();
+                    final T[] yJ     = y[j];
+                    final T[] xBlock = xBlocks[jBlock * cBlocks + kBlock];
+                    int index = (j - jStart) * kWidth;
+                    for (int k = 0; k < kWidth; ++k) {
+                        yJ[k]           =yJ[k].multiply(factor);
+                        xBlock[index++] = yJ[k];
+                    }
+
+                    final T[] qrtJ = qrt[j];
+                    for (int i = 0; i < j; ++i) {
+                        final T rIJ  = qrtJ[i];
+                        final T[] yI = y[i];
+                        for (int k = 0; k < kWidth; ++k) {
+                            yI[k] = yI[k].subtract(yJ[k].multiply(rIJ));
+                        }
+                    }
+                }
+            }
+
+            return new BlockFieldMatrix<T>(n, columns, xBlocks, false);
+        }
+
+        /**
+         * {@inheritDoc}
+         * @throws MathIllegalArgumentException if the decomposed matrix is singular.
+         */
+        @Override
+        public FieldMatrix<T> getInverse() {
+            return solve(MatrixUtils.createFieldIdentityMatrix(threshold.getField(), qrt[0].length));
+        }
+
+        /**
+         * Check singularity.
+         *
+         * @param diag Diagonal elements of the R matrix.
+         * @param min Singularity threshold.
+         * @param raise Whether to raise a {@link MathIllegalArgumentException}
+         * if any element of the diagonal fails the check.
+         * @return {@code true} if any element of the diagonal is smaller
+         * or equal to {@code min}.
+         * @throws MathIllegalArgumentException if the matrix is singular and
+         * {@code raise} is {@code true}.
+         */
+        private boolean checkSingular(T[] diag,
+                                             T min,
+                                             boolean raise) {
+            final int len = diag.length;
+            for (int i = 0; i < len; i++) {
+                final T d = diag[i];
+                if (FastMath.abs(d.getReal()) <= min.getReal()) {
+                    if (raise) {
+                        throw new MathIllegalArgumentException(LocalizedCoreFormats.SINGULAR_MATRIX);
+                    } else {
+                        return true;
+                    }
+                }
+            }
+            return false;
+        }
+    }
+}


Test patch:
diff --git a/hipparchus-core/src/test/java/org/hipparchus/linear/FieldQRDecompositionTest.java b/hipparchus-core/src/test/java/org/hipparchus/linear/FieldQRDecompositionTest.java
new file mode 100644
index 000000000..8d92199e7
--- /dev/null
+++ b/hipparchus-core/src/test/java/org/hipparchus/linear/FieldQRDecompositionTest.java
@@ -0,0 +1,438 @@
+/*
+ * 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.
+ */
+
+package org.hipparchus.linear;
+
+import java.util.Random;
+
+import org.hipparchus.Field;
+import org.hipparchus.RealFieldElement;
+import org.hipparchus.analysis.differentiation.DerivativeStructure;
+import org.hipparchus.exception.MathIllegalArgumentException;
+import org.hipparchus.util.Decimal64Field;
+import org.hipparchus.util.MathArrays;
+import org.junit.Assert;
+import org.junit.Test;
+
+
+public class FieldQRDecompositionTest {
+    private double[][] testData3x3NonSingular = {
+            { 12, -51, 4 },
+            { 6, 167, -68 },
+            { -4, 24, -41 }, };
+
+    private double[][] testData3x3Singular = {
+            { 1, 4, 7, },
+            { 2, 5, 8, },
+            { 3, 6, 9, }, };
+
+    private double[][] testData3x4 = {
+            { 12, -51, 4, 1 },
+            { 6, 167, -68, 2 },
+            { -4, 24, -41, 3 }, };
+
+    private double[][] testData4x3 = {
+            { 12, -51, 4, },
+            { 6, 167, -68, },
+            { -4, 24, -41, },
+            { -5, 34, 7, }, };
+
+    DerivativeStructure zero = new DerivativeStructure(1, 1, 0, 0);
+    Field<DerivativeStructure> DSField = zero.getField();
+    private static final double entryTolerance = 10e-16;
+
+    private static final double normTolerance = 10e-14;
+
+    /**Testing if the dimensions are correct.*/
+    @Test
+    public void testDimensions() {
+        doTestDimensions(DSField);   
+        doTestDimensions(Decimal64Field.getInstance());   
+    }
+
+    /**Testing if is impossible to solve QR.*/
+    @Test(expected=MathIllegalArgumentException.class)
+    public void testQRSingular(){
+        QRSingular(DSField);
+        QRSingular(Decimal64Field.getInstance());
+    }
+
+    /**Testing if Q is orthogonal*/
+    @Test
+    public void testQOrthogonal(){
+        QOrthogonal(DSField);
+        QOrthogonal(Decimal64Field.getInstance());
+    }
+
+    /**Testing if A = Q * R*/
+    @Test
+    public void testAEqualQR(){
+        AEqualQR(DSField);
+        AEqualQR(Decimal64Field.getInstance());
+    }
+
+    /**Testing if R is upper triangular.*/
+    @Test
+    public void testRUpperTriangular(){
+        RUpperTriangular(DSField);
+        RUpperTriangular(Decimal64Field.getInstance());
+    }
+
+    /**Testing if H is trapezoidal.*/
+    @Test
+    public void testHTrapezoidal(){
+        HTrapezoidal(DSField);
+        HTrapezoidal(Decimal64Field.getInstance());
+    }
+    /**Testing the values of the matrices.*/
+    @Test
+    public void testMatricesValues(){
+        MatricesValues(DSField);
+        MatricesValues(Decimal64Field.getInstance());
+    }
+
+    /**Testing if there is an error inverting a non invertible matrix.*/
+    @Test(expected=MathIllegalArgumentException.class)
+    public void testNonInvertible(){
+        NonInvertible(DSField);
+        NonInvertible(Decimal64Field.getInstance());
+    }
+    /**Testing to invert a tall and skinny matrix.*/
+    @Test
+    public void testInvertTallSkinny(){
+        InvertTallSkinny(DSField);
+        InvertTallSkinny(Decimal64Field.getInstance());
+    }
+    /**Testing to invert a short and wide matrix.*/
+    @Test
+    public void testInvertShortWide(){
+        InvertShortWide(DSField);
+        InvertShortWide(Decimal64Field.getInstance());
+    }
+
+    private <T extends RealFieldElement<T>> void doTestDimensions(Field<T> field){
+        T[][] data3x3NS= convert(field, testData3x3NonSingular ); 
+        T[][] data3x4= convert(field, testData3x4            );
+        T[][] data4x3= convert(field, testData4x3            );
+        checkDimension(MatrixUtils.createFieldMatrix( data3x3NS));
+
+        checkDimension(MatrixUtils.createFieldMatrix( data4x3));
+
+        checkDimension(MatrixUtils.createFieldMatrix( data3x4));
+
+        Random r = new Random(643895747384642l);
+        int    p = (5 * BlockFieldMatrix.BLOCK_SIZE) / 4;
+        int    q = (7 * BlockFieldMatrix.BLOCK_SIZE) / 4;
+        checkDimension(createTestMatrix(field,r, p, q));
+        checkDimension(createTestMatrix(field, r, q, p));
+
+    }
+
+    private <T extends RealFieldElement<T>> void checkDimension(FieldMatrix<T> m) {
+        int rows = m.getRowDimension();
+        int columns = m.getColumnDimension();
+        FieldQRDecomposition<T> qr = new FieldQRDecomposition<T>(m);
+        Assert.assertEquals(rows,    qr.getQ().getRowDimension());
+        Assert.assertEquals(rows,    qr.getQ().getColumnDimension());
+        Assert.assertEquals(rows,    qr.getR().getRowDimension());
+        Assert.assertEquals(columns, qr.getR().getColumnDimension());
+    }
+
+    private  <T extends RealFieldElement<T>> void AEqualQR(Field<T> field) {
+        T[][] data3x3NS= convert(field, testData3x3NonSingular ); 
+        T[][] data3x3S= convert(field, testData3x3Singular    );
+        T[][] data3x4= convert(field, testData3x4            );
+        T[][] data4x3= convert(field, testData4x3            );
+        checkAEqualQR(MatrixUtils.createFieldMatrix( data3x3NS));
+
+        checkAEqualQR(MatrixUtils.createFieldMatrix( data3x3S));
+
+        checkAEqualQR(MatrixUtils.createFieldMatrix( data3x4));
+
+        checkAEqualQR(MatrixUtils.createFieldMatrix( data4x3));
+
+        Random r = new Random(643895747384642l);
+        int    p = (5 * BlockFieldMatrix.BLOCK_SIZE) / 4;
+        int    q = (7 * BlockFieldMatrix.BLOCK_SIZE) / 4;
+        checkAEqualQR(createTestMatrix(field, r, p, q));
+
+        checkAEqualQR(createTestMatrix(field, r, q, p));
+
+    }
+
+    private  <T extends RealFieldElement<T>> void checkAEqualQR(FieldMatrix<T> m) {
+        FieldQRDecomposition<T> qr = new FieldQRDecomposition<>(m);
+        T norm = norm(qr.getQ().multiply(qr.getR()).subtract(m));
+        Assert.assertEquals(0, norm.getReal(), normTolerance);
+    }
+
+    private  <T extends RealFieldElement<T>> void QOrthogonal(Field<T> field) {
+        T[][] data3x3NS= convert(field, testData3x3NonSingular ); 
+        T[][] data3x3S= convert(field, testData3x3Singular    );
+        T[][] data3x4= convert(field, testData3x4            );
+        T[][] data4x3= convert(field, testData4x3            );
+        checkQOrthogonal(MatrixUtils.createFieldMatrix( data3x3NS));
+
+        checkQOrthogonal(MatrixUtils.createFieldMatrix( data3x3S));
+
+        checkQOrthogonal(MatrixUtils.createFieldMatrix( data3x4));
+
+        checkQOrthogonal(MatrixUtils.createFieldMatrix( data4x3));
+
+        Random r = new Random(643895747384642l);
+        int    p = (5 * BlockFieldMatrix.BLOCK_SIZE) / 4;
+        int    q = (7 * BlockFieldMatrix.BLOCK_SIZE) / 4;
+        checkQOrthogonal(createTestMatrix(field, r, p, q));
+
+        checkQOrthogonal(createTestMatrix(field, r, q, p));
+
+    }
+
+    private  <T extends RealFieldElement<T>> void checkQOrthogonal(FieldMatrix<T> m) {
+        FieldQRDecomposition<T> qr = new FieldQRDecomposition<T>(m);
+        FieldMatrix<T> eye = MatrixUtils.createFieldIdentityMatrix(m.getField(),m.getRowDimension());
+        T norm = norm(qr.getQT().multiply(qr.getQ()).subtract(eye));
+        Assert.assertEquals(0, norm.getReal(), normTolerance);
+    }
+
+    private  <T extends RealFieldElement<T>> void RUpperTriangular(Field<T> field) {
+        T[][] data3x3NS= convert(field, testData3x3NonSingular ); 
+        T[][] data3x3S= convert(field, testData3x3Singular    );
+        T[][] data3x4= convert(field, testData3x4            );
+        T[][] data4x3= convert(field, testData4x3            );
+        FieldMatrix<T> matrix = MatrixUtils.createFieldMatrix( data3x3NS);
+        checkUpperTriangular(new FieldQRDecomposition<T>(matrix).getR());
+
+        matrix = MatrixUtils.createFieldMatrix( data3x3S);
+        checkUpperTriangular(new FieldQRDecomposition<T>(matrix).getR());
+
+        matrix = MatrixUtils.createFieldMatrix( data3x4);
+        checkUpperTriangular(new FieldQRDecomposition<T>(matrix).getR());
+
+        matrix = MatrixUtils.createFieldMatrix( data4x3);
+        checkUpperTriangular(new FieldQRDecomposition<T>(matrix).getR());
+
+        Random r = new Random(643895747384642l);
+        int    p = (5 * BlockFieldMatrix.BLOCK_SIZE) / 4;
+        int    q = (7 * BlockFieldMatrix.BLOCK_SIZE) / 4;
+        matrix = createTestMatrix(field, r, p, q);
+        checkUpperTriangular(new FieldQRDecomposition<T>(matrix).getR());
+
+        matrix = createTestMatrix(field, r, p, q);
+        checkUpperTriangular(new FieldQRDecomposition<T>(matrix).getR());
+
+    }
+
+    private  <T extends RealFieldElement<T>> void checkUpperTriangular(FieldMatrix<T> m) {
+        m.walkInOptimizedOrder(new DefaultFieldMatrixPreservingVisitor<T>(m.getField().getZero()) {
+            @Override
+            public void visit(int row, int column, T value) {
+                if (column < row) {
+                    Assert.assertEquals(0.0, value.getReal(), entryTolerance);
+                }
+            }
+        });
+    }
+
+    private <T extends RealFieldElement<T>> void HTrapezoidal(Field<T> field) {
+        T[][] data3x3NS= convert(field, testData3x3NonSingular ); 
+        T[][] data3x3S= convert(field, testData3x3Singular    );
+        T[][] data3x4= convert(field, testData3x4            );
+        T[][] data4x3= convert(field, testData4x3            );
+        FieldMatrix<T> matrix = MatrixUtils.createFieldMatrix( data3x3NS);
+        checkTrapezoidal(new FieldQRDecomposition<T>(matrix).getH());
+
+        matrix = MatrixUtils.createFieldMatrix( data3x3S);
+        checkTrapezoidal(new FieldQRDecomposition<T>(matrix).getH());
+
+        matrix = MatrixUtils.createFieldMatrix( data3x4);
+        checkTrapezoidal(new FieldQRDecomposition<T>(matrix).getH());
+
+        matrix = MatrixUtils.createFieldMatrix( data4x3);
+        checkTrapezoidal(new FieldQRDecomposition<T>(matrix).getH());
+
+        Random r = new Random(643895747384642l);
+        int    p = (5 * BlockFieldMatrix.BLOCK_SIZE) / 4;
+        int    q = (7 * BlockFieldMatrix.BLOCK_SIZE) / 4;
+        matrix = createTestMatrix(field, r, p, q);
+        checkTrapezoidal(new FieldQRDecomposition<T>(matrix).getH());
+
+        matrix = createTestMatrix(field, r, p, q);
+        checkTrapezoidal(new FieldQRDecomposition<T>(matrix).getH());
+
+    }
+
+    private  <T extends RealFieldElement<T>> void checkTrapezoidal(FieldMatrix<T> m) {
+        m.walkInOptimizedOrder(new DefaultFieldMatrixPreservingVisitor<T>(m.getField().getZero()) {
+            @Override
+            public void visit(int row, int column, T value) {
+                if (column > row) {
+                    Assert.assertEquals(0.0, value.getReal(), entryTolerance);
+                }
+            }
+        });
+    }
+    
+    private  <T extends RealFieldElement<T>> void MatricesValues(Field<T> field) {
+        T[][] data3x3NS= convert(field, testData3x3NonSingular ); 
+        FieldQRDecomposition<T> qr =
+            new FieldQRDecomposition<T>(MatrixUtils.createFieldMatrix( data3x3NS));
+        FieldMatrix<T> qRef = MatrixUtils.createFieldMatrix( convert(field,new double[][] {
+                { -12.0 / 14.0,   69.0 / 175.0,  -58.0 / 175.0 },
+                {  -6.0 / 14.0, -158.0 / 175.0,    6.0 / 175.0 },
+                {   4.0 / 14.0,  -30.0 / 175.0, -165.0 / 175.0 }
+        }));
+        FieldMatrix<T> rRef = MatrixUtils.createFieldMatrix( convert(field, new double[][] {
+                { -14.0,  -21.0, 14.0 },
+                {   0.0, -175.0, 70.0 },
+                {   0.0,    0.0, 35.0 }
+        }));
+        FieldMatrix<T> hRef = MatrixUtils.createFieldMatrix( convert(field,new double[][] {
+                { 26.0 / 14.0, 0.0, 0.0 },
+                {  6.0 / 14.0, 648.0 / 325.0, 0.0 },
+                { -4.0 / 14.0,  36.0 / 325.0, 2.0 }
+        }));
+
+        // check values against known references
+        FieldMatrix<T> q = qr.getQ();
+        Assert.assertEquals(0, norm(q.subtract(qRef)).getReal(), 1.0e-13);
+        FieldMatrix<T> qT = qr.getQT();
+        Assert.assertEquals(0, norm(qT.subtract(qRef.transpose())).getReal(), 1.0e-13);
+        FieldMatrix<T> r = qr.getR();
+        Assert.assertEquals(0, norm(r.subtract(rRef)).getReal(), 1.0e-13);
+        FieldMatrix<T> h = qr.getH();
+        Assert.assertEquals(0, norm(h.subtract(hRef)).getReal(), 1.0e-13);
+
+        // check the same cached instance is returned the second time
+        Assert.assertTrue(q == qr.getQ());
+        Assert.assertTrue(r == qr.getR());
+        Assert.assertTrue(h == qr.getH());
+
+    }
+
+    private  <T extends RealFieldElement<T>> void NonInvertible(Field<T> field) {
+        T[][] data3x3S= convert(field, testData3x3Singular    );
+        FieldQRDecomposition<T> qr =
+            new FieldQRDecomposition<T>(MatrixUtils.createFieldMatrix( data3x3S));
+        qr.getSolver().getInverse();
+    }
+
+    private  <T extends RealFieldElement<T>> void InvertTallSkinny(Field<T> field) {
+        T[][] data4x3= convert(field, testData4x3            );
+        FieldMatrix<T> a     = MatrixUtils.createFieldMatrix(data4x3);
+        FieldMatrix<T> pinv  = new FieldQRDecomposition<T>(a).getSolver().getInverse();
+        Assert.assertEquals(0, norm(pinv.multiply(a).subtract(MatrixUtils.createFieldIdentityMatrix(field, 3))).getReal(), 1.0e-6);
+    }
+
+    private  <T extends RealFieldElement<T>> void InvertShortWide(Field<T> field) {
+        T[][] data3x4= convert(field, testData3x4            );
+        FieldMatrix<T> a = MatrixUtils.createFieldMatrix( data3x4);
+        FieldMatrix<T> pinv  = new FieldQRDecomposition<T>(a).getSolver().getInverse();
+        Assert.assertEquals(0,norm( a.multiply(pinv).subtract(MatrixUtils.createFieldIdentityMatrix(field, 3))).getReal(), 1.0e-6);
+        Assert.assertEquals(0,norm( pinv.multiply(a).getSubMatrix(0, 2, 0, 2).subtract(MatrixUtils.createFieldIdentityMatrix(field, 3))).getReal(), 1.0e-6);
+    }
+
+    private  <T extends RealFieldElement<T>> FieldMatrix<T> createTestMatrix(Field<T> field, final Random r, final int rows, final int columns) {
+        FieldMatrix<T> m = MatrixUtils.createFieldMatrix(field, rows, columns);
+        m.walkInOptimizedOrder(new DefaultFieldMatrixChangingVisitor<T>(field.getOne()){
+            @Override
+            public T visit(int row, int column,T value) {
+                return field.getZero().add(2.0 * r.nextDouble() - 1.0);
+            }
+        });
+        return m;
+    }
+
+    private  <T extends RealFieldElement<T>> void QRSingular(Field<T> field) {
+        final FieldMatrix<T> a = MatrixUtils.createFieldMatrix(convert( field, new double[][] {
+            { 1, 6, 4 }, { 2, 4, -1 }, { -1, 2, 5 }
+        }));
+        T[] vv = MathArrays.buildArray(field,3);
+        vv[0] = field.getZero().add(5);
+        vv[1] = field.getZero().add(6);
+        vv[2] = field.getZero().add(1);
+        
+        final FieldVector<T> b = new ArrayFieldVector<T>(field, vv);
+        new FieldQRDecomposition<T>(a, field.getZero().add(1.0e-15)).getSolver().solve(b);
+    }
+    
+    private <T extends RealFieldElement<T>> T norm(FieldMatrix<T> FM ){
+        return walkInColumnOrder(FM, new FieldMatrixPreservingVisitor<T>() {
+
+            /** Last row index. */
+            private double endRow;
+
+            /** Sum of absolute values on one column. */
+            private T columnSum;
+
+            /** Maximal sum across all columns. */
+            private T maxColSum;
+
+            /** {@inheritDoc} */
+            @Override
+            public void start(final int rows, final int columns,
+                              final int startRow, final int endRow,
+                              final int startColumn, final int endColumn) {
+                this.endRow = endRow;
+                columnSum   = FM.getField().getZero();
+                maxColSum   = FM.getField().getZero();
+            }
+
+            /** {@inheritDoc} */
+            @Override
+            public void visit(final int row, final int column, final T value) {
+                columnSum = columnSum.add(value).abs();
+                if (row == endRow) {
+                    maxColSum = (maxColSum.getReal() > columnSum.getReal()) ? maxColSum : columnSum ;
+                    columnSum = FM.getField().getZero();
+                }
+            }
+
+            /** {@inheritDoc} */
+            @Override
+            public T end() {
+                return maxColSum;
+            }
+
+        });
+    }
+    
+    private <T extends RealFieldElement<T>> T walkInColumnOrder(FieldMatrix<T> FM, FieldMatrixPreservingVisitor<T> visitor) {
+        final int rows    = FM.getRowDimension();
+        final int columns = FM.getColumnDimension();
+        visitor.start(rows, columns, 0, rows - 1, 0, columns - 1);
+        for (int column = 0; column < columns; ++column) {
+            for (int row = 0; row < rows; ++row) {
+                visitor.visit(row, column, FM.getEntry(row, column));
+            }
+        }
+        return visitor.end();
+    }
+    
+    private <T extends RealFieldElement<T>> T[][] convert(Field<T> field, double[][] value){
+        T[][] res = MathArrays.buildArray(field, value.length, value[0].length);
+        for (int ii = 0; ii < (value.length); ii++){
+            for (int jj = 0; jj < (value[0].length); jj++){
+                res[ii][jj] = field.getZero().add(value[ii][jj]);
+            }
+        }
+        return res;
+        
+    }
+
+
+}