ilu0_subst.c

二维热能方程的求解

#include <stdio.h>
#include <math.h>
typedef  double ATYPE;
typedef  float DTYPE;

void ilu0_subst(nx,ny,level,ailu,x,ierr)
int nx,ny,level;
int *ierr;
ATYPE ailu[][nx][5];
ATYPE x[][nx];
{
   int i,j;

  /********** forward ILU substitution ***************/

    j=0;
       for(i=1; i<nx; i++)
          x[j][i]= (x[j][i]-ailu[j][i][1]*x[j][i-1]);

    for(j=1; j<ny; j++){
       i=0;
          x[j][i]= (x[j][i]-ailu[j][i][0]*x[j-1][i]);
       for(i=1; i<nx; i++)
          x[j][i]= (x[j][i]-ailu[j][i][0]*x[j-1][i]
                           -ailu[j][i][1]*x[j][i-1]);
       }

  /********** back ILU substitution ***************/

    j=ny-1;
       i=nx-1;
          x[j][i] /= ailu[j][i][2];
       for(i=nx-2; i>=0; i--)
          x[j][i]= (x[j][i]-ailu[j][i][3]*x[j][i+1])/ailu[j][i][2];

    for(j=ny-2; j>=0; j--){
       i=nx-1;
          x[j][i]= (x[j][i]-ailu[j][i][4]*x[j+1][i])/ailu[j][i][2];
       for(i=nx-2; i>=0; i--)
          x[j][i]= (x[j][i]-ailu[j][i][4]*x[j+1][i]
                           -ailu[j][i][3]*x[j][i+1])/ailu[j][i][2];
       }
}