/*  SetupMT.c  */

#include "../BridgeMT.h"

#define MYDEBUG 1

#if MYDEBUG > 0
static int count_Setup = 0 ;
static double time_Setup = 0.0 ;
#endif

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

   given InpMtx objects that contain A and B, initialize the bridge
   data structure for the multithreaded factor's, solve's and mvm's.

   data -- pointer to a Bridge object
   pprbtype -- pointer to value containing problem type
     *prbtype = 1 --> A X = B X Lambda, vibration problem
     *prbtype = 2 --> A X = B X Lambda, buckling problem
     *prbtype = 3 --> A X = X Lambda, simple eigenvalue problem
   pneqns   -- pointer to value containing number of equations
   pmxbsz   -- pointer to value containing blocksize
   A        -- pointer to InpMtx object containing A
   B        -- pointer to InpMtx object containing B
   pseed    -- pointer to value containing a random number seed
   pmsglvl  -- pointer to value containing a message level
   msgFile  -- message file pointer
   pnthread -- pointer to value containing # of threads

   return value --
      1 -- normal return
     -1 -- data is NULL
     -2 -- pprbtype is NULL
     -3 -- *pprbtype is invalid
     -4 -- pneqns is NULL
     -5 -- *pneqns is invalid
     -6 -- pmxbsz is NULL
     -7 -- *pmxbsz is invalid
     -8 -- A and B are NULL
     -9 -- pseed is NULL
    -10 -- pmsglvl is NULL
    -11 -- *pmsglvl > 0 and msgFile is NULL
    -12 -- pnthreads is NULL
    -13 -- *pnthreads is invalid

   created -- 98aug10, cca
   ------------------------------------------------------------------
*/
int
SetupMT (
   void     *data,
   int      *pprbtype,
   int      *pneqns,
   int      *pmxbsz,
   InpMtx   *A,
   InpMtx   *B,
   int      *pseed,
   int      *pnthread,
   int      *pmsglvl,
   FILE     *msgFile
) {
BridgeMT   *bridge = (BridgeMT *) data ;
double     sigma[2] ;
DV         *cumopsDV ;
Graph      *graph ;
int        maxdomainsize, maxsize, maxzeros, msglvl, mxbsz,
           nedges, neqns, nthread, prbtype, seed ;
IVL        *adjIVL ;
#if MYDEBUG > 0
double   t1, t2 ;
MARKTIME(t1) ;
count_Setup++ ;
fprintf(stdout, "\n (%d) SetupMT()", count_Setup) ;
fflush(stdout) ;
#endif
/*
   --------------------
   check the input data
   --------------------
*/
if ( data == NULL ) {
   fprintf(stderr, "\n fatal error in SetupMT()"
           "\n data is NULL\n") ;
   return(-1) ;
}
if ( pprbtype == NULL ) {
   fprintf(stderr, "\n fatal error in SetupMT()"
           "\n pprbtype is NULL\n") ;
   return(-2) ;
}
prbtype = *pprbtype ;
if ( prbtype < 1 || prbtype > 3 ) {
   fprintf(stderr, "\n fatal error in SetupMT()"
           "\n prbtype = %d, is invalid\n", prbtype) ;
   return(-3) ;
}
if ( pneqns == NULL ) {
   fprintf(stderr, "\n fatal error in SetupMT()"
           "\n pneqns is NULL\n") ;
   return(-4) ;
}
neqns = *pneqns ;
if ( neqns <= 0 ) {
   fprintf(stderr, "\n fatal error in SetupMT()"
           "\n neqns = %d, is invalid\n", neqns) ;
   return(-5) ;
}
if ( pmxbsz == NULL ) {
   fprintf(stderr, "\n fatal error in SetupMT()"
           "\n pmxbsz is NULL\n") ;
   return(-6) ;
}
mxbsz = *pmxbsz ;
if ( mxbsz <= 0 ) {
   fprintf(stderr, "\n fatal error in SetupMT()"
           "\n *pmxbsz = %d, is invalid\n", mxbsz) ;
   return(-7) ;
}
if ( A == NULL && B == NULL ) {
   fprintf(stderr, "\n fatal error in SetupMT()"
           "\n A and B are NULL\n") ;
   return(-8) ;
}
if ( pseed == NULL ) {
   fprintf(stderr, "\n fatal error in SetupMT()"
           "\n pseed is NULL\n") ;
   return(-9) ;
}
seed = *pseed ;
if ( pmsglvl == NULL ) {
   fprintf(stderr, "\n fatal error in SetupMT()"
           "\n pmsglvl is NULL\n") ;
   return(-10) ;
}
msglvl = *pmsglvl ;
if ( msglvl > 0 && msgFile == NULL ) {
   fprintf(stderr, "\n fatal error in SetupMT()"
           "\n msglvl = %d, msgFile = NULL\n", msglvl) ;
   return(-11) ;
}
if ( pnthread == NULL ) {
   fprintf(stderr, "\n fatal error in SetupMT()"
           "\n pnthread is NULL\n") ;
   return(-12) ;
}
nthread = *pnthread ;
if ( nthread <= 0 ) {
   fprintf(stderr, "\n fatal error in SetupMT()"
           "\n nthread = %d, is invalid\n", nthread) ;
   return(-13) ;
}
bridge->msglvl  = *pmsglvl ;
/*
bridge->msglvl  = 2 ;
*/
bridge->msgFile = msgFile ;
if ( msglvl > 2 ) {
   fprintf(msgFile, "\n\n inside SetupMT()"
     "\n neqns = %d, prbtype = %d, mxbsz = %d, seed = %d, nthread = %d",
     neqns, prbtype, mxbsz, seed, nthread) ;
   if ( A != NULL ) {
      fprintf(msgFile, "\n\n matrix A") ;
      InpMtx_writeForHumanEye(A, msgFile) ;
   }
   if ( B != NULL ) {
      fprintf(msgFile, "\n\n matrix B") ;
      InpMtx_writeForHumanEye(B, msgFile) ;
   }
   fflush(msgFile) ;
}
bridge->prbtype = prbtype ;
bridge->neqns   = neqns   ;
bridge->mxbsz   = mxbsz   ;
bridge->A       = A       ;
bridge->B       = B       ;
bridge->seed    = seed    ;
bridge->nthread = nthread ;
/*
   ----------------------------
   create and initialize pencil
   ----------------------------
*/
sigma[0] = 1.0; sigma[1] = 0.0;
bridge->pencil = Pencil_new() ;
Pencil_setDefaultFields(bridge->pencil) ;
Pencil_init(bridge->pencil, SPOOLES_REAL, SPOOLES_SYMMETRIC,
            A, sigma, B) ;
/*
   --------------------------------
   convert to row or column vectors
   --------------------------------
*/
if ( A != NULL ) {
   if ( ! INPMTX_IS_BY_ROWS(A) && ! INPMTX_IS_BY_COLUMNS(A) ) {
      InpMtx_changeCoordType(A, INPMTX_BY_ROWS) ;
   }
   if ( ! INPMTX_IS_BY_VECTORS(A) ) {
      InpMtx_changeStorageMode(A, INPMTX_BY_VECTORS) ;
   }
}
if ( B != NULL ) {
   if ( ! INPMTX_IS_BY_ROWS(B) && ! INPMTX_IS_BY_COLUMNS(B) ) {
      InpMtx_changeCoordType(B, INPMTX_BY_ROWS) ;
   }
   if ( ! INPMTX_IS_BY_VECTORS(B) ) {
      InpMtx_changeStorageMode(B, INPMTX_BY_VECTORS) ;
   }
}
/*
   -------------------------------
   create a Graph object for A + B
   -------------------------------
*/
graph  = Graph_new() ;
adjIVL = Pencil_fullAdjacency(bridge->pencil);
nedges = IVL_tsize(adjIVL),
Graph_init2(graph, 0, bridge->neqns, 0, nedges,
            bridge->neqns, nedges, adjIVL, NULL, NULL) ;
if ( msglvl > 2 ) {
   fprintf(msgFile, "\n\n graph of the input matrix") ;
   Graph_writeForHumanEye(graph, msgFile) ;
   fflush(msgFile) ;
}
/*
   ---------------
   order the graph
   ---------------
*/
maxdomainsize = neqns / 64 ;
if ( maxdomainsize == 0 ) {
   maxdomainsize = 1 ;
}
maxzeros  = (int) (0.01*neqns) ;
maxsize   = 64 ;
bridge->frontETree = orderViaBestOfNDandMS(graph, maxdomainsize,
                     maxzeros, maxsize, bridge->seed, msglvl, msgFile) ;
if ( msglvl > 2 ) {
   fprintf(msgFile, "\n\n front tree from ordering") ;
   ETree_writeForHumanEye(bridge->frontETree, msgFile) ;
   fflush(msgFile) ;
}
/*
   ----------------------------------------------
   get the old-to-new and new-to-old permutations
   ----------------------------------------------
*/
bridge->oldToNewIV = ETree_oldToNewVtxPerm(bridge->frontETree) ;
bridge->newToOldIV = ETree_newToOldVtxPerm(bridge->frontETree) ;
if ( msglvl > 2 ) {
   fprintf(msgFile, "\n\n old-to-new permutation") ;
   IV_writeForHumanEye(bridge->oldToNewIV, msgFile) ;
   fprintf(msgFile, "\n\n new-to-old permutation") ;
   IV_writeForHumanEye(bridge->newToOldIV, msgFile) ;
   fflush(msgFile) ;
}
/*
   --------------------------------------
   permute the vertices in the front tree
   --------------------------------------
*/
ETree_permuteVertices(bridge->frontETree, bridge->oldToNewIV) ;
if ( msglvl > 2 ) {
   fprintf(msgFile, "\n\n permuted front etree") ;
   ETree_writeForHumanEye(bridge->frontETree, msgFile) ;
   fflush(msgFile) ;
}
/*
   -------------------------------------------
   permute the entries in the pencil. 
   note, after the permutation the 
   entries are mapped into the upper triangle.
   -------------------------------------------
*/
Pencil_permute(bridge->pencil, bridge->oldToNewIV, bridge->oldToNewIV) ;
Pencil_mapToUpperTriangle(bridge->pencil) ;
Pencil_changeCoordType(bridge->pencil, INPMTX_BY_CHEVRONS) ;
Pencil_changeStorageMode(bridge->pencil, INPMTX_BY_VECTORS) ;
if ( msglvl > 2 ) {
   fprintf(msgFile, "\n\n permuted pencil") ;
   Pencil_writeForHumanEye(bridge->pencil, msgFile) ;
   fflush(msgFile) ;
}
/*
   ----------------------------------
   compute the symbolic factorization
   ----------------------------------
*/
bridge->symbfacIVL = SymbFac_initFromPencil(bridge->frontETree, 
                                            bridge->pencil) ;
if ( msglvl > 2 ) {
   fprintf(msgFile, "\n\n symbolic factorization") ;
   IVL_writeForHumanEye(bridge->symbfacIVL, msgFile) ;
   fflush(msgFile) ;
}
/*
   --------------------------------------------------
   create a FrontMtx object to hold the factorization
   --------------------------------------------------
*/
bridge->frontmtx = FrontMtx_new() ;
/*
   ------------------------------------------------------------
   create a SubMtxManager object to hold the factor submatrices
   ------------------------------------------------------------
*/
bridge->mtxmanager = SubMtxManager_new() ;
SubMtxManager_init(bridge->mtxmanager, LOCK_IN_PROCESS, 0) ;
/*
   ------------------------------------------------------------
   allocate the working objects X and Y for the matrix multiply
   ------------------------------------------------------------
*/
bridge->X = DenseMtx_new() ;
DenseMtx_init(bridge->X, SPOOLES_REAL, 0, 0, neqns, mxbsz, 1, neqns) ;
bridge->Y = DenseMtx_new() ;
DenseMtx_init(bridge->Y, SPOOLES_REAL, 0, 0, neqns, mxbsz, 1, neqns) ;
/*
   --------------------------------------------------------
   setup up ownersIV, the map from fronts to owning threads
   --------------------------------------------------------
*/
cumopsDV = DV_new() ;
DV_init(cumopsDV, bridge->nthread, NULL );
bridge->ownersIV = ETree_ddMap(bridge->frontETree, SPOOLES_REAL,
                               SPOOLES_SYMMETRIC, cumopsDV,
                               1./(2.*bridge->nthread) );
if ( msglvl >= 0 ) {
   fprintf(msgFile, "\n\n load balance for the processors") ;
   DV_writeForHumanEye(cumopsDV, msgFile) ;
   fprintf(msgFile, "\n\n map from fronts to threads") ;
   IV_writeForHumanEye(bridge->ownersIV, msgFile) ;
   fflush(msgFile) ;
}
/*
   -------------------------------------------------------------------
   create a SolveMap object, the map from submatries to owning threads
   -------------------------------------------------------------------
*/
bridge->solvemap = SolveMap_new() ;
/*
   ------------------------
   free the working storage
   ------------------------
*/
Graph_free(graph) ;
DV_free(cumopsDV);

#if MYDEBUG > 0
MARKTIME(t2) ;
time_Setup += t2 - t1 ;
fprintf(stdout, ", %8.3f seconds, %8.3f total time",
        t2 - t1, time_Setup) ;
fflush(stdout) ;
#endif

return(1) ; }

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