FrontMtx/src/factorUtil.c

/*  factorUtil.c  */

#include "../FrontMtx.h"
#include "../../timings.h"

/*--------------------------------------------------------------------*/
/*
   -----------------------------
   purpose -- initialize a front

   created -- 98may04, cca
   -----------------------------
*/
void
FrontMtx_initializeFront (
   FrontMtx   *frontmtx,
   Chv        *frontJ,
   int        J
) {
int   ncolJ, nD, nrowJ ;
int   *colindJ, *ivec ;
/*
   -----------------------------------
   get the number of internal vertices
   -----------------------------------
*/
nD = ETree_frontSize(frontmtx->frontETree, J) ;
/*
   --------------------------------------
   get the internal and boundary vertices
   --------------------------------------
*/
IVL_listAndSize(frontmtx->symbfacIVL, J, &ncolJ, &colindJ) ;
#if MYDEBUG > 0
fprintf(stdout, "\n front %d, column indices", J) ;
IVfprintf(stdout, ncolJ, colindJ) ;
fflush(stdout) ;
#endif
/*
   ------------------------------------------------------
   initialize the front's dimensions, indices and entries
   ------------------------------------------------------
*/
#if MYDEBUG > 0
fprintf(stdout, "\n front %d, nD %d, nU %d", J, nD, ncolJ - nD) ;
fflush(stdout) ;
#endif
Chv_init(frontJ, J, nD, ncolJ - nD, ncolJ - nD,
         frontmtx->type, frontmtx->symmetryflag) ;
/*
   -----------------------
   fill the column indices
   -----------------------
*/
Chv_columnIndices(frontJ, &ncolJ, &ivec) ;
IVcopy(ncolJ, ivec, colindJ) ;
#if MYDEBUG > 0
fprintf(stdout, "\n front %d, colind = %p", J, ivec) ;
IVfprintf(stdout, ncolJ, ivec) ;
fflush(stdout) ;
#endif
if ( FRONTMTX_IS_NONSYMMETRIC(frontmtx) ) {
/*
   --------------------
   fill the row indices
   --------------------
*/
   Chv_rowIndices(frontJ, &nrowJ, &ivec) ;
   IVcopy(nrowJ, ivec, colindJ) ;
}
/*
   ------------------------
   zero the front's entries
   ------------------------
*/
Chv_zero(frontJ) ;
#if MYDEBUG > 0
fprintf(stdout, "\n leaving Front_initializeFront()") ;
Chv_writeForHumanEye(frontJ, stdout) ;
fflush(stdout) ;
#endif

return ; }

/*--------------------------------------------------------------------*/
/*
   ----------------
   static functions
   ----------------
*/
static void assembleAggregates ( Chv *frontJ, ChvList *aggList,
   ChvManager *chvmanager, double cpus[], int msglvl, FILE *msgFile ) ;
static char assemblePostponedData ( FrontMtx *frontmtx, Chv *frontJ,
   int *pndelay, Chv *fronts[], ChvList *postList,
   ChvManager *chvmanager, double cpus[], int msglvl, FILE *msgFile ) ;
static int factorFront ( FrontMtx *frontmtx, Chv *frontJ,
   int ndelay, double tau, IV *pivotsizesIV, DV *workDV, int stats[],
   double cpus[], int msglvl, FILE *msgFile ) ;
static void storeEntries ( FrontMtx *frontmtx, Chv *frontJ, int nelim,
   double droptol, IV *pivotsizesIV, ChvList *postList,
   ChvManager *chvmanager, int parent[], int stats[], double cpus[],
   int msglvl, FILE *msgFile ) ;
/*
   ------------------------------------------------------------------
   purpose -- to visit a front during a factorization.
              note: this method is called by the serial,
              multithreaded and MPI factorization codes.

   frontmtx     -- front matrix object
   pencil       -- matrix pencil for A + sigma*B
   J            -- id of front we are working on
   myid         -- id of thread or process
   owners[]     -- map from fronts to owning threads,
                   in a serial environment, owners = NULL
   fronts[]     -- vector of pointers to fronts
   lookahead    -- parameter controls the partial upward visiting
                   of ancestors before returning
   tau          -- used when pivoting enabled,
                   |L_{j,i}| and |U_{i,j}| <= tau
   droptol      -- used for an approximate factorization
                   stored |L_{j,i}| and |U_{i,j}| > droptol
   status[]     -- status vector for the fronts,
                   'W' -- needs to be woke up
                   'A' -- active front
                   'F' -- front is finished
   heads[]      -- vector of pointers to IP objects that store the
                   front-to-front update lists
   pivotsizesIV -- IV object used during the factorization of a front
                   when pivoting is enabled
   workDV       -- DV object used for working storage
   parent       -- front parent vector
   aggList      -- list object used in parallel environment, used
                   to store aggregate fronts
   postList     -- list object used in pivoting and/or parallel
                   environments, used to stored delayed data and/or
                   to synchronize the factorization
   chvmanager   -- used to manage working storage for Chv objects
   stats[]      -- statistics vector
   cpus[]       -- vector to hold breakdown of cpu times
   msglvl       -- message level
   msgFil       -- message file

   created -- 98may04, cca
   ------------------------------------------------------------------
*/
char
FrontMtx_factorVisit (
   FrontMtx     *frontmtx,
   Pencil       *pencil,
   int          J,
   int          myid,
   int          owners[],
   Chv          *fronts[],
   int          lookahead,
   double       tau,
   double       droptol,
   char         status[],
   IP           *heads[],
   IV           *pivotsizesIV,
   DV           *workDV,
   int          parent[],
   ChvList      *aggList,
   ChvList      *postList,
   ChvManager   *chvmanager,
   int          stats[],
   double       cpus[],
   int          msglvl,
   FILE         *msgFile
) {
/*
   ---------------
   local variables
   ---------------
*/
char     allAggregatesHere, allPostponedAssmb, allUpdatesDone ;
Chv     *frontJ ;
double   t1, t2 ;
int      K, ndelay, nelim ;

if ( status[J] == 'F' ) {
   return('F') ;
}
allUpdatesDone    = 'N' ;
allAggregatesHere = 'N' ;
allPostponedAssmb = 'N' ;
frontJ = NULL ;
if ( heads[J] != NULL ) {
/*
   --------------------------------
   internal updates to be performed
   --------------------------------
*/
   if ( (frontJ = fronts[J]) == NULL ) {
/*
      --------------------
      initialize the front
      --------------------
*/
      fronts[J] = FrontMtx_setupFront(frontmtx, pencil, J, myid,
                      owners, chvmanager, cpus, msglvl, msgFile) ;
      frontJ = fronts[J] ;
      status[J] = 'A' ;
   }
/*
   ---------------------------------
   compute updates from owned fronts
   ---------------------------------
*/
   if ( msglvl > 1 ) {
      fprintf(msgFile, "\n performing updates") ;
      fflush(msgFile) ;
   }
   MARKTIME(t1) ;
   FrontMtx_update(frontmtx, frontJ, heads,
                    status, workDV, msglvl, msgFile) ;
   MARKTIME(t2) ;
   cpus[2] += t2 - t1 ;
}
if ( heads[J] == NULL ) {
   allUpdatesDone = 'Y' ;
}
if ( owners == NULL || owners[J] == myid ) {
/*
   --------------
   front is owned
   --------------
*/
   if ( (frontJ = fronts[J]) == NULL ) {
/*
      --------------------
      initialize the front
      --------------------
*/
      fronts[J] = FrontMtx_setupFront(frontmtx, pencil, J, myid,
                            owners, chvmanager, cpus, msglvl, msgFile) ;
      frontJ = fronts[J] ;
      status[J] = 'A' ;
   }
   if ( aggList != NULL ) {
/*
      ------------------------------------------
      we are operating in a parallel environment
      ------------------------------------------
*/
      if ( ChvList_isListNonempty(aggList, J) == 1 ) {
/*
         -------------------------------------------------
         assemble any aggregate updates from other threads
         -------------------------------------------------
*/
         assembleAggregates(frontJ, aggList, chvmanager,
                            cpus, msglvl, msgFile) ;
      }
      if ( ChvList_isCountZero(aggList, J) == 1 ) {
/*
         -------------------------------------------------------
         all aggregates are assembled or in the list to assemble
         -------------------------------------------------------
*/
         if ( ChvList_isListNonempty(aggList, J) == 1 ) {
/*
            -------------------------------------------------
            assemble any aggregate updates from other threads
            -------------------------------------------------
*/
            assembleAggregates(frontJ, aggList, chvmanager,
                               cpus, msglvl, msgFile) ;
         }
         allAggregatesHere = 'Y' ;
      }
   } else {
      allAggregatesHere = 'Y' ;
   }
   if ( msglvl > 1 ) {
      fprintf(msgFile, "\n\n allUpdatesDone %c, allAggregatesHere %c",
              allUpdatesDone, allAggregatesHere) ;
      fflush(msgFile) ;
   }
   if ( allUpdatesDone == 'Y' && allAggregatesHere == 'Y' ) {
/*
      -------------------------------------
      all internal updates and all external
      aggregates have been incorporated
      -------------------------------------
*/
      if ( postList != NULL ) {
/*
         ---------------------------------------------
         front is ready to assemble any postponed data
         ---------------------------------------------
*/
         allPostponedAssmb = assemblePostponedData(frontmtx, frontJ,
                                    &ndelay, fronts, postList,
                                    chvmanager, cpus, msglvl, msgFile) ;
         frontJ = fronts[J] ;
      } else {
         allPostponedAssmb = 'Y' ;
         ndelay = 0 ;
      }
      if ( msglvl > 1 ) {
         fprintf(msgFile,
                 "\n\n allPostponedAssmb %c", allPostponedAssmb) ;
         fflush(msgFile) ;
      }
      if ( allPostponedAssmb == 'Y' ) {
/*
         -----------------------------
         front is ready to be factored
         -----------------------------
*/
         nelim = factorFront(frontmtx, frontJ, ndelay, tau,
                             pivotsizesIV, workDV, stats,
                             cpus, msglvl, msgFile) ;
         if ( msglvl > 1 ) {
            fprintf(msgFile,
                    "\n\n J = %d, nelim = %d", frontJ->id, nelim) ;
            fflush(msgFile) ;
         }
         if ( (! FRONTMTX_IS_PIVOTING(frontmtx)) && nelim < frontJ->nD){
/*
            ------------------------------------------------
            a singular front matrix has been found.
            release the front and set the status to error
            ------------------------------------------------
*/
            ChvManager_releaseObject(chvmanager, frontJ) ;
            fronts[J] = NULL ;
            status[J] = 'E' ;
         } else {
/*
            -------------------------------------------
            store any factor entries and postponed data
            -------------------------------------------
*/
            storeEntries(frontmtx, frontJ, nelim, droptol, pivotsizesIV,
                         postList, chvmanager, parent, stats,
                         cpus, msglvl, msgFile) ;
/*
            ------------------------------------------------
            release the front and set the status to finished
            ------------------------------------------------
*/
            ChvManager_releaseObject(chvmanager, frontJ) ;
            fronts[J] = NULL ;
            status[J] = 'F' ;
         }
      }
   }
} else if ( allUpdatesDone == 'Y' ) {
   if ( frontJ != NULL ) {
/*
      --------------------------------------------------------
      front is not owned and all internal updates are complete
      put aggregate update front on the aggregate list
      --------------------------------------------------------
*/
      if ( msglvl > 1 ) {
         fprintf(msgFile, "\n done with unowned front %3d", J) ;
         fflush(msgFile) ;
      }
      if ( msglvl > 3 ) {
         Chv_writeForHumanEye(frontJ, msgFile) ;
         fflush(msgFile) ;
      }
      MARKTIME(t1) ;
      ChvList_addObjectToList(aggList, frontJ, J) ;
      MARKTIME(t2) ;
#if MYDEBUG > 0
      fprintf(stdout, "\n thread %2d, done with unowned front %3d",
              myid, J) ;
      fflush(stdout) ;
#endif
      cpus[7] += t2 - t1 ;
   }
   status[J] = 'F' ;
}
if ( --lookahead >= 0 && (K = parent[J]) != -1 ) {
/*
   -----------------------------
   visit parent before returning
   -----------------------------
*/
   FrontMtx_factorVisit(frontmtx, pencil, J, myid, owners, fronts,
                   lookahead, tau, droptol, status, heads, pivotsizesIV,
                   workDV, parent, aggList, postList, chvmanager,
                   stats, cpus, msglvl, msgFile) ;
}
return(status[J]) ; }

/*--------------------------------------------------------------------*/
/*
   --------------------------------------------
   purpose -- set up the front's data structure

   created -- 98mar27, cca
   --------------------------------------------
*/
Chv *
FrontMtx_setupFront (
   FrontMtx     *frontmtx,
   Pencil       *pencil,
   int          J,
   int          myid,
   int          owners[],
   ChvManager   *chvmanager,
   double       cpus[],
   int          msglvl,
   FILE         *msgFile
) {
Chv     *frontJ ;
double   t1, t2 ;
int      nbytes, nD, nL, nU ;

if ( msglvl > 4 ) {
   fprintf(msgFile,
           "\n\n inside FrontMtx_setupFront()"
           "\n frontmtx %p, pencil %p, J %d, myid %d"
           "\n owners %p, chvmanager %p, cpus %p"
           "\n msglvl %d, msgFile %p",
           frontmtx, pencil, J, myid, owners, chvmanager,
           cpus, msglvl, msgFile) ;
   fflush(msgFile) ;
}
MARKTIME(t1) ;
FrontMtx_initialFrontDimensions(frontmtx, J,
                                 &nD, &nL, &nU, &nbytes) ;
if ( msglvl > 2 ) {
   fprintf(msgFile, "\n nD %d, nL %d, nU %d, nbytes %d",
           nD, nL, nU, nbytes) ;
   fflush(msgFile) ;
}
frontJ = ChvManager_newObjectOfSizeNbytes(chvmanager, nbytes) ;
if ( msglvl > 2 ) {
   fprintf(msgFile, "\n frontJ = %p", frontJ) ;
   fflush(msgFile) ;
}
Chv_init(frontJ, J, nD, nL, nU, frontmtx->type, frontmtx->symmetryflag);
FrontMtx_initializeFront(frontmtx, frontJ, J) ;
MARKTIME(t2) ;
cpus[0] += t2 - t1 ;
if ( pencil != NULL && (owners == NULL || owners[J] == myid) ) {
/*
   -------------------------------------------------
   thread owns this front, load the original entries
   -------------------------------------------------
*/
   MARKTIME(t1) ;
   FrontMtx_loadEntries(frontJ, pencil, msglvl, msgFile) ;
   MARKTIME(t2) ;
   cpus[1] += t2 - t1 ;
   if ( msglvl > 1 ) {
      fprintf(msgFile, "\n original entries loaded") ;
      fflush(msgFile) ;
   }
/*
   if ( FRONTMTX_IS_HERMITIAN(frontmtx) && J == frontmtx->nfront - 1 ) {
      fprintf(msgFile, "\n last front after entries loaded") ;
      Chv_writeForHumanEye(frontJ, msgFile) ;
      fflush(msgFile) ;
   }
*/
}
if ( msglvl > 1 ) {
   fprintf(msgFile, "\n\n front initialized") ;
   fflush(msgFile) ;
}
if ( msglvl > 3 ) {
   Chv_writeForHumanEye(frontJ, msgFile) ;
   fflush(msgFile) ;
}
return(frontJ) ; }

/*--------------------------------------------------------------------*/
/*
   ------------------------------------------------------------
   purpose -- assemble any aggregate updates from other threads
   ------------------------------------------------------------
*/
static void
assembleAggregates (
   Chv          *frontJ,
   ChvList      *aggList,
   ChvManager   *chvmanager,
   double       cpus[],
   int          msglvl,
   FILE         *msgFile
) {
Chv     *chv, *headchv ;
double   t1, t2 ;

MARKTIME(t1) ;
headchv = ChvList_getList(aggList, frontJ->id) ;
for ( chv = headchv ; chv != NULL ; chv = chv->next ) {
   Chv_assembleChv(frontJ, chv) ;
}
MARKTIME(t2) ;
cpus[8] += t2 - t1 ;
ChvManager_releaseListOfObjects(chvmanager, headchv) ;
if ( msglvl > 3 ) {
   fprintf(msgFile, "\n after assembly") ;
   Chv_writeForHumanEye(frontJ, msgFile) ;
   fflush(msgFile) ;
}
return ; }

/*--------------------------------------------------------------------*/
/*
   --------------------------------------
   purpose -- assemble any postponed data

   created -- 98mar27, cca
   --------------------------------------
*/
static char
assemblePostponedData (
   FrontMtx     *frontmtx,
   Chv          *frontJ,
   int          *pndelay,
   Chv          *fronts[],
   ChvList      *postList,
   ChvManager   *chvmanager,
   double       cpus[],
   int          msglvl,
   FILE         *msgFile
) {
char     rc ;
double   t1, t2 ;
int      J ;

if ( msglvl > 4 ) {
   fprintf(msgFile, "\n\n frontmtx %p, frontJ %p, pndelay %p"
           "\n fronts %p, postList %p, chvmanager %p, cpus %p",
           frontmtx, frontJ, pndelay,
           fronts, postList, chvmanager, cpus) ;
   fflush(msgFile) ;
}


J = frontJ->id ;
if ( msglvl > 1 ) {
   if ( ChvList_isCountZero(postList, J) == 1) {
      fprintf(msgFile, "\n all postponed data is here") ;
      fflush(msgFile) ;
   } else {
      fprintf(msgFile, "\n still waiting on postponed data") ;
      fflush(msgFile) ;
   }
}
if ( ChvList_isCountZero(postList, J) == 1 ) {
   if ( msglvl > 1 ) {
      fprintf(msgFile, "\n assembling postponed data") ;
      fflush(msgFile) ;
   }
/*
   ---------------------------
   assemble any postponed data
   ---------------------------
*/
   MARKTIME(t1) ;
   fronts[J] = FrontMtx_assemblePostponedData(frontmtx, frontJ,
                                   postList, chvmanager, pndelay) ;
   if ( frontJ != fronts[J] ) {
      if ( msglvl > 1 ) {
         fprintf(msgFile, "\n releasing old front") ;
         fflush(msgFile) ;
      }
      ChvManager_releaseObject(chvmanager, frontJ) ;
   }
   MARKTIME(t2) ;
   cpus[3] += t2 - t1 ;
   rc = 'Y' ;
} else {
   rc = 'N' ;
}
return(rc) ; }

/*--------------------------------------------------------------------*/
/*
   ---------------------------
   purpose -- factor the front

   created -- 98mar27, cca
   ---------------------------
*/
static int
factorFront (
   FrontMtx   *frontmtx,
   Chv        *frontJ,
   int        ndelay,
   double     tau,
   IV         *pivotsizesIV,
   DV         *workDV,
   int        stats[],
   double     cpus[],
   int        msglvl,
   FILE       *msgFile
) {
double   t1, t2 ;
int      nelim, npost ;

if ( msglvl > 2 ) {
   fprintf(msgFile, "\n frontJ = %p, ndelay = %d",
           frontJ, ndelay) ;
   fprintf(msgFile, "\n tau = %12.4e", tau) ;
   fprintf(msgFile, "\n stats %p, cpus %p", stats, cpus) ;
   fflush(msgFile) ;
}
if ( msglvl > 2 ) {
   Chv_writeForHumanEye(frontJ, msgFile) ;
   fflush(msgFile) ;
}
MARKTIME(t1) ;
if ( FRONTMTX_IS_PIVOTING(frontmtx) ) {
   nelim = Chv_factorWithPivoting(frontJ, ndelay,
                                  frontmtx->pivotingflag, pivotsizesIV,
                                  workDV, tau, &stats[1]) ;
} else {
   nelim = Chv_factorWithNoPivoting(frontJ, frontmtx->patchinfo) ;
}
npost = frontJ->nD - nelim ;
stats[2] += npost ;
if (  !(FRONTMTX_IS_PIVOTING(frontmtx))
   || FRONTMTX_IS_NONSYMMETRIC(frontmtx) ) {
   stats[0] += nelim ;
} else {
   stats[0] += IV_size(pivotsizesIV) ;
}
MARKTIME(t2) ;
cpus[4] += t2 - t1 ;
if ( msglvl > 1 ) {
   fprintf(msgFile, "\n\n front %d, nelim = %d, npost = %d",
           frontJ->id, nelim, npost) ;
   fflush(msgFile) ;
}
if ( msglvl > 2 ) {
   Chv_writeForHumanEye(frontJ, msgFile) ;
   fflush(msgFile) ;
}
return(nelim) ; }

/*--------------------------------------------------------------------*/
/*
   ---------------------------------------------------------------
   purpose -- store the factor entries in the front matrix object,
              and extract any postponed data and put it in the
              list object

   created -- 98mar27, cca
   ---------------------------------------------------------------
*/
static void
storeEntries (
   FrontMtx     *frontmtx,
   Chv          *frontJ,
   int          nelim,
   double       droptol,
   IV           *pivotsizesIV,
   ChvList      *postList,
   ChvManager   *chvmanager,
   int          parent[],
   int          stats[],
   double       cpus[],
   int          msglvl,
   FILE         *msgFile
) {
double   t1, t2 ;
int      npost ;

npost = frontJ->nD - nelim ;
if ( msglvl > 1 ) {
   fprintf(msgFile, "\n storing factor data, nelim = %d", nelim) ;
   fflush(msgFile) ;
}
MARKTIME(t1) ;
frontJ->nD -= npost ;
frontJ->nL += npost ;
frontJ->nU += npost ;
FrontMtx_storeFront(frontmtx, frontJ, pivotsizesIV, droptol,
                     msglvl, msgFile) ;
frontJ->nD += npost ;
frontJ->nL -= npost ;
frontJ->nU -= npost ;
MARKTIME(t2) ;
cpus[6] += t2 - t1 ;
if ( postList != NULL ) {
   Chv   *postponedchv ;
   if ( npost > 0 ) {
/*
      ---------------------------------
      there is postponed data,
      extract and put on postponed list
      ---------------------------------
*/
      postponedchv = frontJ ;
      if ( msglvl > 2 ) {
         fprintf(msgFile, "\n postponed data for front %d",
                 frontJ->id) ;
         Chv_writeForHumanEye(postponedchv, msgFile) ;
         fflush(msgFile) ;
      }
   } else {
      postponedchv = NULL ;
   }
   if ( msglvl > 1 ) {
      fprintf(msgFile, "\n storing postponed data %p", postponedchv);
      fflush(msgFile) ;
   }
   MARKTIME(t1) ;
   FrontMtx_storePostponedData(frontmtx, postponedchv, npost,
                       parent[frontJ->id], postList, chvmanager) ;
   MARKTIME(t2) ;
   cpus[5] += t2 - t1 ;
}
return ; }

/*--------------------------------------------------------------------*/
/*
   --------------------------------------------------------------------
   purpose -- to set up the link data structures
      for a parallel factorization for process myid

   return value -- pointer to IP *heads[nfront+2], which contains
      the beginning of a list of IP objects that store the remaining
      updates to the fronts.
      note, heads[nfront] is the first IP object in the free list.
      heads[nfront+1] is the base address of the allocated IP objects.

   created -- 98mar07, cca
   --------------------------------------------------------------------
*/
IP **
FrontMtx_factorSetup (
   FrontMtx   *frontmtx,
   IV         *frontOwnersIV,
   int        myid,
   int        msglvl,
   FILE       *msgFile
) {
int   count, ii, J, K, nadj, nfront ;
int   *adj, *mark, *owners, *vtxToFront ;
IP    *ip ;
IP    **heads ;
IVL   *symbfacIVL ;
/*
   ---------------------------------------------------
   count the number of updates this front will perform
   ---------------------------------------------------
*/
nfront = FrontMtx_nfront(frontmtx) ;
if ( frontOwnersIV != NULL ) {
   owners = IV_entries(frontOwnersIV) ;
} else {
   owners = NULL ;
}
symbfacIVL = frontmtx->symbfacIVL ;
vtxToFront = ETree_vtxToFront(frontmtx->frontETree) ;
mark = IVinit(nfront, -1) ;
for ( J = count = 0 ; J < nfront ; J++ ) {
   if ( owners == NULL || owners[J] == myid ) {
      IVL_listAndSize(symbfacIVL, J, &nadj, &adj) ;
      mark[J] = J ;
      for ( ii = 0 ; ii < nadj ; ii++ ) {
         K = vtxToFront[adj[ii]] ;
         if ( mark[K] != J ) {
            mark[K] = J ;
            count++ ;
         }
      }
   }
}
/*
   --------------------------------------------------
   set up the update head/links for the factorization
   --------------------------------------------------
*/
ALLOCATE(heads, struct _IP *, nfront + 2) ;
for ( J = 0 ; J <= nfront + 1 ; J++ ) {
   heads[J] = NULL ;
}
heads[nfront] = heads[nfront+1] = IP_init(count, 1) ;
IVfill(nfront, mark, -1) ;
for ( J = 0 ; J < nfront ; J++ ) {
   if ( owners == NULL || owners[J] == myid ) {
      IVL_listAndSize(symbfacIVL, J, &nadj, &adj) ;
      mark[J] = J ;
      for ( ii = 0 ; ii < nadj ; ii++ ) {
         K = vtxToFront[adj[ii]] ;
         if ( mark[K] != J ) {
            mark[K] = J ;
            ip = heads[nfront] ; heads[nfront] = ip->next ;
            ip->val = J ; ip->next = heads[K] ; heads[K] = ip ;
            if ( msglvl > 3 ) {
               fprintf(msgFile, "\n linking L(%d,%d) to L(%d,%d)",
                     K, J, K, (ip->next == NULL) ? -1 : ip->next->val) ;
               fflush(msgFile) ;
            }
         }
      }
   }
}
IVfree(mark) ;

return(heads) ; }

/*--------------------------------------------------------------------*/
/*
   -----------------------------------------------
   create and return the nactiveChild vector.
   nactiveChild[J] contains the number of children
   of J that belong to an active path

   created -- 97jul03, cca
   -----------------------------------------------
*/
int *
FrontMtx_nactiveChild (
   FrontMtx   *frontmtx,
   char       *status,
   int        myid
) {
int    J, K, nfront ;
int    *nactiveChild, *par ;
/*
   ---------------
   check the input
   ---------------
*/
if ( frontmtx == NULL || status == NULL || myid < 0 ) {
   fprintf(stderr, "\n fatal error in FrontMtx_nativeChild(%p,%p,%d)"
           "\n bad input\n", frontmtx, status, myid) ;
   exit(-1) ;
}
nfront = frontmtx->nfront ;
par    = ETree_par(frontmtx->frontETree) ;
/*
   ---------------------------------------
   create and fill the nactiveChild vector
   ---------------------------------------
*/
nactiveChild = IVinit(nfront, 0) ;
for ( J = 0 ; J < nfront ; J++ ) {
   if ( status[J] == 'W' && (K = par[J]) != -1 ) {
      nactiveChild[K]++ ;
   }
}
return(nactiveChild) ; }

/*--------------------------------------------------------------------*/
/*
   ---------------------------------------------------------
   create, initialize and return a Ideq object
   that will be used for a parallel factorization,
   forward solve, or backward solve.

   the status[] vector will be set as follows:
      status[J] = activeflag   if J is on an active path
      status[J] = inactiveflag if J is not on an active path

   created -- 98mar27, cca
   ---------------------------------------------------------
*/
Ideq *
FrontMtx_setUpDequeue (
   FrontMtx   *frontmtx,
   int        owners[],
   int        myid,
   char       status[],
   IP         *heads[],
   char       activeFlag,
   char       inactiveFlag,
   int        msglvl,
   FILE       *msgFile
) {
Ideq   *dequeue ;
int    J, K, nfront, npath ;
int    *par ;
Tree   *tree ;
/*
   ---------------
   check the input
   ---------------
*/
if (  frontmtx == NULL || owners == NULL
   || status == NULL || myid < 0 ) {
   fprintf(stderr,
           "\n fatal error in FrontMtx_setUpDequeue()"
           "\n frontmtx %p, owners %p, myid %d, status %p"
           "\n heads %p, activeFlag %c, inactiveFlag %c"
           "\n bad input\n",
           frontmtx, owners, myid, status, heads,
           activeFlag, inactiveFlag) ;
   exit(-1) ;
}
nfront = frontmtx->nfront ;
tree   = frontmtx->tree ;
par    = tree->par ;
/*
   ------------------------------------
   count the number of active paths in
   the tree and set the status[] vector
   ------------------------------------
*/
CVfill(nfront, status, inactiveFlag) ;
for ( J = npath = 0 ; J < nfront ; J++ ) {
   if ( status[J] == inactiveFlag ) {
      if ( owners[J] == myid || (heads != NULL && heads[J] != NULL) ) {
         npath++ ;
         for ( K = J ;
               K != -1 && status[K] != activeFlag ;
               K = par[K] ) {
            status[K] = activeFlag ;
         }
      }
   }
}
dequeue = Ideq_new() ;
Ideq_resize(dequeue, npath) ;

return(dequeue) ; }

/*--------------------------------------------------------------------*/
/*
   -----------------------------------------------------------------
   purpose -- load the dequeue with the leaves of the active subtree
              used for a factorization and forward solve

   created -- 98mar27, cca
   -----------------------------------------------------------------
*/
void
FrontMtx_loadActiveLeaves (
   FrontMtx   *frontmtx,
   char       status[],
   char       activeFlag,
   Ideq       *dequeue
) {
int   I, J, nactivechild, nfront ;
int   *fch, *sib ;

nfront = frontmtx->nfront ;
fch    = frontmtx->tree->fch ;
sib    = frontmtx->tree->sib ;
for ( J = 0 ; J < nfront ; J++ ) {
   if ( status[J] == activeFlag ) {
      if ( fch[J] == -1 ) {
         Ideq_insertAtTail(dequeue, J) ;
      } else {
         nactivechild = 0 ;
         for ( I = fch[J] ; I != -1 ; I = sib[I] ) {
            if ( status[I] == activeFlag ) {
               nactivechild++ ;
               break ;
            }
         }
         if ( nactivechild == 0 ) {
            Ideq_insertAtTail(dequeue, J) ;
         }
      }
   }
}
return ; }

/*--------------------------------------------------------------------*/
/*
   -----------------------------------------------
   create, initialize and return a ChvList object
   to deal with postponed chevrons

   created -- 97jul03, cca
   -----------------------------------------------
*/
ChvList *
FrontMtx_postList (
   FrontMtx   *frontmtx,
   IV         *frontOwnersIV,
   int        lockflag
) {
char      *flags ;
ChvList   *postList ;
int       count, I, J, jthread, nchild, nfront, nthread ;
int       *counts, *fch, *frontOwners, *mark, *sib ;
/*
   ---------------
   check the input
   ---------------
*/
if (  frontmtx == NULL || frontOwnersIV == NULL
   || lockflag < 0 || lockflag > 2 ) {
   fprintf(stderr,
           "\n fatal error in FrontMtx_postList(%p,%p,%d)"
           "\n bad input\n", frontmtx, frontOwnersIV, lockflag) ;
   exit(-1) ;
}
fch = ETree_fch(frontmtx->frontETree) ;
sib = ETree_sib(frontmtx->frontETree) ;
IV_sizeAndEntries(frontOwnersIV, &nfront, &frontOwners) ;
counts = IVinit(nfront+1, 0) ;
if ( lockflag > 0 ) {
   flags = CVinit(nfront+1, 'N') ;
} else {
   flags = NULL ;
}
nthread = 1 + IV_max(frontOwnersIV) ;
mark    = IVinit(nthread, -1) ;
/*
   --------------------
   loop over the fronts
   --------------------
*/
for ( J = 0 ; J < nfront ; J++ ) {
   count = nchild = 0 ;
   for ( I = fch[J] ; I != -1 ; I = sib[I] ) {
      nchild++ ;
      jthread = frontOwners[I] ;
      if ( mark[jthread] != J ) {
         mark[jthread] = J ;
         count++ ;
      }
   }
   counts[J] = nchild ;
   if ( flags != NULL ) {
      if ( count > 1 ) {
         flags[J] = 'Y' ;
      } else {
         flags[J] = 'N' ;
      }
   }
}
count = nchild = 0 ;
for ( J = ETree_root(frontmtx->frontETree) ; J != -1 ; J = sib[J] ) {
   nchild++ ;
   jthread = frontOwners[J] ;
   if ( mark[jthread] != J ) {
      mark[jthread] = J ;
      count++ ;
   }
}
counts[nfront] = nchild ;
if ( flags != NULL ) {
   if ( count > 1 ) {
      flags[nfront] = 'Y' ;
   } else {
      flags[nfront] = 'N' ;
   }
}
/*
   -----------------------------------------
   create and initialize the ChvList object
   -----------------------------------------
*/
postList = ChvList_new() ;
ChvList_init(postList, nfront+1, counts, lockflag, flags) ;
/*
   ------------------------
   free the working storage
   ------------------------
*/
IVfree(mark) ;
IVfree(counts) ;
if ( flags != NULL ) {
   CVfree(flags) ;
}

return(postList) ; }

/*--------------------------------------------------------------------*/
/*
   -----------------------------------------------
   create, initialize and return a ChvList object
   to deal with aggregate chevrons

   created -- 97jul03, cca
   -----------------------------------------------
*/
ChvList *
FrontMtx_aggregateList (
   FrontMtx   *frontmtx,
   IV         *frontOwnersIV,
   int        lockflag
) {
char       *flags ;
ChvList   *aggList ;
int        count, ii, I, J, jthread, K, myid, nfront, nthread, size ;
int        *counts, *frontOwners, *head, *indices, *link,
           *mark, *offsets, *vtxToFront ;
IVL        *symbfacIVL ;

/*
   ---------------
   check the input
   ---------------
*/
if (  frontmtx == NULL || frontOwnersIV == NULL
   || lockflag < 0 || lockflag > 2 ) {
   fprintf(stderr,
           "\n fatal error in FrontMtx_aggregateList(%p,%p,%d)"
           "\n bad input\n", frontmtx, frontOwnersIV, lockflag) ;
   exit(-1) ;
}
symbfacIVL = frontmtx->symbfacIVL ;
vtxToFront = ETree_vtxToFront(frontmtx->frontETree) ;
IV_sizeAndEntries(frontOwnersIV, &nfront, &frontOwners) ;
nthread = 1 + IV_max(frontOwnersIV) ;
mark    = IVinit(nthread, -1) ;
head    = IVinit(nfront,  -1) ;
link    = IVinit(nfront,  -1) ;
offsets = IVinit(nfront,   0) ;
counts  = IVinit(nfront,   0) ;
if ( lockflag > 0 ) {
   flags = CVinit(nfront, 'N') ;
} else {
   flags = NULL ;
}
/*
   --------------------
   loop over the fronts
   --------------------
*/
for ( J = 0 ; J < nfront ; J++ ) {
   myid = frontOwners[J] ;
#if MYDEBUG > 0
   fprintf(stdout, "\n\n front %d, owner %d", J, myid) ;
   fflush(stdout) ;
#endif
   mark[myid] = J ;
   count = 0 ;
/*
   ---------------------------------------------------
   loop over all descendent fronts that might update J
   ---------------------------------------------------
*/
   while ( (I = head[J]) != -1 ) {
      head[J] = link[I] ;
      jthread = frontOwners[I] ;
#if MYDEBUG > 0
      fprintf(stdout, "\n descendent front %d, owner %d", I,
jthread) ;
      fflush(stdout) ;
#endif
      if ( mark[jthread] != J ) {
/*
         --------------------------------
         expect an aggregate from jthread
         --------------------------------
*/
#if MYDEBUG > 0
         fprintf(stdout, ", incoming aggregate") ;
         fflush(stdout) ;
#endif
         mark[jthread] = J ;
         count++ ;
      }
/*
      --------------------------------------------------
      link front I to next ancestor that it does not own
      --------------------------------------------------
*/
      IVL_listAndSize(symbfacIVL, I, &size, &indices) ;
      for ( ii = offsets[I] ; ii < size ; ii++ ) {
         if (  (K = vtxToFront[indices[ii]]) > J
            && frontOwners[K] != jthread ) {
#if MYDEBUG > 0
            fprintf(stdout, ", link to %d", K) ;
            fflush(stdout) ;
#endif
            offsets[I] = ii ;
            link[I] = head[K] ;
            head[K] = I ;
            break ;
         }
      }
   }
/*
   -------------------------------------
   set the number of incoming aggregates
   -------------------------------------
*/
   counts[J] = count ;
#if MYDEBUG > 0
   fprintf(stdout, "\n counts[%d] = %d", J, counts[J]) ;
   fflush(stdout) ;
#endif
/*
   ---------------------------------------------------
   set the flags to see if the list needs to be locked
   ---------------------------------------------------
*/
   if ( flags != NULL ) {
      if ( count > 1 ) {
         flags[J] = 'Y' ;
      } else {
         flags[J] = 'N' ;
      }
#if MYDEBUG > 0
      fprintf(stdout, ", flags[%d] = %c", J, flags[J]) ;
      fflush(stdout) ;
#endif
   }
/*
   --------------------------------------------------
   link front J to next ancestor that it does not own
   --------------------------------------------------
*/
   IVL_listAndSize(symbfacIVL, J, &size, &indices) ;
   for ( ii = 0 ; ii < size ; ii++ ) {
      if (  (K = vtxToFront[indices[ii]]) > J && frontOwners[K] != myid
) {
#if MYDEBUG > 0
         fprintf(stdout, ", linking to %d", K) ;
         fflush(stdout) ;
#endif
         offsets[J] = ii ;
         link[J] = head[K] ;
         head[K] = J ;
         break ;
      }
   }
}
#if MYDEBUG > 0
fprintf(stdout, "\n counts") ;
IVfprintf(stdout, nfront, counts) ;
fflush(stdout) ;
#endif
/*
   -----------------------------------------
   create and initialize the ChvList object
   -----------------------------------------
*/
aggList = ChvList_new() ;
ChvList_init(aggList, nfront, counts, lockflag, flags) ;
/*
   ------------------------
   free the working storage
   ------------------------
*/
IVfree(counts) ;
IVfree(head) ;
IVfree(link) ;
IVfree(offsets) ;
IVfree(mark) ;
if ( flags != NULL ) {
   CVfree(flags) ;
}
return(aggList) ; }

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