| Line 190... |
Line 190... |
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "EPSKrylovSchur_Slice"
|
#define __FUNCT__ "EPSKrylovSchur_Slice"
|
static PetscErrorCode EPSKrylovSchur_Slice(EPS eps)
|
static PetscErrorCode EPSKrylovSchur_Slice(EPS eps)
|
{
|
{
|
PetscErrorCode ierr;
|
PetscErrorCode ierr;
|
PetscInt i,conv,k,l,ld,nv,m,*iwork,p,j;
|
PetscInt i,conv,k,l,ld,nv,*iwork,p,j;
|
Vec u=eps->work[0];
|
Vec u=eps->work[0];
|
PetscScalar *Q,*A,nu,rtmp,alpha;
|
PetscScalar *Q,*A,nu,rtmp,alpha;
|
PetscReal *a,*b,beta;
|
PetscReal *a,*b,beta;
|
PetscBool breakdown;
|
PetscBool breakdown;
|
PetscInt count0,count1;
|
PetscInt count0,count1;
|
| Line 242... |
Line 242... |
}
|
}
|
/* Restart loop */
|
/* Restart loop */
|
while (eps->reason == EPS_CONVERGED_ITERATING) {
|
while (eps->reason == EPS_CONVERGED_ITERATING) {
|
eps->its++; sr->itsKs++;
|
eps->its++; sr->itsKs++;
|
/* Compute an nv-step Lanczos factorization */
|
/* Compute an nv-step Lanczos factorization */
|
m = PetscMin(eps->nconv+eps->mpd,eps->ncv);
|
nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
|
ierr = PSGetArrayReal(eps->ps,PS_MAT_T,&a);CHKERRQ(ierr);
|
ierr = PSGetArrayReal(eps->ps,PS_MAT_T,&a);CHKERRQ(ierr);
|
b = a + ld;
|
b = a + ld;
|
ierr = EPSFullLanczos(eps,a+l,b+l,eps->V,eps->nconv+l,&m,u,&breakdown);CHKERRQ(ierr);
|
ierr = EPSFullLanczos(eps,a,b,eps->V,eps->nconv+l,&nv,u,&breakdown);CHKERRQ(ierr);
|
if(breakdown){/* explicit purification*/
|
if(breakdown){/* explicit purification*/
|
sPres->expf = PETSC_TRUE;
|
sPres->expf = PETSC_TRUE;
|
}
|
}
|
nv = m - eps->nconv;
|
|
beta = b[nv-1];
|
beta = b[nv-1];
|
ierr = PSRestoreArrayReal(eps->ps,PS_MAT_T,&a);CHKERRQ(ierr);
|
ierr = PSRestoreArrayReal(eps->ps,PS_MAT_T,&a);CHKERRQ(ierr);
|
ierr = PSSetDimensions(eps->ps,nv,0,l);CHKERRQ(ierr);
|
ierr = PSSetDimensions(eps->ps,nv,eps->nconv,eps->nconv+l);CHKERRQ(ierr);
|
if (l==0) {
|
if (l==0) {
|
ierr = PSSetState(eps->ps,PS_STATE_INTERMEDIATE);CHKERRQ(ierr);
|
ierr = PSSetState(eps->ps,PS_STATE_INTERMEDIATE);CHKERRQ(ierr);
|
} else {
|
} else {
|
ierr = PSSetState(eps->ps,PS_STATE_RAW);CHKERRQ(ierr);
|
ierr = PSSetState(eps->ps,PS_STATE_RAW);CHKERRQ(ierr);
|
}
|
}
|
| Line 264... |
Line 263... |
/* Solve projected problem and compute residual norm estimates */
|
/* Solve projected problem and compute residual norm estimates */
|
if(eps->its == 1 && l > 0){/* After rational update */
|
if(eps->its == 1 && l > 0){/* After rational update */
|
ierr = PSGetArray(eps->ps,PS_MAT_A,&A);CHKERRQ(ierr);
|
ierr = PSGetArray(eps->ps,PS_MAT_A,&A);CHKERRQ(ierr);
|
ierr = PSGetArrayReal(eps->ps,PS_MAT_T,&a);CHKERRQ(ierr);
|
ierr = PSGetArrayReal(eps->ps,PS_MAT_T,&a);CHKERRQ(ierr);
|
b = a + ld;
|
b = a + ld;
|
A[l*ld+l-1] = A[(l-1)*ld+l];
|
k = eps->nconv+l;
|
A[l*ld+l] = a[l];
|
A[k*ld+k-1] = A[(k-1)*ld+k];
|
for(j=l+1; j< nv; j++){
|
A[k*ld+k] = a[k];
|
|
for(j=k+1; j< nv; j++){
|
A[j*ld+j] = a[j];
|
A[j*ld+j] = a[j];
|
A[j*ld+j-1] = b[j-1] ;
|
A[j*ld+j-1] = b[j-1] ;
|
A[(j-1)*ld+j] = b[j-1];
|
A[(j-1)*ld+j] = b[j-1];
|
}
|
}
|
ierr = PSRestoreArray(eps->ps,PS_MAT_A,&A);CHKERRQ(ierr);
|
ierr = PSRestoreArray(eps->ps,PS_MAT_A,&A);CHKERRQ(ierr);
|
ierr = PSRestoreArrayReal(eps->ps,PS_MAT_T,&a);CHKERRQ(ierr);
|
ierr = PSRestoreArrayReal(eps->ps,PS_MAT_T,&a);CHKERRQ(ierr);
|
ierr = PSSolve(eps->ps,eps->eigr+eps->nconv,PETSC_NULL);CHKERRQ(ierr);
|
ierr = PSSolve(eps->ps,eps->eigr,PETSC_NULL);CHKERRQ(ierr);
|
|
ierr = PSSetCompact(eps->ps,PETSC_TRUE);
|
}else{/* Restart */
|
}else{/* Restart */
|
ierr = PSSolve(eps->ps,eps->eigr+eps->nconv,PETSC_NULL);CHKERRQ(ierr);
|
ierr = PSSolve(eps->ps,eps->eigr,PETSC_NULL);CHKERRQ(ierr);
|
}
|
}
|
ierr = PSSort(eps->ps,eps->eigr+eps->nconv,PETSC_NULL,eps->which_func,eps->which_ctx);CHKERRQ(ierr);
|
ierr = PSSort(eps->ps,eps->eigr,PETSC_NULL,eps->which_func,eps->which_ctx);CHKERRQ(ierr);
|
/* Residual */
|
/* Residual */
|
ierr = PSGetArray(eps->ps,PS_MAT_Q,&Q);CHKERRQ(ierr);
|
ierr = EPSKrylovConvergence(eps,PETSC_TRUE,eps->nconv,nv-eps->nconv,eps->V,nv,beta,1.0,&k);CHKERRQ(ierr);
|
ierr = EPSKrylovConvergence(eps,PETSC_TRUE,eps->nconv,nv,eps->V+eps->nconv,nv,beta,1.0,&k);CHKERRQ(ierr);
|
|
/* Check convergence */
|
/* Check convergence */
|
|
|
ierr = PSGetArrayReal(eps->ps,PS_MAT_T,&a);CHKERRQ(ierr);
|
ierr = PSGetArrayReal(eps->ps,PS_MAT_T,&a);CHKERRQ(ierr);
|
b = a + ld;
|
b = a + ld;
|
conv=k=j=0;
|
conv = 0;
|
for(i=0;i<nv;i++)if(eps->errest[eps->nconv+i] < eps->tol)conv++;
|
j = k = eps->nconv;
|
for(i=0;i<nv;i++){
|
for(i=eps->nconv;i<nv;i++)if(eps->errest[i] < eps->tol)conv++;
|
if(eps->errest[eps->nconv+i] < eps->tol){
|
for(i=eps->nconv;i<nv;i++){
|
|
if(eps->errest[i] < eps->tol){
|
iwork[j++]=i;
|
iwork[j++]=i;
|
}else iwork[conv+k++]=i;
|
}else iwork[conv+k++]=i;
|
}
|
}
|
for(i=0;i<nv;i++){
|
for(i=eps->nconv;i<nv;i++){
|
a[i]=PetscRealPart(eps->eigr[eps->nconv+i]);
|
a[i]=PetscRealPart(eps->eigr[i]);
|
b[i]=eps->errest[eps->nconv+i];
|
b[i]=eps->errest[i];
|
}
|
}
|
for(i=0;i<nv;i++){
|
for(i=eps->nconv;i<nv;i++){
|
eps->eigr[eps->nconv+i] = a[iwork[i]];
|
eps->eigr[i] = a[iwork[i]];
|
eps->errest[eps->nconv+i] = b[iwork[i]];
|
eps->errest[i] = b[iwork[i]];
|
}
|
}
|
for( i=0;i<nv;i++){
|
for(i=eps->nconv;i<nv;i++){
|
p=iwork[i];
|
a[i]=PetscRealPart(eps->eigr[i]);
|
if(p!=i){
|
b[i]=eps->errest[i];
|
j=i+1;
|
|
while(iwork[j]!=i)j++;
|
|
iwork[j]=p;iwork[i]=i;
|
|
for(k=0;k<nv;k++){
|
|
rtmp=Q[k+p*ld];Q[k+p*ld]=Q[k+i*ld];Q[k+i*ld]=rtmp;
|
|
}
|
|
}
|
|
}
|
}
|
|
|
k=eps->nconv+conv;
|
k=eps->nconv+conv;
|
ierr = PSRestoreArray(eps->ps,PS_MAT_Q,&Q);CHKERRQ(ierr);
|
ierr = PSPermuteColumns_Private(eps->ps,eps->nconv,nv,PS_MAT_Q,iwork);CHKERRQ(ierr);
|
ierr = PSRestoreArrayReal(eps->ps,PS_MAT_T,&a);CHKERRQ(ierr);
|
ierr = PSRestoreArrayReal(eps->ps,PS_MAT_T,&a);CHKERRQ(ierr);
|
|
|
/* Checking values obtained for completing */
|
/* Checking values obtained for completing */
|
for(i=0;i<k;i++){
|
for(i=0;i<k;i++){
|
sr->back[i]=eps->eigr[i];
|
sr->back[i]=eps->eigr[i];
|
| Line 346... |
Line 342... |
iterCompl = sr->iterCompl;
|
iterCompl = sr->iterCompl;
|
}else eps->reason = EPS_CONVERGED_TOL;
|
}else eps->reason = EPS_CONVERGED_TOL;
|
}
|
}
|
}
|
}
|
/* Update l */
|
/* Update l */
|
if(eps->reason == EPS_CONVERGED_ITERATING )l = (eps->nconv+nv-k)/2;
|
if(eps->reason == EPS_CONVERGED_ITERATING )l = (nv-k)/2;
|
else l=eps->nconv+nv-k;
|
else l=nv-k;
|
if(breakdown)l=0;
|
if(breakdown)l=0;
|
|
|
if (eps->reason == EPS_CONVERGED_ITERATING) {
|
if (eps->reason == EPS_CONVERGED_ITERATING) {
|
if (breakdown) {
|
if (breakdown) {
|
/* Start a new Lanczos factorization */
|
/* Start a new Lanczos factorization */
|
| Line 364... |
Line 360... |
} else {
|
} else {
|
/* Prepare the Rayleigh quotient for restart */
|
/* Prepare the Rayleigh quotient for restart */
|
ierr = PSGetArrayReal(eps->ps,PS_MAT_T,&a);CHKERRQ(ierr);
|
ierr = PSGetArrayReal(eps->ps,PS_MAT_T,&a);CHKERRQ(ierr);
|
ierr = PSGetArray(eps->ps,PS_MAT_Q,&Q);CHKERRQ(ierr);
|
ierr = PSGetArray(eps->ps,PS_MAT_Q,&Q);CHKERRQ(ierr);
|
b = a + ld;
|
b = a + ld;
|
for (i=0;i<l;i++) {
|
for (i=k;i<k+l;i++) {
|
a[i] = PetscRealPart(eps->eigr[i+k]);
|
a[i] = PetscRealPart(eps->eigr[i]);
|
b[i] = PetscRealPart(Q[nv-1+(i+k-eps->nconv)*ld]*beta);
|
b[i] = PetscRealPart(Q[nv-1+i*ld]*beta);
|
}
|
}
|
ierr = PSRestoreArrayReal(eps->ps,PS_MAT_T,&a);CHKERRQ(ierr);
|
ierr = PSRestoreArrayReal(eps->ps,PS_MAT_T,&a);CHKERRQ(ierr);
|
ierr = PSRestoreArray(eps->ps,PS_MAT_Q,&Q);CHKERRQ(ierr);
|
ierr = PSRestoreArray(eps->ps,PS_MAT_Q,&Q);CHKERRQ(ierr);
|
}
|
}
|
}
|
}
|
/* Update the corresponding vectors V(:,idx) = V*Q(:,idx) */
|
/* Update the corresponding vectors V(:,idx) = V*Q(:,idx) */
|
ierr = PSGetArray(eps->ps,PS_MAT_Q,&Q);CHKERRQ(ierr);
|
ierr = PSGetArray(eps->ps,PS_MAT_Q,&Q);CHKERRQ(ierr);
|
ierr = SlepcUpdateVectors(nv,eps->V+eps->nconv,0,k+l-eps->nconv,Q,ld,PETSC_FALSE);CHKERRQ(ierr);
|
ierr = SlepcUpdateVectors(nv,eps->V,eps->nconv,k+l,Q,ld,PETSC_FALSE);CHKERRQ(ierr);
|
/* Purification */
|
/* Purification */
|
|
|
if(!sPres->expf){/* u not saved if breakdown */
|
if(!sPres->expf){/* u not saved if breakdown */
|
for(i=eps->nconv; i<k;i++){
|
for(i=eps->nconv; i<k;i++){
|
alpha = (Q[nv-1+(i-eps->nconv)*ld])/PetscRealPart(eps->eigr[i]);
|
alpha = (Q[nv-1+i*ld])/PetscRealPart(eps->eigr[i]);
|
ierr = VecAXPY(eps->V[i], alpha, u);CHKERRQ(ierr);
|
ierr = VecAXPY(eps->V[i], alpha, u);CHKERRQ(ierr);
|
}
|
}
|
}
|
}
|
|
|
ierr = PSRestoreArray(eps->ps,PS_MAT_Q,&Q);CHKERRQ(ierr);
|
ierr = PSRestoreArray(eps->ps,PS_MAT_Q,&Q);CHKERRQ(ierr);
|
/* Normalize u and append it to V */
|
/* Normalize u and append it to V */
|
if ( eps->reason == EPS_CONVERGED_ITERATING && !breakdown) {
|
if ( eps->reason == EPS_CONVERGED_ITERATING && !breakdown) {
|
ierr = VecAXPBY(eps->V[k+l],1.0/beta,0.0,u);CHKERRQ(ierr);
|
ierr = VecAXPBY(eps->V[k+l],1.0/beta,0.0,u);CHKERRQ(ierr);
|
}
|
}
|
/* Monitor */
|
/* Monitor */
|
if(eps->numbermonitors >0){
|
if(eps->numbermonitors >0){
|
aux = auxc = 0;
|
aux = auxc = 0;
|
for(i=0;i<nv+eps->nconv;i++){
|
for(i=0;i<nv;i++){
|
sr->back[i]=eps->eigr[i];
|
sr->back[i]=eps->eigr[i];
|
}
|
}
|
ierr = STBackTransform(eps->OP,nv+eps->nconv,sr->back,eps->eigi);CHKERRQ(ierr);
|
ierr = STBackTransform(eps->OP,nv,sr->back,eps->eigi);CHKERRQ(ierr);
|
for(i=0;i<nv+eps->nconv;i++){
|
for(i=0;i<nv;i++){
|
lambda = PetscRealPart(sr->back[i]);
|
lambda = PetscRealPart(sr->back[i]);
|
if( ((sr->dir)*(lambda - sPres->ext[0]) > 0)&& ((sr->dir)*(sPres->ext[1] - lambda) > 0)){
|
if( ((sr->dir)*(lambda - sPres->ext[0]) > 0)&& ((sr->dir)*(sPres->ext[1] - lambda) > 0)){
|
sr->monit[sr->indexEig+aux]=eps->eigr[i];
|
sr->monit[sr->indexEig+aux]=eps->eigr[i];
|
sr->errest[sr->indexEig+aux]=eps->errest[i];
|
sr->errest[sr->indexEig+aux]=eps->errest[i];
|
aux++;
|
aux++;
|
if(eps->errest[i] < eps->tol)auxc++;
|
if(eps->errest[i] < eps->tol)auxc++;
|
}
|
}
|
}
|
}
|
ierr = EPSMonitor(eps,eps->its,auxc+sr->indexEig,sr->monit,sr->eigi,sr->errest,sr->indexEig+aux);CHKERRQ(ierr);
|
ierr = EPSMonitor(eps,eps->its,auxc+sr->indexEig,sr->monit,sr->eigi,sr->errest,sr->indexEig+aux);CHKERRQ(ierr);
|
}
|
}
|
conv = k - eps->nconv;
|
|
eps->nconv = k;
|
eps->nconv = k;
|
if(eps->reason != EPS_CONVERGED_ITERATING){
|
if(eps->reason != EPS_CONVERGED_ITERATING){
|
/* Store approximated values for next shift */
|
/* Store approximated values for next shift */
|
ierr = PSGetArray(eps->ps,PS_MAT_Q,&Q);CHKERRQ(ierr);
|
ierr = PSGetArray(eps->ps,PS_MAT_Q,&Q);CHKERRQ(ierr);
|
sr->nS = l;
|
sr->nS = l;
|
for (i=0;i<l;i++) {
|
for (i=0;i<l;i++) {
|
sr->S[i] = eps->eigr[i+k];/* Diagonal elements */
|
sr->S[i] = eps->eigr[i+k];/* Diagonal elements */
|
sr->S[i+l] = Q[nv-1+(i+conv)*ld]*beta; /* Out of diagonal elements */
|
sr->S[i+l] = Q[nv-1+(i+k)*ld]*beta; /* Out of diagonal elements */
|
}
|
}
|
sr->beta = beta;
|
sr->beta = beta;
|
ierr = PSRestoreArray(eps->ps,PS_MAT_Q,&Q);CHKERRQ(ierr);
|
ierr = PSRestoreArray(eps->ps,PS_MAT_Q,&Q);CHKERRQ(ierr);
|
}
|
}
|
}
|
}
|