/* factor.c */
#include "../Chv.h"
#define MYDEBUG 0
#define MYCHECK 0
/*--------------------------------------------------------------------*/
/*
------------------------------------------------
purpose -- to factor the front without pivoting
return value -- # of eliminated rows and columns
created -- 98aug27, cca
------------------------------------------------
*/
int
Chv_factorWithNoPivoting (
Chv *chv,
PatchAndGoInfo *info
) {
Chv wrkChv ;
int ncol, nD, nelim, nrow ;
int *colind, *rowind ;
/*
---------------
check the input
---------------
*/
if ( chv == NULL ) {
fprintf(stderr, "\n fatal error in Chv_factorWithNoPivoting()"
"\n bad input\n") ;
exit(-1) ;
}
nD = chv->nD ;
/*
--------------------------
set up the working chevron
--------------------------
*/
Chv_setDefaultFields(&wrkChv) ;
Chv_rowIndices(chv, &nrow, &rowind) ;
Chv_columnIndices(chv, &ncol, &colind) ;
Chv_initWithPointers(&wrkChv, chv->id, nD, chv->nL, chv->nU, chv->type,
chv->symflag, rowind, colind, Chv_entries(chv)) ;
/*
-------------------------------------
switch over the patch-and-go strategy
-------------------------------------
*/
if ( info == NULL ) {
/*
------------------
simple no pivoting
------------------
*/
nelim = 0 ;
while ( nelim < nD ) {
if ( Chv_r1upd(&wrkChv) == 0 ) {
break ;
}
Chv_shift(&wrkChv, 1) ;
nelim++ ;
}
} else if ( info->strategy == 1 ) {
double colmaxabs, diagabs, offmaxabs, rowmaxabs ;
/*
----------------------------------------
Patch-and-go for optimization matrices.
if |diag| < toosmall * max|offdiag| then
diag = 1.0
offdiag = 0.0
endif
----------------------------------------
*/
for ( nelim = 0 ; nelim < nD ; nelim++ ) {
Chv_maxabsInChevron(&wrkChv, 0, &diagabs, &rowmaxabs, &colmaxabs);
offmaxabs = (rowmaxabs >= colmaxabs) ? rowmaxabs : colmaxabs ;
if ( diagabs <= info->toosmall * offmaxabs ) {
if ( CHV_IS_REAL(chv) ) {
wrkChv.entries[0] = 1.0 ;
} else {
wrkChv.entries[0] = 1.0 ;
wrkChv.entries[1] = 0.0 ;
}
Chv_zeroOffdiagonalOfChevron(chv, 0) ;
if ( info->fudgeIV != NULL ) {
IV_push(info->fudgeIV, chv->colind[0]) ;
}
}
Chv_r1upd(&wrkChv) ;
Chv_shift(&wrkChv, 1) ;
}
} else if ( info->strategy == 2 ) {
double colmaxabs, diagabs, olddiag, newdiag, offmaxabs, rowmaxabs ;
/*
----------------------------------------------
Patch-and-go for structural analysis matrices.
if |diag| < fudge then
diag = fudge * max(max|offdiag|, 1.0)
endif
----------------------------------------------
*/
for ( nelim = 0 ; nelim < nD ; nelim++ ) {
Chv_maxabsInChevron(&wrkChv, 0, &diagabs, &rowmaxabs, &colmaxabs);
offmaxabs = (rowmaxabs >= colmaxabs) ? rowmaxabs : colmaxabs ;
if ( diagabs <= info->fudge ) {
olddiag = diagabs ;
if ( offmaxabs < 1.0 ) {
offmaxabs = 1.0 ;
}
if ( CHV_IS_REAL(chv) ) {
wrkChv.entries[0] = newdiag = info->fudge * offmaxabs ;
} else {
wrkChv.entries[0] = newdiag = info->fudge * offmaxabs ;
wrkChv.entries[1] = 0.0 ;
}
if ( info->fudgeIV != NULL ) {
IV_push(info->fudgeIV, chv->colind[0]) ;
}
if ( info->fudgeDV != NULL ) {
DV_push(info->fudgeDV, fabs(olddiag - newdiag)) ;
}
}
Chv_r1upd(&wrkChv) ;
Chv_shift(&wrkChv, 1) ;
}
} else {
fprintf(stderr, "\n fatal error in Chv_factorWithNoPivoting()"
"\n patch-and-go info != NULL, strategy = %d",
info->strategy) ;
exit(-1) ;
}
return(nelim) ; }
/*--------------------------------------------------------------------*/
/*
------------------------------------------------------------------
purpose -- factor the pivot chevron with pivoting
ndelay -- number of delayed rows and columns
pivotflag -- enable pivoting or not
0 --> no pivoting
1 --> enable pivoting
pivotsizesIV -- IV object that holds the sizes of the pivots,
used only when the front is symmetric or hermitian
and pivoting is enabled
workDV -- DV object used for working storage, resized as necessary
tau -- upper bound on the magnitude of the entries
in the factors, used only when pivoting is enabled
pntest -- pointer to be incremented with the number of pivot tests
return value -- # of eliminated rows and columns
created -- 98aug27, cca
------------------------------------------------------------------
*/
int
Chv_factorWithPivoting (
Chv *chv,
int ndelay,
int pivotflag,
IV *pivotsizesIV,
DV *workDV,
double tau,
int *pntest
) {
Chv wrkChv ;
int irow, jcol, ncol, nD, nelim, nrow, pivotsize ;
int *colind, *rowind ;
/*
---------------
check the input
---------------
*/
if ( chv == NULL || pivotflag != 1 || ndelay < 0 ) {
fprintf(stderr, "\n fatal error in Chv_factorWithPivoting()"
"\n bad input\n") ;
exit(-1) ;
}
if ( workDV == NULL ) {
fprintf(stderr, "\n fatal error in Chv_factorWithPivoting()"
"\n workDV is NULL \n") ;
exit(-1) ;
}
if ( tau < 1.0 ) {
fprintf(stderr, "\n fatal error in Chv_factorWithPivoting()"
"\n tau = %f is invalid \n", tau) ;
exit(-1) ;
}
if ( CHV_IS_REAL(chv) && CHV_IS_SYMMETRIC(chv)
&& pivotsizesIV == NULL ) {
fprintf(stderr, "\n fatal error in Chv_factorWithPivoting()"
"\n real symmetric front"
"\n pivoting enabled, pivotsizesIV is NULL\n") ;
exit(-1) ;
}
if ( CHV_IS_COMPLEX(chv)
&& (CHV_IS_SYMMETRIC(chv) || CHV_IS_HERMITIAN(chv))
&& pivotsizesIV == NULL ) {
fprintf(stderr, "\n fatal error in Chv_factorWithPivoting()"
"\n complex symmetric or hermitian front"
"\n pivoting enabled, pivotsizesIV is NULL\n") ;
exit(-1) ;
}
nD = chv->nD ;
/*
--------------------------
set up the working chevron
--------------------------
*/
Chv_setDefaultFields(&wrkChv) ;
Chv_rowIndices(chv, &nrow, &rowind) ;
Chv_columnIndices(chv, &ncol, &colind) ;
Chv_initWithPointers(&wrkChv, chv->id, nD, chv->nL, chv->nU, chv->type,
chv->symflag, rowind, colind, Chv_entries(chv)) ;
#if MYDEBUG > 0
fprintf(stdout, "\n\n after initializing wrkChv") ;
Chv_writeForHumanEye(&wrkChv, stdout) ;
fflush(stdout) ;
#endif
/*
-----------------------------
switch over the symmetry flag
-----------------------------
*/
if ( CHV_IS_SYMMETRIC(chv) || CHV_IS_HERMITIAN(chv) ) {
#if MYDEBUG > 0
fprintf(stdout,
"\n\n pivoting, hermitian or symmetric front %d", chv->id) ;
fflush(stdout) ;
#endif
/*
-------------------------
symmetric structure front
-------------------------
*/
IV_setSize(pivotsizesIV, 0) ;
nelim = 0 ;
while ( nelim < nD ) {
/*
-------------------------
find the 1x1 or 2x2 pivot
-------------------------
*/
#if MYDEBUG > 0
fprintf(stdout,
"\n trying to find pivot, nelim = %d, nD = %d, ndelay = %d",
nelim, nD, ndelay) ;
fflush(stdout) ;
#endif
pivotsize = Chv_findPivot(&wrkChv, workDV, tau, ndelay,
&irow, &jcol, pntest) ;
if ( irow > jcol ) {
int itemp = irow ;
irow = jcol ;
jcol = itemp ;
}
#if MYDEBUG > 0
fprintf(stdout,
"\n pivotsize = %d, local irow = %d, local jcol = %d",
pivotsize, irow, jcol) ;
fflush(stdout) ;
#endif
irow += nelim ;
jcol += nelim ;
if ( pivotsize == 0 ) {
/*
---------------------------------
no pivot found, break out of loop
---------------------------------
*/
break ;
} else {
ndelay = 0 ;
if ( irow == jcol ) {
/*
------------------------------------------------------
1x1 pivot found, swap row and column, update and shift
------------------------------------------------------
*/
#if MYDEBUG > 0
fprintf(stdout, "\n\n before swaps") ;
Chv_writeForHumanEye(chv, stdout) ;
fflush(stdout) ;
#endif
Chv_swapRowsAndColumns(chv, nelim, irow) ;
#if MYDEBUG > 0
fprintf(stdout, "\n\n after swaps") ;
Chv_writeForHumanEye(chv, stdout) ;
fflush(stdout) ;
#endif
#if MYDEBUG > 0
fprintf(stdout, "\n\n nelim = %d, before update", nelim) ;
Chv_writeForHumanEye(&wrkChv, stdout) ;
fflush(stdout) ;
#endif
if ( Chv_r1upd(&wrkChv) == 0 ) {
break ;
}
#if MYDEBUG > 0
fprintf(stdout, "\n\n nelim = %d, after update", nelim) ;
Chv_writeForHumanEye(&wrkChv, stdout) ;
fflush(stdout) ;
#endif
Chv_shift(&wrkChv, 1) ;
nelim++ ;
IV_push(pivotsizesIV, 1) ;
} else {
/*
--------------------------------------------------------
2x2 pivot found, swap rows and columns, update and shift
--------------------------------------------------------
*/
#if MYDEBUG > 0
fprintf(stdout, "\n\n before swaps") ;
Chv_writeForHumanEye(chv, stdout) ;
fflush(stdout) ;
#endif
Chv_swapRowsAndColumns(chv, nelim, irow) ;
Chv_swapRowsAndColumns(chv, nelim+1, jcol) ;
#if MYDEBUG > 0
fprintf(stdout, "\n\n after swaps") ;
Chv_writeForHumanEye(chv, stdout) ;
fflush(stdout) ;
#endif
#if MYDEBUG > 0
fprintf(stdout, "\n\n irow = %d, jcol = %d", irow, jcol) ;
fprintf(stdout, "\n\n nelim = %d, before update", nelim) ;
#endif
#if MYDEBUG > 0
Chv_writeForHumanEye(&wrkChv, stdout) ;
fflush(stdout) ;
#endif
if ( Chv_r2upd(&wrkChv) == 0 ) {
break ;
}
#if MYDEBUG > 0
fprintf(stdout, "\n\n nelim = %d, after update", nelim) ;
Chv_writeForHumanEye(&wrkChv, stdout) ;
fflush(stdout) ;
#endif
Chv_shift(&wrkChv, 2) ;
nelim += 2 ;
IV_push(pivotsizesIV, 2) ;
}
}
#if MYDEBUG > 0
fprintf(stdout, "\n\n ok, done with this pivot") ;
fflush(stdout) ;
#endif
}
} else {
/*
------------------
nonsymmetric front
------------------
*/
nelim = 0 ;
while ( nelim < nD ) {
/*
------------------
find the 1x1 pivot
------------------
*/
pivotsize = Chv_findPivot(&wrkChv, workDV, tau, ndelay,
&irow, &jcol, pntest) ;
irow += nelim ;
jcol += nelim ;
#if MYDEBUG > 0
fprintf(stdout, "\n\n irow = %d, jcol = %d", irow, jcol) ;
fflush(stdout) ;
#endif
if ( pivotsize == 0 ) {
/*
---------------------------------
no pivot found, break out of loop
---------------------------------
*/
break ;
} else {
ndelay = 0 ;
/*
------------------------------------------------------
1x1 pivot found, swap row and column, update and shift
------------------------------------------------------
*/
#if MYDEBUG > 1
fprintf(stdout, "\n\n before swaps") ;
Chv_writeForHumanEye(chv, stdout) ;
fflush(stdout) ;
#endif
Chv_swapRows(chv, nelim, irow) ;
Chv_swapColumns(chv, nelim, jcol) ;
#if MYDEBUG > 1
fprintf(stdout, "\n\n after swaps") ;
Chv_writeForHumanEye(chv, stdout) ;
fflush(stdout) ;
#endif
#if MYDEBUG > 0
fprintf(stdout, "\n\n nelim = %d, before update", nelim) ;
fflush(stdout) ;
#endif
#if MYDEBUG > 1
Chv_writeForHumanEye(&wrkChv, stdout) ;
fflush(stdout) ;
#endif
if ( Chv_r1upd(&wrkChv) == 0 ) {
break ;
}
#if MYDEBUG > 0
fprintf(stdout, "\n\n nelim = %d, after update", nelim) ;
#endif
#if MYDEBUG > 1
Chv_writeForHumanEye(&wrkChv, stdout) ;
fflush(stdout) ;
#endif
Chv_shift(&wrkChv, 1) ;
nelim++ ;
}
}
}
return(nelim) ; }
/*--------------------------------------------------------------------*/
static int symmetric1x1 ( Chv *chv ) ;
static int nonsym1x1 ( Chv *chv) ;
static int symmetric2x2 ( Chv *chv ) ;
/*--------------------------------------------------------------------*/
/*
---------------------------------------------------------
perform a rank one update using the first row and column.
this is used in the (L + I)D(I + U) factorization
return code ---
0 if the pivot was zero
1 if the pivot was nonzero
created -- 98jan23, cca
---------------------------------------------------------
*/
int
Chv_r1upd (
Chv *chv
) {
int rc = 0 ;
if ( chv == NULL ) {
fprintf(stderr, "\n fatal error in Chv_r1upd(%p)"
"\n bad input\n", chv) ;
exit(-1) ;
}
if ( CHV_IS_NONSYMMETRIC(chv) ) {
rc = nonsym1x1(chv) ;
} else if ( CHV_IS_SYMMETRIC(chv) || CHV_IS_HERMITIAN(chv) ) {
rc = symmetric1x1(chv) ;
} else {
fprintf(stderr, "\n fatal error in Chv_r1upd(%p)"
"\n symflag = %d\n", chv, chv->symflag) ;
exit(-1) ;
}
return(rc) ; }
/*--------------------------------------------------------------------*/
/*
------------------------------------------------------------------
perform a rank two update using the first two rows.
used in the (U^T + I)D(I + U) and (U^H + I)D(I + U) factorizations
return code ---
0 if the pivot was zero
1 if the pivot was nonzero
created -- 98jan23, cca
------------------------------------------------------------------
*/
int
Chv_r2upd (
Chv *chv
) {
int rc = 0 ;
if ( chv == NULL ) {
fprintf(stderr, "\n fatal error in Chv_r2upd(%p)"
"\n bad input\n", chv) ;
exit(-1) ;
}
if ( !(CHV_IS_SYMMETRIC(chv) || CHV_IS_HERMITIAN(chv)) ) {
fprintf(stderr, "\n fatal error in Chv_r2upd(%p)"
"\n symflag = %d\n", chv, chv->symflag) ;
exit(-1) ;
}
rc = symmetric2x2(chv) ;
return(rc) ; }
/*--------------------------------------------------------------------*/
/*
------------------------------------------------------------
perform an internal rank-1 update for a symmetric chevron
return code ---
0 if the pivot was zero
1 if the pivot was nonzero
created -- 98jan23, cca
------------------------------------------------------------
*/
static int
symmetric1x1 (
Chv *chv
) {
double dimag, dreal, fac1, fac2, limag, lreal, uimag, ureal ;
double *entries ;
int dloc, dstride, kchv, nD, nL, nU, uijloc, ukbeg, usize ;
/*
-------------------------------------
get dimensions and pointer to entries
-------------------------------------
*/
Chv_dimensions(chv, &nD, &nL, &nU) ;
entries = Chv_entries(chv) ;
#if MYDEBUG > 0
fprintf(stdout, "\n nD = %d, nL = %d, nU = %d, entries = %p",
nD, nL, nU, entries) ;
#endif
/*
----------------------------------------------
dloc : offset to the first diagonal element
dstride : stride to next diagonal element
usize : size of first row in upper part
----------------------------------------------
*/
dloc = 0 ;
dstride = nD + nU ;
usize = nD + nU - 1 ;
#if MYDEBUG > 0
fprintf(stdout, "\n dloc = %d, dstride = %d, usize = %d",
dloc, dstride, usize) ;
#endif
/*
----------------------
check for a zero pivot
----------------------
*/
if ( CHV_IS_REAL(chv) ) {
dreal = entries[dloc] ;
if ( dreal == 0.0 ) {
return(0) ;
}
fac1 = 1./dreal ;
} else if ( CHV_IS_COMPLEX(chv) ) {
dreal = entries[2*dloc] ;
dimag = entries[2*dloc+1] ;
if ( dreal == 0.0 && dimag == 0.0 ) {
return(0) ;
}
#if MYDEBUG > 0
fprintf(stdout, "\n chv->id = %d : (dreal,dimag) = <%12.4e,%12.4e>",
chv->id, dreal, dimag) ;
fprintf(stdout, "\n (dreal,dimag) = <%12.4e,%12.4e>",
dreal, dimag) ;
#endif
/*
------------------------------
compute (fac1,fac2) = 1/d(0,0)
------------------------------
*/
if ( CHV_IS_HERMITIAN(chv) ) {
fac1 = 1./dreal ; fac2 = 0.0 ;
entries[2*dloc+1] = 0.0 ;
} else {
Zrecip(dreal, dimag, &fac1, &fac2) ;
}
#if MYDEBUG > 0
fprintf(stdout, "\n 1/(dreal,dimag) = <%12.4e,%12.4e>",
fac1, fac2) ;
#endif
}
/*
------------------------
scale the first row of U
------------------------
*/
if ( CHV_IS_REAL(chv) ) {
DVscale(usize, &entries[1], fac1) ;
} else if ( CHV_IS_HERMITIAN(chv) ) {
DVscale(2*usize, &entries[2], fac1) ;
} else {
ZVscale(usize, &entries[2], fac1, fac2) ;
}
/*
------------------------------------
loop over the following chevrons
uijloc -- offset into uij multiplier
------------------------------------
*/
uijloc = dloc + 1 ;
for ( kchv = 1 ; kchv < nD ; kchv++ ) {
/*
--------------------------------------------------
dloc now points to next diagonal location
ukbeg -- offset into start of row in upper part
--------------------------------------------------
*/
dloc += dstride ;
ukbeg = dloc + 1 ;
#if MYDEBUG > 0
fprintf(stdout, "\n kchv = %5d, dloc = %5d"
"\n uijloc = %5d, usize = %5d, ukbeg = %5d",
kchv, dloc, uijloc, usize, ukbeg) ;
#endif
/*
------------------------------------
pull out the multiplier coefficients
------------------------------------
*/
if ( CHV_IS_REAL(chv) ) {
ureal = entries[uijloc] ;
lreal = dreal*ureal ;
} else if ( CHV_IS_COMPLEX(chv) ) {
ureal = entries[2*uijloc] ;
uimag = entries[2*uijloc+1] ;
if (CHV_IS_HERMITIAN(chv) ) {
lreal = dreal*ureal ;
limag = -dreal*uimag ;
} else {
lreal = dreal*ureal - dimag*uimag ;
limag = dreal*uimag + dimag*ureal ;
}
#if MYDEBUG > 0
fprintf(stdout,
"\n (lreal,limag) = <%12.4e,%12.4e>"
"\n (ureal,uimag) = <%12.4e,%12.4e>",
lreal, limag, ureal, uimag) ;
#endif
}
/*
-------------------------------------------------
update the upper row including the diagonal entry
-------------------------------------------------
*/
if ( CHV_IS_REAL(chv) ) {
DVaxpy(usize, &entries[dloc], -lreal, &entries[uijloc]) ;
} else if ( CHV_IS_COMPLEX(chv) ) {
ZVaxpy(usize, &entries[2*dloc], -lreal, -limag,
&entries[2*uijloc]) ;
}
/*
----------------------------------
adjust offsets and diagonal stride
----------------------------------
*/
uijloc++ ; dstride-- ; usize-- ;
}
return(1) ; }
/*--------------------------------------------------------------------*/
/*
------------------------------------------------------------
perform an internal rank-1 update for a nonsymmetric chevron
return code ---
0 if the pivot was zero
1 if the pivot was nonzero
created -- 98jan23, cca
------------------------------------------------------------
*/
static int
nonsym1x1 (
Chv *chv
) {
double dimag, dreal, fac1, fac2, limag, lreal, uimag, ureal ;
double *entries ;
int dloc, dstride, kchv, ljiloc, lkbeg, lsize, nD, nL, nU,
uijloc, ukbeg, usize ;
/*
-------------------------------------
get dimensions and pointer to entries
-------------------------------------
*/
Chv_dimensions(chv, &nD, &nL, &nU) ;
entries = Chv_entries(chv) ;
#if MYDEBUG > 0
fprintf(stdout, "\n nD = %d, nL = %d, nU = %d, entries = %p",
nD, nL, nU, entries) ;
#endif
/*
----------------------------------------------
dloc : offset to the first diagonal element
dstride : stride to next diagonal element
lsize : size of first column in lower part
usize : size of first row in upper part
----------------------------------------------
*/
dloc = nD + nL - 1 ;
dstride = 2*nD + nL + nU - 2 ;
lsize = nD + nL - 1 ;
usize = nD + nU - 1 ;
#if MYDEBUG > 0
fprintf(stdout, "\n dloc = %d, dstride = %d, lsize = %d, usize = %d",
dloc, dstride, lsize, usize) ;
#endif
/*
----------------------
check for a zero pivot
----------------------
*/
if ( CHV_IS_REAL(chv) ) {
dreal = entries[dloc] ;
if ( dreal == 0.0 ) {
return(0) ;
}
} else if ( CHV_IS_COMPLEX(chv) ) {
dreal = entries[2*dloc] ;
dimag = entries[2*dloc+1] ;
if ( dreal == 0.0 && dimag == 0.0 ) {
return(0) ;
}
#if MYDEBUG > 0
fprintf(stdout, "\n (dreal,dimag) = <%12.4e,%12.4e>",
dreal, dimag) ;
#endif
}
/*
-----------------------------------------
compute the inverse of the diagonal pivot
real: fac1 = 1/d(0,0)
complex: (fac1,fac2) = 1/d(0,0)
-----------------------------------------
*/
if ( CHV_IS_REAL(chv) ) {
fac1 = 1./dreal ;
} else if ( CHV_IS_COMPLEX(chv) ) {
Zrecip(dreal, dimag, &fac1, &fac2) ;
#if MYDEBUG > 0
fprintf(stdout, "\n 1/(dreal,dimag) = <%12.4e,%12.4e>",
fac1, fac2) ;
#endif
}
/*
---------------------------
scale the first column of L
(fac1,fac2) = 1/d(0,0)
---------------------------
*/
if ( CHV_IS_REAL(chv) ) {
DVscale(lsize, entries, fac1) ;
} else if ( CHV_IS_COMPLEX(chv) ) {
ZVscale(lsize, entries, fac1, fac2) ;
#if MYDEBUG > 2
{ double real, imag ;
int irow ;
fprintf(stdout, "\n entries in L after scaling") ;
for ( irow = 1 ; irow < nD + nL ; irow++ ) {
Chv_entry(chv, irow, 0, &real, &imag) ;
fprintf(stdout, "\n %% A(%d,%d) = %20.12e + %20.12ei",
irow, 0, real, imag) ;
}
}
#endif
}
/*
------------------------------------
loop over the following chevrons
ljiloc -- offset into lij multiplier
uijloc -- offset into uij multiplier
------------------------------------
*/
ljiloc = dloc - 1 ;
uijloc = dloc + 1 ;
for ( kchv = 1 ; kchv < nD ; kchv++ ) {
/*
--------------------------------------------------
dloc now points to next diagonal location
lsize and usize decremented
lkbeg -- offset into start of column in lower part
ukbeg -- offset into start of row in upper part
--------------------------------------------------
*/
dloc += dstride ;
lsize-- ;
usize-- ;
lkbeg = dloc - lsize ;
ukbeg = dloc + 1 ;
#if MYDEBUG > 0
fprintf(stdout,
"\n kchv = %5d, dloc = %5d"
"\n ljiloc = %5d, uijloc = %5d"
"\n lsize = %5d, usize = %5d"
"\n lkbeg = %5d, ukbeg = %5d",
kchv, dloc, ljiloc, uijloc, lsize, usize, lkbeg, ukbeg) ;
#endif
/*
------------------------------------
pull out the multiplier coefficients
------------------------------------
*/
if ( CHV_IS_REAL(chv) ) {
lreal = entries[ljiloc] ;
ureal = entries[uijloc] ;
} else if ( CHV_IS_COMPLEX(chv) ) {
lreal = entries[2*ljiloc] ;
limag = entries[2*ljiloc+1] ;
ureal = entries[2*uijloc] ;
uimag = entries[2*uijloc+1] ;
#if MYDEBUG > 0
fprintf(stdout,
"\n (lreal,limag) = <%12.4e,%12.4e>"
"\n (ureal,uimag) = <%12.4e,%12.4e>",
lreal, limag, ureal, uimag) ;
#endif
}
/*
-------------------------
update the diagonal entry
-------------------------
*/
if ( CHV_IS_REAL(chv) ) {
entries[dloc] -= lreal*ureal ;
} else if ( CHV_IS_COMPLEX(chv) ) {
entries[2*dloc] -= lreal*ureal - limag*uimag ;
entries[2*dloc+1] -= lreal*uimag + limag*ureal ;
}
/*
-----------------------
update the lower column
-----------------------
*/
if ( CHV_IS_REAL(chv) ) {
DVaxpy(lsize, &entries[lkbeg], -ureal, entries) ;
} else if ( CHV_IS_COMPLEX(chv) ) {
ZVaxpy(lsize, &entries[2*lkbeg], -ureal, -uimag, entries) ;
}
/*
--------------------
update the upper row
--------------------
*/
if ( CHV_IS_REAL(chv) ) {
DVaxpy(usize, &entries[ukbeg], -lreal, &entries[uijloc+1]) ;
} else if ( CHV_IS_COMPLEX(chv) ) {
ZVaxpy(usize, &entries[2*ukbeg],
-lreal, -limag, &entries[2*uijloc+2]) ;
}
/*
----------------------------------
adjust offsets and diagonal stride
----------------------------------
*/
ljiloc-- ; uijloc++ ; dstride -= 2 ;
}
/*
------------------
scale the row of U
------------------
*/
usize = nD + nU - 1 ;
if ( CHV_IS_REAL(chv) ) {
DVscale(usize, &entries[nD+nL], fac1) ;
} else if ( CHV_IS_COMPLEX(chv) ) {
ZVscale(usize, &entries[2*(nD+nL)], fac1, fac2) ;
}
return(1) ; }
/*--------------------------------------------------------------------*/
/*
-------------------------------------
perform an internal rank-2 update for
a hermitian or symmetric chevron
return code ---
0 if the pivot was zero
1 if the pivot was nonzero
created -- 98jan23, cca
-------------------------------------
*/
static int
symmetric2x2 (
Chv *chv
) {
double areal, aimag, breal, bimag, creal, cimag,
ereal, eimag, freal, fimag, greal, gimag,
l0imag, l1imag, l0real, l1real,
u0imag, u1imag, u0real, u1real ;
double *entries ;
int dloc, dstride, kchv, nD, nL, nU, rc,
u0jloc, u1jloc, ukbeg, usize ;
/*
-------------------------------------
get dimensions and pointer to entries
-------------------------------------
*/
Chv_dimensions(chv, &nD, &nL, &nU) ;
entries = Chv_entries(chv) ;
#if MYDEBUG > 0
fprintf(stdout, "\n nD = %d, nL = %d, nU = %d, entries = %p",
nD, nL, nU, entries) ;
#endif
/*
----------------------------------------
check for a zero pivot
D = [ a b ] for hermitian
[ conj(b) c ]
D = [ a b ] for symmetric
[ b c ]
compute the inverse of D
E = inv(D) = [ e f ] for hermitian
[ conj(f) g ]
E = inv(D) = [ e f ] for symmetric
[ f g ]
----------------------------------------
*/
if ( CHV_IS_REAL(chv) ) {
double denom ;
areal = entries[0] ;
breal = entries[1] ;
creal = entries[nD+nU] ;
if ( (denom = areal*creal - breal*breal) == 0.0 ) {
rc = 0 ;
} else {
rc = 1 ;
denom = 1./denom ;
ereal = creal*denom ;
freal = -breal*denom ;
greal = areal*denom ;
}
} else if ( CHV_IS_COMPLEX(chv) ) {
areal = entries[0] ;
aimag = entries[1] ;
breal = entries[2] ;
bimag = entries[3] ;
creal = entries[2*(nD+nU)] ;
cimag = entries[2*(nD+nU)+1] ;
#if MYDEBUG > 0
if ( CHV_IS_HERMITIAN(chv) ) {
fprintf(stdout,
"\n hermitian D = [ <%12.4e,%12.4e> <%12.4e,%12.4e> ]"
"\n [ <%12.4e,%12.4e> <%12.4e,%12.4e> ]",
areal, aimag, breal, bimag,
breal, -bimag, creal, cimag) ;
} else {
fprintf(stdout,
"\n symmetric D = [ <%12.4e,%12.4e> <%12.4e,%12.4e> ]"
"\n [ <%12.4e,%12.4e> <%12.4e,%12.4e> ]",
areal, aimag, breal, bimag,
breal, bimag, creal, cimag) ;
}
#endif
if ( CHV_IS_HERMITIAN(chv) ) {
rc = Zrecip2(areal, 0.0, breal, bimag,
breal, -bimag, creal, 0.0,
&ereal, NULL, &freal, &fimag,
NULL, NULL, &greal, NULL) ;
eimag = gimag = 0.0 ;
} else {
rc = Zrecip2(areal, aimag, breal, bimag,
breal, bimag, creal, cimag,
&ereal, &eimag, &freal, &fimag,
NULL, NULL, &greal, &gimag) ;
}
} else {
fprintf(stderr, "\n fatal error in Chv_symmetric2x2"
"\n chevron must be real or complex") ;
exit(-1) ;
}
if ( rc == 0 ) {
/*
-------------------------
pivot is singular, return
-------------------------
*/
return(0) ;
}
#if MYDEBUG > 0
if ( CHV_IS_HERMITIAN(chv) ) {
fprintf(stdout,
"\n hermitian DINV = [ <%12.4e,%12.4e> <%12.4e,%12.4e> ]"
"\n [ <%12.4e,%12.4e> <%12.4e,%12.4e> ]",
ereal, eimag, freal, fimag,
freal, -fimag, greal, gimag) ;
} else {
fprintf(stdout,
"\n symmetric DINV = [ <%12.4e,%12.4e> <%12.4e,%12.4e> ]"
"\n [ <%12.4e,%12.4e> <%12.4e,%12.4e> ]",
ereal, eimag, freal, fimag,
freal, fimag, greal, gimag) ;
}
#endif
/*
-----------------------------
scale the first two rows of U
-----------------------------
*/
u0jloc = 2 ;
u1jloc = nD + nU + 1 ;
usize = nD + nU - 2 ;
if ( CHV_IS_REAL(chv) ) {
DVscale2(usize, &entries[u0jloc], &entries[u1jloc],
ereal, freal, freal, greal) ;
} else if ( CHV_IS_HERMITIAN(chv) ) {
ZVscale2(usize, &entries[2*u0jloc], &entries[2*u1jloc],
ereal, 0.0, freal, fimag, freal, -fimag, greal, 0.0) ;
} else {
ZVscale2(usize, &entries[2*u0jloc], &entries[2*u1jloc],
ereal, eimag, freal, fimag, freal, fimag, greal, gimag) ;
}
#if MYDEBUG > 2
{ double real, imag ;
int irow, jcol ;
fprintf(stdout, "\n entries in U after scaling") ;
for ( irow = 0 ; irow <= 1 ; irow++ ) {
for ( jcol = 2 ; jcol < nD + nU ; jcol++ ) {
Chv_entry(chv, 0, jcol, &real, &imag) ;
fprintf(stdout, "\n %% A(%d,%d) = %20.12e + %20.12ei",
0, jcol, real, imag) ;
}
}
}
#endif
/*
------------------------------------
loop over the following chevrons
u0jloc -- offset into u0j multiplier
u1jloc -- offset into u1j multiplier
------------------------------------
*/
usize = nD + nU - 2 ;
dloc = nD + nU ;
dstride = nD + nU - 1 ;
for ( kchv = 2 ; kchv < nD ; kchv++ ) {
/*
--------------------------------------------------
dloc now points to next diagonal location
ukbeg -- offset into start of row in upper part
--------------------------------------------------
*/
dloc += dstride ;
ukbeg = dloc + 1 ;
#if MYDEBUG > 0
fprintf(stdout, "\n kchv = %5d, dloc = %5d"
"\n u0jloc = %5d, u1jloc = %d, usize = %5d, ukbeg = %5d",
kchv, dloc, u0jloc, u1jloc, usize, ukbeg) ;
#endif
/*
------------------------------------
pull out the multiplier coefficients
------------------------------------
*/
if ( CHV_IS_REAL(chv) ) {
u0real = entries[u0jloc] ;
u1real = entries[u1jloc] ;
l0real = u0real*areal + u1real*breal ;
l1real = u0real*breal + u1real*creal ;
} else if ( CHV_IS_COMPLEX(chv) ) {
u0real = entries[2*u0jloc] ;
u0imag = entries[2*u0jloc+1] ;
u1real = entries[2*u1jloc] ;
u1imag = entries[2*u1jloc+1] ;
if ( CHV_IS_HERMITIAN(chv) ) {
l0real = u0real*areal + u1real*breal - u1imag*bimag ;
l0imag = -u0imag*areal - u1real*bimag - u1imag*breal ;
l1real = u0real*breal + u0imag*bimag + u1real*creal ;
l1imag = u0real*bimag - u0imag*breal - u1imag*creal ;
} else {
l0real = u0real*areal - u0imag*aimag
+ u1real*breal - u1imag*bimag ;
l0imag = u0real*aimag + u0imag*areal
+ u1real*bimag + u1imag*breal ;
l1real = u0real*breal - u0imag*bimag
+ u1real*creal - u1imag*cimag ;
l1imag = u0real*bimag + u0imag*breal
+ u1real*cimag + u1imag*creal ;
}
#if MYDEBUG > 0
fprintf(stdout,
"\n (l0real,l0imag) = <%12.4e,%12.4e>"
"\n (l1real,l1imag) = <%12.4e,%12.4e>",
l0real, l0imag, l1real, l1imag) ;
#endif
}
/*
-------------------------------------------------
update the upper row including the diagonal entry
-------------------------------------------------
*/
if ( CHV_IS_REAL(chv) ) {
DVaxpy2(usize, &entries[dloc], -l0real, &entries[u0jloc],
-l1real, &entries[u1jloc]) ;
} else if ( CHV_IS_COMPLEX(chv) ) {
ZVaxpy2(usize, &entries[2*dloc],
-l0real, -l0imag, &entries[2*u0jloc],
-l1real, -l1imag, &entries[2*u1jloc]) ;
}
/*
----------------------------------
adjust offsets and diagonal stride
----------------------------------
*/
u0jloc++ ; u1jloc++ ; dstride-- ; usize-- ;
}
return(1) ; }
/*--------------------------------------------------------------------*/
/*
------------------------------------------------------------------
purpose -- looking at just a single chevron inside the Chv object,
find the absolute value of the diagonal element, and
the maximum absolute values of the offdiagonal elements
in the chevron's row and column.
created -- 98aug26, cca
------------------------------------------------------------------
*/
void
Chv_maxabsInChevron (
Chv *chv,
int ichv,
double *pdiagmaxabs,
double *prowmaxabs,
double *pcolmaxabs
) {
double *pdiag ;
int length, loc, nD, nL, nU ;
/*
---------------
check the input
---------------
*/
if ( chv == NULL || ichv < 0 || ichv >= chv->nD
|| pdiagmaxabs == NULL || prowmaxabs == NULL || pcolmaxabs == NULL ) {
fprintf(stderr, "\n fatal error in Chv_maxabsInChevron()"
"\n bad input\n") ;
exit(-1) ;
}
Chv_dimensions(chv, &nD, &nL, &nU) ;
pdiag = Chv_diagLocation(chv, ichv) ;
length = nD - ichv - 1 + nU ;
if ( CHV_IS_REAL(chv) ) {
if ( CHV_IS_SYMMETRIC(chv) ) {
*pdiagmaxabs = fabs(*pdiag) ;
*prowmaxabs = DVmaxabs(length, pdiag + 1, &loc) ;
*pcolmaxabs = *prowmaxabs ;
} else if ( CHV_IS_NONSYMMETRIC(chv) ) {
*pdiagmaxabs = fabs(*pdiag) ;
*prowmaxabs = DVmaxabs(length, pdiag + 1, &loc) ;
*pcolmaxabs = DVmaxabs(length, pdiag - length, &loc) ;
} else {
fprintf(stderr, "\n fatal error in Chv_maxabsInChevron()"
"\n chv is real, chv->symflag = %d"
"\n must be symmetric or nonsymmetric\n", chv->symflag) ;
exit(-1) ;
}
} else if ( CHV_IS_COMPLEX(chv) ) {
if ( CHV_IS_SYMMETRIC(chv) || CHV_IS_HERMITIAN(chv) ) {
*pdiagmaxabs = Zabs(*pdiag, *(pdiag+1)) ;
*prowmaxabs = ZVmaxabs(length, pdiag + 2) ;
*pcolmaxabs = *prowmaxabs ;
} else if ( CHV_IS_NONSYMMETRIC(chv) ) {
*pdiagmaxabs = Zabs(*pdiag, *(pdiag+1)) ;
*prowmaxabs = ZVmaxabs(length, pdiag + 2) ;
*pcolmaxabs = ZVmaxabs(length, pdiag - 2*length) ;
} else {
fprintf(stderr, "\n fatal error in Chv_maxabsInChevron()"
"\n chv is complex, chv->symflag = %d"
"\n must be symmetric or nonsymmetric\n", chv->symflag) ;
exit(-1) ;
}
} else {
fprintf(stderr, "\n fatal error in Chv_maxabsInChevron()"
"\n chv->type = %d, must be real or complex\n", chv->type) ;
exit(-1) ;
}
return ; }
/*--------------------------------------------------------------------*/
/*
-------------------------------------------------------
purpose -- zero the offdiagonal entries of chevron ichv
created -- 98aug26, cca
-------------------------------------------------------
*/
void
Chv_zeroOffdiagonalOfChevron (
Chv *chv,
int ichv
) {
double *pdiag ;
int length, nD, nL, nU ;
/*
---------------
check the input
---------------
*/
if ( chv == NULL || ichv < 0 || ichv >= chv->nD ) {
fprintf(stderr, "\n fatal error in Chv_zeroOffdiagonalOfChevron()"
"\n bad input\n") ;
exit(-1) ;
}
Chv_dimensions(chv, &nD, &nL, &nU) ;
pdiag = Chv_diagLocation(chv, ichv) ;
length = nD - ichv - 1 + nU ;
if ( CHV_IS_REAL(chv) ) {
if ( CHV_IS_SYMMETRIC(chv) ) {
DVzero(length, pdiag+1) ;
} else if ( CHV_IS_NONSYMMETRIC(chv) ) {
DVzero(length, pdiag + 1) ;
DVzero(length, pdiag - length) ;
} else {
fprintf(stderr,
"\n fatal error in Chv_zeroOffdiagonalOfChevron()"
"\n chv is real, chv->symflag = %d"
"\n must be symmetric or nonsymmetric\n", chv->symflag) ;
exit(-1) ;
}
} else if ( CHV_IS_COMPLEX(chv) ) {
if ( CHV_IS_SYMMETRIC(chv) || CHV_IS_HERMITIAN(chv) ) {
ZVzero(length, pdiag+2) ;
} else if ( CHV_IS_NONSYMMETRIC(chv) ) {
ZVzero(length, pdiag+2) ;
ZVzero(length, pdiag-2*length) ;
} else {
fprintf(stderr,
"\n fatal error in Chv_zeroOffdiagonalOfChevron()"
"\n chv is complex, chv->symflag = %d"
"\n must be symmetric or nonsymmetric\n", chv->symflag) ;
exit(-1) ;
}
} else {
fprintf(stderr, "\n fatal error in Chv_zeroOffdiagonalOfChevron()"
"\n chv->type = %d, must be real or complex\n", chv->type) ;
exit(-1) ;
}
return ; }
/*--------------------------------------------------------------------*/
syntax highlighted by Code2HTML, v. 0.9.1