TridiagonalSolver { public static double[] solvTriDag(final TridiagonalMatrix aM, final double[] b) { ArgumentChecker.notNull(aM, "null matrix"); ArgumentChecker.notNull(b, "null vector"); final double[] d = aM.getDiagonal(); final int n = d.length; ArgumentChecker.isTrue(n == b.length, "vector y wrong length for matrix"); final double[] y = Arrays.copyOf(b, n); final double[] l = aM.getLowerSubDiagonalData(); final double[] u = aM.getUpperSubDiagonalData(); final double[] x = new double[n]; for (int i = 1; i < n; i++) { final 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; }  static double[] solvTriDag(final TridiagonalMatrix aM, final double[] b); static DoubleMatrix1D solvTriDag(final TridiagonalMatrix aM, final DoubleMatrix1D b);  }
@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 DoubleMatrix1D xVec = new DoubleMatrix1D(x); final DoubleMatrix1D yVec = (DoubleMatrix1D) MA.multiply(m, xVec); final double[] xSolv = solvTriDag(m, yVec).getData(); for (int i = 0; i < n; i++) { assertEquals(x[i], xSolv[i], 1e-9); } DoubleMatrix1D resi = (DoubleMatrix1D) MA.subtract(MA.multiply(m, new DoubleMatrix1D(xSolv)), yVec); double err = MA.getNorm2(resi); assertEquals(0.0, err, 1e-14); }