////////////////////////////////////////////////////////////////////
// Solve a tridiagonal matrix equation Ax=b.
//
// Be sure to have zero as the first element of "upper"
// and the last element of "lower".
//
// Entry:
//   lower - matrix entries just under diagonal
//   diag  - diagonal matrix entries
//   upper - entries above the diagonal
//   b     - vector of known values (right-hand side)
//   num   - size of input vector (and diagonals)
//      
// Exit:
//   x - solution vector
////////////////////////////////////////////////////////////////////

void tridiagonal_solve(
  float[] lower,
  float[] diag,
  float[] upper,
  float[] b,
  float[] x,
  int num
)
{
  int j;
  float beta;
  float[] gamma = new float[num];

  if (diag[0] == 0.0) {
    println ("Cannot have zero in first element of diagonal.");
    exit();
  }

  beta = diag[0];
  x[0] = b[0] / beta;

  for (j = 1; j < num; j++) {
    gamma[j] = upper[j-1] / beta;
    beta = diag[j] - lower[j] * gamma[j];
    if (beta == 0.0) {
      println ("Error in tridiagonal solve.");
      exit ();
    }
    x[j] = (b[j] - lower[j] * x[j-1]) / beta;
  }

  for (j = num-2; j >= 0; j--)
    x[j] -= gamma[j+1] * x[j+1];
}