ETree/src/util.c

/*  util.c  */

#include "../ETree.h"

#define MYDEBUG 0

/*--------------------------------------------------------------------*/
/*
   ----------------------------------------------
   return the number of bytes taken by the object

   created -- 95nov15, cca
   ----------------------------------------------
*/
int
ETree_sizeOf (
   ETree   *etree
) {
int   nbytes ;
/*
   ---------------
   check the input
   ---------------
*/
if ( etree == NULL ) {
   fprintf(stderr, "\n fatal error in ETree_sizeOf(%p)"
           "\n bad input\n", etree) ;
   exit(-1) ;
}
nbytes = sizeof(struct _ETree) ;
if ( etree->tree != NULL ) {
   nbytes += Tree_sizeOf(etree->tree) ;
}
if ( etree->nodwghtsIV != NULL ) {
   nbytes += IV_sizeOf(etree->nodwghtsIV) ;
}
if ( etree->nodwghtsIV != NULL ) {
   nbytes += IV_sizeOf(etree->bndwghtsIV) ;
}
if ( etree->vtxToFrontIV != NULL ) {
   nbytes += IV_sizeOf(etree->vtxToFrontIV) ;
}
return(nbytes) ; }

/*--------------------------------------------------------------------*/
/*
   ----------------------------------------
   return the number of factor indices

   created  -- 95nov15, cca
   modified -- 96jan11, cca
   ----------------------------------------
*/
int
ETree_nFactorIndices (
   ETree   *etree
) {
int   nb, nfront, nv, nind, nvtx, v ;
int   *bndwghts, *nodwghts ;
/*
   ---------------
   check the input
   ---------------
*/
if ( etree == NULL
   || (nfront = etree->nfront) <= 0
   || (nvtx = etree->nvtx) <= 0 ) {
   fprintf(stderr, "\n fatal error in ETree_nFactorIndices(%p)"
           "\n bad input\n", etree) ;
   exit(-1) ;
}
nodwghts = IV_entries(etree->nodwghtsIV) ;
bndwghts = IV_entries(etree->bndwghtsIV) ;
nind = 0 ;
for ( v = 0 ; v < nfront ; v++ ) {
   nv = nodwghts[v] ;
   nb = bndwghts[v] ;
   nind += nv + nb ;
}
return(nind) ; }

/*--------------------------------------------------------------------*/
/*
   ------------------------------------------
   return the number of factor entries

   symflag -- symmetry flag
     0 (SPOOLES_SYMMETRIC)    -- symmetric
     1 (SPOOLES_HERMITIAN)    -- hermitian
     2 (SPOOLES_NONSYMMETRIC) -- nonsymmetric

   created -- 98jun05, cca
   ------------------------------------------
*/
int
ETree_nFactorEntries (
   ETree   *etree,
   int     symflag
) {
int   J, nfront, nvtx, nzf ;
/*
   ---------------
   check the input
   ---------------
*/
if ( etree == NULL
   || (nfront = etree->nfront) <= 0
   || (nvtx = etree->nvtx) <= 0 ) {
   fprintf(stderr, "\n fatal error in ETree_nFactorEntries(%p,%d)"
           "\n bad input\n", etree, symflag) ;
   exit(-1) ;
}
nzf = 0 ;
for ( J = 0 ; J < nfront ; J++ ) {
   nzf += ETree_nFactorEntriesInFront(etree, symflag, J) ;
}
return(nzf) ; }

/*--------------------------------------------------------------------*/
/*
   ------------------------------------------
   return the number of factor operations

   type -- type of matrix entries
     1 (SPOOLES_REAL)    -- real entries
     2 (SPOOLES_COMPLEX) -- complex entries
   symflag -- symmetry flag
     0 (SPOOLES_SYMMETRIC)    -- symmetric
     1 (SPOOLES_HERMITIAN)    -- hermitian
     2 (SPOOLES_NONSYMMETRIC) -- nonsymmetric

   created -- 98jun05, cca
   ------------------------------------------
*/
double
ETree_nFactorOps (
   ETree   *etree,
   int     type,
   int     symflag
) {
double   ops ;
int      J, nfront, nvtx ;
/*
   ---------------
   check the input
   ---------------
*/
if ( etree == NULL
   || (nfront = etree->nfront) <= 0
   || (nvtx = etree->nvtx) <= 0 ) {
   fprintf(stderr, "\n fatal error in ETree_nFactorOps(%p,%d,%d)"
           "\n bad input\n", etree, type, symflag) ;
   exit(-1) ;
}
ops = 0 ;
for ( J = 0 ; J < nfront ; J++ ) {
   ops += ETree_nInternalOpsInFront(etree, type, symflag, J)
        + ETree_nExternalOpsInFront(etree, type, symflag, J) ;
}
return(ops) ; }

/*--------------------------------------------------------------------*/
/*
   ----------------------------------------
   return the number of entries an LU front

   created -- 96dec04, cca
   ----------------------------------------
*/
double
ETree_nFactorEntriesInFront (
   ETree   *etree,
   int     symflag,
   int     J
) {
int   b, m, nent ;
/*
   ---------------
   check the input
   ---------------
*/
if ( etree == NULL
   || etree->nfront <= 0
   || J < 0 || J >= etree->nfront ) {
   fprintf(stderr,
           "\n fatal error in ETree_nFactorEntriesInFront(%p,%d,%d)"
           "\n bad input\n", etree, symflag, J) ;
   exit(-1) ;
}
b = IV_entry(etree->nodwghtsIV, J) ;
m = IV_entry(etree->bndwghtsIV, J) ;
switch ( symflag ) {
case SPOOLES_SYMMETRIC :
case SPOOLES_HERMITIAN :
   nent = (b*(b+1))/2 + b*m ;
   break ;
case SPOOLES_NONSYMMETRIC :
   nent = b*b + 2*b*m ;
   break ;
default :
   fprintf(stderr,
           "\n fatal error in ETree_nFactorEntriesInFront(%p,%d,%d)"
           "\n bad symflag\n", etree, symflag, J) ;
   break ;
}

return(nent) ; }

/*--------------------------------------------------------------------*/
/*
   -------------------------------------------------------
   return the number of internal LU operations for a front

   created -- 96dec04, cca
   -------------------------------------------------------
*/
double
ETree_nInternalOpsInFront (
   ETree   *etree,
   int     type,
   int     symflag,
   int     J
) {
double   b, m, ops ;
/*
   ---------------
   check the input
   ---------------
*/
if ( etree == NULL
   || etree->nfront <= 0
   || J < 0 || J >= etree->nfront ) {
   fprintf(stderr,
           "\n fatal error in ETree_nInternalOpsInFront(%p,%d,%d,%d)"
           "\n bad input\n", etree, type, symflag, J) ;
   exit(-1) ;
}
b = ETree_frontSize(etree, J) ;
m = ETree_frontBoundarySize(etree, J) ;
switch ( symflag ) {
case SPOOLES_SYMMETRIC :
case SPOOLES_HERMITIAN :
   ops = (b*(b+1)*(2*b+1))/6. + m*b*b ;
   break ;
case SPOOLES_NONSYMMETRIC :
   ops = b*(2*b*b+1)/3. + 2*m*b*b ;
   break ;
default :
   fprintf(stderr,
           "\n fatal error in ETree_nInternalOpsInFront(%p,%d,%d,%d)"
           "\n bad symflag\n", etree, type, symflag, J) ;
   break ;
}
switch ( type ) {
case SPOOLES_REAL :
   break ;
case SPOOLES_COMPLEX :
   ops = 4*ops ;
   break ;
default :
   fprintf(stderr,
           "\n fatal error in ETree_nInternalOpsInFront(%p,%d,%d,%d)"
           "\n bad type\n", etree, type, symflag, J) ;
   break ;
}
return(ops) ; }

/*--------------------------------------------------------------------*/
/*
   -------------------------------------------------------
   return the number of external LU operations for a front

   created -- 96dec04, cca
   -------------------------------------------------------
*/
double
ETree_nExternalOpsInFront (
   ETree   *etree,
   int     type,
   int     symflag,
   int     J
) {
double   b, m, ops ;
/*
   ---------------
   check the input
   ---------------
*/
if ( etree == NULL
   || etree->nfront <= 0
   || J < 0 || J >= etree->nfront ) {
   fprintf(stderr,
           "\n fatal error in ETree_nExternalOpsInFront(%p,%d,%d,%d)"
           "\n bad input\n", etree, J, type, symflag) ;
   exit(-1) ;
}
b = IV_entry(etree->nodwghtsIV, J) ;
m = IV_entry(etree->bndwghtsIV, J) ;
switch ( symflag ) {
case SPOOLES_SYMMETRIC :
case SPOOLES_HERMITIAN :
   ops = m*(m+1)*b ;
   break ;
case SPOOLES_NONSYMMETRIC :
   ops = 2*b*m*m ;
   break ;
default :
   break ;
}
switch ( type ) {
case SPOOLES_REAL :
   break ;
case SPOOLES_COMPLEX :
   ops = 4*ops ;
   break ;
default :
   fprintf(stderr,
           "\n fatal error in ETree_nExternalOpsInFront(%p,%d,%d,%d)"
           "\n bad input\n", etree, J, type, symflag) ;
   break ;
}
return(ops) ; }

/*--------------------------------------------------------------------*/
/*
   ---------------------------------------------------------
   return a DV object that contains the number of operations
   for each front using a backward looking algorithm

   created -- 96dec04, cca
   ---------------------------------------------------------
*/
DV *
ETree_backwardOps (
   ETree   *etree,
   int     type,
   int     symflag,
   int     *vwghts,
   IVL     *symbfacIVL
) {
double   extops, opsKbndK, opsKK ;
double   *ops ;
DV       *opsDV ;
int      bndwghtJ, ii, J, k, K, kwght,
         nadj, nfront, size, wghtJ, wghtK ;
int      *counts, *indices, *list, *mark, *vtxToFront ;
/*
   ---------------
   check the input
   ---------------
*/
if ( etree == NULL || symbfacIVL == NULL ) {
   fprintf(stderr, "\n fatal error in ETree_backwardOps(%p,%p,%p)"
           "\n bad input\n", etree, vwghts, symbfacIVL) ;
   exit(-1) ;
}
nfront     = etree->nfront ;
vtxToFront = IV_entries(etree->vtxToFrontIV) ;
list   = IVinit(nfront, -1) ;
mark   = IVinit(nfront, -1) ;
counts = IVinit(nfront, 0) ;
/*
   ----------------------------
   initialize the ops DV object
   ----------------------------
*/
opsDV = DV_new() ;
DV_init(opsDV, nfront, NULL) ;
ops = DV_entries(opsDV) ;
DV_fill(opsDV, 0.0) ;
/*
   -------------------
   fill the ops vector
   -------------------
*/
for ( J = 0 ; J < nfront ; J++ ) {
   ops[J]  += ETree_nInternalOpsInFront(etree, type, symflag, J) ;
   wghtJ    = ETree_frontSize(etree, J) ;
   bndwghtJ = ETree_frontBoundarySize(etree, J) ;
   IVL_listAndSize(symbfacIVL, J, &size, &indices) ;
   for ( ii = nadj = 0 ; ii < size ; ii++ ) {
      k = indices[ii] ;
      if ( (K = vtxToFront[k]) != J ) {
         kwght = (vwghts == NULL) ? 1 : vwghts[k] ;
         if ( mark[K] != J ) {
            counts[K] = 0 ;
            mark[K] = J ;
            list[nadj++] = K ;
         }
         counts[K] += kwght ;
      }
   }
   IVqsortUp(nadj, list) ;
   extops = 0.0 ;
   for ( ii = 0 ; ii < nadj ; ii++ ) {
      K = list[ii] ;
      wghtK = counts[K] ;
      bndwghtJ -= wghtK ;
      if ( type == SPOOLES_REAL ) {
         opsKbndK = 2*wghtJ*wghtK*bndwghtJ ;
         if ( symflag == SPOOLES_SYMMETRIC ) {
            opsKK = wghtJ*wghtK*(wghtK+1) ;
         } else if ( symflag == SPOOLES_NONSYMMETRIC ) {
            opsKK = 2*wghtJ*wghtK*wghtK ;
         }
      } else if ( type == SPOOLES_COMPLEX ) {
         opsKbndK = 8*wghtJ*wghtK*bndwghtJ ;
         if (  symflag == SPOOLES_SYMMETRIC
            || symflag == SPOOLES_HERMITIAN ) {
            opsKK = 4*wghtJ*wghtK*(wghtK+1) ;
         } else if ( symflag == SPOOLES_NONSYMMETRIC ) {
            opsKK = 8*wghtJ*wghtK*wghtK ;
         }
      }
      extops += opsKK + opsKbndK ;
      ops[K] += opsKK + opsKbndK ;
      if ( symflag == SPOOLES_NONSYMMETRIC ) {
         extops += opsKbndK ;
         ops[K] += opsKbndK ;
      }
   }
/*
   fprintf(stdout,
       "\n front %d, %.0f internal ops, %.0f external ops, %.0f extops",
           J, ETree_nInternalOpsInFront(etree, type, symflag, J),
           ETree_nExternalOpsInFront(etree, type, symflag, J), extops) ;
*/
}
IVfree(list) ;
IVfree(mark) ;
IVfree(counts) ;

return(opsDV) ; }

/*--------------------------------------------------------------------*/
/*
   ------------------------------------
   return an IV object that contains
   the number of entries for each front

   created -- 98jan30, cca
   ------------------------------------
*/
IV *
ETree_factorEntriesIV (
   ETree   *etree,
   int     symflag
) {
int   J, nfront ;
int   *nzf ;
IV    *nzfIV ;
/*
   ---------------
   check the input
   ---------------
*/
if ( etree == NULL ) {
   fprintf(stderr, "\n fatal error in ETree_factorEntriesIV(%p,%d)"
           "\n bad input\n", etree, symflag) ;
   exit(-1) ;
}
nfront = ETree_nfront(etree) ;
/*
   ------------------------
   initialize the IV object
   ------------------------
*/
nzfIV = IV_new() ;
IV_init(nzfIV, nfront, NULL) ;
nzf = IV_entries(nzfIV) ;
IV_fill(nzfIV, 0) ;
/*
   -------------------
   fill the nzf vector
   -------------------
*/
for ( J = 0 ; J < nfront ; J++ ) {
   nzf[J] = ETree_nFactorEntriesInFront(etree, symflag, J) ;
}
return(nzfIV) ; }

/*--------------------------------------------------------------------*/
/*
   ---------------------------------------------------------
   return a DV object that contains the number of operations
   for each front using a forward-looking algorithm

   created -- 96dec04, cca
   ---------------------------------------------------------
*/
DV *
ETree_forwardOps (
   ETree   *etree,
   int     type,
   int     symflag
) {
double   *ops ;
DV       *opsDV ;
int      J, nfront ;
/*
   ---------------
   check the input
   ---------------
*/
if ( etree == NULL ) {
   fprintf(stderr, "\n fatal error in ETree_forwardOps(%p)"
           "\n bad input\n", etree) ;
   exit(-1) ;
}
nfront     = etree->nfront ;
opsDV = DV_new() ;
DV_init(opsDV, nfront, NULL) ;
ops = DV_entries(opsDV) ;
DV_fill(opsDV, 0.0) ;
for ( J = 0 ; J < nfront ; J++ ) {
   ops[J] += ETree_nInternalOpsInFront(etree, type, symflag, J)
          +  ETree_nExternalOpsInFront(etree, type, symflag, J) ;
}
return(opsDV) ; }

/*--------------------------------------------------------------------*/
/*
   ---------------------------------------------------------------
   given an IV object that maps uncompressed vertices to vertices,
   create and return an ETree object that is relative to the
   uncompressed graph.

   created -- 97feb13, cca
   ---------------------------------------------------------------
*/
ETree *
ETree_expand (
   ETree   *etree,
   IV      *eqmapIV
) {
ETree   *etree2 ;
int     ii, ndof, nfront ;
int     *map, *vtxToFront, *vtxToFront2 ;
/*
   ---------------
   check the input
   ---------------
*/
if ( etree == NULL || eqmapIV == NULL ) {
   fprintf(stderr, "\n fatal error in ETree_expand(%p,%p)"
           "\n bad input\n", etree, eqmapIV) ;
   exit(-1) ;
}
nfront = etree->nfront ;
IV_sizeAndEntries(eqmapIV, &ndof, &map) ;
/*
   ---------------------------
   create the new ETree object
   ---------------------------
*/
etree2 = ETree_new() ;
ETree_init1(etree2, nfront, ndof) ;
IV_copy(etree2->nodwghtsIV, etree->nodwghtsIV) ;
IV_copy(etree2->bndwghtsIV, etree->bndwghtsIV) ;
etree2->tree->root = etree->tree->root ;
IVcopy(nfront, etree2->tree->par, etree->tree->par) ;
IVcopy(nfront, etree2->tree->fch, etree->tree->fch) ;
IVcopy(nfront, etree2->tree->sib, etree->tree->sib) ;
vtxToFront  = IV_entries(etree->vtxToFrontIV) ;
vtxToFront2 = IV_entries(etree2->vtxToFrontIV) ;
for ( ii = 0 ; ii < ndof ; ii++ ) {
   vtxToFront2[ii] = vtxToFront[map[ii]] ;
}

return(etree2) ; }

/*--------------------------------------------------------------------*/
/*
   --------------------------------------------------------------
   this method is used to splice together two front trees
   when the domain vertices and schur complement vertices
   have been ordered separately.

   etree0 -- the lower front tree is for vertices in the domains.
   graph0 -- graph for all the vertices
   mapIV  -- IV object that maps vertices to schur complement
             vertices, if IV_entry(mapIV, v) < 0 then v is
             a domain vertex.
   etree1 -- the upper front tree is for vertices in the schur
             complement.

   created -- 97feb01, cca
   --------------------------------------------------------------
*/
ETree *
ETree_spliceTwoETrees (
   ETree   *etree0,
   Graph   *graph0,
   IV      *mapIV,
   ETree   *etree1
) {
ETree   *etree2 ;
int     *bndwghts0, *bndwghts1, *bndwghts2, *fch0, *head0, *link0,
        *map, *mark, *nodwghts0, *nodwghts1, *nodwghts2, *par0,
        *par1, *par2, *sib0, *vadj, *vtxToFront0, *vtxToFront1,
        *vtxToFront2 ;
int     ii, J, K, nfront0, nfront1, nfront2, nvtx, phi, v, vsize, w ;
/*
   ---------------
   check the input
   ---------------
*/
if (  etree0 == NULL || graph0 == NULL
   || mapIV == NULL || etree1 == NULL ) {
   fprintf(stderr,
           "\n fatal error in ETree_spliceTwoETrees(%p,%p,%p,%p)"
           "\n bad input\n",
           etree0, graph0, mapIV, etree1) ;
   exit(-1) ;
}
nfront0     = etree0->nfront ;
nvtx        = etree0->nvtx    ;
par0        = etree0->tree->par ;
fch0        = etree0->tree->fch ;
sib0        = etree0->tree->sib ;
nodwghts0   = IV_entries(etree0->nodwghtsIV) ;
bndwghts0   = IV_entries(etree0->bndwghtsIV) ;
vtxToFront0 = IV_entries(etree0->vtxToFrontIV) ;
nfront1     = etree1->nfront ;
par1        = etree1->tree->par ;
bndwghts1   = IV_entries(etree1->bndwghtsIV) ;
nodwghts1   = IV_entries(etree1->nodwghtsIV) ;
vtxToFront1 = IV_entries(etree1->vtxToFrontIV) ;
map         = IV_entries(mapIV) ;
/*
   -------------------------
   create the new front tree
   -------------------------
*/
nfront2 = nfront0 + nfront1 ;
etree2 = ETree_new() ;
ETree_init1(etree2, nfront2, etree0->nvtx) ;
par2        = etree2->tree->par ;
nodwghts2   = IV_entries(etree2->nodwghtsIV) ;
bndwghts2   = IV_entries(etree2->bndwghtsIV) ;
vtxToFront2 = IV_entries(etree2->vtxToFrontIV) ;
/*
   --------------------------------------------------
   fill the parent fields for fronts in the same tree
   --------------------------------------------------
*/
for ( J = 0 ; J < nfront0 ; J++ ) {
   par2[J]      = par0[J] ;
   nodwghts2[J] = nodwghts0[J] ;
   bndwghts2[J] = bndwghts0[J] ;
}
for ( J = 0 ; J < nfront1 ; J++ ) {
   par2[J+nfront0]      = nfront0 + par1[J] ;
   nodwghts2[J+nfront0] = nodwghts1[J] ;
   bndwghts2[J+nfront0] = bndwghts1[J] ;
}
/*
   ---------------------------
   set the vertex to front map
   ---------------------------
*/
for ( v = 0 ; v < nvtx ; v++ ) {
   if ( (J = vtxToFront0[v]) >= 0 ) {
      vtxToFront2[v] = J ;
   } else {
      vtxToFront2[v] = vtxToFront1[map[v]] + nfront0 ;
   }
}
/*
   ---------------------------------------------
   link the vertices to fronts in the lower tree
   ---------------------------------------------
*/
head0 = IVinit(nfront0, -1) ;
link0 = IVinit(nvtx,    -1) ;
for ( v = 0 ; v < nvtx ; v++ ) {
   if ( (J = vtxToFront0[v]) >= 0 ) {
      link0[v] = head0[J] ;
      head0[J] = v ;
   }
}
/*
   -------------------------------------------------------
   link roots of the lower tree to nodes in the upper tree
   -------------------------------------------------------
*/
mark = IVinit(nvtx, -1) ;
for ( J = etree0->tree->root ; J != -1 ; J = sib0[J] ) {
/*
   ---------------------------------------
   K is the parent front in the upper tree
   ---------------------------------------
*/
   K = nfront1 ;
/*
   ---------------------------------------
   loop over vertices in the lower front J
   ---------------------------------------
*/
   for ( v = head0[J] ; v != -1 ; v = link0[v] ) {
      Graph_adjAndSize(graph0, v, &vsize, &vadj) ;
      for ( ii = 0 ; ii < vsize ; ii++ ) {
         w = vadj[ii] ;
         if ( vtxToFront0[w] < 0 ) {
            phi = map[w] ;
/*
            ---------------------------------------------------
            w is a vertex that belongs to phi in the upper tree
            ---------------------------------------------------
*/
            if ( mark[phi] != J ) {
               mark[phi] = J ;
               if ( K > vtxToFront1[phi] ) {
/*
                  ------------------------------------
                  least numbered adjacent front so far
                  ------------------------------------
*/
                  K = vtxToFront1[phi] ;
               }
            }
         }
      }
   }
   if ( K < nfront1 ) {
   /*
      --------------------
      set the parent field
      --------------------
   */
      par2[J] = nfront0 + K ;
   }
}
/*
   -----------------------------
   set the remaining tree fields
   -----------------------------
*/
Tree_setFchSibRoot(etree2->tree) ;
/*
   ------------------------
   free the working storage
   ------------------------
*/
IVfree(head0) ;
IVfree(link0) ;
IVfree(mark)  ;

return(etree2) ; }

/*--------------------------------------------------------------------*/