| Line 34... |
Line 34... |
#include <slepcblaslapack.h>
|
#include <slepcblaslapack.h>
|
|
|
PetscErrorCode dvd_calcpairs_proj(dvdDashboard *d);
|
PetscErrorCode dvd_calcpairs_proj(dvdDashboard *d);
|
PetscErrorCode dvd_calcpairs_qz_start(dvdDashboard *d);
|
PetscErrorCode dvd_calcpairs_qz_start(dvdDashboard *d);
|
PetscErrorCode dvd_calcpairs_projeig_eig(dvdDashboard *d);
|
PetscErrorCode dvd_calcpairs_projeig_eig(dvdDashboard *d);
|
|
PetscErrorCode dvd_calcpairs_projeig_ind(dvdDashboard *d);
|
PetscErrorCode dvd_calcpairs_projeig_qz_std(dvdDashboard *d);
|
PetscErrorCode dvd_calcpairs_projeig_qz_std(dvdDashboard *d);
|
PetscErrorCode dvd_calcpairs_projeig_qz_gen(dvdDashboard *d);
|
PetscErrorCode dvd_calcpairs_projeig_qz_gen(dvdDashboard *d);
|
PetscErrorCode dvd_calcpairs_selectPairs_qz(dvdDashboard *d, PetscInt n);
|
PetscErrorCode dvd_calcpairs_selectPairs_qz(dvdDashboard *d, PetscInt n);
|
PetscErrorCode dvd_calcpairs_selectPairs_eig(dvdDashboard *d, PetscInt n);
|
PetscErrorCode dvd_calcpairs_selectPairs_eig(dvdDashboard *d, PetscInt n);
|
PetscErrorCode dvd_calcpairs_X(dvdDashboard *d, PetscInt r_s, PetscInt r_e,
|
PetscErrorCode dvd_calcpairs_X(dvdDashboard *d, PetscInt r_s, PetscInt r_e,
|
| Line 70... |
Line 71... |
orthoV_type_t orth, IP ipI,
|
orthoV_type_t orth, IP ipI,
|
PetscInt cX_proj, PetscBool harm)
|
PetscInt cX_proj, PetscBool harm)
|
{
|
{
|
PetscErrorCode ierr;
|
PetscErrorCode ierr;
|
PetscInt i,max_cS;
|
PetscInt i,max_cS;
|
PetscBool std_probl,her_probl;
|
PetscBool std_probl,her_probl,ind_probl;
|
|
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
|
|
std_probl = DVD_IS(d->sEP, DVD_EP_STD)?PETSC_TRUE:PETSC_FALSE;
|
std_probl = DVD_IS(d->sEP, DVD_EP_STD)?PETSC_TRUE:PETSC_FALSE;
|
her_probl = DVD_IS(d->sEP, DVD_EP_HERMITIAN)?PETSC_TRUE:PETSC_FALSE;
|
her_probl = DVD_IS(d->sEP, DVD_EP_HERMITIAN)?PETSC_TRUE:PETSC_FALSE;
|
|
ind_probl = DVD_IS(d->sEP, DVD_EP_INDEFINITE)?PETSC_TRUE:PETSC_FALSE;
|
|
|
/* Setting configuration constrains */
|
/* Setting configuration constrains */
|
#if !defined(PETSC_USE_COMPLEX)
|
#if !defined(PETSC_USE_COMPLEX)
|
/* if the last converged eigenvalue is complex its conjugate pair is also
|
/* if the last converged eigenvalue is complex its conjugate pair is also
|
converged */
|
converged */
|
| Line 91... |
Line 93... |
d->size_real_V = b->max_size_V+b->max_nev;
|
d->size_real_V = b->max_size_V+b->max_nev;
|
d->W_shift = d->B?PETSC_TRUE:PETSC_FALSE;
|
d->W_shift = d->B?PETSC_TRUE:PETSC_FALSE;
|
d->size_real_W = harm?(b->max_size_V+(d->W_shift?b->max_nev:b->max_size_cP)):0;
|
d->size_real_W = harm?(b->max_size_V+(d->W_shift?b->max_nev:b->max_size_cP)):0;
|
d->size_real_AV = b->max_size_V+b->max_size_cP;
|
d->size_real_AV = b->max_size_V+b->max_size_cP;
|
d->size_BDS = 0;
|
d->size_BDS = 0;
|
if (d->B && her_probl && (orth == DVD_ORTHOV_I || orth == DVD_ORTHOV_BOneMV)) {
|
if (d->B && (her_probl || ind_probl) && (orth == DVD_ORTHOV_I || orth == DVD_ORTHOV_BOneMV)) {
|
d->size_real_BV = b->size_V; d->BV_shift = PETSC_TRUE;
|
d->size_real_BV = b->size_V; d->BV_shift = PETSC_TRUE;
|
if (orth == DVD_ORTHOV_BOneMV) d->size_BDS = d->eps->nds;
|
if (orth == DVD_ORTHOV_BOneMV) d->size_BDS = d->eps->nds;
|
} else if (d->B) {
|
} else if (d->B) {
|
d->size_real_BV = b->max_size_V + b->max_size_P; d->BV_shift = PETSC_FALSE;
|
d->size_real_BV = b->max_size_V + b->max_size_P; d->BV_shift = PETSC_FALSE;
|
} else {
|
} else {
|
| Line 105... |
Line 107... |
d->size_real_BV + d->size_BDS;
|
d->size_real_BV + d->size_BDS;
|
b->own_scalars+= b->max_size_proj*b->max_size_proj*2*(std_probl?1:2) +
|
b->own_scalars+= b->max_size_proj*b->max_size_proj*2*(std_probl?1:2) +
|
/* H, G?, S, T? */
|
/* H, G?, S, T? */
|
b->max_size_proj*b->max_size_proj*(std_probl?1:2) +
|
b->max_size_proj*b->max_size_proj*(std_probl?1:2) +
|
/* pX, pY? */
|
/* pX, pY? */
|
b->max_nev*b->max_nev*(her_probl?0:(!d->B?1:2));
|
b->max_nev*b->max_nev*(her_probl?0:(!d->B?1:2)) +
|
/* cS?, cT? */
|
/* cS?, cT? */
|
|
FromRealToScalar(d->size_real_V)*(ind_probl?1:0) + /* nBV */
|
|
FromRealToScalar(b->max_size_proj)*(ind_probl?1:0); /* nBpX */
|
b->max_size_auxV = PetscMax(b->max_size_auxV, b->max_size_X);
|
b->max_size_auxV = PetscMax(b->max_size_auxV, b->max_size_X);
|
/* updateV0 */
|
/* updateV0 */
|
max_cS = PetscMax(b->max_size_X,cX_proj);
|
max_cS = PetscMax(b->max_size_X,cX_proj);
|
b->max_size_auxS = PetscMax(PetscMax(
|
b->max_size_auxS = PetscMax(PetscMax(
|
b->max_size_auxS,
|
b->max_size_auxS,
|
| Line 127... |
Line 131... |
std_probl?0:(b->max_size_proj*11+16) /* projeig */);
|
std_probl?0:(b->max_size_proj*11+16) /* projeig */);
|
#if defined(PETSC_USE_COMPLEX)
|
#if defined(PETSC_USE_COMPLEX)
|
b->max_size_auxS = PetscMax(b->max_size_auxS, b->max_size_V);
|
b->max_size_auxS = PetscMax(b->max_size_auxS, b->max_size_V);
|
/* dvd_calcpairs_projeig_eig */
|
/* dvd_calcpairs_projeig_eig */
|
#endif
|
#endif
|
|
if (ind_probl)
|
|
b->max_size_auxS = PetscMax(b->max_size_auxS,2*b->max_size_proj);
|
|
/* SlepcDenseIOrth */
|
|
|
/* Setup the step */
|
/* Setup the step */
|
if (b->state >= DVD_STATE_CONF) {
|
if (b->state >= DVD_STATE_CONF) {
|
d->max_cX_in_proj = cX_proj;
|
d->max_cX_in_proj = cX_proj;
|
d->max_size_P = b->max_size_P;
|
d->max_size_P = b->max_size_P;
|
| Line 150... |
Line 157... |
d->cS = b->free_scalars; b->free_scalars+= b->max_nev*b->max_nev;
|
d->cS = b->free_scalars; b->free_scalars+= b->max_nev*b->max_nev;
|
d->max_size_cS = d->ldcS = b->max_nev;
|
d->max_size_cS = d->ldcS = b->max_nev;
|
} else {
|
} else {
|
d->cS = PETSC_NULL;
|
d->cS = PETSC_NULL;
|
d->max_size_cS = d->ldcS = 0;
|
d->max_size_cS = d->ldcS = 0;
|
|
d->orthoV_type = orth;
|
|
if (ind_probl) {
|
|
d->real_nBV = b->free_scalars; b->free_scalars+= FromRealToScalar(d->size_real_V);
|
|
d->nBpX = b->free_scalars; b->free_scalars+= FromRealToScalar(d->max_size_proj);
|
|
} else
|
|
d->real_nBV = d->nBDS = d->nBpX = PETSC_NULL;
|
}
|
}
|
d->ipV = ipI;
|
d->ipV = ipI;
|
d->ipW = ipI;
|
d->ipW = ipI;
|
if (orth == DVD_ORTHOV_BOneMV) {
|
if (orth == DVD_ORTHOV_BOneMV) {
|
d->BDS = b->free_vecs; b->free_vecs+= d->eps->nds;
|
d->BDS = b->free_vecs; b->free_vecs+= d->eps->nds;
|
| Line 201... |
Line 214... |
|
|
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "dvd_calcpairs_qz_start"
|
#define __FUNCT__ "dvd_calcpairs_qz_start"
|
PetscErrorCode dvd_calcpairs_qz_start(dvdDashboard *d)
|
PetscErrorCode dvd_calcpairs_qz_start(dvdDashboard *d)
|
{
|
{
|
PetscBool her_probl;
|
PetscBool her_probl,ind_probl;
|
PetscInt i;
|
PetscInt i;
|
|
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
her_probl = DVD_IS(d->sEP, DVD_EP_HERMITIAN)?PETSC_TRUE:PETSC_FALSE;
|
her_probl = DVD_IS(d->sEP, DVD_EP_HERMITIAN)?PETSC_TRUE:PETSC_FALSE;
|
|
ind_probl = DVD_IS(d->sEP, DVD_EP_INDEFINITE)?PETSC_TRUE:PETSC_FALSE;
|
|
|
d->size_V = 0;
|
d->size_V = 0;
|
d->V = d->real_V;
|
d->V = d->real_V;
|
d->cX = d->real_V;
|
d->cX = d->real_V;
|
d->size_cX = 0;
|
d->size_cX = 0;
|
| Line 232... |
Line 246... |
d->cY = d->B && !her_probl ? d->W : PETSC_NULL;
|
d->cY = d->B && !her_probl ? d->W : PETSC_NULL;
|
d->BcX = d->orthoV_type == DVD_ORTHOV_I && d->B && her_probl ? d->BcX : PETSC_NULL;
|
d->BcX = d->orthoV_type == DVD_ORTHOV_I && d->B && her_probl ? d->BcX : PETSC_NULL;
|
d->size_cY = 0;
|
d->size_cY = 0;
|
d->size_BcX = 0;
|
d->size_BcX = 0;
|
d->cX_in_V = d->cX_in_H = d->cX_in_G = d->cX_in_W = d->cX_in_AV = d->cX_in_BV = 0;
|
d->cX_in_V = d->cX_in_H = d->cX_in_G = d->cX_in_W = d->cX_in_AV = d->cX_in_BV = 0;
|
|
d->nBV = d->nBcX = d->real_nBV;
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
}
|
}
|
|
|
|
|
#undef __FUNCT__
|
#undef __FUNCT__
|
| Line 317... |
Line 332... |
d->V_tra_s = d->V_tra_e = 0;
|
d->V_tra_s = d->V_tra_e = 0;
|
d->V_new_s = d->V_new_e;
|
d->V_new_s = d->V_new_e;
|
|
|
/* Solve the projected problem */
|
/* Solve the projected problem */
|
if (DVD_IS(d->sEP, DVD_EP_STD)) {
|
if (DVD_IS(d->sEP, DVD_EP_STD)) {
|
if (DVD_IS(d->sEP, DVD_EP_HERMITIAN)) {
|
if (DVD_IS(d->sEP, DVD_EP_INDEFINITE)) {
|
|
ierr = dvd_calcpairs_projeig_ind(d); CHKERRQ(ierr);
|
|
} else if (DVD_IS(d->sEP, DVD_EP_HERMITIAN)) {
|
ierr = dvd_calcpairs_projeig_eig(d); CHKERRQ(ierr);
|
ierr = dvd_calcpairs_projeig_eig(d); CHKERRQ(ierr);
|
} else {
|
} else {
|
ierr = dvd_calcpairs_projeig_qz_std(d); CHKERRQ(ierr);
|
ierr = dvd_calcpairs_projeig_qz_std(d); CHKERRQ(ierr);
|
}
|
}
|
} else {
|
} else {
|
| Line 350... |
Line 367... |
/* auxV: V_tra_s, DvdMult_copy_func: 1 */
|
/* auxV: V_tra_s, DvdMult_copy_func: 1 */
|
PetscErrorCode dvd_calcpairs_updateV0(dvdDashboard *d, DvdReduction *r,
|
PetscErrorCode dvd_calcpairs_updateV0(dvdDashboard *d, DvdReduction *r,
|
DvdMult_copy_func **sr)
|
DvdMult_copy_func **sr)
|
{
|
{
|
PetscErrorCode ierr;
|
PetscErrorCode ierr;
|
PetscInt rm;
|
PetscInt rm,i;
|
|
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
|
|
if (d->V_tra_s == 0 && d->V_tra_e == 0) { PetscFunctionReturn(0); }
|
if (d->V_tra_s == 0 && d->V_tra_e == 0) { PetscFunctionReturn(0); }
|
|
|
|
/* Update nBcX and nBV */
|
|
if (d->nBcX && d->nBpX && d->nBV) {
|
|
d->nBV+= d->V_tra_s;
|
|
for (i=0; i<d->V_tra_s; i++) d->nBcX[d->size_cX+i] = d->nBpX[i];
|
|
for (i=d->V_tra_s; i<d->V_tra_e; i++) d->nBV[i-d->V_tra_s] = d->nBpX[i];
|
|
}
|
|
|
/* cX <- [cX V*MT(0:V_tra_s-1)], V <- V*MT(V_tra_s:V_tra_e) */
|
/* cX <- [cX V*MT(0:V_tra_s-1)], V <- V*MT(V_tra_s:V_tra_e) */
|
ierr = dvd_calcpairs_updateBV0_gen(d,d->real_V,&d->size_cX,&d->V,&d->size_V,&d->max_size_V,PETSC_TRUE,&d->cX_in_V,d->MTX);CHKERRQ(ierr);
|
ierr = dvd_calcpairs_updateBV0_gen(d,d->real_V,&d->size_cX,&d->V,&d->size_V,&d->max_size_V,PETSC_TRUE,&d->cX_in_V,d->MTX);CHKERRQ(ierr);
|
|
|
/* Udpate cS for standard problems */
|
/* Udpate cS for standard problems */
|
| Line 372... |
Line 396... |
|
|
/* cS(:, size_cS:) <- cX' * auxV */
|
/* cS(:, size_cS:) <- cX' * auxV */
|
ierr = VecsMultS(&d->cS[d->ldcS*d->size_cS], 0, d->ldcS, d->cX, 0, d->size_cX-rm, d->auxV, 0, d->V_tra_s-rm, r, (*sr)++); CHKERRQ(ierr);
|
ierr = VecsMultS(&d->cS[d->ldcS*d->size_cS], 0, d->ldcS, d->cX, 0, d->size_cX-rm, d->auxV, 0, d->V_tra_s-rm, r, (*sr)++); CHKERRQ(ierr);
|
d->size_cS+= d->V_tra_s-rm;
|
d->size_cS+= d->V_tra_s-rm;
|
}
|
}
|
|
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
}
|
}
|
|
|
|
|
#undef __FUNCT__
|
#undef __FUNCT__
|
| Line 394... |
Line 418... |
/* Check consistency */
|
/* Check consistency */
|
if (d->size_V != d->V_new_s) { SETERRQ(PETSC_COMM_SELF,1, "Consistency broken!"); }
|
if (d->size_V != d->V_new_s) { SETERRQ(PETSC_COMM_SELF,1, "Consistency broken!"); }
|
|
|
/* V <- gs([cX V(0:V_new_s-1)], V(V_new_s:V_new_e-1)) */
|
/* V <- gs([cX V(0:V_new_s-1)], V(V_new_s:V_new_e-1)) */
|
if (d->orthoV_type == DVD_ORTHOV_BOneMV) {
|
if (d->orthoV_type == DVD_ORTHOV_BOneMV) {
|
ierr = dvd_BorthV(d->ipV, d->eps->DS, d->BDS, d->eps->nds, d->cX, d->real_BV,
|
ierr = dvd_BorthV(d->ipV, d->eps->DS, d->BDS, d->nBDS, d->eps->nds, d->cX, d->real_BV, d->nBcX,
|
d->size_cX, d->V, d->BV, d->V_new_s, d->V_new_e,
|
d->size_cX, d->V, d->BV, d->nBV, d->V_new_s, d->V_new_e,
|
d->auxS, d->eps->rand); CHKERRQ(ierr);
|
d->auxS, d->eps->rand); CHKERRQ(ierr);
|
d->size_BV = d->V_new_e;
|
d->size_BV = d->V_new_e;
|
} else {
|
} else {
|
ierr = dvd_orthV(d->ipV, d->eps->DS, d->eps->nds, cX, d->size_cX, d->V,
|
ierr = dvd_orthV(d->ipV, d->eps->DS, d->eps->nds, cX, d->size_cX, d->V,
|
d->V_new_s, d->V_new_e, d->auxS, d->eps->rand);
|
d->V_new_s, d->V_new_e, d->auxS, d->eps->rand);
|
| Line 630... |
Line 654... |
d->calcpairs_selectPairs = dvd_calcpairs_selectPairs_eig;
|
d->calcpairs_selectPairs = dvd_calcpairs_selectPairs_eig;
|
|
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
}
|
}
|
|
|
|
|
|
/* in complex, d->size_H real auxiliar values are needed */
|
|
#undef __FUNCT__
|
|
#define __FUNCT__ "dvd_calcpairs_projeig_ind"
|
|
PetscErrorCode dvd_calcpairs_projeig_ind(dvdDashboard *d)
|
|
{
|
|
PetscErrorCode ierr;
|
|
PetscInt i,j;
|
|
|
|
PetscFunctionBegin;
|
|
|
|
/* S <- H */
|
|
d->ldS = d->ldpX = d->size_H;
|
|
ierr = SlepcDenseCopyTriang(d->S, 0, d->size_H, d->H, d->sH, d->ldH,
|
|
d->size_H, d->size_H); CHKERRQ(ierr);
|
|
|
|
/* S <- nBV * S */
|
|
for (i=0; i<d->size_H; i++) {
|
|
if (d->nBV[i] < 0.0)
|
|
for (j=0; j<d->size_H; j++) d->S[i+j*d->ldS]*= d->nBV[i];
|
|
}
|
|
|
|
/* S = pX * L */
|
|
ierr = EPSDenseNHEP(d->size_H, d->S, d->eigr-d->cX_in_H, d->eigi-d->cX_in_H, d->pX, PETSC_NULL); CHKERRQ(ierr);
|
|
|
|
d->calcpairs_selectPairs = dvd_calcpairs_selectPairs_eig;
|
|
|
|
PetscFunctionReturn(0);
|
|
}
|
|
|
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "dvd_calcpairs_projeig_qz_std"
|
#define __FUNCT__ "dvd_calcpairs_projeig_qz_std"
|
PetscErrorCode dvd_calcpairs_projeig_qz_std(dvdDashboard *d)
|
PetscErrorCode dvd_calcpairs_projeig_qz_std(dvdDashboard *d)
|
{
|
{
|
| Line 731... |
Line 784... |
ierr = EPSSortDenseHEP(d->eps, d->size_H, 0, d->eigr-d->cX_in_H, d->pX, d->ldpX);
|
ierr = EPSSortDenseHEP(d->eps, d->size_H, 0, d->eigr-d->cX_in_H, d->pX, d->ldpX);
|
CHKERRQ(ierr);
|
CHKERRQ(ierr);
|
|
|
if (d->calcpairs_eigs_trans) {
|
if (d->calcpairs_eigs_trans) {
|
ierr = d->calcpairs_eigs_trans(d); CHKERRQ(ierr);
|
ierr = d->calcpairs_eigs_trans(d); CHKERRQ(ierr);
|
|
}
|
|
|
|
if (DVD_IS(d->sEP, DVD_EP_INDEFINITE)) {
|
|
ierr = SlepcDenseIOrth(d->pX,d->ldpX,d->size_H,d->size_H,d->nBV-d->cX_in_H,d->nBpX,d->auxS,d->size_auxS);CHKERRQ(ierr);
|
}
|
}
|
|
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
}
|
}
|
|
|
| Line 843... |
Line 900... |
|
|
/* Otherwise, nR[i] <- ||R[i]|| */
|
/* Otherwise, nR[i] <- ||R[i]|| */
|
else
|
else
|
cX = PETSC_NULL;
|
cX = PETSC_NULL;
|
|
|
if (cX) for (i=0; i<r_e-r_s; i++) {
|
if (cX) {
|
ierr = IPOrthogonalize(d->ipI, 0, PETSC_NULL, d->size_cX, PETSC_NULL,
|
Vec auxV;
|
cX, R[i], PETSC_NULL, &d->nR[r_s+i], &lindep);
|
ierr = VecDuplicate(d->auxV[0],&auxV);CHKERRQ(ierr);
|
CHKERRQ(ierr);
|
for (i=0; i<r_e-r_s; i++) {
|
if(lindep || (d->nR[r_s+i] < PETSC_MACHINE_EPSILON)) {
|
if (cX && d->orthoV_type == DVD_ORTHOV_BOneMV) {
|
ierr = PetscInfo2(d->eps,"The computed eigenvector residual %D is too low, %G!\n",r_s+i,d->nR[r_s+i]);CHKERRQ(ierr);
|
ierr = IPBOrthogonalize(d->ipV,d->eps->nds,d->eps->DS,d->BDS,d->nBDS,d->size_cX,PETSC_NULL,d->cX,d->real_BV,d->nBcX,R[i],auxV,PETSC_NULL,&d->nR[r_s+i],&lindep);CHKERRQ(ierr);
|
|
} else {
|
|
ierr = IPOrthogonalize(d->ipI, 0, PETSC_NULL, d->size_cX, PETSC_NULL, cX, R[i], PETSC_NULL, &d->nR[r_s+i], &lindep); CHKERRQ(ierr);
|
|
}
|
|
if(lindep || (PetscAbs(d->nR[r_s+i]) < PETSC_MACHINE_EPSILON)) {
|
|
ierr = PetscInfo2(d->eps,"The computed eigenvector residual %D is too low, %G!\n",r_s+i,d->nR[r_s+i]);CHKERRQ(ierr);
|
|
}
|
}
|
}
|
|
ierr = VecDestroy(&auxV);CHKERRQ(ierr);
|
} else {
|
}
|
|
if (!cX || (cX && d->orthoV_type == DVD_ORTHOV_BOneMV)) {
|
for(i=0; i<r_e-r_s; i++) {
|
for(i=0; i<r_e-r_s; i++) {
|
ierr = VecNormBegin(R[i],NORM_2,&d->nR[r_s+i]); CHKERRQ(ierr);
|
ierr = VecNormBegin(R[i],NORM_2,&d->nR[r_s+i]); CHKERRQ(ierr);
|
}
|
}
|
for(i=0; i<r_e-r_s; i++) {
|
for(i=0; i<r_e-r_s; i++) {
|
ierr = VecNormEnd(R[i],NORM_2,&d->nR[r_s+i]); CHKERRQ(ierr);
|
ierr = VecNormEnd(R[i],NORM_2,&d->nR[r_s+i]); CHKERRQ(ierr);
|