/* util.c */
#include "../Chv.h"
#define MYDEBUG 0
/*--------------------------------------------------------------------*/
/*
-------------------------------------------------------
shift the indices, entries and adjust the nD dimension.
note: shift can be positive or negative
created -- 98apr30, cca
-------------------------------------------------------
*/
void
Chv_shift (
Chv *chv,
int shift
) {
int ii, stride ;
/*
---------------
check the input
---------------
*/
if ( chv == NULL ) {
fprintf(stderr, "\n fatal error in Chv_shift(%p,%d)"
"\n bad input\n", chv, shift) ;
exit(-1) ;
}
if ( shift == 0 ) {
return ;
}
if ( CHV_IS_REAL(chv) ) {
if ( CHV_IS_SYMMETRIC(chv) ) {
/*
--------------------
chevron is symmetric
--------------------
*/
chv->colind += shift ;
if ( shift < 0 ) {
stride = chv->nD + chv->nU + 1 ;
for ( ii = shift ; ii < 0 ; ii++ ) {
chv->entries -= stride ;
stride++ ;
}
} else {
stride = chv->nD + chv->nU ;
for ( ii = 0 ; ii < shift ; ii++ ) {
chv->entries += stride ;
stride-- ;
}
}
chv->nD -= shift ;
} else if ( CHV_IS_NONSYMMETRIC(chv) ) {
/*
-----------------------
chevron is nonsymmetric
-----------------------
*/
chv->rowind += shift ;
chv->colind += shift ;
if ( shift < 0 ) {
stride = 2*chv->nD + chv->nL + chv->nU + 1 ;
for ( ii = shift ; ii < 0 ; ii++ ) {
chv->entries -= stride ;
stride += 2 ;
}
} else {
stride = 2*chv->nD + chv->nL + chv->nU - 1 ;
for ( ii = 0 ; ii < shift ; ii++ ) {
chv->entries += stride ;
stride -= 2 ;
}
}
chv->nD -= shift ;
} else {
fprintf(stderr, "\n fatal error in Chv_shift(%p,%d)"
"\n type is SPOOLES_REAL, symflag = %d"
"\n must be SPOOLES_SYMMETRIC or SPOOLES_NONSYMMETRIC\n",
chv, shift, chv->symflag) ;
exit(-1) ;
}
} else if ( CHV_IS_COMPLEX(chv) ) {
if ( CHV_IS_SYMMETRIC(chv) || CHV_IS_HERMITIAN(chv) ) {
/*
--------------------
chevron is symmetric
--------------------
*/
chv->colind += shift ;
if ( shift < 0 ) {
stride = chv->nD + chv->nU + 1 ;
for ( ii = shift ; ii < 0 ; ii++ ) {
chv->entries -= 2*stride ;
stride++ ;
}
} else {
stride = chv->nD + chv->nU ;
for ( ii = 0 ; ii < shift ; ii++ ) {
chv->entries += 2*stride ;
stride-- ;
}
}
chv->nD -= shift ;
} else if ( CHV_IS_NONSYMMETRIC(chv) ) {
/*
-----------------------
chevron is nonsymmetric
-----------------------
*/
chv->rowind += shift ;
chv->colind += shift ;
if ( shift < 0 ) {
stride = 2*chv->nD + chv->nL + chv->nU + 1 ;
for ( ii = shift ; ii < 0 ; ii++ ) {
chv->entries -= 2*stride ;
stride += 2 ;
}
} else {
stride = 2*chv->nD + chv->nL + chv->nU - 1 ;
for ( ii = 0 ; ii < shift ; ii++ ) {
chv->entries += 2*stride ;
stride -= 2 ;
}
}
chv->nD -= shift ;
} else {
fprintf(stderr, "\n fatal error in Chv_shift(%p,%d)"
"\n type is SPOOLES_COMPLEX, symflag = %d"
"\n must be SPOOLES_SYMMETRIC, SPOOLES_HERMITIAN"
"\n or SPOOLES_NONSYMMETRIC\n",
chv, shift, chv->symflag) ;
exit(-1) ;
}
}
return ; }
/*--------------------------------------------------------------------*/
/*
----------------------------------------------------------
return the maximum magnitude of the entries in the chevron
created -- 98apr30, cca
----------------------------------------------------------
*/
double
Chv_maxabs (
Chv *chv
) {
double maxabs ;
/*
---------------
check the input
---------------
*/
if ( chv == NULL ) {
fprintf(stderr, "\n fatal error in Chv_maxabs(%p)"
"\n bad input\n", chv) ;
exit(-1) ;
}
maxabs = 0.0 ;
if ( CHV_IS_REAL(chv) ) {
int loc ;
maxabs = DVmaxabs(Chv_nent(chv), Chv_entries(chv), &loc) ;
} else if ( CHV_IS_COMPLEX(chv) ) {
maxabs = ZVmaxabs(Chv_nent(chv), Chv_entries(chv)) ;
} else {
fprintf(stderr, "\n fatal error in Chv_maxabs(%p)"
"\n type is %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
chv, chv->type) ;
exit(-1) ;
}
return(maxabs) ; }
/*--------------------------------------------------------------------*/
/*
-------------------------------------------------------
return the frobenius norm of the entries in the chevron
created -- 98apr30, cca
-------------------------------------------------------
*/
double
Chv_frobNorm (
Chv *chv
) {
double sum ;
/*
---------------
check the input
---------------
*/
if ( chv == NULL ) {
fprintf(stderr, "\n fatal error in Chv_frobNorm(%p)"
"\n bad input\n", chv) ;
exit(-1) ;
}
if ( CHV_IS_REAL(chv) ) {
double value, *entries ;
int ii, nent ;
nent = Chv_nent(chv) ;
entries = Chv_entries(chv) ;
for ( ii = 0, sum = 0.0 ; ii < nent ; ii++ ) {
value = entries[ii] ;
sum += value*value ;
}
} else if ( CHV_IS_COMPLEX(chv) ) {
double imag, real, *entries ;
int ii, nent ;
nent = Chv_nent(chv) ;
entries = Chv_entries(chv) ;
for ( ii = 0, sum = 0.0 ; ii < nent ; ii++ ) {
real = entries[2*ii] ;
imag = entries[2*ii+1] ;
sum += real*real + imag*imag ;
}
} else {
fprintf(stderr, "\n fatal error in Chv_frobNorm(%p)"
"\n type is %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
chv, chv->type) ;
exit(-1) ;
}
return(sqrt(sum)) ; }
/*--------------------------------------------------------------------*/
/*
-----------------------
subtract chvI from chvJ
created -- 98apr30, cca
-----------------------
*/
void
Chv_sub (
Chv *chvJ,
Chv *chvI
) {
double *entriesI, *entriesJ ;
int ii, nDI, nDJ, nent, nLI, nLJ, nUI, nUJ ;
/*
---------------
check the input
---------------
*/
if ( chvI == NULL || chvJ == NULL ) {
fprintf(stderr, "\n fatal error in Chv_sub(%p,%p)"
"\n bad input\n", chvI, chvJ) ;
exit(-1) ;
}
Chv_dimensions(chvJ, &nDJ, &nLJ, &nUJ) ;
Chv_dimensions(chvI, &nDI, &nLI, &nUI) ;
if ( nDJ != nDI || nLJ != nLI || nUJ != nUI ) {
fprintf(stderr, "\n fatal error in Chv_sub(%p,%p)"
"\n dimensions do not match\n", chvJ, chvI) ;
exit(-1) ;
}
entriesJ = Chv_entries(chvJ) ;
entriesI = Chv_entries(chvI) ;
if ( entriesJ == NULL || entriesI == NULL ) {
fprintf(stderr, "\n fatal error in Chv_sub(%p,%p)"
"\n entriesJ = %p, entriesI = %p\n",
chvJ, chvI, entriesJ, entriesI) ;
exit(-1) ;
}
if ( CHV_IS_REAL(chvJ) && CHV_IS_REAL(chvI) ) {
nent = Chv_nent(chvJ) ;
for ( ii = 0 ; ii < nent ; ii++ ) {
entriesJ[ii] -= entriesI[ii] ;
}
} else if ( CHV_IS_COMPLEX(chvJ) && CHV_IS_COMPLEX(chvI) ) {
nent = Chv_nent(chvJ) ;
for ( ii = 0 ; ii < nent ; ii++ ) {
entriesJ[2*ii] -= entriesI[2*ii] ;
entriesJ[2*ii+1] -= entriesI[2*ii+1] ;
}
} else {
fprintf(stderr, "\n fatal error in Chv_sub(%p,%p)"
"\n typeJ = %d, typeI = %d"
"\n both must be SPOOLES_REAL or SPOOLES_COMPLEX",
chvJ, chvI, chvJ->type, chvI->type) ;
exit(-1) ;
}
return ; }
/*--------------------------------------------------------------------*/
/*
-------------------------------
zero the entries in the chevron
created -- 98apr30, cca
-------------------------------
*/
void
Chv_zero (
Chv *chv
) {
/*
---------------
check the input
---------------
*/
if ( chv == NULL ) {
fprintf(stderr, "\n fatal error in Chv_zero(%p)"
"\n bad input\n", chv) ;
exit(-1) ;
}
if ( CHV_IS_REAL(chv) ) {
DVzero(Chv_nent(chv), Chv_entries(chv)) ;
} else if ( CHV_IS_COMPLEX(chv) ) {
ZVzero(Chv_nent(chv), Chv_entries(chv)) ;
} else {
fprintf(stderr, "\n fatal error in Chv_zero(%p)"
"\n type = %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
chv, chv->type) ;
exit(-1) ;
}
return ; }
/*--------------------------------------------------------------------*/
/*
-------------------------------
fill A2 object with (1,1) block
created -- 98apr30, cca
-------------------------------
*/
void
Chv_fill11block (
Chv *chv,
A2 *mtx
) {
double *entries ;
int nD, nL, nU ;
/*
---------------
check the input
---------------
*/
if ( chv == NULL || mtx == NULL ) {
fprintf(stderr, "\n fatal error in Chv_fill11block(%p,%p)"
"\n bad input\n", chv, mtx) ;
exit(-1) ;
}
if ( ! (CHV_IS_REAL(chv) || CHV_IS_COMPLEX(chv)) ) {
fprintf(stderr, "\n fatal error in Chv_fill11block(%p,%p)"
"\n type = %d, must be real or complex\n",
chv, mtx, chv->type) ;
exit(-1) ;
}
Chv_dimensions(chv, &nD, &nL, &nU) ;
entries = Chv_entries(chv) ;
if ( CHV_IS_REAL(chv) ) {
A2_init(mtx, SPOOLES_REAL, nD, nD, 1, nD, NULL) ;
A2_zero(mtx) ;
if ( CHV_IS_SYMMETRIC(chv) ) {
int ii, iioff, ioff, istride, jj ;
/*
----------------
chv is symmetric
----------------
*/
ioff = 0 ;
istride = nD + nU ;
for ( ii = 0 ; ii < nD ; ii++ ) {
A2_setRealEntry(mtx, ii, ii, entries[ioff]) ;
for ( jj = ii + 1, iioff = ioff + 1 ;
jj < nD ;
jj++, iioff++ ) {
A2_setRealEntry(mtx, ii, jj, entries[iioff]) ;
A2_setRealEntry(mtx, jj, ii, 0.0) ;
}
ioff += istride ;
istride-- ;
}
} else if ( CHV_IS_NONSYMMETRIC(chv) ) {
int ii, iioff, ioff, istride, jj ;
/*
---------------------
chv is nonsymmetric
---------------------
*/
ioff = nD + nL - 1 ;
istride = 2*nD + nU + nL - 2 ;
for ( ii = 0 ; ii < nD ; ii++ ) {
A2_setRealEntry(mtx, ii, ii, entries[ioff]) ;
for ( jj = ii + 1, iioff = ioff + 1 ;
jj < nD ;
jj++, iioff++ ) {
A2_setRealEntry(mtx, ii, jj, entries[iioff]) ;
}
for ( jj = ii + 1, iioff = ioff - 1 ;
jj < nD ;
jj++, iioff-- ) {
A2_setRealEntry(mtx, jj, ii, entries[iioff]) ;
}
ioff += istride ;
istride -= 2 ;
}
}
} else if ( CHV_IS_COMPLEX(chv) ) {
A2_init(mtx, SPOOLES_COMPLEX, nD, nD, 1, nD, NULL) ;
A2_zero(mtx) ;
if ( CHV_IS_SYMMETRIC(chv) || CHV_IS_HERMITIAN(chv) ) {
int ii, iioff, ioff, istride, jj ;
/*
----------------
chv is symmetric
----------------
*/
ioff = 0 ;
istride = nD + nU ;
for ( ii = 0 ; ii < nD ; ii++ ) {
A2_setComplexEntry(mtx, ii, ii,
entries[2*ioff], entries[2*ioff+1]) ;
for ( jj = ii + 1, iioff = ioff + 1 ;
jj < nD ;
jj++, iioff++ ) {
A2_setComplexEntry(mtx, ii, jj,
entries[2*iioff], entries[2*iioff+1]) ;
A2_setComplexEntry(mtx, jj, ii, 0.0, 0.0) ;
}
ioff += istride ;
istride-- ;
}
} else if ( CHV_IS_NONSYMMETRIC(chv) ) {
int ii, iioff, ioff, istride, jj ;
/*
---------------------
chv is nonsymmetric
---------------------
*/
ioff = nD + nL - 1 ;
istride = 2*nD + nU + nL - 2 ;
for ( ii = 0 ; ii < nD ; ii++ ) {
#if MYDEBUG > 0
fprintf(stdout,
"\n ii %d, ioff %d, setting (%d,%d) = (%20.12e,%20.12e)",
ii, ioff, ii, ii, entries[2*ioff], entries[2*ioff+1]) ;
#endif
A2_setComplexEntry(mtx, ii, ii,
entries[2*ioff], entries[2*ioff+1]) ;
for ( jj = ii + 1, iioff = ioff + 1 ;
jj < nD ;
jj++, iioff++ ) {
#if MYDEBUG > 0
fprintf(stdout,
"\n jj %d, iioff %d, setting (%d,%d) = (%20.12e,%20.12e)",
jj, iioff, ii, ii, entries[2*iioff], entries[2*iioff+1]) ;
#endif
A2_setComplexEntry(mtx, ii, jj,
entries[2*iioff], entries[2*iioff+1]) ;
}
for ( jj = ii + 1, iioff = ioff - 1 ;
jj < nD ;
jj++, iioff-- ) {
#if MYDEBUG > 0
fprintf(stdout,
"\n jj %d, iioff %d, setting (%d,%d) = (%20.12e,%20.12e)",
jj, iioff, ii, ii, entries[2*iioff], entries[2*iioff+1]) ;
#endif
A2_setComplexEntry(mtx, jj, ii,
entries[2*iioff], entries[2*iioff+1]) ;
}
ioff += istride ;
istride -= 2 ;
}
}
}
return ; }
/*--------------------------------------------------------------------*/
/*
-------------------------------
fill A2 object with (1,2) block
created -- 98apr30, cca
-------------------------------
*/
void
Chv_fill12block (
Chv *chv,
A2 *mtx
) {
double *entries ;
int nD, nL, nU ;
/*
---------------
check the input
---------------
*/
if ( chv == NULL || mtx == NULL ) {
fprintf(stderr, "\n fatal error in Chv_fill12block(%p,%p)"
"\n bad input\n", chv, mtx) ;
exit(-1) ;
}
if ( ! (CHV_IS_REAL(chv) || CHV_IS_COMPLEX(chv)) ) {
fprintf(stderr, "\n fatal error in Chv_fill12block(%p,%p)"
"\n type = %d, must be real or complex\n",
chv, mtx, chv->type) ;
exit(-1) ;
}
Chv_dimensions(chv, &nD, &nL, &nU) ;
entries = Chv_entries(chv) ;
/*
---------------------------------
resize the A2 object as necessary
---------------------------------
*/
if ( CHV_IS_REAL(chv) ) {
A2_init(mtx, SPOOLES_REAL, nD, nU, 1, nD, NULL) ;
A2_zero(mtx) ;
if ( CHV_IS_SYMMETRIC(chv) ) {
int ii, iioff, ioff, istride, jj ;
/*
------------------
chv is symmetric
------------------
*/
ioff = 0 ;
istride = nD + nU ;
for ( ii = 0 ; ii < nD ; ii++ ) {
for ( jj = 0, iioff = ioff + nD - ii ;
jj < nU ;
jj++, iioff++ ) {
A2_setRealEntry(mtx, ii, jj, entries[iioff]) ;
}
ioff += istride ;
istride-- ;
}
} else if ( CHV_IS_NONSYMMETRIC(chv) ) {
int ii, iioff, ioff, istride, jj ;
/*
---------------------
chv is nonsymmetric
---------------------
*/
ioff = nD + nL - 1 ;
istride = 2*nD + nU + nL - 2 ;
for ( ii = 0 ; ii < nD ; ii++ ) {
for ( jj = 0, iioff = ioff + nD - ii ;
jj < nU ;
jj++, iioff++ ) {
A2_setRealEntry(mtx, ii, jj, entries[iioff]) ;
}
ioff += istride ;
istride -= 2 ;
}
}
} else if ( CHV_IS_COMPLEX(chv) ) {
A2_init(mtx, SPOOLES_COMPLEX, nD, nU, 1, nD, NULL) ;
A2_zero(mtx) ;
if ( CHV_IS_SYMMETRIC(chv) || CHV_IS_HERMITIAN(chv) ) {
int ii, iioff, ioff, istride, jj ;
/*
------------------
chv is symmetric
------------------
*/
ioff = 0 ;
istride = nD + nU ;
for ( ii = 0 ; ii < nD ; ii++ ) {
for ( jj = 0, iioff = ioff + nD - ii ;
jj < nU ;
jj++, iioff++ ) {
A2_setComplexEntry(mtx, ii, jj,
entries[2*iioff], entries[2*iioff+1]) ;
#if MYDEBUG > 0
fprintf(stdout,
"\n 21: ii %d, jj %d, iioff %d, setting (%d,%d) = (%20.12e,%20.12e)",
ii, jj, iioff, ii, ii, entries[2*iioff], entries[2*iioff+1]) ;
#endif
}
ioff += istride ;
istride-- ;
}
} else if ( CHV_IS_NONSYMMETRIC(chv) ) {
int ii, iioff, ioff, istride, jj ;
/*
---------------------
chv is nonsymmetric
---------------------
*/
ioff = nD + nL - 1 ;
istride = 2*nD + nU + nL - 2 ;
for ( ii = 0 ; ii < nD ; ii++ ) {
for ( jj = 0, iioff = ioff + nD - ii ;
jj < nU ;
jj++, iioff++ ) {
#if MYDEBUG > 0
fprintf(stdout,
"\n 21: ii %d, jj %d, iioff %d, setting (%d,%d) = (%20.12e,%20.12e)",
ii, jj, iioff, ii, ii, entries[2*iioff], entries[2*iioff+1]) ;
#endif
A2_setComplexEntry(mtx, ii, jj,
entries[2*iioff], entries[2*iioff+1]) ;
}
ioff += istride ;
istride -= 2 ;
}
}
}
return ; }
/*--------------------------------------------------------------------*/
/*
-------------------------------
fill A2 object with (2,1) block
created -- 98apr30, cca
-------------------------------
*/
void
Chv_fill21block (
Chv *chv,
A2 *mtx
) {
double *entries ;
int nD, nL, nU ;
/*
---------------
check the input
---------------
*/
if ( chv == NULL || mtx == NULL ) {
fprintf(stderr, "\n fatal error in Chv_fillReal21block(%p,%p)"
"\n bad input\n", chv, mtx) ;
exit(-1) ;
}
if ( ! (CHV_IS_REAL(chv) || CHV_IS_COMPLEX(chv)) ) {
fprintf(stderr, "\n fatal error in Chv_fill21block(%p,%p)"
"\n type = %d, must be real or complex\n",
chv, mtx, chv->type) ;
exit(-1) ;
}
Chv_dimensions(chv, &nD, &nL, &nU) ;
entries = Chv_entries(chv) ;
if ( CHV_IS_REAL(chv) ) {
int ii, iioff, ioff, istride, jj ;
A2_init(mtx, SPOOLES_REAL, nL, nD, nD, 1, NULL) ;
A2_zero(mtx) ;
ioff = nD + nL - 1 ;
istride = 2*nD + nU + nL - 2 ;
for ( ii = 0 ; ii < nD ; ii++ ) {
for ( jj = 0, iioff = ioff - nD + ii ;
jj < nL ;
jj++, iioff-- ) {
A2_setRealEntry(mtx, jj, ii, entries[iioff]) ;
}
ioff += istride ;
istride -= 2 ;
}
} else if ( CHV_IS_COMPLEX(chv) ) {
int ii, iioff, ioff, istride, jj ;
A2_init(mtx, SPOOLES_COMPLEX, nL, nD, nD, 1, NULL) ;
A2_zero(mtx) ;
ioff = nD + nL - 1 ;
istride = 2*nD + nU + nL - 2 ;
for ( ii = 0 ; ii < nD ; ii++ ) {
for ( jj = 0, iioff = ioff - nD + ii ;
jj < nL ;
jj++, iioff-- ) {
#if MYDEBUG > 0
fprintf(stdout,
"\n 12: ii %d, jj %d, iioff %d, setting (%d,%d) = (%20.12e,%20.12e)",
ii, jj, iioff, ii, ii, entries[2*iioff], entries[2*iioff+1]) ;
#endif
A2_setComplexEntry(mtx, jj, ii,
entries[2*iioff], entries[2*iioff+1]) ;
}
ioff += istride ;
istride -= 2 ;
}
}
return ; }
/*--------------------------------------------------------------------*/
syntax highlighted by Code2HTML, v. 0.9.1