main.c

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

    fprintf(fo,"mc64af/PQ|");
#elif defined MC64_AMD_PQ
    perm0=PERMMC64AMD;
    perm =PERMAMD;
    permf=PERMPQ;
    fprintf(fo,"mc64ad/PQ|");
#elif defined MC64_METIS_E_PQ
    perm0=PERMMC64METISE;
    perm =PERMMETISE;
    permf=PERMPQ;
    fprintf(fo,"mc64me/PQ|");
#elif defined MC64_METIS_N_PQ
    perm0=PERMMC64METISN;
    perm =PERMMETISN;
    permf=PERMPQ;
    fprintf(fo,"mc64mn/PQ|");
#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; // DROP_TOL; // 
    droptols[1]=DROP_TOL; 

    // choose the iterative solver, default: 8 (GMRES)
    param.ipar[5]=SOLVER;

    // 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-12;

    // turn of column scaling
    //param.ipar[7]&=~(2+16+128);
    
    // use diagonal compensation
    //flags|=DIAGONAL_COMPENSATION;
    //flags|=DIAGONAL_COMPENSATION|TISMENETSKY_SC;
    //flags|=TISMENETSKY_SC;

    // use Standard Schur-complement
    //flags&=~SIMPLE_SC;

    // turn off aggressive dropping
    flags&=~AGGRESSIVE_DROPPING;

    // turn off destroying the coarse grid matrix
    //flags&=~DISCARD_MATRIX;

    // 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;

    // limit maximum number of entries in each column of L, row of U
    //param.ipar[3]=ELBOW*(nz/n)+5;
    // limit maximum number of entries in each row of S
    //param.ipar[9]=ELBOW*(nz/n)+5;

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

    // manually change the drop tolerance for the approximate Schur complement
    //param.fpar[8]=DROP_TOL;


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

    evaluate_time(&time_start,&systime);
    ierr=AMGFACTOR(&A, &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! */
	            printf("factorization successful completed\n");
	            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
    */
    AT=A;
    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]=3
	*/
	ierr=AMGSOLVER(A, 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
	      printf("iteration successful completed\n");
	      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|");
	      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|");
	      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",param.fpar[23]);

    
	/* 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,&(A.nr),&mynrhs,
		 "solution vectors computed by ILUPACK                                    ",
		 "sol     ",
		 j, 72, 8);

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

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

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