/* 
  Thomas algorithm solving tridiagonal linear system. 
	Description:
   int Sweep(double *u,double *a,double *b,double *c,double *f,
             int n,double *bet,double *gam);
	Parameters:
   u - solution (size n), on output;
   a,b,c,f - arrays coefficients (as in the algorithm's description above, size n)
   n - the solution and coefficients arrays size;
   bet, gam - additional arrays (size n+1).
	   Returns: 
   0 - if the algorithm fails,
   1 - if OK.

  Reduction algorithm to solve tridiagonal system.
	Description:
   void Reduce(double *u,double *a,double *b,double *c,double *f,int n);
	Parameters:
   u - solution (size n), on output;
   a,b,c,f - arrays coefficients (as in the algorithm's description above, size n)
   n - the solution and coefficients arrays size;
   Note:
   The algorithm is implenmented without singularity checkup for better 
   performance. To apply it, one must include checkup for division by 
   zero any time the division takes place.

*/

/* Thomas algorithm */
int Sweep(double *u,double *a,double *b,double *c,double *f,int n,
          double *bet,double *gam) {
 register int n;
 double fi;
 /* calculating the additional arrays (factorization of the system) */
 bet[0]=gam[0]=0.;
 for(i=0;i<n;i++) {
  fi=b[i];
  if(i>0) fi-=a[i]*bet[i];
  if(fi==0.) return(0);
  fi=1./fi; 
  if(i<n-1) bet[i+1]=c[i]*fi; 
  gam[i+1]=(f[i]+a[i]*gam[i])*fi;
 }
 /* back substitution */
 u[n-1]=gam[n];
 for(i=n-1;i&gt;0;i--) u[i-1]=bet[i]*u[i]+gam[i];
 return(1);
}

/* Reduction techniques */
void Reduce(double *u,double *a,double *b,double *c,double *f,int n) {
 register int i,j;int quarter,half,full;
 double a1,c1;
 /* Calculating quarter, central and "right" point numbers */
 i=n;full=1;
 while(i>>=1) full<<=1;
 if(n>full+1) full<<=1;	
 half=full>>1;quarter=half>>1;
 /* Forward pass: coefficients recalculation 
    Note: coefficients at odd positions remain intact! */
 for(i=1;i<half;i<<=1) {
  /* leftmost point case */
  c1=c[0]/b[i]; a[0]=0.; c[0]=c[i]*c1; b[0]-=a[i]*c1; f[0]+=f[i]*c1;
  for(j=i+i;j<n;j+=i+i) {
   if(j+i>=n) {		/* too far points (including the rightmost) */
    a1=a[j]/b[j-i]; c[j]=0.; a[j]=a[j-i]*a1; b[j]-=c[j-i]*a1; f[j]+=f[j-i]*a1;
   }
   else {		/* the main case */
    c1=c[j]/b[j+i]; a1=a[j]/b[j-i];
    c[j]=c[j+i]*c1; a[j]=a[j-i]*a1;
    b[j]-=a[j+i]*c1+c[j-i]*a1; f[j]+=f[j+i]*c1+f[j-i]*a1;
   }
  }
 }
 /* Central coefficients b and f calculation */
 if(full>=n) {		/* case when n!=2^order+1 */
  a1=a[half]/b[0]; b[half]-=c[0]*a1; f[half]+=f[0]*a1;
 }
 else {			/* when it's equal */
  c1=c[half]/b[full]; a1=a[half]/b[0];
  b[half]-=a[full]*c1+c[0]*a1; f[half]+=f[full]*c1+f[0]*a1;
 }
 /* results in the central and leftmost points */
 u[half]=f[half]/b[half];u[0]=(c[0]*u[half]+f[0])/b[0];
 /* result in the rightmost point if n==2^order+1 */
 if(full<n) u[full]=(a[full]*u[half]+f[full])/b[full];
 /* Backward pass: results in all the other points */
 for(i=quarter;i>=1;i>>=1) for(j=i;j<n;j+=i+i) {
  /* too far points case */
  if(j+i>=n) u[j]=(f[j]+a[j]*u[j-i])/b[j];
  /* the main case */
  else u[j]=(f[j]+c[j]*u[j+i]+a[j]*u[j-i])/b[j];
 }
}