mainexpert.c

数学算法的实现库。可以实现常见的线性计算。

#endif
	  }// end for i

	  // uncompress rhs
	  for (i=0; i<rhsptr[1]-rhsptr[0]; i++) {
	      j=rhsind[i]-1;
	      rhs[j]=rhsval[i];
	  } // end for i
       } // end if
       for (i=0; i<n; i++) {
#if defined _DOUBLE_REAL_ || defined _SINGLE_REAL_
	   rhs[i]*=leftscale[i];
#else
	   val=rhs[i].r;
	   rhs[i].r=val*leftscale[i].r-rhs[i].i*leftscale[i].i;
	   rhs[i].i=val*leftscale[i].i+rhs[i].i*leftscale[i].r;
#endif
       } // end for i
    } // end else


    //printf("matrix-vector multiply done\n");fflush(STDOUT);
    /*
    printf("true right hand side\n");
    for (i=0; i<n; i++)
      printf("%12.4e",rhs[i]);
    printf("\n");
    fflush(STDOUT);
    printf("start GMRES\n");fflush(STDOUT);
    */

    /* --------------------   apply preconditioner   --------------------- */
    if (ierr==0)
    {
#ifdef USE_SPARSKIT
       /* initial solution */
       if (rhstyp[1]=='G' || rhstyp[1]=='g') {
	  j=nrhs*n;
	  if (rhstyp[0]=='M' || rhstyp[0]=='m')
	     j=n;
	  for (i=0; i<n; i++,j++)
	      sol[i]=rhs[j];
       }
       else
	  for (i=0; i<n; i++)
#if defined _DOUBLE_REAL_ || defined _SINGLE_REAL_
	      sol[i]=0;
#else
	      sol[i].r=sol[i].i=0;
#endif


       /* init solver */
       param.ipar[20]=0;
       /* right preconditioning */
       param.ipar[21]=2;
       if (param.ipar[21]==0)
	 printf("no preconditioning\n");
       else if (param.ipar[21]==1)
	 printf("left preconditioning\n");
       else if (param.ipar[21]==2)
	 printf("right preconditioning\n");
       else if (param.ipar[21]==3)
	 printf("both sides preconditioning\n");
       
       /* stopping criteria, backward error */
       param.ipar[22]=3;
       /* number of restart length for GMRES,FOM,... */
       param.ipar[24]=nrestart;
       /* maximum number of iteration steps */
       param.ipar[25]=max_it;
       /* init with zeros */
       param.ipar[26]=
	 param.ipar[27]=
	 param.ipar[28]=
	 param.ipar[29]=
	 param.ipar[30]=
	 param.ipar[31]=
	 param.ipar[32]=
	 param.ipar[33]=0;
       /* relative error tolerance */
       param.fpar[20]=restol;
       /* absolute error tolerance */
       param.fpar[21]=0;
       /* init with zeros */
       param.fpar[22]=
	 param.fpar[23]=
	 param.fpar[24]=
	 param.fpar[25]=
	 param.fpar[26]=
	 param.fpar[27]=
	 param.fpar[28]=
	 param.fpar[29]=
	 param.fpar[30]=0;
       
       /* allocate memory for iterative solver depending on our choice */
       i=param.ipar[24];
       switch (SOLVER)
	 {
	 case  1:   /* pcg */
	   i=5*n;
	   printf("use PCG as iterative solver\n");
	   break;
	 case  2:   /* cgnr */
	   i=5*n;
	   printf("use CGNR as iterative solver\n");
	   break;
	 case  3:   /* bcg */
	   printf("use BCG as iterative solver\n");
	   i=7*n;
	   break;
	 case  4:   /* dbcg */
	   i=11*n;
	   printf("use DBCG as iterative solver\n");
	   break;
	 case  5:   /* bcgstab */
	   i=8*n;
	   printf("use BCGSTAB as iterative solver\n");
	   break;
	 case  6:   /* tfqmr */
	   i=11*n;
	   printf("use TFQMR as iterative solver\n");
	   break;
	 case  7:   /* fom */
	   i=(n+3)*(i+2)+(i+1)*i/2;
	   printf("use FOM(%d) as iterative solver\n",param.ipar[24]);
	   break;
	 case  8:   /* gmres */
	   i=(n+3)*(i+2)+(i+1)*i/2;
	   printf("use GMRES(%d) as iterative solver\n",param.ipar[24]);
	   break;
	 case  9:   /* fgmres */
	   i=2*n*(i+1) + ((i+1)*i)/2 + 3*i + 2;
	   printf("use FGMRES(%d) as iterative solver\n",param.ipar[24]);
	   break;
	 case 10:   /* dqgmres */
	   i=n+(i+1)*(2*n+4);
	   printf("use DQGMRES(%d) as iterative solver\n",param.ipar[24]);
	   break;
	 } /* end switch */
       // printf("allocate %d spaces of memory\n",i);
       w=(FLOAT *)MALLOC((size_t)i*sizeof(FLOAT),"main:w"); 
       param.ipar[23]=i;
       


       // init buffer for column sumns
       rbuff=(REALS *)dbuff;
       for (i=0; i<n; i++) 
	 rbuff[i]=0.0;
       // val = maximum row sum
       val=0.0;
       for (i=0; i<n; i++) {
	 j=1;
	 k=A.ia[i+1]-A.ia[i];
	 val=MAX(val,ASUM(&k,A.a+A.ia[i]-1,&j));
	 for (j=A.ia[i]-1; j<A.ia[i+1]-1; j++) {
	   k=A.ja[j]-1;
	   rbuff[k]+=FABS(A.a[j]);
	 } // end for j
       } // end for i
       vb=0;
       for (i=0; i<n; i++)
	 vb=MAX(vb,rbuff[i]);
       param.fpar[11]=sqrt((double)val*vb);



       /* start iterative solver */
       evaluate_time(&time_start,&systime);
       flags=-1;
       while (flags)
	 {
	   /* which solver do we use? */
	   switch (SOLVER)
	     {
	     case  1:   /* pcg */
	       PCG(&n,rhs,sol,param.ipar+20,param.fpar+20,w);
	       break;
	     case  3:   /* bcg */
	       BCG(&n,rhs,sol,param.ipar+20,param.fpar+20,w);
	       break;
	     case  8:   /* gmres */
	       GMRES(&n,rhs,sol,param.ipar+20,param.fpar+20,w);
	       break;
	     case  9:   /* fgmres */
	       FGMRES(&n,rhs,sol,param.ipar+20,param.fpar+20,w);
	       break;
	     default:
	       printf("solver not available\n");
	       exit(0);
	       break;
	     } /* end switch */
	   /* what is needed for the solver? */
	   


	   w--;
	   switch (param.ipar[20])
	     {
	     case  1:   /* apply matrix vector multiplication */
	       /*
		 printf("right hand side\n");
		 for (i=0; i<n;i++)
		     printf("%12.4e",w[param.ipar[27]+i]);
		 printf("\n");
		 fflush(STDOUT);
	       */
	       // printf("apply mat vec\n"); fflush(STDOUT);
	       	 
	       evaluate_time(&timeAx_start,&systime);
	       nAx++;
	       CSRMATVEC(A,w+param.ipar[27],w+param.ipar[28]);
	       evaluate_time(&timeAx_stop,&systime);
	       secndsAx+=timeAx_stop-timeAx_start;
	       
	       //  printf("mat vec applied\n"); fflush(STDOUT);
	       /* 
		 printf("mat vec solution\n");
		 for (i=0; i<n;i++)
		     printf("%12.4e",w[param.ipar[28]+i]);
		 printf("\n");
		 fflush(STDOUT);
	       */
	       break;
	     case  2:   /* apply transposed matrix vector multiplication */
	       //   printf("apply matt vec\n"); fflush(STDOUT);
	       evaluate_time(&timeAx_start,&systime);
	       nAx++;
	       CSRMATTVEC(A,w+param.ipar[27],w+param.ipar[28]);
	       evaluate_time(&timeAx_stop,&systime);
	       secndsAx+=timeAx_stop-timeAx_start;
	       //  printf("matt vec applied\n"); fflush(STDOUT);
	       break;
	     case  3:   /* (no) left preconditioner solver */
	       //   printf("apply prec vec\n"); fflush(STDOUT);
	       i=1;
	       if (param.ipar[21]==1)
		  AMGSOL(PRE, nlev, &param, w+param.ipar[27],w+param.ipar[28], param.dbuff);
	       else
		 if (param.ipar[21]==3)
		    AMGLSOL(PRE, nlev, &param, w+param.ipar[27],w+param.ipar[28], param.dbuff);
		 else
		    COPY(&n,w+param.ipar[27],&i,w+param.ipar[28],&i);
	       //   printf("prec vec applied\n"); fflush(STDOUT);
	       break;
	     case  4:   /* (no) left preconditioner transposed solver */
	       //   printf("apply prect vec\n"); fflush(STDOUT);
	       i=1;
	       if (param.ipar[21]==1)
		  AMGTSOL(PRE, nlev, &param, w+param.ipar[27],w+param.ipar[28], param.dbuff);
	       else
		  if (param.ipar[21]==3)
		     AMGTLSOL(PRE, nlev, &param, w+param.ipar[27],w+param.ipar[28], param.dbuff);
		  else
		     COPY(&n,w+param.ipar[27],&i,w+param.ipar[28],&i);
	       //   printf("prect vec applied\n"); fflush(STDOUT);
	       break;
	     case  5:   /* (no) right preconditioner solver */
	       
	       // printf("apply prec vec\n"); fflush(STDOUT);
	       /*
		 printf("original right hand side\n");
		 for (i=0; i<n;i++)
		     printf("%12.4e",w[param.ipar[27]+i]);
		  printf("\n");
		  fflush(STDOUT);
	       */
	       i=1;
	       if (param.ipar[21]==2)
		 AMGSOL(PRE, nlev, &param, w+param.ipar[27],w+param.ipar[28], param.dbuff);
	       else 
		  if (param.ipar[21]==3)
		     AMGDUSOL(PRE, nlev, &param, w+param.ipar[27],w+param.ipar[28], param.dbuff);
		  else
		     COPY(&n,w+param.ipar[27],&i,w+param.ipar[28],&i);
	       // printf("prec vec applied\n"); fflush(STDOUT);
	       break;
	     case  6:   /* (no) right preconditioner transposed solver */
	       //   printf("apply prect vec\n"); fflush(STDOUT);
	       i=1;
	       if (param.ipar[21]==2)
		  AMGTSOL(PRE, nlev, &param, w+param.ipar[27],w+param.ipar[28], param.dbuff);
	       else 
		 if (param.ipar[21]==3)
		    AMGTDUSOL(PRE, nlev, &param, w+param.ipar[27],w+param.ipar[28], param.dbuff);
		 else
		    COPY(&n,w+param.ipar[27],&i,w+param.ipar[28],&i);
	       // printf("prect vec applied\n"); fflush(STDOUT);
	       break;
	     case 10:   /* call my own stopping test routine */
	       break;
	     default:   /* the iterative solver terminated with code=param.ipar[20] */
	       flags=0;
	       break;
	     } /* end switch */
	   w++;
	 } /* end while */
       for (i=0; i<n; i++) { 
#if defined _DOUBLE_REAL_ || defined _SINGLE_REAL_
	   sol[i]*=rightscale[i];
#else
	   val=sol[i].r;
	   sol[i].r=val*rightscale[i].r-sol[i].i*rightscale[i].i;
	   sol[i].i=val*rightscale[i].i+sol[i].i*rightscale[i].r;
#endif
       } // end for i
       evaluate_time(&time_stop,&systime);
       secnds=time_stop-time_start;



       /* why did the iterative solver stop? */
       switch (param.ipar[20])
	 {
	 case  0:  /* everything is fine */
	   break;
	 case -1:  /* too many iterations */
	   printf("number of iteration steps exceeds its limit\n");
	   break;
	 case -2:  /* not enough work space */
	   printf("not enough work space provided\n");
	   break;
	 case -3:  /* not enough work space */
	   printf("algorithm breaks down ");
	   /* Arnoldi-type solvers */
	   if (SOLVER>=7)
	     {
	       switch (param.ipar[31])
		 {
		 case -1:  
		   printf("due to zero input vector");
		   break;
		 case -2:  
		   printf("since input vector contains abnormal numbers");
		   break;
		 case -3:  
		   printf("since input vector is a linear combination of others");
		   break;
		 case -4:  
		   printf("since triangular system in GMRES/FOM/etc. has null rank");
		   break;
		 } /* end switch */
	     } /* end if */
	   printf("\n");
	   break;
	 default:  /* unknown */
	   printf("solver exited with error code %d\n",param.ipar[20]);
	 } /* end switch */

       printf("number of iteration steps: %d\n",param.ipar[26]);
       printf("time: %8.1le [sec]\n",  (double)secnds);
       printf("      %8.1le [sec]\n\n",(double)ILUPACK_secnds[5]); 
       
       if (param.ipar[26]>=MAX_IT) {
	  fprintf(fo,"  inf  |");
	  fprintf(fo," inf\n");
       }
       else 
	 if (param.ipar[20]!=0) {
	    fprintf(fo,"  NaN  |");
	    fprintf(fo," NaN\n");
	 }
	 else {
	    fprintf(fo,"%7.1le|",(double)secnds);
	    fprintf(fo," %3d\n",param.ipar[26]);
	 }
       fflush(fo);
       
       printf("time matrix-vector multiplication: %8.1le [sec]\n",  (double)secndsAx/((double)nAx));
       printf("                                   %8.1le [sec]\n\n",(double)ILUPACK_secnds[6]);
       printf("residual norms:\n");
       printf("initial: %8.1le\n",param.fpar[22]);
       printf("target:  %8.1le\n",param.fpar[23]);
       printf("current: %8.1le\n",param.fpar[25]);
       printf("convergence rate: %8.1le\n",param.fpar[26]);


       if (nrhs==0) {
	 for (i=0; i<n; i++)
#if defined _DOUBLE_REAL_ || defined _SINGLE_REAL_
	   sol[i]-=1;
#else
	   sol[i].r-=1;
#endif
	 i=1;
	 if ((rhstyp[2]=='X' || rhstyp[2]=='x') || (rhstyp[0]!='M' && rhstyp[0]!='m'))
	    printf("rel. error in the solution: %8.1le\n\n",NRM(&n,sol,&i)/sqrt(1.0*n));
	 else printf("\n");
       }
       else {
	  /* rhs */
	  m=1;
	  if (rhstyp[1]=='G' || rhstyp[1]=='g') m++;

	  if (rhstyp[0]=='M' || rhstyp[0]=='m')
	     rhs+=n+(m-1)*n*nrhs;
	  else
	     rhs+=n*m*nrhs;
	  if (rhstyp[2]=='X' || rhstyp[2]=='x') {
	     for (i=0; i<n; i++) {
#if defined _DOUBLE_REAL_ || defined _SINGLE_REAL_
	         sol[i]-=rhs[i];
#else
		 sol[i].r-=rhs[i].r;
		 sol[i].i-=rhs[i].i;
#endif
	     } // end for i
	  } // end if

	  i=1;
	  if (rhstyp[2]=='X' || rhstyp[2]=='x')
	     printf("rel. error in the solution: %8.1le\n\n",NRM(&n,sol,&i)/NRM(&n,rhs,&i));
	  else printf("\n");

	  if (rhstyp[0]=='M' || rhstyp[0]=='m')
	     rhs-=n+(m-1)*n*nrhs;
	  else
	     rhs-=n*m*nrhs;
       }
       free(w);

       // advance rhs (set up for multiple rhs loops)
       if (rhstyp[0]=='M' || rhstyp[0]=='m')
	  rhs+=n;
#endif /* USE_SPARSKIT */


       printf("\n\n");
    }
    else
       printf("Iterative solver(s) cannot be applied\n\n");


    /* ------------------   END apply preconditioner   ------------------- */








    /* ------------------------------------------------------------------- */
    /* ------------------------   END MAIN part   ------------------------ */
    /* ------------------------------------------------------------------- */
    fclose(fo);
} /* end main */