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