maint.c

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

    permf=PERMPQ;
    fprintf(fo,"mc64ad/PQ|");
    //printf("prescribe MC64 ordering\n");
#elif defined MC64_METIS_E_PQ
    perm0=PERMMC64METISE;
    perm =PERMMETISE;
    permf=PERMPQ;
    fprintf(fo,"mc64me/PQ|");
    //printf("prescribe MC64 ordering\n");
#elif defined MC64_METIS_N_PQ
    perm0=PERMMC64METISN;
    perm =PERMMETISN;
    permf=PERMPQ;
    fprintf(fo,"mc64mn/PQ|");
    //printf("prescribe MC64 ordering\n");
#elif defined MWM_RCM_PQ
    perm0=PERMMWMRCM;
    perm =PERMRCM;
    permf=PERMPQ;
    fprintf(fo,"mwrc/PQs |");
#elif defined MWM_MMD_PQ
    perm0=PERMMWMMMD;
    perm =PERMMMD;
    permf=PERMPQ;
    fprintf(fo,"mwmd/PQs |");
#elif defined MWM_AMF_PQ
    perm0=PERMMWMAMF;
    perm =PERMAMF;
    permf=PERMPQ;
    fprintf(fo,"mwaf/PQs |");
#elif defined MWM_AMD_PQ
    perm0=PERMMWMAMD;
    perm =PERMAMD;
    permf=PERMPQ;
    fprintf(fo,"mwad/PQs |");
#elif defined MWM_METIS_E_PQ
    perm0=PERMMWMMETISE;
    perm =PERMMETISE;
    permf=PERMPQ;
    fprintf(fo,"mwme/PQs |");
#elif defined MWM_METIS_N_PQ
    perm0=PERMMWMMETISN;
    perm =PERMMETISN;
    permf=PERMPQ;
    fprintf(fo,"mwmn/PQs |");
#else
    fprintf(fo,"    /    |");
#endif


    // modify the default settings
    AMGGETPARAMS(param, &flags, &elbow, droptols, &condest,
		 &restol, &max_it, &nrestart);

    // ONLY for mixed reordering strategies it is useful to
    // set the 'FINAL_PIVOTING' flag
#if !defined RCM_PQ && !defined AMF_PQ && !defined AMD_PQ && !defined MMD_PQ && !defined FC_PQ && !defined METIS_E_PQ && !defined METIS_N_PQ && !defined MWM_RCM_PQ && !defined MWM_MMD_PQ && !defined MWM_AMF_PQ && !defined MWM_AMD_PQ && !defined MWM_METIS_E_PQ && !defined MWM_METIS_N_PQ && !defined MC64_RCM_PQ && !defined MC64_MMD_PQ && !defined MC64_AMF_PQ && !defined MC64_AMD_PQ && !defined MC64_METIS_E_PQ && !defined MC64_METIS_N_PQ
    flags&=~FINAL_PIVOTING;
#endif

    // change flags if mpiluc is desired
#ifdef USE_MPILUC
    flags|=MULTI_PILUC;
#endif

    // overwrite the default drop tolerances
    droptols[0]=1.0;
    droptols[1]=DROP_TOL; 

    // overwrite the default value for elbow space
    elbow=ELBOW;

    // overwrite default values for condest
    condest=CONDEST;

    // overwrite the default value for the residual norm
    restol=1e-8;

    // use Standard Schur-complement
    flags|=SIMPLE_SC;

    // keep the (1,2) block and the (1,2) block in the multilevel ILU
    // factorization, this may lead to more fill-in
    // flags&=~COARSE_REDUCE;
 
    // interchange order of scalings for all three cases
    param.ipar[7]|=4+32+256;

    // rewrite the updated parameters
    AMGSETPARAMS(AT, &param, flags, elbow, droptols, condest,
		 restol, max_it, nrestart);


    // print some messages that give information about flags and reorderings
#include "messages.c"


    evaluate_time(&time_start,&systime);
    ierr=AMGFACTOR(&AT, &PRE, &nlev, &param, perm0,perm, permf);
    // update buffer size
    nibuff=param.nibuff;
    ndbuff=param.ndbuff;
    evaluate_time(&time_stop,&systime);
    secnds=time_stop-time_start;

    switch (ierr)
    {
           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 at level %d\n",nlev);
		    fprintf(fo,"input matrix may be wrong\n");fclose(fo); 
		    break;
           case -2: /* The matrix L overflows the array alu */
	            printf("The matrix L overflows the array alu at level %d\n",nlev);
		    fprintf(fo,"out of memory\n");fclose(fo); 
		    break;
           case -3: /* The matrix U overflows the array alu */
	            printf("The matrix U overflows the array alu at level %d\n",nlev);
		    fprintf(fo,"out of memory\n");fclose(fo); 
		    break;
           case -4: /* Illegal value for lfil */
	            printf("Illegal value for lfil at level %d\n",nlev);
		    fprintf(fo,"Illegal value for lfil\n");fclose(fo); 
		    break;
           case -5: /* zero row encountered */
	            printf("zero row encountered at level %d\n",nlev);
		    fprintf(fo,"zero row encountered\n");fclose(fo); 
		    break;
           case -6: /* zero column encountered */
	            printf("zero column encountered at level %d\n",nlev);
		    fprintf(fo,"zero column encountered\n");fclose(fo); 
		    break;
           case -7: /* buffers too small */
	            printf("buffers are too small\n");
		    // This error message would not be necessary if AMGsetup is
		    // called with the correct values of nibuff, ndbuff
		    printf("increase buffer size to at least %ld (float), %ld (int)\n",
			   ndbuff, nibuff);
		    fflush(stdout);
		    fprintf(fo,"buffers are too small\n");fclose(fo); 
           default: /* zero pivot encountered at step number ierr */
	            printf("zero pivot encountered at step number %d of level %d\n",ierr,nlev);
		    fprintf(fo,"zero pivot encountered\n");fclose(fo); 
		    break;
    } /* end switch */
    if (ierr) {
       fprintf(fo,"Iterative solver(s) cannot be applied\n");
       fflush(fo);
       exit(ierr);
    }

    // print some statistics about the levels, their size and the 
    // computation time
#include "printperformance.c"

    
    /* some decisions about the right hand side, the exact solution and the 
       initial guess

       - if a right hand side is provided, then solve the system according
         to this given right hand side.
	 Otherwise set x=(1,...,1)^T and use b=Ax as right hand side
	 --> rhs

       - if an initial guess is provided, then use this as starting vector.
         Otherwise use x=0 as initial guess
	 --> sol
       
       - for statistics also compute the norm of the initial residual --> val
    */
    sumtime=0.0;
    sumit=0;
    for (l=0; l<mynrhs; l++) {
#include "initvectors.c"


        evaluate_time(&time_start,&systime);
	/* left  preconditionig:                   param.ipar[21]=1

	   right preconditionig:                   param.ipar[21]=2

	   double sided preconditionig, 
	   L^{-1} as left preconditioner and 
           (DU)^{-1} as right preconditioner:      param.ipar[21]=2
	*/
	// transposed system is stored
	param.ipar[4]=1; 
	ierr=AMGSOLVER(AT, PRE, nlev, &param, sol+A.nr*l, rhs+A.nr*l);
	evaluate_time(&time_stop,&systime);
	secnds=time_stop-time_start;

	// why did the iterative solver stop?
	switch (ierr) {
	case  0:  // everything is fine
	      sumtime+=ILUPACK_secnds[5];
	      sumit+=param.ipar[26];
	      if (l==mynrhs-1) {
		 sumtime/=mynrhs;
		 sumit   =sumit/((double)mynrhs)+.5;
		 fprintf(fo,"%7.1le|",(double)sumtime);
		 fprintf(fo," %3d|",sumit);
		 fprintf(fo,"%7.1le\n",(double)sumtime+ILUPACK_secnds[7]);
	      }
	      break;
	case -1:  // too many iterations
	      printf("number of iteration steps exceeds its limit\n");
	      fprintf(fo,"  inf  |");
	      fprintf(fo," inf\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 ");
	      printf("\n");
	      fprintf(fo,"  NaN  |");
	      fprintf(fo," NaN\n");
	      break;
	default:  /* unknown */
              printf("solver exited with error code %d\n",ierr);
	} // end switch 
	fflush(fo);

	// stop if necessary
	if (ierr)
	   mynrhs=l;

	printf("%3d. right hand side\n",l+1);
	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]); 
       
	printf("time matrix-vector multiplication: %8.1le [sec]\n",ILUPACK_secnds[6]);

	printf("residual norms:\n");
	printf("initial: %8.1le\n",val);
	printf("target:  %8.1le\n",val*restol);
	fflush(stdout);
    

	/* some final statistics 
	   - about the true current residual of the computed solution and
	   - the relative error in the solution though the exact solution is known
	*/
#include "finalres.c"
    } // end for l
    fclose(fo);
 
    // outputfile name, add ".sol" extension
    i=fnamelen;
    j=-1;
    while (j) {
          i--;
	  if (i<0)
	    j=0;
	  if (fname[i]=='.')
	    j=0;
    }
    fname[++i]='s';
    fname[++i]='o';
    fname[++i]='l';
    i++;
    j=i;
    while (i<100)
      fname[i++]='\0';

    WRITEVECTORS(fname,sol,&(AT.nr),&mynrhs,
		 "solution vectors computed by ILUPACK                                    ",
		 "sol     ",
		 j, 72, 8);

    /*
      // sample code for rereading the vectors from C
      READVECTORS(fname,sol,&(AT.nr),&mynrhs,title,key,j, 72, 8);
      title[72]='\0';
      key[8]='\0';
    */

    // finally release memory of the preconditioner
    AMGDELETE(AT,PRE,nlev,&param);

    // release memory used for the input matrix
    free(AT.ia);
    free(AT.ja);
    free(AT.a );
    free(rhs );
} // end maint