Actual source code: ks-indef.c
slepc-3.11.2 2019-07-30
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: }