Actual source code: ks-symm.c
slepc-3.13.4 2020-09-02
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2020, 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,Op;
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: STGetOperator(eps->st,&Op);
44: BVMatLanczos(eps->V,Op,a,b,eps->nconv+l,&nv,&breakdown);
45: STRestoreOperator(eps->st,&Op);
46: beta = b[nv-1];
47: DSRestoreArrayReal(eps->ds,DS_MAT_T,&a);
48: DSSetDimensions(eps->ds,nv,0,eps->nconv,eps->nconv+l);
49: if (l==0) {
50: DSSetState(eps->ds,DS_STATE_INTERMEDIATE);
51: } else {
52: DSSetState(eps->ds,DS_STATE_RAW);
53: }
54: BVSetActiveColumns(eps->V,eps->nconv,nv);
56: /* Solve projected problem */
57: DSSolve(eps->ds,eps->eigr,NULL);
58: if (eps->arbitrary) { EPSGetArbitraryValues(eps,eps->rr,eps->ri); }
59: DSSort(eps->ds,eps->eigr,NULL,eps->rr,eps->ri,NULL);
60: DSUpdateExtraRow(eps->ds);
61: DSSynchronize(eps->ds,eps->eigr,NULL);
63: /* Check convergence */
64: EPSKrylovConvergence(eps,PETSC_FALSE,eps->nconv,nv-eps->nconv,beta,0.0,1.0,&k);
65: (*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx);
66: nconv = k;
68: /* Update l */
69: if (eps->reason != EPS_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
70: else l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
71: if (!ctx->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
73: if (eps->reason == EPS_CONVERGED_ITERATING) {
74: if (breakdown || k==nv) {
75: /* Start a new Lanczos factorization */
76: PetscInfo2(eps,"Breakdown in Krylov-Schur method (it=%D norm=%g)\n",eps->its,(double)beta);
77: if (k<eps->nev) {
78: EPSGetStartVector(eps,k,&breakdown);
79: if (breakdown) {
80: eps->reason = EPS_DIVERGED_BREAKDOWN;
81: PetscInfo(eps,"Unable to generate more start vectors\n");
82: }
83: }
84: } else {
85: /* Prepare the Rayleigh quotient for restart */
86: DSTruncate(eps->ds,k+l);
87: }
88: }
89: /* Update the corresponding vectors V(:,idx) = V*Q(:,idx) */
90: DSGetMat(eps->ds,DS_MAT_Q,&U);
91: BVMultInPlace(eps->V,U,eps->nconv,k+l);
92: MatDestroy(&U);
94: /* Normalize u and append it to V */
95: if (eps->reason == EPS_CONVERGED_ITERATING && !breakdown) {
96: BVCopyColumn(eps->V,nv,k+l);
97: }
99: eps->nconv = k;
100: EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,nv);
101: }
102: DSSetDimensions(eps->ds,eps->nconv,0,0,0);
103: DSSetState(eps->ds,DS_STATE_RAW);
104: return(0);
105: }