@Override protected double[] matrixEqnSolver(final double[][] doubMat, final double[] doubVec) { final int sizeM1 = doubMat.length - 1; final double[] a = new double[sizeM1]; final double[] b = new double[sizeM1 + 1]; final double[] c = new double[sizeM1]; for (int i = 0; i < sizeM1; ++i) { a[i] = doubMat[i][i + 1]; b[i] = doubMat[i][i]; c[i] = doubMat[i + 1][i]; } b[sizeM1] = doubMat[sizeM1][sizeM1]; final TridiagonalMatrix m = new TridiagonalMatrix(b, a, c); return TridiagonalSolver.solvTriDag(m, doubVec); }
private DoubleArray multiply(TridiagonalMatrix matrix, DoubleArray vector) { double[] a = matrix.getLowerSubDiagonalData(); double[] b = matrix.getDiagonalData(); double[] c = matrix.getUpperSubDiagonalData(); double[] x = vector.toArrayUnsafe(); int n = x.length; ArgChecker.isTrue(b.length == n, "Matrix/vector size mismatch"); double[] res = new double[n]; int i; res[0] = b[0] * x[0] + c[0] * x[1]; res[n - 1] = b[n - 1] * x[n - 1] + a[n - 2] * x[n - 2]; for (i = 1; i < n - 1; i++) { res[i] = a[i - 1] * x[i - 1] + b[i] * x[i] + c[i] * x[i + 1]; } return DoubleArray.ofUnsafe(res); }
/** * @return Returns the tridiagonal matrix as a {@link DoubleMatrix} */ public DoubleMatrix toDoubleMatrix() { if (_matrix == null) { calMatrix(); } return _matrix; }
/** * Solves the system Ax = y for the unknown vector x, where A is a tridiagonal matrix and y is a vector. * This takes order n operations where n is the size of the system * (number of linear equations), as opposed to order n^3 for the general problem. * @param aM tridiagonal matrix * @param b known vector (must be same length as rows/columns of matrix) * @return vector (as an array of doubles) with same length as y */ public static double[] solvTriDag(TridiagonalMatrix aM, double[] b) { ArgChecker.notNull(aM, "null matrix"); ArgChecker.notNull(b, "null vector"); double[] d = aM.getDiagonal(); //b is modified, so get copy of diagonal int n = d.length; ArgChecker.isTrue(n == b.length, "vector y wrong length for matrix"); double[] y = Arrays.copyOf(b, n); double[] l = aM.getLowerSubDiagonalData(); double[] u = aM.getUpperSubDiagonalData(); double[] x = new double[n]; for (int i = 1; i < n; i++) { double m = l[i - 1] / d[i - 1]; d[i] = d[i] - m * u[i - 1]; y[i] = y[i] - m * y[i - 1]; } x[n - 1] = y[n - 1] / d[n - 1]; for (int i = n - 2; i >= 0; i--) { x[i] = (y[i] - u[i] * x[i + 1]) / d[i]; } return x; }
@Test public void testGetters() { assertTrue(Arrays.equals(A, M.getDiagonalData())); assertTrue(Arrays.equals(B, M.getUpperSubDiagonalData())); assertTrue(Arrays.equals(C, M.getLowerSubDiagonalData())); final int n = A.length; final DoubleMatrix matrix = M.toDoubleMatrix(); for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { if (i == j) { assertEquals(matrix.get(i, j), A[i], 0); } else if (j == i + 1) { assertEquals(matrix.get(i, j), B[j - 1], 0); } else if (j == i - 1) { assertEquals(matrix.get(i, j), C[j], 0); } else { assertEquals(matrix.get(i, j), 0, 0); } } } }
@Test public void testHashCodeAndEquals() { final double[] a = Arrays.copyOf(A, A.length); final double[] b = Arrays.copyOf(B, B.length); final double[] c = Arrays.copyOf(C, C.length); TridiagonalMatrix other = new TridiagonalMatrix(a, b, c); assertEquals(other, M); assertEquals(other.hashCode(), M.hashCode()); a[1] = 1000; other = new TridiagonalMatrix(a, B, C); assertFalse(other.equals(M)); b[1] = 1000; other = new TridiagonalMatrix(A, b, C); assertFalse(other.equals(M)); c[1] = 1000; other = new TridiagonalMatrix(A, B, c); assertFalse(other.equals(M)); } }
@Test public void testTridiagonalMultiply() { final int n = 37; double[] l = new double[n - 1]; double[] c = new double[n]; double[] u = new double[n - 1]; double[] x = new double[n]; for (int ii = 0; ii < n; ii++) { c[ii] = RANDOM.nextRandom(); x[ii] = RANDOM.nextRandom(); if (ii < n - 1) { l[ii] = RANDOM.nextRandom(); u[ii] = RANDOM.nextRandom(); } } final TridiagonalMatrix m = new TridiagonalMatrix(c, u, l); final DoubleArray xVec = DoubleArray.copyOf(x); DoubleArray y1 = (DoubleArray) ALGEBRA.multiply(m, xVec); DoubleMatrix full = m.toDoubleMatrix(); DoubleArray y2 = (DoubleArray) ALGEBRA.multiply(full, xVec); for (int i = 0; i < n; i++) { assertEquals(y1.get(i), y2.get(i), 1e-12); } }
private DoubleArray multiply(DoubleArray vector, TridiagonalMatrix matrix) { double[] a = matrix.getLowerSubDiagonalData(); double[] b = matrix.getDiagonalData(); double[] c = matrix.getUpperSubDiagonalData(); double[] x = vector.toArrayUnsafe(); int n = x.length; ArgChecker.isTrue(b.length == n, "Matrix/vector size mismatch"); double[] res = new double[n]; int i; res[0] = b[0] * x[0] + a[0] * x[1]; res[n - 1] = b[n - 1] * x[n - 1] + c[n - 2] * x[n - 2]; for (i = 1; i < n - 1; i++) { res[i] = a[i] * x[i + 1] + b[i] * x[i] + c[i - 1] * x[i - 1]; } return DoubleArray.ofUnsafe(res); }
@Test(expectedExceptions = IllegalArgumentException.class) public void testNullB() { new TridiagonalMatrix(A, null, C); }
@Override public DoubleMatrix apply(TridiagonalMatrix x) { ArgChecker.notNull(x, "x"); double[] a = x.getDiagonalData(); double[] b = x.getUpperSubDiagonalData(); double[] c = x.getLowerSubDiagonalData(); int n = a.length; int i, j, k;
@Test(expectedExceptions = IllegalArgumentException.class) public void testNullC() { new TridiagonalMatrix(A, B, null); }
@Test(expectedExceptions = IllegalArgumentException.class) public void testWrongC() { new TridiagonalMatrix(A, B, new double[] {1, 2, 3, 4, 5, 6, 7 }); }
@Test(expectedExceptions = IllegalArgumentException.class) public void testNullA() { new TridiagonalMatrix(null, B, C); }
@Test(expectedExceptions = IllegalArgumentException.class) public void testWrongB() { new TridiagonalMatrix(A, new double[] {1, 2, 3, 4, 5, 6, 7, 8 }, C); }
private static DoubleMatrix getInverseTridiagonalMatrix(double[] deltaX, boolean leftNatural, boolean rightNatural) { InverseTridiagonalMatrixCalculator invertor = new InverseTridiagonalMatrixCalculator(); int n = deltaX.length + 1; double[] a = new double[n]; double[] b = new double[n - 1]; double[] c = new double[n - 1]; for (int i = 1; i < n - 1; i++) { a[i] = (deltaX[i - 1] + deltaX[i]) / 3.0; b[i] = deltaX[i] / 6.0; c[i - 1] = deltaX[i - 1] / 6.0; } // Boundary condition if (leftNatural) { a[0] = 1.0; b[0] = 0.0; } else { a[0] = -deltaX[0] / 3.0; b[0] = deltaX[0] / 6.0; } if (rightNatural) { a[n - 1] = 1.0; c[n - 2] = 0.0; } else { a[n - 1] = deltaX[n - 2] / 3.0; c[n - 2] = deltaX[n - 2] / 6.0; } TridiagonalMatrix tridiagonal = new TridiagonalMatrix(a, b, c); return invertor.apply(tridiagonal); }
@Test public void testInvertIdentity() { final int n = 11; final double[] a = new double[n]; final double[] b = new double[n - 1]; final double[] c = new double[n - 1]; int i, j; for (i = 0; i < n; i++) { a[i] = 1.0; } final DoubleMatrix res = CALCULATOR.apply(new TridiagonalMatrix(a, b, c)); for (i = 0; i < n; i++) { for (j = 0; j < n; j++) { assertEquals((i == j ? 1.0 : 0.0), res.get(i, j), EPS); } } }
@Override protected DoubleArray[] combinedMatrixEqnSolver(final double[][] doubMat1, final double[] doubVec, final double[][] doubMat2) { final int size = doubVec.length; final DoubleArray[] res = new DoubleArray[size + 1]; final DoubleMatrix doubMat2Matrix = DoubleMatrix.copyOf(doubMat2); final double[] u = new double[size - 1]; final double[] d = new double[size]; final double[] l = new double[size - 1]; for (int i = 0; i < size - 1; ++i) { u[i] = doubMat1[i][i + 1]; d[i] = doubMat1[i][i]; l[i] = doubMat1[i + 1][i]; } d[size - 1] = doubMat1[size - 1][size - 1]; final TridiagonalMatrix m = new TridiagonalMatrix(d, u, l); res[0] = DoubleArray.copyOf(TridiagonalSolver.solvTriDag(m, doubVec)); for (int i = 0; i < size; ++i) { DoubleArray doubMat2Colum = doubMat2Matrix.column(i); res[i + 1] = TridiagonalSolver.solvTriDag(m, doubMat2Colum); } return res; }
@Test public void test() { final int n = 97; double[] a = new double[n - 1]; double[] b = new double[n]; double[] c = new double[n - 1]; double[] x = new double[n]; for (int ii = 0; ii < n; ii++) { b[ii] = RANDOM.nextRandom(); x[ii] = RANDOM.nextRandom(); if (ii < n - 1) { a[ii] = RANDOM.nextRandom(); c[ii] = RANDOM.nextRandom(); } } final TridiagonalMatrix m = new TridiagonalMatrix(b, a, c); final DoubleArray xVec = DoubleArray.copyOf(x); final DoubleArray yVec = (DoubleArray) MA.multiply(m, xVec); final double[] xSolv = solvTriDag(m, yVec).toArray(); for (int i = 0; i < n; i++) { assertEquals(x[i], xSolv[i], 1e-9); } DoubleArray resi = (DoubleArray) MA.subtract(MA.multiply(m, DoubleArray.copyOf(xSolv)), yVec); double err = MA.getNorm2(resi); assertEquals(0.0, err, 1e-14); }