void tridiag

/**********************************************************************\
 * Subroutine to solve a tridiagonal system of the form                *
 *   diags[0][j]*u[j-1] + diags[1][j]*u[j] + diags[2][j]*u[j+1] = b[j] *
 * for j = 0, 1, 2, 3, ... J (where diags[0][0] = diags[2][J] = 0.)    *
 *                                                                     *
 * The vector b is altered by tridiag and contains the values of the   *
 * solution vector u after tridiag.                                    *
 *                                                                     *
 * The matrix diags is also altered:                                   *
 * therefore, to solve the same system, but with different right hand  *
 * side b, again, restore the matrix diags.                            *
 *                                                                     *
 * No pivoting is performed.                                           *
 \*********************************************************************/

(
 int J,                    // index of the last unknown
 float diags[3][J_DIM],    // the coefficients of the equations (destroyed)
 float b[J_DIM]            // in: the right hand sides   out: the solution x
 )

{
  /************************ local variables *************************/

  // equation index
  int j;

  // amount of the previous equation to substract
  float fac;

  /********************* executable statements **********************/

  // reduce the matrix to bidiagonal form:
  for(j=1; j<=J; j++)
    {
      // amount of the previous equation to substract to eliminate diags[0][j]
      fac=diags[0][j]/diags[1][j-1];

      // substract this amount
      diags[1][j]-=fac*diags[2][j-1];
      b[j]-=fac*b[j-1];

    }

  // solve the bidiagonal system and store the solution in b[j]:
  b[J]=b[J]/diags[1][J];
  for(j=J-1; j>=0; j--)
    {
      b[j]=(b[j]-diags[2][j]*b[j+1])/diags[1][j];
    }

}