mainspdexpert.c

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

    SPDAMGSETPARAMS(A, &param, flags, elbow, droptols, condest,
		    restol, max_it);


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


    evaluate_time(&time_start,&systime);
    ierr=SPDAMGSETUP(&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);
    }
#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=⪯
    nnzU=0;
    tmp=0;

    for (i=1; i<=nlev; i++) {
        // fill-in LU
        printf("level %3d, block size %7d\n",i,next->LU.nr); fflush(stdout);
	l=nnzU;
	if (i<nlev || next->LU.ja!=NULL) {
	   nnzU+=next->LU.ja[next->LU.nr-1]-next->LU.ja[0];
	}
	if (i==nlev) {
	   if (next->LU.ja==NULL) {
	      printf("switched to full matrix processing\n");fflush(STDOUT);
	      tmp=-1;
	      j=next->LU.nr;
	      nnzU+=(j*(j-1))/2;
	   }
	}
	printf("  local fill-in %7d(%5.1lfav)\n",
	       nnzU-l+next->LU.nr,(1.0*(nnzU-l+next->LU.nr))/next->LU.nr);

	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;
    }
    printf("\ntotal fill-in sum%8d(%5.1lfav)\n",
	   nnzU+n,(1.0*(nnzU+n))/n);
    printf("fill-in factor:      %5.1lf\n",(1.0*(nnzU+n))/nz);
    
    if (tmp) {
       // nnzU-j*(j+1)/2+n-j memory for sparse data structures
       //                    indices (weight 1/3) and values (weight 2/3)
       // j*(j+1)/2          memory for dense data, no indices (weight 2/3)
       printf("memory usage factor: %5.1lf\n",(1.0*(nnzU-(j*(j+1))/2+n-j))/nz
	                                     +(j*(j+1))/(3.0*nz));
       fprintf(fo,"%5.1f|",(1.0*(nnzU-(j*(j+1))/2+n-j))/nz
                          +(j*(j+1))/(3.0*nz));
    }
    else {
       printf("memory usage factor: %5.1lf\n",(1.0*(nnzU+n))/nz);
       fprintf(fo,"%5.1f|",(1.0*(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 MPILDLC\n"); 
#else
    printf("refined timings for  PILDLC\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("MPILDLC(sum over all levels):%8.1le [sec]\n",ILUPACK_secnds[2]); 
#else
    printf("PILDLC (sum over all levels):%8.1le [sec]\n",ILUPACK_secnds[2]); 
#endif
    printf("ILDLC (if used):             %8.1le [sec]\n",ILUPACK_secnds[3]); 
    printf("LDLP (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;
    }
    */
    

    // 3n words are required by AMGsol
    param.ndbuff=MAX(param.ndbuff,3*(size_t)n);
    ndbuff=param.ndbuff;
    param.dbuff=(FLOAT *)REALLOC(param.dbuff,param.ndbuff*sizeof(FLOAT),
				 "main:ndbuff");
    

    // remember that scaling is explicitly applied to A
    scale =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/scale[i];
#else
	val=1.0/(scale[i].r*scale[i].r
                +scale[i].i*scale[i].i);
	vb=sol[i].r;
	sol[i].r=( vb*scale[i].r+sol[i].i*scale[i].i)*val;
	sol[i].i=(-vb*scale[i].i+sol[i].i*scale[i].r)*val;
        //sol[i].r= scale[i].r*val;
        //sol[i].i=-scale[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
       HERMATVEC(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;
#endif
	  }// end for i

	  // uncompress rhs
	  for (i=0; i<rhsptr[1]-rhsptr[0]; i++) {
	      j=rhsind[i]-1;
	      rhs[j]=rhsval[i];
	  } // end for i
       } // end if
       for (i=0; i<n; i++) {
#if defined _DOUBLE_REAL_ || defined _SINGLE_REAL_
	   rhs[i]*=scale[i];
#else
	   val=rhs[i].r;
	   rhs[i].r=val*scale[i].r-rhs[i].i*scale[i].i;
	   rhs[i].i=val*scale[i].i+rhs[i].i*scale[i].r;
#endif
       } // end for i
    } // end else


    //printf("matrix-vector multiply done\n");fflush(STDOUT);
    /*
    printf("true right hand side\n");
    for (i=0; i<n; i++)
      printf("%12.4e",rhs[i]);
    printf("\n");
    fflush(STDOUT);
    printf("start GMRES\n");fflush(STDOUT);
    */

    /* --------------------   apply preconditioner   --------------------- */
    if (ierr==0)
    {
#ifdef USE_SPARSKIT
       /* initial solution */
       if (rhstyp[1]=='G' || rhstyp[1]=='g') {
	  j=nrhs*n;
	  if (rhstyp[0]=='M' || rhstyp[0]=='m')
	     j=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_