| Line 175... |
Line 175... |
#define __FUNCT__ "EPSSolve_Subspace"
|
#define __FUNCT__ "EPSSolve_Subspace"
|
PetscErrorCode EPSSolve_Subspace(EPS eps)
|
PetscErrorCode EPSSolve_Subspace(EPS eps)
|
{
|
{
|
PetscErrorCode ierr;
|
PetscErrorCode ierr;
|
EPS_SUBSPACE *ctx = (EPS_SUBSPACE*)eps->data;
|
EPS_SUBSPACE *ctx = (EPS_SUBSPACE*)eps->data;
|
PetscInt i,k,ngrp,nogrp,*itrsd,*itrsdold,
|
PetscInt i,k,ld,ngrp,nogrp,*itrsd,*itrsdold,
|
nxtsrr,idsrr,idort,nxtort,nv,ncv = eps->ncv,its;
|
nxtsrr,idsrr,idort,nxtort,nv,ncv = eps->ncv,its;
|
PetscScalar *T,*U;
|
PetscScalar *T,*U;
|
PetscReal arsd,oarsd,ctr,octr,ae,oae,*rsd,norm,tcond=1.0;
|
PetscReal arsd,oarsd,ctr,octr,ae,oae,*rsd,norm,tcond=1.0;
|
PetscBool breakdown;
|
PetscBool breakdown;
|
/* Parameters */
|
/* Parameters */
|
| Line 216... |
Line 216... |
|
|
while (eps->its<eps->max_it) {
|
while (eps->its<eps->max_it) {
|
eps->its++;
|
eps->its++;
|
nv = PetscMin(eps->nconv+eps->mpd,ncv);
|
nv = PetscMin(eps->nconv+eps->mpd,ncv);
|
ierr = PSSetDimensions(eps->ps,nv,eps->nconv,0);CHKERRQ(ierr);
|
ierr = PSSetDimensions(eps->ps,nv,eps->nconv,0);CHKERRQ(ierr);
|
|
ierr = PSGetLeadingDimension(eps->ps,&ld);CHKERRQ(ierr);
|
|
|
/* Find group in previously computed eigenvalues */
|
/* Find group in previously computed eigenvalues */
|
ierr = EPSSubspaceFindGroup(eps->nconv,nv,eps->eigr,eps->eigi,rsd,grptol,&nogrp,&octr,&oae,&oarsd);CHKERRQ(ierr);
|
ierr = EPSSubspaceFindGroup(eps->nconv,nv,eps->eigr,eps->eigi,rsd,grptol,&nogrp,&octr,&oae,&oarsd);CHKERRQ(ierr);
|
|
|
/* AV(:,idx) = OP * V(:,idx) */
|
/* AV(:,idx) = OP * V(:,idx) */
|
| Line 228... |
Line 229... |
}
|
}
|
|
|
/* T(:,idx) = V' * AV(:,idx) */
|
/* T(:,idx) = V' * AV(:,idx) */
|
ierr = PSGetArray(eps->ps,PS_MAT_A,&T);CHKERRQ(ierr);
|
ierr = PSGetArray(eps->ps,PS_MAT_A,&T);CHKERRQ(ierr);
|
for (i=eps->nconv;i<nv;i++) {
|
for (i=eps->nconv;i<nv;i++) {
|
ierr = VecMDot(ctx->AV[i],nv,eps->V,T+i*ncv);CHKERRQ(ierr);
|
ierr = VecMDot(ctx->AV[i],nv,eps->V,T+i*ld);CHKERRQ(ierr);
|
}
|
}
|
ierr = PSRestoreArray(eps->ps,PS_MAT_A,&T);CHKERRQ(ierr);
|
ierr = PSRestoreArray(eps->ps,PS_MAT_A,&T);CHKERRQ(ierr);
|
ierr = PSSetState(eps->ps,PS_STATE_RAW);CHKERRQ(ierr);
|
ierr = PSSetState(eps->ps,PS_STATE_RAW);CHKERRQ(ierr);
|
|
|
/* Solve projected problem */
|
/* Solve projected problem */
|
ierr = PSSolve(eps->ps,eps->eigr,eps->eigi);CHKERRQ(ierr);
|
ierr = PSSolve(eps->ps,eps->eigr,eps->eigi);CHKERRQ(ierr);
|
ierr = PSSort(eps->ps,eps->eigr,eps->eigi,eps->which_func,eps->which_ctx);CHKERRQ(ierr);
|
ierr = PSSort(eps->ps,eps->eigr,eps->eigi,eps->which_func,eps->which_ctx);CHKERRQ(ierr);
|
|
|
/* Update vectors V(:,idx) = V * U(:,idx) */
|
/* Update vectors V(:,idx) = V * U(:,idx) */
|
ierr = PSGetArray(eps->ps,PS_MAT_Q,&U);CHKERRQ(ierr);
|
ierr = PSGetArray(eps->ps,PS_MAT_Q,&U);CHKERRQ(ierr);
|
ierr = SlepcUpdateVectors(nv,ctx->AV,eps->nconv,nv,U,nv,PETSC_FALSE);CHKERRQ(ierr);
|
ierr = SlepcUpdateVectors(nv,ctx->AV,eps->nconv,nv,U,ld,PETSC_FALSE);CHKERRQ(ierr);
|
ierr = SlepcUpdateVectors(nv,eps->V,eps->nconv,nv,U,nv,PETSC_FALSE);CHKERRQ(ierr);
|
ierr = SlepcUpdateVectors(nv,eps->V,eps->nconv,nv,U,ld,PETSC_FALSE);CHKERRQ(ierr);
|
ierr = PSRestoreArray(eps->ps,PS_MAT_Q,&U);CHKERRQ(ierr);
|
ierr = PSRestoreArray(eps->ps,PS_MAT_Q,&U);CHKERRQ(ierr);
|
|
|
/* Convergence check */
|
/* Convergence check */
|
ierr = PSGetArray(eps->ps,PS_MAT_A,&T);CHKERRQ(ierr);
|
ierr = PSGetArray(eps->ps,PS_MAT_A,&T);CHKERRQ(ierr);
|
ierr = EPSSubspaceResidualNorms(eps->V,ctx->AV,T,eps->nconv,nv,ncv,eps->work[0],rsd);CHKERRQ(ierr);
|
ierr = EPSSubspaceResidualNorms(eps->V,ctx->AV,T,eps->nconv,nv,ld,eps->work[0],rsd);CHKERRQ(ierr);
|
ierr = PSRestoreArray(eps->ps,PS_MAT_A,&T);CHKERRQ(ierr);
|
ierr = PSRestoreArray(eps->ps,PS_MAT_A,&T);CHKERRQ(ierr);
|
|
|
for (i=eps->nconv;i<nv;i++) {
|
for (i=eps->nconv;i<nv;i++) {
|
itrsdold[i] = itrsd[i];
|
itrsdold[i] = itrsd[i];
|
itrsd[i] = its;
|
itrsd[i] = its;
|