maingnlsym.c

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

#if !defined RCM_PP && !defined AMF_PP && !defined AMD_PP && !defined MMD_PP && !defined FC_PP && !defined METIS_E_PP && !defined METIS_N_PP 
    flags&=~FINAL_PIVOTING;
#endif
    // change flags if mpiluc is desired
#ifdef USE_MPILUC
    flags|=MULTI_PILUC;
#endif

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

    // choose the iterative solver, default: 8 (GMRES)
    param.ipar[5]=SOLVER;
    // number of restart length for GMRES,FOM,... 
    param.ipar[24]=30;

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

    // overwrite default values for condest
    condest=CONDEST;

    // do not discard matrix entries on an early stage
    //flags&=~AGGRESSIVE_DROPPING;
    //flags&=~DISCARD_MATRIX;

    // rewrite the updated parameters
    MYSYMAMGSETPARAMS(AS, &param, flags, elbow, droptols, condest,
		      restol, max_it);

    // 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 "symmessages.c"
       
    evaluate_time(&time_start,&systime);
    ierr=MYSYMAMGFACTOR(&AS, &PRE, &nlev, &param, perm0,perm, permf);
    // explicitly rescale original matrix
    for (i=0; i<A.nr; i++)
        for (j=A.ia[i]-1; j<A.ia[i+1]-1; j++) {
	    k=A.ja[j]-1;
#if defined _DOUBLE_REAL_ || defined _SINGLE_REAL_
	    A.a[j]*=PRE.colscal[i]*PRE.colscal[k];
#else
	    A.a[j].r*=PRE.colscal[i].r*PRE.colscal[k].r;
	    A.a[j].i*=PRE.colscal[i].r*PRE.colscal[k].r;
#endif
	}
  
    // 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);
    }


#ifdef PRINT_INFO
    next=&PRE;
    for (i=1; i<=nlev; i++) {
        // fill-in LU
        printf("level %3d, block size %7d\n",i,next->LU.nr); fflush(stdout);
	if (i<nlev || next->LU.ja!=NULL) {
	  printf("U-factor");
	  printf("\n");fflush(stdout);
	  for (l=0; l<next->LU.nr; ) {
	    if (next->LU.ja[next->LU.nr+1+l]==0){
	      for (j=next->LU.ja[l];j<next->LU.ja[l+1]; j++) {
		printf("%8d",next->LU.ja[j-1]);
	      }
	      printf("\n");fflush(stdout);
	      for (j=next->LU.ja[l];j<next->LU.ja[l+1]; j++) {
		printf("%8.1le",next->LU.a[j-1]);
	      }
	      l++;
	    }
	    else {
	      for (j=next->LU.ja[l];j<next->LU.ja[l+1]; j++) {
		printf("%8d",next->LU.ja[j-1]);
	      }
	      printf("\n");fflush(stdout);
	      for (j=next->LU.ja[l];j<next->LU.ja[l+1]; j++) {
		printf("%8.1le",next->LU.a[next->LU.ja[l]+2*(j-next->LU.ja[l])-1]);
	      }
	      printf("\n");fflush(stdout);
	      for (j=next->LU.ja[l];j<next->LU.ja[l+1]; j++) {
		printf("%8.1le",next->LU.a[next->LU.ja[l]+2*(j-next->LU.ja[l])]);
	      }
	      l+=2;
	    }
	    printf("\n");fflush(stdout);
	  }

	  printf("Block diagonal factor\n");
	  for (l=0; l<next->LU.nr;) {
	    if (next->LU.ja[next->LU.nr+1+l]==0){
	      printf("%8.1le",next->LU.a[l]);
	      l++;
	    }
	    else {
	      printf("%8.1le%8.1le",next->LU.a[l],next->LU.a[next->LU.nr+1+l]);
	      l+=2;
	    }
	  }
	  printf("\n");fflush(stdout);
	  for (l=0; l<next->LU.nr; ) {
	    if (next->LU.ja[next->LU.nr+1+l]==0) {
	      printf("        ");
	      l++;
	    }
	    else {
	      printf("%8.1le%8.1le",next->LU.a[next->LU.nr+1+l],next->LU.a[l+1]);
	      l+=2;
	    }
	  }
	  printf("\n");fflush(stdout);



	}
	if (i==nlev) {
	   if (next->LU.ja==NULL) {
	      printf("switched to full matrix processing\n");fflush(STDOUT);
	   }
	}

	if (i<nlev) {
	   // fill-in F
	   nnzU+=next->F.ia[next->F.nr]-1;
	   printf("level %3d->%3d, block size (%7d,%7d)\n",i,i+1,next->LU.nr,next->F.nc);
	   printf("  local fill-in F %7d(%5.1lfav pr)\n",
		  next->F.ia[next->F.nr]-1,(1.0*(next->F.ia[next->F.nr]-1))/next->LU.nr);
	}
	next=next->next;
    }
#endif



    // print some statistics about the levels, their size and the 
    // computation time
#include "symprintperformance.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=MYGNLSYMAMGSOLVER(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
	      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;
    }
    extension[0]=fname[++i]; fname[i]='s';
    extension[1]=fname[++i]; fname[i]='o';
    extension[2]=fname[++i]; 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';
    */

    // outputfile name, rewrite ".rsa,..." extension
    i=fnamelen;
    j=-1;
    while (j) {
          i--;
	  if (i<0)
	    j=0;
	  if (fname[i]=='.')
	    j=0;
    }
    fname[++i]=extension[0];
    fname[++i]=extension[1];
    fname[++i]=extension[2];
    i++;
    j=i;
    while (i<100)
      fname[i++]='\0';


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

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

}