org.apache.commons.math.linear.MatrixUtils Java Examples

/** test rank */
public void testRank() {
    DecompositionSolver solver =
        new QRDecompositionImpl(MatrixUtils.createRealMatrix(testData3x3NonSingular)).getSolver();
    assertTrue(solver.isNonSingular());

    solver = new QRDecompositionImpl(MatrixUtils.createRealMatrix(testData3x3Singular)).getSolver();
    assertFalse(solver.isNonSingular());

    solver = new QRDecompositionImpl(MatrixUtils.createRealMatrix(testData3x4)).getSolver();
    assertTrue(solver.isNonSingular());

    solver = new QRDecompositionImpl(MatrixUtils.createRealMatrix(testData4x3)).getSolver();
    assertTrue(solver.isNonSingular());

}
 

/** test that R is upper triangular */
public void testRUpperTriangular() {
    RealMatrix matrix = MatrixUtils.createRealMatrix(testData3x3NonSingular);
    checkUpperTriangular(new QRDecompositionImpl(matrix).getR());

    matrix = MatrixUtils.createRealMatrix(testData3x3Singular);
    checkUpperTriangular(new QRDecompositionImpl(matrix).getR());

    matrix = MatrixUtils.createRealMatrix(testData3x4);
    checkUpperTriangular(new QRDecompositionImpl(matrix).getR());

    matrix = MatrixUtils.createRealMatrix(testData4x3);
    checkUpperTriangular(new QRDecompositionImpl(matrix).getR());

    Random r = new Random(643895747384642l);
    int    p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
    int    q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
    matrix = createTestMatrix(r, p, q);
    checkUpperTriangular(new QRDecompositionImpl(matrix).getR());

    matrix = createTestMatrix(r, p, q);
    checkUpperTriangular(new QRDecompositionImpl(matrix).getR());

}
 

public void testMath226()
    throws DimensionMismatchException, NotPositiveDefiniteMatrixException {
    double[] mean = { 1, 1, 10, 1 };
    double[][] cov = {
            { 1, 3, 2, 6 },
            { 3, 13, 16, 2 },
            { 2, 16, 38, -1 },
            { 6, 2, -1, 197 }
    };
    RealMatrix covRM = MatrixUtils.createRealMatrix(cov);
    JDKRandomGenerator jg = new JDKRandomGenerator();
    jg.setSeed(5322145245211l);
    NormalizedRandomGenerator rg = new GaussianRandomGenerator(jg);
    CorrelatedRandomVectorGenerator sg =
        new CorrelatedRandomVectorGenerator(mean, covRM, 0.00001, rg);

    for (int i = 0; i < 10; i++) {
        double[] generated = sg.nextVector();
        assertTrue(Math.abs(generated[0] - 1) > 0.1);
    }

}
 

private RealMatrix createTestMatrix(final Random r, final int rows, final int columns) {
    RealMatrix m = MatrixUtils.createRealMatrix(rows, columns);
    m.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor(){
        @Override
        public double visit(int row, int column, double value)
            throws MatrixVisitorException {
            return 2.0 * r.nextDouble() - 1.0;
        }
    });
    return m;
}
 

/** test non-symmetric matrix */
@Test(expected = NotSymmetricMatrixException.class)
public void testNotSymmetricMatrixException() throws MathException {
    double[][] changed = testData.clone();
    changed[0][changed[0].length - 1] += 1.0e-5;
    new CholeskyDecompositionImpl(MatrixUtils.createRealMatrix(changed));
}
 

/** test that V is orthogonal */
public void testVOrthogonal() {
    RealMatrix v = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN).getV();
    RealMatrix vTv = v.transpose().multiply(v);
    RealMatrix id  = MatrixUtils.createRealIdentityMatrix(vTv.getRowDimension());
    assertEquals(0, vTv.subtract(id).getNorm(), 2.0e-13);
}
 

/** test dimensions */
public void testDimensions() {
    checkDimension(MatrixUtils.createRealMatrix(testData3x3NonSingular));

    checkDimension(MatrixUtils.createRealMatrix(testData4x3));

    checkDimension(MatrixUtils.createRealMatrix(testData3x4));

    Random r = new Random(643895747384642l);
    int    p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
    int    q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
    checkDimension(createTestMatrix(r, p, q));
    checkDimension(createTestMatrix(r, q, p));

}
 

/** test dimensions */
@Test
public void testDimensions() throws MathException {
    CholeskyDecomposition llt =
        new CholeskyDecompositionImpl(MatrixUtils.createRealMatrix(testData));
    assertEquals(testData.length, llt.getL().getRowDimension());
    assertEquals(testData.length, llt.getL().getColumnDimension());
    assertEquals(testData.length, llt.getLT().getRowDimension());
    assertEquals(testData.length, llt.getLT().getColumnDimension());
}
 

/** test matrices values */
public void testMatricesValues1() {
   LUDecomposition lu =
        new LUDecompositionImpl(MatrixUtils.createRealMatrix(testData));
    RealMatrix lRef = MatrixUtils.createRealMatrix(new double[][] {
            { 1.0, 0.0, 0.0 },
            { 0.5, 1.0, 0.0 },
            { 0.5, 0.2, 1.0 }
    });
    RealMatrix uRef = MatrixUtils.createRealMatrix(new double[][] {
            { 2.0,  5.0, 3.0 },
            { 0.0, -2.5, 6.5 },
            { 0.0,  0.0, 0.2 }
    });
    RealMatrix pRef = MatrixUtils.createRealMatrix(new double[][] {
            { 0.0, 1.0, 0.0 },
            { 0.0, 0.0, 1.0 },
            { 1.0, 0.0, 0.0 }
    });
    int[] pivotRef = { 1, 2, 0 };

    // check values against known references
    RealMatrix l = lu.getL();
    assertEquals(0, l.subtract(lRef).getNorm(), 1.0e-13);
    RealMatrix u = lu.getU();
    assertEquals(0, u.subtract(uRef).getNorm(), 1.0e-13);
    RealMatrix p = lu.getP();
    assertEquals(0, p.subtract(pRef).getNorm(), 1.0e-13);
    int[] pivot = lu.getPivot();
    for (int i = 0; i < pivotRef.length; ++i) {
        assertEquals(pivotRef[i], pivot[i]);
    }

    // check the same cached instance is returned the second time
    assertTrue(l == lu.getL());
    assertTrue(u == lu.getU());
    assertTrue(p == lu.getP());

}
 

public static RealMatrix createDiagonalMatrix(final double[] diagonal,
                                              final int rows, final int columns) {
    final double[][] dData = new double[rows][columns];
    for (int i = 0; i < Math.min(rows, columns); ++i) {
        dData[i][i] = diagonal[i];
    }
    return MatrixUtils.createRealMatrix(dData);
}
 

/**
 * Correct the current state estimate with an actual measurement.
 *
 * @param z
 *            the measurement vector
 * @throws DimensionMismatchException
 *             if the dimension of the measurement vector does not fit
 * @throws org.apache.commons.math.linear.SingularMatrixException
 *             if the covariance matrix could not be inverted
 */
public void correct(final RealVector z) {
    // sanity checks
    MathUtils.checkNotNull(z);
    if (z.getDimension() != measurementMatrix.getRowDimension()) {
        throw new DimensionMismatchException(z.getDimension(),
                                             measurementMatrix.getRowDimension());
    }

    // S = H * P(k) - * H' + R
    RealMatrix s = measurementMatrix.multiply(errorCovariance)
        .multiply(measurementMatrixT)
        .add(measurementModel.getMeasurementNoise());

    // invert S
    // as the error covariance matrix is a symmetric positive
    // semi-definite matrix, we can use the cholesky decomposition
    DecompositionSolver solver = new CholeskyDecompositionImpl(s).getSolver();
    RealMatrix invertedS = solver.getInverse();

    // Inn = z(k) - H * xHat(k)-
    RealVector innovation = z.subtract(measurementMatrix.operate(stateEstimation));

    // calculate gain matrix
    // K(k) = P(k)- * H' * (H * P(k)- * H' + R)^-1
    // K(k) = P(k)- * H' * S^-1
    RealMatrix kalmanGain = errorCovariance.multiply(measurementMatrixT).multiply(invertedS);

    // update estimate with measurement z(k)
    // xHat(k) = xHat(k)- + K * Inn
    stateEstimation = stateEstimation.add(kalmanGain.operate(innovation));

    // update covariance of prediction error
    // P(k) = (I - K * H) * P(k)-
    RealMatrix identity = MatrixUtils.createRealIdentityMatrix(kalmanGain.getRowDimension());
    errorCovariance = identity.subtract(kalmanGain.multiply(measurementMatrix)).multiply(errorCovariance);
}
 

/** test that L is lower triangular with unit diagonal */
public void testLLowerTriangular() {
    RealMatrix matrix = MatrixUtils.createRealMatrix(testData);
    RealMatrix l = new LUDecompositionImpl(matrix).getL();
    for (int i = 0; i < l.getRowDimension(); i++) {
        assertEquals(l.getEntry(i, i), 1, entryTolerance);
        for (int j = i + 1; j < l.getColumnDimension(); j++) {
            assertEquals(l.getEntry(i, j), 0, entryTolerance);
        }
    }
}
 

/** test matrices values */
public void testMatricesValues1() {
   SingularValueDecomposition svd =
        new SingularValueDecompositionImpl(MatrixUtils.createRealMatrix(testSquare));
    RealMatrix uRef = MatrixUtils.createRealMatrix(new double[][] {
            { 3.0 / 5.0, -4.0 / 5.0 },
            { 4.0 / 5.0,  3.0 / 5.0 }
    });
    RealMatrix sRef = MatrixUtils.createRealMatrix(new double[][] {
            { 3.0, 0.0 },
            { 0.0, 1.0 }
    });
    RealMatrix vRef = MatrixUtils.createRealMatrix(new double[][] {
            { 4.0 / 5.0,  3.0 / 5.0 },
            { 3.0 / 5.0, -4.0 / 5.0 }
    });

    // check values against known references
    RealMatrix u = svd.getU();
    assertEquals(0, u.subtract(uRef).getNorm(), normTolerance);
    RealMatrix s = svd.getS();
    assertEquals(0, s.subtract(sRef).getNorm(), normTolerance);
    RealMatrix v = svd.getV();
    assertEquals(0, v.subtract(vRef).getNorm(), normTolerance);

    // check the same cached instance is returned the second time
    assertTrue(u == svd.getU());
    assertTrue(s == svd.getS());
    assertTrue(v == svd.getV());

}
 

/** test that L is lower triangular */
@Test
public void testLLowerTriangular() throws MathException {
    RealMatrix matrix = MatrixUtils.createRealMatrix(testData);
    RealMatrix l = new CholeskyDecompositionImpl(matrix).getL();
    for (int i = 0; i < l.getRowDimension(); i++) {
        for (int j = i + 1; j < l.getColumnDimension(); j++) {
            assertEquals(0.0, l.getEntry(i, j), 0.0);
        }
    }
}
 

/** test matrices values */
public void testMatricesValues2() {

    RealMatrix uRef = MatrixUtils.createRealMatrix(new double[][] {
        {  0.0 / 5.0,  3.0 / 5.0,  0.0 / 5.0 },
        { -4.0 / 5.0,  0.0 / 5.0, -3.0 / 5.0 },
        {  0.0 / 5.0,  4.0 / 5.0,  0.0 / 5.0 },
        { -3.0 / 5.0,  0.0 / 5.0,  4.0 / 5.0 }
    });
    RealMatrix sRef = MatrixUtils.createRealMatrix(new double[][] {
        { 4.0, 0.0, 0.0 },
        { 0.0, 3.0, 0.0 },
        { 0.0, 0.0, 2.0 }
    });
    RealMatrix vRef = MatrixUtils.createRealMatrix(new double[][] {
        {  80.0 / 125.0,  -60.0 / 125.0, 75.0 / 125.0 },
        {  24.0 / 125.0,  107.0 / 125.0, 60.0 / 125.0 },
        { -93.0 / 125.0,  -24.0 / 125.0, 80.0 / 125.0 }
    });

    // check values against known references
    SingularValueDecomposition svd =
        new SingularValueDecompositionImpl(MatrixUtils.createRealMatrix(testNonSquare));
    RealMatrix u = svd.getU();
    assertEquals(0, u.subtract(uRef).getNorm(), normTolerance);
    RealMatrix s = svd.getS();
    assertEquals(0, s.subtract(sRef).getNorm(), normTolerance);
    RealMatrix v = svd.getV();
    assertEquals(0, v.subtract(vRef).getNorm(), normTolerance);

    // check the same cached instance is returned the second time
    assertTrue(u == svd.getU());
    assertTrue(s == svd.getS());
    assertTrue(v == svd.getV());

}
 

@Test
public void testMatricesValues() {
   BiDiagonalTransformer transformer =
        new BiDiagonalTransformer(MatrixUtils.createRealMatrix(testSquare));
   final double s17 = FastMath.sqrt(17.0);
    RealMatrix uRef = MatrixUtils.createRealMatrix(new double[][] {
            {  -8 / (5 * s17), 19 / (5 * s17) },
            { -19 / (5 * s17), -8 / (5 * s17) }
    });
    RealMatrix bRef = MatrixUtils.createRealMatrix(new double[][] {
            { -3 * s17 / 5, 32 * s17 / 85 },
            {      0.0,     -5 * s17 / 17 }
    });
    RealMatrix vRef = MatrixUtils.createRealMatrix(new double[][] {
            { 1.0,  0.0 },
            { 0.0, -1.0 }
    });

    // check values against known references
    RealMatrix u = transformer.getU();
    Assert.assertEquals(0, u.subtract(uRef).getNorm(), 1.0e-14);
    RealMatrix b = transformer.getB();
    Assert.assertEquals(0, b.subtract(bRef).getNorm(), 1.0e-14);
    RealMatrix v = transformer.getV();
    Assert.assertEquals(0, v.subtract(vRef).getNorm(), 1.0e-14);

    // check the same cached instance is returned the second time
    Assert.assertTrue(u == transformer.getU());
    Assert.assertTrue(b == transformer.getB());
    Assert.assertTrue(v == transformer.getV());

}
 

private void checkMatricesValues(double[][] matrix, double[][] qRef,
                                 double[] mainDiagnonal,
                                 double[] secondaryDiagonal) {
    TriDiagonalTransformer transformer =
        new TriDiagonalTransformer(MatrixUtils.createRealMatrix(matrix));

    // check values against known references
    RealMatrix q = transformer.getQ();
    assertEquals(0, q.subtract(MatrixUtils.createRealMatrix(qRef)).getNorm(), 1.0e-14);

    RealMatrix t = transformer.getT();
    double[][] tData = new double[mainDiagnonal.length][mainDiagnonal.length];
    for (int i = 0; i < mainDiagnonal.length; ++i) {
        tData[i][i] = mainDiagnonal[i];
        if (i > 0) {
            tData[i][i - 1] = secondaryDiagonal[i - 1];
        }
        if (i < secondaryDiagonal.length) {
            tData[i][i + 1] = secondaryDiagonal[i];
        }
    }
    assertEquals(0, t.subtract(MatrixUtils.createRealMatrix(tData)).getNorm(), 1.0e-14);

    // check the same cached instance is returned the second time
    assertTrue(q == transformer.getQ());
    assertTrue(t == transformer.getT());

}
 

/** test invertible matrix */
public void testInvertible() {
    Random r = new Random(9994100315209l);
    RealMatrix m =
        EigenDecompositionImplTest.createTestMatrix(r, new double[] { 1.0, 0.5, -1.0, -2.0, -3.0 });
    DecompositionSolver es = new EigenDecompositionImpl(m, MathUtils.SAFE_MIN).getSolver();
    assertTrue(es.isNonSingular());
    RealMatrix inverse = es.getInverse();
    RealMatrix error =
        m.multiply(inverse).subtract(MatrixUtils.createRealIdentityMatrix(m.getRowDimension()));
    assertEquals(0, error.getNorm(), 4.0e-15);
}
 

/** test that L is lower triangular */
@Test
public void testLLowerTriangular() throws MathException {
    RealMatrix matrix = MatrixUtils.createRealMatrix(testData);
    RealMatrix l = new CholeskyDecompositionImpl(matrix).getL();
    for (int i = 0; i < l.getRowDimension(); i++) {
        for (int j = i + 1; j < l.getColumnDimension(); j++) {
            assertEquals(0.0, l.getEntry(i, j), 0.0);
        }
    }
}
 

/** test invertible matrix */
public void testInvertible() {
    Random r = new Random(9994100315209l);
    RealMatrix m =
        EigenDecompositionImplTest.createTestMatrix(r, new double[] { 1.0, 0.5, -1.0, -2.0, -3.0 });
    DecompositionSolver es = new EigenDecompositionImpl(m, MathUtils.SAFE_MIN).getSolver();
    assertTrue(es.isNonSingular());
    RealMatrix inverse = es.getInverse();
    RealMatrix error =
        m.multiply(inverse).subtract(MatrixUtils.createRealIdentityMatrix(m.getRowDimension()));
    assertEquals(0, error.getNorm(), 4.0e-15);
}
 

/** test invertible matrix */
public void testInvertible() {
    Random r = new Random(9994100315209l);
    RealMatrix m =
        EigenDecompositionImplTest.createTestMatrix(r, new double[] { 1.0, 0.5, -1.0, -2.0, -3.0 });
    DecompositionSolver es = new EigenDecompositionImpl(m, MathUtils.SAFE_MIN).getSolver();
    assertTrue(es.isNonSingular());
    RealMatrix inverse = es.getInverse();
    RealMatrix error =
        m.multiply(inverse).subtract(MatrixUtils.createRealIdentityMatrix(m.getRowDimension()));
    assertEquals(0, error.getNorm(), 4.0e-15);
}
 

public void testNonSquare() {
    try {
        new TriDiagonalTransformer(MatrixUtils.createRealMatrix(new double[3][2]));
        fail("an exception should have been thrown");
    } catch (InvalidMatrixException ime) {
        // expected behavior
    } catch (Exception e) {
        fail("wrong exception caught");
    }
}
 

/** test matrices values */
public void testMatricesValues1() {
   LUDecomposition lu =
        new LUDecompositionImpl(MatrixUtils.createRealMatrix(testData));
    RealMatrix lRef = MatrixUtils.createRealMatrix(new double[][] {
            { 1.0, 0.0, 0.0 },
            { 0.5, 1.0, 0.0 },
            { 0.5, 0.2, 1.0 }
    });
    RealMatrix uRef = MatrixUtils.createRealMatrix(new double[][] {
            { 2.0,  5.0, 3.0 },
            { 0.0, -2.5, 6.5 },
            { 0.0,  0.0, 0.2 }
    });
    RealMatrix pRef = MatrixUtils.createRealMatrix(new double[][] {
            { 0.0, 1.0, 0.0 },
            { 0.0, 0.0, 1.0 },
            { 1.0, 0.0, 0.0 }
    });
    int[] pivotRef = { 1, 2, 0 };

    // check values against known references
    RealMatrix l = lu.getL();
    assertEquals(0, l.subtract(lRef).getNorm(), 1.0e-13);
    RealMatrix u = lu.getU();
    assertEquals(0, u.subtract(uRef).getNorm(), 1.0e-13);
    RealMatrix p = lu.getP();
    assertEquals(0, p.subtract(pRef).getNorm(), 1.0e-13);
    int[] pivot = lu.getPivot();
    for (int i = 0; i < pivotRef.length; ++i) {
        assertEquals(pivotRef[i], pivot[i]);
    }

    // check the same cached instance is returned the second time
    assertTrue(l == lu.getL());
    assertTrue(u == lu.getU());
    assertTrue(p == lu.getP());
    
}
 

/** test dimensions */
@Test
public void testDimensions() throws MathException {
    CholeskyDecomposition llt =
        new CholeskyDecompositionImpl(MatrixUtils.createRealMatrix(testData));
    assertEquals(testData.length, llt.getL().getRowDimension());
    assertEquals(testData.length, llt.getL().getColumnDimension());
    assertEquals(testData.length, llt.getLT().getRowDimension());
    assertEquals(testData.length, llt.getLT().getColumnDimension());
}
 

/**
 * Get the covariance matrix of unbound estimated parameters.
 * @param problem estimation problem
 * @return covariance matrix
 * @exception EstimationException if the covariance matrix
 * cannot be computed (singular problem)
 */
public double[][] getCovariances(EstimationProblem problem)
  throws EstimationException {

    // set up the jacobian
    updateJacobian();

    // compute transpose(J).J, avoiding building big intermediate matrices
    final int n = problem.getMeasurements().length;
    final int m = problem.getUnboundParameters().length;
    final int max  = m * n;
    double[][] jTj = new double[m][m];
    for (int i = 0; i < m; ++i) {
        for (int j = i; j < m; ++j) {
            double sum = 0;
            for (int k = 0; k < max; k += m) {
                sum += jacobian[k + i] * jacobian[k + j];
            }
            jTj[i][j] = sum;
            jTj[j][i] = sum;
        }
    }

    try {
        // compute the covariances matrix
        RealMatrix inverse =
            new LUDecompositionImpl(MatrixUtils.createRealMatrix(jTj)).getSolver().getInverse();
        return inverse.getData();
    } catch (InvalidMatrixException ime) {
        throw new EstimationException("unable to compute covariances: singular problem");
    }

}
 

@Test
public void testMatricesValues() {
   BiDiagonalTransformer transformer =
        new BiDiagonalTransformer(MatrixUtils.createRealMatrix(testSquare));
   final double s17 = FastMath.sqrt(17.0);
    RealMatrix uRef = MatrixUtils.createRealMatrix(new double[][] {
            {  -8 / (5 * s17), 19 / (5 * s17) },
            { -19 / (5 * s17), -8 / (5 * s17) }
    });
    RealMatrix bRef = MatrixUtils.createRealMatrix(new double[][] {
            { -3 * s17 / 5, 32 * s17 / 85 },
            {      0.0,     -5 * s17 / 17 }
    });
    RealMatrix vRef = MatrixUtils.createRealMatrix(new double[][] {
            { 1.0,  0.0 },
            { 0.0, -1.0 }
    });

    // check values against known references
    RealMatrix u = transformer.getU();
    Assert.assertEquals(0, u.subtract(uRef).getNorm(), 1.0e-14);
    RealMatrix b = transformer.getB();
    Assert.assertEquals(0, b.subtract(bRef).getNorm(), 1.0e-14);
    RealMatrix v = transformer.getV();
    Assert.assertEquals(0, v.subtract(vRef).getNorm(), 1.0e-14);

    // check the same cached instance is returned the second time
    Assert.assertTrue(u == transformer.getU());
    Assert.assertTrue(b == transformer.getB());
    Assert.assertTrue(v == transformer.getV());

}
 

private RealMatrix createTestMatrix(final Random r, final int rows, final int columns) {
    RealMatrix m = MatrixUtils.createRealMatrix(rows, columns);
    m.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor(){
        @Override
        public double visit(int row, int column, double value)
            throws MatrixVisitorException {
            return 2.0 * r.nextDouble() - 1.0;
        }
    });
    return m;
}
 

/** test matrices values */
public void testMatricesValues1() {
   SingularValueDecomposition svd =
        new SingularValueDecompositionImpl(MatrixUtils.createRealMatrix(testSquare));
    RealMatrix uRef = MatrixUtils.createRealMatrix(new double[][] {
            { 3.0 / 5.0, -4.0 / 5.0 },
            { 4.0 / 5.0,  3.0 / 5.0 }
    });
    RealMatrix sRef = MatrixUtils.createRealMatrix(new double[][] {
            { 3.0, 0.0 },
            { 0.0, 1.0 }
    });
    RealMatrix vRef = MatrixUtils.createRealMatrix(new double[][] {
            { 4.0 / 5.0,  3.0 / 5.0 },
            { 3.0 / 5.0, -4.0 / 5.0 }
    });

    // check values against known references
    RealMatrix u = svd.getU();
    assertEquals(0, u.subtract(uRef).getNorm(), normTolerance);
    RealMatrix s = svd.getS();
    assertEquals(0, s.subtract(sRef).getNorm(), normTolerance);
    RealMatrix v = svd.getV();
    assertEquals(0, v.subtract(vRef).getNorm(), normTolerance);

    // check the same cached instance is returned the second time
    assertTrue(u == svd.getU());
    assertTrue(s == svd.getS());
    assertTrue(v == svd.getV());
    
}
 

/** test matrices values */
public void testMatricesValues1() {
   LUDecomposition lu =
        new LUDecompositionImpl(MatrixUtils.createRealMatrix(testData));
    RealMatrix lRef = MatrixUtils.createRealMatrix(new double[][] {
            { 1.0, 0.0, 0.0 },
            { 0.5, 1.0, 0.0 },
            { 0.5, 0.2, 1.0 }
    });
    RealMatrix uRef = MatrixUtils.createRealMatrix(new double[][] {
            { 2.0,  5.0, 3.0 },
            { 0.0, -2.5, 6.5 },
            { 0.0,  0.0, 0.2 }
    });
    RealMatrix pRef = MatrixUtils.createRealMatrix(new double[][] {
            { 0.0, 1.0, 0.0 },
            { 0.0, 0.0, 1.0 },
            { 1.0, 0.0, 0.0 }
    });
    int[] pivotRef = { 1, 2, 0 };

    // check values against known references
    RealMatrix l = lu.getL();
    assertEquals(0, l.subtract(lRef).getNorm(), 1.0e-13);
    RealMatrix u = lu.getU();
    assertEquals(0, u.subtract(uRef).getNorm(), 1.0e-13);
    RealMatrix p = lu.getP();
    assertEquals(0, p.subtract(pRef).getNorm(), 1.0e-13);
    int[] pivot = lu.getPivot();
    for (int i = 0; i < pivotRef.length; ++i) {
        assertEquals(pivotRef[i], pivot[i]);
    }

    // check the same cached instance is returned the second time
    assertTrue(l == lu.getL());
    assertTrue(u == lu.getU());
    assertTrue(p == lu.getP());
    
}
 

@Test
public void testMatricesValues() {
   BiDiagonalTransformer transformer =
        new BiDiagonalTransformer(MatrixUtils.createRealMatrix(testSquare));
   final double s17 = Math.sqrt(17.0);
    RealMatrix uRef = MatrixUtils.createRealMatrix(new double[][] {
            {  -8 / (5 * s17), 19 / (5 * s17) },
            { -19 / (5 * s17), -8 / (5 * s17) }
    });
    RealMatrix bRef = MatrixUtils.createRealMatrix(new double[][] {
            { -3 * s17 / 5, 32 * s17 / 85 },
            {      0.0,     -5 * s17 / 17 }
    });
    RealMatrix vRef = MatrixUtils.createRealMatrix(new double[][] {
            { 1.0,  0.0 },
            { 0.0, -1.0 }
    });

    // check values against known references
    RealMatrix u = transformer.getU();
    Assert.assertEquals(0, u.subtract(uRef).getNorm(), 1.0e-14);
    RealMatrix b = transformer.getB();
    Assert.assertEquals(0, b.subtract(bRef).getNorm(), 1.0e-14);
    RealMatrix v = transformer.getV();
    Assert.assertEquals(0, v.subtract(vRef).getNorm(), 1.0e-14);

    // check the same cached instance is returned the second time
    Assert.assertTrue(u == transformer.getU());
    Assert.assertTrue(b == transformer.getB());
    Assert.assertTrue(v == transformer.getV());

}