pcg_ilu0.c

二维热能方程的求解

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

void pcg_ilu0(nx,ny,itmax,idpcond,level,idpsrc,iterPCG,
              a,ailu,x,b,wksp,tol,ierr)
int nx,ny;
ATYPE a[];
ATYPE x[];
ATYPE ailu[];
ATYPE b[];
ATYPE wksp[][nx*ny];
VTYPE *tol;
int itmax,idpcond,level,idpsrc;
int *iterPCG,*ierr;
{
   int i,j,iter,nxny;
   int idp, idq, idr,idz;
   VTYPE rho0,rho1,alpha,beta,r0norm;
   VTYPE errnorm,product;

   if(level>=1) printf("PCG_ILU0:\n");
   if(level>=2) printf(" nx=%d ny=%d\n itmax=%d level=%d\n ierr=%d tol=%g\n",
             nx,ny,itmax,level,*ierr,*tol);

   nxny = nx*ny;

   idp=0;  /* wksp[idp] saves vector p  */
   idq=1;  /* wksp[idq] saves vector Ap */
   idr=2;  /* wksp[idr] saves vector r  */
   idz=3;  /* wksp[idz] saves vector z  */

 /************** initialization for iteration *******************/

   #if defined (_IBMR2) || defined (AIX) || defined (hpux)
      mtxvec(&nx,&ny,a,x,wksp[idq],&level,ierr);
   #else
      mtxvec_(&nx,&ny,a,x,wksp[idq],&level,ierr);
   #endif

   for(i=0; i<nxny; i++)
      wksp[idr][i]=b[i]-wksp[idq][i];   /* r=b-Ax */

   /* dot_prod(nxny,level,wksp[idr],wksp[idr],&product,ierr); */

   dot_prod(nxny,level,b,b,&product,ierr);
   r0norm=sqrt(product);

   for(i=0; i<nxny; i++)
      wksp[idz][i]=wksp[idr][i];        /* Mz=r */

   if(idpcond) ilu0_subst(nx,ny,level,ailu,wksp[idz],ierr);

   for(i=0; i<nxny; i++)
      wksp[idp][i]=wksp[idz][i];        /* p=z */

   dot_prod(nxny,level,wksp[idz],wksp[idr],&product,ierr); /* rho=z.r */
   rho0=product;

   if (*ierr!=0){printf("PCG_ILU: error in initialization.\n"); return;}

 /************** now, start the iteration ***********************/

   for(iter=0; iter<=itmax; iter++){

      #if defined (_IBMR2) || defined (AIX) || defined (hpux)
         mtxvec(&nx,&ny,a,wksp[idp],wksp[idq],&level,ierr);
      #else
         mtxvec_(&nx,&ny,a,wksp[idp],wksp[idq],&level,ierr);
      #endif

      dot_prod(nxny,level,wksp[idp],wksp[idq],&product,ierr);
      alpha=rho0/product;

      for(i=0; i<nxny; i++)
         x[i] += alpha*wksp[idp][i];

      for(i=0; i<nxny; i++)
         wksp[idr][i] -= alpha*wksp[idq][i];

      dot_prod(nxny,level,wksp[idr],wksp[idr],&product,ierr);
      errnorm=sqrt(product)/r0norm;
      if(level>=3) printf(" errnorm[%d]=%g\n",iter,errnorm);

      if(errnorm<*tol)
         {if(level>0) printf(" CONVERGENCE in %d iterations.\n",iter);
          *iterPCG=iter; break; }

      for(i=0; i<nxny; i++)
         wksp[idz][i]=wksp[idr][i];        /* Mz=r */

      if(idpcond) ilu0_subst(nx,ny,level,ailu,wksp[idz],ierr);

      dot_prod(nxny,level,wksp[idz],wksp[idr],&product,ierr);
      rho1=product;

      beta=rho1/rho0;

      for(i=0; i<nxny; i++)
         wksp[idp][i] = wksp[idz][i]+ beta*wksp[idp][i];

      if(level>=5){
         printf(" rho  [%d]=%f\n",iter,rho0);
         printf(" alpha[%d]=%f\n",iter,alpha);
         printf(" beta [%d]=%f\n",iter,beta);
         }

      rho0=rho1;

      if (*ierr!=0){printf("PCG_ILU: error in iteration %d\n.",iter); 
          return;}

   }  /* end of the iteration loop */

   if(iter>=itmax)
      printf(" PCG-ILU did NOT converge in %d iterations.\n",iter-1);

}