/* storeFront.c */
#include "../FrontMtx.h"
/*--------------------------------------------------------------------*/
/*
----------------------------------------------------------------
purpose -- to store the factor indicies and entries from front J
pivotsizesIV -- used for symmetric or hermitian and pivoting
droptol -- used for drop tolerance factorization,
an entry is stored if its magnitude > droptol
created -- 98may04, cca
----------------------------------------------------------------
*/
void
FrontMtx_storeFront (
FrontMtx *frontmtx,
Chv *frontJ,
IV *pivotsizesIV,
double droptol,
int msglvl,
FILE *msgFile
) {
SubMtx *mtx ;
double *entries ;
int inc1, inc2, ipivot, irow, jcol, J, m, nbytes,
ncol, nD, nent, nfront, nL, npivot, nrow, nU ;
int *colind, *colindJ, *firstlocs, *indices, *pivots,
*pivotsizes, *rowind, *rowindJ, *sizes ;
SubMtxManager *manager ;
/*
-------------------------------------
extract information from front matrix
-------------------------------------
*/
nfront = frontmtx->nfront ;
manager = frontmtx->manager ;
if ( ( FRONTMTX_IS_SYMMETRIC(frontmtx)
|| FRONTMTX_IS_HERMITIAN(frontmtx) )
&& FRONTMTX_IS_PIVOTING(frontmtx) ) {
IV_sizeAndEntries(pivotsizesIV, &npivot, &pivotsizes) ;
} else {
npivot = 0, pivotsizes = NULL ;
}
/*
--------------------------------
extract information from chevron
--------------------------------
*/
J = frontJ->id ;
Chv_dimensions(frontJ, &nD, &nL, &nU) ;
Chv_columnIndices(frontJ, &ncol, &colindJ) ;
Chv_rowIndices(frontJ, &nrow, &rowindJ) ;
if ( msglvl > 2 ) {
fprintf(msgFile, "\n inside FrontMtx_storeFront, J = %d", J) ;
fflush(msgFile) ;
}
if ( FRONTMTX_IS_PIVOTING(frontmtx) ) {
if ( frontmtx->lock != NULL ) {
if ( msglvl > 2 ) {
fprintf(msgFile, "\n locking lock") ;
fflush(msgFile) ;
}
Lock_lock(frontmtx->lock) ;
frontmtx->nlocks++ ;
}
/*
--------------------
store the front size
--------------------
*/
FrontMtx_setFrontSize(frontmtx, J, nD) ;
/*
------------------------
store the column indices
------------------------
*/
if ( msglvl > 2 ) {
fprintf(msgFile, "\n setting column indices, J = %d", J) ;
fflush(msgFile) ;
}
IVL_setList(frontmtx->coladjIVL, J, ncol, colindJ) ;
if ( FRONTMTX_IS_NONSYMMETRIC(frontmtx) ) {
/*
---------------------
store the row indices
---------------------
*/
if ( msglvl > 2 ) {
fprintf(msgFile, "\n setting row indices, J = %d", J) ;
fflush(msgFile) ;
}
IVL_setList(frontmtx->rowadjIVL, J, nrow, rowindJ) ;
}
if ( frontmtx->lock != NULL ) {
if ( msglvl > 2 ) {
fprintf(msgFile, "\n unlocking lock") ;
fflush(msgFile) ;
}
Lock_unlock(frontmtx->lock) ;
}
}
if ( nD == 0 ) {
return ;
}
/*
------------------------
store the D_{J,J} matrix
------------------------
*/
if ( pivotsizes != NULL ) {
/*
---------------------------------------------------
symmetric and pivoting, store block diagonal matrix
---------------------------------------------------
*/
nent = Chv_countEntries(frontJ, npivot, pivotsizes, CHV_DIAGONAL) ;
nbytes = SubMtx_nbytesNeeded(frontJ->type, SUBMTX_BLOCK_DIAGONAL_SYM,
nD, nD, nent) ;
mtx = SubMtxManager_newObjectOfSizeNbytes(manager, nbytes) ;
if ( FRONTMTX_IS_SYMMETRIC(frontmtx) ) {
SubMtx_init(mtx, frontJ->type, SUBMTX_BLOCK_DIAGONAL_SYM,
J, J, nD, nD, nent) ;
} else if ( FRONTMTX_IS_HERMITIAN(frontmtx) ) {
SubMtx_init(mtx, frontJ->type, SUBMTX_BLOCK_DIAGONAL_HERM,
J, J, nD, nD, nent) ;
}
SubMtx_blockDiagonalInfo(mtx, &nrow, &nent, &pivots, &entries) ;
IVzero(nD, pivots) ;
IVcopy(npivot, pivots, pivotsizes) ;
Chv_copyEntriesToVector(frontJ, npivot, pivotsizes, nent, entries,
CHV_DIAGONAL, CHV_BY_ROWS) ;
if ( msglvl > 2 ) {
fprintf(msgFile, "\n\n block diagonal matrix") ;
SubMtx_writeForHumanEye(mtx, msgFile) ;
fflush(msgFile) ;
}
frontmtx->nentD += nent ;
} else {
/*
---------------
diagonal matrix
---------------
*/
mtx = frontmtx->p_mtxDJJ[J] ;
if ( mtx == NULL ) {
nbytes = SubMtx_nbytesNeeded(frontJ->type, SUBMTX_DIAGONAL,
nD, nD, nD) ;
mtx = SubMtxManager_newObjectOfSizeNbytes(manager, nbytes) ;
SubMtx_init(mtx, frontJ->type, SUBMTX_DIAGONAL, J, J, nD, nD, nD);
}
SubMtx_diagonalInfo(mtx, &nent, &entries) ;
Chv_copyEntriesToVector(frontJ, npivot, pivotsizes, nent, entries,
CHV_DIAGONAL, CHV_BY_ROWS) ;
frontmtx->nentD += nD ;
if ( msglvl > 2 ) {
fprintf(msgFile, "\n\n entries in SubMtx") ;
DVfprintf(msgFile, nent, entries) ;
fprintf(msgFile, "\n\n block diagonal matrix") ;
SubMtx_writeForHumanEye(mtx, msgFile) ;
fflush(msgFile) ;
}
}
frontmtx->p_mtxDJJ[J] = mtx ;
SubMtx_columnIndices(mtx, &ncol, &colind) ;
IVcopy(ncol, colind, colindJ) ;
SubMtx_rowIndices(mtx, &nrow, &rowind) ;
IVcopy(nrow, rowind, rowindJ) ;
/*
------------------------------
store the upper U_{J,J} matrix
------------------------------
*/
if ( FRONTMTX_IS_DENSE_FRONTS(frontmtx) ) {
/*
-----------
dense front
-----------
*/
nent = Chv_countEntries(frontJ, npivot, pivotsizes,
CHV_STRICT_UPPER_11) ;
if ( nent > 0 ) {
mtx = frontmtx->p_mtxUJJ[J] ;
if ( mtx == NULL ) {
nbytes = SubMtx_nbytesNeeded(frontJ->type,
SUBMTX_DENSE_SUBCOLUMNS, nD, nD,nent);
mtx = SubMtxManager_newObjectOfSizeNbytes(manager, nbytes) ;
SubMtx_init(mtx, frontJ->type, SUBMTX_DENSE_SUBCOLUMNS,
J, J, nD, nD, nent) ;
}
SubMtx_denseSubcolumnsInfo(mtx, &ncol, &nent,
&firstlocs, &sizes, &entries) ;
if ( pivotsizes != NULL ) {
for ( jcol = ipivot = 0 ; jcol < nD ; ipivot++ ) {
m = pivotsizes[ipivot] ;
if ( m == 1 ) {
firstlocs[jcol] = 0 ;
sizes[jcol] = jcol ;
jcol++ ;
} else if ( m == 2 ) {
firstlocs[jcol] = firstlocs[jcol+1] = 0 ;
sizes[jcol] = sizes[jcol+1] = jcol ;
jcol += 2 ;
}
}
} else {
for ( jcol = 0 ; jcol < nD ; jcol++ ) {
firstlocs[jcol] = 0 ;
sizes[jcol] = jcol ;
}
}
Chv_copyEntriesToVector(frontJ, npivot, pivotsizes, nent, entries,
CHV_STRICT_UPPER_11, CHV_BY_COLUMNS) ;
frontmtx->p_mtxUJJ[J] = mtx ;
SubMtx_columnIndices(mtx, &ncol, &colind) ;
IVcopy(ncol, colind, colindJ) ;
SubMtx_rowIndices(mtx, &nrow, &rowind) ;
IVcopy(nrow, rowind, rowindJ) ;
if ( FRONTMTX_IS_PIVOTING(frontmtx) ) {
frontmtx->nentU += nent ;
}
if ( msglvl > 2 ) {
fprintf(msgFile, "\n\n UJJ matrix") ;
SubMtx_writeForHumanEye(mtx, msgFile) ;
fflush(msgFile) ;
}
}
} else {
/*
------------
sparse front
------------
*/
nent = Chv_countBigEntries(frontJ, npivot, pivotsizes,
CHV_STRICT_UPPER_11, droptol) ;
if ( nent > 0 ) {
nbytes = SubMtx_nbytesNeeded(frontJ->type, SUBMTX_SPARSE_COLUMNS,
nD, nD, nent) ;
if ( msglvl > 2 ) {
fprintf(msgFile,
"\n U_{%d,%d}, nD %d, nent %d, nbytes %d",
J, J, nD, nent, nbytes) ;
fflush(msgFile) ;
}
mtx = SubMtxManager_newObjectOfSizeNbytes(manager, nbytes) ;
SubMtx_init(mtx, frontJ->type, SUBMTX_SPARSE_COLUMNS,
J, J, nD, nD, nent) ;
SubMtx_sparseColumnsInfo(mtx, &ncol, &nent,
&sizes, &indices, &entries) ;
Chv_copyBigEntriesToVector(frontJ, npivot, pivotsizes, sizes,
indices, entries, CHV_STRICT_UPPER_11,
CHV_BY_COLUMNS, droptol) ;
frontmtx->p_mtxUJJ[J] = mtx ;
SubMtx_columnIndices(mtx, &ncol, &colind) ;
IVcopy(ncol, colind, colindJ) ;
SubMtx_rowIndices(mtx, &nrow, &rowind) ;
IVcopy(nrow, rowind, rowindJ) ;
frontmtx->nentU += nent ;
if ( msglvl > 2 ) {
fprintf(msgFile, "\n\n UJJ matrix") ;
SubMtx_writeForHumanEye(mtx, msgFile) ;
fflush(msgFile) ;
}
}
}
/*
-----------------------------------
store the upper U_{J,bnd{J}} matrix
-----------------------------------
*/
if ( FRONTMTX_IS_DENSE_FRONTS(frontmtx) ) {
/*
-----------
dense front
-----------
*/
nent = Chv_countEntries(frontJ, npivot, pivotsizes, CHV_UPPER_12) ;
if ( nent > 0 ) {
mtx = frontmtx->p_mtxUJN[J] ;
if ( mtx == NULL ) {
nbytes = SubMtx_nbytesNeeded(frontJ->type,
SUBMTX_DENSE_COLUMNS, nD, nU, nent) ;
mtx = SubMtxManager_newObjectOfSizeNbytes(manager, nbytes) ;
SubMtx_init(mtx, frontJ->type, SUBMTX_DENSE_COLUMNS,
J, nfront, nD, nU, nent) ;
}
SubMtx_denseInfo(mtx, &nrow, &ncol, &inc1, &inc2, &entries) ;
Chv_copyEntriesToVector(frontJ, npivot, pivotsizes, nent, entries,
CHV_UPPER_12, CHV_BY_COLUMNS) ;
frontmtx->p_mtxUJN[J] = mtx ;
SubMtx_columnIndices(mtx, &ncol, &colind) ;
IVcopy(ncol, colind, colindJ + nD) ;
SubMtx_rowIndices(mtx, &nrow, &rowind) ;
IVcopy(nrow, rowind, rowindJ) ;
if ( FRONTMTX_IS_PIVOTING(frontmtx) ) {
frontmtx->nentU += nent ;
}
if ( msglvl > 2 ) {
fprintf(msgFile, "\n\n UJJN matrix") ;
SubMtx_writeForHumanEye(mtx, msgFile) ;
fflush(msgFile) ;
}
}
} else {
/*
------------
sparse front
------------
*/
nent = Chv_countBigEntries(frontJ, npivot, pivotsizes,
CHV_UPPER_12, droptol) ;
if ( nent > 0 ) {
nbytes = SubMtx_nbytesNeeded(frontJ->type, SUBMTX_SPARSE_COLUMNS,
nD, nU, nent) ;
if ( msglvl > 2 ) {
fprintf(msgFile,
"\n U_{%d,%d}, nD %d, nU %d, nent %d, nbytes %d",
J, nfront, nD, nU, nent, nbytes) ;
fflush(msgFile) ;
}
mtx = SubMtxManager_newObjectOfSizeNbytes(manager, nbytes) ;
SubMtx_init(mtx, frontJ->type, SUBMTX_SPARSE_COLUMNS,
J, nfront, nD, nU, nent) ;
SubMtx_sparseColumnsInfo(mtx,
&nrow, &nent, &sizes, &indices, &entries) ;
Chv_copyBigEntriesToVector(frontJ, npivot, pivotsizes, sizes,
indices, entries, CHV_UPPER_12,
CHV_BY_COLUMNS, droptol) ;
frontmtx->p_mtxUJN[J] = mtx ;
SubMtx_columnIndices(mtx, &ncol, &colind) ;
IVcopy(ncol, colind, colindJ + nD) ;
SubMtx_rowIndices(mtx, &nrow, &rowind) ;
IVcopy(nrow, rowind, rowindJ) ;
frontmtx->nentU += nent ;
if ( msglvl > 2 ) {
fprintf(msgFile, "\n\n UJJN matrix") ;
SubMtx_writeForHumanEye(mtx, msgFile) ;
fflush(msgFile) ;
}
}
}
if ( FRONTMTX_IS_NONSYMMETRIC(frontmtx) ) {
/*
------------------------------
store the lower L_{J,J} matrix
------------------------------
*/
if ( FRONTMTX_IS_DENSE_FRONTS(frontmtx) ) {
/*
-----------
dense front
-----------
*/
nent = Chv_countEntries(frontJ, npivot, pivotsizes,
CHV_STRICT_LOWER_11) ;
if ( nent > 0 ) {
mtx = frontmtx->p_mtxLJJ[J] ;
if ( mtx == NULL ) {
nbytes = SubMtx_nbytesNeeded(frontJ->type,
SUBMTX_DENSE_SUBROWS, nD, nD,nent);
mtx = SubMtxManager_newObjectOfSizeNbytes(manager, nbytes) ;
SubMtx_init(mtx, frontJ->type, SUBMTX_DENSE_SUBROWS,
J, J, nD, nD, nent) ;
}
SubMtx_denseSubrowsInfo(mtx, &nrow, &nent,
&firstlocs, &sizes, &entries) ;
for ( irow = 0 ; irow < nD ; irow++ ) {
firstlocs[irow] = 0 ;
sizes[irow] = irow ;
}
Chv_copyEntriesToVector(frontJ, npivot, pivotsizes, nent,
entries, CHV_STRICT_LOWER_11, CHV_BY_ROWS) ;
frontmtx->p_mtxLJJ[J] = mtx ;
SubMtx_columnIndices(mtx, &ncol, &colind) ;
IVcopy(ncol, colind, colindJ) ;
SubMtx_rowIndices(mtx, &nrow, &rowind) ;
IVcopy(nrow, rowind, rowindJ) ;
if ( FRONTMTX_IS_PIVOTING(frontmtx) ) {
frontmtx->nentL += nent ;
}
if ( msglvl > 2 ) {
fprintf(msgFile, "\n\n LJJ matrix") ;
SubMtx_writeForHumanEye(mtx, msgFile) ;
fflush(msgFile) ;
}
}
} else {
/*
------------
sparse front
------------
*/
nent = Chv_countBigEntries(frontJ, npivot, pivotsizes,
CHV_STRICT_LOWER_11, droptol);
if ( nent > 0 ) {
nbytes = SubMtx_nbytesNeeded(frontJ->type, SUBMTX_SPARSE_ROWS,
nD, nD, nent) ;
if ( msglvl > 2 ) {
fprintf(msgFile,
"\n L_{%d,%d}, nD %d, nent %d, nbytes %d",
J, J, nD, nent, nbytes) ;
fflush(msgFile) ;
}
mtx = SubMtxManager_newObjectOfSizeNbytes(manager, nbytes) ;
SubMtx_init(mtx, frontJ->type, SUBMTX_SPARSE_ROWS,
J, J, nD, nD, nent) ;
SubMtx_sparseRowsInfo(mtx,
&nrow, &nent, &sizes, &indices, &entries);
Chv_copyBigEntriesToVector(frontJ, npivot, pivotsizes, sizes,
indices, entries, CHV_STRICT_LOWER_11,
CHV_BY_ROWS, droptol) ;
frontmtx->p_mtxLJJ[J] = mtx ;
SubMtx_columnIndices(mtx, &ncol, &colind) ;
IVcopy(ncol, colind, colindJ) ;
SubMtx_rowIndices(mtx, &nrow, &rowind) ;
IVcopy(nrow, rowind, rowindJ) ;
frontmtx->nentL += nent ;
if ( msglvl > 2 ) {
fprintf(msgFile, "\n\n LJJ matrix") ;
SubMtx_writeForHumanEye(mtx, msgFile) ;
fflush(msgFile) ;
}
}
}
/*
-----------------------------------
store the lower L_{bnd{J},J} matrix
-----------------------------------
*/
if ( FRONTMTX_IS_DENSE_FRONTS(frontmtx) ) {
/*
-----------
dense front
-----------
*/
nent = Chv_countEntries(frontJ, npivot, pivotsizes, CHV_LOWER_21);
if ( nent > 0 ) {
mtx = frontmtx->p_mtxLNJ[J] ;
if ( mtx == NULL ) {
nbytes = SubMtx_nbytesNeeded(frontJ->type,
SUBMTX_DENSE_ROWS, nL, nD, nent) ;
mtx = SubMtxManager_newObjectOfSizeNbytes(manager, nbytes) ;
SubMtx_init(mtx, frontJ->type, SUBMTX_DENSE_ROWS,
nfront, J, nL, nD, nent) ;
}
SubMtx_denseInfo(mtx, &nrow, &ncol, &inc1, &inc2, &entries) ;
Chv_copyEntriesToVector(frontJ, npivot, pivotsizes, nent,
entries, CHV_LOWER_21, CHV_BY_ROWS) ;
frontmtx->p_mtxLNJ[J] = mtx ;
SubMtx_columnIndices(mtx, &ncol, &colind) ;
IVcopy(ncol, colind, colindJ) ;
SubMtx_rowIndices(mtx, &nrow, &rowind) ;
IVcopy(nrow, rowind, rowindJ + nD) ;
if ( FRONTMTX_IS_PIVOTING(frontmtx) ) {
frontmtx->nentL += nent ;
}
if ( msglvl > 2 ) {
fprintf(msgFile, "\n\n LNJ matrix") ;
SubMtx_writeForHumanEye(mtx, msgFile) ;
fflush(msgFile) ;
}
}
} else {
/*
------------
sparse front
------------
*/
nent = Chv_countBigEntries(frontJ, npivot, pivotsizes,
CHV_LOWER_21, droptol) ;
if ( nent > 0 ) {
nbytes = SubMtx_nbytesNeeded(frontJ->type, SUBMTX_SPARSE_ROWS,
nL, nD, nent) ;
if ( msglvl > 2 ) {
fprintf(msgFile,
"\n L_{%d,%d}, nL %d, nD %d, nent %d, nbytes %d",
nfront, J, nL, nD, nent, nbytes) ;
fflush(msgFile) ;
}
mtx = SubMtxManager_newObjectOfSizeNbytes(manager, nbytes) ;
SubMtx_init(mtx, frontJ->type, SUBMTX_SPARSE_ROWS,
nfront, J, nL, nD, nent) ;
SubMtx_sparseRowsInfo(mtx,
&nrow, &nent, &sizes, &indices, &entries) ;
Chv_copyBigEntriesToVector(frontJ, npivot, pivotsizes, sizes,
indices, entries, CHV_LOWER_21,
CHV_BY_ROWS, droptol) ;
frontmtx->p_mtxLNJ[J] = mtx ;
SubMtx_columnIndices(mtx, &ncol, &colind) ;
IVcopy(ncol, colind, colindJ) ;
SubMtx_rowIndices(mtx, &nrow, &rowind) ;
IVcopy(nrow, rowind, rowindJ + nD) ;
frontmtx->nentL += nent ;
if ( msglvl > 2 ) {
fprintf(msgFile, "\n\n LNJ matrix") ;
SubMtx_writeForHumanEye(mtx, msgFile) ;
fflush(msgFile) ;
}
}
}
}
return ; }
/*--------------------------------------------------------------------*/
syntax highlighted by Code2HTML, v. 0.9.1