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 : SLEPc nonlinear eigensolver: "nleigs"
12 :
13 : Method: NLEIGS
14 :
15 : Algorithm:
16 :
17 : Fully rational Krylov method for nonlinear eigenvalue problems.
18 :
19 : References:
20 :
21 : [1] S. Guttel et al., "NLEIGS: A class of robust fully rational Krylov
22 : method for nonlinear eigenvalue problems", SIAM J. Sci. Comput.
23 : 36(6):A2842-A2864, 2014.
24 : */
25 :
26 : #include <slepc/private/nepimpl.h> /*I "slepcnep.h" I*/
27 : #include <slepcblaslapack.h>
28 : #include "nleigs.h"
29 :
30 11380 : PetscErrorCode NEPNLEIGSBackTransform(PetscObject ob,PetscInt n,PetscScalar *valr,PetscScalar *vali)
31 : {
32 11380 : NEP nep;
33 11380 : PetscInt j;
34 : #if !defined(PETSC_USE_COMPLEX)
35 11380 : PetscScalar t;
36 : #endif
37 :
38 11380 : PetscFunctionBegin;
39 11380 : nep = (NEP)ob;
40 : #if !defined(PETSC_USE_COMPLEX)
41 35318 : for (j=0;j<n;j++) {
42 23938 : if (vali[j] == 0) valr[j] = 1.0 / valr[j] + nep->target;
43 : else {
44 5631 : t = valr[j] * valr[j] + vali[j] * vali[j];
45 5631 : valr[j] = valr[j] / t + nep->target;
46 5631 : vali[j] = - vali[j] / t;
47 : }
48 : }
49 : #else
50 : for (j=0;j<n;j++) {
51 : valr[j] = 1.0 / valr[j] + nep->target;
52 : }
53 : #endif
54 11380 : PetscFunctionReturn(PETSC_SUCCESS);
55 : }
56 :
57 : /* Computes the roots of a polynomial */
58 3 : static PetscErrorCode NEPNLEIGSAuxiliarPRootFinder(PetscInt deg,PetscScalar *polcoeffs,PetscScalar *wr,PetscScalar *wi,PetscBool *avail)
59 : {
60 3 : PetscScalar *C;
61 3 : PetscBLASInt n_,lwork;
62 3 : PetscInt i;
63 : #if defined(PETSC_USE_COMPLEX)
64 : PetscReal *rwork=NULL;
65 : #endif
66 3 : PetscScalar *work;
67 3 : PetscBLASInt info;
68 :
69 3 : PetscFunctionBegin;
70 3 : *avail = PETSC_TRUE;
71 3 : if (deg>0) {
72 3 : PetscCall(PetscCalloc1(deg*deg,&C));
73 3 : PetscCall(PetscBLASIntCast(deg,&n_));
74 3 : for (i=0;i<deg-1;i++) {
75 0 : C[(deg+1)*i+1] = 1.0;
76 0 : C[(deg-1)*deg+i] = -polcoeffs[deg-i]/polcoeffs[0];
77 : }
78 3 : C[deg*deg+-1] = -polcoeffs[1]/polcoeffs[0];
79 3 : PetscCall(PetscBLASIntCast(3*deg,&lwork));
80 :
81 3 : PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
82 : #if !defined(PETSC_USE_COMPLEX)
83 3 : PetscCall(PetscMalloc1(lwork,&work));
84 3 : PetscCallBLAS("LAPACKgeev",LAPACKgeev_("N","N",&n_,C,&n_,wr,wi,NULL,&n_,NULL,&n_,work,&lwork,&info));
85 3 : if (info) *avail = PETSC_FALSE;
86 3 : PetscCall(PetscFree(work));
87 : #else
88 : PetscCall(PetscMalloc2(2*deg,&rwork,lwork,&work));
89 : PetscCallBLAS("LAPACKgeev",LAPACKgeev_("N","N",&n_,C,&n_,wr,NULL,&n_,NULL,&n_,work,&lwork,rwork,&info));
90 : if (info) *avail = PETSC_FALSE;
91 : PetscCall(PetscFree2(rwork,work));
92 : #endif
93 3 : PetscCall(PetscFPTrapPop());
94 3 : PetscCall(PetscFree(C));
95 : }
96 3 : PetscFunctionReturn(PETSC_SUCCESS);
97 : }
98 :
99 13 : static PetscErrorCode NEPNLEIGSAuxiliarRmDuplicates(PetscInt nin,PetscScalar *pin,PetscInt *nout,PetscScalar *pout,PetscInt max)
100 : {
101 13 : PetscInt i,j;
102 :
103 13 : PetscFunctionBegin;
104 29 : for (i=0;i<nin;i++) {
105 16 : if (max && *nout>=max) break;
106 16 : pout[(*nout)++] = pin[i];
107 20 : for (j=0;j<*nout-1;j++)
108 4 : if (PetscAbsScalar(pin[i]-pout[j])<PETSC_MACHINE_EPSILON*100) {
109 0 : (*nout)--;
110 0 : break;
111 : }
112 : }
113 13 : PetscFunctionReturn(PETSC_SUCCESS);
114 : }
115 :
116 11 : static PetscErrorCode NEPNLEIGSFNSingularities(FN f,PetscInt *nisol,PetscScalar **isol,PetscBool *rational)
117 : {
118 11 : FNCombineType ctype;
119 11 : FN f1,f2;
120 11 : PetscInt i,nq,nisol1,nisol2;
121 11 : PetscScalar *qcoeff,*wr,*wi,*isol1,*isol2;
122 11 : PetscBool flg,avail,rat1,rat2;
123 :
124 11 : PetscFunctionBegin;
125 11 : *rational = PETSC_FALSE;
126 11 : PetscCall(PetscObjectTypeCompare((PetscObject)f,FNRATIONAL,&flg));
127 11 : if (flg) {
128 10 : *rational = PETSC_TRUE;
129 10 : PetscCall(FNRationalGetDenominator(f,&nq,&qcoeff));
130 10 : if (nq>1) {
131 3 : PetscCall(PetscMalloc2(nq-1,&wr,nq-1,&wi));
132 3 : PetscCall(NEPNLEIGSAuxiliarPRootFinder(nq-1,qcoeff,wr,wi,&avail));
133 3 : if (avail) {
134 3 : PetscCall(PetscCalloc1(nq-1,isol));
135 3 : *nisol = 0;
136 6 : for (i=0;i<nq-1;i++)
137 : #if !defined(PETSC_USE_COMPLEX)
138 3 : if (wi[i]==0)
139 : #endif
140 3 : (*isol)[(*nisol)++] = wr[i];
141 3 : nq = *nisol; *nisol = 0;
142 6 : for (i=0;i<nq;i++) wr[i] = (*isol)[i];
143 3 : PetscCall(NEPNLEIGSAuxiliarRmDuplicates(nq,wr,nisol,*isol,0));
144 3 : PetscCall(PetscFree2(wr,wi));
145 0 : } else { *nisol=0; *isol = NULL; }
146 7 : } else { *nisol = 0; *isol = NULL; }
147 10 : PetscCall(PetscFree(qcoeff));
148 : }
149 11 : PetscCall(PetscObjectTypeCompare((PetscObject)f,FNCOMBINE,&flg));
150 11 : if (flg) {
151 1 : PetscCall(FNCombineGetChildren(f,&ctype,&f1,&f2));
152 1 : if (ctype != FN_COMBINE_COMPOSE && ctype != FN_COMBINE_DIVIDE) {
153 1 : PetscCall(NEPNLEIGSFNSingularities(f1,&nisol1,&isol1,&rat1));
154 1 : PetscCall(NEPNLEIGSFNSingularities(f2,&nisol2,&isol2,&rat2));
155 1 : if (nisol1+nisol2>0) {
156 1 : PetscCall(PetscCalloc1(nisol1+nisol2,isol));
157 1 : *nisol = 0;
158 1 : PetscCall(NEPNLEIGSAuxiliarRmDuplicates(nisol1,isol1,nisol,*isol,0));
159 1 : PetscCall(NEPNLEIGSAuxiliarRmDuplicates(nisol2,isol2,nisol,*isol,0));
160 : }
161 1 : *rational = (rat1&&rat2)?PETSC_TRUE:PETSC_FALSE;
162 1 : PetscCall(PetscFree(isol1));
163 1 : PetscCall(PetscFree(isol2));
164 : }
165 : }
166 11 : PetscFunctionReturn(PETSC_SUCCESS);
167 : }
168 :
169 3 : static PetscErrorCode NEPNLEIGSRationalSingularities(NEP nep,PetscInt *ndptx,PetscScalar *dxi,PetscBool *rational)
170 : {
171 3 : PetscInt nt,i,nisol;
172 3 : FN f;
173 3 : PetscScalar *isol;
174 3 : PetscBool rat;
175 :
176 3 : PetscFunctionBegin;
177 3 : *rational = PETSC_TRUE;
178 3 : *ndptx = 0;
179 3 : PetscCall(NEPGetSplitOperatorInfo(nep,&nt,NULL));
180 12 : for (i=0;i<nt;i++) {
181 9 : PetscCall(NEPGetSplitOperatorTerm(nep,i,NULL,&f));
182 9 : PetscCall(NEPNLEIGSFNSingularities(f,&nisol,&isol,&rat));
183 9 : if (nisol) {
184 3 : PetscCall(NEPNLEIGSAuxiliarRmDuplicates(nisol,isol,ndptx,dxi,0));
185 3 : PetscCall(PetscFree(isol));
186 : }
187 9 : *rational = ((*rational)&&rat)?PETSC_TRUE:PETSC_FALSE;
188 : }
189 3 : PetscFunctionReturn(PETSC_SUCCESS);
190 : }
191 :
192 : #if defined(SLEPC_MISSING_LAPACK_GGEV3)
193 : #define LAPGEEV "ggev"
194 : #else
195 : #define LAPGEEV "ggev3"
196 : #endif
197 :
198 : /* Adaptive Anderson-Antoulas algorithm */
199 22 : static PetscErrorCode NEPNLEIGSAAAComputation(NEP nep,PetscInt ndpt,PetscScalar *ds,PetscScalar *F,PetscInt *ndptx,PetscScalar *dxi)
200 : {
201 22 : NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
202 22 : PetscScalar mean=0.0,*z,*f,*C,*A,*VT,*work,*ww,szero=0.0,sone=1.0;
203 22 : PetscScalar *N,*D;
204 22 : PetscReal *S,norm,err,*R;
205 22 : PetscInt i,k,j,idx=0,cont;
206 22 : PetscBLASInt n_,m_,lda_,lwork,info,one=1;
207 : #if defined(PETSC_USE_COMPLEX)
208 : PetscReal *rwork;
209 : #endif
210 :
211 22 : PetscFunctionBegin;
212 22 : PetscCall(PetscBLASIntCast(8*ndpt,&lwork));
213 22 : PetscCall(PetscMalloc5(ndpt,&R,ndpt,&z,ndpt,&f,ndpt*ndpt,&C,ndpt,&ww));
214 22 : PetscCall(PetscMalloc6(ndpt*ndpt,&A,ndpt,&S,ndpt*ndpt,&VT,lwork,&work,ndpt,&D,ndpt,&N));
215 : #if defined(PETSC_USE_COMPLEX)
216 : PetscCall(PetscMalloc1(8*ndpt,&rwork));
217 : #endif
218 22 : PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
219 : norm = 0.0;
220 2222 : for (i=0;i<ndpt;i++) {
221 2200 : mean += F[i];
222 3509 : norm = PetscMax(PetscAbsScalar(F[i]),norm);
223 : }
224 22 : mean /= ndpt;
225 22 : PetscCall(PetscBLASIntCast(ndpt,&lda_));
226 2222 : for (i=0;i<ndpt;i++) R[i] = PetscAbsScalar(F[i]-mean);
227 : /* next support point */
228 : err = 0.0;
229 2222 : for (i=0;i<ndpt;i++) if (R[i]>=err) {idx = i; err = R[i];}
230 48 : for (k=0;k<ndpt-1;k++) {
231 48 : z[k] = ds[idx]; f[k] = F[idx]; R[idx] = -1.0;
232 : /* next column of Cauchy matrix */
233 4848 : for (i=0;i<ndpt;i++) {
234 4800 : C[i+k*ndpt] = 1.0/(ds[i]-ds[idx]);
235 : }
236 :
237 48 : PetscCall(PetscArrayzero(A,ndpt*ndpt));
238 : cont = 0;
239 4848 : for (i=0;i<ndpt;i++) {
240 4800 : if (R[i]!=-1.0) {
241 12846 : for (j=0;j<=k;j++)A[cont+j*ndpt] = C[i+j*ndpt]*F[i]-C[i+j*ndpt]*f[j];
242 4717 : cont++;
243 : }
244 : }
245 48 : PetscCall(PetscBLASIntCast(cont,&m_));
246 48 : PetscCall(PetscBLASIntCast(k+1,&n_));
247 : #if defined(PETSC_USE_COMPLEX)
248 : PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","A",&m_,&n_,A,&lda_,S,NULL,&lda_,VT,&lda_,work,&lwork,rwork,&info));
249 : #else
250 48 : PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","A",&m_,&n_,A,&lda_,S,NULL,&lda_,VT,&lda_,work,&lwork,&info));
251 : #endif
252 48 : SlepcCheckLapackInfo("gesvd",info);
253 131 : for (i=0;i<=k;i++) {
254 83 : ww[i] = PetscConj(VT[i*ndpt+k]);
255 83 : D[i] = ww[i]*f[i];
256 : }
257 48 : PetscCallBLAS("BLASgemv",BLASgemv_("N",&lda_,&n_,&sone,C,&lda_,D,&one,&szero,N,&one));
258 48 : PetscCallBLAS("BLASgemv",BLASgemv_("N",&lda_,&n_,&sone,C,&lda_,ww,&one,&szero,D,&one));
259 4848 : for (i=0;i<ndpt;i++) if (R[i]>=0) R[i] = PetscAbsScalar(F[i]-N[i]/D[i]);
260 : /* next support point */
261 : err = 0.0;
262 4848 : for (i=0;i<ndpt;i++) if (R[i]>=err) {idx = i; err = R[i];}
263 48 : if (err <= ctx->ddtol*norm) break;
264 : }
265 :
266 22 : PetscCheck(k<ndpt-1,PetscObjectComm((PetscObject)nep),PETSC_ERR_CONV_FAILED,"Failed to determine singularities automatically in general problem");
267 : /* poles */
268 22 : PetscCall(PetscArrayzero(C,ndpt*ndpt));
269 22 : PetscCall(PetscArrayzero(A,ndpt*ndpt));
270 70 : for (i=0;i<=k;i++) {
271 48 : C[i+ndpt*i] = 1.0;
272 48 : A[(i+1)*ndpt] = ww[i];
273 48 : A[i+1] = 1.0;
274 48 : A[i+1+(i+1)*ndpt] = z[i];
275 : }
276 22 : C[0] = 0.0; C[k+1+(k+1)*ndpt] = 1.0;
277 22 : n_++;
278 : #if defined(PETSC_USE_COMPLEX)
279 : PetscCallBLAS("LAPACK" LAPGEEV,LAPACKggevalt_("N","N",&n_,A,&lda_,C,&lda_,D,N,NULL,&lda_,NULL,&lda_,work,&lwork,rwork,&info));
280 : #else
281 22 : PetscCallBLAS("LAPACK" LAPGEEV,LAPACKggevalt_("N","N",&n_,A,&lda_,C,&lda_,D,VT,N,NULL,&lda_,NULL,&lda_,work,&lwork,&info));
282 : #endif
283 22 : SlepcCheckLapackInfo(LAPGEEV,info);
284 : cont = 0.0;
285 92 : for (i=0;i<n_;i++) if (N[i]!=0.0) {
286 20 : dxi[cont++] = D[i]/N[i];
287 : }
288 22 : *ndptx = cont;
289 22 : PetscCall(PetscFPTrapPop());
290 22 : PetscCall(PetscFree5(R,z,f,C,ww));
291 22 : PetscCall(PetscFree6(A,S,VT,work,D,N));
292 : #if defined(PETSC_USE_COMPLEX)
293 : PetscCall(PetscFree(rwork));
294 : #endif
295 22 : PetscFunctionReturn(PETSC_SUCCESS);
296 : }
297 :
298 : /* Singularities using Adaptive Anderson-Antoulas algorithm */
299 12 : static PetscErrorCode NEPNLEIGSAAASingularities(NEP nep,PetscInt ndpt,PetscScalar *ds,PetscInt *ndptx,PetscScalar *dxi)
300 : {
301 12 : Vec u,v,w;
302 12 : PetscRandom rand=NULL;
303 12 : PetscScalar *F,*isol;
304 12 : PetscInt i,k,nisol,nt;
305 12 : Mat T;
306 12 : FN f;
307 :
308 12 : PetscFunctionBegin;
309 12 : PetscCall(PetscMalloc1(ndpt,&F));
310 12 : if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
311 5 : PetscCall(PetscMalloc1(ndpt,&isol));
312 5 : *ndptx = 0;
313 5 : PetscCall(NEPGetSplitOperatorInfo(nep,&nt,NULL));
314 5 : nisol = *ndptx;
315 20 : for (k=0;k<nt;k++) {
316 15 : PetscCall(NEPGetSplitOperatorTerm(nep,k,NULL,&f));
317 1515 : for (i=0;i<ndpt;i++) PetscCall(FNEvaluateFunction(f,ds[i],&F[i]));
318 15 : PetscCall(NEPNLEIGSAAAComputation(nep,ndpt,ds,F,&nisol,isol));
319 15 : if (nisol) PetscCall(NEPNLEIGSAuxiliarRmDuplicates(nisol,isol,ndptx,dxi,ndpt));
320 : }
321 5 : PetscCall(PetscFree(isol));
322 : } else {
323 7 : PetscCall(MatCreateVecs(nep->function,&u,NULL));
324 7 : PetscCall(VecDuplicate(u,&v));
325 7 : PetscCall(VecDuplicate(u,&w));
326 7 : if (nep->V) PetscCall(BVGetRandomContext(nep->V,&rand));
327 7 : PetscCall(VecSetRandom(u,rand));
328 7 : PetscCall(VecNormalize(u,NULL));
329 7 : PetscCall(VecSetRandom(v,rand));
330 7 : PetscCall(VecNormalize(v,NULL));
331 7 : T = nep->function;
332 707 : for (i=0;i<ndpt;i++) {
333 700 : PetscCall(NEPComputeFunction(nep,ds[i],T,T));
334 700 : PetscCall(MatMult(T,v,w));
335 700 : PetscCall(VecDot(w,u,&F[i]));
336 : }
337 7 : PetscCall(NEPNLEIGSAAAComputation(nep,ndpt,ds,F,ndptx,dxi));
338 7 : PetscCall(VecDestroy(&u));
339 7 : PetscCall(VecDestroy(&v));
340 7 : PetscCall(VecDestroy(&w));
341 : }
342 12 : PetscCall(PetscFree(F));
343 12 : PetscFunctionReturn(PETSC_SUCCESS);
344 : }
345 :
346 27 : static PetscErrorCode NEPNLEIGSLejaBagbyPoints(NEP nep)
347 : {
348 27 : NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
349 27 : PetscInt i,k,ndpt=NDPOINTS,ndptx=NDPOINTS;
350 27 : PetscScalar *ds,*dsi,*dxi,*nrs,*nrxi,*s=ctx->s,*xi=ctx->xi,*beta=ctx->beta;
351 27 : PetscReal maxnrs,minnrxi;
352 27 : PetscBool rational;
353 : #if !defined(PETSC_USE_COMPLEX)
354 27 : PetscReal a,b,h;
355 : #endif
356 :
357 27 : PetscFunctionBegin;
358 27 : if (!ctx->computesingularities && nep->problem_type!=NEP_RATIONAL) ndpt = ndptx = LBPOINTS;
359 27 : PetscCall(PetscMalloc5(ndpt+1,&ds,ndpt+1,&dsi,ndpt,&dxi,ndpt+1,&nrs,ndpt,&nrxi));
360 :
361 : /* Discretize the target region boundary */
362 27 : PetscCall(RGComputeContour(nep->rg,ndpt,ds,dsi));
363 : #if !defined(PETSC_USE_COMPLEX)
364 131227 : for (i=0;i<ndpt;i++) if (dsi[i]!=0.0) break;
365 27 : if (i<ndpt) {
366 2 : PetscCheck(nep->problem_type==NEP_RATIONAL,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"NLEIGS with real arithmetic requires the target set to be included in the real axis");
367 : /* Select a segment in the real axis */
368 2 : PetscCall(RGComputeBoundingBox(nep->rg,&a,&b,NULL,NULL));
369 2 : PetscCheck(a>-PETSC_MAX_REAL && b<PETSC_MAX_REAL,PetscObjectComm((PetscObject)nep),PETSC_ERR_USER_INPUT,"NLEIGS requires a bounded target set");
370 2 : h = (b-a)/ndpt;
371 20002 : for (i=0;i<ndpt;i++) {ds[i] = a+h*i; dsi[i] = 0.0;}
372 : }
373 : #endif
374 : /* Discretize the singularity region */
375 27 : if (ctx->computesingularities) PetscCall(ctx->computesingularities(nep,&ndptx,dxi,ctx->singularitiesctx));
376 : else {
377 15 : if (nep->problem_type==NEP_RATIONAL) {
378 3 : PetscCall(NEPNLEIGSRationalSingularities(nep,&ndptx,dxi,&rational));
379 3 : PetscCheck(rational,PetscObjectComm((PetscObject)nep),PETSC_ERR_CONV_FAILED,"Failed to determine singularities automatically in rational problem; consider solving the problem as general");
380 : } else {
381 : /* AAA algorithm */
382 12 : PetscCall(NEPNLEIGSAAASingularities(nep,ndpt,ds,&ndptx,dxi));
383 : }
384 : }
385 : /* Look for Leja-Bagby points in the discretization sets */
386 27 : s[0] = ds[0];
387 27 : xi[0] = (ndptx>0)?dxi[0]:PETSC_INFINITY;
388 27 : PetscCheck(PetscAbsScalar(xi[0])>=10*PETSC_MACHINE_EPSILON,PetscObjectComm((PetscObject)nep),PETSC_ERR_USER_INPUT,"Singularity point 0 is nearly zero: %g; consider removing the singularity or shifting the problem",(double)PetscAbsScalar(xi[0]));
389 27 : beta[0] = 1.0; /* scaling factors are also computed here */
390 151227 : for (i=0;i<ndpt;i++) {
391 151200 : nrs[i] = 1.0;
392 151200 : nrxi[i] = 1.0;
393 : }
394 2325 : for (k=1;k<ctx->ddmaxit;k++) {
395 : maxnrs = 0.0;
396 12062598 : minnrxi = PETSC_MAX_REAL;
397 12062598 : for (i=0;i<ndpt;i++) {
398 12060300 : nrs[i] *= ((ds[i]-s[k-1])/(1.0-ds[i]/xi[k-1]))/beta[k-1];
399 12060300 : if (PetscAbsScalar(nrs[i])>maxnrs) {maxnrs = PetscAbsScalar(nrs[i]); s[k] = ds[i];}
400 : }
401 2298 : if (ndptx>k) {
402 8980016 : for (i=1;i<ndptx;i++) {
403 8979110 : nrxi[i] *= ((dxi[i]-s[k-1])/(1.0-dxi[i]/xi[k-1]))/beta[k-1];
404 8979110 : if (PetscAbsScalar(nrxi[i])<minnrxi) {minnrxi = PetscAbsScalar(nrxi[i]); xi[k] = dxi[i];}
405 : }
406 906 : PetscCheck(PetscAbsScalar(xi[k])>=10*PETSC_MACHINE_EPSILON,PetscObjectComm((PetscObject)nep),PETSC_ERR_USER_INPUT,"Singularity point %" PetscInt_FMT " is nearly zero: %g; consider removing the singularity or shifting the problem",k,(double)PetscAbsScalar(xi[k]));
407 1392 : } else xi[k] = PETSC_INFINITY;
408 2298 : beta[k] = maxnrs;
409 : }
410 27 : PetscCall(PetscFree5(ds,dsi,dxi,nrs,nrxi));
411 27 : PetscFunctionReturn(PETSC_SUCCESS);
412 : }
413 :
414 1206 : PetscErrorCode NEPNLEIGSEvalNRTFunct(NEP nep,PetscInt k,PetscScalar sigma,PetscScalar *b)
415 : {
416 1206 : NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
417 1206 : PetscInt i;
418 1206 : PetscScalar *beta=ctx->beta,*s=ctx->s,*xi=ctx->xi;
419 :
420 1206 : PetscFunctionBegin;
421 1206 : b[0] = 1.0/beta[0];
422 11295 : for (i=0;i<k;i++) {
423 10089 : b[i+1] = ((sigma-s[i])*b[i])/(beta[i+1]*(1.0-sigma/xi[i]));
424 : }
425 1206 : PetscFunctionReturn(PETSC_SUCCESS);
426 : }
427 :
428 18022 : static PetscErrorCode MatMult_Fun(Mat A,Vec x,Vec y)
429 : {
430 18022 : NEP_NLEIGS_MATSHELL *ctx;
431 18022 : PetscInt i;
432 :
433 18022 : PetscFunctionBeginUser;
434 18022 : PetscCall(MatShellGetContext(A,&ctx));
435 18022 : PetscCall(MatMult(ctx->A[0],x,y));
436 18022 : if (ctx->coeff[0]!=1.0) PetscCall(VecScale(y,ctx->coeff[0]));
437 217565 : for (i=1;i<ctx->nmat;i++) {
438 199543 : PetscCall(MatMult(ctx->A[i],x,ctx->t));
439 199543 : PetscCall(VecAXPY(y,ctx->coeff[i],ctx->t));
440 : }
441 18022 : PetscFunctionReturn(PETSC_SUCCESS);
442 : }
443 :
444 16237 : static PetscErrorCode MatMultTranspose_Fun(Mat A,Vec x,Vec y)
445 : {
446 16237 : NEP_NLEIGS_MATSHELL *ctx;
447 16237 : PetscInt i;
448 :
449 16237 : PetscFunctionBeginUser;
450 16237 : PetscCall(MatShellGetContext(A,&ctx));
451 16237 : PetscCall(MatMultTranspose(ctx->A[0],x,y));
452 16237 : if (ctx->coeff[0]!=1.0) PetscCall(VecScale(y,ctx->coeff[0]));
453 196925 : for (i=1;i<ctx->nmat;i++) {
454 180688 : PetscCall(MatMultTranspose(ctx->A[i],x,ctx->t));
455 180688 : PetscCall(VecAXPY(y,ctx->coeff[i],ctx->t));
456 : }
457 16237 : PetscFunctionReturn(PETSC_SUCCESS);
458 : }
459 :
460 6 : static PetscErrorCode MatGetDiagonal_Fun(Mat A,Vec diag)
461 : {
462 6 : NEP_NLEIGS_MATSHELL *ctx;
463 6 : PetscInt i;
464 :
465 6 : PetscFunctionBeginUser;
466 6 : PetscCall(MatShellGetContext(A,&ctx));
467 6 : PetscCall(MatGetDiagonal(ctx->A[0],diag));
468 6 : if (ctx->coeff[0]!=1.0) PetscCall(VecScale(diag,ctx->coeff[0]));
469 75 : for (i=1;i<ctx->nmat;i++) {
470 69 : PetscCall(MatGetDiagonal(ctx->A[i],ctx->t));
471 69 : PetscCall(VecAXPY(diag,ctx->coeff[i],ctx->t));
472 : }
473 6 : PetscFunctionReturn(PETSC_SUCCESS);
474 : }
475 :
476 3 : static PetscErrorCode MatDuplicate_Fun(Mat A,MatDuplicateOption op,Mat *B)
477 : {
478 3 : PetscInt m,n,M,N,i;
479 3 : NEP_NLEIGS_MATSHELL *ctxnew,*ctx;
480 3 : void (*fun)(void);
481 :
482 3 : PetscFunctionBeginUser;
483 3 : PetscCall(MatShellGetContext(A,&ctx));
484 3 : PetscCall(PetscNew(&ctxnew));
485 3 : ctxnew->nmat = ctx->nmat;
486 3 : ctxnew->maxnmat = ctx->maxnmat;
487 3 : PetscCall(PetscMalloc2(ctxnew->maxnmat,&ctxnew->A,ctxnew->maxnmat,&ctxnew->coeff));
488 6 : for (i=0;i<ctx->nmat;i++) {
489 3 : PetscCall(PetscObjectReference((PetscObject)ctx->A[i]));
490 3 : ctxnew->A[i] = ctx->A[i];
491 3 : ctxnew->coeff[i] = ctx->coeff[i];
492 : }
493 3 : PetscCall(MatGetSize(ctx->A[0],&M,&N));
494 3 : PetscCall(MatGetLocalSize(ctx->A[0],&m,&n));
495 3 : PetscCall(VecDuplicate(ctx->t,&ctxnew->t));
496 3 : PetscCall(MatCreateShell(PetscObjectComm((PetscObject)A),m,n,M,N,(void*)ctxnew,B));
497 3 : PetscCall(MatShellSetManageScalingShifts(*B));
498 3 : PetscCall(MatShellGetOperation(A,MATOP_MULT,&fun));
499 3 : PetscCall(MatShellSetOperation(*B,MATOP_MULT,fun));
500 3 : PetscCall(MatShellGetOperation(A,MATOP_MULT_TRANSPOSE,&fun));
501 3 : PetscCall(MatShellSetOperation(*B,MATOP_MULT_TRANSPOSE,fun));
502 3 : PetscCall(MatShellGetOperation(A,MATOP_GET_DIAGONAL,&fun));
503 3 : PetscCall(MatShellSetOperation(*B,MATOP_GET_DIAGONAL,fun));
504 3 : PetscCall(MatShellGetOperation(A,MATOP_DUPLICATE,&fun));
505 3 : PetscCall(MatShellSetOperation(*B,MATOP_DUPLICATE,fun));
506 3 : PetscCall(MatShellGetOperation(A,MATOP_DESTROY,&fun));
507 3 : PetscCall(MatShellSetOperation(*B,MATOP_DESTROY,fun));
508 3 : PetscCall(MatShellGetOperation(A,MATOP_AXPY,&fun));
509 3 : PetscCall(MatShellSetOperation(*B,MATOP_AXPY,fun));
510 3 : PetscFunctionReturn(PETSC_SUCCESS);
511 : }
512 :
513 78 : static PetscErrorCode MatDestroy_Fun(Mat A)
514 : {
515 78 : NEP_NLEIGS_MATSHELL *ctx;
516 78 : PetscInt i;
517 :
518 78 : PetscFunctionBeginUser;
519 78 : if (A) {
520 78 : PetscCall(MatShellGetContext(A,&ctx));
521 984 : for (i=0;i<ctx->nmat;i++) PetscCall(MatDestroy(&ctx->A[i]));
522 78 : PetscCall(VecDestroy(&ctx->t));
523 78 : PetscCall(PetscFree2(ctx->A,ctx->coeff));
524 78 : PetscCall(PetscFree(ctx));
525 : }
526 78 : PetscFunctionReturn(PETSC_SUCCESS);
527 : }
528 :
529 828 : static PetscErrorCode MatAXPY_Fun(Mat Y,PetscScalar a,Mat X,MatStructure str)
530 : {
531 828 : NEP_NLEIGS_MATSHELL *ctxY,*ctxX;
532 828 : PetscInt i,j;
533 828 : PetscBool found;
534 :
535 828 : PetscFunctionBeginUser;
536 828 : PetscCall(MatShellGetContext(Y,&ctxY));
537 828 : PetscCall(MatShellGetContext(X,&ctxX));
538 7728 : for (i=0;i<ctxX->nmat;i++) {
539 : found = PETSC_FALSE;
540 56997 : for (j=0;!found&&j<ctxY->nmat;j++) {
541 50097 : if (ctxX->A[i]==ctxY->A[j]) {
542 6072 : found = PETSC_TRUE;
543 6072 : ctxY->coeff[j] += a*ctxX->coeff[i];
544 : }
545 : }
546 6900 : if (!found) {
547 828 : ctxY->coeff[ctxY->nmat] = a*ctxX->coeff[i];
548 828 : ctxY->A[ctxY->nmat++] = ctxX->A[i];
549 6900 : PetscCall(PetscObjectReference((PetscObject)ctxX->A[i]));
550 : }
551 : }
552 828 : PetscFunctionReturn(PETSC_SUCCESS);
553 : }
554 :
555 57 : static PetscErrorCode MatScale_Fun(Mat M,PetscScalar a)
556 : {
557 57 : NEP_NLEIGS_MATSHELL *ctx;
558 57 : PetscInt i;
559 :
560 57 : PetscFunctionBeginUser;
561 57 : PetscCall(MatShellGetContext(M,&ctx));
562 792 : for (i=0;i<ctx->nmat;i++) ctx->coeff[i] *= a;
563 57 : PetscFunctionReturn(PETSC_SUCCESS);
564 : }
565 :
566 75 : static PetscErrorCode NLEIGSMatToMatShellArray(Mat A,Mat *Ms,PetscInt maxnmat)
567 : {
568 75 : NEP_NLEIGS_MATSHELL *ctx;
569 75 : PetscInt m,n,M,N;
570 75 : PetscBool has;
571 :
572 75 : PetscFunctionBegin;
573 75 : PetscCall(MatHasOperation(A,MATOP_DUPLICATE,&has));
574 75 : PetscCheck(has,PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"MatDuplicate operation required");
575 75 : PetscCall(PetscNew(&ctx));
576 75 : ctx->maxnmat = maxnmat;
577 75 : PetscCall(PetscMalloc2(ctx->maxnmat,&ctx->A,ctx->maxnmat,&ctx->coeff));
578 75 : PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&ctx->A[0]));
579 75 : ctx->nmat = 1;
580 75 : ctx->coeff[0] = 1.0;
581 75 : PetscCall(MatCreateVecs(A,&ctx->t,NULL));
582 75 : PetscCall(MatGetSize(A,&M,&N));
583 75 : PetscCall(MatGetLocalSize(A,&m,&n));
584 75 : PetscCall(MatCreateShell(PetscObjectComm((PetscObject)A),m,n,M,N,(void*)ctx,Ms));
585 75 : PetscCall(MatShellSetManageScalingShifts(*Ms));
586 75 : PetscCall(MatShellSetOperation(*Ms,MATOP_MULT,(void(*)(void))MatMult_Fun));
587 75 : PetscCall(MatShellSetOperation(*Ms,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_Fun));
588 75 : PetscCall(MatShellSetOperation(*Ms,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Fun));
589 75 : PetscCall(MatShellSetOperation(*Ms,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_Fun));
590 75 : PetscCall(MatShellSetOperation(*Ms,MATOP_DESTROY,(void(*)(void))MatDestroy_Fun));
591 75 : PetscCall(MatShellSetOperation(*Ms,MATOP_AXPY,(void(*)(void))MatAXPY_Fun));
592 75 : PetscCall(MatShellSetOperation(*Ms,MATOP_SCALE,(void(*)(void))MatScale_Fun));
593 75 : PetscFunctionReturn(PETSC_SUCCESS);
594 : }
595 :
596 : /*
597 : MatIsShellAny - returns true if any of the n matrices is a shell matrix
598 : */
599 15 : static PetscErrorCode MatIsShellAny(Mat *A,PetscInt n,PetscBool *shell)
600 : {
601 15 : PetscInt i;
602 15 : PetscBool flg;
603 :
604 15 : PetscFunctionBegin;
605 15 : *shell = PETSC_FALSE;
606 47 : for (i=0;i<n;i++) {
607 35 : PetscCall(MatIsShell(A[i],&flg));
608 35 : if (flg) { *shell = PETSC_TRUE; break; }
609 : }
610 15 : PetscFunctionReturn(PETSC_SUCCESS);
611 : }
612 :
613 15 : static PetscErrorCode NEPNLEIGSDividedDifferences_split(NEP nep)
614 : {
615 15 : PetscErrorCode ierr;
616 15 : NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
617 15 : PetscInt k,j,i,maxnmat,nmax;
618 15 : PetscReal norm0,norm,*matnorm;
619 15 : PetscScalar *s=ctx->s,*beta=ctx->beta,*xi=ctx->xi,*b,alpha,*coeffs,*pK,*pH,sone=1.0;
620 15 : Mat T,P,Ts,K,H;
621 15 : PetscBool shell,hasmnorm=PETSC_FALSE,matrix=PETSC_TRUE;
622 15 : PetscBLASInt n_;
623 :
624 15 : PetscFunctionBegin;
625 15 : nmax = ctx->ddmaxit;
626 15 : PetscCall(PetscMalloc1(nep->nt*nmax,&ctx->coeffD));
627 15 : PetscCall(PetscMalloc3(nmax+1,&b,nmax+1,&coeffs,nep->nt,&matnorm));
628 47 : for (j=0;j<nep->nt;j++) {
629 35 : PetscCall(MatHasOperation(nep->A[j],MATOP_NORM,&hasmnorm));
630 35 : if (!hasmnorm) break;
631 32 : PetscCall(MatNorm(nep->A[j],NORM_INFINITY,matnorm+j));
632 : }
633 : /* Try matrix functions scheme */
634 15 : PetscCall(PetscCalloc2(nmax*nmax,&pK,nmax*nmax,&pH));
635 1215 : for (i=0;i<nmax-1;i++) {
636 1200 : pK[(nmax+1)*i] = 1.0;
637 1200 : pK[(nmax+1)*i+1] = beta[i+1]/xi[i];
638 1200 : pH[(nmax+1)*i] = s[i];
639 1200 : pH[(nmax+1)*i+1] = beta[i+1];
640 : }
641 15 : pH[nmax*nmax-1] = s[nmax-1];
642 15 : pK[nmax*nmax-1] = 1.0;
643 15 : PetscCall(PetscBLASIntCast(nmax,&n_));
644 15 : PetscCallBLAS("BLAStrsm",BLAStrsm_("R","L","N","U",&n_,&n_,&sone,pK,&n_,pH,&n_));
645 : /* The matrix to be used is in H. K will be a work-space matrix */
646 15 : PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,nmax,nmax,pH,&H));
647 15 : PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,nmax,nmax,pK,&K));
648 53 : for (j=0;matrix&&j<nep->nt;j++) {
649 38 : PetscCall(PetscPushErrorHandler(PetscReturnErrorHandler,NULL));
650 38 : ierr = FNEvaluateFunctionMat(nep->f[j],H,K);
651 38 : PetscCall(PetscPopErrorHandler());
652 38 : if (!ierr) {
653 3183 : for (i=0;i<nmax;i++) ctx->coeffD[j+i*nep->nt] = pK[i]*beta[0];
654 : } else {
655 0 : matrix = PETSC_FALSE;
656 38 : PetscCall(PetscFPTrapPop());
657 : }
658 : }
659 15 : PetscCall(MatDestroy(&H));
660 15 : PetscCall(MatDestroy(&K));
661 15 : if (!matrix) {
662 0 : for (j=0;j<nep->nt;j++) {
663 0 : PetscCall(FNEvaluateFunction(nep->f[j],s[0],ctx->coeffD+j));
664 0 : ctx->coeffD[j] *= beta[0];
665 : }
666 : }
667 15 : if (hasmnorm) {
668 : norm0 = 0.0;
669 44 : for (j=0;j<nep->nt;j++) norm0 += matnorm[j]*PetscAbsScalar(ctx->coeffD[j]);
670 : } else {
671 : norm0 = 0.0;
672 9 : for (j=0;j<nep->nt;j++) norm0 = PetscMax(PetscAbsScalar(ctx->coeffD[j]),norm0);
673 : }
674 15 : ctx->nmat = ctx->ddmaxit;
675 167 : for (k=1;k<ctx->ddmaxit;k++) {
676 165 : if (!matrix) {
677 0 : PetscCall(NEPNLEIGSEvalNRTFunct(nep,k,s[k],b));
678 0 : for (i=0;i<nep->nt;i++) {
679 0 : PetscCall(FNEvaluateFunction(nep->f[i],s[k],ctx->coeffD+k*nep->nt+i));
680 0 : for (j=0;j<k;j++) {
681 0 : ctx->coeffD[k*nep->nt+i] -= b[j]*ctx->coeffD[i+nep->nt*j];
682 : }
683 0 : ctx->coeffD[k*nep->nt+i] /= b[k];
684 : }
685 : }
686 165 : if (hasmnorm) {
687 : norm = 0.0;
688 342 : for (j=0;j<nep->nt;j++) norm += matnorm[j]*PetscAbsScalar(ctx->coeffD[k*nep->nt+j]);
689 : } else {
690 : norm = 0.0;
691 180 : for (j=0;j<nep->nt;j++) norm = PetscMax(PetscAbsScalar(ctx->coeffD[k*nep->nt+j]),norm);
692 : }
693 165 : if (k>1 && norm/norm0 < ctx->ddtol) {
694 13 : ctx->nmat = k+1;
695 13 : break;
696 : }
697 : }
698 15 : if (!ctx->ksp) PetscCall(NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp));
699 15 : PetscCall(MatIsShellAny(nep->A,nep->nt,&shell));
700 15 : maxnmat = PetscMax(ctx->ddmaxit,nep->nt);
701 36 : for (i=0;i<ctx->nshiftsw;i++) {
702 21 : PetscCall(NEPNLEIGSEvalNRTFunct(nep,ctx->nmat-1,ctx->shifts[i],coeffs));
703 21 : if (!shell) PetscCall(MatDuplicate(nep->A[0],MAT_COPY_VALUES,&T));
704 3 : else PetscCall(NLEIGSMatToMatShellArray(nep->A[0],&T,maxnmat));
705 21 : if (nep->P) { /* user-defined preconditioner */
706 1 : PetscCall(MatDuplicate(nep->P[0],MAT_COPY_VALUES,&P));
707 20 : } else P=T;
708 21 : alpha = 0.0;
709 321 : for (j=0;j<ctx->nmat;j++) alpha += coeffs[j]*ctx->coeffD[j*nep->nt];
710 21 : PetscCall(MatScale(T,alpha));
711 21 : if (nep->P) PetscCall(MatScale(P,alpha));
712 50 : for (k=1;k<nep->nt;k++) {
713 : alpha = 0.0;
714 364 : for (j=0;j<ctx->nmat;j++) alpha += coeffs[j]*ctx->coeffD[j*nep->nt+k];
715 29 : if (shell) PetscCall(NLEIGSMatToMatShellArray(nep->A[k],&Ts,maxnmat));
716 29 : PetscCall(MatAXPY(T,alpha,shell?Ts:nep->A[k],nep->mstr));
717 29 : if (nep->P) PetscCall(MatAXPY(P,alpha,nep->P[k],nep->mstrp));
718 29 : if (shell) PetscCall(MatDestroy(&Ts));
719 : }
720 21 : PetscCall(NEP_KSPSetOperators(ctx->ksp[i],T,P));
721 21 : PetscCall(KSPSetUp(ctx->ksp[i]));
722 21 : PetscCall(MatDestroy(&T));
723 21 : if (nep->P) PetscCall(MatDestroy(&P));
724 : }
725 15 : PetscCall(PetscFree3(b,coeffs,matnorm));
726 15 : PetscCall(PetscFree2(pK,pH));
727 15 : PetscFunctionReturn(PETSC_SUCCESS);
728 : }
729 :
730 12 : static PetscErrorCode NEPNLEIGSDividedDifferences_callback(NEP nep)
731 : {
732 12 : NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
733 12 : PetscInt k,j,i,maxnmat;
734 12 : PetscReal norm0,norm;
735 12 : PetscScalar *s=ctx->s,*beta=ctx->beta,*b,*coeffs;
736 12 : Mat *D=ctx->D,*DP,T,P;
737 12 : PetscBool shell,has,vec=PETSC_FALSE,precond=(nep->function_pre!=nep->function)?PETSC_TRUE:PETSC_FALSE;
738 12 : PetscRandom rand=NULL;
739 12 : Vec w[2];
740 :
741 12 : PetscFunctionBegin;
742 12 : PetscCall(PetscMalloc2(ctx->ddmaxit+1,&b,ctx->ddmaxit+1,&coeffs));
743 12 : if (nep->V) PetscCall(BVGetRandomContext(nep->V,&rand));
744 12 : T = nep->function;
745 12 : P = nep->function_pre;
746 12 : PetscCall(NEPComputeFunction(nep,s[0],T,P));
747 12 : PetscCall(MatIsShell(T,&shell));
748 12 : maxnmat = PetscMax(ctx->ddmaxit,nep->nt);
749 12 : if (!shell) PetscCall(MatDuplicate(T,MAT_COPY_VALUES,&D[0]));
750 3 : else PetscCall(NLEIGSMatToMatShellArray(T,&D[0],maxnmat));
751 12 : if (beta[0]!=1.0) PetscCall(MatScale(D[0],1.0/beta[0]));
752 12 : PetscCall(MatHasOperation(D[0],MATOP_NORM,&has));
753 12 : if (has) PetscCall(MatNorm(D[0],NORM_FROBENIUS,&norm0));
754 : else {
755 3 : PetscCall(MatCreateVecs(D[0],NULL,&w[0]));
756 3 : PetscCall(VecDuplicate(w[0],&w[1]));
757 3 : PetscCall(VecDuplicate(w[0],&ctx->vrn));
758 3 : PetscCall(VecSetRandomNormal(ctx->vrn,rand,w[0],w[1]));
759 3 : PetscCall(VecNormalize(ctx->vrn,NULL));
760 3 : vec = PETSC_TRUE;
761 3 : PetscCall(MatNormEstimate(D[0],ctx->vrn,w[0],&norm0));
762 : }
763 12 : if (precond) {
764 1 : PetscCall(PetscMalloc1(ctx->ddmaxit,&DP));
765 1 : PetscCall(MatDuplicate(P,MAT_COPY_VALUES,&DP[0]));
766 : }
767 12 : ctx->nmat = ctx->ddmaxit;
768 129 : for (k=1;k<ctx->ddmaxit;k++) {
769 128 : PetscCall(NEPNLEIGSEvalNRTFunct(nep,k,s[k],b));
770 128 : PetscCall(NEPComputeFunction(nep,s[k],T,P));
771 128 : if (!shell) PetscCall(MatDuplicate(T,MAT_COPY_VALUES,&D[k]));
772 66 : else PetscCall(NLEIGSMatToMatShellArray(T,&D[k],maxnmat));
773 1372 : for (j=0;j<k;j++) PetscCall(MatAXPY(D[k],-b[j],D[j],nep->mstr));
774 128 : PetscCall(MatScale(D[k],1.0/b[k]));
775 128 : PetscCall(MatHasOperation(D[k],MATOP_NORM,&has));
776 128 : if (has) PetscCall(MatNorm(D[k],NORM_FROBENIUS,&norm));
777 : else {
778 66 : if (!vec) {
779 0 : PetscCall(MatCreateVecs(D[k],NULL,&w[0]));
780 0 : PetscCall(VecDuplicate(w[0],&w[1]));
781 0 : PetscCall(VecDuplicate(w[0],&ctx->vrn));
782 0 : PetscCall(VecSetRandomNormal(ctx->vrn,rand,w[0],w[1]));
783 0 : PetscCall(VecNormalize(ctx->vrn,NULL));
784 : vec = PETSC_TRUE;
785 : }
786 66 : PetscCall(MatNormEstimate(D[k],ctx->vrn,w[0],&norm));
787 : }
788 128 : if (precond) {
789 3 : PetscCall(MatDuplicate(P,MAT_COPY_VALUES,&DP[k]));
790 9 : for (j=0;j<k;j++) PetscCall(MatAXPY(DP[k],-b[j],DP[j],nep->mstrp));
791 3 : PetscCall(MatScale(DP[k],1.0/b[k]));
792 : }
793 128 : if (k>1 && norm/norm0 < ctx->ddtol && k>1) {
794 11 : ctx->nmat = k+1;
795 11 : break;
796 : }
797 : }
798 12 : if (!ctx->ksp) PetscCall(NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp));
799 27 : for (i=0;i<ctx->nshiftsw;i++) {
800 15 : PetscCall(NEPNLEIGSEvalNRTFunct(nep,ctx->nmat-1,ctx->shifts[i],coeffs));
801 15 : PetscCall(MatDuplicate(D[0],MAT_COPY_VALUES,&T));
802 15 : if (coeffs[0]!=1.0) PetscCall(MatScale(T,coeffs[0]));
803 200 : for (j=1;j<ctx->nmat;j++) PetscCall(MatAXPY(T,coeffs[j],D[j],nep->mstr));
804 15 : if (precond) {
805 1 : PetscCall(MatDuplicate(DP[0],MAT_COPY_VALUES,&P));
806 1 : if (coeffs[0]!=1.0) PetscCall(MatScale(P,coeffs[0]));
807 4 : for (j=1;j<ctx->nmat;j++) PetscCall(MatAXPY(P,coeffs[j],DP[j],nep->mstrp));
808 14 : } else P=T;
809 15 : PetscCall(NEP_KSPSetOperators(ctx->ksp[i],T,P));
810 15 : PetscCall(KSPSetUp(ctx->ksp[i]));
811 15 : PetscCall(MatDestroy(&T));
812 : }
813 12 : PetscCall(PetscFree2(b,coeffs));
814 12 : if (vec) {
815 3 : PetscCall(VecDestroy(&w[0]));
816 3 : PetscCall(VecDestroy(&w[1]));
817 : }
818 12 : if (precond) {
819 1 : PetscCall(MatDestroy(&P));
820 1 : PetscCall(MatDestroyMatrices(ctx->nmat,&DP));
821 : }
822 12 : PetscFunctionReturn(PETSC_SUCCESS);
823 : }
824 :
825 : /*
826 : NEPKrylovConvergence - This is the analogue to EPSKrylovConvergence.
827 : */
828 72 : static PetscErrorCode NEPNLEIGSKrylovConvergence(NEP nep,PetscBool getall,PetscInt kini,PetscInt nits,PetscReal betah,PetscScalar betak,PetscInt *kout,Vec *w)
829 : {
830 72 : PetscInt k,newk,marker,inside;
831 72 : PetscScalar re,im;
832 72 : PetscReal resnorm,tt;
833 72 : PetscBool istrivial;
834 72 : NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
835 :
836 72 : PetscFunctionBegin;
837 72 : PetscCall(RGIsTrivial(nep->rg,&istrivial));
838 72 : marker = -1;
839 72 : if (nep->trackall) getall = PETSC_TRUE;
840 177 : for (k=kini;k<kini+nits;k++) {
841 : /* eigenvalue */
842 177 : re = nep->eigr[k];
843 177 : im = nep->eigi[k];
844 177 : if (!istrivial) {
845 177 : if (!ctx->nshifts) PetscCall(NEPNLEIGSBackTransform((PetscObject)nep,1,&re,&im));
846 177 : PetscCall(RGCheckInside(nep->rg,1,&re,&im,&inside));
847 177 : if (marker==-1 && inside<0) marker = k;
848 : }
849 177 : newk = k;
850 177 : PetscCall(DSVectors(nep->ds,DS_MAT_X,&newk,&resnorm));
851 177 : tt = (ctx->nshifts)?SlepcAbsEigenvalue(betak-nep->eigr[k]*betah,nep->eigi[k]*betah):betah;
852 177 : resnorm *= PetscAbsReal(tt);
853 : /* error estimate */
854 177 : PetscCall((*nep->converged)(nep,nep->eigr[k],nep->eigi[k],resnorm,&nep->errest[k],nep->convergedctx));
855 177 : if (marker==-1 && nep->errest[k] >= nep->tol) marker = k;
856 177 : if (newk==k+1) {
857 0 : nep->errest[k+1] = nep->errest[k];
858 0 : k++;
859 : }
860 177 : if (marker!=-1 && !getall) break;
861 : }
862 72 : if (marker!=-1) k = marker;
863 72 : *kout = k;
864 72 : PetscFunctionReturn(PETSC_SUCCESS);
865 : }
866 :
867 27 : static PetscErrorCode NEPSetUp_NLEIGS(NEP nep)
868 : {
869 27 : PetscInt k,in;
870 27 : PetscScalar zero=0.0;
871 27 : NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
872 27 : SlepcSC sc;
873 27 : PetscBool istrivial;
874 :
875 27 : PetscFunctionBegin;
876 27 : PetscCall(NEPSetDimensions_Default(nep,nep->nev,&nep->ncv,&nep->mpd));
877 27 : PetscCheck(nep->ncv<=nep->nev+nep->mpd,PetscObjectComm((PetscObject)nep),PETSC_ERR_USER_INPUT,"The value of ncv must not be larger than nev+mpd");
878 27 : if (nep->max_it==PETSC_DETERMINE) nep->max_it = PetscMax(5000,2*nep->n/nep->ncv);
879 27 : if (!ctx->ddmaxit) ctx->ddmaxit = LBPOINTS;
880 27 : PetscCall(RGIsTrivial(nep->rg,&istrivial));
881 27 : PetscCheck(!istrivial,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"NEPNLEIGS requires a nontrivial region defining the target set");
882 27 : if (!nep->which) nep->which = NEP_TARGET_MAGNITUDE;
883 27 : PetscCheck(nep->which==NEP_TARGET_MAGNITUDE || nep->which==NEP_TARGET_REAL || nep->which==NEP_TARGET_IMAGINARY || nep->which==NEP_WHICH_USER,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"This solver supports only target selection of eigenvalues");
884 :
885 : /* Initialize the NLEIGS context structure */
886 27 : k = ctx->ddmaxit;
887 27 : PetscCall(PetscMalloc4(k,&ctx->s,k,&ctx->xi,k,&ctx->beta,k,&ctx->D));
888 27 : nep->data = ctx;
889 27 : if (nep->tol==(PetscReal)PETSC_DETERMINE) nep->tol = SLEPC_DEFAULT_TOL;
890 27 : if (ctx->ddtol==(PetscReal)PETSC_DETERMINE) ctx->ddtol = nep->tol/10.0;
891 27 : if (!ctx->keep) ctx->keep = 0.5;
892 :
893 : /* Compute Leja-Bagby points and scaling values */
894 27 : PetscCall(NEPNLEIGSLejaBagbyPoints(nep));
895 27 : if (nep->problem_type!=NEP_RATIONAL) {
896 24 : PetscCall(RGCheckInside(nep->rg,1,&nep->target,&zero,&in));
897 24 : PetscCheck(in>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"The target is not inside the target set");
898 : }
899 :
900 : /* Compute the divided difference matrices */
901 27 : if (nep->fui==NEP_USER_INTERFACE_SPLIT) PetscCall(NEPNLEIGSDividedDifferences_split(nep));
902 12 : else PetscCall(NEPNLEIGSDividedDifferences_callback(nep));
903 27 : PetscCall(NEPAllocateSolution(nep,ctx->nmat-1));
904 27 : PetscCall(NEPSetWorkVecs(nep,4));
905 27 : if (!ctx->fullbasis) {
906 24 : PetscCheck(!nep->twosided,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Two-sided variant requires the full-basis option, rerun with -nep_nleigs_full_basis");
907 : /* set-up DS and transfer split operator functions */
908 45 : PetscCall(DSSetType(nep->ds,ctx->nshifts?DSGNHEP:DSNHEP));
909 24 : PetscCall(DSAllocate(nep->ds,nep->ncv+1));
910 24 : PetscCall(DSGetSlepcSC(nep->ds,&sc));
911 24 : if (!ctx->nshifts) sc->map = NEPNLEIGSBackTransform;
912 24 : PetscCall(DSSetExtraRow(nep->ds,PETSC_TRUE));
913 24 : sc->mapobj = (PetscObject)nep;
914 24 : sc->rg = nep->rg;
915 24 : sc->comparison = nep->sc->comparison;
916 24 : sc->comparisonctx = nep->sc->comparisonctx;
917 24 : PetscCall(BVDestroy(&ctx->V));
918 24 : PetscCall(BVCreateTensor(nep->V,ctx->nmat-1,&ctx->V));
919 24 : nep->ops->solve = NEPSolve_NLEIGS;
920 24 : nep->ops->computevectors = NEPComputeVectors_Schur;
921 : } else {
922 3 : PetscCall(NEPSetUp_NLEIGS_FullBasis(nep));
923 3 : nep->ops->solve = NEPSolve_NLEIGS_FullBasis;
924 3 : nep->ops->computevectors = NULL;
925 : }
926 27 : PetscFunctionReturn(PETSC_SUCCESS);
927 : }
928 :
929 : /*
930 : Extend the TOAR basis by applying the matrix operator
931 : over a vector which is decomposed on the TOAR way
932 : Input:
933 : - S,V: define the latest Arnoldi vector (nv vectors in V)
934 : Output:
935 : - t: new vector extending the TOAR basis
936 : - r: temporally coefficients to compute the TOAR coefficients
937 : for the new Arnoldi vector
938 : Workspace: t_ (two vectors)
939 : */
940 890 : static PetscErrorCode NEPTOARExtendBasis(NEP nep,PetscInt idxrktg,PetscScalar *S,PetscInt ls,PetscInt nv,BV W,BV V,Vec t,PetscScalar *r,PetscInt lr,Vec *t_)
941 : {
942 890 : NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
943 890 : PetscInt deg=ctx->nmat-1,k,j;
944 890 : Vec v=t_[0],q=t_[1],w;
945 890 : PetscScalar *beta=ctx->beta,*s=ctx->s,*xi=ctx->xi,*coeffs,sigma;
946 :
947 890 : PetscFunctionBegin;
948 890 : if (!ctx->ksp) PetscCall(NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp));
949 890 : sigma = ctx->shifts[idxrktg];
950 890 : PetscCall(BVSetActiveColumns(nep->V,0,nv));
951 890 : PetscCall(PetscMalloc1(ctx->nmat,&coeffs));
952 890 : PetscCheck(PetscAbsScalar(s[deg-2]-sigma)>100*PETSC_MACHINE_EPSILON,PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"Breakdown in NLEIGS");
953 : /* i-part stored in (i-1) position */
954 23084 : for (j=0;j<nv;j++) {
955 22194 : r[(deg-2)*lr+j] = (S[(deg-2)*ls+j]+(beta[deg-1]/xi[deg-2])*S[(deg-1)*ls+j])/(s[deg-2]-sigma);
956 : }
957 890 : PetscCall(BVSetActiveColumns(W,0,deg));
958 890 : PetscCall(BVGetColumn(W,deg-1,&w));
959 890 : PetscCall(BVMultVec(V,1.0/beta[deg],0,w,S+(deg-1)*ls));
960 890 : PetscCall(BVRestoreColumn(W,deg-1,&w));
961 890 : PetscCall(BVGetColumn(W,deg-2,&w));
962 890 : PetscCall(BVMultVec(V,1.0,0.0,w,r+(deg-2)*lr));
963 890 : PetscCall(BVRestoreColumn(W,deg-2,&w));
964 8009 : for (k=deg-2;k>0;k--) {
965 7119 : PetscCheck(PetscAbsScalar(s[k-1]-sigma)>100*PETSC_MACHINE_EPSILON,PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"Breakdown in NLEIGS");
966 209473 : for (j=0;j<nv;j++) r[(k-1)*lr+j] = (S[(k-1)*ls+j]+(beta[k]/xi[k-1])*S[k*ls+j]-beta[k]*(1.0-sigma/xi[k-1])*r[(k)*lr+j])/(s[k-1]-sigma);
967 7119 : PetscCall(BVGetColumn(W,k-1,&w));
968 7119 : PetscCall(BVMultVec(V,1.0,0.0,w,r+(k-1)*lr));
969 7119 : PetscCall(BVRestoreColumn(W,k-1,&w));
970 : }
971 890 : if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
972 5682 : for (j=0;j<ctx->nmat-2;j++) coeffs[j] = ctx->coeffD[nep->nt*j];
973 615 : coeffs[ctx->nmat-2] = ctx->coeffD[nep->nt*(ctx->nmat-1)];
974 615 : PetscCall(BVMultVec(W,1.0,0.0,v,coeffs));
975 615 : PetscCall(MatMult(nep->A[0],v,q));
976 1634 : for (k=1;k<nep->nt;k++) {
977 7188 : for (j=0;j<ctx->nmat-2;j++) coeffs[j] = ctx->coeffD[nep->nt*j+k];
978 1019 : coeffs[ctx->nmat-2] = ctx->coeffD[nep->nt*(ctx->nmat-1)+k];
979 1019 : PetscCall(BVMultVec(W,1.0,0,v,coeffs));
980 1019 : PetscCall(MatMult(nep->A[k],v,t));
981 1019 : PetscCall(VecAXPY(q,1.0,t));
982 : }
983 615 : PetscCall(KSPSolve(ctx->ksp[idxrktg],q,t));
984 615 : PetscCall(VecScale(t,-1.0));
985 : } else {
986 3217 : for (k=0;k<deg-1;k++) {
987 2942 : PetscCall(BVGetColumn(W,k,&w));
988 2942 : PetscCall(MatMult(ctx->D[k],w,q));
989 2942 : PetscCall(BVRestoreColumn(W,k,&w));
990 2942 : PetscCall(BVInsertVec(W,k,q));
991 : }
992 275 : PetscCall(BVGetColumn(W,deg-1,&w));
993 275 : PetscCall(MatMult(ctx->D[deg],w,q));
994 275 : PetscCall(BVRestoreColumn(W,k,&w));
995 275 : PetscCall(BVInsertVec(W,k,q));
996 3492 : for (j=0;j<ctx->nmat-1;j++) coeffs[j] = 1.0;
997 275 : PetscCall(BVMultVec(W,1.0,0.0,q,coeffs));
998 275 : PetscCall(KSPSolve(ctx->ksp[idxrktg],q,t));
999 275 : PetscCall(VecScale(t,-1.0));
1000 : }
1001 890 : PetscCall(PetscFree(coeffs));
1002 890 : PetscFunctionReturn(PETSC_SUCCESS);
1003 : }
1004 :
1005 : /*
1006 : Compute TOAR coefficients of the blocks of the new Arnoldi vector computed
1007 : */
1008 890 : static PetscErrorCode NEPTOARCoefficients(NEP nep,PetscScalar sigma,PetscInt nv,PetscScalar *S,PetscInt ls,PetscScalar *r,PetscInt lr,PetscScalar *x,PetscScalar *work)
1009 : {
1010 890 : NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1011 890 : PetscInt k,j,d=ctx->nmat-1;
1012 890 : PetscScalar *t=work;
1013 :
1014 890 : PetscFunctionBegin;
1015 890 : PetscCall(NEPNLEIGSEvalNRTFunct(nep,d-1,sigma,t));
1016 8899 : for (k=0;k<d-1;k++) {
1017 239483 : for (j=0;j<=nv;j++) r[k*lr+j] += t[k]*x[j];
1018 : }
1019 23917 : for (j=0;j<=nv;j++) r[(d-1)*lr+j] = t[d-1]*x[j];
1020 890 : PetscFunctionReturn(PETSC_SUCCESS);
1021 : }
1022 :
1023 : /*
1024 : Compute continuation vector coefficients for the Rational-Krylov run.
1025 : dim(work) >= (end-ini)*(end-ini+1) + end+1 + 2*(end-ini+1), dim(t) = end.
1026 : */
1027 890 : static PetscErrorCode NEPNLEIGS_RKcontinuation(NEP nep,PetscInt ini,PetscInt end,PetscScalar *K,PetscScalar *H,PetscInt ld,PetscScalar sigma,PetscScalar *S,PetscInt lds,PetscScalar *cont,PetscScalar *t,PetscScalar *work)
1028 : {
1029 890 : PetscScalar *x,*W,*tau,sone=1.0,szero=0.0;
1030 890 : PetscInt i,j,n1,n,nwu=0;
1031 890 : PetscBLASInt info,n_,n1_,one=1,dim,lds_;
1032 890 : NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
1033 :
1034 890 : PetscFunctionBegin;
1035 890 : if (!ctx->nshifts || !end) {
1036 823 : t[0] = 1;
1037 823 : PetscCall(PetscArraycpy(cont,S+end*lds,lds));
1038 : } else {
1039 67 : n = end-ini;
1040 67 : n1 = n+1;
1041 67 : x = work+nwu;
1042 67 : nwu += end+1;
1043 67 : tau = work+nwu;
1044 67 : nwu += n;
1045 67 : W = work+nwu;
1046 67 : nwu += n1*n;
1047 742 : for (j=ini;j<end;j++) {
1048 9705 : for (i=ini;i<=end;i++) W[(j-ini)*n1+i-ini] = K[j*ld+i] -H[j*ld+i]*sigma;
1049 : }
1050 67 : PetscCall(PetscBLASIntCast(n,&n_));
1051 67 : PetscCall(PetscBLASIntCast(n1,&n1_));
1052 67 : PetscCall(PetscBLASIntCast(end+1,&dim));
1053 67 : PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
1054 67 : PetscCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&n1_,&n_,W,&n1_,tau,work+nwu,&n1_,&info));
1055 67 : SlepcCheckLapackInfo("geqrf",info);
1056 742 : for (i=0;i<end;i++) t[i] = 0.0;
1057 67 : t[end] = 1.0;
1058 742 : for (j=n-1;j>=0;j--) {
1059 4515 : for (i=0;i<ini+j;i++) x[i] = 0.0;
1060 675 : x[ini+j] = 1.0;
1061 5190 : for (i=j+1;i<n1;i++) x[i+ini] = W[i+n1*j];
1062 675 : tau[j] = PetscConj(tau[j]);
1063 675 : PetscCallBLAS("LAPACKlarf",LAPACKlarf_("L",&dim,&one,x,&one,tau+j,t,&dim,work+nwu));
1064 : }
1065 67 : PetscCall(PetscBLASIntCast(lds,&lds_));
1066 67 : PetscCallBLAS("BLASgemv",BLASgemv_("N",&lds_,&n1_,&sone,S,&lds_,t,&one,&szero,cont,&one));
1067 67 : PetscCall(PetscFPTrapPop());
1068 : }
1069 890 : PetscFunctionReturn(PETSC_SUCCESS);
1070 : }
1071 :
1072 : /*
1073 : Compute a run of Arnoldi iterations
1074 : */
1075 72 : static PetscErrorCode NEPNLEIGSTOARrun(NEP nep,Mat MK,Mat MH,BV W,PetscInt k,PetscInt *M,PetscReal *betah,PetscScalar *betak,PetscBool *breakdown,Vec *t_)
1076 : {
1077 72 : NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
1078 72 : PetscInt i,j,m=*M,lwa,deg=ctx->nmat-1,lds,nqt,ld,l,ldh;
1079 72 : Vec t;
1080 72 : PetscReal norm=0.0;
1081 72 : PetscScalar *x,*work,*tt,sigma=1.0,*cont,*S,*K=NULL,*H;
1082 72 : PetscBool lindep;
1083 72 : Mat MS;
1084 :
1085 72 : PetscFunctionBegin;
1086 72 : *betah = 0.0; *betak = 0.0;
1087 72 : PetscCall(MatDenseGetArray(MH,&H));
1088 72 : if (MK) PetscCall(MatDenseGetArray(MK,&K));
1089 72 : PetscCall(MatDenseGetLDA(MH,&ldh));
1090 72 : PetscCall(BVTensorGetFactors(ctx->V,NULL,&MS));
1091 72 : PetscCall(MatDenseGetArray(MS,&S));
1092 72 : PetscCall(BVGetSizes(nep->V,NULL,NULL,&ld));
1093 72 : lds = ld*deg;
1094 72 : PetscCall(BVGetActiveColumns(nep->V,&l,&nqt));
1095 72 : lwa = PetscMax(ld,deg)+(m+1)*(m+1)+4*(m+1);
1096 72 : PetscCall(PetscMalloc4(ld,&x,lwa,&work,m+1,&tt,lds,&cont));
1097 72 : PetscCall(BVSetActiveColumns(ctx->V,0,m));
1098 962 : for (j=k;j<m;j++) {
1099 890 : sigma = ctx->shifts[(++(ctx->idxrk))%ctx->nshiftsw];
1100 :
1101 : /* Continuation vector */
1102 890 : PetscCall(NEPNLEIGS_RKcontinuation(nep,0,j,K,H,ldh,sigma,S,lds,cont,tt,work));
1103 :
1104 : /* apply operator */
1105 890 : PetscCall(BVGetColumn(nep->V,nqt,&t));
1106 890 : PetscCall(NEPTOARExtendBasis(nep,(ctx->idxrk)%ctx->nshiftsw,cont,ld,nqt,W,nep->V,t,S+(j+1)*lds,ld,t_));
1107 890 : PetscCall(BVRestoreColumn(nep->V,nqt,&t));
1108 :
1109 : /* orthogonalize */
1110 890 : PetscCall(BVOrthogonalizeColumn(nep->V,nqt,x,&norm,&lindep));
1111 890 : if (!lindep) {
1112 833 : x[nqt] = norm;
1113 833 : PetscCall(BVScaleColumn(nep->V,nqt,1.0/norm));
1114 833 : nqt++;
1115 57 : } else x[nqt] = 0.0;
1116 :
1117 890 : PetscCall(NEPTOARCoefficients(nep,sigma,nqt-1,cont,ld,S+(j+1)*lds,ld,x,work));
1118 :
1119 : /* Level-2 orthogonalization */
1120 890 : PetscCall(BVOrthogonalizeColumn(ctx->V,j+1,H+j*ldh,&norm,breakdown));
1121 890 : H[j+1+ldh*j] = norm;
1122 890 : if (ctx->nshifts && MK) {
1123 815 : for (i=0;i<=j;i++) K[i+ldh*j] = sigma*H[i+ldh*j] + tt[i];
1124 70 : K[j+1+ldh*j] = sigma*H[j+1+ldh*j];
1125 : }
1126 890 : if (*breakdown) {
1127 0 : *M = j+1;
1128 0 : break;
1129 : }
1130 890 : PetscCall(BVScaleColumn(ctx->V,j+1,1.0/norm));
1131 890 : PetscCall(BVSetActiveColumns(nep->V,l,nqt));
1132 : }
1133 72 : *betah = norm;
1134 72 : if (ctx->nshifts) *betak = norm*sigma;
1135 72 : PetscCall(PetscFree4(x,work,tt,cont));
1136 72 : PetscCall(MatDenseRestoreArray(MS,&S));
1137 72 : PetscCall(MatDenseRestoreArray(MH,&H));
1138 72 : if (MK) PetscCall(MatDenseRestoreArray(MK,&K));
1139 72 : PetscCall(BVTensorRestoreFactors(ctx->V,NULL,&MS));
1140 72 : PetscFunctionReturn(PETSC_SUCCESS);
1141 : }
1142 :
1143 24 : PetscErrorCode NEPSolve_NLEIGS(NEP nep)
1144 : {
1145 24 : NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
1146 24 : PetscInt i,k=0,l,nv=0,ld,lds,nq;
1147 24 : PetscInt deg=ctx->nmat-1,nconv=0,dsn,dsk;
1148 24 : PetscScalar *pU,betak=0,*eigr,*eigi;
1149 24 : const PetscScalar *S;
1150 24 : PetscReal betah;
1151 24 : PetscBool falselock=PETSC_FALSE,breakdown=PETSC_FALSE;
1152 24 : BV W;
1153 24 : Mat H,K=NULL,MS,MQ,U;
1154 :
1155 24 : PetscFunctionBegin;
1156 24 : if (ctx->lock) {
1157 : /* undocumented option to use a cheaper locking instead of the true locking */
1158 21 : PetscCall(PetscOptionsGetBool(NULL,NULL,"-nep_nleigs_falselocking",&falselock,NULL));
1159 : }
1160 :
1161 24 : PetscCall(BVGetSizes(nep->V,NULL,NULL,&ld));
1162 24 : lds = deg*ld;
1163 24 : if (!ctx->nshifts) PetscCall(PetscMalloc2(nep->ncv,&eigr,nep->ncv,&eigi));
1164 3 : else { eigr = nep->eigr; eigi = nep->eigi; }
1165 24 : PetscCall(BVDuplicateResize(nep->V,PetscMax(nep->nt-1,ctx->nmat-1),&W));
1166 :
1167 : /* clean projected matrix (including the extra-arrow) */
1168 24 : PetscCall(DSSetDimensions(nep->ds,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_DETERMINE));
1169 24 : PetscCall(DSGetMat(nep->ds,DS_MAT_A,&H));
1170 24 : PetscCall(MatZeroEntries(H));
1171 24 : PetscCall(DSRestoreMat(nep->ds,DS_MAT_A,&H));
1172 24 : if (ctx->nshifts) {
1173 3 : PetscCall(DSGetMat(nep->ds,DS_MAT_B,&H));
1174 3 : PetscCall(MatZeroEntries(H));
1175 3 : PetscCall(DSRestoreMat(nep->ds,DS_MAT_B,&H));
1176 : }
1177 :
1178 : /* Get the starting Arnoldi vector */
1179 24 : PetscCall(BVTensorBuildFirstColumn(ctx->V,nep->nini));
1180 :
1181 : /* Restart loop */
1182 24 : l = 0;
1183 24 : while (nep->reason == NEP_CONVERGED_ITERATING) {
1184 72 : nep->its++;
1185 :
1186 : /* Compute an nv-step Krylov relation */
1187 72 : nv = PetscMin(nep->nconv+nep->mpd,nep->ncv);
1188 72 : if (ctx->nshifts) PetscCall(DSGetMat(nep->ds,DS_MAT_A,&K));
1189 72 : PetscCall(DSGetMat(nep->ds,ctx->nshifts?DS_MAT_B:DS_MAT_A,&H));
1190 72 : PetscCall(NEPNLEIGSTOARrun(nep,K,H,W,nep->nconv+l,&nv,&betah,&betak,&breakdown,nep->work));
1191 72 : PetscCall(DSRestoreMat(nep->ds,ctx->nshifts?DS_MAT_B:DS_MAT_A,&H));
1192 72 : if (ctx->nshifts) PetscCall(DSRestoreMat(nep->ds,DS_MAT_A,&K));
1193 72 : PetscCall(DSSetDimensions(nep->ds,nv,nep->nconv,nep->nconv+l));
1194 72 : if (l==0) PetscCall(DSSetState(nep->ds,DS_STATE_INTERMEDIATE));
1195 48 : else PetscCall(DSSetState(nep->ds,DS_STATE_RAW));
1196 :
1197 : /* Solve projected problem */
1198 72 : PetscCall(DSSolve(nep->ds,nep->eigr,nep->eigi));
1199 72 : PetscCall(DSSort(nep->ds,nep->eigr,nep->eigi,NULL,NULL,NULL));
1200 72 : PetscCall(DSUpdateExtraRow(nep->ds));
1201 72 : PetscCall(DSSynchronize(nep->ds,nep->eigr,nep->eigi));
1202 :
1203 : /* Check convergence */
1204 72 : PetscCall(NEPNLEIGSKrylovConvergence(nep,PETSC_FALSE,nep->nconv,nv-nep->nconv,betah,betak,&k,nep->work));
1205 72 : PetscCall((*nep->stopping)(nep,nep->its,nep->max_it,k,nep->nev,&nep->reason,nep->stoppingctx));
1206 :
1207 : /* Update l */
1208 72 : if (nep->reason != NEP_CONVERGED_ITERATING || breakdown) l = 0;
1209 : else {
1210 48 : l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
1211 48 : PetscCall(DSGetTruncateSize(nep->ds,k,nv,&l));
1212 48 : if (!breakdown) {
1213 : /* Prepare the Rayleigh quotient for restart */
1214 48 : PetscCall(DSGetDimensions(nep->ds,&dsn,NULL,&dsk,NULL));
1215 48 : PetscCall(DSSetDimensions(nep->ds,dsn,k,dsk));
1216 48 : PetscCall(DSTruncate(nep->ds,k+l,PETSC_FALSE));
1217 : }
1218 : }
1219 72 : nconv = k;
1220 72 : if (!ctx->lock && nep->reason == NEP_CONVERGED_ITERATING && !breakdown) { l += k; k = 0; }
1221 72 : if (l) PetscCall(PetscInfo(nep,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
1222 :
1223 : /* Update S */
1224 139 : PetscCall(DSGetMat(nep->ds,ctx->nshifts?DS_MAT_Z:DS_MAT_Q,&MQ));
1225 72 : PetscCall(BVMultInPlace(ctx->V,MQ,nep->nconv,k+l));
1226 139 : PetscCall(DSRestoreMat(nep->ds,ctx->nshifts?DS_MAT_Z:DS_MAT_Q,&MQ));
1227 :
1228 : /* Copy last column of S */
1229 72 : PetscCall(BVCopyColumn(ctx->V,nv,k+l));
1230 :
1231 72 : if (breakdown && nep->reason == NEP_CONVERGED_ITERATING) {
1232 : /* Stop if breakdown */
1233 0 : PetscCall(PetscInfo(nep,"Breakdown (it=%" PetscInt_FMT " norm=%g)\n",nep->its,(double)betah));
1234 0 : nep->reason = NEP_DIVERGED_BREAKDOWN;
1235 : }
1236 72 : if (nep->reason != NEP_CONVERGED_ITERATING) l--;
1237 : /* truncate S */
1238 72 : PetscCall(BVGetActiveColumns(nep->V,NULL,&nq));
1239 72 : if (k+l+deg<=nq) {
1240 66 : PetscCall(BVSetActiveColumns(ctx->V,nep->nconv,k+l+1));
1241 66 : if (!falselock && ctx->lock) PetscCall(BVTensorCompress(ctx->V,k-nep->nconv));
1242 6 : else PetscCall(BVTensorCompress(ctx->V,0));
1243 : }
1244 72 : nep->nconv = k;
1245 72 : if (!ctx->nshifts) {
1246 1395 : for (i=0;i<nv;i++) { eigr[i] = nep->eigr[i]; eigi[i] = nep->eigi[i]; }
1247 67 : PetscCall(NEPNLEIGSBackTransform((PetscObject)nep,nv,eigr,eigi));
1248 : }
1249 96 : PetscCall(NEPMonitor(nep,nep->its,nconv,eigr,eigi,nep->errest,nv));
1250 : }
1251 24 : nep->nconv = nconv;
1252 24 : if (nep->nconv>0) {
1253 24 : PetscCall(BVSetActiveColumns(ctx->V,0,nep->nconv));
1254 24 : PetscCall(BVGetActiveColumns(nep->V,NULL,&nq));
1255 24 : PetscCall(BVSetActiveColumns(nep->V,0,nq));
1256 24 : if (nq>nep->nconv) {
1257 4 : PetscCall(BVTensorCompress(ctx->V,nep->nconv));
1258 4 : PetscCall(BVSetActiveColumns(nep->V,0,nep->nconv));
1259 4 : nq = nep->nconv;
1260 : }
1261 24 : if (ctx->nshifts) {
1262 3 : PetscCall(DSGetMat(nep->ds,DS_MAT_B,&MQ));
1263 3 : PetscCall(BVMultInPlace(ctx->V,MQ,0,nep->nconv));
1264 3 : PetscCall(DSRestoreMat(nep->ds,DS_MAT_B,&MQ));
1265 : }
1266 24 : PetscCall(BVTensorGetFactors(ctx->V,NULL,&MS));
1267 24 : PetscCall(MatDenseGetArrayRead(MS,&S));
1268 24 : PetscCall(PetscMalloc1(nq*nep->nconv,&pU));
1269 123 : for (i=0;i<nep->nconv;i++) PetscCall(PetscArraycpy(pU+i*nq,S+i*lds,nq));
1270 24 : PetscCall(MatDenseRestoreArrayRead(MS,&S));
1271 24 : PetscCall(BVTensorRestoreFactors(ctx->V,NULL,&MS));
1272 24 : PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,nq,nep->nconv,pU,&U));
1273 24 : PetscCall(BVSetActiveColumns(nep->V,0,nq));
1274 24 : PetscCall(BVMultInPlace(nep->V,U,0,nep->nconv));
1275 24 : PetscCall(BVSetActiveColumns(nep->V,0,nep->nconv));
1276 24 : PetscCall(MatDestroy(&U));
1277 24 : PetscCall(PetscFree(pU));
1278 24 : PetscCall(DSTruncate(nep->ds,nep->nconv,PETSC_TRUE));
1279 : }
1280 :
1281 : /* Map eigenvalues back to the original problem */
1282 24 : if (!ctx->nshifts) {
1283 21 : PetscCall(NEPNLEIGSBackTransform((PetscObject)nep,nep->nconv,nep->eigr,nep->eigi));
1284 21 : PetscCall(PetscFree2(eigr,eigi));
1285 : }
1286 24 : PetscCall(BVDestroy(&W));
1287 24 : PetscFunctionReturn(PETSC_SUCCESS);
1288 : }
1289 :
1290 12 : static PetscErrorCode NEPNLEIGSSetSingularitiesFunction_NLEIGS(NEP nep,PetscErrorCode (*fun)(NEP,PetscInt*,PetscScalar*,void*),void *ctx)
1291 : {
1292 12 : NEP_NLEIGS *nepctx=(NEP_NLEIGS*)nep->data;
1293 :
1294 12 : PetscFunctionBegin;
1295 12 : if (fun) nepctx->computesingularities = fun;
1296 12 : if (ctx) nepctx->singularitiesctx = ctx;
1297 12 : nep->state = NEP_STATE_INITIAL;
1298 12 : PetscFunctionReturn(PETSC_SUCCESS);
1299 : }
1300 :
1301 : /*@C
1302 : NEPNLEIGSSetSingularitiesFunction - Sets a user function to compute a discretization
1303 : of the singularity set (where T(.) is not analytic).
1304 :
1305 : Logically Collective
1306 :
1307 : Input Parameters:
1308 : + nep - the NEP context
1309 : . fun - user function (if NULL then NEP retains any previously set value)
1310 : - ctx - [optional] user-defined context for private data for the function
1311 : (may be NULL, in which case NEP retains any previously set value)
1312 :
1313 : Calling sequence of fun:
1314 : $ PetscErrorCode fun(NEP nep,PetscInt *maxnp,PetscScalar *xi,void *ctx)
1315 : + nep - the NEP context
1316 : . maxnp - on input number of requested points in the discretization (can be set)
1317 : . xi - computed values of the discretization
1318 : - ctx - optional context, as set by NEPNLEIGSSetSingularitiesFunction()
1319 :
1320 : Notes:
1321 : The user-defined function can set a smaller value of maxnp if necessary.
1322 : It is wrong to return a larger value.
1323 :
1324 : If the problem type has been set to rational with NEPSetProblemType(),
1325 : then it is not necessary to set the singularities explicitly since the
1326 : solver will try to determine them automatically.
1327 :
1328 : Level: intermediate
1329 :
1330 : .seealso: NEPNLEIGSGetSingularitiesFunction(), NEPSetProblemType()
1331 : @*/
1332 12 : PetscErrorCode NEPNLEIGSSetSingularitiesFunction(NEP nep,PetscErrorCode (*fun)(NEP nep,PetscInt *maxnp,PetscScalar *xi,void *ctx),void *ctx)
1333 : {
1334 12 : PetscFunctionBegin;
1335 12 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
1336 12 : PetscTryMethod(nep,"NEPNLEIGSSetSingularitiesFunction_C",(NEP,PetscErrorCode(*)(NEP,PetscInt*,PetscScalar*,void*),void*),(nep,fun,ctx));
1337 12 : PetscFunctionReturn(PETSC_SUCCESS);
1338 : }
1339 :
1340 1 : static PetscErrorCode NEPNLEIGSGetSingularitiesFunction_NLEIGS(NEP nep,PetscErrorCode (**fun)(NEP,PetscInt*,PetscScalar*,void*),void **ctx)
1341 : {
1342 1 : NEP_NLEIGS *nepctx=(NEP_NLEIGS*)nep->data;
1343 :
1344 1 : PetscFunctionBegin;
1345 1 : if (fun) *fun = nepctx->computesingularities;
1346 1 : if (ctx) *ctx = nepctx->singularitiesctx;
1347 1 : PetscFunctionReturn(PETSC_SUCCESS);
1348 : }
1349 :
1350 : /*@C
1351 : NEPNLEIGSGetSingularitiesFunction - Returns the Function and optionally the user
1352 : provided context for computing a discretization of the singularity set.
1353 :
1354 : Not Collective
1355 :
1356 : Input Parameter:
1357 : . nep - the nonlinear eigensolver context
1358 :
1359 : Output Parameters:
1360 : + fun - location to put the function (or NULL)
1361 : - ctx - location to stash the function context (or NULL)
1362 :
1363 : Level: advanced
1364 :
1365 : .seealso: NEPNLEIGSSetSingularitiesFunction()
1366 : @*/
1367 1 : PetscErrorCode NEPNLEIGSGetSingularitiesFunction(NEP nep,PetscErrorCode (**fun)(NEP,PetscInt*,PetscScalar*,void*),void **ctx)
1368 : {
1369 1 : PetscFunctionBegin;
1370 1 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
1371 1 : PetscUseMethod(nep,"NEPNLEIGSGetSingularitiesFunction_C",(NEP,PetscErrorCode(**)(NEP,PetscInt*,PetscScalar*,void*),void**),(nep,fun,ctx));
1372 1 : PetscFunctionReturn(PETSC_SUCCESS);
1373 : }
1374 :
1375 3 : static PetscErrorCode NEPNLEIGSSetRestart_NLEIGS(NEP nep,PetscReal keep)
1376 : {
1377 3 : NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1378 :
1379 3 : PetscFunctionBegin;
1380 3 : if (keep==(PetscReal)PETSC_DEFAULT || keep==(PetscReal)PETSC_DECIDE) ctx->keep = 0.5;
1381 : else {
1382 3 : PetscCheck(keep>=0.1 && keep<=0.9,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The keep argument must be in the range [0.1,0.9]");
1383 3 : ctx->keep = keep;
1384 : }
1385 3 : PetscFunctionReturn(PETSC_SUCCESS);
1386 : }
1387 :
1388 : /*@
1389 : NEPNLEIGSSetRestart - Sets the restart parameter for the NLEIGS
1390 : method, in particular the proportion of basis vectors that must be kept
1391 : after restart.
1392 :
1393 : Logically Collective
1394 :
1395 : Input Parameters:
1396 : + nep - the nonlinear eigensolver context
1397 : - keep - the number of vectors to be kept at restart
1398 :
1399 : Options Database Key:
1400 : . -nep_nleigs_restart - Sets the restart parameter
1401 :
1402 : Notes:
1403 : Allowed values are in the range [0.1,0.9]. The default is 0.5.
1404 :
1405 : Level: advanced
1406 :
1407 : .seealso: NEPNLEIGSGetRestart()
1408 : @*/
1409 3 : PetscErrorCode NEPNLEIGSSetRestart(NEP nep,PetscReal keep)
1410 : {
1411 3 : PetscFunctionBegin;
1412 3 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
1413 9 : PetscValidLogicalCollectiveReal(nep,keep,2);
1414 3 : PetscTryMethod(nep,"NEPNLEIGSSetRestart_C",(NEP,PetscReal),(nep,keep));
1415 3 : PetscFunctionReturn(PETSC_SUCCESS);
1416 : }
1417 :
1418 6 : static PetscErrorCode NEPNLEIGSGetRestart_NLEIGS(NEP nep,PetscReal *keep)
1419 : {
1420 6 : NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1421 :
1422 6 : PetscFunctionBegin;
1423 6 : *keep = ctx->keep;
1424 6 : PetscFunctionReturn(PETSC_SUCCESS);
1425 : }
1426 :
1427 : /*@
1428 : NEPNLEIGSGetRestart - Gets the restart parameter used in the NLEIGS method.
1429 :
1430 : Not Collective
1431 :
1432 : Input Parameter:
1433 : . nep - the nonlinear eigensolver context
1434 :
1435 : Output Parameter:
1436 : . keep - the restart parameter
1437 :
1438 : Level: advanced
1439 :
1440 : .seealso: NEPNLEIGSSetRestart()
1441 : @*/
1442 6 : PetscErrorCode NEPNLEIGSGetRestart(NEP nep,PetscReal *keep)
1443 : {
1444 6 : PetscFunctionBegin;
1445 6 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
1446 6 : PetscAssertPointer(keep,2);
1447 6 : PetscUseMethod(nep,"NEPNLEIGSGetRestart_C",(NEP,PetscReal*),(nep,keep));
1448 6 : PetscFunctionReturn(PETSC_SUCCESS);
1449 : }
1450 :
1451 3 : static PetscErrorCode NEPNLEIGSSetLocking_NLEIGS(NEP nep,PetscBool lock)
1452 : {
1453 3 : NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1454 :
1455 3 : PetscFunctionBegin;
1456 3 : ctx->lock = lock;
1457 3 : PetscFunctionReturn(PETSC_SUCCESS);
1458 : }
1459 :
1460 : /*@
1461 : NEPNLEIGSSetLocking - Choose between locking and non-locking variants of
1462 : the NLEIGS method.
1463 :
1464 : Logically Collective
1465 :
1466 : Input Parameters:
1467 : + nep - the nonlinear eigensolver context
1468 : - lock - true if the locking variant must be selected
1469 :
1470 : Options Database Key:
1471 : . -nep_nleigs_locking - Sets the locking flag
1472 :
1473 : Notes:
1474 : The default is to lock converged eigenpairs when the method restarts.
1475 : This behaviour can be changed so that all directions are kept in the
1476 : working subspace even if already converged to working accuracy (the
1477 : non-locking variant).
1478 :
1479 : Level: advanced
1480 :
1481 : .seealso: NEPNLEIGSGetLocking()
1482 : @*/
1483 3 : PetscErrorCode NEPNLEIGSSetLocking(NEP nep,PetscBool lock)
1484 : {
1485 3 : PetscFunctionBegin;
1486 3 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
1487 9 : PetscValidLogicalCollectiveBool(nep,lock,2);
1488 3 : PetscTryMethod(nep,"NEPNLEIGSSetLocking_C",(NEP,PetscBool),(nep,lock));
1489 3 : PetscFunctionReturn(PETSC_SUCCESS);
1490 : }
1491 :
1492 6 : static PetscErrorCode NEPNLEIGSGetLocking_NLEIGS(NEP nep,PetscBool *lock)
1493 : {
1494 6 : NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1495 :
1496 6 : PetscFunctionBegin;
1497 6 : *lock = ctx->lock;
1498 6 : PetscFunctionReturn(PETSC_SUCCESS);
1499 : }
1500 :
1501 : /*@
1502 : NEPNLEIGSGetLocking - Gets the locking flag used in the NLEIGS method.
1503 :
1504 : Not Collective
1505 :
1506 : Input Parameter:
1507 : . nep - the nonlinear eigensolver context
1508 :
1509 : Output Parameter:
1510 : . lock - the locking flag
1511 :
1512 : Level: advanced
1513 :
1514 : .seealso: NEPNLEIGSSetLocking()
1515 : @*/
1516 6 : PetscErrorCode NEPNLEIGSGetLocking(NEP nep,PetscBool *lock)
1517 : {
1518 6 : PetscFunctionBegin;
1519 6 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
1520 6 : PetscAssertPointer(lock,2);
1521 6 : PetscUseMethod(nep,"NEPNLEIGSGetLocking_C",(NEP,PetscBool*),(nep,lock));
1522 6 : PetscFunctionReturn(PETSC_SUCCESS);
1523 : }
1524 :
1525 9 : static PetscErrorCode NEPNLEIGSSetInterpolation_NLEIGS(NEP nep,PetscReal tol,PetscInt degree)
1526 : {
1527 9 : NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1528 :
1529 9 : PetscFunctionBegin;
1530 9 : if (tol == (PetscReal)PETSC_DETERMINE) {
1531 6 : ctx->ddtol = PETSC_DETERMINE;
1532 6 : nep->state = NEP_STATE_INITIAL;
1533 3 : } else if (tol != (PetscReal)PETSC_CURRENT) {
1534 3 : PetscCheck(tol>0.0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of tol. Must be > 0");
1535 3 : ctx->ddtol = tol;
1536 : }
1537 9 : if (degree == PETSC_DETERMINE) {
1538 0 : ctx->ddmaxit = 0;
1539 0 : if (nep->state) PetscCall(NEPReset(nep));
1540 0 : nep->state = NEP_STATE_INITIAL;
1541 9 : } else if (degree != PETSC_CURRENT) {
1542 9 : PetscCheck(degree>0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of degree. Must be > 0");
1543 9 : if (ctx->ddmaxit != degree) {
1544 9 : ctx->ddmaxit = degree;
1545 9 : if (nep->state) PetscCall(NEPReset(nep));
1546 9 : nep->state = NEP_STATE_INITIAL;
1547 : }
1548 : }
1549 9 : PetscFunctionReturn(PETSC_SUCCESS);
1550 : }
1551 :
1552 : /*@
1553 : NEPNLEIGSSetInterpolation - Sets the tolerance and maximum degree
1554 : when building the interpolation via divided differences.
1555 :
1556 : Collective
1557 :
1558 : Input Parameters:
1559 : + nep - the nonlinear eigensolver context
1560 : . tol - tolerance to stop computing divided differences
1561 : - degree - maximum degree of interpolation
1562 :
1563 : Options Database Key:
1564 : + -nep_nleigs_interpolation_tol <tol> - Sets the tolerance to stop computing divided differences
1565 : - -nep_nleigs_interpolation_degree <degree> - Sets the maximum degree of interpolation
1566 :
1567 : Notes:
1568 : PETSC_CURRENT can be used to preserve the current value of any of the
1569 : arguments, and PETSC_DETERMINE to set them to a default value.
1570 :
1571 : Level: advanced
1572 :
1573 : .seealso: NEPNLEIGSGetInterpolation()
1574 : @*/
1575 9 : PetscErrorCode NEPNLEIGSSetInterpolation(NEP nep,PetscReal tol,PetscInt degree)
1576 : {
1577 9 : PetscFunctionBegin;
1578 9 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
1579 27 : PetscValidLogicalCollectiveReal(nep,tol,2);
1580 27 : PetscValidLogicalCollectiveInt(nep,degree,3);
1581 9 : PetscTryMethod(nep,"NEPNLEIGSSetInterpolation_C",(NEP,PetscReal,PetscInt),(nep,tol,degree));
1582 9 : PetscFunctionReturn(PETSC_SUCCESS);
1583 : }
1584 :
1585 31 : static PetscErrorCode NEPNLEIGSGetInterpolation_NLEIGS(NEP nep,PetscReal *tol,PetscInt *degree)
1586 : {
1587 31 : NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1588 :
1589 31 : PetscFunctionBegin;
1590 31 : if (tol) *tol = ctx->ddtol;
1591 31 : if (degree) *degree = ctx->ddmaxit;
1592 31 : PetscFunctionReturn(PETSC_SUCCESS);
1593 : }
1594 :
1595 : /*@
1596 : NEPNLEIGSGetInterpolation - Gets the tolerance and maximum degree
1597 : when building the interpolation via divided differences.
1598 :
1599 : Not Collective
1600 :
1601 : Input Parameter:
1602 : . nep - the nonlinear eigensolver context
1603 :
1604 : Output Parameters:
1605 : + tol - tolerance to stop computing divided differences
1606 : - degree - maximum degree of interpolation
1607 :
1608 : Level: advanced
1609 :
1610 : .seealso: NEPNLEIGSSetInterpolation()
1611 : @*/
1612 31 : PetscErrorCode NEPNLEIGSGetInterpolation(NEP nep,PetscReal *tol,PetscInt *degree)
1613 : {
1614 31 : PetscFunctionBegin;
1615 31 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
1616 31 : PetscTryMethod(nep,"NEPNLEIGSGetInterpolation_C",(NEP,PetscReal*,PetscInt*),(nep,tol,degree));
1617 31 : PetscFunctionReturn(PETSC_SUCCESS);
1618 : }
1619 :
1620 3 : static PetscErrorCode NEPNLEIGSSetRKShifts_NLEIGS(NEP nep,PetscInt ns,PetscScalar *shifts)
1621 : {
1622 3 : NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1623 3 : PetscInt i;
1624 :
1625 3 : PetscFunctionBegin;
1626 3 : PetscCheck(ns>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"Number of shifts must be non-negative");
1627 3 : if (ctx->nshifts) PetscCall(PetscFree(ctx->shifts));
1628 3 : for (i=0;i<ctx->nshiftsw;i++) PetscCall(KSPDestroy(&ctx->ksp[i]));
1629 3 : PetscCall(PetscFree(ctx->ksp));
1630 3 : ctx->ksp = NULL;
1631 3 : if (ns) {
1632 3 : PetscCall(PetscMalloc1(ns,&ctx->shifts));
1633 15 : for (i=0;i<ns;i++) ctx->shifts[i] = shifts[i];
1634 : }
1635 3 : ctx->nshifts = ns;
1636 3 : nep->state = NEP_STATE_INITIAL;
1637 3 : PetscFunctionReturn(PETSC_SUCCESS);
1638 : }
1639 :
1640 : /*@
1641 : NEPNLEIGSSetRKShifts - Sets a list of shifts to be used in the Rational
1642 : Krylov method.
1643 :
1644 : Collective
1645 :
1646 : Input Parameters:
1647 : + nep - the nonlinear eigensolver context
1648 : . ns - number of shifts
1649 : - shifts - array of scalar values specifying the shifts
1650 :
1651 : Options Database Key:
1652 : . -nep_nleigs_rk_shifts - Sets the list of shifts
1653 :
1654 : Notes:
1655 : If only one shift is provided, the built subspace built is equivalent to
1656 : shift-and-invert Krylov-Schur (provided that the absolute convergence
1657 : criterion is used).
1658 :
1659 : In the case of real scalars, complex shifts are not allowed. In the
1660 : command line, a comma-separated list of complex values can be provided with
1661 : the format [+/-][realnumber][+/-]realnumberi with no spaces, e.g.
1662 : -nep_nleigs_rk_shifts 1.0+2.0i,1.5+2.0i,1.0+1.5i
1663 :
1664 : Use ns=0 to remove previously set shifts.
1665 :
1666 : Level: advanced
1667 :
1668 : .seealso: NEPNLEIGSGetRKShifts()
1669 : @*/
1670 3 : PetscErrorCode NEPNLEIGSSetRKShifts(NEP nep,PetscInt ns,PetscScalar shifts[])
1671 : {
1672 3 : PetscFunctionBegin;
1673 3 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
1674 9 : PetscValidLogicalCollectiveInt(nep,ns,2);
1675 3 : if (ns) PetscAssertPointer(shifts,3);
1676 3 : PetscTryMethod(nep,"NEPNLEIGSSetRKShifts_C",(NEP,PetscInt,PetscScalar*),(nep,ns,shifts));
1677 3 : PetscFunctionReturn(PETSC_SUCCESS);
1678 : }
1679 :
1680 1 : static PetscErrorCode NEPNLEIGSGetRKShifts_NLEIGS(NEP nep,PetscInt *ns,PetscScalar **shifts)
1681 : {
1682 1 : NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1683 1 : PetscInt i;
1684 :
1685 1 : PetscFunctionBegin;
1686 1 : *ns = ctx->nshifts;
1687 1 : if (ctx->nshifts) {
1688 1 : PetscCall(PetscMalloc1(ctx->nshifts,shifts));
1689 5 : for (i=0;i<ctx->nshifts;i++) (*shifts)[i] = ctx->shifts[i];
1690 : }
1691 1 : PetscFunctionReturn(PETSC_SUCCESS);
1692 : }
1693 :
1694 : /*@C
1695 : NEPNLEIGSGetRKShifts - Gets the list of shifts used in the Rational
1696 : Krylov method.
1697 :
1698 : Not Collective
1699 :
1700 : Input Parameter:
1701 : . nep - the nonlinear eigensolver context
1702 :
1703 : Output Parameters:
1704 : + ns - number of shifts
1705 : - shifts - array of shifts
1706 :
1707 : Note:
1708 : The user is responsible for deallocating the returned array.
1709 :
1710 : Level: advanced
1711 :
1712 : .seealso: NEPNLEIGSSetRKShifts()
1713 : @*/
1714 1 : PetscErrorCode NEPNLEIGSGetRKShifts(NEP nep,PetscInt *ns,PetscScalar *shifts[])
1715 : {
1716 1 : PetscFunctionBegin;
1717 1 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
1718 1 : PetscAssertPointer(ns,2);
1719 1 : PetscAssertPointer(shifts,3);
1720 1 : PetscTryMethod(nep,"NEPNLEIGSGetRKShifts_C",(NEP,PetscInt*,PetscScalar**),(nep,ns,shifts));
1721 1 : PetscFunctionReturn(PETSC_SUCCESS);
1722 : }
1723 :
1724 25 : static PetscErrorCode NEPNLEIGSGetKSPs_NLEIGS(NEP nep,PetscInt *nsolve,KSP **ksp)
1725 : {
1726 25 : NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
1727 25 : PetscInt i;
1728 25 : PC pc;
1729 :
1730 25 : PetscFunctionBegin;
1731 25 : if (!ctx->ksp) {
1732 25 : PetscCall(NEPNLEIGSSetShifts(nep,&ctx->nshiftsw));
1733 25 : PetscCall(PetscMalloc1(ctx->nshiftsw,&ctx->ksp));
1734 59 : for (i=0;i<ctx->nshiftsw;i++) {
1735 34 : PetscCall(KSPCreate(PetscObjectComm((PetscObject)nep),&ctx->ksp[i]));
1736 34 : PetscCall(PetscObjectIncrementTabLevel((PetscObject)ctx->ksp[i],(PetscObject)nep,1));
1737 34 : PetscCall(KSPSetOptionsPrefix(ctx->ksp[i],((PetscObject)nep)->prefix));
1738 34 : PetscCall(KSPAppendOptionsPrefix(ctx->ksp[i],"nep_nleigs_"));
1739 34 : PetscCall(PetscObjectSetOptions((PetscObject)ctx->ksp[i],((PetscObject)nep)->options));
1740 34 : PetscCall(KSPSetErrorIfNotConverged(ctx->ksp[i],PETSC_TRUE));
1741 57 : PetscCall(KSPSetTolerances(ctx->ksp[i],1e-3*SlepcDefaultTol(nep->tol),PETSC_CURRENT,PETSC_CURRENT,PETSC_CURRENT));
1742 34 : PetscCall(KSPGetPC(ctx->ksp[i],&pc));
1743 34 : if ((nep->fui==NEP_USER_INTERFACE_SPLIT && nep->P) || (nep->fui==NEP_USER_INTERFACE_CALLBACK && nep->function_pre!=nep->function)) {
1744 2 : PetscCall(KSPSetType(ctx->ksp[i],KSPBCGS));
1745 2 : PetscCall(PCSetType(pc,PCBJACOBI));
1746 : } else {
1747 32 : PetscCall(KSPSetType(ctx->ksp[i],KSPPREONLY));
1748 34 : PetscCall(PCSetType(pc,PCLU));
1749 : }
1750 : }
1751 : }
1752 25 : if (nsolve) *nsolve = ctx->nshiftsw;
1753 25 : if (ksp) *ksp = ctx->ksp;
1754 25 : PetscFunctionReturn(PETSC_SUCCESS);
1755 : }
1756 :
1757 : /*@C
1758 : NEPNLEIGSGetKSPs - Retrieve the array of linear solver objects associated with
1759 : the nonlinear eigenvalue solver.
1760 :
1761 : Collective
1762 :
1763 : Input Parameter:
1764 : . nep - nonlinear eigenvalue solver
1765 :
1766 : Output Parameters:
1767 : + nsolve - number of returned KSP objects
1768 : - ksp - array of linear solver object
1769 :
1770 : Notes:
1771 : The number of KSP objects is equal to the number of shifts provided by the user,
1772 : or 1 if the user did not provide shifts.
1773 :
1774 : Level: advanced
1775 :
1776 : .seealso: NEPNLEIGSSetRKShifts()
1777 : @*/
1778 25 : PetscErrorCode NEPNLEIGSGetKSPs(NEP nep,PetscInt *nsolve,KSP **ksp)
1779 : {
1780 25 : PetscFunctionBegin;
1781 25 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
1782 25 : PetscUseMethod(nep,"NEPNLEIGSGetKSPs_C",(NEP,PetscInt*,KSP**),(nep,nsolve,ksp));
1783 25 : PetscFunctionReturn(PETSC_SUCCESS);
1784 : }
1785 :
1786 4 : static PetscErrorCode NEPNLEIGSSetFullBasis_NLEIGS(NEP nep,PetscBool fullbasis)
1787 : {
1788 4 : NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1789 :
1790 4 : PetscFunctionBegin;
1791 4 : if (fullbasis!=ctx->fullbasis) {
1792 3 : ctx->fullbasis = fullbasis;
1793 3 : nep->state = NEP_STATE_INITIAL;
1794 3 : nep->useds = PetscNot(fullbasis);
1795 : }
1796 4 : PetscFunctionReturn(PETSC_SUCCESS);
1797 : }
1798 :
1799 : /*@
1800 : NEPNLEIGSSetFullBasis - Choose between TOAR-basis (default) and full-basis
1801 : variants of the NLEIGS method.
1802 :
1803 : Logically Collective
1804 :
1805 : Input Parameters:
1806 : + nep - the nonlinear eigensolver context
1807 : - fullbasis - true if the full-basis variant must be selected
1808 :
1809 : Options Database Key:
1810 : . -nep_nleigs_full_basis - Sets the full-basis flag
1811 :
1812 : Notes:
1813 : The default is to use a compact representation of the Krylov basis, that is,
1814 : V = (I otimes U) S, with a tensor BV. This behaviour can be changed so that
1815 : the full basis V is explicitly stored and operated with. This variant is more
1816 : expensive in terms of memory and computation, but is necessary in some cases,
1817 : particularly for two-sided computations, see NEPSetTwoSided().
1818 :
1819 : In the full-basis variant, the NLEIGS solver uses an EPS object to explicitly
1820 : solve the linearized eigenproblem, see NEPNLEIGSGetEPS().
1821 :
1822 : Level: advanced
1823 :
1824 : .seealso: NEPNLEIGSGetFullBasis(), NEPNLEIGSGetEPS(), NEPSetTwoSided(), BVCreateTensor()
1825 : @*/
1826 4 : PetscErrorCode NEPNLEIGSSetFullBasis(NEP nep,PetscBool fullbasis)
1827 : {
1828 4 : PetscFunctionBegin;
1829 4 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
1830 12 : PetscValidLogicalCollectiveBool(nep,fullbasis,2);
1831 4 : PetscTryMethod(nep,"NEPNLEIGSSetFullBasis_C",(NEP,PetscBool),(nep,fullbasis));
1832 4 : PetscFunctionReturn(PETSC_SUCCESS);
1833 : }
1834 :
1835 1 : static PetscErrorCode NEPNLEIGSGetFullBasis_NLEIGS(NEP nep,PetscBool *fullbasis)
1836 : {
1837 1 : NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1838 :
1839 1 : PetscFunctionBegin;
1840 1 : *fullbasis = ctx->fullbasis;
1841 1 : PetscFunctionReturn(PETSC_SUCCESS);
1842 : }
1843 :
1844 : /*@
1845 : NEPNLEIGSGetFullBasis - Gets the flag that indicates if NLEIGS is using the
1846 : full-basis variant.
1847 :
1848 : Not Collective
1849 :
1850 : Input Parameter:
1851 : . nep - the nonlinear eigensolver context
1852 :
1853 : Output Parameter:
1854 : . fullbasis - the flag
1855 :
1856 : Level: advanced
1857 :
1858 : .seealso: NEPNLEIGSSetFullBasis()
1859 : @*/
1860 1 : PetscErrorCode NEPNLEIGSGetFullBasis(NEP nep,PetscBool *fullbasis)
1861 : {
1862 1 : PetscFunctionBegin;
1863 1 : PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
1864 1 : PetscAssertPointer(fullbasis,2);
1865 1 : PetscUseMethod(nep,"NEPNLEIGSGetFullBasis_C",(NEP,PetscBool*),(nep,fullbasis));
1866 1 : PetscFunctionReturn(PETSC_SUCCESS);
1867 : }
1868 :
1869 : #define SHIFTMAX 30
1870 :
1871 25 : static PetscErrorCode NEPSetFromOptions_NLEIGS(NEP nep,PetscOptionItems *PetscOptionsObject)
1872 : {
1873 25 : NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
1874 25 : PetscInt i=0,k;
1875 25 : PetscBool flg1,flg2,b;
1876 25 : PetscReal r;
1877 25 : PetscScalar array[SHIFTMAX];
1878 :
1879 25 : PetscFunctionBegin;
1880 25 : PetscOptionsHeadBegin(PetscOptionsObject,"NEP NLEIGS Options");
1881 :
1882 25 : PetscCall(PetscOptionsReal("-nep_nleigs_restart","Proportion of vectors kept after restart","NEPNLEIGSSetRestart",0.5,&r,&flg1));
1883 25 : if (flg1) PetscCall(NEPNLEIGSSetRestart(nep,r));
1884 :
1885 25 : PetscCall(PetscOptionsBool("-nep_nleigs_locking","Choose between locking and non-locking variants","NEPNLEIGSSetLocking",PETSC_FALSE,&b,&flg1));
1886 25 : if (flg1) PetscCall(NEPNLEIGSSetLocking(nep,b));
1887 :
1888 25 : PetscCall(PetscOptionsBool("-nep_nleigs_full_basis","Choose between TOAR and full-basis variants","NEPNLEIGSSetFullBasis",PETSC_FALSE,&b,&flg1));
1889 25 : if (flg1) PetscCall(NEPNLEIGSSetFullBasis(nep,b));
1890 :
1891 25 : PetscCall(NEPNLEIGSGetInterpolation(nep,&r,&i));
1892 25 : if (!i) i = PETSC_DETERMINE;
1893 25 : PetscCall(PetscOptionsInt("-nep_nleigs_interpolation_degree","Maximum number of terms for interpolation via divided differences","NEPNLEIGSSetInterpolation",i,&i,&flg1));
1894 25 : PetscCall(PetscOptionsReal("-nep_nleigs_interpolation_tol","Tolerance for interpolation via divided differences","NEPNLEIGSSetInterpolation",r,&r,&flg2));
1895 25 : if (flg1 || flg2) PetscCall(NEPNLEIGSSetInterpolation(nep,r,i));
1896 :
1897 25 : k = SHIFTMAX;
1898 775 : for (i=0;i<k;i++) array[i] = 0;
1899 25 : PetscCall(PetscOptionsScalarArray("-nep_nleigs_rk_shifts","Shifts for Rational Krylov","NEPNLEIGSSetRKShifts",array,&k,&flg1));
1900 25 : if (flg1) PetscCall(NEPNLEIGSSetRKShifts(nep,k,array));
1901 :
1902 25 : PetscOptionsHeadEnd();
1903 :
1904 25 : if (!ctx->ksp) PetscCall(NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp));
1905 59 : for (i=0;i<ctx->nshiftsw;i++) PetscCall(KSPSetFromOptions(ctx->ksp[i]));
1906 :
1907 25 : if (ctx->fullbasis) {
1908 3 : if (!ctx->eps) PetscCall(NEPNLEIGSGetEPS(nep,&ctx->eps));
1909 3 : PetscCall(EPSSetFromOptions(ctx->eps));
1910 : }
1911 25 : PetscFunctionReturn(PETSC_SUCCESS);
1912 : }
1913 :
1914 1 : static PetscErrorCode NEPView_NLEIGS(NEP nep,PetscViewer viewer)
1915 : {
1916 1 : NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1917 1 : PetscBool isascii;
1918 1 : PetscInt i;
1919 1 : char str[50];
1920 :
1921 1 : PetscFunctionBegin;
1922 1 : PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
1923 1 : if (isascii) {
1924 1 : PetscCall(PetscViewerASCIIPrintf(viewer," %d%% of basis vectors kept after restart\n",(int)(100*ctx->keep)));
1925 1 : if (ctx->fullbasis) PetscCall(PetscViewerASCIIPrintf(viewer," using the full-basis variant\n"));
1926 1 : else PetscCall(PetscViewerASCIIPrintf(viewer," using the %slocking variant\n",ctx->lock?"":"non-"));
1927 1 : PetscCall(PetscViewerASCIIPrintf(viewer," divided difference terms: used=%" PetscInt_FMT ", max=%" PetscInt_FMT "\n",ctx->nmat,ctx->ddmaxit));
1928 1 : PetscCall(PetscViewerASCIIPrintf(viewer," tolerance for divided difference convergence: %g\n",(double)ctx->ddtol));
1929 1 : if (ctx->nshifts) {
1930 1 : PetscCall(PetscViewerASCIIPrintf(viewer," RK shifts: "));
1931 1 : PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE));
1932 5 : for (i=0;i<ctx->nshifts;i++) {
1933 4 : PetscCall(SlepcSNPrintfScalar(str,sizeof(str),ctx->shifts[i],PETSC_FALSE));
1934 5 : PetscCall(PetscViewerASCIIPrintf(viewer,"%s%s",str,(i<ctx->nshifts-1)?",":""));
1935 : }
1936 1 : PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
1937 1 : PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE));
1938 : }
1939 1 : if (!ctx->ksp) PetscCall(NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp));
1940 1 : PetscCall(PetscViewerASCIIPushTab(viewer));
1941 1 : PetscCall(KSPView(ctx->ksp[0],viewer));
1942 1 : PetscCall(PetscViewerASCIIPopTab(viewer));
1943 1 : if (ctx->fullbasis) {
1944 0 : if (!ctx->eps) PetscCall(NEPNLEIGSGetEPS(nep,&ctx->eps));
1945 0 : PetscCall(PetscViewerASCIIPushTab(viewer));
1946 0 : PetscCall(EPSView(ctx->eps,viewer));
1947 0 : PetscCall(PetscViewerASCIIPopTab(viewer));
1948 : }
1949 : }
1950 1 : PetscFunctionReturn(PETSC_SUCCESS);
1951 : }
1952 :
1953 27 : static PetscErrorCode NEPReset_NLEIGS(NEP nep)
1954 : {
1955 27 : PetscInt k;
1956 27 : NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1957 :
1958 27 : PetscFunctionBegin;
1959 27 : if (nep->fui==NEP_USER_INTERFACE_SPLIT) PetscCall(PetscFree(ctx->coeffD));
1960 : else {
1961 152 : for (k=0;k<ctx->nmat;k++) PetscCall(MatDestroy(&ctx->D[k]));
1962 : }
1963 27 : PetscCall(PetscFree4(ctx->s,ctx->xi,ctx->beta,ctx->D));
1964 63 : for (k=0;k<ctx->nshiftsw;k++) PetscCall(KSPReset(ctx->ksp[k]));
1965 27 : PetscCall(VecDestroy(&ctx->vrn));
1966 27 : if (ctx->fullbasis) {
1967 3 : PetscCall(MatDestroy(&ctx->A));
1968 3 : PetscCall(EPSReset(ctx->eps));
1969 15 : for (k=0;k<4;k++) PetscCall(VecDestroy(&ctx->w[k]));
1970 : }
1971 27 : PetscFunctionReturn(PETSC_SUCCESS);
1972 : }
1973 :
1974 25 : static PetscErrorCode NEPDestroy_NLEIGS(NEP nep)
1975 : {
1976 25 : PetscInt k;
1977 25 : NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
1978 :
1979 25 : PetscFunctionBegin;
1980 25 : PetscCall(BVDestroy(&ctx->V));
1981 59 : for (k=0;k<ctx->nshiftsw;k++) PetscCall(KSPDestroy(&ctx->ksp[k]));
1982 25 : PetscCall(PetscFree(ctx->ksp));
1983 25 : if (ctx->nshifts) PetscCall(PetscFree(ctx->shifts));
1984 25 : if (ctx->fullbasis) PetscCall(EPSDestroy(&ctx->eps));
1985 25 : PetscCall(PetscFree(nep->data));
1986 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetSingularitiesFunction_C",NULL));
1987 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetSingularitiesFunction_C",NULL));
1988 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRestart_C",NULL));
1989 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRestart_C",NULL));
1990 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetLocking_C",NULL));
1991 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetLocking_C",NULL));
1992 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetInterpolation_C",NULL));
1993 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetInterpolation_C",NULL));
1994 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRKShifts_C",NULL));
1995 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRKShifts_C",NULL));
1996 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetKSPs_C",NULL));
1997 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetFullBasis_C",NULL));
1998 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetFullBasis_C",NULL));
1999 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetEPS_C",NULL));
2000 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetEPS_C",NULL));
2001 25 : PetscFunctionReturn(PETSC_SUCCESS);
2002 : }
2003 :
2004 25 : SLEPC_EXTERN PetscErrorCode NEPCreate_NLEIGS(NEP nep)
2005 : {
2006 25 : NEP_NLEIGS *ctx;
2007 :
2008 25 : PetscFunctionBegin;
2009 25 : PetscCall(PetscNew(&ctx));
2010 25 : nep->data = (void*)ctx;
2011 25 : ctx->lock = PETSC_TRUE;
2012 25 : ctx->ddtol = PETSC_DETERMINE;
2013 :
2014 25 : nep->useds = PETSC_TRUE;
2015 :
2016 25 : nep->ops->setup = NEPSetUp_NLEIGS;
2017 25 : nep->ops->setfromoptions = NEPSetFromOptions_NLEIGS;
2018 25 : nep->ops->view = NEPView_NLEIGS;
2019 25 : nep->ops->destroy = NEPDestroy_NLEIGS;
2020 25 : nep->ops->reset = NEPReset_NLEIGS;
2021 :
2022 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetSingularitiesFunction_C",NEPNLEIGSSetSingularitiesFunction_NLEIGS));
2023 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetSingularitiesFunction_C",NEPNLEIGSGetSingularitiesFunction_NLEIGS));
2024 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRestart_C",NEPNLEIGSSetRestart_NLEIGS));
2025 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRestart_C",NEPNLEIGSGetRestart_NLEIGS));
2026 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetLocking_C",NEPNLEIGSSetLocking_NLEIGS));
2027 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetLocking_C",NEPNLEIGSGetLocking_NLEIGS));
2028 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetInterpolation_C",NEPNLEIGSSetInterpolation_NLEIGS));
2029 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetInterpolation_C",NEPNLEIGSGetInterpolation_NLEIGS));
2030 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRKShifts_C",NEPNLEIGSSetRKShifts_NLEIGS));
2031 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRKShifts_C",NEPNLEIGSGetRKShifts_NLEIGS));
2032 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetKSPs_C",NEPNLEIGSGetKSPs_NLEIGS));
2033 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetFullBasis_C",NEPNLEIGSSetFullBasis_NLEIGS));
2034 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetFullBasis_C",NEPNLEIGSGetFullBasis_NLEIGS));
2035 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetEPS_C",NEPNLEIGSSetEPS_NLEIGS));
2036 25 : PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetEPS_C",NEPNLEIGSGetEPS_NLEIGS));
2037 25 : PetscFunctionReturn(PETSC_SUCCESS);
2038 : }
|
