mainildlc.c

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

    else
       printf("   dual threshold dropping\n");
    if (param & TISMENETSKY_SC)
       printf("   Tismenetsky-like Schur complement\n");
    else
       printf("   simple Schur complement\n");
    if (param & NO_SHIFT)
       printf("   zero pivots are kept\n");
    else
       printf("   zero pivots are shifted away\n");
    if (param & REPEAT_FACT)
       printf("   second pass\n");
    else
       printf("   first pass\n");
    if (param & IMPROVED_ESTIMATE)
       printf("   improved estimate of the inverse\n");
    else
       printf("   simple estimate of the inverse\n");
    if (param & DIAGONAL_COMPENSATION)
       printf("   diagonal compensation\n");
    else
       printf("   diagonal entries unmodified\n");
    if (param & ENSURE_SPD)
       printf("   SPD property preserved\n");
    else
       printf("   no attempt to preserve SPD\n");
    fflush(stdout);

    evaluate_time(&time_start,&systime);
    ILDLC(&n,A.a,A.ja,A.ia,ipar,fpar,&param,
	 ilutmat.a,ilutmat.ja,ipar+1,w,ws,ipar+2);
    evaluate_time(&time_stop,&systime);
    secnds=time_stop-time_start;
    fprintf(fo,"   |");
    free(w);
    free(ws);

    /*
    printf("diagonal entries\n");
    for (i=0; i<n; i++)
      printf("%8.1e",ilutmat.a[i]);
      // printf("%8.1e%8.1e",ilutmat.a[i].r,ilutmat.a[i].i);
    printf("\n");fflush(stdout);
    printf("U entries\n");
    for (i=0; i<n; i++) {
      printf("%8d:\n",i+1);
      for (j=ilutmat.ja[i]; j<ilutmat.ja[i+1]; j++)
	printf("%16d",ilutmat.ja[j-1]);
      printf("\n");
      for (j=ilutmat.ja[i]; j<ilutmat.ja[i+1]; j++)
	printf("%8.1e",ilutmat.a[j-1]);
        // printf("%8.1e%8.1e",ilutmat.a[j-1].r,ilutmat.a[j-1].i);
      printf("\n");
    }
    fflush(stdout);
    */

    nc=0;
    for (i=0; i<n; i++)
        nc+=ilutmat.ja[i+1]-ilutmat.ja[i];
    printf("\nfill-in U: %5d(%4.1fav.), totally %5d(%4.1fav.)\n",
	   nc,(1.0*nc)/n,nc+n,(1.0*(nc+n))/n);
    printf("fill-in factor: %5.1f\n",(1.0*(nc+n))/nz);
    fprintf(fo,"%5.1f|",(1.0*(tmp3+nc+n))/nz);
    printf("time: %8.1e [sec]\n\n",(double)secnds);
    fprintf(fo,"%7.1le|",(double)secnds+secndsprep);
    switch (ipar[2]) {
           case  0: /* perfect! */
	            break;
           case -1: /* Error. input matrix may be wrong.
                       (The elimination process has generated a
			row in L or U whose length is .gt.  n.) */
	            printf("Error. input matrix may be wrong\n");
		    break;
           case -2: /* The matrix L overflows the array alu */
	            printf("The matrix L overflows the array alu\n");
		    break;
           case -3: /* The matrix U overflows the array alu */
	            printf("The matrix U overflows the array alu\n");
		    break;
           case -4: /* Illegal value for lfil */
	            printf("Illegal value for lfil\n");
		    break;
           case -5: /* zero row encountered */
	            printf("zero row encountered\n");
		    break;
           default: /* zero pivot encountered at step number ierr */
	            printf("zero pivot encountered at step number %d\n",ipar[2]);
		    break;
    } // end switch


    if (ipar[2]!=0) {
	printf("Iterative solver(s) cannot be applied\n\n");
	fprintf(fo,"Iterative solver(s) cannot be applied\n");
	fflush(fo);
	exit(0);
    } // end if

    /* --------------------   apply preconditioner   --------------------- */
    // 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
    
    // rescale initial solution
    for (i=0; i<n; i++) {  
#if defined _DOUBLE_REAL_ || defined _SINGLE_REAL_
        sol[i]/=scale[i];
#else
	sol[i].r/=scale[i].r;
	sol[i].i/=scale[i].r;
#endif
    } // end for i


    // init solver
    ipar[0]=0;
    // right preconditioning
    ipar[1]=2;
    if (ipar[1]==0)
      printf("no preconditioning\n");
    else if (ipar[1]==1)
      printf("left preconditioning\n");
    else if (ipar[1]==2)
      printf("right preconditioning\n");
    else if (ipar[1]==3)
      printf("both sides preconditioning\n");
       
    // stopping criteria, energy norm
    ipar[2]=3;
    // number of restart length for GMRES,FOM,...
    ipar[4]=30;
    // maximum number of iteration steps
    ipar[5]=1.1*n+10;
    ipar[5]=MIN(ipar[5],MAX_IT);
    // init with zeros
    ipar[6]=ipar[7]=ipar[8]=ipar[9]=ipar[10]=ipar[11]=ipar[12]=ipar[13]=0;
    // relative error tolerance
    fpar[0]=1e-8;
    // absolute error tolerance
    fpar[1]=0;
    // init with zeros
    fpar[2]=fpar[3]=fpar[4]=fpar[5]=fpar[6]=fpar[7]=fpar[8]=fpar[9]=fpar[10]=0;
       
    // allocate memory for iterative solver depending on our choice */
    i=ipar[4];
    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",ipar[4]);
      break;
    case  8:   // gmres
      i=(n+3)*(i+2)+(i+1)*i/2;
      printf("use GMRES(%d) as iterative solver\n",ipar[4]);
      break;
    case  9:   // fgmres
      i=2*n*(i+1)+(i+1)*i/2;
      printf("use FGMRES(%d) as iterative solver\n",ipar[4]);
      break;
    case 10:   // dqgmres
      i=n+(i+1)*(2*n+4);
      printf("use DQGMRES(%d) as iterative solver\n",ipar[4]);
      break;
    } // end switch
    w=(FLOAT *)MALLOC((size_t)i*sizeof(FLOAT),"mainildlc:w");
    ipar[3]=i;
       
    // start iterative solver
    evaluate_time(&time_start,&systime);
    flag=-1;
    while (flag) {
          // which solver do we use?
	  switch (SOLVER) {
	  case  1:   // pcg
	    PCG(&n,rhs,sol,ipar,fpar,w);
	    break;
	    /*
	     case  2:   // cgnr
	       cgnr(&n,rhs,sol,ipar,fpar,w);
	       break;
	     case  3:   // bcg
	       bcg(&n,rhs,sol,ipar,fpar,w);
	       break;
	     case  4:   // dbcg
	       dbcg(&n,rhs,sol,ipar,fpar,w);
	       break;
	     case  5:   // bcgstab
	       bcgstab(&n,rhs,sol,ipar,fpar,w);
	       break;
	     case  6:   // tfqmr
	       tfqmr(&n,rhs,sol,ipar,fpar,w);
	       break;
	     case  7:   // fom 
	       fom(&n,rhs,sol,ipar,fpar,w);
	       break;
	     case  8:   // gmres
	       GMRES(&n,rhs,sol,ipar,fpar,w);
	       break;
	     case  9:   // fgmres
	       fgmres(&n,rhs,sol,ipar,fpar,w);
	       break;
	     case 10:   // dqgmres
	       dqgmres(&n,rhs,sol,ipar,fpar,w);
	       break;
	    */
	  } // end switch
	  // what is needed for the solver?
	  w--;
	  switch (ipar[0]) {
	  case  1:   // apply matrix vector multiplication
	    HERMATVEC(A,w+ipar[7],w+ipar[8]);
	    break;
	  case  2:   // apply transposed matrix vector multiplication
	    HERMATVEC(A,w+ipar[7],w+ipar[8]);
	    break;
	  case  3:   // (no) left preconditioner solver
	    i=1;
	    if (ipar[1]==1)
	       ILDLCSOL(&n, w+ipar[7],w+ipar[8], ilutmat.a,ilutmat.ja);
	    else
	       COPY(&n,w+ipar[7],&i,w+ipar[8],&i);
	    break;
	  case  4:   // (no) left preconditioner transposed solver
	    i=1;
	    if (ipar[1]==1)
	       ILDLCSOL(&n, w+ipar[7],w+ipar[8], ilutmat.a,ilutmat.ja);
	    else
	       COPY(&n,w+ipar[7],&i,w+ipar[8],&i);
	    break;
	  case  5:   // (no) right preconditioner solver
	    i=1;
	    if (ipar[1]==2)
	       ILDLCSOL(&n, w+ipar[7],w+ipar[8], ilutmat.a,ilutmat.ja);
	    else
	       COPY(&n,w+ipar[7],&i,w+ipar[8],&i);
	    break;
	  case  6:   // (no) right preconditioner transposed solver
	    i=1;
	    if (ipar[1]==2)
	       ILDLCSOL(&n, w+ipar[7],w+ipar[8], ilutmat.a,ilutmat.ja);
	    else
	       COPY(&n,w+ipar[7],&i,w+ipar[8],&i);
	    break;
	  case 10:   /* call my own stopping test routine */
	    break;
	  default:   /* the iterative solver terminated with code=ipar[0] */
	    flag=0;
	    break;
	  } // end switch
	  w++;
    } // end while

    // rescale solution
    for (i=0; i<n; i++) {  
#if defined _DOUBLE_REAL_ || defined _SINGLE_REAL_
        sol[i]*=scale[i];
#else
	sol[i].r*=scale[i].r;
	sol[i].i*=scale[i].r;
#endif
    } // end for i
    evaluate_time(&time_stop,&systime);
    secnds=time_stop-time_start;
    // why did the iterative solver stop?
    switch (ipar[0]) {
    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 (ipar[11]) {
	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",ipar[0]);
    } // end switch

    printf("number of iteration steps: %d\n",ipar[6]);
    printf("time: %8.1le [sec]\n\n", (double)secnds);
    if (ipar[6]>=MAX_IT) {
       fprintf(fo,"  inf  |");
       fprintf(fo," inf\n");
    } // end if
    else 
       if (ipar[0]!=0) {
	  fprintf(fo,"  NaN  |"); 
	  fprintf(fo," NaN\n");
       } // end if 
       else {
	  fprintf(fo,"%7.1le|",(double)secnds);
	  fprintf(fo," %3d\n",ipar[6]);
       } // if-else

    fflush(fo);
    
    printf("residual norms:\n");
    printf("initial: %8.1e\n",fpar[2]);
    printf("target:  %8.1e\n",fpar[3]);
    printf("current: %8.1e\n",fpar[5]);

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

    free(w);



    /* ------------------------------------------------------------------- */
    /* ------------------------   END MAIN part   ------------------------ */
    /* ------------------------------------------------------------------- */
} // end main