/* scalevec.c */
#include "../SubMtx.h"
/*--------------------------------------------------------------------*/
static void diagonal_scale3vec ( SubMtx *mtxA, double y0[], double y1[],
double y2[], double x0[], double x1[], double x2[] ) ;
static void diagonal_scale2vec ( SubMtx *mtxA, double y0[], double y1[],
double x0[], double x1[] ) ;
static void diagonal_scale1vec ( SubMtx *mtxA,
double y0[], double x0[] ) ;
static void block_diagonal_sym_scale3vec ( SubMtx *mtxA, double y0[],
double y1[], double y2[], double x0[], double x1[], double x2[] ) ;
static void block_diagonal_sym_scale2vec ( SubMtx *mtxA, double y0[],
double y1[], double x0[], double x1[] ) ;
static void block_diagonal_sym_scale1vec ( SubMtx *mtxA, double y0[],
double x0[] ) ;
static void block_diagonal_herm_scale3vec ( SubMtx *mtxA, double y0[],
double y1[], double y2[], double x0[], double x1[], double x2[] ) ;
static void block_diagonal_herm_scale2vec ( SubMtx *mtxA, double y0[],
double y1[], double x0[], double x1[] ) ;
static void block_diagonal_herm_scale1vec ( SubMtx *mtxA, double y0[],
double x0[] ) ;
/*--------------------------------------------------------------------*/
/*
-----------------------------------
purpose -- compute [y0] := A * [x0]
created -- 98apr17, cca
-----------------------------------
*/
void
SubMtx_scale1vec (
SubMtx *mtxA,
double y0[],
double x0[]
) {
/*
---------------
check the input
---------------
*/
if ( mtxA == NULL || y0 == NULL || x0 == NULL ) {
fprintf(stderr, "\n fatal error in SubMtx_scale1vec(%p,%p,%p)"
"\n bad input\n", mtxA, y0, x0) ;
exit(-1) ;
}
if ( ! (SUBMTX_IS_REAL(mtxA) || SUBMTX_IS_COMPLEX(mtxA)) ) {
fprintf(stderr, "\n fatal error in SubMtx_scale1vec(%p,%p,%p)"
"\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
mtxA, y0, x0, mtxA->type) ;
exit(-1) ;
}
switch ( mtxA->mode ) {
case SUBMTX_DIAGONAL :
diagonal_scale1vec(mtxA, y0, x0) ;
break ;
case SUBMTX_BLOCK_DIAGONAL_SYM :
block_diagonal_sym_scale1vec(mtxA, y0, x0) ;
break ;
case SUBMTX_BLOCK_DIAGONAL_HERM :
if ( ! SUBMTX_IS_COMPLEX(mtxA) ) {
fprintf(stderr, "\n fatal error in SubMtx_scale1vec(%p,%p,%p)"
"\n hermitian matrix, type %d is not SPOOLES_COMPLEX\n",
mtxA, y0, x0, mtxA->type) ;
exit(-1) ;
}
block_diagonal_herm_scale1vec(mtxA, y0, x0) ;
break ;
default :
fprintf(stderr, "\n fatal error in SubMtx_scale1vec()"
"\n matrix mode not supported"
"\n must be SUBMTX_DIAGONAL,"
"\n or SUBMTX_BLOCK_DIAGONAL_SYM"
"\n or SUBMTX_BLOCK_DIAGONAL_HERM\n") ;
exit(-1) ;
}
return ; }
/*--------------------------------------------------------------------*/
/*
-------------------------------------
purpose -- compute [y0 y1] := [x0 x1]
created -- 98apr17, cca
-------------------------------------
*/
void
SubMtx_scale2vec (
SubMtx *mtxA,
double y0[],
double y1[],
double x0[],
double x1[]
) {
/*
---------------
check the input
---------------
*/
if ( mtxA == NULL || y0 == NULL || y1 == NULL
|| x0 == NULL || x1 == NULL ) {
fprintf(stderr, "\n fatal error in SubMtx_scale2vec(%p,%p,%p,%p,%p)"
"\n bad input\n", mtxA, y0, y1, x0, x1) ;
exit(-1) ;
}
if ( ! (SUBMTX_IS_REAL(mtxA) || SUBMTX_IS_COMPLEX(mtxA)) ) {
fprintf(stderr, "\n fatal error in SubMtx_scale2vec(%p,%p,%p,%p,%p)"
"\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
mtxA, y0, y1, x0, x1, mtxA->type) ;
exit(-1) ;
}
switch ( mtxA->mode ) {
case SUBMTX_DIAGONAL :
diagonal_scale2vec(mtxA, y0, y1, x0, x1) ;
break ;
case SUBMTX_BLOCK_DIAGONAL_SYM :
block_diagonal_sym_scale2vec(mtxA, y0, y1, x0, x1) ;
break ;
case SUBMTX_BLOCK_DIAGONAL_HERM :
if ( ! SUBMTX_IS_COMPLEX(mtxA) ) {
fprintf(stderr,
"\n fatal error in SubMtx_scale2vec(%p,%p,%p,%p,%p)"
"\n hermitian matrix, type %d is not SPOOLES_COMPLEX\n",
mtxA, y0, y1, x0, x1, mtxA->type) ;
exit(-1) ;
}
block_diagonal_herm_scale2vec(mtxA, y0, y1, x0, x1) ;
break ;
default :
fprintf(stderr, "\n fatal error in SubMtx_scale2vec()"
"\n matrix type not supported"
"\n must be SUBMTX_DIAGONAL,"
"\n or SUBMTX_BLOCK_DIAGONAL_SYM"
"\n or SUBMTX_BLOCK_DIAGONAL_HERM\n") ;
exit(-1) ;
}
return ; }
/*--------------------------------------------------------------------*/
/*
-----------------------------------------------
purpose -- compute [y0 y1 y2] := A * [x0 x1 x2]
created -- 98apr17, cca
-----------------------------------------------
*/
void
SubMtx_scale3vec (
SubMtx *mtxA,
double y0[],
double y1[],
double y2[],
double x0[],
double x1[],
double x2[]
) {
/*
---------------
check the input
---------------
*/
if ( mtxA == NULL || y0 == NULL || y1 == NULL || y2 == NULL
|| x0 == NULL || x1 == NULL || x2 == NULL ) {
fprintf(stderr,
"\n fatal error in SubMtx_scale3vec(%p,%p,%p,%p,%p,%p,%p)"
"\n bad input\n", mtxA, y0, y1, y2, x0, x1, x2) ;
exit(-1) ;
}
if ( ! (SUBMTX_IS_REAL(mtxA) || SUBMTX_IS_COMPLEX(mtxA)) ) {
fprintf(stderr,
"\n fatal error in SubMtx_scale3vec(%p,%p,%p,%p,%p,%p,%p)"
"\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
mtxA, y0, y1, y2, x0, x1, x2, mtxA->type) ;
exit(-1) ;
}
switch ( mtxA->mode ) {
case SUBMTX_DIAGONAL :
diagonal_scale3vec(mtxA, y0, y1, y2, x0, x1, x2) ;
break ;
case SUBMTX_BLOCK_DIAGONAL_SYM :
block_diagonal_sym_scale3vec(mtxA, y0, y1, y2, x0, x1, x2) ;
break ;
case SUBMTX_BLOCK_DIAGONAL_HERM :
if ( ! SUBMTX_IS_COMPLEX(mtxA) ) {
fprintf(stderr,
"\n fatal error in SubMtx_scale3vec(%p,%p,%p,%p,%p,%p,%p)"
"\n hermitian matrix, type %d is not SPOOLES_COMPLEX\n",
mtxA, y0, y1, y2, x0, x1, x2, mtxA->type) ;
exit(-1) ;
}
block_diagonal_herm_scale3vec(mtxA, y0, y1, y2, x0, x1, x2) ;
break ;
default :
fprintf(stderr, "\n fatal error in SubMtx_scale3vec()"
"\n matrix type not supported"
"\n must be SUBMTX_DIAGONAL,"
"\n or SUBMTX_BLOCK_DIAGONAL_SYM"
"\n or SUBMTX_BLOCK_DIAGONAL_HERM\n") ;
exit(-1) ;
}
return ; }
/*--------------------------------------------------------------------*/
/*
-----------------------------------------------
purpose -- compute [y0 y1 y2] := A * [x0 x1 x2]
where A is diagonal
created -- 98apr17, cca
-----------------------------------------------
*/
static void
diagonal_scale3vec (
SubMtx *mtxA,
double y0[],
double y1[],
double y2[],
double x0[],
double x1[],
double x2[]
) {
double *entriesA ;
int nrowA ;
SubMtx_diagonalInfo(mtxA, &nrowA, &entriesA) ;
if ( SUBMTX_IS_REAL(mtxA) ) {
double a ;
int irowA ;
for ( irowA = 0 ; irowA < nrowA ; irowA++ ) {
a = entriesA[irowA] ;
y0[irowA] = a*x0[irowA] ;
y1[irowA] = a*x1[irowA] ;
y2[irowA] = a*x2[irowA] ;
}
} else if ( SUBMTX_IS_COMPLEX(mtxA) ) {
double ai, ar, xi0, xi1, xi2, xr0, xr1, xr2 ;
int iloc, irowA, rloc ;
for ( irowA = rloc = 0, iloc = 1 ;
irowA < nrowA ;
irowA++, rloc += 2, iloc += 2 ) {
xr0 = x0[rloc] ; xi0 = x0[iloc] ;
xr1 = x1[rloc] ; xi1 = x1[iloc] ;
xr2 = x2[rloc] ; xi2 = x2[iloc] ;
ar = entriesA[rloc] ; ai = entriesA[iloc] ;
y0[rloc] = ar*xr0 - ai*xi0 ; y0[iloc] = ar*xi0 + ai*xr0 ;
y1[rloc] = ar*xr1 - ai*xi1 ; y1[iloc] = ar*xi1 + ai*xr1 ;
y2[rloc] = ar*xr2 - ai*xi2 ; y2[iloc] = ar*xi2 + ai*xr2 ;
}
}
return ; }
/*--------------------------------------------------------------------*/
/*
-----------------------------------------
purpose -- compute [y0 y1] := A * [x0 x1]
where A is diagonal
created -- 98apr17, cca
-----------------------------------------
*/
static void
diagonal_scale2vec (
SubMtx *mtxA,
double y0[],
double y1[],
double x0[],
double x1[]
) {
double *entriesA ;
int nrowA ;
SubMtx_diagonalInfo(mtxA, &nrowA, &entriesA) ;
if ( SUBMTX_IS_REAL(mtxA) ) {
double a ;
int irowA ;
for ( irowA = 0 ; irowA < nrowA ; irowA++ ) {
a = entriesA[irowA] ;
y0[irowA] = a*x0[irowA] ;
y1[irowA] = a*x1[irowA] ;
}
} else if ( SUBMTX_IS_COMPLEX(mtxA) ) {
double ai, ar, xi0, xi1, xr0, xr1 ;
int iloc, irowA, rloc ;
for ( irowA = rloc = 0, iloc = 1 ;
irowA < nrowA ;
irowA++, rloc += 2, iloc += 2 ) {
xr0 = x0[rloc] ; xi0 = x0[iloc] ;
xr1 = x1[rloc] ; xi1 = x1[iloc] ;
ar = entriesA[rloc] ; ai = entriesA[iloc] ;
y0[rloc] = ar*xr0 - ai*xi0 ; y0[iloc] = ar*xi0 + ai*xr0 ;
y1[rloc] = ar*xr1 - ai*xi1 ; y1[iloc] = ar*xi1 + ai*xr1 ;
}
}
return ; }
/*--------------------------------------------------------------------*/
/*
-----------------------------------
purpose -- compute [y0] := A * [x0]
where A is diagonal
created -- 98apr17, cca
-----------------------------------
*/
static void
diagonal_scale1vec (
SubMtx *mtxA,
double y0[],
double x0[]
) {
double *entriesA ;
int nrowA ;
SubMtx_diagonalInfo(mtxA, &nrowA, &entriesA) ;
if ( SUBMTX_IS_REAL(mtxA) ) {
double a ;
int irowA ;
for ( irowA = 0 ; irowA < nrowA ; irowA++ ) {
a = entriesA[irowA] ;
y0[irowA] = a*x0[irowA] ;
}
} else if ( SUBMTX_IS_COMPLEX(mtxA) ) {
double ai, ar, xi0, xr0 ;
int iloc, irowA, rloc ;
for ( irowA = rloc = 0, iloc = 1 ;
irowA < nrowA ;
irowA++, rloc += 2, iloc += 2 ) {
xr0 = x0[rloc] ; xi0 = x0[iloc] ;
ar = entriesA[rloc] ; ai = entriesA[iloc] ;
y0[rloc] = ar*xr0 - ai*xi0 ; y0[iloc] = ar*xi0 + ai*xr0 ;
}
}
return ; }
/*--------------------------------------------------------------------*/
/*
--------------------------------------------------
purpose -- compute [y0 y1 y2] := A * [x0 x1 x2]
where A is block diagonal and symmetric
created -- 98apr17, cca
--------------------------------------------------
*/
static void
block_diagonal_sym_scale3vec (
SubMtx *mtxA,
double y0[],
double y1[],
double y2[],
double x0[],
double x1[],
double x2[]
) {
double *entries ;
int nentA, nrowA ;
int *pivotsizes ;
SubMtx_blockDiagonalInfo(mtxA, &nrowA, &nentA, &pivotsizes, &entries) ;
if ( SUBMTX_IS_REAL(mtxA) ) {
int ipivot, irowA, kk ;
for ( irowA = ipivot = kk = 0 ; irowA < nrowA ; ipivot++ ) {
if ( pivotsizes[ipivot] == 1 ) {
double a00, x00, x01, x02 ;
a00 = entries[kk] ;
x00 = x0[irowA] ;
x01 = x1[irowA] ;
x02 = x2[irowA] ;
y0[irowA] = a00*x00 ;
y1[irowA] = a00*x01 ;
y2[irowA] = a00*x02 ;
kk++, irowA++ ;
} else if ( pivotsizes[ipivot] == 2 ) {
double a00, a01, a11,
x00, x01, x02, x10, x11, x12 ;
a00 = entries[kk] ;
a01 = entries[kk+1] ;
a11 = entries[kk+2] ;
x00 = x0[irowA] ;
x01 = x1[irowA] ;
x02 = x2[irowA] ;
x10 = x0[irowA+1] ;
x11 = x1[irowA+1] ;
x12 = x2[irowA+1] ;
y0[irowA] = a00*x00 + a01*x10 ;
y1[irowA] = a00*x01 + a01*x11 ;
y2[irowA] = a00*x02 + a01*x12 ;
y0[irowA+1] = a01*x00 + a11*x10 ;
y1[irowA+1] = a01*x01 + a11*x11 ;
y2[irowA+1] = a01*x02 + a11*x12 ;
kk += 3, irowA += 2 ;
} else {
fprintf(stderr, "\n fatal error in SubMtx_scale3vec()"
"\n pivotsizes[%d] = %d", ipivot, pivotsizes[ipivot]) ;
exit(-1) ;
}
}
} else if ( SUBMTX_IS_COMPLEX(mtxA) ) {
int iloc, ipivot, irowA, kk, rloc ;
for ( irowA = ipivot = kk = rloc = 0, iloc = 1 ;
irowA < nrowA ;
ipivot++ ) {
if ( pivotsizes[ipivot] == 1 ) {
double ai00, ar00, xi0, xi1, xi2, xr0, xr1, xr2 ;
ar00 = entries[kk] ; ai00 = entries[kk+1] ;
xr0 = x0[rloc] ; xi0 = x0[iloc] ;
xr1 = x1[rloc] ; xi1 = x1[iloc] ;
xr2 = x2[rloc] ; xi2 = x2[iloc] ;
y0[rloc] = ar00*xr0 - ai00*xi0; y0[iloc] = ar00*xi0 + ai00*xr0;
y1[rloc] = ar00*xr1 - ai00*xi1; y1[iloc] = ar00*xi1 + ai00*xr1;
y2[rloc] = ar00*xr2 - ai00*xi2; y2[iloc] = ar00*xi2 + ai00*xr2;
kk += 2, irowA++, rloc += 2, iloc += 2 ;
} else if ( pivotsizes[ipivot] == 2 ) {
double ai00, ai01, ai11, ar00, ar01, ar11,
xi00, xi01, xi02, xi10, xi11, xi12,
xr00, xr01, xr02, xr10, xr11, xr12 ;
int iloc0, iloc1, rloc0, rloc1 ;
ar00 = entries[kk] ; ai00 = entries[kk+1] ;
ar01 = entries[kk+2] ; ai01 = entries[kk+3] ;
ar11 = entries[kk+4] ; ai11 = entries[kk+5] ;
rloc0 = rloc ; iloc0 = iloc ;
rloc1 = rloc + 2 ; iloc1 = iloc + 2 ;
xr00 = x0[rloc0] ; xi00 = x0[iloc0] ;
xr01 = x1[rloc0] ; xi01 = x1[iloc0] ;
xr02 = x2[rloc0] ; xi02 = x2[iloc0] ;
xr10 = x0[rloc1] ; xi10 = x0[iloc1] ;
xr11 = x1[rloc1] ; xi11 = x1[iloc1] ;
xr12 = x2[rloc1] ; xi12 = x2[iloc1] ;
y0[rloc0] = ar00*xr00 - ai00*xi00 + ar01*xr10 - ai01*xi10 ;
y0[iloc0] = ar00*xi00 + ai00*xr00 + ar01*xi10 + ai01*xr10 ;
y1[rloc0] = ar00*xr01 - ai00*xi01 + ar01*xr11 - ai01*xi11 ;
y1[iloc0] = ar00*xi01 + ai00*xr01 + ar01*xi11 + ai01*xr11 ;
y2[rloc0] = ar00*xr02 - ai00*xi02 + ar01*xr12 - ai01*xi12 ;
y2[iloc0] = ar00*xi02 + ai00*xr02 + ar01*xi12 + ai01*xr12 ;
y0[rloc1] = ar01*xr00 - ai01*xi00 + ar11*xr10 - ai11*xi10 ;
y0[iloc1] = ar01*xi00 + ai01*xr00 + ar11*xi10 + ai11*xr10 ;
y1[rloc1] = ar01*xr01 - ai01*xi01 + ar11*xr11 - ai11*xi11 ;
y1[iloc1] = ar01*xi01 + ai01*xr01 + ar11*xi11 + ai11*xr11 ;
y2[rloc1] = ar01*xr02 - ai01*xi02 + ar11*xr12 - ai11*xi12 ;
y2[iloc1] = ar01*xi02 + ai01*xr02 + ar11*xi12 + ai11*xr12 ;
kk += 6, irowA += 2 ; rloc += 4 ; iloc += 4 ;
} else {
fprintf(stderr, "\n fatal error in SubMtx_scale3vec()"
"\n pivotsizes[%d] = %d", ipivot, pivotsizes[ipivot]) ;
exit(-1) ;
}
}
}
return ; }
/*--------------------------------------------------------------------*/
/*
--------------------------------------------------
purpose -- compute [y0 y1] := A * [x0 x1]
where A is block diagonal and symmetric
created -- 98apr17, cca
--------------------------------------------------
*/
static void
block_diagonal_sym_scale2vec (
SubMtx *mtxA,
double y0[],
double y1[],
double x0[],
double x1[]
) {
double *entries ;
int nentA, nrowA ;
int *pivotsizes ;
SubMtx_blockDiagonalInfo(mtxA, &nrowA, &nentA, &pivotsizes, &entries) ;
if ( SUBMTX_IS_REAL(mtxA) ) {
int ipivot, irowA, kk ;
for ( irowA = ipivot = kk = 0 ; irowA < nrowA ; ipivot++ ) {
if ( pivotsizes[ipivot] == 1 ) {
double a00, x00, x01 ;
a00 = entries[kk] ;
x00 = x0[irowA] ;
x01 = x1[irowA] ;
y0[irowA] = a00*x00 ;
y1[irowA] = a00*x01 ;
kk++, irowA++ ;
} else if ( pivotsizes[ipivot] == 2 ) {
double a00, a01, a11, x00, x01, x10, x11 ;
a00 = entries[kk] ;
a01 = entries[kk+1] ;
a11 = entries[kk+2] ;
x00 = x0[irowA] ;
x01 = x1[irowA] ;
x10 = x0[irowA+1] ;
x11 = x1[irowA+1] ;
y0[irowA] = a00*x00 + a01*x10 ;
y1[irowA] = a00*x01 + a01*x11 ;
y0[irowA+1] = a01*x00 + a11*x10 ;
y1[irowA+1] = a01*x01 + a11*x11 ;
kk += 3, irowA += 2 ;
} else {
fprintf(stderr, "\n fatal error in SubMtx_scale3vec()"
"\n pivotsizes[%d] = %d", ipivot, pivotsizes[ipivot]) ;
exit(-1) ;
}
}
} else if ( SUBMTX_IS_COMPLEX(mtxA) ) {
int iloc, ipivot, irowA, kk, rloc ;
for ( irowA = ipivot = kk = rloc = 0, iloc = 1 ;
irowA < nrowA ;
ipivot++ ) {
if ( pivotsizes[ipivot] == 1 ) {
double ai00, ar00, xi0, xi1, xr0, xr1 ;
ar00 = entries[kk] ; ai00 = entries[kk+1] ;
xr0 = x0[rloc] ; xi0 = x0[iloc] ;
xr1 = x1[rloc] ; xi1 = x1[iloc] ;
y0[rloc] = ar00*xr0 - ai00*xi0; y0[iloc] = ar00*xi0 + ai00*xr0;
y1[rloc] = ar00*xr1 - ai00*xi1; y1[iloc] = ar00*xi1 + ai00*xr1;
kk += 2, irowA++, rloc += 2, iloc += 2 ;
} else if ( pivotsizes[ipivot] == 2 ) {
double ai00, ai01, ai11, ar00, ar01, ar11,
xi00, xi01, xi10, xi11, xr00, xr01, xr10, xr11 ;
int iloc0, iloc1, rloc0, rloc1 ;
ar00 = entries[kk] ; ai00 = entries[kk+1] ;
ar01 = entries[kk+2] ; ai01 = entries[kk+3] ;
ar11 = entries[kk+4] ; ai11 = entries[kk+5] ;
rloc0 = rloc ; iloc0 = iloc ;
rloc1 = rloc + 2 ; iloc1 = iloc + 2 ;
xr00 = x0[rloc0] ; xi00 = x0[iloc0] ;
xr01 = x1[rloc0] ; xi01 = x1[iloc0] ;
xr10 = x0[rloc1] ; xi10 = x0[iloc1] ;
xr11 = x1[rloc1] ; xi11 = x1[iloc1] ;
y0[rloc0] = ar00*xr00 - ai00*xi00 + ar01*xr10 - ai01*xi10 ;
y0[iloc0] = ar00*xi00 + ai00*xr00 + ar01*xi10 + ai01*xr10 ;
y1[rloc0] = ar00*xr01 - ai00*xi01 + ar01*xr11 - ai01*xi11 ;
y1[iloc0] = ar00*xi01 + ai00*xr01 + ar01*xi11 + ai01*xr11 ;
y0[rloc1] = ar01*xr00 - ai01*xi00 + ar11*xr10 - ai11*xi10 ;
y0[iloc1] = ar01*xi00 + ai01*xr00 + ar11*xi10 + ai11*xr10 ;
y1[rloc1] = ar01*xr01 - ai01*xi01 + ar11*xr11 - ai11*xi11 ;
y1[iloc1] = ar01*xi01 + ai01*xr01 + ar11*xi11 + ai11*xr11 ;
kk += 6, irowA += 2 ; rloc += 4 ; iloc += 4 ;
} else {
fprintf(stderr, "\n fatal error in SubMtx_scale2vec()"
"\n pivotsizes[%d] = %d", ipivot, pivotsizes[ipivot]) ;
exit(-1) ;
}
}
}
return ; }
/*--------------------------------------------------------------------*/
/*
--------------------------------------------------
purpose -- compute [y0] := A * [x0]
where A is block diagonal and symmetric
created -- 98apr17, cca
--------------------------------------------------
*/
static void
block_diagonal_sym_scale1vec (
SubMtx *mtxA,
double y0[],
double x0[]
) {
double *entries ;
int nentA, nrowA ;
int *pivotsizes ;
SubMtx_blockDiagonalInfo(mtxA, &nrowA, &nentA, &pivotsizes, &entries) ;
if ( SUBMTX_IS_REAL(mtxA) ) {
int ipivot, irowA, kk ;
for ( irowA = ipivot = kk = 0 ; irowA < nrowA ; ipivot++ ) {
if ( pivotsizes[ipivot] == 1 ) {
double a00, x00 ;
a00 = entries[kk] ;
x00 = x0[irowA] ;
y0[irowA] = a00*x00 ;
kk++, irowA++ ;
} else if ( pivotsizes[ipivot] == 2 ) {
double a00, a01, a11, x00, x10 ;
a00 = entries[kk] ;
a01 = entries[kk+1] ;
a11 = entries[kk+2] ;
x00 = x0[irowA] ;
x10 = x0[irowA+1] ;
y0[irowA] = a00*x00 + a01*x10 ;
y0[irowA+1] = a01*x00 + a11*x10 ;
kk += 3, irowA += 2 ;
} else {
fprintf(stderr, "\n fatal error in SubMtx_scale3vec()"
"\n pivotsizes[%d] = %d", ipivot, pivotsizes[ipivot]) ;
exit(-1) ;
}
}
} else if ( SUBMTX_IS_COMPLEX(mtxA) ) {
int iloc, ipivot, irowA, kk, rloc ;
for ( irowA = ipivot = kk = rloc = 0, iloc = 1 ;
irowA < nrowA ;
ipivot++ ) {
if ( pivotsizes[ipivot] == 1 ) {
double ai00, ar00, xi0, xr0 ;
ar00 = entries[kk] ; ai00 = entries[kk+1] ;
xr0 = x0[rloc] ; xi0 = x0[iloc] ;
y0[rloc] = ar00*xr0 - ai00*xi0; y0[iloc] = ar00*xi0 + ai00*xr0;
kk += 2, irowA++, rloc += 2, iloc += 2 ;
} else if ( pivotsizes[ipivot] == 2 ) {
double ai00, ai01, ai11, ar00, ar01, ar11,
xi00, xi10, xr00, xr10 ;
int iloc0, iloc1, rloc0, rloc1 ;
ar00 = entries[kk] ; ai00 = entries[kk+1] ;
ar01 = entries[kk+2] ; ai01 = entries[kk+3] ;
ar11 = entries[kk+4] ; ai11 = entries[kk+5] ;
rloc0 = rloc ; iloc0 = iloc ;
rloc1 = rloc + 2 ; iloc1 = iloc + 2 ;
xr00 = x0[rloc0] ; xi00 = x0[iloc0] ;
xr10 = x0[rloc1] ; xi10 = x0[iloc1] ;
y0[rloc0] = ar00*xr00 - ai00*xi00 + ar01*xr10 - ai01*xi10 ;
y0[iloc0] = ar00*xi00 + ai00*xr00 + ar01*xi10 + ai01*xr10 ;
y0[rloc1] = ar01*xr00 - ai01*xi00 + ar11*xr10 - ai11*xi10 ;
y0[iloc1] = ar01*xi00 + ai01*xr00 + ar11*xi10 + ai11*xr10 ;
kk += 6, irowA += 2 ; rloc += 4 ; iloc += 4 ;
} else {
fprintf(stderr, "\n fatal error in SubMtx_scale1vec()"
"\n pivotsizes[%d] = %d", ipivot, pivotsizes[ipivot]) ;
exit(-1) ;
}
}
}
return ; }
/*--------------------------------------------------------------------*/
/*
--------------------------------------------------
purpose -- compute [y0 y1 y2] := A * [x0 x1 x2]
where A is block diagonal and hermitian
created -- 98apr17, cca
--------------------------------------------------
*/
static void
block_diagonal_herm_scale3vec (
SubMtx *mtxA,
double y0[],
double y1[],
double y2[],
double x0[],
double x1[],
double x2[]
) {
double *entries ;
int iloc, ipivot, irowA, kk, nentA, nrowA, rloc ;
int *pivotsizes ;
SubMtx_blockDiagonalInfo(mtxA, &nrowA, &nentA, &pivotsizes, &entries) ;
for ( irowA = ipivot = kk = rloc = 0, iloc = 1 ;
irowA < nrowA ;
ipivot++ ) {
if ( pivotsizes[ipivot] == 1 ) {
double ar00, ai00, xi0, xi1, xi2, xr0, xr1, xr2 ;
/*
ar00 = entries[kk++] ; ai00 = entries[kk++] ;
*/
ar00 = entries[kk++] ; ai00 = 0.0 ; kk++ ;
xr0 = x0[rloc] ; xi0 = x0[iloc] ;
xr1 = x1[rloc] ; xi1 = x1[iloc] ;
xr2 = x2[rloc] ; xi2 = x2[iloc] ;
y0[rloc] = ar00*xr0 ; y0[iloc] = ar00*xi0 ;
y1[rloc] = ar00*xr1 ; y1[iloc] = ar00*xi1 ;
y2[rloc] = ar00*xr2 ; y2[iloc] = ar00*xi2 ;
irowA++, rloc += 2, iloc += 2 ;
} else if ( pivotsizes[ipivot] == 2 ) {
double ai00, ai01, ai11, ar00, ar01, ar11,
xi00, xi01, xi02, xi10, xi11, xi12,
xr00, xr01, xr02, xr10, xr11, xr12 ;
int iloc0, iloc1, rloc0, rloc1 ;
rloc0 = rloc ; iloc0 = iloc ;
rloc1 = rloc + 2 ; iloc1 = iloc + 2 ;
/*
ar00 = entries[kk++] ; ai00 = entries[kk++] ;
ar01 = entries[kk++] ; ai01 = entries[kk++] ;
ar11 = entries[kk++] ; ai11 = entries[kk++] ;
*/
ar00 = entries[kk++] ; ai00 = 0.0 ; kk++ ;
ar01 = entries[kk++] ; ai01 = entries[kk++] ;
ar11 = entries[kk++] ; ai11 = 0.0 ; kk++ ;
xr00 = x0[rloc0] ; xi00 = x0[iloc0] ;
xr01 = x1[rloc0] ; xi01 = x1[iloc0] ;
xr02 = x2[rloc0] ; xi02 = x2[iloc0] ;
xr10 = x0[rloc1] ; xi10 = x0[iloc1] ;
xr11 = x1[rloc1] ; xi11 = x1[iloc1] ;
xr12 = x2[rloc1] ; xi12 = x2[iloc1] ;
y0[rloc0] = ar00*xr00 + ar01*xr10 - ai01*xi10 ;
y0[iloc0] = ar00*xi00 + ar01*xi10 + ai01*xr10 ;
y1[rloc0] = ar00*xr01 + ar01*xr11 - ai01*xi11 ;
y1[iloc0] = ar00*xi01 + ar01*xi11 + ai01*xr11 ;
y2[rloc0] = ar00*xr02 + ar01*xr12 - ai01*xi12 ;
y2[iloc0] = ar00*xi02 + ar01*xi12 + ai01*xr12 ;
y0[rloc1] = ar01*xr00 + ai01*xi00 + ar11*xr10 ;
y0[iloc1] = ar01*xi00 - ai01*xr00 + ar11*xi10 ;
y1[rloc1] = ar01*xr01 + ai01*xi01 + ar11*xr11 ;
y1[iloc1] = ar01*xi01 - ai01*xr01 + ar11*xi11 ;
y2[rloc1] = ar01*xr02 + ai01*xi02 + ar11*xr12 ;
y2[iloc1] = ar01*xi02 - ai01*xr02 + ar11*xi12 ;
irowA += 2 ; rloc += 4 ; iloc += 4 ;
} else {
fprintf(stderr, "\n fatal error in SubMtx_scale3vec()"
"\n pivotsizes[%d] = %d", ipivot, pivotsizes[ipivot]) ;
exit(-1) ;
}
}
return ; }
/*--------------------------------------------------------------------*/
/*
--------------------------------------------------
purpose -- compute [y0 y1] := A * [x0 x1]
where A is block diagonal and hermitian
created -- 98apr17, cca
--------------------------------------------------
*/
static void
block_diagonal_herm_scale2vec (
SubMtx *mtxA,
double y0[],
double y1[],
double x0[],
double x1[]
) {
double *entries ;
int iloc, ipivot, irowA, kk, nentA, nrowA, rloc ;
int *pivotsizes ;
SubMtx_blockDiagonalInfo(mtxA, &nrowA, &nentA, &pivotsizes, &entries) ;
for ( irowA = ipivot = kk = rloc = 0, iloc = 1 ;
irowA < nrowA ;
ipivot++ ) {
if ( pivotsizes[ipivot] == 1 ) {
double ai00, ar00, xi0, xi1, xr0, xr1 ;
/*
ar00 = entries[kk++] ; ai00 = entries[kk++] ;
*/
ar00 = entries[kk++] ; ai00 = 0.0 ; kk++ ;
xr0 = x0[rloc] ; xi0 = x0[iloc] ;
xr1 = x1[rloc] ; xi1 = x1[iloc] ;
y0[rloc] = ar00*xr0 - ai00*xi0 ; y0[iloc] = ar00*xi0 + ai00*xr0 ;
y1[rloc] = ar00*xr1 - ai00*xi1 ; y1[iloc] = ar00*xi1 + ai00*xr1 ;
irowA++, rloc += 2, iloc += 2 ;
} else if ( pivotsizes[ipivot] == 2 ) {
double ai00, ai01, ai11, ar00, ar01, ar11,
xi00, xi01, xi10, xi11, xr00, xr01, xr10, xr11 ;
int iloc0, iloc1, rloc0, rloc1 ;
rloc0 = rloc ; iloc0 = iloc ;
rloc1 = rloc + 2 ; iloc1 = iloc + 2 ;
/*
ar00 = entries[kk++] ; ai00 = entries[kk++] ;
ar01 = entries[kk++] ; ai01 = entries[kk++] ;
ar11 = entries[kk++] ; ai11 = entries[kk++] ;
*/
ar00 = entries[kk++] ; ai00 = 0.0 ; kk++ ;
ar01 = entries[kk++] ; ai01 = entries[kk++] ;
ar11 = entries[kk++] ; ai11 = 0.0 ; kk++ ;
xr00 = x0[rloc0] ; xi00 = x0[iloc0] ;
xr01 = x1[rloc0] ; xi01 = x1[iloc0] ;
xr10 = x0[rloc1] ; xi10 = x0[iloc1] ;
xr11 = x1[rloc1] ; xi11 = x1[iloc1] ;
y0[rloc0] = ar00*xr00 + ar01*xr10 - ai01*xi10 ;
y0[iloc0] = ar00*xi00 + ar01*xi10 + ai01*xr10 ;
y1[rloc0] = ar00*xr01 + ar01*xr11 - ai01*xi11 ;
y1[iloc0] = ar00*xi01 + ar01*xi11 + ai01*xr11 ;
y0[rloc1] = ar01*xr00 + ai01*xi00 + ar11*xr10 ;
y0[iloc1] = ar01*xi00 - ai01*xr00 + ar11*xi10 ;
y1[rloc1] = ar01*xr01 + ai01*xi01 + ar11*xr11 ;
y1[iloc1] = ar01*xi01 - ai01*xr01 + ar11*xi11 ;
irowA += 2 ; rloc += 4 ; iloc += 4 ;
} else {
fprintf(stderr, "\n fatal error in SubMtx_scale2vec()"
"\n pivotsizes[%d] = %d", ipivot, pivotsizes[ipivot]) ;
exit(-1) ;
}
}
return ; }
/*--------------------------------------------------------------------*/
/*
--------------------------------------------------
purpose -- compute [y0] := A * [x0]
where A is block diagonal and hermitian
created -- 98apr17, cca
--------------------------------------------------
*/
static void
block_diagonal_herm_scale1vec (
SubMtx *mtxA,
double y0[],
double x0[]
) {
double *entries ;
int iloc, ipivot, irowA, kk, nentA, nrowA, rloc ;
int *pivotsizes ;
SubMtx_blockDiagonalInfo(mtxA, &nrowA, &nentA, &pivotsizes, &entries) ;
for ( irowA = ipivot = kk = rloc = 0, iloc = 1 ;
irowA < nrowA ;
ipivot++ ) {
if ( pivotsizes[ipivot] == 1 ) {
double ai00, ar00, xi0, xr0 ;
/*
ar00 = entries[kk++] ; ai00 = entries[kk++] ;
*/
ar00 = entries[kk++] ; ai00 = 0.0 ; kk++ ;
xr0 = x0[rloc] ; xi0 = x0[iloc] ;
y0[rloc] = ar00*xr0 - ai00*xi0 ; y0[iloc] = ar00*xi0 + ai00*xr0 ;
irowA++, rloc += 2, iloc += 2 ;
} else if ( pivotsizes[ipivot] == 2 ) {
double ai00, ai01, ai11, ar00, ar01, ar11,
xi00, xi10, xr00, xr10 ;
int iloc0, iloc1, rloc0, rloc1 ;
rloc0 = rloc ; iloc0 = iloc ;
rloc1 = rloc + 2 ; iloc1 = iloc + 2 ;
/*
ar00 = entries[kk++] ; ai00 = entries[kk++] ;
ar01 = entries[kk++] ; ai01 = entries[kk++] ;
ar11 = entries[kk++] ; ai11 = entries[kk++] ;
*/
ar00 = entries[kk++] ; ai00 = 0.0 ; kk++ ;
ar01 = entries[kk++] ; ai01 = entries[kk++] ;
ar11 = entries[kk++] ; ai11 = 0.0 ; kk++ ;
xr00 = x0[rloc0] ; xi00 = x0[iloc0] ;
xr10 = x0[rloc1] ; xi10 = x0[iloc1] ;
y0[rloc0] = ar00*xr00 + ar01*xr10 - ai01*xi10 ;
y0[iloc0] = ar00*xi00 + ar01*xi10 + ai01*xr10 ;
y0[rloc1] = ar01*xr00 + ai01*xi00 + ar11*xr10 ;
y0[iloc1] = ar01*xi00 - ai01*xr00 + ar11*xi10 ;
irowA += 2 ; rloc += 4 ; iloc += 4 ;
} else {
fprintf(stderr, "\n fatal error in SubMtx_scale1vec()"
"\n pivotsizes[%d] = %d", ipivot, pivotsizes[ipivot]) ;
exit(-1) ;
}
}
return ; }
/*--------------------------------------------------------------------*/
syntax highlighted by Code2HTML, v. 0.9.1