symamgsetup.c

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

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

#ifdef PRINT_INLINE
	  printf("%7d\n",current->nB);fflush(STDOUT);
#endif


#ifdef PRINT_INFO
	  printf("1.computed U-factor\n");fflush(STDOUT);
	  for (i=0; i<current->nB; ) {
	    if (jlu[n+1+i]==0){
	      for (k=jlu[i];k<jlu[i+1]; k++) 
		printf("%8d",jlu[k-1]);
	      printf("\n");fflush(STDOUT);
	      for (k=jlu[i];k<jlu[i+1]; k++) 
		printf("%8.1le",alu[k-1]);
	      printf("\n");fflush(STDOUT);
	      i++;
	    }
	    else {
	      for (k=jlu[i];k<jlu[i+1]; k++) 
		printf("%8d",jlu[k-1]);
	      printf("\n");fflush(STDOUT);
	      for (k=jlu[i];k<jlu[i+1]; k++) 
		printf("%8.1le",alu[jlu[i]+2*(k-jlu[i])-1]);
	      printf("\n");fflush(STDOUT);
	      for (k=jlu[i];k<jlu[i+1]; k++) 
		printf("%8.1le",alu[jlu[i]+2*(k-jlu[i])]);
	      printf("\n");fflush(STDOUT);
	      i+=2;
	    }
	  }
	  printf("\n");fflush(STDOUT);
	  printf("Block diagonal factor\n");
	  for (k=0; k<current->nB;) {
	    if (jlu[n+1+k]==0){
	      printf("%8.1le",alu[k]);
	      k++;
	    }
	    else {
	      printf("%8.1le%8.1le",alu[k],alu[n+1+k]);
	      k+=2;
	    }
	  }
	  printf("\n");fflush(stdout);
	  for (k=0; k<current->nB; ) {
	    if (jlu[n+1+k]==0) {
	      printf("        ");
	      k++;
	    }
	    else {
	      printf("%8.1le%8.1le",alu[n+1+k],alu[k+1]);
	      k+=2;
	    }
	  }
	  printf("\n");fflush(stdout);
	  printf("computed Schur complement\n");fflush(STDOUT);
	  for (i=current->nB; i<n; i++) {
	    for (k=jlu[i];k<jlu[i+1]; k++) 
	      printf("%8d",jlu[k-1]);
	    printf("\n");fflush(STDOUT);
	    for (k=jlu[i];k<jlu[i+1]; k++) 
	      printf("%8.1le",alu[k-1]);
	    printf("\n");fflush(STDOUT);
	  }

	  for (i=0; i<=current->nB; i++)
	    printf("%8d",i+1);
	  printf("\n");fflush(stdout);
	  for (i=0; i<=current->nB; i++)
	    printf("%8d",jlu[i]);
	  printf("\n");fflush(stdout);
	  for (i=0; i<current->nB; i++)
	    printf("%8d",jlu[n+1+i]);
	  printf("\n");fflush(stdout);
#endif


	  
	  // COMMENT: Maybe in a later version the extraction of E and F and discarding the old
	  //          S should be done in one shot. This is reasonable, since (for *nlev>1) the old
	  //          matrix S is not needed any longer and to create some new space to store E and
	  //          F is not really necessary. But this is very technical, since the rows of S
	  //          have to be reordered in order to construct E and F :(


	  // appearently piluc failed
	  // BUT we set the final pivoting flag
	  // AND we haven't switched to final pivoting yet
	  if ((nB-current->nB > amgcancel*nB || 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);
	     }

	     nB=current->nB;
	     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);
		imem=(LONG_INT)mem;
		// if DISCARD_MATRIX is set, then at any level>1 the system
		// matrix is discarded as soon as possible to save memory
		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);
		printf("myimem=%8ld,condest=%8.1le\n",myimem,condest);
		mymaximem=MAX(mymaximem,myimem);

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

#ifdef PRINT_INLINE
	     printf("%7d\n",current->nB);fflush(STDOUT);
#endif



#ifdef PRINT_INFO
	     printf("2.computed U-factor\n");fflush(STDOUT);
	     for (i=0; i<current->nB; ) {
	       if (jlu[n+1+i]==0){
		 for (k=jlu[i];k<jlu[i+1]; k++) 
		   printf("%8d",jlu[k-1]);
		 printf("\n");fflush(STDOUT);
		 for (k=jlu[i];k<jlu[i+1]; k++) 
		   printf("%8.1le",alu[k-1]);
		 printf("\n");fflush(STDOUT);
		 i++;
	       }
	       else {
		 for (k=jlu[i];k<jlu[i+1]; k++) 
		   printf("%8d",jlu[k-1]);
		 printf("\n");fflush(STDOUT);
		 for (k=jlu[i];k<jlu[i+1]; k++) 
		   printf("%8.1le",alu[jlu[i]+2*(k-jlu[i])-1]);
		 printf("\n");fflush(STDOUT);
		 for (k=jlu[i];k<jlu[i+1]; k++) 
		   printf("%8.1le",alu[jlu[i]+2*(k-jlu[i])]);
		 printf("\n");fflush(STDOUT);
		 i+=2;
	       }
	     }
	     printf("\n");fflush(stdout);
	     printf("Block diagonal factor\n");
	     for (k=0; k<current->nB;) {
	       if (jlu[n+1+k]==0){
		 printf("%8.1le",alu[k]);
		 k++;
	       }
	       else {
		 printf("%8.1le%8.1le",alu[k],alu[n+1+k]);
		 k+=2;
	       }
	     }
	     printf("\n");fflush(stdout);
	     for (k=0; k<current->nB; ) {
	       if (jlu[n+1+k]==0) {
		 printf("        ");
		 k++;
	       }
	       else {
		 printf("%8.1le%8.1le",alu[n+1+k],alu[k+1]);
		 k+=2;
	       }
	     }
	     printf("\n");fflush(stdout);
	     printf("computed Schur complement\n");fflush(STDOUT);
	     for (i=current->nB; i<n; i++) {
	       for (k=jlu[i];k<jlu[i+1]; k++) 
		 printf("%8d",jlu[k-1]);
	       printf("\n");fflush(STDOUT);
	       for (k=jlu[i];k<jlu[i+1]; k++) 
		 printf("%8.1le",alu[k-1]);
	       printf("\n");fflush(STDOUT);
	     }
	     for (i=0; i<=current->nB; i++)
	       printf("%8d",i+1);
	     for (i=0; i<=current->nB; i++)
	       printf("%8d",jlu[i]);
	     printf("\n");fflush(stdout);
	     printf("\n");fflush(stdout);
	     for (i=0; i<current->nB; i++)
	       printf("%8d",jlu[n+1+i]);
	     printf("\n");fflush(stdout);
#endif
	    
	  } // end if

	  
	  // maximum (peek) memory observed during sympiluc
	  ILUPACK_mem[4]=MAX(ILUPACK_mem[4],ILUPACK_mem[2]+mymaximem);
	  ILUPACK_mem[5]=MAX(ILUPACK_mem[5],ILUPACK_mem[3]+mymaximem);

	  /* if piluc reduction was successful */
	  if (nB-current->nB<=amgcancel*nB && nB>0) {
	    if (n-current->nB) {
	       MYSYMAMGEXTRACT(&current->F,S, current->p,current->invq,
			       current->nB);
	       ILUPACK_mem[6]+=current->F.nr+1
		              +current->F.ia[current->F.nr]-1;
	       ILUPACK_mem[7]+=current->F.ia[current->F.nr]-1;
	       current->E=current->F;
#ifdef PRINT_INFO2
	       printf("lev %d, F allocated\n",*nlev);fflush(stdout);
#endif
	    }
	    else {
	       current->E.ia=current->F.ia=NULL;
	       current->E.ja=current->F.ja=NULL;
	       current->E.a =current->F.a =NULL;
	       current->E.nr=current->F.nc=0;
	       current->E.nc=current->F.nr=0;
	    }

	    /* discard old S by simply shifting the new S and the L and U factors,
	       because they immediately follow S in the memory */

	    if (*nlev>1) {
	       /* number of nonzeros in the old S */
	       k=S.ia[n]-1;
	       /* number of nonzeros in U and the new S */
	       nnz=jlu[n]-1;
	       alu-=k;
	       jlu-=k;
	       delta[*nlev-1]-=k;
	       for (i=0; i<nnz; i++) {
		   alu[i]=alu[i+k];
		   jlu[i]=jlu[i+k];
	       } /* end for */
	       mem+=k;
	       ILUPACK_mem[2]-=k;
	       ILUPACK_mem[3]-=k;
	    }
	   
	    current->LU.nr=current->nB;
	    current->LU.nc=current->LU.nr;
	    current->LU.a =alu;
	    current->LU.ja=jlu;
	    current->LU.ia=NULL;
	    
	    nnz=jlu[current->nB]-1;
	    S.nr=n-current->nB; 
	    S.nc=S.nr;
	    S.a=alu+nnz;
	    S.ja=jlu+nnz;
	    S.ia=jlu+current->nB;

	    for (i=0; i<=S.nr; i++) {
#ifdef PRINT_INFO
	        printf("!%4d,%4d,%4d\n",i+1,S.ia[i],S.ia[i]-nnz);
#endif
	        S.ia[i]-=nnz;
	    }
	     
	    nnz+=jlu[n]-1;
	    alu+=nnz;
	    jlu+=nnz;
	    delta[*nlev]=delta[*nlev-1]+nnz;
	    mem-=(size_t)nnz;
	    ILUPACK_mem[2]+=(size_t)nnz;
	    ILUPACK_mem[3]+=(size_t)nnz;
	    // maximum (peek) memory observed during sympiluc