symamgsetup.c

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

	    ILUPACK_mem[4]=MAX(ILUPACK_mem[4],ILUPACK_mem[2]);
	    ILUPACK_mem[5]=MAX(ILUPACK_mem[5],ILUPACK_mem[3]);


#ifdef PRINT_INFO
	    printf("3.computed U-factor\n");fflush(STDOUT);
	    for (i=0; i<current->nB; ) {
	      if (current->LU.ja[n+1+i]==0){
		for (k=current->LU.ja[i];k<current->LU.ja[i+1]; k++) 
		  printf("%8d",current->LU.ja[k-1]);
		printf("\n");fflush(STDOUT);
		for (k=current->LU.ja[i];k<current->LU.ja[i+1]; k++) 
		  printf("%8.1le",current->LU.a[k-1]);
		printf("\n");fflush(STDOUT);
		i++;
	      }
	      else {
		for (k=current->LU.ja[i];k<current->LU.ja[i+1]; k++) 
		  printf("%8d",current->LU.ja[k-1]);
		printf("\n");fflush(STDOUT);
		for (k=current->LU.ja[i];k<current->LU.ja[i+1]; k++) 
		  printf("%8.1le",current->LU.a[current->LU.ja[i]+2*(k-current->LU.ja[i])-1]);
		printf("\n");fflush(STDOUT);
		for (k=current->LU.ja[i];k<current->LU.ja[i+1]; k++) 
		  printf("%8.1le",current->LU.a[current->LU.ja[i]+2*(k-current->LU.ja[i])]);
		printf("\n");fflush(STDOUT);
		i+=2;
	      }
	    }
	    printf("\n");fflush(stdout);
	    printf("Block diagonal factor\n");
	    for (k=0; k<current->nB;) {
	      if (current->LU.ja[n+1+k]==0){
		printf("%8.1le",current->LU.a[k]);
		k++;
	      }
	      else {
		printf("%8.1le%8.1le",current->LU.a[k],current->LU.a[n+1+k]);
		k+=2;
	      }
	    }
	    printf("\n");fflush(stdout);
	    for (k=0; k<current->nB; ) {
	      if (current->LU.ja[n+1+k]==0) {
		printf("        ");
		k++;
	      }
	      else {
		printf("%8.1le%8.1le",current->LU.a[n+1+k],current->LU.a[k+1]);
		k+=2;
	      }
	    }
	    printf("\n");fflush(stdout);
	    printf("computed Schur complement\n");fflush(STDOUT);
	    for (i=0; i<S.nr; i++) {
	      for (k=S.ia[i];k<S.ia[i+1]; k++) 
		printf("%8d",S.ja[k-1]);
	      printf("\n");fflush(STDOUT);
	      for (k=S.ia[i];k<S.ia[i+1]; k++) 
		printf("%8.1le",S.a[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",current->LU.ja[i]);
	    printf("\n");fflush(stdout);
	    for (i=0; i<current->nB; i++)
	      printf("%8d",current->LU.ja[n+1+i]);
	    printf("\n");fflush(stdout);
#endif


	  } /* end if (nB-current->nB<=amgcancel*nB) */
	  else { /* sympiluc failed to produce a sensible reduction */

#ifdef PRINT_INLINE
	     printf("level %d,  piluc failed (nB=%d), switch to SYMILUC\n",*nlev,current->nB);fflush(STDOUT);
#endif

	     /* switch to ILDLC */

	     // start counter for ILDLC
	     evaluate_time(&time_start,&systime);
	     current->LUperm=NULL;
	     for (i=0; i<n; i++)
	         current->p[i]=current->invq[i]=i+1;
	     
	     imem=(LONG_INT)mem;
	     myimem=imem;
	     SYMILUC(&n,S.a,S.ja,S.ia,&n,IPparam->fpar+7,&PILUCparam,
		     current->p,current->invq,alu,jlu,&myimem,
		     IPparam->dbuff,IPparam->ibuff,&ierr);
	     ILUPACK_mem[4]=MAX(ILUPACK_mem[4],ILUPACK_mem[2]+myimem);
	     ILUPACK_mem[5]=MAX(ILUPACK_mem[5],ILUPACK_mem[3]+myimem);

	     /*
	     for (i=0; i<n; i++) 
	       printf("%8.1le\n",alu[i]);
	     printf("\n");fflush(stdout);
	     for (i=0; i<n; i++) {
	       printf("%8d:\n",i+1);fflush(stdout);
	       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);
	     }
	     */

	     // computation time for ILDLC
	     evaluate_time(&time_stop,&systime);
	     ILUP_sec[3]=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;
	     
	     /* discard old S by simply shifting the new U factor, because
		it immediately follows the old S in the memory */
	     if (*nlev>1) {
	        /* number of nonzeros in the old S */
	        k=S.ia[n]-1;
		/* number of nonzeros in U */
		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=(size_t)jlu[n]-1;
	     ILUPACK_mem[2]+=(size_t)nnz;
	     ILUPACK_mem[3]+=(size_t)nnz;
	     ILUPACK_mem[4]=MAX(ILUPACK_mem[4],ILUPACK_mem[2]);
	     ILUPACK_mem[5]=MAX(ILUPACK_mem[5],ILUPACK_mem[3]);

	     n=-1;
	     S.nr=0;
	  } /* end if-else */
	} /* end if-else (*nlev>1 && 3*(S.ia[n]-1)>S.nr*S.nc) */
  } /* end while */
  
  // reduce memory to the size that was really needed
  if (!(param&RE_FACTOR)) {
    // printf("matrix will not be re-factored again\n");fflush(stdout);
     IPparam->njlu=ILUPACK_mem[2];
     IPparam->jlu=(integer *)  REALLOC(IPparam->jlu,
				   IPparam->njlu*sizeof(integer),  
				   "symAMGsetup:jlu");
     IPparam->nalu=ILUPACK_mem[3];
     IPparam->alu=(FLOAT *)REALLOC(IPparam->alu,
				   IPparam->nalu*sizeof(FLOAT),
				   "symAMGsetup:alu");

     /*
     current=PRE;
     j=0;
     for (i=0; i<*nlev; i++) {
       if (i==0)
	 printf("%8d,%8ld\n", delta[i],(long)0);
       else {
	 j+=(integer)(((unsigned long)(current->LU.a)-(unsigned long)(current->prev->LU.a))/sizeof(FLOAT));
	 printf("%8d,%8ld\n", delta[i],j);
       }
       current=current->next;
     }
     if ((unsigned long)IPparam->alu!=(unsigned long)PRE->LU.a)
       printf("pointers will be remapped\n");
     */

     // remap pointers
     current=PRE;
     for (i=0; i<*nlev; i++) {

       // sparse case
       if (current->LU.ja!=NULL)
	  current->LU.ja=IPparam->jlu+delta[i];
       else // dense case
	  current->LU.ia=current->prev->LU.ja+current->prev->LU.nr;

       current->LU.a =IPparam->alu+delta[i];
       
       current=current->next;
     }

     /*
     current=PRE;
     for (i=1; i<*nlev; i++) {
       current=current->next;
     }
     for (i=0; i<current->LU.nr; i++) 
       printf("%8.1le\n",current->LU.a[i]);
     printf("\n");fflush(stdout);
     for (i=0; i<current->LU.nr; i++) {
       printf("%8d:\n",i+1);fflush(stdout);
       for (k=current->LU.ja[i]; k<current->LU.ja[i+1];k++)
	 printf("%8d",current->LU.ja[k-1]);
       printf("\n");fflush(stdout);
       for (k=current->LU.ja[i]; k<current->LU.ja[i+1];k++)
	 printf("%8.1le",current->LU.a[k-1]);
       printf("\n");fflush(stdout);
     }
     */

#ifdef PRINT_MEM
     printf("jlu(%8.2lfMB), alu(%8.2lfMB) finally re-allocated\n",ILUPACK_mem[2]/IMEGA,ILUPACK_mem[3]/DMEGA);fflush(stdout);
#endif
  }
  else {
    // printf("matrix may be re-factored\n");fflush(stdout);
  }

  current=PRE;
  n=A->nr;
  for (j=1; j<=*nlev; j++) {
      // level i
    
      // case of a sparse inverse-based block ILU
      if (j<*nlev || current->LU.ja!=NULL) {
	 // shift block diagonal parts properly
	 if (j==*nlev-1 && current->next->LU.ja==NULL) {
	    /* exceptional case, if the next level is the final level
	       and the final level uses a dense factorization 
	       In this case move the space in ja at nB+1,...,n to 
               n+1+nB+1,...n+1+n
	    */
	    for (i=0; i<n-current->nB; i++) {
		current->LU.ja[n+1+current->nB+i]=current->LU.ja[current->nB+i];
	    }
	    // re-adjust pointer
	    current->next->LU.ia=current->LU.ja+n+1+current->nB;
	 }
	 for (i=0; i<current->nB; i++) {
	     current->LU.ja[current->nB+1+i]=current->LU.ja[n+1+i];
	     current->LU.a[current->nB+1+i] =current->LU.a[n+1+i];
	 } // end for i
	 // now we do not need to refer to the dense Schur complement anymore
	 current->LU.ja[current->nB]=current->LU.ja[current->nB-1];

#ifdef PRINT_INFO
	    printf("4.computed U-factor\n");fflush(STDOUT);
	    for (i=0; i<current->nB; ) {
	      if (current->LU.ja[current->nB+1+i]==0){
		for (k=current->LU.ja[i];k<current->LU.ja[i+1]; k++) 
		  printf("%8d",current->LU.ja[k-1]);
		printf("\n");fflush(STDOUT);
		for (k=current->LU.ja[i];k<current->LU.ja[i+1]; k++) 
		  printf("%8.1le",current->LU.a[k-1]);
		printf("\n");fflush(STDOUT);
		i++;
	      }
	      else {
		for (k=current->LU.ja[i];k<current->LU.ja[i+1]; k++) 
		  printf("%8d",current->LU.ja[k-1]);
		printf("\n");fflush(STDOUT);
		for (k=current->LU.ja[i];k<current->LU.ja[i+1]; k++) 
		  printf("%8.1le",current->LU.a[current->LU.ja[i]+2*(k-current->LU.ja[i])-1]);
		printf("\n");fflush(STDOUT);
		for (k=current->LU.ja[i];k<current->LU.ja[i+1]; k++) 
		  printf("%8.1le",current->LU.a[current->LU.ja[i]+2*(k-current->LU.ja[i])]);
		printf("\n");fflush(STDOUT);
		i+=2;
	      }
	    }
	    printf("\n");fflush(stdout);
	    printf("Block diagonal factor\n");
	    for (k=0; k<current->nB;) {
	      if (current->LU.ja[current->nB+1+k]==0){
		printf("%8.1le",current->LU.a[k]);
		k++;
	      }
	      else {
		printf("%8.1le%8.1le",current->LU.a[k],current->LU.a[current->nB+1+k]);
		k+=2;
	      }
	    }
	    printf("\n");fflush(stdout);
	    for (k=0; k<current->nB; ) {
	      if (current->LU.ja[current->nB+1+k]==0) {
		printf("        ");
		k++;
	      }
	      else {
		printf("%8.1le%8.1le",current->LU.a[current->nB+1+k],current->LU.a[k+1]);
		k+=2;
	      }
	    }
	    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",current->LU.ja[i]);
	    printf("\n");fflush(stdout);
	    for (i=0; i<current->nB; i++)
	      printf("%8d",current->LU.ja[current->nB+1+i]);
	    printf("\n");fflush(stdout);
#endif


      }
      else {
	  // full matrix processing in the final step
#ifdef PRINT_INFO
	    printf("dense U-factor\n");fflush(STDOUT);
	    for (i=0; i<current->nB; i++)
	      printf("%8d",current->LU.ia[i]);
	    printf("\n");fflush(stdout);
	    j=0;
	    for (k=0; k<current->nB; k++) {
	        for (i=k; i<current->nB; i++) {
		    printf("%8.1le",current->LU.a[j++]);
		}
		printf("\n");fflush(stdout);
	    }
	    printf("\n");fflush(stdout);
#endif

      }
      n-=current->nB;
      current=current->next;
  } // end for j

  evaluate_time(&time_stop,&systime);
  // compute total setup time
  ILUP_sec[7]=time_stop-ILUP_sec[7];
  for (i=0; i<ILUPACK_secnds_length; i++)
      ILUPACK_secnds[i]=ILUP_sec[i];

  // export final condest
  IPparam->fpar[2]=condest;

  ILUPACK_mem[0]=IPparam->nibuff;
  ILUPACK_mem[1]=IPparam->ndbuff;

  /*
  if (ierr==0) {
     printf("maximum elbow required during factorization %8.1f\n",
	    MAX(ILUPACK_mem[4],ILUPACK_mem[5])/(A->ia[A->nc]-1.0)+.05);
     printf("maximum buffer size  %8.1f\n",
	    MAX(ILUPACK_mem[0]*sizeof(int)*1.0/sizeof(FLOAT),ILUPACK_mem[1])/(A->ia[A->nc]-1.0)+.05);
     printf("maximum size F       %8.1f\n",
	    MAX(ILUPACK_mem[6]*sizeof(int)*1.0/sizeof(FLOAT),ILUPACK_mem[7])/(A->ia[A->nc]-1.0)+.05);
     printf("maximum size scaling %8.1f\n",
	    MAX(ILUPACK_mem[8]*sizeof(int)*1.0/sizeof(FLOAT),ILUPACK_mem[9])/(A->ia[A->nc]-1.0)+.05);
     printf("maximum size drivers %8.1f\n",
	    MAX(ILUPACK_mem[10]*sizeof(int)*1.0/sizeof(FLOAT),ILUPACK_mem[11])/(A->ia[A->nc]-1.0)+.05);
     fflush(stdout);
  }
  */

  free(delta);
  return (ierr);
} /* end symAMGsetup */