Subversion Repositories slepc-dev

Compare Revisions

Ignore whitespace Rev 2588 → Rev 2589

/branches/slepc-3_2-proj/src/eps/impls/davidson/common/dvd_blas.c
314,6 → 314,11
SETERRQ(PETSC_COMM_SELF,1, "Leading dimension error");
}
 
if (sY == 0 && sX == 0) {
ierr = SlepcDenseCopy(Y, ldY, X, ldX, rX, cX);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
 
if (rX != cX) {
SETERRQ(PETSC_COMM_SELF,1, "Rectangular matrices not supported");
}
594,9 → 599,12
/* Full M matrix */
} else if (DVD_ISNOT(sM,DVD_MAT_UTRIANG) &&
DVD_ISNOT(sM,DVD_MAT_LTRIANG)) {
if (workS1)
if (workS1) {
Wr = workS1;
else {
if (PetscMin(W - workS1, workS1 - W) < ((eU-sU)*sV+(eV-sV)*eU)) {
SETERRQ(PETSC_COMM_SELF,1, "Consistency broken!");
}
} else {
ierr = PetscMalloc(sizeof(PetscScalar)*((eU-sU)*sV+(eV-sV)*eU), &Wr);
CHKERRQ(ierr);
}
630,9 → 638,12
/* Upper triangular M matrix */
} else if (DVD_IS(sM,DVD_MAT_UTRIANG) &&
DVD_ISNOT(sM,DVD_MAT_LTRIANG)) {
if (workS1)
if (workS1) {
Wr = workS1;
else {
if (PetscMin(W - workS1, workS1 - W) < (eV-sV)*eU) {
SETERRQ(PETSC_COMM_SELF,1, "Consistency broken!");
}
} else {
ierr = PetscMalloc(sizeof(PetscScalar)*(eV-sV)*eU, &Wr);
CHKERRQ(ierr);
}
658,9 → 669,12
/* Lower triangular M matrix */
} else if (DVD_ISNOT(sM,DVD_MAT_UTRIANG) &&
DVD_IS(sM,DVD_MAT_LTRIANG)) {
if (workS1)
if (workS1) {
Wr = workS1;
else {
if (PetscMin(W - workS1, workS1 - W) < (eU-sU)*eV) {
SETERRQ(PETSC_COMM_SELF,1, "Consistency broken!");
}
} else {
ierr = PetscMalloc(sizeof(PetscScalar)*(eU-sU)*eV, &Wr);
CHKERRQ(ierr);
}
/branches/slepc-3_2-proj/src/eps/impls/davidson/common/dvd_calcpairs.c
90,7 → 90,7
d->size_real_V = b->size_V;
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_AV = b->size_V;
d->size_real_AV = b->max_size_V+b->max_size_P+b->max_size_cP;
d->size_BDS = 0;
if (d->B && her_probl && (orth == DVD_ORTHOV_I || orth == DVD_ORTHOV_BOneMV)) {
d->size_real_BV = b->size_V; d->BV_shift = PETSC_TRUE;
108,6 → 108,8
/* pX, pY? */
b->max_nev*b->max_nev*(her_probl?0:(std_probl?1:2));
/* cS?, cT? */
b->max_size_auxV = PetscMax(b->max_size_auxV, b->max_size_X);
/* updateV0 */
b->max_size_auxS = PetscMax(PetscMax(
b->max_size_auxS,
b->max_size_X*b->max_size_proj*2*(std_probl?1:2) + /* updateAV1,BV1 */
129,20 → 131,19
/* Setup the step */
if (b->state >= DVD_STATE_CONF) {
d->max_cX_in_proj = cX_proj;
d->real_V = b->free_vecs; b->free_vecs+= b->size_V;
d->max_size_P = b->max_size_P;
d->real_V = b->free_vecs; b->free_vecs+= d->size_real_V;
if (harm) {
d->size_real_W = d->max_size_V+(d->B?b->max_nev:0);
d->real_W = b->free_vecs; b->free_vecs+= d->size_real_W;
} else {
d->real_W = PETSC_NULL;
}
d->real_AV = d->AV = b->free_vecs; b->free_vecs+= b->size_V;
d->max_size_AV = b->size_V;
d->real_AV = d->AV = b->free_vecs; b->free_vecs+= d->size_real_AV;
d->max_size_proj = b->max_size_proj;
d->real_H = b->free_scalars; b->free_scalars+= b->max_size_proj*b->max_size_proj;
d->ldH = b->max_size_proj;
d->pX = b->free_scalars; b->free_scalars+= b->max_size_V*b->max_size_proj;
d->S = b->free_scalars; b->free_scalars+= b->max_size_V*b->max_size_proj;
d->pX = b->free_scalars; b->free_scalars+= b->max_size_proj*b->max_size_proj;
d->S = b->free_scalars; b->free_scalars+= b->max_size_proj*b->max_size_proj;
if (!her_probl) {
d->cS = b->free_scalars; b->free_scalars+= b->max_nev*b->max_nev;
d->max_size_cS = b->max_nev;
209,14 → 210,14
 
d->size_V = 0;
d->tX = d->real_V;
d->V = d->real_V+d->max_size_cX_in_impr;
d->cX= d->real_V;
d->V = d->real_V+d->max_size_P;
d->cX = d->real_V;
d->size_cX = 0;
d->max_size_V = d->size_real_V-d->max_size_cX_in_impr;
d->max_size_V = d->size_real_V-d->max_size_P;
d->size_AV = 0;
d->KZ = d->tKZ = d->real_AV;
d->AV = d->real_AV+d->max_size_cX_in_impr;
d->max_size_AV = d->size_real_AV-d->max_size_cX_in_impr;
d->AV = d->real_AV+d->max_size_cX_in_impr+d->max_size_P;
d->max_size_AV = d->size_real_AV-d->max_size_cX_in_impr-d->max_size_P;
d->size_H = 0;
d->H = d->real_H;
for (i=0; i<d->max_size_cS*d->max_size_cS; i++) d->cS[i] = 0.0;
249,21 → 250,26
{
PetscErrorCode ierr;
DvdReduction r;
static const PetscInt MAX_OPS = 7;
DvdReductionChunk
ops[4];
ops[MAX_OPS];
DvdMult_copy_func
sr[4], *sr0 = sr;
PetscInt t=d->V_tra_e-d->V_tra_s,n=d->size_V+d->V_new_e-d->V_new_s,
size_in = 2*((d->size_cX+t)*t+n*n), i;
PetscScalar *in = d->auxS, *out = in+size_in;
sr[MAX_OPS], *sr0 = sr;
PetscInt size_in, i;
PetscScalar *in = d->auxS, *out;
 
PetscFunctionBegin;
 
size_in = (d->size_cX+d->V_tra_s)*d->V_tra_s*(d->cT?2:(d->cS?1:0)) +
d->size_H*d->V_tra_e*(d->BV?2:1) +
((d->size_H+d->V_new_e-d->V_new_s)*2+(d->V_new_e-d->V_new_s)*(d->V_new_e-d->V_new_s))*(d->BV?2:1);
out = in+size_in;
 
/* Check consistency */
if (2*size_in > d->size_auxS) { SETERRQ(PETSC_COMM_SELF,1, "Consistency broken!"); }
 
/* Prepare reductions */
ierr = SlepcAllReduceSumBegin(ops, 4, in, out, size_in, &r,
ierr = SlepcAllReduceSumBegin(ops, MAX_OPS, in, out, size_in, &r,
((PetscObject)d->V[0])->comm); CHKERRQ(ierr);
/* Allocate size_in */
d->auxS+= size_in;
334,6 → 340,7
 
#undef __FUNCT__
#define __FUNCT__ "dvd_calcpairs_updateV0"
/* auxV: V_tra_s, DvdMult_copy_func: 1 */
PetscErrorCode dvd_calcpairs_updateV0(dvdDashboard *d, DvdReduction *r,
DvdMult_copy_func **sr)
{
344,11 → 351,11
if (d->V_tra_s == 0 && d->V_tra_e == 0) { PetscFunctionReturn(0); }
 
/* cX <- [cX V*MT(0:V_tra_s-1)], V <- V*MT(V_tra_s:V_tra_e) */
if (d->V_tra_s > 0 && &d->cX[d->nconv] != d->V) {
if (d->V_tra_s > 0 && &d->cX[d->size_cX] != d->V) {
/* if cX and V are not contiguous */
ierr = SlepcUpdateVectorsZ(&d->cX[d->nconv], 0.0, 1.0, d->V, d->size_V,
ierr = SlepcUpdateVectorsZ(&d->cX[d->size_cX], 0.0, 1.0, d->V, d->size_V,
d->MTX, d->ldMTX, d->size_MT, d->V_tra_s); CHKERRQ(ierr);
ierr = SlepcUpdateVectorsZ(d->V, 0.0, 1.0, d->V, d->size_V,
ierr = SlepcUpdateVectorsZ(&d->V[d->V_tra_s], 0.0, 1.0, d->V, d->size_V,
&d->MTX[d->ldMTX*d->V_tra_s], d->ldMTX, d->size_MT,
d->V_tra_e-d->V_tra_s); CHKERRQ(ierr);
} else {
359,6 → 366,7
d->tX+= d->V_tra_s;
d->max_size_V-= d->V_tra_s;
d->size_V = d->V_tra_e-d->V_tra_s;
d->size_cX+= d->V_tra_s;
 
/* Udpate cS for standard problems */
if (d->cS && !d->cT && d->V_tra_s > 0 && !d->cY) {
395,7 → 403,7
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)) */
if (d->ipV_oneMV) {
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,
d->size_cX, d->V, d->BV, d->V_new_s, d->V_new_e,
d->auxS, d->eps->rand); CHKERRQ(ierr);
411,6 → 419,7
 
#undef __FUNCT__
#define __FUNCT__ "dvd_calcpairs_updateW0"
/* auxV: V_tra_s, DvdMult_copy_func: 2 */
PetscErrorCode dvd_calcpairs_updateW0(dvdDashboard *d, DvdReduction *r,
DvdMult_copy_func **sr)
{
454,6 → 463,7
 
#undef __FUNCT__
#define __FUNCT__ "dvd_calcpairs_updateW1"
/* auxS: size_cX+V_new_e+1 */
PetscErrorCode dvd_calcpairs_updateW1(dvdDashboard *d)
{
PetscErrorCode ierr;
479,6 → 489,7
 
#undef __FUNCT__
#define __FUNCT__ "dvd_calcpairs_updateAV0"
/* auxS: size_H*(V_tra_e-V_tra_s) */
PetscErrorCode dvd_calcpairs_updateAV0(dvdDashboard *d)
{
PetscErrorCode ierr;
520,7 → 531,7
MTY(tra_s:)'*auxS ] */
ierr = SlepcDenseCopyTriang(d->H, d->sH, d->ldH, d->auxS, 0, d->ldH,
cp, cMT); CHKERRQ(ierr);
ierr = SlepcDenseMatProdTriang(d->H, d->sH, d->ldH,
ierr = SlepcDenseMatProdTriang(&d->H[cp], d->sH, d->ldH,
&MTY[d->ldMTX*tra_s], 0, d->ldMTX, d->size_MT, cMT, PETSC_TRUE,
d->auxS, 0, d->ldH, d->size_H, cp+cMT, PETSC_FALSE);
CHKERRQ(ierr);
533,6 → 544,7
 
#undef __FUNCT__
#define __FUNCT__ "dvd_calcpairs_updateAV1"
/* DvdMult_copy_func: 2 */
PetscErrorCode dvd_calcpairs_updateAV1(dvdDashboard *d, DvdReduction *r,
DvdMult_copy_func **sr)
{
554,7 → 566,7
being old=0:V_new_s-1, new=V_new_s:V_new_e-1 */
V_new_s = d->V_new_s+d->cX_in_H;
V_tra_s = d->V_tra_s+d->cX_in_H;
if (&d->cX[d->nconv] != d->V && W == d->V) {
if (&d->cX[d->size_cX] != d->V && W == d->V) {
/* if W=[cX V] is not contiguous */
ierr = VecsMultS(d->H, d->sH, d->ldH,
&d->cX[d->size_cX-d->cX_in_H], 0, d->cX_in_H,
564,7 → 576,7
&d->AV[-d->size_cX], V_new_s, d->V_new_e, r, (*sr)++); CHKERRQ(ierr);
} else {
ierr = VecsMultS(d->H, d->sH, d->ldH, W-d->cX_in_H, V_new_s, d->V_new_e,
d->V-d->cX_in_H, V_new_s, d->V_new_e, r, (*sr)++); CHKERRQ(ierr);
d->AV-d->cX_in_H, V_new_s, d->V_new_e, r, (*sr)++); CHKERRQ(ierr);
}
d->size_H = d->V_new_e;
 
573,6 → 585,7
 
#undef __FUNCT__
#define __FUNCT__ "dvd_calcpairs_updateBV0"
/* auxS: max(BcX*(size_cX+V_new_e+1), size_G*(V_tra_e-V_tra_s)) */
PetscErrorCode dvd_calcpairs_updateBV0(dvdDashboard *d)
{
PetscErrorCode ierr;
634,6 → 647,7
 
#undef __FUNCT__
#define __FUNCT__ "dvd_calcpairs_updateBV1"
/* DvdMult_copy_func: 2 */
PetscErrorCode dvd_calcpairs_updateBV1(dvdDashboard *d, DvdReduction *r,
DvdMult_copy_func **sr)
{
650,7 → 664,7
being old=0:V_new_s-1, new=V_new_s:V_new_e-1 */
V_new_s = d->V_new_s+d->cX_in_G;
V_tra_s = d->V_tra_s+d->cX_in_G;
if (&d->cX[d->nconv] != d->V && W == d->V) {
if (&d->cX[d->size_cX] != d->V && W == d->V) {
/* if W=[cX V] is not contiguous */
ierr = VecsMultS(d->G, d->sG, d->ldH,
&d->cX[d->size_cX-d->cX_in_G], 0, d->cX_in_G,
660,7 → 674,7
&d->BV[-d->size_cX], V_new_s, d->V_new_e, r, (*sr)++); CHKERRQ(ierr);
} else {
ierr = VecsMultS(d->G, d->sG, d->ldH, W-d->cX_in_G, V_new_s, d->V_new_e,
d->V-d->cX_in_G, V_new_s, d->V_new_e, r, (*sr)++); CHKERRQ(ierr);
d->BV-d->cX_in_G, V_new_s, d->V_new_e, r, (*sr)++); CHKERRQ(ierr);
}
d->size_G = d->V_new_e;
 
796,7 → 810,7
 
PetscFunctionBegin;
 
ierr = EPSSortDenseHEP(d->eps, d->size_H, d->cX_in_H, d->eigr, d->pX, d->ldpX);
ierr = EPSSortDenseHEP(d->eps, d->size_H, d->cX_in_H, d->eigr-d->cX_in_H, d->pX, d->ldpX);
CHKERRQ(ierr);
 
if (d->calcpairs_eigs_trans) {
/branches/slepc-3_2-proj/src/eps/impls/davidson/common/dvd_updatev.c
52,7 → 52,6
*oldU, /* previous projected right igenvectors */
*oldV; /* previous projected left eigenvectors */
PetscInt
size_real_V, /* original size of real_* */
ldoldU, /* leading dimension of oldU */
size_oldU; /* size of oldU */
PetscBool
102,11 → 101,11
/* Setup the step */
if (b->state >= DVD_STATE_CONF) {
ierr = PetscMalloc(sizeof(dvdManagV_basic), &data); CHKERRQ(ierr);
data->mpd = b->max_size_V;
data->min_size_V = min_size_V;
d->bs = bs;
d->max_size_X = b->max_size_X;
data->plusk = plusk;
data->size_real_V = b->size_V;
data->allResiduals = allResiduals;
 
d->real_eigr = b->free_scalars; b->free_scalars+= b->size_V;
158,12 → 157,13
d->eigr = d->ceigr = d->real_eigr;
d->eigi = d->ceigi = d->real_eigi;
#if defined(PETSC_USE_COMPLEX)
for(i=0; i<data->size_real_V; i++) d->eigi[i] = 0.0;
for(i=0; i<d->size_real_V; i++) d->eigi[i] = 0.0;
#endif
d->nR = d->real_nR;
for(i=0; i<data->size_real_V; i++) d->nR[i] = PETSC_MAX_REAL;
for(i=0; i<d->size_real_V; i++) d->nR[i] = PETSC_MAX_REAL;
d->nX = d->real_nX;
d->errest = d->real_errest;
for(i=0; i<d->size_real_V; i++) d->errest[i] = PETSC_MAX_REAL;
data->ldoldU = 0;
data->oldV = PETSC_NULL;
d->ldMTY = 0;
189,8 → 189,6
restart = (d->size_V + d->max_size_X > PetscMin(data->mpd,d->max_size_V))?
PETSC_TRUE:PETSC_FALSE;
 
#undef ONE
 
/* Check old isRestarting function */
if (!restart && data->old_isRestarting)
restart = data->old_isRestarting(d);
230,6 → 228,7
 
/* If the subspaces doesn't need restart, add new vector */
if (!d->isRestarting(d)) {
d->size_D = 0;
ierr = dvd_updateV_update_gen(d); CHKERRQ(ierr);
 
/* If some vector were add, exit */
458,6 → 457,7
&size_D); CHKERRQ(ierr);
 
/* If D is empty, exit */
d->size_D = size_D;
if (size_D == 0) { PetscFunctionReturn(0); }
 
/* Get the converged pairs */
468,7 → 468,6
/* Notify the changes in V */
d->V_tra_s = 0; d->V_tra_e = 0;
d->V_new_s = d->size_V; d->V_new_e = d->size_V+size_D;
d->size_D = size_D;
 
/* Save the projected eigenvectors */
if (data->plusk > 0) {
/branches/slepc-3_2-proj/src/eps/impls/davidson/common/dvd_improvex.c
115,11 → 115,11
#if !defined(PETSC_USE_COMPLEX)
if (!herm) {
max_bs++;
b->max_size_X = PetscMax(b->max_size_X, max_bs);
b->max_size_P = PetscMax(b->max_size_P, 2);
} else
#endif
b->max_size_P = PetscMax(b->max_size_P, 1);
b->max_size_X = PetscMax(b->max_size_X, max_bs);
size_P = b->max_size_P+cX_impr;
b->max_size_auxV = PetscMax(b->max_size_auxV, b->max_size_X*(ksp?2:1));
/* kr?, auxV */
566,7 → 566,7
/* XKZ <- X'*KZ */
size_cX = PetscMin(d->size_cX, d->max_size_cX_in_impr);
data->s_cX = d->size_cX - size_cX;
wS0 = *auxS; wS1 = wS0+2*n*d->size_cX;
wS0 = *auxS; wS1 = wS0+2*n*d->size_cX+n*n;
ierr = VecsMult(&data->XKZ[size_cX*data->ldXKZ+size_cX], 0, data->ldXKZ,
&d->cX[data->s_cX], size_cX, size_cX+n,
&d->KZ[data->s_cX], size_cX, size_cX+n, wS0, wS1);CHKERRQ(ierr);
/branches/slepc-3_2-proj/src/eps/impls/davidson/common/davidson.c
284,10 → 284,10
CHKERRQ(ierr);
 
/* Associate the eigenvalues to the EPS */
eps->eigr = dvd->eigr;
eps->eigi = dvd->eigi;
eps->errest = dvd->errest;
eps->V = dvd->V;
eps->eigr = dvd->real_eigr;
eps->eigi = dvd->real_eigi;
eps->errest = dvd->real_errest;
eps->V = dvd->real_V;
 
PetscFunctionReturn(0);
/branches/slepc-3_2-proj/src/eps/impls/davidson/common/davidson.h
176,6 → 176,7
/* max vectors from cX in the projected problem */
max_size_cX_in_impr,
/* max vectros from cX in the projector */
max_size_P, /* max unconverged vectors in the projector */
bs; /* max vectors that expands the subspace every iteration */
EPS eps; /* Connection to SLEPc */