Actual source code: ks-indef.c

slepc-3.18.1 2022-11-02
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-, 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:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
 21:   PetscInt        k,l,ld,nv,t,nconv=0;
 22:   Mat             U,D;
 23:   Vec             vomega,w=eps->work[0];
 24:   PetscReal       *a,*b,beta,beta1=1.0,*omega;
 25:   PetscBool       breakdown=PETSC_FALSE,symmlost=PETSC_FALSE;

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

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

 32:   /* Extract sigma[0] from BV, computed during normalization */
 33:   DSSetDimensions(eps->ds,1,PETSC_DEFAULT,PETSC_DEFAULT);
 34:   BVSetActiveColumns(eps->V,0,1);
 35:   DSGetMatAndColumn(eps->ds,DS_MAT_D,0,&D,&vomega);
 36:   BVGetSignature(eps->V,vomega);
 37:   DSRestoreMatAndColumn(eps->ds,DS_MAT_D,0,&D,&vomega);
 38:   l = 0;

 40:   /* Restart loop */
 41:   while (eps->reason == EPS_CONVERGED_ITERATING) {
 42:     eps->its++;

 44:     /* Compute an nv-step Lanczos factorization */
 45:     nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
 46:     DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l);
 47:     DSGetArrayReal(eps->ds,DS_MAT_T,&a);
 48:     b = a + ld;
 49:     DSGetArrayReal(eps->ds,DS_MAT_D,&omega);
 50:     EPSPseudoLanczos(eps,a,b,omega,eps->nconv+l,&nv,&breakdown,&symmlost,NULL,w);
 51:     if (symmlost) {
 52:       eps->reason = EPS_DIVERGED_SYMMETRY_LOST;
 53:       if (nv==eps->nconv+l+1) { eps->nconv = nconv; break; }
 54:     }
 55:     beta = b[nv-1];
 56:     DSRestoreArrayReal(eps->ds,DS_MAT_T,&a);
 57:     DSRestoreArrayReal(eps->ds,DS_MAT_D,&omega);
 58:     DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l);
 59:     DSSetState(eps->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE);
 60:     BVSetActiveColumns(eps->V,eps->nconv,nv);

 62:     /* Solve projected problem */
 63:     DSSolve(eps->ds,eps->eigr,eps->eigi);
 64:     DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL);
 65:     DSUpdateExtraRow(eps->ds);
 66:     DSSynchronize(eps->ds,eps->eigr,eps->eigi);

 68:     /* Check convergence */
 69:     DSGetDimensions(eps->ds,NULL,NULL,NULL,&t);
 70: #if 0
 71:     /* take into account also left residual */
 72:     BVGetColumn(eps->V,nv,&u);
 73:     VecNorm(u,NORM_2,&beta1);
 74:     BVRestoreColumn(eps->V,nv,&u);
 75:     VecNorm(w,NORM_2,&beta2);  /* w contains B*V[nv] */
 76:     beta1 = PetscMax(beta1,beta2);
 77: #endif
 78:     EPSKrylovConvergence(eps,PETSC_FALSE,eps->nconv,t-eps->nconv,beta*beta1,0.0,1.0,&k);
 79:     (*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx);
 80:     nconv = k;

 82:     /* Update l */
 83:     if (eps->reason != EPS_CONVERGED_ITERATING || breakdown) l = 0;
 84:     else {
 85:       l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
 86:       l = PetscMin(l,t);
 87:       DSGetTruncateSize(eps->ds,k,t,&l);
 88:     }
 89:     if (!ctx->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
 90:     if (l) PetscInfo(eps,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l);

 92:     if (eps->reason == EPS_CONVERGED_ITERATING) {
 94:       /* Prepare the Rayleigh quotient for restart */
 95:       DSTruncate(eps->ds,k+l,PETSC_FALSE);
 96:     }
 97:     /* Update the corresponding vectors V(:,idx) = V*Q(:,idx) */
 98:     DSGetMat(eps->ds,DS_MAT_Q,&U);
 99:     BVMultInPlace(eps->V,U,eps->nconv,k+l);
100:     DSRestoreMat(eps->ds,DS_MAT_Q,&U);

102:     /* Move restart vector and update signature */
103:     if (eps->reason == EPS_CONVERGED_ITERATING && !breakdown) {
104:       BVCopyColumn(eps->V,nv,k+l);
105:       BVSetActiveColumns(eps->V,0,k+l);
106:       DSGetMatAndColumn(eps->ds,DS_MAT_D,0,&D,&vomega);
107:       BVSetSignature(eps->V,vomega);
108:       DSRestoreMatAndColumn(eps->ds,DS_MAT_D,0,&D,&vomega);
109:     }

111:     eps->nconv = k;
112:     EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,nv);
113:   }
114:   DSTruncate(eps->ds,eps->nconv,PETSC_TRUE);
115:   return 0;
116: }