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