/* solveUtil.c */
#include "../FrontMtx.h"
#include "../../SubMtxList.h"
/*--------------------------------------------------------------------*/
static SubMtx * initBJ ( int type, int J, int nJ, int nrhs,
SubMtxManager *mtxmanager, int msglvl, FILE *msgFile ) ;
static void computeForwardUpdates ( FrontMtx *frontmtx, SubMtx *BJ,
int J, IP *heads[], char frontIsDone[], SubMtx *p_mtx[],
int msglvl, FILE *msgFile ) ;
static void assembleAggregates ( int J, SubMtx *BJ, SubMtxList *aggList,
SubMtxManager *mtxmanager, int msglvl, FILE *msgFile ) ;
static void computeBackwardUpdates ( FrontMtx *frontmtx, SubMtx *ZJ,
int J, IP *heads[], char frontIsDone[], SubMtx *p_mtx[],
int msglvl, FILE *msgFile ) ;
/*--------------------------------------------------------------------*/
/*
---------------------------------------------
load the right hand side for the owned fronts
created -- 98mar19, cca
---------------------------------------------
*/
SubMtx **
FrontMtx_loadRightHandSide (
FrontMtx *frontmtx,
DenseMtx *rhsmtx,
int owners[],
int myid,
SubMtxManager *mtxmanager,
int msglvl,
FILE *msgFile
) {
char localrhs ;
SubMtx *BJ ;
SubMtx **p_agg ;
double *bJ, *rhs ;
int inc1, inc2, irow, jrhs, J, kk, nbytes, ncolJ,
neqns, nfront, nJ, nrhs, nrowInRhs, nrowJ ;
int *rowind, *rowindJ, *rowmap ;
nrowInRhs = rhsmtx->nrow ;
neqns = frontmtx->neqns ;
if ( nrowInRhs != neqns ) {
/*
----------------------------------------------------
the rhs matrix is only part of the total rhs 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 rhs matrix.
----------------------------------------------------
*/
rowmap = IVinit(neqns, -1) ;
rowind = rhsmtx->rowind ;
if ( msglvl > 2 ) {
fprintf(msgFile, "\n rhsmtx->rowind") ;
IVfprintf(msgFile, rhsmtx->nrow, rowind) ;
fflush(msgFile) ;
}
for ( irow = 0 ; irow < nrowInRhs ; irow++ ) {
rowmap[rowind[irow]] = irow ;
}
if ( msglvl > 1 ) {
fprintf(msgFile, "\n rowmap") ;
IVfprintf(msgFile, neqns, rowmap) ;
fflush(msgFile) ;
}
localrhs = 'T' ;
} else {
localrhs = 'F' ;
}
/*
--------------------------------------
allocate the vector of SubMtx pointers
--------------------------------------
*/
nfront = FrontMtx_nfront(frontmtx) ;
ALLOCATE(p_agg, struct _SubMtx *, nfront) ;
for ( J = 0 ; J < nfront ; J++ ) {
p_agg[J] = NULL ;
}
DenseMtx_dimensions(rhsmtx, &neqns, &nrhs) ;
for ( J = 0 ; J < nfront ; J++ ) {
if ( (owners == NULL || owners[J] == myid)
&& (nJ = FrontMtx_frontSize(frontmtx, J)) > 0 ) {
FrontMtx_rowIndices(frontmtx, J, &nrowJ, &rowindJ) ;
if ( localrhs == 'T' ) {
/*
------------------------------------------------------
map the global row indices into the local row indices
------------------------------------------------------
*/
for ( irow = 0 ; irow < nJ ; irow++ ) {
rowindJ[irow] = rowmap[rowindJ[irow]] ;
}
}
/*
------------------------------------------------
initialize bmtx, a SubMtx object to hold b_{J,*}
------------------------------------------------
*/
nbytes = SubMtx_nbytesNeeded(frontmtx->type, SUBMTX_DENSE_COLUMNS,
nJ, nrhs, nJ*nrhs);
BJ = SubMtxManager_newObjectOfSizeNbytes(mtxmanager, nbytes) ;
SubMtx_init(BJ, frontmtx->type, SUBMTX_DENSE_COLUMNS,
J, 0, nJ, nrhs, nJ*nrhs) ;
p_agg[J] = BJ ;
if ( BJ == NULL ) {
fprintf(stderr,
"\n fatal error in load rhs(%d), BJ = NULL", J) ;
exit(-1) ;
}
/*
-----------------------
load BJ with b_{J,*}
-----------------------
*/
rhs = DenseMtx_entries(rhsmtx) ;
SubMtx_denseInfo(BJ, &nrowJ, &ncolJ, &inc1, &inc2, &bJ) ;
if ( FRONTMTX_IS_REAL(frontmtx) ) {
for ( jrhs = 0 ; jrhs < nrhs ; jrhs++ ) {
for ( irow = 0 ; irow < nJ ; irow++ ) {
kk = rowindJ[irow] ;
bJ[irow] = rhs[kk] ;
}
bJ += nJ ;
rhs += nrowInRhs ;
}
} else if ( FRONTMTX_IS_COMPLEX(frontmtx) ) {
for ( jrhs = 0 ; jrhs < nrhs ; jrhs++ ) {
for ( irow = 0 ; irow < nJ ; irow++ ) {
kk = rowindJ[irow] ;
bJ[2*irow] = rhs[2*kk] ;
bJ[2*irow+1] = rhs[2*kk+1] ;
}
bJ += 2*nJ ;
rhs += 2*nrowInRhs ;
}
}
if ( msglvl > 3 ) {
fprintf(msgFile, "\n\n rhs for J = %d, BJ = %p", J, BJ) ;
fflush(msgFile) ;
SubMtx_writeForHumanEye(BJ, msgFile) ;
fflush(msgFile) ;
}
if ( localrhs == 'T' ) {
/*
-----------------------------------------------------------
map the local row indices back into the global row indices
-----------------------------------------------------------
*/
for ( irow = 0 ; irow < nJ ; irow++ ) {
rowindJ[irow] = rowind[rowindJ[irow]] ;
}
}
}
}
if ( localrhs == 'T' ) {
IVfree(rowmap) ;
}
return(p_agg) ; }
/*--------------------------------------------------------------------*/
/*
--------------------------------------
visit front J during the forward solve
created -- 98mar27, cca
--------------------------------------
*/
void
FrontMtx_forwardVisit (
FrontMtx *frontmtx,
int J,
int nrhs,
int *owners,
int myid,
SubMtxManager *mtxmanager,
SubMtxList *aggList,
SubMtx *p_mtx[],
char frontIsDone[],
IP *heads[],
SubMtx *p_agg[],
char status[],
int msglvl,
FILE *msgFile
) {
char aggDone, updDone ;
SubMtx *BJ, *LJJ, *UJJ ;
int nJ ;
if ( (nJ = FrontMtx_frontSize(frontmtx, J)) == 0 ) {
/*
-----------------------------------------------------
front has no eliminated rows or columns, quick return
-----------------------------------------------------
*/
if ( owners == NULL || owners[J] == myid ) {
frontIsDone[J] = 'Y' ;
}
status[J] = 'F' ;
return ;
}
if ( heads[J] != NULL ) {
/*
-------------------------------------
there are internal updates to perform
-------------------------------------
*/
if ( (BJ = p_agg[J]) == NULL ) {
/*
---------------------------
create the aggregate object
---------------------------
*/
BJ = p_agg[J] = initBJ(frontmtx->type, J, nJ, nrhs,
mtxmanager, msglvl, msgFile) ;
}
if ( msglvl > 3 ) {
fprintf(msgFile, "\n\n BJ = %p", BJ) ;
fflush(msgFile) ;
SubMtx_writeForHumanEye(BJ, msgFile) ;
fflush(msgFile) ;
}
/*
---------------------------
compute any waiting updates
---------------------------
*/
computeForwardUpdates(frontmtx, BJ, J, heads, frontIsDone, p_mtx,
msglvl, msgFile) ;
}
if ( heads[J] == NULL ) {
updDone = 'Y' ;
} else {
updDone = 'N' ;
}
if ( aggList != NULL && owners[J] == myid ) {
/*
-----------------------
assemble any aggregates
-----------------------
*/
aggDone = 'N' ;
if ( (BJ = p_agg[J]) == NULL ) {
fprintf(stderr,
"\n 2. fatal error in forwardVisit(%d), BJ = NULL", J) ;
exit(-1) ;
}
assembleAggregates(J, BJ, aggList, mtxmanager,
msglvl, msgFile) ;
if ( SubMtxList_isCountZero(aggList, J) == 1 ) {
if ( msglvl > 3 ) {
fprintf(msgFile, "\n\n aggregate count is zero") ;
fflush(msgFile) ;
}
assembleAggregates(J, BJ, aggList, mtxmanager,
msglvl, msgFile) ;
aggDone = 'Y' ;
}
} else {
aggDone = 'Y' ;
}
if ( msglvl > 3 ) {
fprintf(msgFile, "\n\n updDone = %c, aggDone = %c",
updDone, aggDone) ;
fflush(msgFile) ;
}
if ( updDone == 'Y' && aggDone == 'Y' ) {
BJ = p_agg[J] ;
if ( owners == NULL || owners[J] == myid ) {
/*
-------------------------------------
owned front, ready for interior solve
-------------------------------------
*/
if ( FRONTMTX_IS_NONSYMMETRIC(frontmtx) ) {
LJJ = FrontMtx_lowerMtx(frontmtx, J, J) ;
if ( LJJ != NULL ) {
SubMtx_solve(LJJ, BJ) ;
}
} else {
UJJ = FrontMtx_upperMtx(frontmtx, J, J) ;
if ( UJJ != NULL ) {
if ( FRONTMTX_IS_SYMMETRIC(frontmtx) ) {
SubMtx_solveT(UJJ, BJ) ;
} else if ( FRONTMTX_IS_HERMITIAN(frontmtx) ) {
SubMtx_solveH(UJJ, BJ) ;
}
}
}
if ( msglvl > 1 ) {
fprintf(msgFile, "\n\n after forward solve") ;
SubMtx_writeForHumanEye(BJ, msgFile) ;
fflush(msgFile) ;
}
/*
------------------------------------------------
move YJ (stored in BJ) into p_mtx[],
signal front as done, and set status to finished
------------------------------------------------
*/
p_agg[J] = NULL ;
p_mtx[J] = BJ ;
frontIsDone[J] = 'Y' ;
} else if ( BJ != NULL ) {
/*
--------------------------------------
unowned front, put into aggregate list
--------------------------------------
*/
if ( msglvl > 3 ) {
fprintf(msgFile, "\n\n putting BJ into aggregate list") ;
fflush(msgFile) ;
}
SubMtxList_addObjectToList(aggList, BJ, J) ;
p_agg[J] = NULL ;
}
status[J] = 'F' ;
}
return ; }
/*--------------------------------------------------------------------*/
/*
--------------------------------------------------
purpose -- visit front J during the diagonal solve
created -- 98feb20, cca
--------------------------------------------------
*/
void
FrontMtx_diagonalVisit (
FrontMtx *frontmtx,
int J,
int owners[],
int myid,
SubMtx *p_mtx[],
char frontIsDone[],
SubMtx *p_agg[],
int msglvl,
FILE *msgFile
) {
if ( owners == NULL || owners[J] == myid ) {
SubMtx *BJ, *DJJ ;
if ( (BJ = p_mtx[J]) != NULL ) {
if ( msglvl > 3 ) {
fprintf(msgFile, "\n\n BJ = %p", BJ) ;
SubMtx_writeForHumanEye(BJ, msgFile) ;
fflush(msgFile) ;
}
DJJ = FrontMtx_diagMtx(frontmtx, J) ;
if ( msglvl > 3 ) {
fprintf(msgFile, "\n\n DJJ = %p", DJJ) ;
SubMtx_writeForHumanEye(DJJ, msgFile) ;
fflush(msgFile) ;
}
SubMtx_solve(DJJ, BJ) ;
if ( msglvl > 1 ) {
fprintf(msgFile, "\n\n b_{%d,*} after diagonal solve", J) ;
SubMtx_writeForHumanEye(BJ, msgFile) ;
fflush(msgFile) ;
}
p_mtx[J] = NULL ;
p_agg[J] = BJ ;
}
frontIsDone[J] = 'Y' ;
}
return ; }
/*--------------------------------------------------------------------*/
/*
---------------------------------------
visit front J during the backward solve
created -- 98mar27, cca
---------------------------------------
*/
void
FrontMtx_backwardVisit (
FrontMtx *frontmtx,
int J,
int nrhs,
int *owners,
int myid,
SubMtxManager *mtxmanager,
SubMtxList *aggList,
SubMtx *p_mtx[],
char frontIsDone[],
IP *heads[],
SubMtx *p_agg[],
char status[],
int msglvl,
FILE *msgFile
) {
char aggDone, updDone ;
SubMtx *UJJ, *ZJ ;
int nJ ;
if ( msglvl > 1 ) {
fprintf(msgFile, "\n inside FrontMtx_backwardVisit(%d), nJ = %d",
J, FrontMtx_frontSize(frontmtx, J)) ;
fflush(msgFile) ;
}
if ( (nJ = FrontMtx_frontSize(frontmtx, J)) == 0 ) {
/*
-----------------------------------------------------
front has no eliminated rows or columns, quick return
-----------------------------------------------------
*/
if ( owners == NULL || owners[J] == myid ) {
frontIsDone[J] = 'Y' ;
}
status[J] = 'F' ;
return ;
}
if ( msglvl > 1 ) {
fprintf(msgFile, "\n heads[%d] = %p", J, heads[J]) ;
fflush(msgFile) ;
}
if ( heads[J] != NULL ) {
/*
-------------------------------------
there are internal updates to perform
-------------------------------------
*/
if ( (ZJ = p_agg[J]) == NULL ) {
/*
---------------------------
create the aggregate object
---------------------------
*/
ZJ = p_agg[J] = initBJ(frontmtx->type, J, nJ, nrhs,
mtxmanager, msglvl, msgFile) ;
}
if ( msglvl > 3 ) {
fprintf(msgFile, "\n\n ZJ = %p", ZJ) ;
SubMtx_writeForHumanEye(ZJ, msgFile) ;
fflush(msgFile) ;
}
/*
---------------------------
compute any waiting updates
---------------------------
*/
computeBackwardUpdates(frontmtx, ZJ, J, heads, frontIsDone, p_mtx,
msglvl, msgFile) ;
}
if ( heads[J] == NULL ) {
updDone = 'Y' ;
} else {
updDone = 'N' ;
}
if ( msglvl > 1 ) {
fprintf(msgFile, "\n updDone = %c", updDone) ;
fflush(msgFile) ;
}
if ( aggList != NULL && owners[J] == myid ) {
/*
-----------------------
assemble any aggregates
-----------------------
*/
aggDone = 'N' ;
if ( (ZJ = p_agg[J]) == NULL ) {
fprintf(stderr,
"\n 2. fatal error in backwardVisit(%d), ZJ = NULL", J) ;
exit(-1) ;
}
assembleAggregates(J, ZJ, aggList, mtxmanager,
msglvl, msgFile) ;
if ( SubMtxList_isCountZero(aggList, J) == 1 ) {
if ( msglvl > 3 ) {
fprintf(msgFile, "\n\n aggregate count is zero") ;
fflush(msgFile) ;
}
assembleAggregates(J, ZJ, aggList, mtxmanager,
msglvl, msgFile) ;
aggDone = 'Y' ;
}
} else {
aggDone = 'Y' ;
}
if ( msglvl > 1 ) {
fprintf(msgFile, "\n aggDone = %c", aggDone) ;
fflush(msgFile) ;
}
if ( updDone == 'Y' && aggDone == 'Y' ) {
ZJ = p_agg[J] ;
if ( owners == NULL || owners[J] == myid ) {
/*
-------------------------------------
owned front, ready for interior solve
-------------------------------------
*/
UJJ = FrontMtx_upperMtx(frontmtx, J, J) ;
if ( UJJ != NULL ) {
SubMtx_solve(UJJ, ZJ) ;
}
if ( msglvl > 1 ) {
fprintf(msgFile, "\n\n after backward solve") ;
SubMtx_writeForHumanEye(ZJ, msgFile) ;
fflush(msgFile) ;
}
/*
------------------------------------------------
move YJ (stored in BJ) into p_mtx[],
signal front as done, and set status to finished
------------------------------------------------
*/
p_agg[J] = NULL ;
p_mtx[J] = ZJ ;
frontIsDone[J] = 'Y' ;
} else if ( ZJ != NULL ) {
/*
--------------------------------------
unowned front, put into aggregate list
--------------------------------------
*/
SubMtxList_addObjectToList(aggList, ZJ, J) ;
p_agg[J] = NULL ;
}
status[J] = 'F' ;
}
if ( msglvl > 1 ) {
fprintf(msgFile, "\n status[%d] = %c", J, status[J]) ;
fflush(msgFile) ;
}
return ; }
/*--------------------------------------------------------------------*/
/*
---------------------------------------------------
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") ;
IVfprintf(msgFile, nJ, colindJ) ;
fflush(msgFile) ;
}
for ( irow = 0 ; irow < nJ ; irow++ ) {
colindJ[irow] = colmap[colindJ[irow]] ;
}
if ( msglvl > 1 ) {
fprintf(msgFile, "\n local row indices") ;
IVfprintf(msgFile, nJ, colindJ) ;
fflush(msgFile) ;
}
}
/*
----------------------------------
store x_{J,*} into solution matrix
----------------------------------
*/
sol = DenseMtx_entries(solmtx) ;
SubMtx_denseInfo(xmtxJ, &nrowJ, &ncolJ, &inc1, &inc2, &xJ) ;
if ( FRONTMTX_IS_REAL(frontmtx) ) {
for ( jrhs = 0 ; jrhs < nrhs ; jrhs++ ) {
for ( irow = 0 ; irow < nJ ; irow++ ) {
kk = colindJ[irow] ;
sol[kk] = xJ[irow] ;
}
sol += neqns ;
xJ += nJ ;
}
} else if ( FRONTMTX_IS_COMPLEX(frontmtx) ) {
for ( jrhs = 0 ; jrhs < nrhs ; jrhs++ ) {
for ( irow = 0 ; irow < nJ ; irow++ ) {
kk = colindJ[irow] ;
sol[2*kk] = xJ[2*irow] ;
sol[2*kk+1] = xJ[2*irow+1] ;
}
sol += 2*neqns ;
xJ += 2*nJ ;
}
}
/*
fprintf(msgFile, "\n solution for front %d stored", J) ;
*/
SubMtxManager_releaseObject(manager, xmtxJ) ;
if ( localsol == 'T' ) {
/*
-----------------------------------------------------------
map the local row indices back into the global row indices
-----------------------------------------------------------
*/
for ( irow = 0 ; irow < nJ ; irow++ ) {
colindJ[irow] = rowind[colindJ[irow]] ;
}
}
}
}
if ( localsol == 'T' ) {
IVfree(colmap) ;
}
/*
fprintf(msgFile, "\n\n SOLUTION") ;
DenseMtx_writeForHumanEye(solmtx, msgFile) ;
*/
return ; }
/*--------------------------------------------------------------------*/
/*
----------------------------------------------
purpose -- initialize an aggregate SubMtx object
created -- 98mar27, cca
----------------------------------------------
*/
static SubMtx *
initBJ (
int type,
int J,
int nJ,
int nrhs,
SubMtxManager *mtxmanager,
int msglvl,
FILE *msgFile
) {
SubMtx *BJ ;
double *entries ;
int inc1, inc2, nbytes ;
/*
------------------------------------------
B_J not yet allocated (must not be owned),
create and zero the entries
------------------------------------------
*/
nbytes = SubMtx_nbytesNeeded(type, SUBMTX_DENSE_COLUMNS,
nJ, nrhs, nJ*nrhs);
BJ = SubMtxManager_newObjectOfSizeNbytes(mtxmanager, nbytes) ;
if ( BJ == NULL ) {
fprintf(stderr,
"\n 1. fatal error in forwardVisit(%d), BJ = NULL", J) ;
exit(-1) ;
}
SubMtx_init(BJ, type, SUBMTX_DENSE_COLUMNS, J, 0, nJ, nrhs, nJ*nrhs) ;
SubMtx_denseInfo(BJ, &nJ, &nrhs, &inc1, &inc2, &entries) ;
if ( type == SPOOLES_REAL ) {
DVzero(nJ*nrhs, entries) ;
} else if ( type == SPOOLES_COMPLEX ) {
DVzero(2*nJ*nrhs, entries) ;
}
return(BJ) ; }
/*--------------------------------------------------------------------*/
/*
------------------------------------
purpose -- compute any updates to BJ
created -- 98mar26, cca
------------------------------------
*/
static void
computeForwardUpdates (
FrontMtx *frontmtx,
SubMtx *BJ,
int J,
IP *heads[],
char frontIsDone[],
SubMtx *p_mtx[],
int msglvl,
FILE *msgFile
) {
SubMtx *LJI, *UIJ, *YI ;
int I ;
IP *ip, *nextip ;
/*
-------------------------------
loop over the remaining updates
-------------------------------
*/
for ( ip = heads[J], heads[J] = NULL ;
ip != NULL ;
ip = nextip ) {
I = ip->val ; nextip = ip->next ;
if ( msglvl > 3 ) {
fprintf(msgFile,
"\n\n frontIsDone[%d] = %c, p_mtx[%d] = %p",
I, frontIsDone[I], I, p_mtx[I]) ;
fflush(msgFile) ;
}
if ( frontIsDone[I] == 'Y' ) {
if ( (YI = p_mtx[I]) != NULL ) {
/*
--------------------------------
Y_I exists and has been computed
--------------------------------
*/
if ( msglvl > 3 ) {
fprintf(msgFile, "\n\n before solve: YI = %p", YI) ;
SubMtx_writeForHumanEye(YI, msgFile) ;
fflush(msgFile) ;
}
if ( FRONTMTX_IS_NONSYMMETRIC(frontmtx) ) {
if ( (LJI = FrontMtx_lowerMtx(frontmtx, J, I)) != NULL ) {
if ( msglvl > 3 ) {
fprintf(msgFile, "\n\n LJI = %p", LJI) ;
SubMtx_writeForHumanEye(LJI, msgFile) ;
fflush(msgFile) ;
}
SubMtx_solveupd(BJ, LJI, YI) ;
}
} else {
if ( (UIJ = FrontMtx_upperMtx(frontmtx, I, J)) != NULL ) {
if ( msglvl > 3 ) {
fprintf(msgFile, "\n\n UIJ = %p", UIJ) ;
SubMtx_writeForHumanEye(UIJ, msgFile) ;
fflush(msgFile) ;
}
if ( FRONTMTX_IS_SYMMETRIC(frontmtx) ) {
SubMtx_solveupdT(BJ, UIJ, YI) ;
} else if ( FRONTMTX_IS_HERMITIAN(frontmtx) ) {
SubMtx_solveupdH(BJ, UIJ, YI) ;
}
}
}
if ( msglvl > 3 ) {
fprintf(msgFile, "\n\n after update: BJ = %p", BJ) ;
SubMtx_writeForHumanEye(BJ, msgFile) ;
fflush(msgFile) ;
}
}
} else {
/*
------------------------
Y_I is not yet available
------------------------
*/
ip->next = heads[J] ;
heads[J] = ip ;
}
}
return ; }
/*--------------------------------------------------------------------*/
/*
--------------------------------------------------------
purpose -- assemble any aggregates in the aggregate list
created -- 98mar26, cca
--------------------------------------------------------
*/
static void
assembleAggregates (
int J,
SubMtx *BJ,
SubMtxList *aggList,
SubMtxManager *mtxmanager,
int msglvl,
FILE *msgFile
) {
SubMtx *BJhat, *BJhead ;
double *entBJ, *entBJhat ;
int inc1, inc1hat, inc2, inc2hat, ncol, ncolhat, nrow, nrowhat ;
if ( BJ == NULL || aggList == NULL ) {
fprintf(stderr,
"\n fatal error in assembleAggregates()"
"\n BJ = %p, aggList = %p", BJ, aggList) ;
exit(-1) ;
}
if ( SubMtxList_isListNonempty(aggList, BJ->rowid) ) {
if ( msglvl > 3 ) {
fprintf(msgFile, "\n\n aggregate list is not-empty") ;
fflush(msgFile) ;
}
SubMtx_denseInfo(BJ, &nrow, &ncol, &inc1, &inc2, &entBJ) ;
if ( msglvl > 2 ) {
fprintf(msgFile,
"\n\n BJ(%d,%d) : nrow %d, ncol %d, inc1 %d, inc2 %d, ent %p",
BJ->rowid, BJ->colid, nrow, ncol, inc1, inc2, entBJ) ;
fflush(msgFile) ;
}
BJhead = SubMtxList_getList(aggList, J) ;
for ( BJhat = BJhead ; BJhat != NULL ; BJhat = BJhat->next ) {
if ( BJhat == NULL ) {
fprintf(stderr,
"\n 3. fatal error in forwardVisit(%d)"
"\n BJhat = NULL", J) ;
exit(-1) ;
}
SubMtx_denseInfo(BJhat, &nrowhat, &ncolhat, &inc1hat, &inc2hat,
&entBJhat) ;
if ( msglvl > 2 ) {
fprintf(msgFile,
"\n BJhat(%d,%d) : nrow %d, ncol %d, inc1 %d, inc2 %d, ent %p",
BJhat->rowid, BJhat->colid,
nrowhat, ncolhat, inc1hat, inc2hat, entBJhat) ;
fflush(msgFile) ;
}
if ( nrow != nrowhat || ncol != ncolhat
|| inc1 != inc1hat || inc2 != inc2hat || entBJhat == NULL ) {
fprintf(stderr, "\n fatal error") ;
exit(-1) ;
}
if ( msglvl > 3 ) {
fprintf(msgFile, "\n\n BJ") ;
SubMtx_writeForHumanEye(BJ, msgFile) ;
fprintf(msgFile, "\n\n BJhat") ;
SubMtx_writeForHumanEye(BJhat, msgFile) ;
fflush(msgFile) ;
}
if ( SUBMTX_IS_REAL(BJhat) ) {
DVadd(nrow*ncol, entBJ, entBJhat) ;
} else if ( SUBMTX_IS_COMPLEX(BJhat) ) {
DVadd(2*nrow*ncol, entBJ, entBJhat) ;
}
}
SubMtxManager_releaseListOfObjects(mtxmanager, BJhead) ;
if ( msglvl > 3 ) {
fprintf(msgFile, "\n\n BJ after assembly") ;
SubMtx_writeForHumanEye(BJ, msgFile) ;
fflush(msgFile) ;
}
}
return ; }
/*--------------------------------------------------------------------*/
/*
------------------------------------
purpose -- compute any updates to ZJ
created -- 98mar26, cca
------------------------------------
*/
static void
computeBackwardUpdates (
FrontMtx *frontmtx,
SubMtx *ZJ,
int J,
IP *heads[],
char frontIsDone[],
SubMtx *p_mtx[],
int msglvl,
FILE *msgFile
) {
SubMtx *UJK, *XK ;
int K ;
IP *ip, *nextip ;
/*
-------------------------------
loop over the remaining updates
-------------------------------
*/
for ( ip = heads[J], heads[J] = NULL ;
ip != NULL ;
ip = nextip ) {
K = ip->val ; nextip = ip->next ;
if ( msglvl > 3 ) {
fprintf(msgFile,
"\n\n frontIsDone[%d] = %c", K, frontIsDone[K]) ;
fflush(msgFile) ;
}
if ( frontIsDone[K] == 'Y' ) {
if ( (XK = p_mtx[K]) != NULL ) {
/*
--------------------------------
X_K exists and has been computed
--------------------------------
*/
if ( msglvl > 3 ) {
fprintf(msgFile, "\n\n before solve: XK = %p", XK) ;
SubMtx_writeForHumanEye(XK, msgFile) ;
fflush(msgFile) ;
}
if ( (UJK = FrontMtx_upperMtx(frontmtx, J, K)) != NULL ) {
if ( msglvl > 3 ) {
fprintf(msgFile, "\n\n UJK = %p", UJK) ;
SubMtx_writeForHumanEye(UJK, msgFile) ;
fflush(msgFile) ;
}
SubMtx_solveupd(ZJ, UJK, XK) ;
}
if ( msglvl > 3 ) {
fprintf(msgFile, "\n\n after update: ZJ = %p", ZJ) ;
SubMtx_writeForHumanEye(ZJ, msgFile) ;
fflush(msgFile) ;
}
}
} else {
/*
------------------------
X_K is not yet available
------------------------
*/
ip->next = heads[J] ;
heads[J] = ip ;
}
}
return ; }
/*--------------------------------------------------------------------*/
/*
---------------------------------------------------
purpose -- load the dequeue with the roots of the
active subtree used for a backward solve
created -- 98mar27, cca
---------------------------------------------------
*/
void
FrontMtx_loadActiveRoots (
FrontMtx *frontmtx,
char status[],
char activeFlag,
Ideq *dequeue
) {
int J ;
int *sib ;
sib = frontmtx->tree->sib ;
for ( J = frontmtx->tree->root ; J != -1 ; J = sib[J] ) {
if ( status[J] == activeFlag ) {
Ideq_insertAtTail(dequeue, J) ;
}
}
return ; }
/*--------------------------------------------------------------------*/
/*
--------------------------------------------------------------------
purpose -- to set up the link data structures for a forward solve
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 -- 98feb20, cca
--------------------------------------------------------------------
*/
IP **
FrontMtx_forwardSetup (
FrontMtx *frontmtx,
int msglvl,
FILE *msgFile
) {
int ii, J, K, nadj, nblock, nfront ;
int *adj ;
IP *ip ;
IP **heads ;
/*
--------------------------------------------------
set up the update head/links for the forward solve
--------------------------------------------------
*/
nblock = FrontMtx_nLowerBlocks(frontmtx) ;
nfront = FrontMtx_nfront(frontmtx) ;
ALLOCATE(heads, struct _IP *, nfront + 2) ;
for ( J = 0 ; J <= nfront + 1 ; J++ ) {
heads[J] = NULL ;
}
heads[nfront] = heads[nfront+1] = IP_init(nblock, 1) ;
for ( J = 0 ; J < nfront ; J++ ) {
FrontMtx_lowerAdjFronts(frontmtx, J, &nadj, &adj) ;
for ( ii = 0 ; ii < nadj ; ii++ ) {
if ( (K = adj[ii]) > 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) ;
}
}
}
}
return(heads) ; }
/*--------------------------------------------------------------------*/
/*
---------------------------------------------------------
purpose -- set up the linked lists for the backward solve
created -- 98feb20, cca
---------------------------------------------------------
*/
IP **
FrontMtx_backwardSetup (
FrontMtx *frontmtx,
int msglvl,
FILE *msgFile
) {
IP *ip ;
IP **heads ;
int ii, J, K, nadj, nblock, nfront ;
int *adj ;
nfront = FrontMtx_nfront(frontmtx) ;
nblock = FrontMtx_nLowerBlocks(frontmtx) ;
ALLOCATE(heads, struct _IP *, nfront + 2) ;
for ( J = 0 ; J <= nfront + 1 ; J++ ) {
heads[J] = NULL ;
}
heads[nfront] = heads[nfront+1] = IP_init(nblock, 1) ;
for ( J = 0 ; J < nfront ; J++ ) {
FrontMtx_upperAdjFronts(frontmtx, J, &nadj, &adj) ;
for ( ii = 0 ; ii < nadj ; ii++ ) {
if ( (K = adj[ii]) > J ) {
if ( heads[nfront] == NULL ) {
fprintf(stderr,
"\n WHOA, heads[nfront] = %p", heads[nfront]) ;
exit(-1) ;
}
ip = heads[nfront] ; heads[nfront] = ip->next ;
ip->val = K ; ip->next = heads[J] ; heads[J] = ip ;
if ( msglvl > 3 ) {
fprintf(msgFile, "\n linking U(%d,%d) to U(%d,%d)",
J, K, J, (ip->next == NULL) ? -1 : ip->next->val) ;
fflush(msgFile) ;
}
}
}
}
return(heads) ; }
/*--------------------------------------------------------------------*/