Actual source code: slp-twosided.c
slepc-3.21.1 2024-04-26
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: Two-sided variant of the NEPSLP solver.
12: */
14: #include <slepc/private/nepimpl.h>
15: #include "slp.h"
17: typedef struct _n_nep_def_ctx *NEP_NEDEF_CTX;
19: struct _n_nep_def_ctx {
20: PetscInt n;
21: PetscBool ref;
22: PetscScalar *eig;
23: BV V,W;
24: };
26: typedef struct { /* context for two-sided solver */
27: Mat Ft;
28: Mat Jt;
29: Vec w;
30: } NEP_SLPTS_MATSHELL;
32: typedef struct { /* context for non-equivalence deflation */
33: NEP_NEDEF_CTX defctx;
34: Mat F;
35: Mat P;
36: Mat J;
37: KSP ksp;
38: PetscBool isJ;
39: PetscScalar lambda;
40: Vec w[2];
41: } NEP_NEDEF_MATSHELL;
43: static PetscErrorCode MatMult_SLPTS_Right(Mat M,Vec x,Vec y)
44: {
45: NEP_SLPTS_MATSHELL *ctx;
47: PetscFunctionBegin;
48: PetscCall(MatShellGetContext(M,&ctx));
49: PetscCall(MatMult(ctx->Jt,x,ctx->w));
50: PetscCall(MatSolve(ctx->Ft,ctx->w,y));
51: PetscFunctionReturn(PETSC_SUCCESS);
52: }
54: static PetscErrorCode MatMult_SLPTS_Left(Mat M,Vec x,Vec y)
55: {
56: NEP_SLPTS_MATSHELL *ctx;
58: PetscFunctionBegin;
59: PetscCall(MatShellGetContext(M,&ctx));
60: PetscCall(MatMultTranspose(ctx->Jt,x,ctx->w));
61: PetscCall(MatSolveTranspose(ctx->Ft,ctx->w,y));
62: PetscFunctionReturn(PETSC_SUCCESS);
63: }
65: static PetscErrorCode MatDestroy_SLPTS(Mat M)
66: {
67: NEP_SLPTS_MATSHELL *ctx;
69: PetscFunctionBegin;
70: PetscCall(MatShellGetContext(M,&ctx));
71: PetscCall(VecDestroy(&ctx->w));
72: PetscCall(PetscFree(ctx));
73: PetscFunctionReturn(PETSC_SUCCESS);
74: }
76: #if defined(PETSC_HAVE_CUDA)
77: static PetscErrorCode MatCreateVecs_SLPTS(Mat M,Vec *left,Vec *right)
78: {
79: NEP_SLPTS_MATSHELL *ctx;
81: PetscFunctionBegin;
82: PetscCall(MatShellGetContext(M,&ctx));
83: if (right) PetscCall(VecDuplicate(ctx->w,right));
84: if (left) PetscCall(VecDuplicate(ctx->w,left));
85: PetscFunctionReturn(PETSC_SUCCESS);
86: }
87: #endif
89: static PetscErrorCode NEPSLPSetUpEPSMat(NEP nep,Mat F,Mat J,PetscBool left,Mat *M)
90: {
91: Mat Mshell;
92: PetscInt nloc,mloc;
93: NEP_SLPTS_MATSHELL *shellctx;
95: PetscFunctionBegin;
96: /* Create mat shell */
97: PetscCall(PetscNew(&shellctx));
98: shellctx->Ft = F;
99: shellctx->Jt = J;
100: PetscCall(MatGetLocalSize(nep->function,&mloc,&nloc));
101: PetscCall(MatCreateShell(PetscObjectComm((PetscObject)nep),nloc,mloc,PETSC_DETERMINE,PETSC_DETERMINE,shellctx,&Mshell));
102: if (left) PetscCall(MatShellSetOperation(Mshell,MATOP_MULT,(void(*)(void))MatMult_SLPTS_Left));
103: else PetscCall(MatShellSetOperation(Mshell,MATOP_MULT,(void(*)(void))MatMult_SLPTS_Right));
104: PetscCall(MatShellSetOperation(Mshell,MATOP_DESTROY,(void(*)(void))MatDestroy_SLPTS));
105: #if defined(PETSC_HAVE_CUDA)
106: PetscCall(MatShellSetOperation(Mshell,MATOP_CREATE_VECS,(void(*)(void))MatCreateVecs_SLPTS));
107: #endif
108: *M = Mshell;
109: PetscCall(MatCreateVecs(nep->function,&shellctx->w,NULL));
110: PetscFunctionReturn(PETSC_SUCCESS);
111: }
113: /* Functions for deflation */
114: static PetscErrorCode NEPDeflationNEDestroy(NEP_NEDEF_CTX defctx)
115: {
116: PetscFunctionBegin;
117: if (!defctx) PetscFunctionReturn(PETSC_SUCCESS);
118: PetscCall(PetscFree(defctx->eig));
119: PetscCall(PetscFree(defctx));
120: PetscFunctionReturn(PETSC_SUCCESS);
121: }
123: static PetscErrorCode NEPDeflationNECreate(NEP nep,BV V,BV W,PetscInt sz,NEP_NEDEF_CTX *defctx)
124: {
125: NEP_NEDEF_CTX op;
127: PetscFunctionBegin;
128: PetscCall(PetscNew(&op));
129: *defctx = op;
130: op->n = 0;
131: op->ref = PETSC_FALSE;
132: PetscCall(PetscCalloc1(sz,&op->eig));
133: PetscCall(PetscObjectStateIncrease((PetscObject)V));
134: PetscCall(PetscObjectStateIncrease((PetscObject)W));
135: op->V = V;
136: op->W = W;
137: PetscFunctionReturn(PETSC_SUCCESS);
138: }
140: static PetscErrorCode NEPDeflationNEComputeFunction(NEP nep,Mat M,PetscScalar lambda)
141: {
142: NEP_NEDEF_MATSHELL *matctx;
144: PetscFunctionBegin;
145: PetscCall(MatShellGetContext(M,&matctx));
146: if (lambda==matctx->lambda) PetscFunctionReturn(PETSC_SUCCESS);
147: PetscCall(NEPComputeFunction(nep,lambda,matctx->F,matctx->P));
148: if (matctx->isJ) PetscCall(NEPComputeJacobian(nep,lambda,matctx->J));
149: if (matctx->ksp) PetscCall(NEP_KSPSetOperators(matctx->ksp,matctx->F,matctx->P));
150: matctx->lambda = lambda;
151: PetscFunctionReturn(PETSC_SUCCESS);
152: }
154: static PetscErrorCode MatMult_NEPDeflationNE(Mat M,Vec x,Vec r)
155: {
156: Vec t,tt;
157: PetscScalar *h,*alpha,lambda,*eig;
158: PetscInt i,k;
159: NEP_NEDEF_MATSHELL *matctx;
161: PetscFunctionBegin;
162: PetscCall(MatShellGetContext(M,&matctx));
163: if (matctx->defctx->n && !matctx->defctx->ref) {
164: k = matctx->defctx->n;
165: lambda = matctx->lambda;
166: eig = matctx->defctx->eig;
167: t = matctx->w[0];
168: PetscCall(VecCopy(x,t));
169: PetscCall(PetscMalloc2(k,&h,k,&alpha));
170: for (i=0;i<k;i++) alpha[i] = (lambda-eig[i]-1.0)/(lambda-eig[i]);
171: PetscCall(BVDotVec(matctx->defctx->V,t,h));
172: for (i=0;i<k;i++) h[i] *= alpha[i];
173: PetscCall(BVMultVec(matctx->defctx->W,-1.0,1.0,t,h));
174: PetscCall(MatMult(matctx->isJ?matctx->J:matctx->F,t,r));
175: if (matctx->isJ) {
176: for (i=0;i<k;i++) h[i] *= (1.0/((lambda-eig[i])*(lambda-eig[i])))/alpha[i];
177: tt = matctx->w[1];
178: PetscCall(BVMultVec(matctx->defctx->W,-1.0,0.0,tt,h));
179: PetscCall(MatMult(matctx->F,tt,t));
180: PetscCall(VecAXPY(r,1.0,t));
181: }
182: PetscCall(PetscFree2(h,alpha));
183: } else PetscCall(MatMult(matctx->isJ?matctx->J:matctx->F,x,r));
184: PetscFunctionReturn(PETSC_SUCCESS);
185: }
187: static PetscErrorCode MatMultTranspose_NEPDeflationNE(Mat M,Vec x,Vec r)
188: {
189: Vec t,tt;
190: PetscScalar *h,*alphaC,lambda,*eig;
191: PetscInt i,k;
192: NEP_NEDEF_MATSHELL *matctx;
194: PetscFunctionBegin;
195: PetscCall(MatShellGetContext(M,&matctx));
196: t = matctx->w[0];
197: PetscCall(VecCopy(x,t));
198: if (matctx->defctx->n && !matctx->defctx->ref) {
199: PetscCall(VecConjugate(t));
200: k = matctx->defctx->n;
201: lambda = matctx->lambda;
202: eig = matctx->defctx->eig;
203: PetscCall(PetscMalloc2(k,&h,k,&alphaC));
204: for (i=0;i<k;i++) alphaC[i] = PetscConj((lambda-eig[i]-1.0)/(lambda-eig[i]));
205: PetscCall(BVDotVec(matctx->defctx->W,t,h));
206: for (i=0;i<k;i++) h[i] *= alphaC[i];
207: PetscCall(BVMultVec(matctx->defctx->V,-1.0,1.0,t,h));
208: PetscCall(VecConjugate(t));
209: PetscCall(MatMultTranspose(matctx->isJ?matctx->J:matctx->F,t,r));
210: if (matctx->isJ) {
211: for (i=0;i<k;i++) h[i] *= PetscConj(1.0/((lambda-eig[i])*(lambda-eig[i])))/alphaC[i];
212: tt = matctx->w[1];
213: PetscCall(BVMultVec(matctx->defctx->V,-1.0,0.0,tt,h));
214: PetscCall(VecConjugate(tt));
215: PetscCall(MatMultTranspose(matctx->F,tt,t));
216: PetscCall(VecAXPY(r,1.0,t));
217: }
218: PetscCall(PetscFree2(h,alphaC));
219: } else PetscCall(MatMultTranspose(matctx->isJ?matctx->J:matctx->F,t,r));
220: PetscFunctionReturn(PETSC_SUCCESS);
221: }
223: static PetscErrorCode MatSolve_NEPDeflationNE(Mat M,Vec b,Vec x)
224: {
225: PetscScalar *h,*alpha,lambda,*eig;
226: PetscInt i,k;
227: NEP_NEDEF_MATSHELL *matctx;
229: PetscFunctionBegin;
230: PetscCall(MatShellGetContext(M,&matctx));
231: if (!matctx->ksp) {
232: PetscCall(VecCopy(b,x));
233: PetscFunctionReturn(PETSC_SUCCESS);
234: }
235: PetscCall(KSPSolve(matctx->ksp,b,x));
236: if (matctx->defctx->n && !matctx->defctx->ref) {
237: k = matctx->defctx->n;
238: lambda = matctx->lambda;
239: eig = matctx->defctx->eig;
240: PetscCall(PetscMalloc2(k,&h,k,&alpha));
241: PetscCall(BVDotVec(matctx->defctx->V,x,h));
242: for (i=0;i<k;i++) alpha[i] = (lambda-eig[i]-1.0)/(lambda-eig[i]);
243: for (i=0;i<k;i++) h[i] *= alpha[i]/(1.0-alpha[i]);
244: PetscCall(BVMultVec(matctx->defctx->W,1.0,1.0,x,h));
245: PetscCall(PetscFree2(h,alpha));
246: }
247: PetscFunctionReturn(PETSC_SUCCESS);
248: }
250: static PetscErrorCode MatSolveTranspose_NEPDeflationNE(Mat M,Vec b,Vec x)
251: {
252: PetscScalar *h,*alphaC,lambda,*eig;
253: PetscInt i,k;
254: NEP_NEDEF_MATSHELL *matctx;
256: PetscFunctionBegin;
257: PetscCall(MatShellGetContext(M,&matctx));
258: if (!matctx->ksp) {
259: PetscCall(VecCopy(b,x));
260: PetscFunctionReturn(PETSC_SUCCESS);
261: }
262: PetscCall(KSPSolveTranspose(matctx->ksp,b,x));
263: if (matctx->defctx->n && !matctx->defctx->ref) {
264: PetscCall(VecConjugate(x));
265: k = matctx->defctx->n;
266: lambda = matctx->lambda;
267: eig = matctx->defctx->eig;
268: PetscCall(PetscMalloc2(k,&h,k,&alphaC));
269: PetscCall(BVDotVec(matctx->defctx->W,x,h));
270: for (i=0;i<k;i++) alphaC[i] = PetscConj((lambda-eig[i]-1.0)/(lambda-eig[i]));
271: for (i=0;i<k;i++) h[i] *= alphaC[i]/(1.0-alphaC[i]);
272: PetscCall(BVMultVec(matctx->defctx->V,1.0,1.0,x,h));
273: PetscCall(PetscFree2(h,alphaC));
274: PetscCall(VecConjugate(x));
275: }
276: PetscFunctionReturn(PETSC_SUCCESS);
277: }
279: static PetscErrorCode MatDestroy_NEPDeflationNE(Mat M)
280: {
281: NEP_NEDEF_MATSHELL *matctx;
283: PetscFunctionBegin;
284: PetscCall(MatShellGetContext(M,&matctx));
285: PetscCall(VecDestroy(&matctx->w[0]));
286: PetscCall(VecDestroy(&matctx->w[1]));
287: PetscCall(PetscFree(matctx));
288: PetscFunctionReturn(PETSC_SUCCESS);
289: }
291: static PetscErrorCode MatCreateVecs_NEPDeflationNE(Mat M,Vec *right,Vec *left)
292: {
293: NEP_NEDEF_MATSHELL *matctx;
295: PetscFunctionBegin;
296: PetscCall(MatShellGetContext(M,&matctx));
297: PetscCall(MatCreateVecs(matctx->F,right,left));
298: PetscFunctionReturn(PETSC_SUCCESS);
299: }
301: static PetscErrorCode NEPDeflationNEFunctionCreate(NEP_NEDEF_CTX defctx,NEP nep,Mat F,Mat P,Mat J,KSP ksp,PetscBool isJ,Mat *Mshell)
302: {
303: NEP_NEDEF_MATSHELL *matctx;
304: PetscInt nloc,mloc;
306: PetscFunctionBegin;
307: /* Create mat shell */
308: PetscCall(PetscNew(&matctx));
309: PetscCall(MatGetLocalSize(nep->function,&mloc,&nloc));
310: PetscCall(MatCreateShell(PetscObjectComm((PetscObject)nep),nloc,mloc,PETSC_DETERMINE,PETSC_DETERMINE,matctx,Mshell));
311: matctx->F = F;
312: matctx->P = P;
313: matctx->J = J;
314: matctx->isJ = isJ;
315: matctx->ksp = ksp;
316: matctx->defctx = defctx;
317: matctx->lambda = PETSC_MAX_REAL;
318: PetscCall(MatCreateVecs(F,&matctx->w[0],NULL));
319: PetscCall(VecDuplicate(matctx->w[0],&matctx->w[1]));
320: PetscCall(MatShellSetOperation(*Mshell,MATOP_MULT,(void(*)(void))MatMult_NEPDeflationNE));
321: PetscCall(MatShellSetOperation(*Mshell,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_NEPDeflationNE));
322: PetscCall(MatShellSetOperation(*Mshell,MATOP_SOLVE,(void(*)(void))MatSolve_NEPDeflationNE));
323: PetscCall(MatShellSetOperation(*Mshell,MATOP_SOLVE_TRANSPOSE,(void(*)(void))MatSolveTranspose_NEPDeflationNE));
324: PetscCall(MatShellSetOperation(*Mshell,MATOP_DESTROY,(void(*)(void))MatDestroy_NEPDeflationNE));
325: PetscCall(MatShellSetOperation(*Mshell,MATOP_CREATE_VECS,(void(*)(void))MatCreateVecs_NEPDeflationNE));
326: PetscFunctionReturn(PETSC_SUCCESS);
327: }
329: static PetscErrorCode NEPDeflationNERecoverEigenvectors(NEP_NEDEF_CTX defctx,Vec u,Vec w,PetscScalar lambda)
330: {
331: PetscScalar *h,*alpha,*eig;
332: PetscInt i,k;
334: PetscFunctionBegin;
335: if (w) PetscCall(VecConjugate(w));
336: if (defctx->n && !defctx->ref) {
337: eig = defctx->eig;
338: k = defctx->n;
339: PetscCall(PetscMalloc2(k,&h,k,&alpha));
340: for (i=0;i<k;i++) alpha[i] = (lambda-eig[i]-1.0)/(lambda-eig[i]);
341: PetscCall(BVDotVec(defctx->V,u,h));
342: for (i=0;i<k;i++) h[i] *= alpha[i];
343: PetscCall(BVMultVec(defctx->W,-1.0,1.0,u,h));
344: PetscCall(VecNormalize(u,NULL));
345: if (w) {
346: PetscCall(BVDotVec(defctx->W,w,h));
347: for (i=0;i<k;i++) alpha[i] = PetscConj((lambda-eig[i]-1.0)/(lambda-eig[i]));
348: for (i=0;i<k;i++) h[i] *= alpha[i];
349: PetscCall(BVMultVec(defctx->V,-1.0,1.0,w,h));
350: PetscCall(VecNormalize(w,NULL));
351: }
352: PetscCall(PetscFree2(h,alpha));
353: }
354: PetscFunctionReturn(PETSC_SUCCESS);
355: }
357: static PetscErrorCode NEPDeflationNELocking(NEP_NEDEF_CTX defctx,Vec u,Vec w,PetscScalar lambda)
358: {
359: PetscInt n;
361: PetscFunctionBegin;
362: n = defctx->n++;
363: defctx->eig[n] = lambda;
364: PetscCall(BVInsertVec(defctx->V,n,u));
365: PetscCall(BVInsertVec(defctx->W,n,w));
366: PetscCall(BVSetActiveColumns(defctx->V,0,defctx->n));
367: PetscCall(BVSetActiveColumns(defctx->W,0,defctx->n));
368: PetscCall(BVBiorthonormalizeColumn(defctx->V,defctx->W,n,NULL));
369: PetscFunctionReturn(PETSC_SUCCESS);
370: }
372: static PetscErrorCode NEPDeflationNESetRefine(NEP_NEDEF_CTX defctx,PetscBool ref)
373: {
374: PetscFunctionBegin;
375: defctx->ref = ref;
376: PetscFunctionReturn(PETSC_SUCCESS);
377: }
379: PetscErrorCode NEPSolve_SLP_Twosided(NEP nep)
380: {
381: NEP_SLP *ctx = (NEP_SLP*)nep->data;
382: Mat mF,mJ,M,Mt;
383: Vec u,r,t,w;
384: BV X,Y;
385: PetscScalar sigma,lambda,mu,im=0.0,mu2,im2;
386: PetscReal resnorm,resl;
387: PetscInt nconv,nconv2,i;
388: PetscBool skip=PETSC_FALSE,lock=PETSC_FALSE;
389: NEP_NEDEF_CTX defctx=NULL; /* Extended operator for deflation */
391: PetscFunctionBegin;
392: /* get initial approximation of eigenvalue and eigenvector */
393: PetscCall(NEPGetDefaultShift(nep,&sigma));
394: if (!nep->nini) PetscCall(BVSetRandomColumn(nep->V,0));
395: PetscCall(BVSetRandomColumn(nep->W,0));
396: lambda = sigma;
397: if (!ctx->ksp) PetscCall(NEPSLPGetKSP(nep,&ctx->ksp));
398: PetscCall(BVDuplicate(nep->V,&X));
399: PetscCall(BVDuplicate(nep->W,&Y));
400: PetscCall(NEPDeflationNECreate(nep,X,Y,nep->nev,&defctx));
401: PetscCall(BVGetColumn(nep->V,0,&t));
402: PetscCall(VecDuplicate(t,&u));
403: PetscCall(VecDuplicate(t,&w));
404: PetscCall(BVRestoreColumn(nep->V,0,&t));
405: PetscCall(BVCopyVec(nep->V,0,u));
406: PetscCall(BVCopyVec(nep->W,0,w));
407: PetscCall(VecDuplicate(u,&r));
408: PetscCall(NEPDeflationNEFunctionCreate(defctx,nep,nep->function,nep->function_pre?nep->function_pre:nep->function,NULL,ctx->ksp,PETSC_FALSE,&mF));
409: PetscCall(NEPDeflationNEFunctionCreate(defctx,nep,nep->function,nep->function,nep->jacobian,NULL,PETSC_TRUE,&mJ));
410: PetscCall(NEPSLPSetUpEPSMat(nep,mF,mJ,PETSC_FALSE,&M));
411: PetscCall(NEPSLPSetUpEPSMat(nep,mF,mJ,PETSC_TRUE,&Mt));
412: PetscCall(EPSSetOperators(ctx->eps,M,NULL));
413: PetscCall(MatDestroy(&M));
414: PetscCall(EPSSetOperators(ctx->epsts,Mt,NULL));
415: PetscCall(MatDestroy(&Mt));
417: /* Restart loop */
418: while (nep->reason == NEP_CONVERGED_ITERATING) {
419: nep->its++;
421: /* form residual, r = T(lambda)*u (used in convergence test only) */
422: PetscCall(NEPDeflationNEComputeFunction(nep,mF,lambda));
423: PetscCall(MatMultTranspose(mF,w,r));
424: PetscCall(VecNorm(r,NORM_2,&resl));
425: PetscCall(MatMult(mF,u,r));
427: /* convergence test */
428: PetscCall(VecNorm(r,NORM_2,&resnorm));
429: resnorm = PetscMax(resnorm,resl);
430: PetscCall((*nep->converged)(nep,lambda,0,resnorm,&nep->errest[nep->nconv],nep->convergedctx));
431: nep->eigr[nep->nconv] = lambda;
432: if (nep->errest[nep->nconv]<=nep->tol || nep->errest[nep->nconv]<=ctx->deftol) {
433: if (nep->errest[nep->nconv]<=ctx->deftol && !defctx->ref && nep->nconv) {
434: PetscCall(NEPDeflationNERecoverEigenvectors(defctx,u,w,lambda));
435: PetscCall(VecConjugate(w));
436: PetscCall(NEPDeflationNESetRefine(defctx,PETSC_TRUE));
437: PetscCall(MatMultTranspose(mF,w,r));
438: PetscCall(VecNorm(r,NORM_2,&resl));
439: PetscCall(MatMult(mF,u,r));
440: PetscCall(VecNorm(r,NORM_2,&resnorm));
441: resnorm = PetscMax(resnorm,resl);
442: PetscCall((*nep->converged)(nep,lambda,0,resnorm,&nep->errest[nep->nconv],nep->convergedctx));
443: if (nep->errest[nep->nconv]<=nep->tol) lock = PETSC_TRUE;
444: } else if (nep->errest[nep->nconv]<=nep->tol) lock = PETSC_TRUE;
445: }
446: if (lock) {
447: lock = PETSC_FALSE;
448: skip = PETSC_TRUE;
449: PetscCall(NEPDeflationNERecoverEigenvectors(defctx,u,w,lambda));
450: PetscCall(NEPDeflationNELocking(defctx,u,w,lambda));
451: PetscCall(NEPDeflationNESetRefine(defctx,PETSC_FALSE));
452: PetscCall(BVInsertVec(nep->V,nep->nconv,u));
453: PetscCall(BVInsertVec(nep->W,nep->nconv,w));
454: PetscCall(VecConjugate(w));
455: nep->nconv = nep->nconv + 1;
456: }
457: PetscCall((*nep->stopping)(nep,nep->its,nep->max_it,nep->nconv,nep->nev,&nep->reason,nep->stoppingctx));
458: if (!skip || nep->reason>0) PetscCall(NEPMonitor(nep,nep->its,nep->nconv,nep->eigr,nep->eigi,nep->errest,(nep->reason>0)?nep->nconv:nep->nconv+1));
460: if (nep->reason == NEP_CONVERGED_ITERATING) {
461: if (!skip) {
462: /* evaluate T(lambda) and T'(lambda) */
463: PetscCall(NEPDeflationNEComputeFunction(nep,mF,lambda));
464: PetscCall(NEPDeflationNEComputeFunction(nep,mJ,lambda));
465: PetscCall(EPSSetInitialSpace(ctx->eps,1,&u));
466: PetscCall(EPSSetInitialSpace(ctx->epsts,1,&w));
468: /* compute new eigenvalue correction mu and eigenvector approximation u */
469: PetscCall(EPSSolve(ctx->eps));
470: PetscCall(EPSSolve(ctx->epsts));
471: PetscCall(EPSGetConverged(ctx->eps,&nconv));
472: PetscCall(EPSGetConverged(ctx->epsts,&nconv2));
473: if (!nconv||!nconv2) {
474: PetscCall(PetscInfo(nep,"iter=%" PetscInt_FMT ", inner iteration failed, stopping solve\n",nep->its));
475: nep->reason = NEP_DIVERGED_LINEAR_SOLVE;
476: break;
477: }
478: PetscCall(EPSGetEigenpair(ctx->eps,0,&mu,&im,u,NULL));
479: for (i=0;i<nconv2;i++) {
480: PetscCall(EPSGetEigenpair(ctx->epsts,i,&mu2,&im2,w,NULL));
481: if (SlepcAbsEigenvalue(mu-mu2,im-im2)/SlepcAbsEigenvalue(mu,im)<nep->tol*1000) break;
482: }
483: if (i==nconv2) {
484: PetscCall(PetscInfo(nep,"iter=%" PetscInt_FMT ", inner iteration failed, stopping solve\n",nep->its));
485: nep->reason = NEP_DIVERGED_LINEAR_SOLVE;
486: break;
487: }
489: mu = 1.0/mu;
490: PetscCheck(PetscAbsScalar(im)<PETSC_MACHINE_EPSILON,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Complex eigenvalue approximation - not implemented in real scalars");
491: } else {
492: nep->its--; /* do not count this as a full iteration */
493: /* use second eigenpair computed in previous iteration */
494: PetscCall(EPSGetConverged(ctx->eps,&nconv));
495: if (nconv>=2 && nconv2>=2) {
496: PetscCall(EPSGetEigenpair(ctx->eps,1,&mu,&im,u,NULL));
497: PetscCall(EPSGetEigenpair(ctx->epsts,1,&mu2,&im2,w,NULL));
498: mu = 1.0/mu;
499: } else {
500: PetscCall(BVSetRandomColumn(nep->V,nep->nconv));
501: PetscCall(BVSetRandomColumn(nep->W,nep->nconv));
502: PetscCall(BVCopyVec(nep->V,nep->nconv,u));
503: PetscCall(BVCopyVec(nep->W,nep->nconv,w));
504: mu = lambda-sigma;
505: }
506: skip = PETSC_FALSE;
507: }
508: /* correct eigenvalue */
509: lambda = lambda - mu;
510: }
511: }
512: PetscCall(VecDestroy(&u));
513: PetscCall(VecDestroy(&w));
514: PetscCall(VecDestroy(&r));
515: PetscCall(MatDestroy(&mF));
516: PetscCall(MatDestroy(&mJ));
517: PetscCall(BVDestroy(&X));
518: PetscCall(BVDestroy(&Y));
519: PetscCall(NEPDeflationNEDestroy(defctx));
520: PetscFunctionReturn(PETSC_SUCCESS);
521: }