/*
   ----------------------------------------
   if the new-to-old vector is not present, 
   create it and fill its entries

   created -- 96mar16, cca
   ----------------------------------------
*/
void
Perm_fillNewToOld (
   Perm   *perm
) {
int   size ;
/*
   ---------------
   check the input
   ---------------
*/
if (  perm == NULL 
   || perm->isPresent < 1 || perm->isPresent > 3
   || (size = perm->size) <= 0 ) {
   fprintf(stderr, "\n fatal error in Perm_fillNewToOld(%p)"
           "\n bad input\n", perm) ;
   spoolesFatal();
}
if ( perm->isPresent == 2 ) {
   int   iold ;
   int   *newToOld = perm->newToOld = IVinit(size, -1) ;
   int   *oldToNew = perm->oldToNew ;

   for ( iold = 0 ; iold < size ; iold++ ) {
      newToOld[oldToNew[iold]] = iold ;
   }
   perm->isPresent = 3 ;
}
return ; }

/*
   ----------------------------------------
   if the old-to-new vector is not present, 
   create it and fill its entries

   created -- 96mar16, cca
   ----------------------------------------
*/
void
Perm_fillOldToNew (
   Perm   *perm
) {
int   size ;
/*
   ---------------
   check the input
   ---------------
*/
if (  perm == NULL 
   || perm->isPresent < 1 || perm->isPresent > 3
   || (size = perm->size) <= 0 ) {
   fprintf(stderr, "\n fatal error in Perm_fillOldToNew(%p)"
           "\n bad input\n", perm) ;
   spoolesFatal();
}
if ( perm->isPresent == 1 ) {
   int   inew ;
   int   *newToOld = perm->newToOld ;
   int   *oldToNew = perm->oldToNew = IVinit(size, -1) ;

   for ( inew = 0 ; inew < size ; inew++ ) {
      oldToNew[newToOld[inew]] = inew ;
   }
   perm->isPresent = 3 ;
}
return ; }

/*
   -------------------------------------------------------------
   purpose -- allocate an int array and fill with the inverse of
              y[]. note, y[] must be a permutation vector.

   created : 95sep22, cca
   -------------------------------------------------------------
*/
int *
IVinverse ( 
   int   size, 
   int   y[] 
) {
int   *x = NULL ;
if ( size > 0 ) {
   if ( y == NULL ) {
      fprintf(stderr, "\n fatal error in IVinverse, invalid data"
              "\n size = %d, y = %p\n", size, y) ;
      exit(-1) ;
   } else {
      int   i, j ;
      x = IVinit(size, -1) ;
      for ( i = 0 ; i < size ; i++ ) {
         j = y[i] ;
         if ( j < 0 || j >= size || x[j] != -1 ) {
            fprintf(stderr, 
                    "\n fatal error in IVinverse"
                    "\n y[%d] = %d, value out-of-range or repeated",
                    i, j) ;
            exit(-1) ;
         }
         x[j] = i ; 
      }
   }
}
return(x) ; }

/*
   ---------------------------------------------------------
   set the compids[] vector using a global map from vertices
   to domains and interface nodes.

   DDmapIV -- IV object that contains the map from vertices
              to domains and interface nodes

   created -- 96mar17, cca
   ---------------------------------------------------------
*/
void
GPart_DDviaProjection (
   GPart   *gpart,
   IV      *DDmapIV
) {
int   *compids, *domainMap, *map, *vtxMap ;
int   dom, domloc, ndom, ndomloc, nvtx, vglob, vloc ;
/*
   ---------------
   check the input
   ---------------
*/
if ( gpart == NULL || DDmapIV == NULL ) {
   fprintf(stderr, "\n fatal error in GPart_DDviaProjection(%p,%p)"
           "\n bad input\n", gpart, DDmapIV) ;
   exit(-1) ;
}
nvtx    = gpart->nvtx ;
compids = IV_entries(&gpart->compidsIV) ;
/*
   --------------------------
   find the number of domains
   --------------------------
*/
vtxMap = IV_entries(&gpart->vtxMapIV) ;
map    = IV_entries(DDmapIV) ;
ndom   = IV_max(DDmapIV) ;
/*
   ------------------------
   check for a quick return
   ------------------------
*/
if ( gpart->par == NULL ) {
   IVcopy(nvtx, compids, map) ;
   gpart->ncomp = ndom ;
   return ;
}
/*
   ----------------------------------------
   fill compids[] with the local domain ids
   ----------------------------------------
*/
domainMap = IVinit(ndom+1, -1) ;
ndomloc = 0 ;
for ( vloc = 0 ; vloc < nvtx ; vloc++ ) {
   vglob = vtxMap[vloc] ;
   if ( (dom = map[vglob]) > 0 ) {
      if ( (domloc = domainMap[dom]) == -1 ) {
         domloc = domainMap[dom] = ++ndomloc ;
      }
      compids[vloc] = domloc ;
   } else {
      compids[vloc] = 0 ;
   }
}
gpart->ncomp = ndomloc ;
IVfree(domainMap) ; 

return ; }

/*
   ----------------------------------------------------------
   check that the permutation object does house a permutation

   return value --
      1 if a true permutation
      0 otherwise
   ----------------------------------------------------------
*/
int
Perm_checkPerm (
   Perm   *perm
) {
int   inew, iold, rc, size ;
int   *counts, *newToOld, *oldToNew ;
/*
   ---------------
   check the input
   ---------------
*/
if (  perm == NULL 
   || perm->isPresent < 1 || perm->isPresent > 3
   || (size = perm->size) <= 0 ) {
   fprintf(stderr, "\n fatal error in Perm_checkPerm(%p)"
           "\n bad input\n", perm) ;
   spoolesFatal();
}
rc = 1 ;
counts = IVinit(size, 0) ;
if ( (newToOld = perm->newToOld) != NULL ) {
   for ( inew = 0 ; inew < size ; inew++ ) {
      if ( 0 <= (iold = newToOld[inew]) && iold < size ) {
         counts[iold]++ ;
      } else {
         IVfree(counts) ;
         return(0) ;
      }
   }
   for ( iold = 0 ; iold < size ; iold++ ) {
      if ( counts[iold] != 1 ) {
         IVfree(counts) ;
         return(0) ;
      }
   }
}
if ( (oldToNew = perm->oldToNew) != NULL ) {
   IVzero(size, counts) ;
   for ( iold = 0 ; iold < size ; iold++ ) {
      if ( 0 <= (inew = oldToNew[iold]) && inew < size ) {
         counts[inew]++ ;
      } else {
         IVfree(counts) ;
         return(0) ;
      }
   }
   for ( inew = 0 ; inew < size ; inew++ ) {
      if ( counts[inew] != 1 ) {
         IVfree(counts) ;
         return(0) ;
      }
   }
}
IVfree(counts) ;

return(rc) ; }

/*
   -----------------------
   initialize the object

   created -- 95oct07, cca
   -----------------------
*/
void
BKL_init (
   BKL     *bkl,
   BPG     *bpg,
   float   alpha
) {
/*
   ---------------
   check the input
   ---------------
*/
if ( bkl == NULL || bpg == NULL ) {
   fprintf(stderr, "\n fatal error in BKL_init(%p,%p,%f)"
           "\n bad input\n", bkl, bpg, alpha) ;
   exit(-1) ;
}
/*
   --------------
   clear the data
   --------------
*/
BKL_clearData(bkl) ;
/*
   ---------------------
   initialize the fields
   ---------------------
*/
bkl->bpg  = bpg ;
bkl->ndom = bpg->nX ;
bkl->nseg = bpg->nY ;
bkl->nreg = bpg->nX + bpg->nY ;
if ( bpg->graph->vwghts == NULL ) {
   bkl->totweight = bkl->nreg ;
   bkl->regwghts  = IVinit(bkl->nreg, 1) ;
} else {
   bkl->regwghts  = bpg->graph->vwghts ;
   bkl->totweight = IVsum(bkl->nreg, bkl->regwghts) ; 
}
bkl->colors = IVinit(bkl->nreg, 0) ;
bkl->alpha  = alpha ;

return ; }

/*
   ------------------------------------------------------------------
   purpose -- set the scalars and allocate the vectors for the object
 
   created -- 98mar19, cca
   ------------------------------------------------------------------
*/
void
SolveMap_init (
   SolveMap   *solvemap,
   int        symmetryflag,
   int        nfront,
   int        nproc,
   int        nblockUpper,
   int        nblockLower
) {
/*
   ---------------
   check the input
   ---------------
*/
if ( solvemap == NULL || symmetryflag < 0 || nfront <= 0 
   || nproc < 0 || nblockUpper < 0 || nblockLower < 0 ) {
   fprintf(stderr, "\n fatal error in SolveMap_init(%p,%d,%d,%d,%d,%d)"
           "\n bad input\n", solvemap, symmetryflag, nfront, 
           nproc, nblockUpper, nblockLower) ;
   spoolesFatal();
}
/*
   ----------------
   clear the object
   ----------------
*/
SolveMap_clearData(solvemap) ;
/*
   ---------------
   set the scalars
   ---------------
*/
solvemap->symmetryflag = symmetryflag ;
solvemap->nfront       = nfront       ;
solvemap->nproc        = nproc        ;
solvemap->nblockUpper  = nblockUpper  ;
solvemap->nblockLower  = nblockLower  ;
/*
   ------------------------
   allocate the data arrays
   ------------------------
*/
solvemap->owners      = IVinit(nfront, -1) ;
solvemap->rowidsUpper = IVinit(nblockUpper, -1) ;
solvemap->colidsUpper = IVinit(nblockUpper, -1) ;
solvemap->mapUpper    = IVinit(nblockUpper, -1) ;
if (  symmetryflag == SPOOLES_NONSYMMETRIC && nblockLower > 0 ) {
   solvemap->rowidsLower = IVinit(nblockLower, -1) ;
   solvemap->colidsLower = IVinit(nblockLower, -1) ;
   solvemap->mapLower    = IVinit(nblockLower, -1) ;
}
return ; }

/*
   -------------------------------------
   return an IV object filled with the
   weights of the component's boundaries

   created -- 96oct21, cca
   -------------------------------------
*/
IV *
GPart_bndWeightsIV (
   GPart   *gpart 
) {
Graph   *graph ;
int     icomp, ii, ncomp, nvtx, v, vsize, vwght, w ;
int     *bnd, *compids, *cweights, *mark, *vadj, *vwghts ;
IV      *bndIV ;
/*
   ---------------
   check the input
   ---------------
*/
if ( gpart == NULL || (graph = gpart->g) == NULL ) {
   fprintf(stderr, "\n fatal error in GPart_bndWeightsIV(%p)"
           "\n bad input\n", gpart) ;
   exit(-1) ;
}
nvtx     = gpart->nvtx  ;
ncomp    = gpart->ncomp ;
compids  = IV_entries(&gpart->compidsIV)  ;
cweights = IV_entries(&gpart->cweightsIV) ;
vwghts   = graph->vwghts ;
bndIV    = IV_new() ;
IV_init(bndIV, 1 + ncomp, NULL) ;
IV_fill(bndIV, 0) ;
bnd  = IV_entries(bndIV) ;
mark = IVinit(ncomp+1, -1) ;
for ( v = 0 ; v < nvtx ; v++ ) {
   if ( compids[v] == 0 ) {
      vwght = (vwghts == NULL) ? 1 : vwghts[v] ;
      Graph_adjAndSize(graph, v, &vsize, &vadj) ;
      for ( ii = 0 ; ii < vsize ; ii++ ) {
         w = vadj[ii] ;
         if ( (icomp = compids[w]) != 0 && mark[icomp] != v ) {
            mark[icomp] = v ;
            bnd[icomp] += vwght ;
         }
      }
   }
}
IVfree(mark) ;

return(bndIV) ; }

/*
   ---------------------------------------------
   simplest initialization method

   if entries != NULL
      the object does not own the entries,
      it just points to the entries base address
   else if size > 0
      the object will own the entries, 
      it allocates a vector of size int's.
   else 
      nothing happens
   endif

   created -- 96aug28, cca
   ---------------------------------------------
*/
void
IV_init (
   IV    *iv,
   int   size,
   int   *entries 
) {
if ( iv == NULL || size < 0 ) {
   fprintf(stderr, "\n fatal error in IV_init(%p,%d,%p)"
           "\n bad input\n", iv, size, entries) ;
   exit(-1) ;
}
/*
   --------------
   clear any data
   --------------
*/
IV_clearData(iv) ;
/*
   -----------------------------
   set the size and maximum size
   -----------------------------
*/
iv->maxsize = iv->size = size ;
/*
   -------------------------
   set vector and owner flag
   -------------------------
*/
if ( entries != NULL ) {
   iv->owned = 0 ;
   iv->vec   = entries ; 
} else if ( size > 0 ) {
   iv->owned = 1 ;
   iv->vec   = IVinit(size, -1) ;
}
/*
fprintf(stdout, 
        "\n %% leaving IV_init, iv %p, size %d, maxsize %d, entries %p",
        iv, iv->size, iv->maxsize, iv->vec) ;
fflush(stdout) ;
*/

return ; }

/*
   ---------------------------------------
   purpose -- initialize the EGraph object

   created -- 96oct24, cca
   ---------------------------------------
*/
void
EGraph_init (
   EGraph   *egraph,
   int      type,
   int      nelem,
   int      nvtx,
   int      IVL_type
) {
/*
   ---------------
   check the input
   ---------------
*/
if (  egraph == NULL || type < 0 || type > 1 
   || nelem <= 0 || nvtx <= 0 ) {
   fprintf(stderr, "\n fatal error in EGraph_init(%p,%d,%d,%d,%d)"
           "\n bad input\n", egraph, type, nelem, nvtx, IVL_type) ;
   exit(-1) ;
}
/*
   ----------------------------
   clear the data in the object
   ----------------------------
*/
EGraph_clearData(egraph) ;
/*
   ---------------------
   initialize the object
   ---------------------
*/
egraph->type  = type ;
egraph->nelem = nelem ;
egraph->nvtx  = nvtx  ;
egraph->adjIVL = IVL_new() ;
IVL_init1(egraph->adjIVL, IVL_type, nelem) ;
if ( type == 1 ) {
   egraph->vwghts = IVinit(nvtx, 0) ;
}

return ; }

/*
   -----------------------------------------
   permute the columns of the matrix
   A(*,*) = A(*,index(*))
   this method calls A2_sortColumnsUp
   but does not overwrite the index[] vector

   created -- 98apr15, cca
   -----------------------------------------
*/
void
A2_permuteColumns (
   A2   *mtx,
   int   ncol,
   int   index[]
) {
int   *colids ;
/*
   ---------------
   check the input
   ---------------
*/
if ( mtx == NULL || ncol < 0 || ncol > mtx->n2 || index == NULL ) {
   fprintf(stderr, "\n fatal error in A2_permuteColumns(%p,%d,%p)"
           "\n bad input\n", mtx, ncol, index) ;
   exit(-1) ;
}
colids = IVinit(ncol, -1) ;
IVcopy(ncol, colids, index) ;
A2_sortColumnsUp(mtx, ncol, colids) ;
IVfree(colids) ;

return ; }

/*
   -------------------------------------------------------------
   purpose -- 

   the IVL object ivl and IV object ownersIV are both found on 
   each process.  the ownersIV object is identical over all the 
   processes, and owners[ii] tells which processes owns list ii 
   of the ivl object. on return from this method, the ivl object 
   is replicated over all the processes. each process sends 
   the lists that it owns to all the other processes.

   created -- 98apr03, cca
   -------------------------------------------------------------
*/
void
IVL_MPI_allgather (
   IVL        *ivl,
   IV         *ownersIV,
   int        stats[],
   int        msglvl,
   FILE       *msgFile,
   int        firsttag,
   MPI_Comm   comm
) {
int          count, destination, ii, ilist, incount, jlist, 
             jproc, left, maxcount, myid, nlist, nmylists, 
             notherlists, nowners, nproc, offset, outcount, 
             right, size, source, tag ;
int          *counts, *inbuffer, *list, *outbuffer, *owners ;
MPI_Status   status ;
/*
   ---------------
   check the input
   ---------------
*/
if ( ivl == NULL || ownersIV == NULL ) {
   fprintf(stderr, "\n fatal error in IVL_MPI_allgather()"
           "\n ivl = %p, ownersIV = %p\n",
           ivl, ownersIV) ;
   exit(-1) ;
}
/*
   ----------------------------------------------
   get id of self, # of processes and # of fronts
   ----------------------------------------------
*/
MPI_Comm_rank(comm, &myid) ;
MPI_Comm_size(comm, &nproc) ;
nlist = ivl->nlist ;
IV_sizeAndEntries(ownersIV, &nowners, &owners) ;
if ( msglvl > 1 ) {
   fprintf(msgFile, "\n\n inside IVL_MPI_allgather()"
           "\n nproc = %d, myid = %d, nlist = %d, nowners = %d",
           nproc, myid, nlist, nowners) ;
   fflush(msgFile) ;
}
if ( nlist != nowners || owners == NULL ) {
   fprintf(stderr, "\n fatal error in IVL_MPI_allgather()"
           "\n nlist = %d, nowners = %d, owners = %p\n",
           nlist, nowners, owners) ;
   exit(-1) ;
}
if ( msglvl > 2 ) {
   fprintf(msgFile, "\n\n ivl") ;
   IVL_writeForHumanEye(ivl, msgFile) ;
   fprintf(msgFile, "\n\n ownersIV") ;
   IV_writeForHumanEye(ownersIV, msgFile) ;
   fflush(msgFile) ;
}
/*
   -----------------------------------------------
   step 1 : determine the size of the message that
            this process will send to the others
   -----------------------------------------------
*/
for ( ilist = 0, outcount = 1 ; ilist < nlist ; ilist++ ) {
   if ( owners[ilist] < 0 || owners[ilist] >= nproc ) {
      fprintf(stderr, "\n owners[%d] = %d", ilist, owners[ilist]) ;
      exit(-1) ;
   }
   if ( owners[ilist] == myid ) {
      outcount += 2 ;
      IVL_listAndSize(ivl, ilist, &size, &list) ;
      outcount += size ;
   }
}
if ( msglvl > 2 ) {
   fprintf(msgFile, "\n\n outcount = %d", outcount) ;
   fflush(msgFile) ;
}
/*
   ----------------------------------------------------
   do an all-to-all gather/scatter
   counts[jproc] = # of int's in the message from jproc
   ----------------------------------------------------
*/
counts = IVinit(nproc, 0) ;
counts[myid] = outcount ;
MPI_Allgather((void *) &counts[myid], 1, MPI_INT,
              (void *) counts,  1, MPI_INT, comm) ;
//.........这里部分代码省略.........

/*--------------------------------------------------------------------*/
int
main ( int argc, char *argv[] )
/*
   ---------------------------------------------------------------
   read BPG from file and get the Dulmage-Mendelsohn decomposition

   created -- 96mar08, cca
   ---------------------------------------------------------------
*/
{
char     *inBPGFileName ;
double   t1, t2 ;
int      ierr, msglvl, rc ;
int      *dmflags, *stats ;
BPG      *bpg ;
FILE     *msgFile ;

if ( argc != 4 ) {
   fprintf(stdout, 
      "\n\n usage : %s msglvl msgFile inFile "
      "\n    msglvl   -- message level"
      "\n    msgFile  -- message file"
      "\n    inFile   -- input file, must be *.bpgf or *.bpgb"
      "\n", argv[0]) ;
   return(0) ;
}
msglvl = atoi(argv[1]) ;
if ( strcmp(argv[2], "stdout") == 0 ) {
   msgFile = stdout ;
} else if ( (msgFile = fopen(argv[2], "a")) == NULL ) {
   fprintf(stderr, "\n fatal error in %s"
           "\n unable to open file %s\n",
           argv[0], argv[2]) ;
   return(-1) ;
}
inBPGFileName  = argv[3] ;
fprintf(msgFile, 
        "\n %s "
        "\n msglvl   -- %d" 
        "\n msgFile  -- %s" 
        "\n inFile   -- %s" 
        "\n",
        argv[0], msglvl, argv[2], inBPGFileName) ;
fflush(msgFile) ;
/*
   ----------------------
   read in the BPG object
   ----------------------
*/
if ( strcmp(inBPGFileName, "none") == 0 ) {
   fprintf(msgFile, "\n no file to read from") ;
   exit(0) ;
}
bpg = BPG_new() ;
MARKTIME(t1) ;
rc = BPG_readFromFile(bpg, inBPGFileName) ;
MARKTIME(t2) ;
fprintf(msgFile, "\n CPU %9.5f : read in graph from file %s",
        t2 - t1, inBPGFileName) ;
if ( rc != 1 ) {
   fprintf(msgFile, "\n return value %d from BPG_readFromFile(%p,%s)",
           rc, bpg, inBPGFileName) ;
   exit(-1) ;
}
fprintf(msgFile, "\n\n after reading BPG object from file %s",
        inBPGFileName) ;
if ( msglvl > 2 ) {
   BPG_writeForHumanEye(bpg, msgFile) ;
} else {
   BPG_writeStats(bpg, msgFile) ;
}
fflush(msgFile) ;
/*
   --------------------------------------------
   test out the max flow DMdecomposition method
   --------------------------------------------
*/
dmflags = IVinit(bpg->nX + bpg->nY, -1) ;
stats   = IVinit(6, 0) ;
MARKTIME(t1) ;
BPG_DMviaMaxFlow(bpg, dmflags, stats, msglvl, msgFile) ;
MARKTIME(t2) ;
fprintf(msgFile, "\n\n CPU %9.5f : find DM via maxflow", t2 - t1) ;
if ( msglvl > 0 ) {
   fprintf(msgFile, 
           "\n\n BPG_DMviaMaxFlow"
           "\n |X_I| = %6d, |X_E| = %6d, |X_R| = %6d"
           "\n |Y_I| = %6d, |Y_E| = %6d, |Y_R| = %6d",
           stats[0], stats[1], stats[2],
           stats[3], stats[4], stats[5]) ;
}
if ( msglvl > 1 ) {
   fprintf(msgFile, "\n dmflags") ;
   IVfp80(msgFile, bpg->nX + bpg->nY, dmflags, 80, &ierr) ;
   fflush(msgFile) ;
}
/*
   ------------------------------------------
   test out the matching DMcomposition method
//.........这里部分代码省略.........

//.........这里部分代码省略.........
}
if ( symflag == 1 && dataType != 2 ) {
   fprintf(stderr, 
           "\n symflag = 1 (hermitian), dataType != 2 (complex)") ;
   spoolesFatal();
}
if ( nrowA <= 0 || ncolA <= 0 || nitem <= 0 ) {
   fprintf(stderr, 
           "\n invalid value: nrow = %d, ncol = %d, nitem = %d",
           nrowA, ncolA, nitem) ;
   spoolesFatal();
}
if ( symflag < 2 && nrowA != ncolA ) {
   fprintf(stderr,
           "\n invalid data: symflag = %d, nrow = %d, ncol = %d",
           symflag, nrowA, ncolA) ;
   spoolesFatal();
}
alpha[0] = alphaReal ;
alpha[1] = alphaImag ;
/*
   ----------------------------
   initialize the matrix object
   ----------------------------
*/
A = InpMtx_new() ;
InpMtx_init(A, storageMode, dataType, 0, 0) ;
drand = Drand_new() ;
/*
   ----------------------------------
   generate a vector of nitem triples
   ----------------------------------
*/
rowids = IVinit(nitem,   -1) ;
Drand_setUniform(drand, 0, nrowA) ;
Drand_fillIvector(drand, nitem, rowids) ;
colids = IVinit(nitem,   -1) ;
Drand_setUniform(drand, 0, ncolA) ;
Drand_fillIvector(drand, nitem, colids) ;
Drand_setUniform(drand, 0.0, 1.0) ;
if ( INPMTX_IS_REAL_ENTRIES(A) ) {
   zvec = DVinit(nitem, 0.0) ;
   Drand_fillDvector(drand, nitem, zvec) ;
} else if ( INPMTX_IS_COMPLEX_ENTRIES(A) ) {
   zvec = ZVinit(nitem, 0.0, 0.0) ;
   Drand_fillDvector(drand, 2*nitem, zvec) ;
}
/*
   -----------------------------------
   assemble the entries entry by entry
   -----------------------------------
*/
if ( msglvl > 1 ) {
   fprintf(msgFile, "\n\n A = zeros(%d,%d) ;", nrowA, ncolA) ;
}
if ( symflag == 1 ) {
/*
   ----------------
   hermitian matrix
   ----------------
*/
   for ( ii = 0 ; ii < nitem ; ii++ ) {
      if ( rowids[ii] == colids[ii] ) {
         zvec[2*ii+1] = 0.0 ;
      }
      if ( rowids[ii] <= colids[ii] ) {

//.........这里部分代码省略.........
   Chv_rowIndices(chvJ, &nrow, &rowind) ;
   IVramp(nrow, rowind, 0, 1) ;
}
if ( msglvl > 3 ) {
   fprintf(msgFile, "\n %% chevron a") ;
   Chv_writeForMatlab(chvJ, "a", msgFile) ;
   fflush(msgFile) ;
}
/*
   --------------------------
   initialize the chvI object
   --------------------------
*/
MARKTIME(t1) ;
chvI = Chv_new() ;
Chv_init(chvI, 0, nD, nL, nU, type, symflag) ;
MARKTIME(t2) ;
fprintf(msgFile, "\n %% CPU : %.3f to initialize matrix objects",
        t2 - t1) ;
Chv_zero(chvI) ;
Chv_columnIndices(chvI, &ncol, &colind) ;
IVramp(ncol, colind, 0, 1) ;
if ( CHV_IS_NONSYMMETRIC(chvI) ) {
   Chv_rowIndices(chvI, &nrow, &rowind) ;
   IVramp(nrow, rowind, 0, 1) ;
}
if ( symflag == 0 && pivotingflag == 1 ) {
/*
   ------------------------------
   create the pivotsizes[] vector
   ------------------------------
*/
   Drand_setUniform(drand, 1, 2.999) ;
   pivotsizes = IVinit(nD, 0) ;
   Drand_fillIvector(drand, nD, pivotsizes) ;
/*
   fprintf(msgFile, "\n initial pivotsizes[] : ") ;
   IVfp80(msgFile, nD, pivotsizes, 80, &ierr) ;
*/
   for ( npivot = count = 0 ; npivot < nD ; npivot++ ) {
      count += pivotsizes[npivot] ;
      if ( count > nD ) {
         pivotsizes[npivot]-- ;
         count-- ;
      } 
      if ( count == nD ) {
         break ;
      }
   }
   npivot++ ;
/*
   fprintf(msgFile, "\n final pivotsizes[] : ") ;
   IVfp80(msgFile, npivot, pivotsizes, 80, &ierr) ;
*/
} else {
   npivot = 0 ;
   pivotsizes = NULL ;
}
/*
   --------------------------------------------------
   first test: copy lower, diagonal and upper entries
   --------------------------------------------------
*/
if ( CHV_IS_NONSYMMETRIC(chvJ) ) {
   nentL = Chv_countEntries(chvJ, npivot, pivotsizes, CHV_STRICT_LOWER);
} else {

/*
   ------------------------------------------------------------
   create and return an ETree object that holds the front tree.

   created  -- 96jun23, cca
   ------------------------------------------------------------
*/
ETree *
MSMD_frontETree (
   MSMD   *msmd
) {
ETree     *etree ;
int       front, iv, nfront, nvtx, root ;
int       *bndwghts, *fch, *nodwghts, *par, *sib, *vtxToFront ;
MSMDvtx   *v, *w ;
/*
   ---------------
   check the input
   ---------------
*/
if ( msmd == NULL ) {
   fprintf(stderr, "\n fatal error in MSMD_frontETree(%p)"
           "\n bad input\n", msmd) ;
   exit(-1) ;
}
nvtx = msmd->nvtx ;
/*
   --------------------------
   count the number of fronts
   --------------------------
*/
nfront = 0 ;
fch = IVinit(nvtx, -1) ;
sib = IVinit(nvtx, -1) ;
root = -1 ;
for ( iv = 0, v = msmd->vertices ; iv < nvtx ; iv++, v++ ) {
#if MYDEBUG > 0
   fprintf(stdout, "\n vertex %d, status %c, wght %d",
           v->id, v->status, v->wght) ;
/*
   MSMDvtx_print(v, stdout) ;
*/
   fflush(stdout) ;
#endif
   switch ( v->status ) {
   case 'L' :
   case 'E' :
      if ( (w = v->par) != NULL ) {
         sib[v->id] = fch[w->id] ;
         fch[w->id] = v->id ;
      } else {
         sib[v->id] = root ;
         root = v->id ;
      }
#if MYDEBUG > 0
   fprintf(stdout, ", new front %d", nfront) ;
   fflush(stdout) ;
#endif
      nfront++ ;
      break ;
   default :
      break ;
   }
}
#if MYDEBUG > 0
fprintf(stdout, "\n %d fronts", nfront) ;
fflush(stdout) ;
#endif
/*
   ---------------------------
   initialize the ETree object
   ---------------------------
*/
etree = ETree_new() ;
ETree_init1(etree, nfront, nvtx) ;
nodwghts   = IV_entries(etree->nodwghtsIV) ;
bndwghts   = IV_entries(etree->bndwghtsIV) ;
vtxToFront = IV_entries(etree->vtxToFrontIV) ;
/*
   ----------------------------------------------
   fill the vtxToFront[] vector so representative 
   vertices are mapped in a post-order traversal
   ----------------------------------------------
*/
nfront = 0 ;
iv = root ;
while ( iv != -1 ) {
   while ( fch[iv] != -1 ) {
       iv = fch[iv] ;
   }
   v = msmd->vertices + iv ;
   vtxToFront[iv] = nfront++ ;
#if MYDEBUG > 0
   fprintf(stdout, "\n v = %d, vwght = %d, vtxToFront[%d] = %d", 
           v->id, v->wght, iv, vtxToFront[iv]) ;
   fflush(stdout) ;
#endif
   while ( sib[iv] == -1 && v->par != NULL ) {
      v = v->par ;
      iv = v->id ;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
if ( myid == root ) {
   InpMtx   *inpmtx ;
   int      nedges, totewght, totvwght, v ;
   int      *adj, *vwghts ;
   IVL      *adjIVL, *ewghtIVL ;
/*
   -----------------------
   generate a random graph
   -----------------------
*/
   inpmtx = InpMtx_new() ;
   InpMtx_init(inpmtx, INPMTX_BY_ROWS, INPMTX_INDICES_ONLY, nitem, 0) ;
   Drand_setDefaultFields(&drand) ;
   Drand_setSeed(&drand, seed) ;
   Drand_setUniform(&drand, 0, nvtx) ;
   Drand_fillIvector(&drand, nitem, InpMtx_ivec1(inpmtx)) ;
   Drand_fillIvector(&drand, nitem, InpMtx_ivec2(inpmtx)) ;
   InpMtx_setNent(inpmtx, nitem) ;
   InpMtx_changeStorageMode(inpmtx, INPMTX_BY_VECTORS) ;
   if ( msglvl > 2 ) {
      fprintf(msgFile, "\n\n inpmtx mtx filled with raw entries") ;
      InpMtx_writeForHumanEye(inpmtx, msgFile) ;
      fflush(msgFile) ;
   }
   adjIVL = InpMtx_fullAdjacency(inpmtx) ;
   if ( msglvl > 2 ) {
      fprintf(msgFile, "\n\n full adjacency structure") ;
      IVL_writeForHumanEye(adjIVL, msgFile) ;
      fflush(msgFile) ;
   }
   nedges = adjIVL->tsize ;
   if ( type == 1 || type == 3 ) {
      Drand_setUniform(&drand, 1, 10) ;
      vwghts = IVinit(nvtx, 0) ;
      Drand_fillIvector(&drand, nvtx, vwghts) ;
      totvwght = IVsum(nvtx, vwghts) ;
      if ( msglvl > 2 ) {
         fprintf(msgFile, "\n\n vertex weights") ;
         IVfprintf(msgFile, nvtx, vwghts) ;
         fflush(msgFile) ;
      }
   } else {
      vwghts = NULL ;
      totvwght = nvtx ;
   }
   if ( msglvl > 2 ) {
      fprintf(msgFile, "\n\n totvwght %d", totvwght) ;
      fflush(msgFile) ;
   }
   if ( type == 2 || type == 3 ) {
      ewghtIVL = IVL_new() ;
      IVL_init1(ewghtIVL, IVL_CHUNKED, nvtx) ;
      Drand_setUniform(&drand, 1, 100) ;
      totewght = 0 ;
      for ( v = 0 ; v < nvtx ; v++ ) {
         IVL_listAndSize(adjIVL, v, &size, &adj) ;
         IVL_setList(ewghtIVL, v, size, NULL) ;
         IVL_listAndSize(ewghtIVL, v, &size, &adj) ;
         Drand_fillIvector(&drand, size, adj) ;
         totewght += IVsum(size, adj) ;
      }
      if ( msglvl > 2 ) {
         fprintf(msgFile, "\n\n ewghtIVL") ;
         IVL_writeForHumanEye(ewghtIVL, msgFile) ;
         fflush(msgFile) ;
      }

/*
   -------------------------------------------------------
   purpose -- to read in a Graph object from a CHACO file
 
   input --
 
      fn -- filename
 
   return value -- 1 if success, 0 if failure
 
   created -- 98sep20, jjs
   --------------------------------------------------------
*/
int
Graph_readFromChacoFile (
   Graph   *graph,
   char    *fn
) {
char    *rc ;
FILE    *fp;
int     nvtx, nedges, format;
char    string[BUFLEN], *s1, *s2;
int     k, v, vsize, w, vwghts, ewghts;
int     *adjncy, *weights, *vwghtsINT;
IVL     *adjIVL, *ewghtIVL;
/*
   ---------------
   check the input
   ---------------
*/
if ((graph == NULL) || (fn == NULL)) {
   fprintf(stderr, "\n error in Graph_readFromFile(%p,%s)"
           "\n bad input\n", graph, fn);
   return(0);
}
/*
   ---------------------
   clear the data fields
   ---------------------
*/
Graph_clearData(graph);
/*
   ----------------------------------------------
   open file and read in nvtx, nedges, and format
   ----------------------------------------------
*/
if ((fp = fopen(fn, "r")) == (FILE*)NULL) {
   fprintf(stderr, "\n error in Graph_readFromChacoFile(%p,%s)"
           "\n unable to open file %s", graph, fn, fn);
   return(0);
}
/*
   -------------
   skip comments
   -------------
*/
do {
   rc = fgets(string, BUFLEN, fp) ;
   if ( rc == NULL ) {
      fprintf(stderr, "\n error in Graph_readFromChacoFile()"
             "\n error skipping comments in file %s\n", fn) ;
      return(0) ;
   }
} while ( string[0] == '%');
/*
   -------------------------------------------------
   read in # vertices, # edges and (optional) format
   -------------------------------------------------
*/
format = 0;
if (sscanf(string, "%d %d %d", &nvtx, &nedges, &format) < 2) {
   fprintf(stderr, "\n error in Graph_readFromChacoFile(%p,%s)"
           "\n unable to read header of file %s", graph, fn, fn);
   return(0);
}
ewghts = ((format % 10) > 0);
vwghts = (((format / 10) % 10) > 0);
if (format >= 100) {
   fprintf(stderr, "\n error in Graph_readFromChacoFile(%p,%s)"
           "\n unknown format", graph, fn);
   return(0);
}
/*
   ------------------------------------------------------------------
   initialize vector(s) to hold adjacency and (optional) edge weights
   ------------------------------------------------------------------
*/
adjncy = IVinit(nvtx, -1) ;
if ( ewghts ) {
   weights = IVinit(nvtx, -1) ;
} else {
   weights = NULL ;
}
/*
   ---------------------------
   initialize the Graph object
   ---------------------------
*/
nedges *= 2;
nedges += nvtx;
//.........这里部分代码省略.........

/*
   -------------------------------------------------
   sort the columns of the matrix in ascending order
   of the colids[] vector. on return, colids is
   in asending order. return value is the number
   of column swaps made.

   created -- 98apr15, cca
   -------------------------------------------------
*/
int
A2_sortColumnsUp (
   A2   *mtx,
   int   ncol,
   int   colids[]
) {
int   ii, mincol, mincolid, nswap, target ;
/*
   ---------------
   check the input
   ---------------
*/
if ( mtx == NULL || mtx->n2 < ncol || ncol < 0 || colids == NULL ) {
   fprintf(stderr, "\n fatal error in A2_sortColumnsUp(%p,%d,%p)"
           "\n bad input\n", mtx, ncol, colids) ;
   if ( mtx != NULL ) {
      A2_writeStats(mtx, stderr) ;
   }
   exit(-1) ;
}
if ( ! (A2_IS_REAL(mtx) || A2_IS_COMPLEX(mtx)) ) {
   fprintf(stderr, "\n fatal error in A2_sortColumnsUp(%p,%d,%p)"
           "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n", 
           mtx, ncol, colids, mtx->type) ;
   exit(-1) ;
}
nswap = 0 ;
if ( mtx->inc2 == 1 ) {
   double   *dvtmp ;
   int      irow, nrow ;
   int      *ivtmp ;
/*
   ---------------------------------------------------
   matrix is stored by rows, so permute each row
   ---------------------------------------------------
*/
   ivtmp = IVinit(ncol, -1) ;
   if ( A2_IS_REAL(mtx) ) {
      dvtmp = DVinit(ncol, 0.0) ;
   } else if ( A2_IS_COMPLEX(mtx) ) {
      dvtmp = DVinit(2*ncol, 0.0) ;
   }
   IVramp(ncol, ivtmp, 0, 1) ;
   IV2qsortUp(ncol, colids, ivtmp) ;
   nrow = mtx->n1 ;
   for ( irow = 0 ; irow < nrow ; irow++ ) {
      if ( A2_IS_REAL(mtx) ) {
         DVcopy(ncol, dvtmp, A2_row(mtx, irow)) ;
         DVgather(ncol, A2_row(mtx, irow), dvtmp, ivtmp) ;
      } else if ( A2_IS_COMPLEX(mtx) ) {
         ZVcopy(ncol, dvtmp, A2_row(mtx, irow)) ;
         ZVgather(ncol, A2_row(mtx, irow), dvtmp, ivtmp) ;
      }
   }
   IVfree(ivtmp) ;
   DVfree(dvtmp) ;
} else {
/*
   ----------------------------------------
   use a simple insertion sort to swap cols
   ----------------------------------------
*/
   for ( target = 0 ; target < ncol ; target++ ) {
      mincol   = target ;
      mincolid = colids[target] ;
      for ( ii = target + 1 ; ii < ncol ; ii++ ) {
         if ( mincolid > colids[ii] ) {
            mincol   = ii ;
            mincolid = colids[ii] ;
         }
      }
      if ( mincol != target ) {
         colids[mincol] = colids[target] ;
         colids[target] = mincolid ;
         A2_swapColumns(mtx, target, mincol) ;
         nswap++ ;
      }
   }
}
return(nswap) ; }

/*
   ---------------------------------------------------
   purpose -- move the solution from the individual
     SubMtx objects into the global solution SubMtx object
 
   created -- 98feb20
   ---------------------------------------------------
*/
void
FrontMtx_storeSolution (
   FrontMtx        *frontmtx,
   int             owners[],
   int             myid,
   SubMtxManager   *manager,
   SubMtx          *p_mtx[],
   DenseMtx        *solmtx,
   int             msglvl,
   FILE            *msgFile
) {
char     localsol ;
SubMtx   *xmtxJ ;
double   *sol, *xJ ;
int      inc1, inc2, irow, jrhs, J, kk,
         ncolJ, neqns, nfront, nJ, nrhs, nrowInSol, nrowJ ;
int      *colindJ, *colmap, *rowind ;

if ( (nrowInSol = solmtx->nrow) != (neqns = frontmtx->neqns) ) {
/*
   --------------------------------------------------------------
   the solution matrix is only part of the total solution matrix.
   (this happens in an MPI environment where the rhs
   is partitioned among the processors.)
   create a map from the global row indices to the
   indices local to this solution matrix.
   --------------------------------------------------------------
*/
   colmap = IVinit(neqns, -1) ;
   rowind = solmtx->rowind ;
   if ( msglvl > 1 ) {
      fprintf(msgFile, "\n solmtx->rowind") ;
      IVfprintf(msgFile, solmtx->nrow, rowind) ;
      fflush(msgFile) ;
   }
   for ( irow = 0 ; irow < nrowInSol ; irow++ ) {
      colmap[rowind[irow]] = irow ;
   }
   localsol = 'T' ;
   if ( msglvl > 1 ) {
      fprintf(msgFile, "\n colmap") ;
      IVfprintf(msgFile, neqns, colmap) ;
      fflush(msgFile) ;
   }
} else {
   localsol = 'F' ;
}
DenseMtx_dimensions(solmtx, &neqns, &nrhs) ;
nfront = FrontMtx_nfront(frontmtx) ;
for ( J = 0 ; J < nfront ; J++ ) {
   if (  (owners == NULL || owners[J] == myid)
      && (nJ = FrontMtx_frontSize(frontmtx, J)) > 0 ) {
      FrontMtx_columnIndices(frontmtx, J, &ncolJ, &colindJ) ;
      xmtxJ = p_mtx[J] ;
      if ( xmtxJ == NULL ) {
         fprintf(stderr,
            "\n fatal error in storeSolution(%d)"
            "\n thread %d, xmtxJ = NULL", J, myid) ;
         exit(-1) ;
      }
      if ( msglvl > 1 ) {
         fprintf(msgFile, "\n storing solution for front %d", J) ;
         SubMtx_writeForHumanEye(xmtxJ, msgFile) ;
         fflush(msgFile) ;
      }
      if ( localsol == 'T' ) {
/*
        ------------------------------------------------------
         map the global row indices into the local row indices
        ------------------------------------------------------
*/
         if ( msglvl > 1 ) {
            fprintf(msgFile, "\n global row indices") ;