symamgsetup.c

a software code for computing selected eigenvalues of large sparse symmetric matrices

#ifdef PRINT_INFO2
	printf("lev %d, colscal allocated\n",*nlev);fflush(stdout);
#endif

	for (i=0; i<n; i++) {
#if defined _DOUBLE_REAL_ || defined _SINGLE_REAL_
	    current->colscal[i]=1;
#else
	    current->colscal[i].r=1;
	    current->colscal[i].i=0;
#endif
	}

	current->p   =(integer *)MALLOC((size_t)n*sizeof(integer),"symAMGsetup:current->p");
	current->invq=(integer *)MALLOC((size_t)n*sizeof(integer),"symAMGsetup:current->invq");
	ILUPACK_mem[8]+=2*n;
	for (i=0; i<n; i++)
	  current->p[i]=current->invq[i]=i+1;
	current->nB=current->n;
#ifdef PRINT_INFO2
	printf("lev %d, p,invq allocated\n",*nlev);fflush(stdout);
#endif


	// auxilliary memory that is currently left over
	IPparam->niaux=IPparam->ndaux=mem;
	IPparam->iaux=jlu;
	IPparam->daux=alu;


	/* check whether the matrix is almost dense */
	/* We define a matrix to be dense if nnz(S)>1/3 n(n+1)/2 */
	nnz=S.ia[n]-1;

#ifdef PRINT_INLINE
	printf("n=%7d, nnz=%9d, %6.1fav",n,nnz,((double)nnz)/n);
#endif

	if (*nlev>1 && ilucancel*nnz>((REALS)n)*(n+1)/2.0) {

#ifdef PRINT_INLINE
	   printf("\nswitched to full matrix processing\n");fflush(STDOUT);
#endif

	   // start counter for SPTRF
	   evaluate_time(&time_start,&systime);
	   /* uncompress S if possible */
	   if (mem<(n*(size_t)(n+1))/2-nnz) {
	      evaluate_time(&time_stop,&systime);
	      // compute total setup time
	      ILUP_sec[7]=time_stop-ILUP_sec[7];
	      // compute TOTAL time
	      ILUP_sec[8]=time_stop-ILUP_sec[8];
	      for (i=0; i<ILUPACK_secnds_length; i++)
		  ILUPACK_secnds[i]=ILUP_sec[i];
	      ierr=-2;
	      current->E.ia=current->F.ia=NULL;
	      current->E.ja=current->F.ja=NULL;
	      current->E.a =current->F.a =NULL;
	      current->LU.ja=jlu;
	      current->absdiag=NULL;

	      // undo scaling
	      r=PRE->rowscal;
	      c=PRE->colscal;
	      // adjust arrays
	      ia=A->ia;
	      ja=A->ja;
	      a =A->a;
	      ja--;
	      a--;
	      c--;
	      for (i=0; i<A->nr; i++) {
		for (j=ia[i]; j<ia[i+1]; j++) {
		    k=ja[j];
#if defined _SINGLE_REAL_ || defined _DOUBLE_REAL_
		    a[j]/=r[i]*c[k];
#else
		    a[j].r/=r[i].r*c[k].r;
		    a[j].i/=r[i].r*c[k].r;
#endif	  
		}
	      }

	      // MYSYMAMGDELETE(*A, *PRE, *nlev, IPparam);
	      return (ierr);
	   }
	   
	   /* position in alu behind the uncompressed matrix */
	   alu+=(n*(size_t)(n+1))/2-nnz;
	   delta[*nlev-1]-=nnz;
	   // additional memory needed
	   ILUPACK_mem[2]+=0;
	   ILUPACK_mem[3]+=(n*(size_t)(n+1))/2-nnz;
	   ILUPACK_mem[5]=MAX(ILUPACK_mem[5],ILUPACK_mem[3]);

	   S.ja--;
	   S.a--;
	   for (i=n-1; i>=0; i--) {
	       for (j=i; j<n; j++) 
#if defined _DOUBLE_REAL_ || defined _SINGLE_REAL_
		   IPparam->dbuff[j]=0;
#else
	           IPparam->dbuff[j].r=IPparam->dbuff[j].i=0;
#endif
	       IPparam->dbuff--;
	       for (j=S.ia[i]; j<S.ia[i+1]; j++)
		   IPparam->dbuff[S.ja[j]]=S.a[j];
	       IPparam->dbuff++;
	       alu-=(n-i);
	       for (j=i; j<n; j++)
		   alu[j-i]=IPparam->dbuff[j];
	   }
	   S.ja++;
	   S.a++;
	   current->LUperm=NULL;
	   
#ifdef PRINT_INFO
	   printf("dense A\n");fflush(STDOUT);
	   for (i=0; i<n; i++)
	     printf("%12d",S.ia[i]);
	   printf("\n");fflush(stdout);
	   j=0;
	   for (k=0; k<n; k++) {
	     for (i=k; i<n; i++) {
	       printf("%12.4le",S.a[j++]);
	     }
	     printf("\n");fflush(stdout);
	   }
	   printf("\n");fflush(stdout);
#endif
	   MYSPTRF("L",&n,S.a,S.ia,&ierr,1); 
#ifdef PRINT_INFO
	   printf("dense U\n");fflush(STDOUT);
	   for (i=0; i<n; i++)
	     printf("%12d",S.ia[i]);
	   printf("\n");fflush(stdout);
	   j=0;
	   for (k=0; k<n; k++) {
	     for (i=k; i<n; i++) {
	       printf("%12.4le",S.a[j++]);
	     }
	     printf("\n");fflush(stdout);
	   }
	   printf("\n");fflush(stdout);
#endif
	   // store SPTRF computation time
	   evaluate_time(&time_stop,&systime);
	   ILUP_sec[4]=time_stop-time_start;

	   if (ierr) { 
	      // compute total setup time
	      ILUP_sec[7]=time_stop-ILUP_sec[7];
	      // compute TOTAL time
	      ILUP_sec[8]=time_stop-ILUP_sec[8];
	      for (i=0; i<ILUPACK_secnds_length; i++)
	          ILUPACK_secnds[i]=ILUP_sec[i];
	      current->E.ia=current->F.ia=NULL;
	      current->E.ja=current->F.ja=NULL;
	      current->E.a =current->F.a =NULL;
	      current->LU.ja=jlu;
	      current->absdiag=NULL;

	      // undo scaling
	      r=PRE->rowscal;
	      c=PRE->colscal;
	      // adjust arrays
	      ia=A->ia;
	      ja=A->ja;
	      a=A->a;
	      ja--;
	      a--;
	      c--;
	      for (i=0; i<A->nr; i++) {
		for (j=ia[i]; j<ia[i+1]; j++) {
		    k=ja[j];
#if defined _SINGLE_REAL_ || defined _DOUBLE_REAL_
		    a[j]/=r[i]*c[k];
#else
		    a[j].r/=r[i].r*c[k].r;
		    a[j].i/=r[i].r*c[k].r;
#endif	  
		}
	      }

	      // MYSYMAMGDELETE(*A, *PRE, *nlev, IPparam);
	      return (ierr);
	   }
	   current->nB=n;
	     
	   current->E.nr=current->F.nc=0;
	   current->E.nc=current->F.nr=0;
	   
	   current->LU.nr=current->nB;
	   current->LU.nc=current->LU.nr;
	   current->LU.a =alu;
	   current->LU.ja=NULL;
	   current->LU.ia=S.ia;

	   n=-1;
	   S.nr=0;
	}
	else { /* the matrix is still sparse */

	  // apply reordering techniques 
	  
	  current->absdiag=NULL;
	  // initial reordering
	  if (*nlev==1 && (param&PREPROCESS_INITIAL_SYSTEM)) {
	     IPparam->ipar[7]&=~1024;
	     IPparam->ipar[7]|=512;
	     ierr=(*perm0)(S, current->rowscal, current->colscal,current->p,
			   current->invq, &current->nB, IPparam);
	     IPparam->ipar[7]&=~(512+1024);
	  }
	  else 
	    // final pivoting
	    if (*nlev>1 && (param&FINAL_PIVOTING) && !regular_pivoting) { 
	       IPparam->ipar[7]|=512+1024;
	       ierr=(*permf)(S, current->rowscal, current->colscal, current->p,
			     current->invq, &current->nB, IPparam);
	       IPparam->ipar[7]&=~(512+1024);
	    } // end if-elseif-if
	    else
	       // regular pivoting
	       if (*nlev>1 && (param&PREPROCESS_SUBSYSTEMS)) { 
		  IPparam->ipar[7]&=~512;
		  IPparam->ipar[7]|=1024;
		  ierr=(*perm)(S, current->rowscal, current->colscal, current->p,
			       current->invq, &current->nB, IPparam);
		  IPparam->ipar[7]&=~(512+1024);
	       } // end if-elseif-elseif

	  // mymaximem=MAX(ILUPACK_mem[12],ILUPACK_mem[13]);

	  evaluate_time(&time_stop,&systime);
	  if (*nlev==1)
	     // time for initial preprocessing + scaling
	     ILUP_sec[0] =time_stop-time_start;
	  else
	     // accumulated time for reordering/pivoting and scaling 
	     // the remaining systems
	     ILUP_sec[1]+=time_stop-time_start;

	  if (ierr) {
	     // compute total setup time
	     ILUP_sec[7]=time_stop-ILUP_sec[7];
	     // compute TOTAL time
	     ILUP_sec[8]=time_stop-ILUP_sec[8];
	     for (i=0; i<ILUPACK_secnds_length; i++)
		 ILUPACK_secnds[i]=ILUP_sec[i];
	     current->E.ia=current->F.ia=NULL;
	     current->E.ja=current->F.ja=NULL;
	     current->E.a =current->F.a =NULL;
	     current->LU.ja=jlu;
	     current->absdiag=NULL;

	     // undo scaling
	     r=PRE->rowscal;
	     c=PRE->colscal;
	     // adjust arrays
	     ia=A->ia;
	     ja=A->ja;
	     a=A->a;
	     ja--;
	     a--;
	     c--;
	     for (i=0; i<A->nr; i++) {
	         for (j=ia[i]; j<ia[i+1]; j++) {
		     k=ja[j];
#if defined _SINGLE_REAL_ || defined _DOUBLE_REAL_
		     a[j]/=r[i]*c[k];
#else
		     a[j].r/=r[i].r*c[k].r;
		     a[j].i/=r[i].r*c[k].r;
#endif	  
		 }
	     }

	     // MYSYMAMGDELETE(*A, *PRE, *nlev, IPparam);
	     return (ierr);
	  }

	  if ((n-current->nB > amgcancel*n || current->nB==0) && 
	      (param&FINAL_PIVOTING) && regular_pivoting) {
	     // switch to final pivoting
#ifdef PRINT_INLINE
  	     printf("\nswitched to final pivoting\n");fflush(STDOUT);
#endif
	     regular_pivoting=0;
	     amgcancel  =IPparam->fpar[5]; 
	     current->nB=n;
	     evaluate_time(&time_start,&systime);
	     IPparam->ipar[7]|=512+1024;
	     ierr=(*permf)(S, current->rowscal, current->colscal, current->p,
			   current->invq, &current->nB, IPparam);
	     IPparam->ipar[7]&=~(512+1024);
	     // mymaximem=MAX(mymaximem,ILUPACK_mem[12]);
	     // mymaximem=MAX(mymaximem,ILUPACK_mem[13]);

	     evaluate_time(&time_stop,&systime);
	     if (*nlev==1)
	        // time for initial preprocessing + scaling
	        ILUP_sec[0] =time_stop-time_start;
	     else
	        // accumulated time for reordering/pivoting and scaling 
	        // the remaining systems
	        ILUP_sec[1]+=time_stop-time_start;

	     if (ierr) {
	        // compute total setup time
	        ILUP_sec[7]=time_stop-ILUP_sec[7];
		// compute TOTAL time
		ILUP_sec[8]=time_stop-ILUP_sec[8];
		for (i=0; i<ILUPACK_secnds_length; i++)
		    ILUPACK_secnds[i]=ILUP_sec[i];
		current->E.ia=current->F.ia=NULL;
		current->E.ja=current->F.ja=NULL;
		current->E.a =current->F.a =NULL;
		current->LU.ja=jlu;
		current->absdiag=NULL;

		// undo scaling
		r=PRE->rowscal;
		c=PRE->colscal;
		// adjust arrays
		ia=A->ia;
		ja=A->ja;
		a=A->a;
		ja--;
		a--;
		c--;
		for (i=0; i<A->nr; i++) {
		    for (j=ia[i]; j<ia[i+1]; j++) {
		        k=ja[j];
#if defined _SINGLE_REAL_ || defined _DOUBLE_REAL_
			a[j]/=r[i]*c[k];
#else
			a[j].r/=r[i].r*c[k].r;
			a[j].i/=r[i].r*c[k].r;
#endif	  
		    }
		}

		// MYSYMAMGDELETE(*A, *PRE, *nlev, IPparam);
		return (ierr);
	     }
	  } // end if


	  nB=current->nB;
	  /* perform PILUC decomposition */

	  current->LUperm=NULL;

#ifdef PRINT_INLINE
	  printf(", nB=%7d->",current->nB);fflush(STDOUT);
#endif
	  
	  if (current->nB>0) {
	     // start counter for (m)piluc
	     evaluate_time(&time_start,&systime);
	     condest1=condest;
	     condest/=2.0;
	     mymaximem=0;
	     do {
	        current->nB=nB;
		ierr=0;
	        condest*=2.0;
		// printf("level %2d, condest=%8.1le\n",*nlev,condest);fflush(stdout);
	        imem=(LONG_INT)mem;
		if (*nlev>1)
		   PILUCparam|=discardA;
		shiftA=0;
		myimem=imem;
		if (param&MULTI_PILUC){
		   SYMPILUC(&n,S.a,S.ja,S.ia,lfil,droptols,&condest,
			    // IPparam->fpar+6,
			    &current->nB,&PILUCparam,current->p,current->invq,
			    alu,jlu,&myimem,IPparam->dbuff,IPparam->ibuff,
			    &shiftA,&amgcancel,&ierr);
		}
		else {
		   SYMPILUC(&n,S.a,S.ja,S.ia,lfil,droptols,&condest,
			    &current->nB,&PILUCparam,current->p,current->invq,
			    alu,jlu,&myimem,IPparam->dbuff,IPparam->ibuff,
			    &shiftA,&amgcancel,&ierr);
		}
		mymaximem=MAX(mymaximem,myimem);
	     }
	     while (nB-current->nB>amgcancel*nB && 
		    condest<1024*condest0 && condest<16*condest1);

	     alu-=shiftA;
	     jlu-=shiftA;
	     delta[*nlev-1]-=shiftA;
	     // stop counter for (m)piluc
	     evaluate_time(&time_stop,&systime);
	     // accumulate collective time for (m)piluc
	     ILUP_sec[2]+=time_stop-time_start;
	     if (ierr) {
	        // compute total setup time
	        ILUP_sec[7]=time_stop-ILUP_sec[7];
	        // compute TOTAL time
		ILUP_sec[8]=time_stop-ILUP_sec[8];
		for (i=0; i<ILUPACK_secnds_length; i++)
		    ILUPACK_secnds[i]=ILUP_sec[i];
		current->E.ia=current->F.ia=NULL;
		current->E.ja=current->F.ja=NULL;
		current->E.a =current->F.a =NULL;
		current->LU.ja=jlu;