// solves triadiagonal linear systems
// using LU decomposition.
//
// Data:
//
// N - size of the system (NxN)
// DM - (i, i-1) diagonal
// D  - (i, i)   diagonal
// DP - (i, i+1) diagonal
//
// Author:      Alexander Bogomolny, CTK Software, Inc.
// URL:         http://www.cut-the-knot.com
// Date:        November 30, 2000
// Copyright:   A. Bogomolny
//              Permission to use and modify the file is therefore granted
//              as long as this comment remains unchanged. Do this at your
//              own risk.
//
class Tridiagonal
{
    int N;

    double[] DM, D, DP; // three diagonals

    // LU decomposition
    double[] A, C; // L has two non zero diagonals A & DM, U has units on main diagonal and C

    public Tridiagonal(int n, double[] m, double[] d, double[] p)
    {
        N = n;
        DM = m;
        D = d;
        DP = p;

        FindLU();
    }

    // the same object may be used repeatedly
    //
    public void SetMatrix(int n, double[] m, double[] d, double[] p)
    {
        N = n;
        DM = m;
        D = d;
        DP = p;

        FindLU();
    }

    public void SetMatrix(double[] m, double[] d, double[] p)
    {
        DM = m;
        D = d;
        DP = p;

        FindLU();
    }

    // standard LU decomposition.
    // See K. Atkinson "An Introduction
    // to Numerical Analysis", p 455 (Chapter 8)
    //
    public void FindLU()
    {
        A = new double[N];
        C = new double[N-1];

        A[0] = D[0];
        C[0] = DP[0]/D[0];

        for (int i = 1; i < N-1; i++)
        {
            A[i] = D[i] - DM[i-1]*C[i-1];
            C[i] = DP[i]/A[i];
        }

        A[N-1] = D[N-1] - DM[N-2]*C[N-2];
    }

    // solve the system for a given
    // right-hand side r
    //
    public double[] Solve(double[] r)
    {
        double[] z = new double[N];

        z[0] = r[0]/A[0];

        for (int i = 1; i < N; i++)
            z[i] = (r[i] - DM[i-1]*z[i-1])/A[i];

        double[] x = new double[N];

        x[N-1] = z[N-1];

        for (int i = N-2; i >= 0; i--)
            x[i] = z[i] - C[i]*x[i+1];

        return x;
    }
}

// solves triadiagonal linear systems
// using LU decomposition.
//
// Data:
//
// N - size of the system (NxN)
// DM - (i, i-1) diagonal
// D  - (i, i)   diagonal
// DP - (i, i+1) diagonal
//
// Author:      Alexander Bogomolny, CTK Software, Inc.
// URL:         http://www.cut-the-knot.com
// Date:        November 30, 2000
// Copyright:   A. Bogomolny
//              Permission to use and modify the file is therefore granted
//              as long as this comment remains unchanged. Do this at your
//              own risk.
//
class Tridiagonal
{
    int N;

    double[] DM, D, DP; // three diagonals

    // LU decomposition
    double[] A, C; // L has two non zero diagonals A & DM, U has units on main diagonal and C

    public Tridiagonal(int n, double[] m, double[] d, double[] p)
    {
        N = n;
        DM = m;
        D = d;
        DP = p;

        FindLU();
    }

    // the same object may be used repeatedly
    //
    public void SetMatrix(int n, double[] m, double[] d, double[] p)
    {
        N = n;
        DM = m;
        D = d;
        DP = p;

        FindLU();
    }

    public void SetMatrix(double[] m, double[] d, double[] p)
    {
        DM = m;
        D = d;
        DP = p;

        FindLU();
    }

    // standard LU decomposition.
    // See K. Atkinson "An Introduction
    // to Numerical Analysis", p 455 (Chapter 8)
    //
    public void FindLU()
    {
        A = new double[N];
        C = new double[N-1];

        A[0] = D[0];
        C[0] = DP[0]/D[0];

        for (int i = 1; i < N-1; i++)
        {
            A[i] = D[i] - DM[i-1]*C[i-1];
            C[i] = DP[i]/A[i];
        }

        A[N-1] = D[N-1] - DM[N-2]*C[N-2];
    }

    // solve the system for a given
    // right-hand side r
    //
    public double[] Solve(double[] r)
    {
        double[] z = new double[N];

        z[0] = r[0]/A[0];

        for (int i = 1; i < N; i++)
            z[i] = (r[i] - DM[i-1]*z[i-1])/A[i];

        double[] x = new double[N];

        x[N-1] = z[N-1];

        for (int i = N-2; i >= 0; i--)
            x[i] = z[i] - C[i]*x[i+1];

        return x;
    }
}