#define xm(p,q) *(x+q+p*n)
#define bm(p,q) *(b+q+p*n)
void multidiag(int n,int L,float *d,float *e,
    float *f,float *m,float *u,float *x,float *b)
{
/*
  solves tridiagonal system according to
  G.E. Forsythe and C.B. Moler " Computer
  solutions to linear algebraic systems,"
  Prentice-Hall publ., 1967, with multiple
  right hand sides: one T matrix, L solutions,

      TX = B

  where T = (d0, f0, 0,          ....)
            (e1, d1, f1, 0       ....)
            ( 0, e2, d2, f2, 0,  ....)
            (        ....            )
            ( 0, ...    0, en-1, dn-1)
  and
        X = X[0..n-1][0...L-1], L x-vectors
        B = B[0..n-1][0...L-1], L b-vectors

                              wpp 5/11/2002
*/
   int i,j;
   u[0] = e[0];
   for(j=0;j<L;j++) xm(0,j) = bm(0,j);
   for(i=1;i<n;i++){
       m[i] = d[i]/u[i-1];
       u[i] = e[i] - m[i]*f[i-1];
       for(j=0;j<L;j++){
           xm(i,j) = bm(i,j) - m[i]*xm(i-1,j);
       }
   }
   for(j=1;j<L;j++) xm(n-1,j) = xm(n-1,j)/u[n-1];
   for(i=n-2;i>=0;i--){
       for(j=0;j<L;j++){
           xm(i,j) = (xm(i,j)-f[i]*xm(i+1,j))/u[i];
       }
   }
}
#undef xm
#undef bm