printperformance.c

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

#ifdef USE_MPILUC
    printf("M-ILUPACK, multilevel structure\n");
#else
    printf("ILUPACK,   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) {
 	   if (flags&COARSE_REDUCE) {
	      // 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);
	   } // end if
	   else {
	      printf("level %3d->%3d, block size (%7d,%7d)\n",i,i+1,next->LU.nr,next->E.nr);
	      fflush(stdout);
	   } // end else
	}
	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,(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 M-ILUPACK multilevel factorization\n"); 
#else
    printf("refined timings for   ILUPACK multilevel factorization\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("ILUTP (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;
    }
    */