Actual source code: ks-symm.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 eigenproblems
 14: */

 16: #include <slepc/private/epsimpl.h>
 17: #include "krylovschur.h"

 19: PetscErrorCode EPSSolve_KrylovSchur_Symm(EPS eps)
 20: {
 21:   PetscErrorCode  ierr;
 22:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
 23:   PetscInt        k,l,ld,nv,nconv;
 24:   Mat             U;
 25:   PetscReal       *a,*b,beta;
 26:   PetscBool       breakdown;

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

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

 35:   /* Restart loop */
 36:   while (eps->reason == EPS_CONVERGED_ITERATING) {
 37:     eps->its++;

 39:     /* Compute an nv-step Lanczos factorization */
 40:     nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
 41:     DSGetArrayReal(eps->ds,DS_MAT_T,&a);
 42:     b = a + ld;
 43:     EPSFullLanczos(eps,a,b,eps->nconv+l,&nv,&breakdown);
 44:     beta = b[nv-1];
 45:     DSRestoreArrayReal(eps->ds,DS_MAT_T,&a);
 46:     DSSetDimensions(eps->ds,nv,0,eps->nconv,eps->nconv+l);
 47:     if (l==0) {
 48:       DSSetState(eps->ds,DS_STATE_INTERMEDIATE);
 49:     } else {
 50:       DSSetState(eps->ds,DS_STATE_RAW);
 51:     }
 52:     BVSetActiveColumns(eps->V,eps->nconv,nv);

 54:     /* Solve projected problem */
 55:     DSSolve(eps->ds,eps->eigr,NULL);
 56:     if (eps->arbitrary) { EPSGetArbitraryValues(eps,eps->rr,eps->ri); }
 57:     DSSort(eps->ds,eps->eigr,NULL,eps->rr,eps->ri,NULL);
 58:     DSUpdateExtraRow(eps->ds);
 59:     DSSynchronize(eps->ds,eps->eigr,NULL);

 61:     /* Check convergence */
 62:     EPSKrylovConvergence(eps,PETSC_FALSE,eps->nconv,nv-eps->nconv,beta,0.0,1.0,&k);
 63:     (*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx);
 64:     nconv = k;

 66:     /* Update l */
 67:     if (eps->reason != EPS_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
 68:     else l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
 69:     if (!ctx->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */

 71:     if (eps->reason == EPS_CONVERGED_ITERATING) {
 72:       if (breakdown || k==nv) {
 73:         /* Start a new Lanczos factorization */
 74:         PetscInfo2(eps,"Breakdown in Krylov-Schur method (it=%D norm=%g)\n",eps->its,(double)beta);
 75:         if (k<eps->nev) {
 76:           EPSGetStartVector(eps,k,&breakdown);
 77:           if (breakdown) {
 78:             eps->reason = EPS_DIVERGED_BREAKDOWN;
 79:             PetscInfo(eps,"Unable to generate more start vectors\n");
 80:           }
 81:         }
 82:       } else {
 83:         /* Prepare the Rayleigh quotient for restart */
 84:         DSTruncate(eps->ds,k+l);
 85:       }
 86:     }
 87:     /* Update the corresponding vectors V(:,idx) = V*Q(:,idx) */
 88:     DSGetMat(eps->ds,DS_MAT_Q,&U);
 89:     BVMultInPlace(eps->V,U,eps->nconv,k+l);
 90:     MatDestroy(&U);

 92:     /* Normalize u and append it to V */
 93:     if (eps->reason == EPS_CONVERGED_ITERATING && !breakdown) {
 94:       BVCopyColumn(eps->V,nv,k+l);
 95:     }

 97:     eps->nconv = k;
 98:     EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,nv);
 99:   }
100:   DSSetDimensions(eps->ds,eps->nconv,0,0,0);
101:   DSSetState(eps->ds,DS_STATE_RAW);
102:   return(0);
103: }