mainexpert.c

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

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

    // overwrite the default value for the maximum number of iterations
    // max_it=MIN(1.1*n+10,MAX_IT);

    // overwrite the default value for the restart
    // nrestart=30;

    // use diagonal compensation
    //flags|=DIAGONAL_COMPENSATION;
    //flags|=DIAGONAL_COMPENSATION|TISMENETSKY_SC;
    //flags|=TISMENETSKY_SC;
    flags|=SIMPLE_SC;

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


    // 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! */
	            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 %d (float), %d (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 USE_MPILUC
    // printf("ARMSMPILUC finished\n");fflush(STDOUT);

    printf("ARMS oo MPILUC, multilevel structure\n");
#else
    // printf(" ARMSPILUC finished\n");fflush(STDOUT);

    printf("ARMS oo  PILUC, multilevel structure\n");
#endif

    fprintf(fo,"%3d|",nlev);

    next=⪯
    nnzL=0;
    nnzU=0;
    tmp=0;
    for (i=1; i<=nlev; i++) {
        // fill-in LU
        printf("level %3d, block size %7d\n",i,next->LU.nr);
	k=nnzL;
	l=nnzU;
	if (i<nlev || next->LU.ja!=NULL) {
	    for (j=0; j<next->LU.nr-1; j++) { 
	        nnzL+=next->LU.ia[j]  -next->LU.ja[j];
		nnzU+=next->LU.ja[j+1]-next->LU.ia[j];
	    }
	}
	if (i==nlev) {
	  if (next->LU.ja==NULL) {
	     printf("switched to full matrix processing\n");fflush(STDOUT);
	     tmp=-1;
	     j=next->LU.nr;
	     nnzL+=(j*j-j)/2;
	     nnzU+=(j*j-j)/2;
	  }
	  else {
	    /* ILUTP case */
	    if (next->LUperm!=NULL) {
#ifdef USE_MPILUC
	      printf("  MPILUC failed, switched to ILUTP\n");
#else
	      printf("   PILUC failed, switched to ILUTP\n");
#endif
	      nnzL+=next->LU.ia[j]  -next->LU.ja[j];
	    }
	  } /* end if-else (i==nlev) */
	}
	printf("  local fill-in L %7d(%5.1lfav),   U %7d(%5.1lfav),   sum %7d(%5.1lfav)\n",
	       nnzL-k,(1.0*(nnzL-k))/next->LU.nr,
	       nnzU-l+next->LU.nr,(1.0*(nnzU-l+next->LU.nr))/next->LU.nr,
	       nnzL-k+nnzU-l+next->LU.nr,((1.0)*(nnzL-l+nnzU-k+next->LU.nr))/next->LU.nr);

	if (i<nlev) {
	   // fill-in E
	   nnzL+=next->E.ia[next->E.nr]-1;
	   // 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->E.nr);
	   printf("  local fill-in E %7d(%5.1lfav pc),F %7d(%5.1lfav pr),sum %7d(%5.1lfav)\n",
		  next->E.ia[next->E.nr]-1,(1.0*(next->E.ia[next->E.nr]-1))/next->LU.nr,
		  next->F.ia[next->F.nr]-1,(1.0*(next->F.ia[next->F.nr]-1))/next->LU.nr,
	          next->E.ia[next->E.nr]-1+next->F.ia[next->F.nr]-1,
		  ((1.0)*(next->E.ia[next->E.nr]-1+next->F.ia[next->F.nr]-1))/next->LU.nr);
	}
	next=next->next;
    }
    printf("\ntotal fill-in L&E%8d(%5.1lfav), U&F%8d(%5.1lfav),   sum%8d(%5.1lfav)\n",
	   nnzL,(1.0*nnzL)/n,nnzU+n,(1.0*(nnzU+n))/n,nnzL+nnzU+n,(1.0*(nnzL+nnzU+n))/n);
    printf("fill-in factor:      %5.1lf\n",(1.0*(nnzL+nnzU+n))/nz);
    
    if (tmp) {
       // nnzL+nnzU-j*j+n-j memory for sparse data structures
       //                   indices (weight 1/3) and values (weight 2/3)
       // j*j               memory for dense data, no indices (weight 2/3)
       printf("memory usage factor: %5.1lf\n",(1.0*(nnzL+nnzU-j*j+n-j))/nz
	                                     +(2.0*(j*j))/(3.0*nz));
       fprintf(fo,"%5.1f|",(1.0*(nnzL+nnzU-j*j+n-j))/nz
	                  +(2.0*(j*j))/(3.0*nz));
    }
    else {
       printf("memory usage factor: %5.1lf\n",(1.0*(nnzL+nnzU+n))/nz);
       fprintf(fo,"%5.1f|",(1.0*(nnzL+nnzU+n))/nz);
    }

    printf("total time: %8.1le [sec]\n",  (double)secnds); 
    printf("            %8.1le [sec]\n\n",(double)ILUPACK_secnds[7]); 
    fprintf(fo,"%7.1le|",(double)secnds);

#ifdef USE_MPILUC
    printf("refined timings for MPILUC\n"); 
#else
    printf("refined timings for  PILUC\n"); 
#endif
    printf("initial preprocessing:       %8.1le [sec]\n",ILUPACK_secnds[0]); 
    printf("reorderings remaining levels:%8.1le [sec]\n",ILUPACK_secnds[1]); 
#ifdef USE_MPILUC
    printf("MPILUC (sum over all levels):%8.1le [sec]\n",ILUPACK_secnds[2]); 
#else
    printf("PILUC (sum over all levels): %8.1le [sec]\n",ILUPACK_secnds[2]); 
#endif
    printf("ILUTP2 (if used):            %8.1le [sec]\n",ILUPACK_secnds[3]); 
    printf("LUPQ (if used):              %8.1le [sec]\n",ILUPACK_secnds[4]); 
    printf("remaining parts:             %8.1le [sec]\n\n",MAX(0.0,(double)secnds
	                                                   -ILUPACK_secnds[0]
							   -ILUPACK_secnds[1]
							   -ILUPACK_secnds[2]
							   -ILUPACK_secnds[3]
							   -ILUPACK_secnds[4]));

    fflush(STDOUT);

    /*
    next=&PRE;
    for (i=1; i<=nlev; i++) {
      printf("%d,%d\n",next->n,next->nB);
      for (j=0; j<next->n; j++)
	printf("%8d",next->p[j]);
      printf("\n");
      for (j=0; j<next->n; j++)
	printf("%8d",next->invq[j]);
      printf("\n");
        next=next->next;
    }
    fflush(STDOUT);

    printf("multilevel structure\n");
    fflush(STDOUT);
    printf("number of levels: %d\n",nlev);
    fflush(STDOUT);
    next=&PRE;
    for (i=1; i<=nlev; i++) {
      printf("total size %d\n",next->n);
      fflush(STDOUT);
      printf("leading block %d\n",next->nB);
      fflush(STDOUT);

      printf("row permutation\n");
      for (j=0; j<next->n; j++)
	printf("%4d",next->p[j]);
      printf("\n");
      fflush(STDOUT);
      printf("inverse column permutation\n");
      for (j=0; j<next->n; j++)
	printf("%4d",next->invq[j]);
      printf("\n");
      fflush(STDOUT);

      if (nlev==1) {
	printf("row scaling\n");
	for (j=0; j<next->n; j++)
	  printf("%12.4e",next->rowscal[j]);
	printf("\n");
	fflush(STDOUT);
	printf("column scaling\n");
	for (j=0; j<next->n; j++)
	  printf("%12.4e",next->colscal[j]);
	printf("\n");
	fflush(STDOUT);
      }

      printf("(2,1) block E (%d,%d)\n",next->E.nr,next->E.nc);
      fflush(STDOUT);
      for (k=0; k<next->E.nr; k++) {
	printf("%3d: ",k+1);
	for (j=next->E.ia[k]-1; j<next->E.ia[k+1]-1;j++)
	  fprintf(STDOUT,"%12d",next->E.ja[j]);
	printf("\n");
	fflush(STDOUT);
	printf("     ");
	for (j=next->E.ia[k]-1; j<next->E.ia[k+1]-1;j++)
	  fprintf(STDOUT,"%12.4e",next->E.a[j]);
	printf("\n");
	fflush(STDOUT);
      }
      printf("(1,2) block F (%d,%d)\n",next->F.nr,next->F.nc);
      fflush(STDOUT);
      for (k=0; k<next->F.nr; k++) {
	printf("%3d: ",k+1);
	for (j=next->F.ia[k]-1; j<next->F.ia[k+1]-1;j++)
	  fprintf(STDOUT,"%12d",next->F.ja[j]);
	printf("\n");
	fflush(STDOUT);
	printf("     ");
	for (j=next->F.ia[k]-1; j<next->F.ia[k+1]-1;j++)
	  fprintf(STDOUT,"%12.4e",next->F.a[j]);
	printf("\n");
	fflush(STDOUT);
      }

      printf("ILU...\n");
      printf("Diagonal part\n");
      for (k=0; k<next->LU.nr; k++) {
	fprintf(STDOUT,"%12.4e",next->LU.a[k]);
      }
      printf("\n");
      printf("L part\n");
      for (k=0; k<next->LU.nr; k++) {
	printf("col %3d: ",k+1);
	for (j=next->LU.ja[k]-1; j<next->LU.ia[k]-1;j++)
	  fprintf(STDOUT,"%12d",next->LU.ja[j]);
	printf("\n");
	fflush(STDOUT);
	printf("         ");
	for (j=next->LU.ja[k]-1; j<next->LU.ia[k]-1;j++)
	  fprintf(STDOUT,"%12.4e",next->LU.a[j]);
	printf("\n");
	fflush(STDOUT);
      }
      printf("U part\n");
      for (k=0; k<next->LU.nr; k++) {
	printf("row %3d: ",k+1);
	for (j=next->LU.ia[k]-1; j<next->LU.ja[k+1]-1;j++)
	  fprintf(STDOUT,"%12d",next->LU.ja[j]);
	printf("\n");
	fflush(STDOUT);
	printf("         ");
	for (j=next->LU.ia[k]-1; j<next->LU.ja[k+1]-1;j++)
	  fprintf(STDOUT,"%12.4e",next->LU.a[j]);
	printf("\n");
	fflush(STDOUT);
      }

      next=next->next;
    }
    */
    

    // remember that scaling is explicitly applied to A
    leftscale =PRE.rowscal;
    rightscale=PRE.colscal;


    if (rhstyp[2]=='X' || rhstyp[2]=='x') {
       j=nrhs*n;
       if (rhstyp[1]=='G' || rhstyp[1]=='g') 
	 j*=2;
       if (rhstyp[0]=='M' || rhstyp[0]=='m')
	 j-=nrhs*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]=1.0;
#else
	   sol[i].r=1.0;
	   sol[i].i=0.0;
#endif
      } // end for i



    /* rescale solution */
    for (i=0; i<n; i++) {
#if defined _DOUBLE_REAL_ || defined _SINGLE_REAL_
        sol[i]=1.0/rightscale[i];
#else
	val=1.0/(rightscale[i].r*rightscale[i].r
                +rightscale[i].i*rightscale[i].i);
	vb=sol[i].r;
	sol[i].r=( vb*rightscale[i].r+sol[i].i*rightscale[i].i)*val;
	sol[i].i=(-vb*rightscale[i].i+sol[i].i*rightscale[i].r)*val;
        //sol[i].r= rightscale[i].r*val;
        //sol[i].i=-rightscale[i].i*val;
#endif
    } // end for i
    /*
    for (i=0; i<n; i++)
      printf("%12.4e",sol[i]);
    printf("\n");
    fflush(STDOUT);
    for (i=0; i<n; i++) {
      printf("%3d:",i+1);
      for (j=A.ia[i]-1; j<A.ia[i+1]-1; j++)
	printf("%12d", A.ja[j]);
      printf("\n");
      fflush(STDOUT);
      printf("    ");
      for (j=A.ia[i]-1; j<A.ia[i+1]-1; j++)
	printf("%12.4e", A.a[j]);
      printf("\n");
      fflush(STDOUT);
    }
    */

    // release part of rhs that may store the uncompressed rhs
    if (rhstyp[0]=='M' || rhstyp[0]=='m')
       rhs-=n;

    if (nrhs==0) {
      /* right hand side rhs=A*sol*/
      CSRMATVEC(A,sol,rhs);
    }
    else {
       if (rhstyp[0]=='M' || rhstyp[0]=='m') {
	  for (i=0; i<n; i++) {
#if defined _DOUBLE_REAL_ || defined _SINGLE_REAL_ 
	      rhs[i]=0;
#else
	      rhs[i].r=rhs[i].i=0;