mainsym3.c

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

	      }
	      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
    */
    sumtime=0.0;
    sumit=0;
    for (l=0; l<mynrhs; l++) {
#include "syminitvectors.c"


        evaluate_time(&time_start,&systime);
	ierr=MYSYMAMGSOLVER(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 "symfinalres.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 );



    // ------------------------------------------------------------------------
    // read in another matrix, here: same matrix for the second time
    // note that we keep the RE_FACTOR flag to indicate that we want
    // to solve a third system
    // ------------------------------------------------------------------------

    
    // save data from old structures.
    // This is only necessary if the size of the matrix changes
    nibuff=param.nibuff; ibuff=param.ibuff;
    ndbuff=param.ndbuff; dbuff=param.dbuff;
    nju   =param.nju;    ju   =param.ju;
    njlu  =param.njlu;   jlu  =param.jlu;
    nalu  =param.nalu;   alu  =param.alu;

#include "spdreadmatrix.c"

    
    // if right hand sides are provided, then run AMGSOLVER for any of these
    // right hand sides. Otherwise use own set of right hand sides

    // allocate memory for the solution vector and some buffer
    sol  =(FLOAT *) MALLOC(mynrhs*(size_t)n*sizeof(FLOAT), "mainspd:sol");


    // set parameters to the default settings. This is only necessary, if the
    // size of the input matrix has changed.
    // Otherwise one can skip AMGINIT and one does not have to save the buffer
    // pointers
    MYSYMAMGINIT(A, &param);
    param.nibuff=nibuff; param.ibuff=ibuff;
    param.ndbuff=ndbuff; param.dbuff=dbuff;
    param.nju=nju;       param.ju =ju;
    param.njlu=njlu;     param.jlu=jlu;
    param.nalu=nalu;     param.alu=alu;
    
    // display selected reordering functions
#ifdef MINIMUM_DEGREE
    fprintf(fo,"mmds/mmds|");
#elif defined REVERSE_CM
    fprintf(fo,"rcms/rcms|");
#elif defined NESTED_DISSECTION
    fprintf(fo,"nds /nds |");
#elif defined IND_SET
    fprintf(fo,"inds/inds|");
#elif defined AMF
    fprintf(fo,"amfs/amfs|");
#elif defined AMD
    fprintf(fo,"amds/amds|");
#elif defined PP_PERM
    fprintf(fo,"PPs /PPs |");
#elif defined MMD_PP
    fprintf(fo,"mmds/PPs |");
#elif defined AMF_PP
    fprintf(fo,"amfs/PPs |");
#elif defined AMD_PP
    fprintf(fo,"amds/PPs |");
#elif defined RCM_PP
    fprintf(fo,"rcms/PPs |");
#elif defined FC_PP
    fprintf(fo,"FCs /PPs |");
#elif defined METIS_E_PP
    fprintf(fo,"mes /PQs |");
#elif defined METIS_N_PP
    fprintf(fo,"mns /PPs |");
#elif defined SMWM_MMD_PP
    fprintf(fo,"mwmd/PPs |");
#elif defined SMWM_AMF_PP
    fprintf(fo,"mwaf/PPs |");
#elif defined SMWM_AMD_PP
    fprintf(fo,"mwad/PPs |");
#elif defined SMWM_RCM_PP
    fprintf(fo,"mwrc/PPs |");
#elif defined SMWM_METIS_E_PP
    fprintf(fo,"mwme/PQs |");
#elif defined SMWM_METIS_N_PP
    fprintf(fo,"mwmn/PPs |");
#elif defined SMC64_MMD_PP
    fprintf(fo,"mc64md/PP|");
#elif defined SMC64_AMF_PP
    fprintf(fo,"mc64af/PP|");
#elif defined SMC64_AMD_PP
    fprintf(fo,"mc64ad/PP|");
#elif defined SMC64_RCM_PP
    fprintf(fo,"mc64rc/PP|");
#elif defined SMC64_METIS_E_PP
    fprintf(fo,"mc64me/PP|");
#elif defined SMC64_METIS_N_PP
    fprintf(fo,"mc64mn/PP|");
#else
    fprintf(fo,"    /    |");
#endif


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

    // ONLY for mixed reordering strategies it is useful to
    // set the 'FINAL_PIVOTING' flag
#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[0]=1.0; // DROP_TOL; //
    droptols[1]=DROP_TOL; 

    // choose the iterative solver, default: 4 (symmetric QMR)
    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-8;

    // turn of scaling
    // param.ipar[7]&=~(3+24+192);
    // param.ipar[8]&=~((1+32+1024)*4);
    // param.ipar[8]|=((1+32+1024)*5);

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

    // indicate that we want to compute more than one ILU
    flags|=RE_FACTOR;

    // rewrite the updated parameters
    MYSYMAMGSETPARAMS(A, &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(&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! */
	            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 "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
    */
    sumtime=0.0;
    sumit=0;
    for (l=0; l<mynrhs; l++) {
#include "syminitvectors.c"


        evaluate_time(&time_start,&systime);
	ierr=MYSYMAMGSOLVER(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|");