/* util.c */
#include "../SubMtx.h"
#define MYDEBUG 0
/*--------------------------------------------------------------------*/
/*
-----------------------------------------
purpose -- to fill rowDV with the entries
in row irow of the matrix
created -- 98may01, cca
-----------------------------------------
*/
void
SubMtx_fillRowDV (
SubMtx *mtx,
int irow,
DV *rowDV
) {
double *rowvec ;
/*
---------------
check the input
---------------
*/
if ( mtx == NULL || irow < 0 || rowDV == NULL ) {
fprintf(stderr, "\n fatal error in SubMtx_fillRowDV(%p,%d,%p)"
"\n bad input\n", mtx, irow, rowDV) ;
exit(-1) ;
}
if ( ! SUBMTX_IS_REAL(mtx) ) {
fprintf(stderr, "\n fatal error in SubMtx_fillRowDV(%p,%d,%p)"
"\n type = %d, must be SPOOLES_REAL\n",
mtx, irow, rowDV, mtx->type) ;
exit(-1) ;
}
DV_setSize(rowDV, mtx->ncol) ;
rowvec = DV_entries(rowDV) ;
DVzero(mtx->ncol, rowvec) ;
/*
--------------------------------------
switch over the different matrix types
--------------------------------------
*/
switch ( mtx->mode ) {
case SUBMTX_DENSE_ROWS :
case SUBMTX_DENSE_COLUMNS : {
double *entries ;
int inc1, inc2, jcol, loc, ncol, nrow ;
SubMtx_denseInfo(mtx, &nrow, &ncol, &inc1, &inc2, &entries) ;
for ( jcol = 0 ; jcol < ncol ; jcol++ ) {
loc = irow*inc1 + jcol*inc2 ;
rowvec[jcol] = entries[loc] ;
}
} break ;
case SUBMTX_SPARSE_ROWS : {
double *entries ;
int ii, jrow, kk, nent, nrow, offset ;
int *indices, *sizes ;
SubMtx_sparseRowsInfo(mtx,
&nrow, &nent, &sizes, &indices, &entries) ;
for ( jrow = offset = 0 ; jrow < irow ; jrow++ ) {
offset += sizes[jrow] ;
}
for ( ii = 0, kk = offset ; ii < sizes[irow] ; ii++, kk++ ) {
rowvec[indices[kk]] = entries[kk] ;
}
} break ;
case SUBMTX_SPARSE_COLUMNS : {
double *entries ;
int ii, jcol, kk, nent, ncol, offset ;
int *indices, *sizes ;
SubMtx_sparseColumnsInfo(mtx, &ncol, &nent,
&sizes, &indices, &entries) ;
for ( jcol = offset = 0 ; jcol < ncol ; jcol++ ) {
for ( ii = 0, kk = offset ; ii < sizes[jcol] ; ii++, kk++ ) {
if ( indices[kk] == irow ) {
rowvec[jcol] = entries[kk] ;
break ;
}
}
offset += sizes[jcol] ;
}
} break ;
case SUBMTX_SPARSE_TRIPLES : {
double *entries ;
int ii, nent ;
int *colids, *rowids ;
SubMtx_sparseTriplesInfo(mtx, &nent, &rowids, &colids, &entries) ;
for ( ii = 0 ; ii < nent ; ii++ ) {
if ( rowids[ii] == irow ) {
rowvec[colids[ii]] = entries[ii] ;
}
}
} break ;
case SUBMTX_DENSE_SUBROWS : {
double *entries ;
int first, ii, jrow, kk, last, nent, nrow, offset ;
int *firstlocs, *sizes ;
SubMtx_denseSubrowsInfo(mtx, &nrow, &nent,
&firstlocs, &sizes, &entries) ;
#if MYDEBUG > 0
fprintf(stdout, "\n %% entries(%d) :", nent) ;
DVfprintf(stdout, nent, entries) ;
#endif
for ( jrow = offset = 0 ; jrow < irow ; jrow++ ) {
offset += sizes[jrow] ;
}
#if MYDEBUG > 0
fprintf(stdout, "\n %% irow = %d, offset = %d", irow, offset) ;
fprintf(stdout, "\n %% first = %d, size = %d",
firstlocs[irow], sizes[irow]) ;
#endif
if ( sizes[irow] > 0 ) {
first = firstlocs[irow] ;
last = first + sizes[irow] - 1 ;
for ( kk = first, ii = offset ; kk <= last ; kk++, ii++ ) {
rowvec[kk] = entries[ii] ;
#if MYDEBUG > 0
fprintf(stdout, "\n %% rowvec[%d] = entries[%d]", kk, ii) ;
#endif
}
}
} break ;
case SUBMTX_DENSE_SUBCOLUMNS : {
double *entries ;
int first, jcol, last, loc, nent, ncol, offset ;
int *firstlocs, *sizes ;
SubMtx_denseSubcolumnsInfo(mtx, &ncol, &nent,
&firstlocs, &sizes, &entries) ;
for ( jcol = offset = 0 ; jcol < ncol ; jcol++ ) {
if ( sizes[jcol] > 0 ) {
first = firstlocs[jcol] ;
last = first + sizes[jcol] - 1 ;
if ( first <= irow && irow <= last ) {
loc = offset + irow - first ;
rowvec[jcol] = entries[loc] ;
}
offset += sizes[jcol] ;
}
}
} break ;
case SUBMTX_DIAGONAL : {
double *entries ;
int nent ;
SubMtx_diagonalInfo(mtx, &nent, &entries) ;
rowvec[irow] = entries[irow] ;
} break ;
case SUBMTX_BLOCK_DIAGONAL_SYM : {
double *entries ;
int ii, ipivot, jrow, kk, m, nrow, nent, stride ;
int *pivotsizes ;
SubMtx_blockDiagonalInfo(mtx, &nrow, &nent, &pivotsizes, &entries) ;
for ( jrow = ipivot = kk = 0 ; jrow <= irow ; ipivot++ ) {
m = pivotsizes[ipivot] ;
/*
fprintf(stdout, "\n jrow %d, m %d, kk %d", jrow, m, kk) ;
fflush(stdout) ;
*/
if ( jrow <= irow && irow < jrow + m ) {
stride = m - 1 ;
kk += irow - jrow ;
/*
fprintf(stdout, "\n 0. kk %d, stride %d", kk, stride) ;
fflush(stdout) ;
*/
for ( ii = jrow ; ii <= irow ; ii++ ) {
/*
fprintf(stdout, "\n 1. ii %d, kk %d", ii, kk) ;
fflush(stdout) ;
*/
rowvec[ii] = entries[kk] ;
kk += stride-- ;
}
for ( ; ii < jrow + m ; ii++ ) {
/*
fprintf(stdout, "\n 2. ii %d, kk %d", ii, kk) ;
fflush(stdout) ;
*/
rowvec[ii] = entries[kk] ;
kk++ ;
}
} else {
kk += (m*(m+1))/2 ;
}
jrow += m ;
}
} break ;
case SUBMTX_BLOCK_DIAGONAL_HERM : {
double *entries ;
int ii, ipivot, jrow, kk, m, nrow, nent, stride ;
int *pivotsizes ;
SubMtx_blockDiagonalInfo(mtx, &nrow, &nent, &pivotsizes, &entries) ;
for ( jrow = ipivot = kk = 0 ; jrow <= irow ; ipivot++ ) {
m = pivotsizes[ipivot] ;
/*
fprintf(stdout, "\n jrow %d, m %d, kk %d", jrow, m, kk) ;
fflush(stdout) ;
*/
if ( jrow <= irow && irow < jrow + m ) {
stride = m - 1 ;
kk += irow - jrow ;
/*
fprintf(stdout, "\n 0. kk %d, stride %d", kk, stride) ;
fflush(stdout) ;
*/
for ( ii = jrow ; ii < irow ; ii++ ) {
/*
fprintf(stdout, "\n 1. ii %d, kk %d", ii, kk) ;
fflush(stdout) ;
*/
rowvec[ii] = entries[kk] ;
kk += stride-- ;
}
for ( ; ii < jrow + m ; ii++ ) {
/*
fprintf(stdout, "\n 2. ii %d, kk %d", ii, kk) ;
fflush(stdout) ;
*/
rowvec[ii] = entries[kk] ;
kk++ ;
}
} else {
kk += (m*(m+1))/2 ;
}
jrow += m ;
}
} break ;
default :
break ;
}
return ; }
/*--------------------------------------------------------------------*/
/*
-----------------------------------------
purpose -- to fill rowZV with the entries
in row irow of the matrix
created -- 98may01, cca
-----------------------------------------
*/
void
SubMtx_fillRowZV (
SubMtx *mtx,
int irow,
ZV *rowZV
) {
double *rowvec ;
/*
---------------
check the input
---------------
*/
if ( mtx == NULL || irow < 0 || rowZV == NULL ) {
fprintf(stderr, "\n fatal error in SubMtx_fillRowZV(%p,%d,%p)"
"\n bad input\n", mtx, irow, rowZV) ;
exit(-1) ;
}
if ( ! SUBMTX_IS_COMPLEX(mtx) ) {
fprintf(stderr, "\n fatal error in SubMtx_fillRowZV(%p,%d,%p)"
"\n type = %d, must be SPOOLES_COMPLEX\n",
mtx, irow, rowZV, mtx->type) ;
exit(-1) ;
}
ZV_setSize(rowZV, mtx->ncol) ;
rowvec = ZV_entries(rowZV) ;
DVzero(2*mtx->ncol, rowvec) ;
/*
--------------------------------------
switch over the different matrix types
--------------------------------------
*/
switch ( mtx->mode ) {
case SUBMTX_DENSE_ROWS :
case SUBMTX_DENSE_COLUMNS : {
double *entries ;
int inc1, inc2, jcol, loc, ncol, nrow ;
SubMtx_denseInfo(mtx, &nrow, &ncol, &inc1, &inc2, &entries) ;
for ( jcol = 0 ; jcol < ncol ; jcol++ ) {
loc = irow*inc1 + jcol*inc2 ;
rowvec[2*jcol] = entries[2*loc] ;
rowvec[2*jcol+1] = entries[2*loc+1] ;
}
} break ;
case SUBMTX_SPARSE_ROWS : {
double *entries ;
int ii, jrow, kk, nent, nrow, offset ;
int *indices, *sizes ;
SubMtx_sparseRowsInfo(mtx, &nrow, &nent, &sizes, &indices, &entries) ;
for ( jrow = offset = 0 ; jrow < irow ; jrow++ ) {
offset += sizes[jrow] ;
}
for ( ii = 0, kk = offset ; ii < sizes[irow] ; ii++, kk++ ) {
rowvec[2*indices[kk]] = entries[2*kk] ;
rowvec[2*indices[kk]+1] = entries[2*kk+1] ;
}
} break ;
case SUBMTX_SPARSE_COLUMNS : {
double *entries ;
int ii, jcol, kk, nent, ncol, offset ;
int *indices, *sizes ;
SubMtx_sparseColumnsInfo(mtx, &ncol, &nent,
&sizes, &indices, &entries) ;
for ( jcol = offset = 0 ; jcol < ncol ; jcol++ ) {
for ( ii = 0, kk = offset ; ii < sizes[jcol] ; ii++, kk++ ) {
if ( indices[kk] == irow ) {
rowvec[2*jcol] = entries[2*kk] ;
rowvec[2*jcol+1] = entries[2*kk+1] ;
break ;
}
}
offset += sizes[jcol] ;
}
} break ;
case SUBMTX_SPARSE_TRIPLES : {
double *entries ;
int ii, nent ;
int *colids, *rowids ;
SubMtx_sparseTriplesInfo(mtx, &nent, &rowids, &colids, &entries) ;
for ( ii = 0 ; ii < nent ; ii++ ) {
if ( rowids[ii] == irow ) {
rowvec[2*colids[ii]] = entries[2*ii] ;
rowvec[2*colids[ii]+1] = entries[2*ii+1] ;
}
}
} break ;
case SUBMTX_DENSE_SUBROWS : {
double *entries ;
int first, ii, jrow, kk, last, nent, nrow, offset ;
int *firstlocs, *sizes ;
SubMtx_denseSubrowsInfo(mtx, &nrow, &nent,
&firstlocs, &sizes, &entries) ;
#if MYDEBUG > 0
fprintf(stdout, "\n %% entries(%d) :", nent) ;
ZVfprintf(stdout, nent, entries) ;
#endif
for ( jrow = offset = 0 ; jrow < irow ; jrow++ ) {
offset += sizes[jrow] ;
}
#if MYDEBUG > 0
fprintf(stdout, "\n %% irow = %d, offset = %d", irow, offset) ;
fprintf(stdout, "\n %% first = %d, size = %d",
firstlocs[irow], sizes[irow]) ;
#endif
if ( sizes[irow] > 0 ) {
first = firstlocs[irow] ;
last = first + sizes[irow] - 1 ;
for ( kk = first, ii = offset ; kk <= last ; kk++, ii++ ) {
rowvec[2*kk] = entries[2*ii] ;
rowvec[2*kk+1] = entries[2*ii+1] ;
#if MYDEBUG > 0
fprintf(stdout, "\n %% rowvec[%d] = entries[%d]", kk, ii) ;
#endif
}
}
} break ;
case SUBMTX_DENSE_SUBCOLUMNS : {
double *entries ;
int first, jcol, last, loc, nent, ncol, offset ;
int *firstlocs, *sizes ;
SubMtx_denseSubcolumnsInfo(mtx, &ncol, &nent,
&firstlocs, &sizes, &entries) ;
for ( jcol = offset = 0 ; jcol < ncol ; jcol++ ) {
if ( sizes[jcol] > 0 ) {
first = firstlocs[jcol] ;
last = first + sizes[jcol] - 1 ;
if ( first <= irow && irow <= last ) {
loc = offset + irow - first ;
rowvec[2*jcol] = entries[2*loc] ;
rowvec[2*jcol+1] = entries[2*loc+1] ;
}
offset += sizes[jcol] ;
}
}
} break ;
case SUBMTX_DIAGONAL : {
double *entries ;
int nent ;
SubMtx_diagonalInfo(mtx, &nent, &entries) ;
rowvec[2*irow] = entries[2*irow] ;
rowvec[2*irow+1] = entries[2*irow+1] ;
} break ;
case SUBMTX_BLOCK_DIAGONAL_SYM : {
double *entries ;
int ii, ipivot, jrow, kk, m, nrow, nent, stride ;
int *pivotsizes ;
SubMtx_blockDiagonalInfo(mtx, &nrow, &nent, &pivotsizes, &entries) ;
for ( jrow = ipivot = kk = 0 ; jrow <= irow ; ipivot++ ) {
m = pivotsizes[ipivot] ;
/*
fprintf(stdout, "\n jrow %d, m %d, kk %d", jrow, m, kk) ;
fflush(stdout) ;
*/
if ( jrow <= irow && irow < jrow + m ) {
stride = m - 1 ;
kk += irow - jrow ;
/*
fprintf(stdout, "\n 0. kk %d, stride %d", kk, stride) ;
fflush(stdout) ;
*/
for ( ii = jrow ; ii <= irow ; ii++ ) {
/*
fprintf(stdout, "\n 1. ii %d, kk %d", ii, kk) ;
fflush(stdout) ;
*/
rowvec[2*ii] = entries[2*kk] ;
rowvec[2*ii+1] = entries[2*kk+1] ;
kk += stride-- ;
}
for ( ; ii < jrow + m ; ii++ ) {
/*
fprintf(stdout, "\n 2. ii %d, kk %d", ii, kk) ;
fflush(stdout) ;
*/
rowvec[2*ii] = entries[2*kk] ;
rowvec[2*ii+1] = entries[2*kk+1] ;
kk++ ;
}
} else {
kk += (m*(m+1))/2 ;
}
jrow += m ;
}
} break ;
case SUBMTX_BLOCK_DIAGONAL_HERM : {
double *entries ;
int ii, ipivot, jrow, kk, m, nrow, nent, stride ;
int *pivotsizes ;
SubMtx_blockDiagonalInfo(mtx, &nrow, &nent, &pivotsizes, &entries) ;
for ( jrow = ipivot = kk = 0 ; jrow <= irow ; ipivot++ ) {
m = pivotsizes[ipivot] ;
/*
fprintf(stdout, "\n jrow %d, m %d, kk %d", jrow, m, kk) ;
fflush(stdout) ;
*/
if ( jrow <= irow && irow < jrow + m ) {
stride = m - 1 ;
kk += irow - jrow ;
/*
fprintf(stdout, "\n 0. kk %d, stride %d", kk, stride) ;
fflush(stdout) ;
*/
for ( ii = jrow ; ii < irow ; ii++ ) {
/*
fprintf(stdout, "\n 1. ii %d, kk %d", ii, kk) ;
fflush(stdout) ;
*/
rowvec[2*ii] = entries[2*kk] ;
rowvec[2*ii+1] = -entries[2*kk+1] ;
kk += stride-- ;
}
for ( ; ii < jrow + m ; ii++ ) {
/*
fprintf(stdout, "\n 2. ii %d, kk %d", ii, kk) ;
fflush(stdout) ;
*/
rowvec[2*ii] = entries[2*kk] ;
rowvec[2*ii+1] = entries[2*kk+1] ;
kk++ ;
}
} else {
kk += (m*(m+1))/2 ;
}
jrow += m ;
}
} break ;
default :
break ;
}
return ; }
/*--------------------------------------------------------------------*/
/*
-----------------------------------------
purpose -- to fill colDV with the entries
in column icol of the matrix
created -- 98may01, cca
-----------------------------------------
*/
void
SubMtx_fillColumnDV (
SubMtx *mtx,
int icol,
DV *colDV
) {
double *colvec ;
/*
---------------
check the input
---------------
*/
if ( mtx == NULL || icol < 0 || colDV == NULL ) {
fprintf(stderr, "\n fatal error in SubMtx_fillColumnDV(%p,%d,%p)"
"\n bad input\n", mtx, icol, colDV) ;
exit(-1) ;
}
if ( ! SUBMTX_IS_REAL(mtx) ) {
fprintf(stderr, "\n fatal error in SubMtx_fillColumnDV(%p,%d,%p)"
"\n bad type %d, must be SPOOLES_REAL\n",
mtx, icol, colDV, mtx->type) ;
exit(-1) ;
}
DV_setSize(colDV, mtx->nrow) ;
colvec = DV_entries(colDV) ;
DVzero(mtx->nrow, colvec) ;
/*
--------------------------------------
switch over the different matrix types
--------------------------------------
*/
switch ( mtx->mode ) {
case SUBMTX_DENSE_ROWS :
case SUBMTX_DENSE_COLUMNS : {
double *entries ;
int inc1, inc2, jrow, loc, ncol, nrow ;
SubMtx_denseInfo(mtx, &nrow, &ncol, &inc1, &inc2, &entries) ;
for ( jrow = 0 ; jrow < nrow ; jrow++ ) {
loc = jrow*inc1 + icol*inc2 ;
colvec[jrow] = entries[loc] ;
}
} break ;
case SUBMTX_SPARSE_COLUMNS : {
double *entries ;
int ii, jcol, kk, ncol, nent, offset ;
int *indices, *sizes ;
SubMtx_sparseColumnsInfo(mtx, &ncol, &nent,
&sizes, &indices, &entries) ;
for ( jcol = offset = 0 ; jcol < icol ; jcol++ ) {
offset += sizes[jcol] ;
}
for ( ii = 0, kk = offset ; ii < sizes[icol] ; ii++, kk++ ) {
colvec[indices[kk]] = entries[kk] ;
}
} break ;
case SUBMTX_SPARSE_ROWS : {
double *entries ;
int ii, jrow, kk, nent, nrow, offset ;
int *indices, *sizes ;
SubMtx_sparseRowsInfo(mtx, &nrow, &nent, &sizes, &indices, &entries) ;
for ( jrow = offset = 0 ; jrow < nrow ; jrow++ ) {
for ( ii = 0, kk = offset ; ii < sizes[jrow] ; ii++, kk++ ) {
if ( indices[kk] == icol ) {
colvec[jrow] = entries[kk] ;
break ;
}
}
offset += sizes[jrow] ;
}
} break ;
case SUBMTX_SPARSE_TRIPLES : {
double *entries ;
int ii, nent ;
int *colids, *rowids ;
SubMtx_sparseTriplesInfo(mtx, &nent, &rowids, &colids, &entries) ;
for ( ii = 0 ; ii < nent ; ii++ ) {
if ( colids[ii] == icol ) {
colvec[rowids[ii]] = entries[ii] ;
}
}
} break ;
case SUBMTX_DENSE_SUBCOLUMNS : {
double *entries ;
int first, ii, jcol, kk, last, ncol, nent, offset ;
int *firstlocs, *sizes ;
SubMtx_denseSubcolumnsInfo(mtx, &ncol, &nent,
&firstlocs, &sizes, &entries) ;
#if MYDEBUG > 0
fprintf(stdout, "\n %% entries(%d) :", nent) ;
DVfprintf(stdout, nent, entries) ;
#endif
for ( jcol = offset = 0 ; jcol < icol ; jcol++ ) {
offset += sizes[jcol] ;
}
#if MYDEBUG > 0
fprintf(stdout, "\n %% icol = %d, offset = %d", icol, offset) ;
fprintf(stdout, "\n %% first = %d, size = %d",
firstlocs[icol], sizes[icol]) ;
#endif
if ( sizes[icol] > 0 ) {
first = firstlocs[icol] ;
last = first + sizes[icol] - 1 ;
for ( kk = first, ii = offset ; kk <= last ; kk++, ii++ ) {
colvec[kk] = entries[ii] ;
#if MYDEBUG > 0
fprintf(stdout, "\n %% colvec[%d] = entries[%d]", kk, ii) ;
#endif
}
}
} break ;
case SUBMTX_DENSE_SUBROWS : {
double *entries ;
int first, jrow, last, loc, nent, nrow, offset ;
int *firstlocs, *sizes ;
SubMtx_denseSubrowsInfo(mtx, &nrow, &nent,
&firstlocs, &sizes, &entries) ;
for ( jrow = offset = 0 ; jrow < nrow ; jrow++ ) {
if ( sizes[jrow] > 0 ) {
first = firstlocs[jrow] ;
last = first + sizes[jrow] - 1 ;
if ( first <= icol && icol <= last ) {
loc = offset + icol - first ;
colvec[jrow] = entries[loc] ;
}
offset += sizes[jrow] ;
}
}
} break ;
case SUBMTX_DIAGONAL : {
double *entries ;
int nent ;
SubMtx_diagonalInfo(mtx, &nent, &entries) ;
colvec[icol] = entries[icol] ;
} break ;
case SUBMTX_BLOCK_DIAGONAL_SYM : {
double *entries ;
int ii, ipivot, jcol, kk, m, nrow, nent, stride ;
int *pivotsizes ;
SubMtx_blockDiagonalInfo(mtx, &nrow, &nent, &pivotsizes, &entries) ;
for ( jcol = ipivot = kk = 0 ; jcol <= icol ; ipivot++ ) {
m = pivotsizes[ipivot] ;
/*
fprintf(stdout, "\n jrow %d, m %d, kk %d", jrow, m, kk) ;
fflush(stdout) ;
*/
if ( jcol <= icol && icol < jcol + m ) {
stride = m - 1 ;
kk += icol - jcol ;
/*
fprintf(stdout, "\n 0. kk %d, stride %d", kk, stride) ;
fflush(stdout) ;
*/
for ( ii = jcol ; ii <= icol ; ii++ ) {
/*
fprintf(stdout, "\n 1. ii %d, kk %d", ii, kk) ;
fflush(stdout) ;
*/
colvec[ii] = entries[kk] ;
kk += stride-- ;
}
for ( ; ii < jcol + m ; ii++ ) {
/*
fprintf(stdout, "\n 2. ii %d, kk %d", ii, kk) ;
fflush(stdout) ;
*/
colvec[ii] = entries[kk] ;
kk++ ;
}
} else {
kk += (m*(m+1))/2 ;
}
jcol += m ;
}
} break ;
case SUBMTX_BLOCK_DIAGONAL_HERM : {
double *entries ;
int ii, ipivot, jcol, kk, m, nrow, nent, stride ;
int *pivotsizes ;
SubMtx_blockDiagonalInfo(mtx, &nrow, &nent, &pivotsizes, &entries) ;
for ( jcol = ipivot = kk = 0 ; jcol <= icol ; ipivot++ ) {
m = pivotsizes[ipivot] ;
/*
fprintf(stdout, "\n jrow %d, m %d, kk %d", jrow, m, kk) ;
fflush(stdout) ;
*/
if ( jcol <= icol && icol < jcol + m ) {
stride = m - 1 ;
kk += icol - jcol ;
/*
fprintf(stdout, "\n 0. kk %d, stride %d", kk, stride) ;
fflush(stdout) ;
*/
for ( ii = jcol ; ii <= icol ; ii++ ) {
/*
fprintf(stdout, "\n 1. ii %d, kk %d", ii, kk) ;
fflush(stdout) ;
*/
colvec[ii] = entries[kk] ;
kk += stride-- ;
}
for ( ; ii < jcol + m ; ii++ ) {
/*
fprintf(stdout, "\n 2. ii %d, kk %d", ii, kk) ;
fflush(stdout) ;
*/
colvec[ii] = entries[kk] ;
kk++ ;
}
} else {
kk += (m*(m+1))/2 ;
}
jcol += m ;
}
} break ;
default :
break ;
}
return ; }
/*--------------------------------------------------------------------*/
/*
-----------------------------------------
purpose -- to fill colZV with the entries
in column icol of the matrix
created -- 98may01, cca
-----------------------------------------
*/
void
SubMtx_fillColumnZV (
SubMtx *mtx,
int icol,
ZV *colZV
) {
double *colvec ;
/*
---------------
check the input
---------------
*/
if ( mtx == NULL || icol < 0 || colZV == NULL ) {
fprintf(stderr, "\n fatal error in SubMtx_fillColumnZV(%p,%d,%p)"
"\n bad input\n", mtx, icol, colZV) ;
exit(-1) ;
}
if ( ! SUBMTX_IS_COMPLEX(mtx) ) {
fprintf(stderr, "\n fatal error in SubMtx_fillColumnZV(%p,%d,%p)"
"\n bad type %d, must be SPOOLES_COMPLEX\n",
mtx, icol, colZV, mtx->type) ;
exit(-1) ;
}
ZV_setSize(colZV, mtx->nrow) ;
colvec = ZV_entries(colZV) ;
DVzero(2*mtx->nrow, colvec) ;
/*
--------------------------------------
switch over the different matrix types
--------------------------------------
*/
switch ( mtx->mode ) {
case SUBMTX_DENSE_ROWS :
case SUBMTX_DENSE_COLUMNS : {
double *entries ;
int inc1, inc2, jrow, loc, ncol, nrow ;
SubMtx_denseInfo(mtx, &nrow, &ncol, &inc1, &inc2, &entries) ;
for ( jrow = 0 ; jrow < nrow ; jrow++ ) {
loc = jrow*inc1 + icol*inc2 ;
colvec[2*jrow] = entries[2*loc] ;
colvec[2*jrow+1] = entries[2*loc+1] ;
}
} break ;
case SUBMTX_SPARSE_COLUMNS : {
double *entries ;
int ii, jcol, kk, ncol, nent, offset ;
int *indices, *sizes ;
SubMtx_sparseColumnsInfo(mtx, &ncol, &nent,
&sizes, &indices, &entries) ;
for ( jcol = offset = 0 ; jcol < icol ; jcol++ ) {
offset += sizes[jcol] ;
}
for ( ii = 0, kk = offset ; ii < sizes[icol] ; ii++, kk++ ) {
colvec[2*indices[kk]] = entries[2*kk] ;
colvec[2*indices[kk]+1] = entries[2*kk+1] ;
}
} break ;
case SUBMTX_SPARSE_ROWS : {
double *entries ;
int ii, jrow, kk, nent, nrow, offset ;
int *indices, *sizes ;
SubMtx_sparseRowsInfo(mtx, &nrow, &nent, &sizes, &indices, &entries) ;
for ( jrow = offset = 0 ; jrow < nrow ; jrow++ ) {
for ( ii = 0, kk = offset ; ii < sizes[jrow] ; ii++, kk++ ) {
if ( indices[kk] == icol ) {
colvec[2*jrow] = entries[2*kk] ;
colvec[2*jrow+1] = entries[2*kk+1] ;
break ;
}
}
offset += sizes[jrow] ;
}
} break ;
case SUBMTX_SPARSE_TRIPLES : {
double *entries ;
int ii, nent ;
int *colids, *rowids ;
SubMtx_sparseTriplesInfo(mtx, &nent, &rowids, &colids, &entries) ;
for ( ii = 0 ; ii < nent ; ii++ ) {
if ( colids[ii] == icol ) {
colvec[2*rowids[ii]] = entries[2*ii] ;
colvec[2*rowids[ii]+1] = entries[2*ii+1] ;
}
}
} break ;
case SUBMTX_DENSE_SUBCOLUMNS : {
double *entries ;
int first, ii, jcol, kk, last, ncol, nent, offset ;
int *firstlocs, *sizes ;
SubMtx_denseSubcolumnsInfo(mtx, &ncol, &nent,
&firstlocs, &sizes, &entries) ;
#if MYDEBUG > 0
fprintf(stdout, "\n %% entries(%d) :", nent) ;
ZVfprintf(stdout, nent, entries) ;
#endif
for ( jcol = offset = 0 ; jcol < icol ; jcol++ ) {
offset += sizes[jcol] ;
}
#if MYDEBUG > 0
fprintf(stdout, "\n %% icol = %d, offset = %d", icol, offset) ;
fprintf(stdout, "\n %% first = %d, size = %d",
firstlocs[icol], sizes[icol]) ;
#endif
if ( sizes[icol] > 0 ) {
first = firstlocs[icol] ;
last = first + sizes[icol] - 1 ;
for ( kk = first, ii = offset ; kk <= last ; kk++, ii++ ) {
colvec[2*kk] = entries[2*ii] ;
colvec[2*kk+1] = entries[2*ii+1] ;
#if MYDEBUG > 0
fprintf(stdout, "\n %% colvec[%d] = entries[%d]", kk, ii) ;
#endif
}
}
} break ;
case SUBMTX_DENSE_SUBROWS : {
double *entries ;
int first, jrow, last, loc, nent, nrow, offset ;
int *firstlocs, *sizes ;
SubMtx_denseSubrowsInfo(mtx, &nrow, &nent,
&firstlocs, &sizes, &entries) ;
for ( jrow = offset = 0 ; jrow < nrow ; jrow++ ) {
if ( sizes[jrow] > 0 ) {
first = firstlocs[jrow] ;
last = first + sizes[jrow] - 1 ;
if ( first <= icol && icol <= last ) {
loc = offset + icol - first ;
colvec[2*jrow] = entries[2*loc] ;
colvec[2*jrow+1] = entries[2*loc+1] ;
}
offset += sizes[jrow] ;
}
}
} break ;
case SUBMTX_DIAGONAL : {
double *entries ;
int nent ;
SubMtx_diagonalInfo(mtx, &nent, &entries) ;
colvec[2*icol] = entries[2*icol] ;
colvec[2*icol+1] = entries[2*icol+1] ;
} break ;
case SUBMTX_BLOCK_DIAGONAL_SYM : {
double *entries ;
int ii, ipivot, jcol, kk, m, nrow, nent, stride ;
int *pivotsizes ;
SubMtx_blockDiagonalInfo(mtx, &nrow, &nent, &pivotsizes, &entries) ;
for ( jcol = ipivot = kk = 0 ; jcol <= icol ; ipivot++ ) {
m = pivotsizes[ipivot] ;
/*
fprintf(stdout, "\n jrow %d, m %d, kk %d", jrow, m, kk) ;
fflush(stdout) ;
*/
if ( jcol <= icol && icol < jcol + m ) {
stride = m - 1 ;
kk += icol - jcol ;
/*
fprintf(stdout, "\n 0. kk %d, stride %d", kk, stride) ;
fflush(stdout) ;
*/
for ( ii = jcol ; ii <= icol ; ii++ ) {
/*
fprintf(stdout, "\n 1. ii %d, kk %d", ii, kk) ;
fflush(stdout) ;
*/
colvec[2*ii] = entries[2*kk] ;
colvec[2*ii+1] = entries[2*kk+1] ;
kk += stride-- ;
}
for ( ; ii < jcol + m ; ii++ ) {
/*
fprintf(stdout, "\n 2. ii %d, kk %d", ii, kk) ;
fflush(stdout) ;
*/
colvec[2*ii] = entries[2*kk] ;
colvec[2*ii+1] = entries[2*kk+1] ;
kk++ ;
}
} else {
kk += (m*(m+1))/2 ;
}
jcol += m ;
}
} break ;
case SUBMTX_BLOCK_DIAGONAL_HERM : {
double *entries ;
int ii, ipivot, jcol, kk, m, nrow, nent, stride ;
int *pivotsizes ;
SubMtx_blockDiagonalInfo(mtx, &nrow, &nent, &pivotsizes, &entries) ;
for ( jcol = ipivot = kk = 0 ; jcol <= icol ; ipivot++ ) {
m = pivotsizes[ipivot] ;
/*
fprintf(stdout, "\n jrow %d, m %d, kk %d", jrow, m, kk) ;
fflush(stdout) ;
*/
if ( jcol <= icol && icol < jcol + m ) {
stride = m - 1 ;
kk += icol - jcol ;
/*
fprintf(stdout, "\n 0. kk %d, stride %d", kk, stride) ;
fflush(stdout) ;
*/
for ( ii = jcol ; ii <= icol ; ii++ ) {
/*
fprintf(stdout, "\n 1. ii %d, kk %d", ii, kk) ;
fflush(stdout) ;
*/
colvec[2*ii] = entries[2*kk] ;
colvec[2*ii+1] = entries[2*kk+1] ;
kk += stride-- ;
}
for ( ; ii < jcol + m ; ii++ ) {
/*
fprintf(stdout, "\n 2. ii %d, kk %d", ii, kk) ;
fflush(stdout) ;
*/
colvec[2*ii] = entries[2*kk] ;
colvec[2*ii+1] = -entries[2*kk+1] ;
kk++ ;
}
} else {
kk += (m*(m+1))/2 ;
}
jcol += m ;
}
} break ;
default :
break ;
}
return ; }
/*--------------------------------------------------------------------*/
/*
------------------------------------------------------
purpose -- return the magnitude of the largest element
created -- 98may01, cca
------------------------------------------------------
*/
double
SubMtx_maxabs (
SubMtx *mtx
) {
double maxabs ;
double *entries ;
int loc, nent ;
/*
---------------
check the input
---------------
*/
if ( mtx == NULL ) {
fprintf(stderr, "\n fatal error in SubMtx_maxabs(%p)"
"\n bad input\n", mtx) ;
exit(-1) ;
}
if ( ! (SUBMTX_IS_REAL(mtx) || SUBMTX_IS_COMPLEX(mtx)) ) {
fprintf(stderr, "\n fatal error in SubMtx_maxabs(%p)"
"\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
mtx, mtx->type) ;
exit(-1) ;
}
/*
--------------------------------------
switch over the different matrix types
--------------------------------------
*/
switch ( mtx->mode ) {
case SUBMTX_DENSE_ROWS :
case SUBMTX_DENSE_COLUMNS : {
int inc1, inc2, ncol, nrow ;
SubMtx_denseInfo(mtx, &nrow, &ncol, &inc1, &inc2, &entries) ;
nent = nrow*ncol ;
} break ;
case SUBMTX_SPARSE_COLUMNS : {
int ncol ;
int *indices, *sizes ;
SubMtx_sparseColumnsInfo(mtx, &ncol, &nent,
&sizes, &indices, &entries) ;
} break ;
case SUBMTX_SPARSE_ROWS : {
int nrow ;
int *indices, *sizes ;
SubMtx_sparseRowsInfo(mtx, &nrow, &nent, &sizes, &indices, &entries) ;
} break ;
case SUBMTX_SPARSE_TRIPLES : {
int *colids, *rowids ;
SubMtx_sparseTriplesInfo(mtx, &nent, &rowids, &colids, &entries) ;
} break ;
case SUBMTX_DENSE_SUBCOLUMNS : {
int ncol ;
int *firstlocs, *sizes ;
SubMtx_denseSubcolumnsInfo(mtx, &ncol, &nent,
&firstlocs, &sizes, &entries) ;
} break ;
case SUBMTX_DENSE_SUBROWS : {
int nrow ;
int *firstlocs, *sizes ;
SubMtx_denseSubrowsInfo(mtx, &nrow, &nent,
&firstlocs, &sizes, &entries) ;
} break ;
case SUBMTX_DIAGONAL : {
SubMtx_diagonalInfo(mtx, &nent, &entries) ;
} break ;
case SUBMTX_BLOCK_DIAGONAL_SYM : {
int nrow ;
int *pivotsizes ;
SubMtx_blockDiagonalInfo(mtx, &nrow, &nent, &pivotsizes, &entries) ;
} break ;
case SUBMTX_BLOCK_DIAGONAL_HERM : {
int nrow ;
int *pivotsizes ;
SubMtx_blockDiagonalInfo(mtx, &nrow, &nent, &pivotsizes, &entries) ;
} break ;
default :
break ;
}
if ( SUBMTX_IS_REAL(mtx) ) {
maxabs = DVmaxabs(nent, entries, &loc) ;
} else if ( SUBMTX_IS_COMPLEX(mtx) ) {
maxabs = ZVmaxabs(nent, entries) ;
}
return(maxabs) ; }
/*--------------------------------------------------------------------*/
/*
-----------------------------------------
purpose -- zero the entries in the matrix
created -- 98may04, cca
-----------------------------------------
*/
void
SubMtx_zero (
SubMtx *mtx
) {
/*
---------------
check the input
---------------
*/
if ( mtx == NULL ) {
fprintf(stderr, "\n fatal error in SubMtx_zero(%p)"
"\n bad input\n", mtx) ;
exit(-1) ;
}
if ( SUBMTX_IS_REAL(mtx) ) {
DVzero(mtx->nent, mtx->entries) ;
} else if ( SUBMTX_IS_COMPLEX(mtx) ) {
DVzero(2*mtx->nent, mtx->entries) ;
}
return ; }
/*--------------------------------------------------------------------*/
syntax highlighted by Code2HTML, v. 0.9.1