mainildlc.c

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

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#include <blas.h>
#include <sparspak.h>
#include <ilupack.h>

#include <ilupackmacros.h>

#define MAX_LINE        255
#define STDERR          stdout
#define STDOUT          stdout
#define PRINT_INFO
#define MAX(A,B)        (((A)>(B))?(A):(B))
#define MIN(A,B)        (((A)<(B))?(A):(B))

// maximum number of iterations independent on n
#define MAX_IT          500

// measure for terminating iterative solver
#define RESTOL_FUNC(A)  sqrt(A)
//#define RESTOL_FUNC(A)  (A)

// use SPD scaling
#define USE_SCALE


// use an iterative solver from SPARSKIT defined by variable SOLVER
#define USE_SPARSKIT







int main(int argc, char **argv)
{
    /* SOLVER choice:   1  pcg
                        2  cgnr
                        3  bcg
                        4  dbcg
                        5  bcgstab
                        6  tfqmr
                        7  fom
                        8  gmres
                        9  fgmres
                       10  dqgmres */
    int SOLVER=1; /* gmres */

    CSRMAT A, ilutmat;

    FILE *fp, *fo; 
    char rhstyp[3], title[72], key[8], type[3], fname[100], foname[100];
    char line[MAX_LINE], *tmpstring, timeout[7], residout[7];
    int  i,j,k,m,fnamelen,n,nc,nz,nrhs,tmp0,tmp,tmp2,tmp3,ierr,ipar[16],flag,
         *invperm, *buff, *ws, *symmmd,
         nrhsix, *rhsptr, *rhsind;
    REALS eps,fpar[16], DROP_TOL;
    FLOAT *rhs,*sol,*w, *dbuff, *scale, *rhsval;
    float  systime, time_start,   time_stop,   secnds, secndsprep;
    int ELBOW, param;
    size_t nibuff, ndbuff;

    /* the last argument passed serves as file name */
    if (argc!=4) {
      printf("usage '%s <drop tol.> <elbow space> <matrix>'\n",argv[0]);
       exit(0);
    }
    i=0;
    while (argv[argc-1][i]!='\0')
    {
          fname[i]=argv[argc-1][i];
          i++;
    } /* end while */
    fname[i]='\0';
    fnamelen=i;
    while (i>=0 && fname[i]!='/')
          i--;
    i++;
    j=0;
    while (i<fnamelen && fname[i]!='.')
          foname[j++]=fname[i++];
    while (j<16)
          foname[j++]=' ';
    foname[j]='\0';

    ELBOW=atoi(argv[argc-2]);
    DROP_TOL=atof(argv[argc-3]);


/*----------------------------------------------------------------------
|     Read a Harwell-Boeing matrix.
|     Use readmtc first time to determine sizes of arrays.
|     Read in values on the second call.
|---------------------------------------------------------------------*/

    nrhs = 0;
    tmp0 = 0;

    if ((fp=fopen(fname,"r"))==NULL) {
        fprintf(STDERR," file %s not found\n",fname);
        exit(0);
    }
    fclose(fp);

    READMTC(&tmp0,&tmp0,&tmp0,fname,A.a,A.ja,A.ia,
	    rhs,&nrhs,rhstyp,&n,&nc,&nz,title,key,type,
	    &nrhsix,rhsptr,rhsind,rhsval,&ierr,fnamelen,2,72,8,3);
    if (ierr) {
        fprintf(STDERR," ierr = %d\n",ierr);
        fprintf(STDERR," error in reading the matrix, stop.\n");
	switch(ierr) {
	case 1:
	  fprintf(STDERR,"too many columns\n");
	  break;  
	case 2:
	  fprintf(STDERR,"too many nonzeros\n");
	  break;  
	case 3:
	  fprintf(STDERR,"too many columns and nonzeros\n");
	  break;  
	case 4:
	  fprintf(STDERR,"right hand side has incompatible type\n");
	  break;  
	case 5:
	  fprintf(STDERR,"too many right hand side entries\n");
	  break;  
	case 6:
	  fprintf(STDERR,"wrong type (real/complex)\n");
	  break;  
	}
        exit(ierr);
    }
    printf("Matrix: %s: size (%d,%d), nnz=%d(%4.1lfav.)\n", fname, n,nc,
	   nz,((double)nz)/n);

    if (fname[fnamelen-1]=='5')
      fo = fopen("out_mc64","aw");
    else
      fo = fopen("out_normal","aw");

    fprintf(fo,"%s|%7.1e|       |",foname,DROP_TOL);

    m=1;
    if (nrhs>0) {
       printf ("Number of right hand sides supplied: %d \n", nrhs) ;
       if (rhstyp[1]=='G' || rhstyp[1]=='g') {
	 m++;
	 printf ("Initial solution(s) offered\n") ;
       }
       else
	 printf ("\n") ;

       if (rhstyp[2]=='X' || rhstyp[2]=='x') {
	 m++;
	 printf ("Exact solution(s) provided\n") ;
       }
       else
	 printf ("\n") ;
    }
    else 
      printf("\n\n\n");
    printf("\n");

    rhsptr=NULL;
    rhsind=NULL;
    rhsval=NULL;
    if (rhstyp[0]=='M' || rhstyp[0]=='m') {
       rhsptr=(int *)  MALLOC((size_t)(nrhs+1)*sizeof(int),"main:rhsptr");
       rhsind=(int *)  MALLOC((size_t)nrhsix*sizeof(int),  "main:rhsind");
       rhsval=(FLOAT *)MALLOC((size_t)nrhsix*sizeof(FLOAT),"main:rhsval");
       // no dense right hand side
       m--;
       m*=n*MAX(1,nrhs);
       // in any case we need at least one vector for the r.h.s.
       m+=n;
    }
    else
       m*=n*MAX(1,nrhs);
    A.ia=(int *)  MALLOC((size_t)(n+1)*sizeof(int),"main:A.ia");
    A.ja=(int *)  MALLOC((size_t)nz*sizeof(int),   "main:A.ja");
    A.a =(FLOAT *)MALLOC((size_t)nz*sizeof(FLOAT), "main:A.a");
    A.nr=n;
    A.nc=n;
    rhs       =(FLOAT *) MALLOC((size_t)m*sizeof(FLOAT),  "main:rhs");
    // advance pointer to reserve space when uncompressing the right hand side
    if (rhstyp[0]=='M' || rhstyp[0]=='m')
       rhs+=n;
    
    sol       =(FLOAT *) MALLOC((size_t)n*sizeof(FLOAT),  "main:sol");
    dbuff     =(FLOAT *) MALLOC(3*(size_t)n*sizeof(FLOAT),"main:dbuff");
    ndbuff=3*(size_t)n;

    tmp = 3;
    tmp2 = n;
    tmp3 = nz;

    if (rhstyp[0]=='M' || rhstyp[0]=='m')
       m-=n;
    READMTC(&tmp2,&tmp3,&tmp,fname,A.a,A.ja,A.ia,
	    rhs,&m,rhstyp,&n,&nc,&nz,title,key,type,
	    &nrhsix,rhsptr,rhsind,rhsval,&ierr,fnamelen,2,72,8,3);
    if (rhstyp[0]=='M' || rhstyp[0]=='m')
       m+=n;



    if (ierr) {
        fprintf(STDERR," ierr = %d\n",ierr);
        fprintf(STDERR," error in reading the matrix, stop.\n");
	fprintf(fo,"out of memory\n");fclose(fo); 
        exit(ierr);
    }
    /*
    for (i=0; i<n;i++) {
      printf("%4d:\n",i+1);
      for (j=A.ia[i];j<A.ia[i+1]; j++)
	printf("%16d",A.ja[j-1]);
      printf("\n");
      fflush(stdout);
      for (j=A.ia[i];j<A.ia[i+1]; j++) 
#if defined _DOUBLE_REAL_ || defined _SINGLE_REAL_
	printf("%16.1e",A.a[j-1]);
#else
	printf("%8.1e%8.1e",A.a[j-1].r,A.a[j-1].i);
#endif
      printf("\n");
      fflush(stdout);
    }// end for i
    */

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

    // if no right hand side is provided, then set sol=1 and rhs=A*sol
    if (nrhs==0) { 
       for (i=0; i<n; i++) {
#if defined _DOUBLE_REAL_ || defined _SINGLE_REAL_
	   sol[i]=1;
#else
	   sol[i].r=1;
	   sol[i].i=0;
#endif
       } // end for i
       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
    } // end if-else
    
    
    evaluate_time(&time_start,&systime);
    scale=(FLOAT *)MALLOC((size_t)n*sizeof(FLOAT),"mainildlc:scale");
    for (i=0; i<n; i++) {
#if defined _DOUBLE_REAL_ || defined _SINGLE_REAL_
        scale[i]=1.0;
#else
        scale[i].r=1.0;
        scale[i].i=0.0;
#endif
    } // end for i

#ifdef USE_SCALE
    i=0;
    SPDSCALE(&(A.nc),A.a, A.ja, A.ia, scale, &i);
    
    if(i!=0) {
      fprintf(STDERR, "SPDSCALE: a zero row/column...\n" );
      return;
    }
#endif /* USE_SCALE */

    // rescale right hand side
    for (i=0; i<n; i++) {
#if defined _DOUBLE_REAL_ || defined _SINGLE_REAL_
        rhs[i]*=scale[i];
#else
	rhs[i].r*=scale[i].r;
	rhs[i].i*=scale[i].r;
#endif
    } // end for i


    /*
    for (i=0; i<n; i++) {
        printf("%4d:\n",i+1);
        for (j=A.ia[i]; j<A.ia[i+1]; j++) {
	  printf("%8d",A.ja[j-1]);
	} // end for j
	printf("\n");fflush(stdout);
        for (j=A.ia[i]; j<A.ia[i+1]; j++) {
#if defined _DOUBLE_REAL_ || defined _SINGLE_REAL_
	    printf("%16.1e",A.a[j-1]);
#else
	    printf("%8.1e%8.1e",A.a[j-1].r,A.a[j-1].i);
#endif
	} // end for j
	printf("\n");fflush(stdout);
    } // end for i
    */


    fprintf(fo,"   s/   s|");
    evaluate_time(&time_stop,&systime);
    secndsprep=time_stop-time_start;
    printf("time: %8.1e [sec]\n",(double)secndsprep);



    /* ------------------------------------------------------------------- */
    /* -----------------------   Start MAIN part   ----------------------- */
    /* ------------------------------------------------------------------- */

    fprintf(fo,"IB-ILDLC  |");










    /* set up paramaters for ILUT */
    /* n       = integer. The row dimension of the matrix A. The matrix      */
    /* a,ja,ia = matrix stored in Compressed Sparse Row format.              */
    /* lfil    = integer. The fill-in parameter. Each row of L and each row
                 of U will have a maximum of lfil elements (excluding the 
		 diagonal element). lfil must be .ge. 0.
		 ** WARNING: THE MEANING OF LFIL HAS CHANGED WITH RESPECT TO
		 EARLIER VERSIONS.                                           */
    ipar[0]=n+1;
    /* droptol = real*8. Sets the threshold for dropping small terms in the
                 factorization. See below for details on dropping strategy.  */
    fpar[0]=DROP_TOL;
    /* iwk     = integer. The lengths of arrays alu and jlu. If the arrays
                 are not big enough to store the ILU factorizations, ilut
		 will stop with an error message.                            */
    ipar[1]=ELBOW*nz;
    /* On return:
       ===========

       alu,jlu = matrix stored in Modified Sparse Row (MSR) format containing
                 the L and U factors together. The diagonal (stored in
		 alu(1:n) ) is inverted. Each i-th row of the alu,jlu matrix
		 contains the i-th row of L (excluding the diagonal entry=1)
		 followed by the i-th row of U.                              */
    ilutmat.ja=(int *)  MALLOC((size_t)ipar[1]*sizeof(int),  "mainildlc:ilutmat.ja");
    ilutmat.a =(FLOAT *)MALLOC((size_t)ipar[1]*sizeof(FLOAT),"mainildlc:ilutmat.a");
    /* ju      = integer array of length n containing the pointers to
                 the beginning of each row of U in the matrix alu,jlu.       */
    ilutmat.ia=(int *)  MALLOC((size_t)n*sizeof(int),"mainildlc:ilutmat.ia");   
    /* ierr    = integer. Error message with the following meaning.
                 ierr  = 0    --> successful return.
		 ierr .gt. 0  --> zero pivot encountered at step number ierr.
		 ierr  = -1   --> Error. input matrix may be wrong.
                                  (The elimination process has generated a
				  row in L or U whose length is .gt.  n.)
                 ierr  = -2   --> The matrix L overflows the array alu.
		 ierr  = -3   --> The matrix U overflows the array alu.
		 ierr  = -4   --> Illegal value for lfil.
		 ierr  = -5   --> zero row encountered.                      */
    ipar[2]=0;
    /* work arrays:
       =============
       jw      = integer work array of length 2*n.
       w       = real work array of length n                                 */
    w =(FLOAT *)MALLOC(3*(size_t)n*sizeof(FLOAT),"mainildlc:w");
    ws=(int *)  MALLOC(7*(size_t)n*sizeof(int),"mainildlc:ws");

    // perform ILUT decomposition
    // DROP_INVERSE|TISMENETSKY_SC|NO_SHIFT|REPEAT_FACT|IMPROVED_ESTIMATE
    // param=DROP_INVERSE|TISMENETSKY_SC|IMPROVED_ESTIMATE|DIAGONAL_COMPENSATION;
    param=DROP_INVERSE|IMPROVED_ESTIMATE|ENSURE_SPD;

    printf("ILDLC PARAMETERS:\n");
    printf("   droptol=%8.1e\n",DROP_TOL);
    printf("   lfil   =%d\n",  ipar[0]);
    printf("   elbow space factor=%4d\n",  ELBOW);
    if (param & DROP_INVERSE)
       printf("   inverse-based dropping\n");