Actual source code: ks-indef.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2019, Universitat Politecnica de Valencia, Spain

  6:    This file is part of SLEPc.
  7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
  8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  9: */
 10: /*
 11:    SLEPc eigensolver: "krylovschur"

 13:    Method: Krylov-Schur for symmetric-indefinite eigenproblems
 14: */
 15: #include <slepc/private/epsimpl.h>
 16: #include "krylovschur.h"

 18: PetscErrorCode EPSSolve_KrylovSchur_Indefinite(EPS eps)
 19: {
 20:   PetscErrorCode  ierr;
 21:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
 22:   PetscInt        i,k,l,ld,nv,t,nconv=0;
 23:   Mat             U;
 24:   Vec             vomega,w=eps->work[0];
 25:   PetscScalar     *Q,*aux;
 26:   PetscReal       *a,*b,*r,beta,beta1=1.0,*omega;
 27:   PetscBool       breakdown=PETSC_FALSE,symmlost=PETSC_FALSE;

 30:   DSGetLeadingDimension(eps->ds,&ld);

 32:   /* Get the starting Lanczos vector */
 33:   EPSGetStartVector(eps,0,NULL);

 35:   /* Extract sigma[0] from BV, computed during normalization */
 36:   BVSetActiveColumns(eps->V,0,1);
 37:   VecCreateSeq(PETSC_COMM_SELF,1,&vomega);
 38:   BVGetSignature(eps->V,vomega);
 39:   VecGetArray(vomega,&aux);
 40:   DSGetArrayReal(eps->ds,DS_MAT_D,&omega);
 41:   omega[0] = PetscRealPart(aux[0]);
 42:   DSRestoreArrayReal(eps->ds,DS_MAT_D,&omega);
 43:   VecRestoreArray(vomega,&aux);
 44:   VecDestroy(&vomega);
 45:   l = 0;

 47:   /* Restart loop */
 48:   while (eps->reason == EPS_CONVERGED_ITERATING) {
 49:     eps->its++;

 51:     /* Compute an nv-step Lanczos factorization */
 52:     nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
 53:     DSGetArrayReal(eps->ds,DS_MAT_T,&a);
 54:     b = a + ld;
 55:     DSGetArrayReal(eps->ds,DS_MAT_D,&omega);
 56:     EPSPseudoLanczos(eps,a,b,omega,eps->nconv+l,&nv,&breakdown,&symmlost,NULL,w);
 57:     if (symmlost) {
 58:       eps->reason = EPS_DIVERGED_SYMMETRY_LOST;
 59:       if (nv==eps->nconv+l+1) { eps->nconv = nconv; break; }
 60:     }
 61:     beta = b[nv-1];
 62:     DSRestoreArrayReal(eps->ds,DS_MAT_T,&a);
 63:     DSRestoreArrayReal(eps->ds,DS_MAT_D,&omega);
 64:     DSSetDimensions(eps->ds,nv,0,eps->nconv,eps->nconv+l);
 65:     if (l==0) {
 66:       DSSetState(eps->ds,DS_STATE_INTERMEDIATE);
 67:     } else {
 68:       DSSetState(eps->ds,DS_STATE_RAW);
 69:     }
 70:     BVSetActiveColumns(eps->V,eps->nconv,nv);

 72:     /* Solve projected problem */
 73:     DSSolve(eps->ds,eps->eigr,eps->eigi);
 74:     DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL);
 75:     DSSynchronize(eps->ds,eps->eigr,eps->eigi);

 77:     /* Check convergence */
 78:     DSGetDimensions(eps->ds,NULL,NULL,NULL,NULL,&t);
 79: #if 0
 80:     /* take into account also left residual */
 81:     BVGetColumn(eps->V,nv,&u);
 82:     VecNorm(u,NORM_2,&beta1);
 83:     BVRestoreColumn(eps->V,nv,&u);
 84:     VecNorm(w,NORM_2,&beta2);  /* w contains B*V[nv] */
 85:     beta1 = PetscMax(beta1,beta2);
 86: #endif
 87:     EPSKrylovConvergence(eps,PETSC_FALSE,eps->nconv,t-eps->nconv,beta*beta1,0.0,1.0,&k);
 88:     (*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx);
 89:     nconv = k;

 91:     /* Update l */
 92:     if (eps->reason != EPS_CONVERGED_ITERATING || breakdown) l = 0;
 93:     else {
 94:       l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
 95:       l = PetscMin(l,t);
 96:       DSGetArrayReal(eps->ds,DS_MAT_T,&a);
 97:       if (*(a+ld+k+l-1)!=0) {
 98:         if (k+l<t-1) l = l+1;
 99:         else l = l-1;
100:       }
101:       DSRestoreArrayReal(eps->ds,DS_MAT_T,&a);
102:     }
103:     if (!ctx->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */

105:     if (eps->reason == EPS_CONVERGED_ITERATING) {
106:       if (breakdown) SETERRQ1(PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Breakdown in Indefinite Krylov-Schur (beta=%g)",beta);
107:       else {
108:         /* Prepare the Rayleigh quotient for restart */
109:         DSGetArray(eps->ds,DS_MAT_Q,&Q);
110:         DSGetArrayReal(eps->ds,DS_MAT_T,&a);
111:         DSGetArrayReal(eps->ds,DS_MAT_D,&omega);
112:         b = a + ld;
113:         r = a + 2*ld;
114:         for (i=k;i<k+l;i++) {
115:           r[i] = PetscRealPart(Q[nv-1+i*ld]*beta);
116:         }
117:         b[k+l-1] = r[k+l-1];
118:         omega[k+l] = omega[nv];
119:         DSRestoreArrayReal(eps->ds,DS_MAT_T,&a);
120:         DSRestoreArray(eps->ds,DS_MAT_Q,&Q);
121:         DSRestoreArrayReal(eps->ds,DS_MAT_D,&omega);
122:       }
123:     }
124:     /* Update the corresponding vectors V(:,idx) = V*Q(:,idx) */
125:     DSGetMat(eps->ds,DS_MAT_Q,&U);
126:     BVMultInPlace(eps->V,U,eps->nconv,k+l);
127:     MatDestroy(&U);

129:     /* Move restart vector and update signature */
130:     if (eps->reason == EPS_CONVERGED_ITERATING && !breakdown) {
131:       BVCopyColumn(eps->V,nv,k+l);
132:       DSGetArrayReal(eps->ds,DS_MAT_D,&omega);
133:       VecCreateSeq(PETSC_COMM_SELF,k+l,&vomega);
134:       VecGetArray(vomega,&aux);
135:       for (i=0;i<k+l;i++) aux[i] = omega[i];
136:       VecRestoreArray(vomega,&aux);
137:       BVSetActiveColumns(eps->V,0,k+l);
138:       BVSetSignature(eps->V,vomega);
139:       VecDestroy(&vomega);
140:       DSRestoreArrayReal(eps->ds,DS_MAT_D,&omega);
141:     }

143:     eps->nconv = k;
144:     EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,nv);
145:   }
146:   DSSetDimensions(eps->ds,eps->nconv,0,0,0);
147:   DSSetState(eps->ds,DS_STATE_RAW);
148:   return(0);
149: }