symamgsetup.c

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

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


#include <lapack.h>
#include <ilupack.h>
#include <ilupackmacros.h>


//#define PRINT_MEM
#define MEGA      1048576.0
#define IMEGA      262144.0
#define DMEGA      131072.0

#if defined _DOUBLE_REAL_ || defined _SINGLE_REAL_
#define CONJG(A)      (A)

#ifdef _DOUBLE_REAL_
#define SYMPILUC             DSYMpiluc
#define SYMILUC              DSYMiluc
#define MYSMATVEC            DSYMmatvec
#define MYSYMAMGFACTOR       DSYMAMGfactor
#define MYSYMAMGDELETE       DSYMAMGdelete
#define MYSYMAMGEXTRACT      DSYMAMGextract
#define MYSPTRF              dsptrf

#else
#define SYMPILUC             SSYMpiluc
#define SYMILUC              SSYMiluc
#define MYSMATVEC            SSYMmatvec
#define MYSYMAMGFACTOR       SSYMAMGfactor
#define MYSYMAMGDELETE       SSYMAMGdelete
#define MYSYMAMGEXTRACT      SSYMAMGextract
#define MYSPTRF              ssptrf
#endif



#else



#ifdef _COMPLEX_SYMMETRIC_
#define CONJG(A)     (A)

#ifdef _SINGLE_COMPLEX_
#define SYMPILUC             CSYMpiluc
#define SYMILUC              CSYMiluc
#define MYSMATVEC            CSYMmatvec
#define MYSYMAMGFACTOR       CSYMAMGfactor
#define MYSYMAMGDELETE       CSYMAMGdelete
#define MYSYMAMGEXTRACT      CSYMAMGextract
#define MYSPTRF              csptrf
#else
#define SYMPILUC             ZSYMpiluc
#define SYMILUC              ZSYMiluc
#define MYSMATVEC            ZSYMmatvec
#define MYSYMAMGFACTOR       ZSYMAMGfactor
#define MYSYMAMGDELETE       ZSYMAMGdelete
#define MYSYMAMGEXTRACT      ZSYMAMGextract
#define MYSPTRF              zsptrf
#endif

#else

#ifdef _SINGLE_COMPLEX_
#define CONJG(A)     (conjg(A))
#define SYMPILUC             CHERpiluc
#define SYMILUC              CHERiluc
#define MYSMATVEC            CHERmatvec
#define MYSYMAMGFACTOR       CHERAMGfactor
#define MYSYMAMGDELETE       CHERAMGdelete
#define MYSYMAMGEXTRACT      CHERAMGextract
#define MYSPTRF              chptrf
#else
#define CONJG(A)     (dconjg(A))
#define SYMPILUC             ZHERpiluc
#define SYMILUC              ZHERiluc
#define MYSMATVEC            ZHERmatvec
#define MYSYMAMGFACTOR       ZHERAMGfactor
#define MYSYMAMGDELETE       ZHERAMGdelete
#define MYSYMAMGEXTRACT      ZHERAMGextract
#define MYSPTRF              zhptrf
#endif

#endif



#endif



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

//#define PRINT_INLINE
//#define PRINT_INFO2





/*
  main driver for the multilevel ILU

  given a symmetrically structured A, a sequence of permutations and 
  incomplete LU factorizations is performed. Partially we obtain an
  approximate partial factorization
          /  B  F \   /U^T   0\ /D  0\ /U UF\
  P^TAP = |       | ~ |       | |    | |    |
          \ F^T C /   \UF^T  I/ \0  S/ \0 I /
  The D,U as well as p,p^{-1} and F are kept and the multilevel strategy
  is repeated for S
*/
integer MYSYMAMGFACTOR(CSRMAT *A, AMGLEVELMAT *PRE, integer *nlev, 
		  ILUPACKPARAM *IPparam,
		  integer (*perm0)(CSRMAT, FLOAT *, FLOAT *, integer *, integer *, integer *, 
			       ILUPACKPARAM *),
		  integer (*perm) (CSRMAT, FLOAT *, FLOAT *, integer *, integer *, integer *, 
			       ILUPACKPARAM *),
		  integer (*permf)(CSRMAT, FLOAT *, FLOAT *, integer *, integer *, integer *, 
			       ILUPACKPARAM *))
{
/*
  A        input matrix, matrix is altered by the function
           diagonal scaling is directly applied to A
           only the diagonal part of A and its upper triangular part
           are stored
  PRE      Linked list of partial preconditioners
  nlev     number of levels(blocks)
  IPparam  ILUPACK parameter list
  perm0    function for initial scaling and preprocessing, applied to
           the initial system only
  perm     regular reordering/scaling applied to all subsequent levels
  permf    final scaling/pivoting strategy, applied when perm has failed
    
           perm...(A,rowscale,colscale,p,invq,nB,IPparam)
           - rescales A, returns the results in rowscale and colscale
           - permutes the rows of A and returns the permutation vector in p
           - permutes similarly the columns of A and returns the inverse 
           - permutation vector of p in invq
           - return the size nB of the leading block B
                     /  B  F \
             P^TAP = |       |
                     \ F^T C /
             that hopefully can be safely factored

  Code written by Matthias Bollhoefer November 2003
*/


  integer i,j,k,n=A->nr, PILUCparam, nnz, nB;
  integer  *jlu, ierr=0, regular_pivoting, *ia, *ja;
  integer lfil[2], param, discardA, shiftA;
  REALS droptols[3], condest, condest0, condest1, amgcancel, ilucancel; 
  REALS seps, ptol, elbow;
  CSRMAT S;
  AMGLEVELMAT *current;
  FLOAT *alu, *a, *r,*c;
  float  systime, time_start, time_stop, secnds, 
         ILUP_sec[ILUPACK_secnds_length];
  LONG_INT imem, myimem,mymaximem;
  size_t   nibuff, ndbuff, mem, *delta;




  // init time counters
  evaluate_time(&time_start,&systime);
  // counter for the total setup time
  for (i=0; i<ILUPACK_secnds_length; i++)
      ILUP_sec[i]=0;
  ILUP_sec[7]=time_start;

  // init memory counters
  for (i=0; i<ILUPACK_mem_length; i++)
      ILUPACK_mem[i]=0;


  // for pildlc the integer buffer is required to have at least
  //    11*n for the simple Schur complement (S version)
  //    14*n for the Tismenetsky-like Schur complement (T version)
  // for mpiluc the integer buffer is required to have at least
  //    12*n for the simple Schur complement (S version)
  //    15*n for the Tismenetsky-like Schur complement (T version)

  // analogous settings for the floating point buffer
  //    2*n  (not improved version)
  //    3*n  (improved version or medium Schur complement version)


  elbow  =IPparam->elbow;
  nibuff =IPparam->nibuff;
  ndbuff =IPparam->ndbuff;
  lfil[0]=IPparam->ipar[3];
  lfil[1]=IPparam->ipar[9];
  param  =IPparam->ipar[6];

  droptols[0]=IPparam->fpar[0];
  droptols[1]=IPparam->fpar[1];
  droptols[2]=IPparam->fpar[8];
  condest    =IPparam->fpar[2];
  condest0=condest;
  ilucancel  =IPparam->fpar[3];
  amgcancel  =IPparam->fpar[4]; 

  //  if (param&RE_FACTOR)
  //   printf("entering ILU, memory should be kept\n");fflush(stdout);


  // analyze how the size of the buffer has to be chosen
  if (param&TISMENETSKY_SC) {
     if (param&MULTI_PILUC)
        nibuff=MAX(nibuff,17*(size_t)n);
     nibuff=MAX(nibuff,16*(size_t)n);
  }
  else {
     if (param&MULTI_PILUC)
        nibuff=MAX(nibuff,14*(size_t)n);
     nibuff=MAX(nibuff,13*(size_t)n); // recommended version!
  }

  if (param&DROP_INVERSE) {
     // if a more reliable estimate for the inverse is required
     if (param&IMPROVED_ESTIMATE || !(param&(SIMPLE_SC|TISMENETSKY_SC)))
        ndbuff=MAX(ndbuff,4*(size_t)n); // recommended version!
     else 
        ndbuff=MAX(ndbuff,3*(size_t)n);
  }
  else {
     // currently the norm of the inverse factors are computed
     // even if classical dropping is desired
     // The comments in (m)piluc still refer to their predecessor
     // iluc, which allows classical dropping in a strict sense

     // if a more reliable estimate for the inverse is required
     if (param&IMPROVED_ESTIMATE || !(param&(SIMPLE_SC|TISMENETSKY_SC)))
        ndbuff=MAX(ndbuff,4*(size_t)n);
     else 
        ndbuff=MAX(ndbuff,3*(size_t)n);
  }


  /* integer and floating point buffer */
#ifdef PRINT_MEM
  printf("ibuffmem=%8.2fMB, dbuffmem=%8.2fMB\n",IPparam->nibuff/IMEGA,IPparam->ndbuff/DMEGA);fflush(stdout);
#endif
  if (IPparam->nibuff==0) {
     // printf("no ibuff, allocating memory for the preconditioner\n");fflush(stdout);
     IPparam->nibuff=nibuff;
     IPparam->ibuff=(integer *)  MALLOC(IPparam->nibuff*sizeof(integer),
				    "symAMGsetup:ibuff");
#ifdef PRINT_MEM
     printf("ibuff(%8.2fMB) allocated\n",nibuff/IMEGA);fflush(stdout);
#endif
  }     
  else if (nibuff>IPparam->nibuff) {
     IPparam->nibuff=nibuff;
     IPparam->ibuff=(integer *)  REALLOC(IPparam->ibuff,
				     IPparam->nibuff*sizeof(integer),
				     "symAMGsetup:ibuff");
#ifdef PRINT_MEM
     printf("ibuff(%8.2fMB) re-allocated\n",nibuff/IMEGA);fflush(stdout);
#endif
  }
  if (IPparam->ndbuff==0) {
     IPparam->ndbuff=ndbuff;
     // printf("no dbuff, allocating memory for the preconditioner\n");fflush(stdout);
     IPparam->dbuff=(FLOAT *)MALLOC(IPparam->ndbuff*sizeof(FLOAT),
				    "symAMGsetup:dbuff");
#ifdef PRINT_MEM
     printf("dbuff(%8.2fMB) allocated\n",ndbuff/DMEGA);fflush(stdout);
#endif
  } 
  else if (ndbuff>IPparam->ndbuff) {
     IPparam->ndbuff=ndbuff;
     IPparam->dbuff=(FLOAT *)REALLOC(IPparam->dbuff,
				     IPparam->ndbuff*sizeof(FLOAT),
				     "symAMGsetup:dbuff");
#ifdef PRINT_MEM
     printf("dbuff(%8.2fMB) re-allocated\n",ndbuff/DMEGA);fflush(stdout);
#endif
  }

  // keep track on the amount of memory consumed by the buffers
  ILUPACK_mem[0]=nibuff;
  ILUPACK_mem[1]=ndbuff;



  // total amount of memory requested by the parameters
  mem=elbow*(size_t)A->ia[A->nr];
#ifdef PRINT_MEM
  printf("jlumem=%8.2fMB, alumem=%8.2fMB\n",IPparam->njlu/IMEGA,IPparam->nalu/DMEGA);fflush(stdout);
#endif
  if (IPparam->njlu==0) {
     IPparam->njlu=mem;
     // printf("no jlu, allocating memory for the preconditioner\n");fflush(stdout);
     IPparam->jlu=(integer *)  MALLOC(IPparam->njlu*sizeof(integer),
				  "symAMGsetup:jlu");
#ifdef PRINT_MEM
     printf("jlu(%8.2lfMB) allocated\n",mem/IMEGA);fflush(stdout);
#endif
  }
  else if (mem>IPparam->njlu) {
     IPparam->njlu=mem;
     IPparam->jlu=(integer *) REALLOC(IPparam->jlu,
				  IPparam->njlu*sizeof(integer),
				  "symAMGsetup:jlu");
#ifdef PRINT_MEM
     printf("jlu(%8.2fMB) re-allocated\n",mem/IMEGA);fflush(stdout);
#endif
  }
  jlu=IPparam->jlu;
  if (IPparam->nalu==0) {
     IPparam->nalu=mem;
     // printf("no alu, allocating memory for the preconditioner\n");fflush(stdout);
     IPparam->alu=(FLOAT *) MALLOC(IPparam->nalu*sizeof(FLOAT),
				   "symAMGsetup:alu");
#ifdef PRINT_MEM
     printf("alu(%8.2fMB) allocated\n",mem/DMEGA);fflush(stdout);
#endif
  }
  else if (mem>IPparam->nalu) {
     IPparam->nalu=mem;
     IPparam->alu=(FLOAT *)REALLOC(IPparam->alu,
				   IPparam->nalu*sizeof(FLOAT),
				   "symAMGsetup:alu");
#ifdef PRINT_MEM
     printf("alu(%8.2fMB) re-allocated\n",mem/DMEGA);fflush(stdout);
#endif
  }
  alu=IPparam->alu;



  if (param&FINAL_PIVOTING)
     regular_pivoting=-1;
  else 
     regular_pivoting=0;
  // default settings for PILUC within ARMS oo PILUC
  PILUCparam=COARSE_REDUCE; //DROP_INVERSE|
  if (param&IMPROVED_ESTIMATE)
     PILUCparam|=IMPROVED_ESTIMATE;
  if (param&DROP_INVERSE)
     PILUCparam|=DROP_INVERSE;
  // additional legal options
  if (param&TISMENETSKY_SC)
     PILUCparam|=TISMENETSKY_SC;
  else if (param&SIMPLE_SC)
    PILUCparam|=SIMPLE_SC;
  if (param&DIAGONAL_COMPENSATION)
     PILUCparam|=DIAGONAL_COMPENSATION;

  if (param&AGGRESSIVE_DROPPING)
     PILUCparam|=AGGRESSIVE_DROPPING;

  discardA=0;
  if (param&DISCARD_MATRIX) {
     discardA=DISCARD_MATRIX;
  }

  // store the logical gap between the allocated memory areas alu,jlu
  // and the pointers to the subsystems
  delta=(size_t *)MALLOC(A->nr*sizeof(size_t),"symAMGsetup:delta");
  for (i=0; i<A->nr; i++)
    delta[i]=0;
  /* block LU decomposition */
  S=*A; 

  n=(1+1/amgcancel)*S.nr;
  current=PRE;
  current->prev=NULL;
  *nlev=0;

  
  // while (0<S.nr & S.nr<=amgcancel*n) {
  while (0<S.nr) {
	(*nlev)++;
	
	// printf("lev %d\n",*nlev);fflush(stdout);

        n=S.nr;
	if (*nlev>1)
	{
	   current->next=(AMGLEVELMAT *)MALLOC((size_t)1*sizeof(AMGLEVELMAT),
					       "symAMGsetup:current->next");
#ifdef PRINT_INFO2
	   printf("lev %d, next level allocated\n",*nlev);fflush(stdout);
#endif
	   current->next->prev=current;
	   current=current->next;
	}
	current->n=n;
	current->next=NULL;
	current->absdiag=NULL;
	

	// start counter for preprocessing/reordering/pivoting  and  scaling
	evaluate_time(&time_start,&systime);
	/* row and column scaling if desired */
	current->colscal=(FLOAT *)MALLOC((size_t)n*sizeof(FLOAT),
					 "symAMGsetup:current->colscal");
	current->rowscal=current->colscal;
	ILUPACK_mem[9]+=n;