/* util.c */
#include "../InpMtx.h"
#define MYDEBUG 0
/*--------------------------------------------------------------------*/
/*
---------------------------------------
given the data is in raw triples,
sort and compress the data
created -- 98jan28, cca
modified -- 98sep04, cca
test to see if the sort is necessary
---------------------------------------
*/
void
InpMtx_sortAndCompress (
InpMtx *inpmtx
) {
int ient, nent, sortMustBeDone ;
int *ivec1, *ivec2 ;
/*
---------------
check the input
---------------
*/
if ( inpmtx == NULL ) {
fprintf(stderr, "\n fatal error in InpMtx_sortAndCompress(%p)"
"\n bad input\n", inpmtx) ;
exit(-1) ;
}
if ( INPMTX_IS_SORTED(inpmtx)
|| INPMTX_IS_BY_VECTORS(inpmtx)
|| (nent = inpmtx->nent) == 0 ) {
inpmtx->storageMode = INPMTX_SORTED ;
return ;
}
ivec1 = InpMtx_ivec1(inpmtx) ;
ivec2 = InpMtx_ivec2(inpmtx) ;
sortMustBeDone = 0 ;
for ( ient = 1 ; ient < nent ; ient++ ) {
if ( ivec1[ient-1] > ivec1[ient]
|| ( ivec1[ient-1] == ivec1[ient]
&& ivec2[ient-1] > ivec2[ient] ) ) {
sortMustBeDone = 1 ;
break ;
}
}
if ( sortMustBeDone == 1 ) {
if ( INPMTX_IS_INDICES_ONLY(inpmtx) ) {
inpmtx->nent = IV2sortUpAndCompress(nent, ivec1, ivec2) ;
} else if ( INPMTX_IS_REAL_ENTRIES(inpmtx) ) {
double *dvec = InpMtx_dvec(inpmtx) ;
inpmtx->nent = IV2DVsortUpAndCompress(nent, ivec1, ivec2, dvec) ;
} else if ( INPMTX_IS_COMPLEX_ENTRIES(inpmtx) ) {
double *dvec = InpMtx_dvec(inpmtx) ;
inpmtx->nent = IV2ZVsortUpAndCompress(nent, ivec1, ivec2, dvec) ;
}
}
inpmtx->storageMode = INPMTX_SORTED ;
return ; }
/*--------------------------------------------------------------------*/
/*
----------------------------------------------------
convert from sorted and compressed triples to vector
created -- 98jan28, cca
----------------------------------------------------
*/
void
InpMtx_convertToVectors (
InpMtx *inpmtx
) {
int *ivec1, *ivec2, *offsets, *sizes, *vecids ;
int first, id, ient, jj, nent, nvector, value ;
/*
---------------
check the input
---------------
*/
if ( inpmtx == NULL ) {
fprintf(stderr, "\n fatal error in InpMtx_convertToVectors(%p)"
"\n bad input\n", inpmtx) ;
exit(-1) ;
}
if ( INPMTX_IS_BY_VECTORS(inpmtx) || (nent = inpmtx->nent) == 0 ) {
inpmtx->storageMode = INPMTX_BY_VECTORS ;
return ;
}
if ( INPMTX_IS_RAW_DATA(inpmtx) ) {
InpMtx_sortAndCompress(inpmtx) ;
}
ivec1 = InpMtx_ivec1(inpmtx) ;
ivec2 = InpMtx_ivec2(inpmtx) ;
/*
-----------------------------------
find the number of distinct vectors
-----------------------------------
*/
value = -1 ;
nvector = 0 ;
for ( ient = 0 ; ient < nent ; ient++ ) {
if ( ivec1[ient] >= 0 && value != ivec1[ient] ) {
value = ivec1[ient] ;
nvector++ ;
#if MYDEBUG > 0
fprintf(stdout, "\n ient %d, value %d, nvector %d",
ient, value, nvector) ;
fflush(stdout) ;
#endif
}
}
#if MYDEBUG > 0
fprintf(stdout, "\n %d vectors", nvector) ;
fflush(stdout) ;
#endif
/*
-----------------------------------------------------------------
adjust the sizes of the sizes[] and offsets[] arrays if necessary
-----------------------------------------------------------------
*/
InpMtx_setNvector(inpmtx, nvector) ;
if ( nvector <= 0 ) {
/*
-----------------------------
matrix has no entries, return
-----------------------------
*/
inpmtx->storageMode = INPMTX_RAW_DATA ;
InpMtx_setNent(inpmtx, 0) ;
return ;
}
vecids = InpMtx_vecids(inpmtx) ;
sizes = InpMtx_sizes(inpmtx) ;
offsets = InpMtx_offsets(inpmtx) ;
/*
------------------------------------------------------------
set the vector sizes and offsets
note: we skip all entries whose first coordinate is negative
------------------------------------------------------------
*/
for ( first = 0 ; first < nent ; first++ ) {
if ( ivec1[first] >= 0 ) {
break ;
}
}
id = 0 ;
if ( first < nent ) {
value = ivec1[first] ;
for ( jj = first + 1 ; jj < nent ; jj++ ) {
if ( ivec1[jj] != value ) {
vecids[id] = value ;
sizes[id] = jj - first ;
offsets[id] = first ;
#if MYDEBUG > 0
fprintf(stdout,
"\n vecids[%d] = %d, sizes[%d] = %d, offsets[%d] = %d",
id, vecids[id], id, sizes[id], id, offsets[id]) ;
fflush(stdout) ;
#endif
first = jj ;
value = ivec1[jj] ;
id++ ;
}
}
vecids[id] = value ;
sizes[id] = jj - first ;
offsets[id] = first ;
#if MYDEBUG > 0
fprintf(stdout,
"\n vecids[%d] = %d, sizes[%d] = %d, offsets[%d] = %d",
id, vecids[id], id, sizes[id], id, offsets[id]) ;
fflush(stdout) ;
#endif
}
inpmtx->storageMode = INPMTX_BY_VECTORS ;
#if MYDEBUG > 0
fprintf(stdout, "\n vecidsIV") ;
IV_writeForHumanEye(&inpmtx->vecidsIV, stdout) ;
fprintf(stdout, "\n sizesIV") ;
IV_writeForHumanEye(&inpmtx->sizesIV, stdout) ;
fprintf(stdout, "\n offsetsIV") ;
IV_writeForHumanEye(&inpmtx->offsetsIV, stdout) ;
fflush(stdout) ;
#endif
return ; }
/*--------------------------------------------------------------------*/
/*
-------------------------
drop off-diagonal entries
created -- 98jan28, cca
-------------------------
*/
void
InpMtx_dropOffdiagonalEntries (
InpMtx *inpmtx
) {
double *dvec ;
int count, ii, nent ;
int *ivec1, *ivec2 ;
/*
---------------
check the input
---------------
*/
if ( inpmtx == NULL ) {
fprintf(stderr,
"\n fatal error in InpMtx_dropOffdiagonalEntries(%p)"
"\n bad input\n", inpmtx) ;
exit(-1) ;
}
if ( !( INPMTX_IS_BY_ROWS(inpmtx)
|| INPMTX_IS_BY_COLUMNS(inpmtx)
|| INPMTX_IS_BY_CHEVRONS(inpmtx) ) ) {
fprintf(stderr,
"\n fatal error in InpMtx_dropOffdiagonalEntries(%p)"
"\n bad coordType = %d\n", inpmtx, inpmtx->coordType) ;
exit(-1) ;
}
nent = inpmtx->nent ;
ivec1 = InpMtx_ivec1(inpmtx) ;
ivec2 = InpMtx_ivec2(inpmtx) ;
count = 0 ;
if ( INPMTX_IS_REAL_ENTRIES(inpmtx)
|| INPMTX_IS_COMPLEX_ENTRIES(inpmtx) ) {
dvec = InpMtx_dvec(inpmtx) ;
}
if ( INPMTX_IS_BY_ROWS(inpmtx) ) {
for ( ii = 0 ; ii < nent ; ii++ ) {
if ( ivec1[ii] == ivec2[ii] ) {
ivec1[count] = ivec1[ii] ;
ivec2[count] = ivec2[ii] ;
if ( INPMTX_IS_REAL_ENTRIES(inpmtx) ) {
dvec[count] = dvec[ii] ;
} else if ( INPMTX_IS_COMPLEX_ENTRIES(inpmtx) ) {
dvec[2*count] = dvec[2*ii] ;
dvec[2*count+1] = dvec[2*ii+1] ;
}
count++ ;
}
}
} else if ( INPMTX_IS_BY_COLUMNS(inpmtx) ) {
for ( ii = 0 ; ii < nent ; ii++ ) {
if ( ivec1[ii] == ivec2[ii] ) {
ivec1[count] = ivec1[ii] ;
ivec2[count] = ivec2[ii] ;
if ( INPMTX_IS_REAL_ENTRIES(inpmtx) ) {
dvec[count] = dvec[ii] ;
} else if ( INPMTX_IS_COMPLEX_ENTRIES(inpmtx) ) {
dvec[2*count] = dvec[2*ii] ;
dvec[2*count+1] = dvec[2*ii+1] ;
}
count++ ;
}
}
} else if ( INPMTX_IS_BY_CHEVRONS(inpmtx) ) {
for ( ii = 0 ; ii < nent ; ii++ ) {
if ( ivec2[ii] == 0 ) {
ivec1[count] = ivec1[ii] ;
ivec2[count] = ivec2[ii] ;
if ( INPMTX_IS_REAL_ENTRIES(inpmtx) ) {
dvec[count] = dvec[ii] ;
} else if ( INPMTX_IS_COMPLEX_ENTRIES(inpmtx) ) {
dvec[2*count] = dvec[2*ii] ;
dvec[2*count+1] = dvec[2*ii+1] ;
}
count++ ;
}
}
}
inpmtx->nent = count ;
IV_setSize(&inpmtx->ivec1IV, count) ;
IV_setSize(&inpmtx->ivec2IV, count) ;
if ( INPMTX_IS_REAL_ENTRIES(inpmtx)
|| INPMTX_IS_COMPLEX_ENTRIES(inpmtx) ) {
DV_setSize(&inpmtx->dvecDV, count) ;
}
return ; }
/*--------------------------------------------------------------------*/
/*
----------------------------------
drop entries in the lower triangle
created -- 98jan28, cca
----------------------------------
*/
void
InpMtx_dropLowerTriangle (
InpMtx *inpmtx
) {
double *dvec ;
int count, ii, nent ;
int *ivec1, *ivec2 ;
/*
---------------
check the input
---------------
*/
if ( inpmtx == NULL ) {
fprintf(stderr, "\n fatal error in InpMtx_dropLowerTriangle(%p)"
"\n bad input\n", inpmtx) ;
exit(-1) ;
}
if ( !( INPMTX_IS_BY_ROWS(inpmtx)
|| INPMTX_IS_BY_COLUMNS(inpmtx)
|| INPMTX_IS_BY_CHEVRONS(inpmtx) ) ) {
fprintf(stderr, "\n fatal error in InpMtx_dropLowerTriangle(%p)"
"\n bad coordType \n", inpmtx) ;
exit(-1) ;
}
nent = inpmtx->nent ;
ivec1 = InpMtx_ivec1(inpmtx) ;
ivec2 = InpMtx_ivec2(inpmtx) ;
count = 0 ;
if ( INPMTX_IS_REAL_ENTRIES(inpmtx)
|| INPMTX_IS_COMPLEX_ENTRIES(inpmtx) ) {
dvec = InpMtx_dvec(inpmtx) ;
}
if ( INPMTX_IS_BY_ROWS(inpmtx) ) {
for ( ii = 0 ; ii < nent ; ii++ ) {
/*
fprintf(stdout, "\n ii = %d, (%d,%d)", ii, ivec1[ii], ivec2[ii]) ;
*/
if ( ivec1[ii] <= ivec2[ii] ) {
ivec1[count] = ivec1[ii] ;
ivec2[count] = ivec2[ii] ;
/*
fprintf(stdout, "\n ivec1[%d] = %d, ivec2[%d] = %d",
count, ivec1[count], count, ivec2[count]) ;
*/
if ( INPMTX_IS_REAL_ENTRIES(inpmtx) ) {
dvec[count] = dvec[ii] ;
} else if ( INPMTX_IS_COMPLEX_ENTRIES(inpmtx) ) {
dvec[2*count] = dvec[2*ii] ;
dvec[2*count+1] = dvec[2*ii+1] ;
}
count++ ;
}
}
} else if ( INPMTX_IS_BY_COLUMNS(inpmtx) ) {
for ( ii = 0 ; ii < nent ; ii++ ) {
if ( ivec1[ii] >= ivec2[ii] ) {
ivec1[count] = ivec1[ii] ;
ivec2[count] = ivec2[ii] ;
if ( INPMTX_IS_REAL_ENTRIES(inpmtx) ) {
dvec[count] = dvec[ii] ;
} else if ( INPMTX_IS_COMPLEX_ENTRIES(inpmtx) ) {
dvec[2*count] = dvec[2*ii] ;
dvec[2*count+1] = dvec[2*ii+1] ;
}
count++ ;
}
}
} else if ( INPMTX_IS_BY_CHEVRONS(inpmtx) ) {
for ( ii = 0 ; ii < nent ; ii++ ) {
if ( ivec2[ii] >= 0 ) {
ivec1[count] = ivec1[ii] ;
ivec2[count] = ivec2[ii] ;
if ( INPMTX_IS_REAL_ENTRIES(inpmtx) ) {
dvec[count] = dvec[ii] ;
} else if ( INPMTX_IS_COMPLEX_ENTRIES(inpmtx) ) {
dvec[2*count] = dvec[2*ii] ;
dvec[2*count+1] = dvec[2*ii+1] ;
}
count++ ;
}
}
}
inpmtx->nent = count ;
IV_setSize(&inpmtx->ivec1IV, count) ;
IV_setSize(&inpmtx->ivec2IV, count) ;
if ( INPMTX_IS_REAL_ENTRIES(inpmtx) ) {
DV_setSize(&inpmtx->dvecDV, count) ;
} else if ( INPMTX_IS_COMPLEX_ENTRIES(inpmtx) ) {
DV_setSize(&inpmtx->dvecDV, 2*count) ;
}
return ; }
/*--------------------------------------------------------------------*/
/*
----------------------------------
drop entries in the upper triangle
created -- 98jan28, cca
----------------------------------
*/
void
InpMtx_dropUpperTriangle (
InpMtx *inpmtx
) {
double *dvec ;
int count, ii, nent ;
int *ivec1, *ivec2 ;
/*
---------------
check the input
---------------
*/
if ( inpmtx == NULL ) {
fprintf(stderr, "\n fatal error in InpMtx_dropUpperTriangle(%p)"
"\n bad input\n", inpmtx) ;
exit(-1) ;
}
if ( !( INPMTX_IS_BY_ROWS(inpmtx)
|| INPMTX_IS_BY_COLUMNS(inpmtx)
|| INPMTX_IS_BY_CHEVRONS(inpmtx) ) ) {
fprintf(stderr, "\n fatal error in InpMtx_dropUpperTriangle(%p)"
"\n bad coordType \n", inpmtx) ;
exit(-1) ;
}
nent = inpmtx->nent ;
ivec1 = InpMtx_ivec1(inpmtx) ;
ivec2 = InpMtx_ivec2(inpmtx) ;
count = 0 ;
if ( INPMTX_IS_REAL_ENTRIES(inpmtx)
|| INPMTX_IS_COMPLEX_ENTRIES(inpmtx) ) {
dvec = InpMtx_dvec(inpmtx) ;
}
if ( INPMTX_IS_BY_ROWS(inpmtx) ) {
for ( ii = 0 ; ii < nent ; ii++ ) {
if ( ivec1[ii] >= ivec2[ii] ) {
ivec1[count] = ivec1[ii] ;
ivec2[count] = ivec2[ii] ;
if ( INPMTX_IS_REAL_ENTRIES(inpmtx) ) {
dvec[count] = dvec[ii] ;
} else if ( INPMTX_IS_COMPLEX_ENTRIES(inpmtx) ) {
dvec[2*count] = dvec[2*ii] ;
dvec[2*count+1] = dvec[2*ii+1] ;
}
count++ ;
}
}
} else if ( INPMTX_IS_BY_COLUMNS(inpmtx) ) {
for ( ii = 0 ; ii < nent ; ii++ ) {
if ( ivec1[ii] <= ivec2[ii] ) {
ivec1[count] = ivec1[ii] ;
ivec2[count] = ivec2[ii] ;
if ( INPMTX_IS_REAL_ENTRIES(inpmtx) ) {
dvec[count] = dvec[ii] ;
} else if ( INPMTX_IS_COMPLEX_ENTRIES(inpmtx) ) {
dvec[2*count] = dvec[2*ii] ;
dvec[2*count+1] = dvec[2*ii+1] ;
}
count++ ;
}
}
} else if ( INPMTX_IS_BY_CHEVRONS(inpmtx) ) {
for ( ii = 0 ; ii < nent ; ii++ ) {
if ( ivec2[ii] <= 0 ) {
ivec1[count] = ivec1[ii] ;
ivec2[count] = ivec2[ii] ;
if ( INPMTX_IS_REAL_ENTRIES(inpmtx) ) {
dvec[count] = dvec[ii] ;
} else if ( INPMTX_IS_COMPLEX_ENTRIES(inpmtx) ) {
dvec[2*count] = dvec[2*ii] ;
dvec[2*count+1] = dvec[2*ii+1] ;
}
count++ ;
}
}
}
inpmtx->nent = count ;
IV_setSize(&inpmtx->ivec1IV, count) ;
IV_setSize(&inpmtx->ivec2IV, count) ;
if ( INPMTX_IS_REAL_ENTRIES(inpmtx)
|| INPMTX_IS_COMPLEX_ENTRIES(inpmtx) ) {
DV_setSize(&inpmtx->dvecDV, count) ;
}
return ; }
/*--------------------------------------------------------------------*/
/*
-----------------------------------
map entries into the lower triangle
created -- 98jan28, cca
-----------------------------------
*/
void
InpMtx_mapToLowerTriangle (
InpMtx *inpmtx
) {
int col, ii, nent, row ;
int *ivec1, *ivec2 ;
/*
---------------
check the input
---------------
*/
if ( inpmtx == NULL ) {
fprintf(stderr, "\n fatal error in InpMtx_mapToLowerTriangle(%p)"
"\n bad input\n", inpmtx) ;
exit(-1) ;
}
if ( !( INPMTX_IS_BY_ROWS(inpmtx)
|| INPMTX_IS_BY_COLUMNS(inpmtx)
|| INPMTX_IS_BY_CHEVRONS(inpmtx) ) ) {
fprintf(stderr, "\n fatal error in InpMtx_mapToLowerTriangle(%p)"
"\n bad coordType\n", inpmtx) ;
exit(-1) ;
}
nent = inpmtx->nent ;
ivec1 = InpMtx_ivec1(inpmtx) ;
ivec2 = InpMtx_ivec2(inpmtx) ;
if ( INPMTX_IS_BY_ROWS(inpmtx) ) {
for ( ii = 0 ; ii < nent ; ii++ ) {
if ( (row = ivec1[ii]) < (col = ivec2[ii]) ) {
ivec1[ii] = col ;
ivec2[ii] = row ;
}
}
} else if ( INPMTX_IS_BY_COLUMNS(inpmtx) ) {
for ( ii = 0 ; ii < nent ; ii++ ) {
if ( (row = ivec2[ii]) < (col = ivec1[ii]) ) {
ivec1[ii] = row ;
ivec2[ii] = col ;
}
}
} else if ( INPMTX_IS_BY_CHEVRONS(inpmtx) ) {
for ( ii = 0 ; ii < nent ; ii++ ) {
if ( ivec2[ii] > 0 ) {
ivec2[ii] = -ivec2[ii] ;
}
}
}
return ; }
/*--------------------------------------------------------------------*/
/*
-----------------------------------
map entries into the upper triangle
created -- 98jan28, cca
-----------------------------------
*/
void
InpMtx_mapToUpperTriangle (
InpMtx *inpmtx
) {
int col, ii, nent, row ;
int *ivec1, *ivec2 ;
/*
---------------
check the input
---------------
*/
if ( inpmtx == NULL ) {
fprintf(stderr, "\n fatal error in InpMtx_mapToUpperTriangle(%p)"
"\n bad input\n", inpmtx) ;
exit(-1) ;
}
if ( !( INPMTX_IS_BY_ROWS(inpmtx)
|| INPMTX_IS_BY_COLUMNS(inpmtx)
|| INPMTX_IS_BY_CHEVRONS(inpmtx) ) ) {
fprintf(stderr, "\n fatal error in InpMtx_mapToUpperTriangle(%p)"
"\n bad coordType = %d, must be 1, 2, or 3\n",
inpmtx, inpmtx->coordType) ;
exit(-1) ;
}
nent = inpmtx->nent ;
ivec1 = InpMtx_ivec1(inpmtx) ;
ivec2 = InpMtx_ivec2(inpmtx) ;
if ( INPMTX_IS_BY_ROWS(inpmtx) ) {
for ( ii = 0 ; ii < nent ; ii++ ) {
if ( (row = ivec1[ii]) > (col = ivec2[ii]) ) {
ivec1[ii] = col ;
ivec2[ii] = row ;
}
}
} else if ( INPMTX_IS_BY_COLUMNS(inpmtx) ) {
for ( ii = 0 ; ii < nent ; ii++ ) {
if ( (row = ivec2[ii]) > (col = ivec1[ii]) ) {
ivec1[ii] = row ;
ivec2[ii] = col ;
}
}
} else if ( INPMTX_IS_BY_CHEVRONS(inpmtx) ) {
for ( ii = 0 ; ii < nent ; ii++ ) {
if ( ivec2[ii] < 0 ) {
ivec2[ii] = -ivec2[ii] ;
}
}
}
inpmtx->storageMode = INPMTX_RAW_DATA ;
return ; }
/*--------------------------------------------------------------------*/
/*
-----------------------------------
map entries into the upper triangle
for a hermitian matrix
created -- 98jan28, cca
-----------------------------------
*/
void
InpMtx_mapToUpperTriangleH (
InpMtx *inpmtx
) {
double *dvec ;
int col, ii, nent, row ;
int *ivec1, *ivec2 ;
/*
---------------
check the input
---------------
*/
if ( inpmtx == NULL ) {
fprintf(stderr, "\n fatal error in InpMtx_mapToUpperTriangle(%p)"
"\n bad input\n", inpmtx) ;
exit(-1) ;
}
if ( !( INPMTX_IS_BY_ROWS(inpmtx)
|| INPMTX_IS_BY_COLUMNS(inpmtx)
|| INPMTX_IS_BY_CHEVRONS(inpmtx) ) ) {
fprintf(stderr, "\n fatal error in InpMtx_mapToUpperTriangle(%p)"
"\n bad coordType = %d, must be 1, 2, or 3\n",
inpmtx, inpmtx->coordType) ;
exit(-1) ;
}
if ( ! INPMTX_IS_COMPLEX_ENTRIES(inpmtx) ) {
fprintf(stderr, "\n fatal error in InpMtx_mapToUpperTriangleH(%p)"
"\n input mode is not complex\n", inpmtx) ;
exit(-1) ;
}
nent = inpmtx->nent ;
ivec1 = InpMtx_ivec1(inpmtx) ;
ivec2 = InpMtx_ivec2(inpmtx) ;
dvec = InpMtx_dvec(inpmtx) ;
if ( INPMTX_IS_BY_ROWS(inpmtx) ) {
for ( ii = 0 ; ii < nent ; ii++ ) {
if ( (row = ivec1[ii]) > (col = ivec2[ii]) ) {
ivec1[ii] = col ;
ivec2[ii] = row ;
dvec[2*ii+1] = -dvec[2*ii+1] ;
}
}
} else if ( INPMTX_IS_BY_COLUMNS(inpmtx) ) {
for ( ii = 0 ; ii < nent ; ii++ ) {
if ( (row = ivec2[ii]) > (col = ivec1[ii]) ) {
ivec1[ii] = row ;
ivec2[ii] = col ;
dvec[2*ii+1] = -dvec[2*ii+1] ;
}
}
} else if ( INPMTX_IS_BY_CHEVRONS(inpmtx) ) {
for ( ii = 0 ; ii < nent ; ii++ ) {
if ( ivec2[ii] < 0 ) {
ivec2[ii] = -ivec2[ii] ;
dvec[2*ii+1] = -dvec[2*ii+1] ;
}
}
}
inpmtx->storageMode = INPMTX_RAW_DATA ;
return ; }
/*--------------------------------------------------------------------*/
/*
-----------------------------------------------------------
purpose -- compute the checksums of the indices and entries
sums[0] = sum_{ii=0}^{nent} abs(ivec1[ii])
sums[1] = sum_{ii=0}^{nent} abs(ivec2[ii])
if real or complex entries then
sums[2] = sum_{ii=0}^{nent} magnitudes of entries
endif
created -- 98may16, cca
-----------------------------------------------------------
*/
void
InpMtx_checksums (
InpMtx *inpmtx,
double sums[]
) {
int ient, nent ;
int *ivec1, *ivec2 ;
/*
---------------
check the input
---------------
*/
if ( inpmtx == NULL ) {
fprintf(stderr, "\n fatal error in InpMtx_checksums(%p,%p)"
"\n bad input\n", inpmtx, sums) ;
exit(-1) ;
}
switch ( inpmtx->inputMode ) {
case INPMTX_INDICES_ONLY :
case SPOOLES_REAL :
case SPOOLES_COMPLEX :
break ;
default :
fprintf(stderr, "\n fatal error in InpMtx_checksums(%p,%p)"
"\n bad inputMode\n", inpmtx, sums) ;
exit(-1) ;
}
sums[0] = sums[1] = sums[2] = 0.0 ;
if ( (nent = InpMtx_nent(inpmtx)) <= 0 ) {
return ;
}
ivec1 = InpMtx_ivec1(inpmtx) ;
ivec2 = InpMtx_ivec2(inpmtx) ;
for ( ient = 0 ; ient < nent ; ient++ ) {
sums[0] += abs(ivec1[ient]) ;
sums[1] += abs(ivec2[ient]) ;
}
switch ( inpmtx->inputMode ) {
case INPMTX_INDICES_ONLY :
break ;
case SPOOLES_REAL : {
double *dvec = InpMtx_dvec(inpmtx) ;
for ( ient = 0 ; ient < nent ; ient++ ) {
sums[2] += fabs(dvec[ient]) ;
}
} break ;
case SPOOLES_COMPLEX : {
double *dvec = InpMtx_dvec(inpmtx) ;
for ( ient = 0 ; ient < nent ; ient++ ) {
sums[2] += Zabs(dvec[2*ient], dvec[2*ient+1]) ;
}
} break ;
}
return ; }
/*--------------------------------------------------------------------*/
/*
----------------------------------------------------------------
purpose -- to create an InpMtx object filled with random entries
input --
mtx -- matrix object, if NULL, it is created
inputMode -- input mode for the object,
indices only, real or complex entries
coordType -- coordinate type for the object,
by rows, by columns or by chevrons
storageMode -- storage mode for the object,
raw data, sorted or by vectors
nrow -- # of rows
ncol -- # of columns
symflag -- symmetry flag for the matrix,
symmetric, hermitian or nonsymmetric
nonzerodiag -- if 1, entries are placed on the diagonal
nitem -- # of items to be placed into the matrix
seed -- random number seed
return value ---
1 -- normal return
-1 -- mtx is NULL
-2 -- bad input mode
-3 -- bad coordinate type
-4 -- bad storage mode
-5 -- nrow or ncol <= 0
-6 -- bad symmetry flag
-7 -- hermitian matrix but not complex
-8 -- symmetric or hermitian matrix but nrow != ncol
-9 -- nitem < 0
----------------------------------------------------------------
*/
int
InpMtx_randomMatrix (
InpMtx *mtx,
int inputMode,
int coordType,
int storageMode,
int nrow,
int ncol,
int symflag,
int nonzerodiag,
int nitem,
int seed
) {
double *dvec ;
Drand *drand ;
int col, ii, neqns, row ;
int *colids, *rowids ;
/*
---------------
check the input
---------------
*/
if ( mtx == NULL ) {
fprintf(stderr, "\n fatal error in InpMtx_randomMatrix"
"\n mtx is NULL\n") ;
return(-1) ;
}
switch ( inputMode ) {
case INPMTX_INDICES_ONLY :
case SPOOLES_REAL :
case SPOOLES_COMPLEX :
break ;
default :
fprintf(stderr, "\n fatal error in InpMtx_randomMatrix"
"\n bad input mode %d\n", inputMode) ;
return(-2) ;
break ;
}
switch ( coordType ) {
case INPMTX_BY_ROWS :
case INPMTX_BY_COLUMNS :
case INPMTX_BY_CHEVRONS :
break ;
default :
fprintf(stderr, "\n fatal error in InpMtx_randomMatrix"
"\n bad coordinate type %d\n", coordType) ;
return(-3) ;
break ;
}
switch ( storageMode ) {
case INPMTX_RAW_DATA :
case INPMTX_SORTED :
case INPMTX_BY_VECTORS :
break ;
default :
fprintf(stderr, "\n fatal error in InpMtx_randomMatrix"
"\n bad storage mode%d\n", storageMode) ;
return(-4) ;
break ;
}
if ( nrow <= 0 || ncol <= 0 ) {
fprintf(stderr, "\n fatal error in InpMtx_randomMatrix"
"\n nrow = %d, ncol = %d\n", nrow, ncol) ;
return(-5) ;
}
switch ( symflag ) {
case SPOOLES_SYMMETRIC :
case SPOOLES_HERMITIAN :
case SPOOLES_NONSYMMETRIC :
break ;
default :
fprintf(stderr, "\n fatal error in InpMtx_randomMatrix"
"\n bad symmetry flag%d\n", symflag) ;
return(-6) ;
break ;
}
if ( symflag == SPOOLES_HERMITIAN && inputMode != SPOOLES_COMPLEX ) {
fprintf(stderr, "\n fatal error in InpMtx_randomMatrix"
"\n symmetryflag is Hermitian, requires complex type\n") ;
return(-7) ;
}
if ( (symflag == SPOOLES_SYMMETRIC || symflag == SPOOLES_HERMITIAN)
&& nrow != ncol ) {
fprintf(stderr, "\n fatal error in InpMtx_randomMatrix"
"\n symmetric or hermitian matrix, nrow %d, ncol%d\n",
nrow, ncol) ;
return(-8) ;
}
if ( nitem < 0 ) {
fprintf(stderr, "\n fatal error in InpMtx_randomMatrix"
"\n nitem = %d\n", nitem) ;
return(-9) ;
}
/*--------------------------------------------------------------------*/
neqns = (nrow <= ncol) ? nrow : ncol ;
if ( nonzerodiag == 1 ) {
nitem += neqns ;
}
/*
---------------------
initialize the object
---------------------
*/
InpMtx_init(mtx, INPMTX_BY_ROWS, inputMode, nitem, 0) ;
/*
----------------
fill the triples
----------------
*/
drand = Drand_new() ;
Drand_setSeed(drand, seed) ;
rowids = IVinit(nitem, -1) ;
colids = IVinit(nitem, -1) ;
if ( nonzerodiag == 1 ) {
IVramp(neqns, rowids, 0, 1) ;
Drand_setUniform(drand, 0, nrow) ;
Drand_fillIvector(drand, nitem - neqns, rowids + neqns) ;
Drand_setUniform(drand, 0, ncol) ;
IVramp(neqns, colids, 0, 1) ;
Drand_fillIvector(drand, nitem - neqns, colids + neqns) ;
} else {
Drand_setUniform(drand, 0, nrow) ;
Drand_fillIvector(drand, nitem, rowids) ;
Drand_setUniform(drand, 0, ncol) ;
Drand_fillIvector(drand, nitem, colids) ;
}
if ( symflag == SPOOLES_SYMMETRIC || symflag == SPOOLES_HERMITIAN ) {
for ( ii = 0 ; ii < nitem ; ii++ ) {
if ( (row = rowids[ii]) > (col = colids[ii]) ) {
rowids[ii] = col ;
colids[ii] = row ;
}
}
}
if ( inputMode == SPOOLES_REAL ) {
dvec = DVinit(nitem, 0.0) ;
Drand_setUniform(drand, 0.0, 1.0) ;
Drand_fillDvector(drand, nitem, dvec) ;
} else if ( inputMode == SPOOLES_COMPLEX ) {
dvec = DVinit(2*nitem, 0.0) ;
Drand_setUniform(drand, 0.0, 1.0) ;
Drand_fillDvector(drand, 2*nitem, dvec) ;
if ( symflag == SPOOLES_HERMITIAN ) {
for ( ii = 0 ; ii < nitem ; ii++ ) {
if ( rowids[ii] == colids[ii] ) {
dvec[2*ii+1] = 0.0 ;
}
}
}
} else {
dvec = NULL ;
}
/*
----------------
load the triples
----------------
*/
switch ( inputMode ) {
case INPMTX_INDICES_ONLY :
InpMtx_inputTriples(mtx, nitem, rowids, colids) ;
break ;
case SPOOLES_REAL :
InpMtx_inputRealTriples(mtx, nitem, rowids, colids, dvec) ;
break ;
case SPOOLES_COMPLEX :
InpMtx_inputComplexTriples(mtx, nitem, rowids, colids, dvec) ;
break ;
}
/*
----------------------------------------
set the coordinate type and storage mode
----------------------------------------
*/
InpMtx_changeCoordType(mtx, coordType) ;
InpMtx_changeStorageMode(mtx, storageMode) ;
/*
------------------------
free the working storage
------------------------
*/
Drand_free(drand) ;
IVfree(rowids) ;
IVfree(colids) ;
if ( dvec != NULL ) {
DVfree(dvec) ;
}
return(1) ; }
/*--------------------------------------------------------------------*/
/*
----------------------------------------------------
determine the range of the matrix,
i.e., the minimum and maximum rows and columns
if pmincol != NULL then *pmincol = minimum column id
if pmaxcol != NULL then *pmaxcol = maximum column id
if pminrow != NULL then *pminrow = minimum row id
if pmaxrow != NULL then *pmaxrow = maximum row id
return value ---
1 -- normal return
-1 -- mtx is NULL
-2 -- no entries in the matrix
-3 -- invalid coordinate type
created -- 98oct15, cca
----------------------------------------------------
*/
int
InpMtx_range (
InpMtx *mtx,
int *pmincol,
int *pmaxcol,
int *pminrow,
int *pmaxrow
) {
int maxcol, maxrow, mincol, minrow, nent ;
/*
---------------
check the input
---------------
*/
if ( mtx == NULL ) {
fprintf(stderr, "\n fatal error in InpMtx_range()"
"\n mtx is NULL\n") ;
return(-1) ;
}
if ( (nent = mtx->nent) <= 0 ) {
fprintf(stderr, "\n fatal error in InpMtx_range()"
"\n %d entries\n", nent) ;
return(-2) ;
}
if ( INPMTX_IS_BY_ROWS(mtx) ) {
int *rowind = InpMtx_ivec1(mtx) ;
int *colind = InpMtx_ivec2(mtx) ;
int col, ii, row ;
minrow = maxrow = rowind[0] ;
mincol = maxcol = colind[0] ;
for ( ii = 1 ; ii < nent ; ii++ ) {
row = rowind[ii] ;
col = colind[ii] ;
if ( minrow > row ) {
minrow = row ;
} else if ( maxrow < row ) {
maxrow = row ;
}
if ( mincol > col ) {
mincol = col ;
} else if ( maxcol < col ) {
maxcol = col ;
}
}
} else if ( INPMTX_IS_BY_COLUMNS(mtx) ) {
int *rowind = InpMtx_ivec2(mtx) ;
int *colind = InpMtx_ivec1(mtx) ;
int col, ii, row ;
minrow = maxrow = rowind[0] ;
mincol = maxcol = colind[0] ;
for ( ii = 1 ; ii < nent ; ii++ ) {
row = rowind[ii] ;
col = colind[ii] ;
if ( minrow > row ) {
minrow = row ;
} else if ( maxrow < row ) {
maxrow = row ;
}
if ( mincol > col ) {
mincol = col ;
} else if ( maxcol < col ) {
maxcol = col ;
}
}
} else if ( INPMTX_IS_BY_CHEVRONS(mtx) ) {
int *chvind = InpMtx_ivec1(mtx) ;
int *offsets = InpMtx_ivec2(mtx) ;
int col, ii, off, row ;
if ( (off = offsets[0]) >= 0 ) {
row = chvind[0], col = row + off ;
} else {
col = chvind[0], row = col - off ;
}
minrow = maxrow = row ;
mincol = maxcol = col ;
for ( ii = 1 ; ii < nent ; ii++ ) {
if ( (off = offsets[ii]) >= 0 ) {
row = chvind[ii], col = row + off ;
} else {
col = chvind[ii], row = col - off ;
}
if ( minrow > row ) {
minrow = row ;
} else if ( maxrow < row ) {
maxrow = row ;
}
if ( mincol > col ) {
mincol = col ;
} else if ( maxcol < col ) {
maxcol = col ;
}
}
} else {
fprintf(stderr, "\n fatal error in InpMtx_range()"
"\n invalid coordType %d\n", mtx->coordType) ;
return(-3) ;
}
if ( pmincol != NULL ) { *pmincol = mincol ; }
if ( pmaxcol != NULL ) { *pmaxcol = maxcol ; }
if ( pminrow != NULL ) { *pminrow = minrow ; }
if ( pmaxrow != NULL ) { *pmaxrow = maxrow ; }
return(1) ; }
/*--------------------------------------------------------------------*/
syntax highlighted by Code2HTML, v. 0.9.1