ilu0.c

来自「二维热能方程的求解」· C语言 代码 · 共 95 行

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

void ilu0(nx,ny,level,a,ailu,ierr)

int nx,ny,level;
int *ierr;
ATYPE a[][nx][5];
ATYPE ailu[][nx][5];

{
   int   i,j,nxm1,nym1,jp1;
   int   notdone,itcount,maxiter;
   VTYPE pivot,mul,pivot_min;
   VTYPE diag_add,diag_more;

   if(level>=1) printf("ILU0: Incomplete LU Factorization\n");
   if(level>=2){
      printf(" :Stabilized Outer-Product Form\n");
      printf(" nx=%d ny=%d\n level=%d ierr=%d\n",nx,ny,level,*ierr);
     }

   if(nx>2 && ny>2)
      pivot_min=a[2][2][2]*1.e-3;
   else
      pivot_min=a[1][1][2]*1.e-3;
   if(pivot_min<0.0){
      printf(" WARNING from ILU0: diagonal shoul be positive\n");
     }
   itcount=0;
   maxiter=10;
   diag_add=0.0;
   diag_more=0.1;

   nxm1=nx-1;
   nym1=ny-1;

   do{ notdone=0;
       for(j=0; j<ny; j++)
       for(i=0; i<nx; i++){
          ailu[j][i][0]=a[j][i][0];
          ailu[j][i][1]=a[j][i][1];
          ailu[j][i][2]=a[j][i][2]+diag_add;
          ailu[j][i][3]=a[j][i][3];
          ailu[j][i][4]=a[j][i][4];
          }

       for(j=0; j<ny; j++){
          jp1=j+1;
          for(i=0; i<nx; i++){
             pivot=ailu[j][i][2];
             if(pivot<pivot_min){
                printf("ILU0: pivot < pivot_min: Let us try again.\n");
                notdone=1; break;
                }
             if(i<nxm1){
                mul=ailu[j][i+1][1]/pivot;
                ailu[j][i+1][1] =mul;
                ailu[j][i+1][2]-=mul*ailu[j][i][3];
                }

             if(j<nym1){
                mul=ailu[jp1][i][0]/pivot;
                ailu[jp1][i][0] =mul;
                ailu[jp1][i][2]-=mul*ailu[j][i][4];
                }
             }
             if(notdone) break;
          }  /* end of the nested loop */

       if( (++itcount>=maxiter) && (notdone) ) {
           printf("ERROR from ILU0:\n");
           printf(" after %d iterations, still pivot < pivot_min.\n",itcount);
           *ierr=44; return;
           }

       diag_add+=diag_more;

    } while(notdone);


   /*************** print out ***************/

   if(level>=3){
      for(j=0; j<ny; j++)
         for(i=0; i<nx; i++)
             printf(" ailu[%d][%d][]=%f  %f  %f  %f  %f\n", j,i,
                ailu[j][i][0],ailu[j][i][1],ailu[j][i][2],
                ailu[j][i][3],ailu[j][i][4]);
      }
}