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 !PetscDefined(USE_COMPLEX)
35: PetscScalar t;
36: #endif
38: PetscFunctionBegin;
39: nep = (NEP)ob;
40: #if !PetscDefined(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 PetscDefined(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 !PetscDefined(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 !PetscDefined(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,one=1;
207: #if PetscDefined(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 PetscDefined(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 PetscDefined(USE_COMPLEX)
248: PetscCallLAPACKInfo("LAPACKgesvd",LAPACKgesvd_("N","A",&m_,&n_,A,&lda_,S,NULL,&lda_,VT,&lda_,work,&lwork,rwork,&info));
249: #else
250: PetscCallLAPACKInfo("LAPACKgesvd",LAPACKgesvd_("N","A",&m_,&n_,A,&lda_,S,NULL,&lda_,VT,&lda_,work,&lwork,&info));
251: #endif
252: for (i=0;i<=k;i++) {
253: ww[i] = PetscConj(VT[i*ndpt+k]);
254: D[i] = ww[i]*f[i];
255: }
256: PetscCallBLAS("BLASgemv",BLASgemv_("N",&lda_,&n_,&sone,C,&lda_,D,&one,&szero,N,&one));
257: PetscCallBLAS("BLASgemv",BLASgemv_("N",&lda_,&n_,&sone,C,&lda_,ww,&one,&szero,D,&one));
258: for (i=0;i<ndpt;i++) if (R[i]>=0) R[i] = PetscAbsScalar(F[i]-N[i]/D[i]);
259: /* next support point */
260: err = 0.0;
261: for (i=0;i<ndpt;i++) if (R[i]>=err) {idx = i; err = R[i];}
262: if (err <= ctx->ddtol*norm) break;
263: }
265: PetscCheck(k<ndpt-1,PetscObjectComm((PetscObject)nep),PETSC_ERR_CONV_FAILED,"Failed to determine singularities automatically in general problem");
266: /* poles */
267: PetscCall(PetscArrayzero(C,ndpt*ndpt));
268: PetscCall(PetscArrayzero(A,ndpt*ndpt));
269: for (i=0;i<=k;i++) {
270: C[i+ndpt*i] = 1.0;
271: A[(i+1)*ndpt] = ww[i];
272: A[i+1] = 1.0;
273: A[i+1+(i+1)*ndpt] = z[i];
274: }
275: C[0] = 0.0; C[k+1+(k+1)*ndpt] = 1.0;
276: n_++;
277: #if PetscDefined(USE_COMPLEX)
278: PetscCallLAPACKInfo("LAPACK" LAPGEEV,LAPACKggevalt_("N","N",&n_,A,&lda_,C,&lda_,D,N,NULL,&lda_,NULL,&lda_,work,&lwork,rwork,&info));
279: #else
280: PetscCallLAPACKInfo("LAPACK" LAPGEEV,LAPACKggevalt_("N","N",&n_,A,&lda_,C,&lda_,D,VT,N,NULL,&lda_,NULL,&lda_,work,&lwork,&info));
281: #endif
282: cont = 0.0;
283: for (i=0;i<n_;i++) if (N[i]!=0.0) {
284: dxi[cont++] = D[i]/N[i];
285: }
286: *ndptx = cont;
287: PetscCall(PetscFPTrapPop());
288: PetscCall(PetscFree5(R,z,f,C,ww));
289: PetscCall(PetscFree6(A,S,VT,work,D,N));
290: #if PetscDefined(USE_COMPLEX)
291: PetscCall(PetscFree(rwork));
292: #endif
293: PetscFunctionReturn(PETSC_SUCCESS);
294: }
296: /* Singularities using Adaptive Anderson-Antoulas algorithm */
297: static PetscErrorCode NEPNLEIGSAAASingularities(NEP nep,PetscInt ndpt,PetscScalar *ds,PetscInt *ndptx,PetscScalar *dxi)
298: {
299: Vec u,v,w;
300: PetscRandom rand=NULL;
301: PetscScalar *F,*isol;
302: PetscInt i,k,nisol,nt;
303: Mat T;
304: FN f;
306: PetscFunctionBegin;
307: PetscCall(PetscMalloc1(ndpt,&F));
308: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
309: PetscCall(PetscMalloc1(ndpt,&isol));
310: *ndptx = 0;
311: PetscCall(NEPGetSplitOperatorInfo(nep,&nt,NULL));
312: nisol = *ndptx;
313: for (k=0;k<nt;k++) {
314: PetscCall(NEPGetSplitOperatorTerm(nep,k,NULL,&f));
315: for (i=0;i<ndpt;i++) PetscCall(FNEvaluateFunction(f,ds[i],&F[i]));
316: PetscCall(NEPNLEIGSAAAComputation(nep,ndpt,ds,F,&nisol,isol));
317: if (nisol) PetscCall(NEPNLEIGSAuxiliarRmDuplicates(nisol,isol,ndptx,dxi,ndpt));
318: }
319: PetscCall(PetscFree(isol));
320: } else {
321: PetscCall(MatCreateVecs(nep->function,&u,NULL));
322: PetscCall(VecDuplicate(u,&v));
323: PetscCall(VecDuplicate(u,&w));
324: if (nep->V) PetscCall(BVGetRandomContext(nep->V,&rand));
325: PetscCall(VecSetRandom(u,rand));
326: PetscCall(VecNormalize(u,NULL));
327: PetscCall(VecSetRandom(v,rand));
328: PetscCall(VecNormalize(v,NULL));
329: T = nep->function;
330: for (i=0;i<ndpt;i++) {
331: PetscCall(NEPComputeFunction(nep,ds[i],T,T));
332: PetscCall(MatMult(T,v,w));
333: PetscCall(VecDot(w,u,&F[i]));
334: }
335: PetscCall(NEPNLEIGSAAAComputation(nep,ndpt,ds,F,ndptx,dxi));
336: PetscCall(VecDestroy(&u));
337: PetscCall(VecDestroy(&v));
338: PetscCall(VecDestroy(&w));
339: }
340: PetscCall(PetscFree(F));
341: PetscFunctionReturn(PETSC_SUCCESS);
342: }
344: static PetscErrorCode NEPNLEIGSLejaBagbyPoints(NEP nep)
345: {
346: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
347: PetscInt i,k,ndpt=NDPOINTS,ndptx=NDPOINTS;
348: PetscScalar *ds,*dsi,*dxi,*nrs,*nrxi,*s=ctx->s,*xi=ctx->xi,*beta=ctx->beta;
349: PetscReal maxnrs,minnrxi;
350: PetscBool rational;
351: #if !PetscDefined(USE_COMPLEX)
352: PetscReal a,b,h;
353: #endif
355: PetscFunctionBegin;
356: if (!ctx->computesingularities && nep->problem_type!=NEP_RATIONAL) ndpt = ndptx = LBPOINTS;
357: PetscCall(PetscMalloc5(ndpt+1,&ds,ndpt+1,&dsi,ndpt,&dxi,ndpt+1,&nrs,ndpt,&nrxi));
359: /* Discretize the target region boundary */
360: PetscCall(RGComputeContour(nep->rg,ndpt,ds,dsi));
361: #if !PetscDefined(USE_COMPLEX)
362: for (i=0;i<ndpt;i++) if (dsi[i]!=0.0) break;
363: if (i<ndpt) {
364: 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");
365: /* Select a segment in the real axis */
366: PetscCall(RGComputeBoundingBox(nep->rg,&a,&b,NULL,NULL));
367: PetscCheck(a>-PETSC_MAX_REAL && b<PETSC_MAX_REAL,PetscObjectComm((PetscObject)nep),PETSC_ERR_USER_INPUT,"NLEIGS requires a bounded target set");
368: h = (b-a)/ndpt;
369: for (i=0;i<ndpt;i++) {ds[i] = a+h*i; dsi[i] = 0.0;}
370: }
371: #endif
372: /* Discretize the singularity region */
373: if (ctx->computesingularities) PetscCall(ctx->computesingularities(nep,&ndptx,dxi,ctx->singularitiesctx));
374: else {
375: if (nep->problem_type==NEP_RATIONAL) {
376: PetscCall(NEPNLEIGSRationalSingularities(nep,&ndptx,dxi,&rational));
377: PetscCheck(rational,PetscObjectComm((PetscObject)nep),PETSC_ERR_CONV_FAILED,"Failed to determine singularities automatically in rational problem; consider solving the problem as general");
378: } else {
379: /* AAA algorithm */
380: PetscCall(NEPNLEIGSAAASingularities(nep,ndpt,ds,&ndptx,dxi));
381: }
382: }
383: /* Look for Leja-Bagby points in the discretization sets */
384: s[0] = ds[0];
385: xi[0] = (ndptx>0)?dxi[0]:PETSC_INFINITY;
386: 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]));
387: beta[0] = 1.0; /* scaling factors are also computed here */
388: for (i=0;i<ndpt;i++) {
389: nrs[i] = 1.0;
390: nrxi[i] = 1.0;
391: }
392: for (k=1;k<ctx->ddmaxit;k++) {
393: maxnrs = 0.0;
394: minnrxi = PETSC_MAX_REAL;
395: for (i=0;i<ndpt;i++) {
396: nrs[i] *= ((ds[i]-s[k-1])/(1.0-ds[i]/xi[k-1]))/beta[k-1];
397: if (PetscAbsScalar(nrs[i])>maxnrs) {maxnrs = PetscAbsScalar(nrs[i]); s[k] = ds[i];}
398: }
399: if (ndptx>k) {
400: for (i=1;i<ndptx;i++) {
401: nrxi[i] *= ((dxi[i]-s[k-1])/(1.0-dxi[i]/xi[k-1]))/beta[k-1];
402: if (PetscAbsScalar(nrxi[i])<minnrxi) {minnrxi = PetscAbsScalar(nrxi[i]); xi[k] = dxi[i];}
403: }
404: 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]));
405: } else xi[k] = PETSC_INFINITY;
406: beta[k] = maxnrs;
407: }
408: PetscCall(PetscFree5(ds,dsi,dxi,nrs,nrxi));
409: PetscFunctionReturn(PETSC_SUCCESS);
410: }
412: PetscErrorCode NEPNLEIGSEvalNRTFunct(NEP nep,PetscInt k,PetscScalar sigma,PetscScalar *b)
413: {
414: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
415: PetscInt i;
416: PetscScalar *beta=ctx->beta,*s=ctx->s,*xi=ctx->xi;
418: PetscFunctionBegin;
419: b[0] = 1.0/beta[0];
420: for (i=0;i<k;i++) {
421: b[i+1] = ((sigma-s[i])*b[i])/(beta[i+1]*(1.0-sigma/xi[i]));
422: }
423: PetscFunctionReturn(PETSC_SUCCESS);
424: }
426: static PetscErrorCode MatMult_Fun(Mat A,Vec x,Vec y)
427: {
428: NEP_NLEIGS_MATSHELL *ctx;
429: PetscInt i;
431: PetscFunctionBeginUser;
432: PetscCall(MatShellGetContext(A,&ctx));
433: PetscCall(MatMult(ctx->A[0],x,y));
434: if (ctx->coeff[0]!=1.0) PetscCall(VecScale(y,ctx->coeff[0]));
435: for (i=1;i<ctx->nmat;i++) {
436: PetscCall(MatMult(ctx->A[i],x,ctx->t));
437: PetscCall(VecAXPY(y,ctx->coeff[i],ctx->t));
438: }
439: PetscFunctionReturn(PETSC_SUCCESS);
440: }
442: static PetscErrorCode MatMultTranspose_Fun(Mat A,Vec x,Vec y)
443: {
444: NEP_NLEIGS_MATSHELL *ctx;
445: PetscInt i;
447: PetscFunctionBeginUser;
448: PetscCall(MatShellGetContext(A,&ctx));
449: PetscCall(MatMultTranspose(ctx->A[0],x,y));
450: if (ctx->coeff[0]!=1.0) PetscCall(VecScale(y,ctx->coeff[0]));
451: for (i=1;i<ctx->nmat;i++) {
452: PetscCall(MatMultTranspose(ctx->A[i],x,ctx->t));
453: PetscCall(VecAXPY(y,ctx->coeff[i],ctx->t));
454: }
455: PetscFunctionReturn(PETSC_SUCCESS);
456: }
458: static PetscErrorCode MatGetDiagonal_Fun(Mat A,Vec diag)
459: {
460: NEP_NLEIGS_MATSHELL *ctx;
461: PetscInt i;
463: PetscFunctionBeginUser;
464: PetscCall(MatShellGetContext(A,&ctx));
465: PetscCall(MatGetDiagonal(ctx->A[0],diag));
466: if (ctx->coeff[0]!=1.0) PetscCall(VecScale(diag,ctx->coeff[0]));
467: for (i=1;i<ctx->nmat;i++) {
468: PetscCall(MatGetDiagonal(ctx->A[i],ctx->t));
469: PetscCall(VecAXPY(diag,ctx->coeff[i],ctx->t));
470: }
471: PetscFunctionReturn(PETSC_SUCCESS);
472: }
474: static PetscErrorCode MatDuplicate_Fun(Mat A,MatDuplicateOption op,Mat *B)
475: {
476: PetscInt m,n,M,N,i;
477: NEP_NLEIGS_MATSHELL *ctxnew,*ctx;
478: PetscErrorCodeFn *fun;
480: PetscFunctionBeginUser;
481: PetscCall(MatShellGetContext(A,&ctx));
482: PetscCall(PetscNew(&ctxnew));
483: ctxnew->nmat = ctx->nmat;
484: ctxnew->maxnmat = ctx->maxnmat;
485: PetscCall(PetscMalloc2(ctxnew->maxnmat,&ctxnew->A,ctxnew->maxnmat,&ctxnew->coeff));
486: for (i=0;i<ctx->nmat;i++) {
487: PetscCall(PetscObjectReference((PetscObject)ctx->A[i]));
488: ctxnew->A[i] = ctx->A[i];
489: ctxnew->coeff[i] = ctx->coeff[i];
490: }
491: PetscCall(MatGetSize(ctx->A[0],&M,&N));
492: PetscCall(MatGetLocalSize(ctx->A[0],&m,&n));
493: PetscCall(VecDuplicate(ctx->t,&ctxnew->t));
494: PetscCall(MatCreateShell(PetscObjectComm((PetscObject)A),m,n,M,N,(void*)ctxnew,B));
495: PetscCall(MatShellSetManageScalingShifts(*B));
496: PetscCall(MatShellGetOperation(A,MATOP_MULT,&fun));
497: PetscCall(MatShellSetOperation(*B,MATOP_MULT,fun));
498: PetscCall(MatShellGetOperation(A,MATOP_MULT_TRANSPOSE,&fun));
499: PetscCall(MatShellSetOperation(*B,MATOP_MULT_TRANSPOSE,fun));
500: PetscCall(MatShellGetOperation(A,MATOP_GET_DIAGONAL,&fun));
501: PetscCall(MatShellSetOperation(*B,MATOP_GET_DIAGONAL,fun));
502: PetscCall(MatShellGetOperation(A,MATOP_DUPLICATE,&fun));
503: PetscCall(MatShellSetOperation(*B,MATOP_DUPLICATE,fun));
504: PetscCall(MatShellGetOperation(A,MATOP_DESTROY,&fun));
505: PetscCall(MatShellSetOperation(*B,MATOP_DESTROY,fun));
506: PetscCall(MatShellGetOperation(A,MATOP_AXPY,&fun));
507: PetscCall(MatShellSetOperation(*B,MATOP_AXPY,fun));
508: PetscFunctionReturn(PETSC_SUCCESS);
509: }
511: static PetscErrorCode MatDestroy_Fun(Mat A)
512: {
513: NEP_NLEIGS_MATSHELL *ctx;
514: PetscInt i;
516: PetscFunctionBeginUser;
517: if (A) {
518: PetscCall(MatShellGetContext(A,&ctx));
519: for (i=0;i<ctx->nmat;i++) PetscCall(MatDestroy(&ctx->A[i]));
520: PetscCall(VecDestroy(&ctx->t));
521: PetscCall(PetscFree2(ctx->A,ctx->coeff));
522: PetscCall(PetscFree(ctx));
523: }
524: PetscFunctionReturn(PETSC_SUCCESS);
525: }
527: static PetscErrorCode MatAXPY_Fun(Mat Y,PetscScalar a,Mat X,MatStructure str)
528: {
529: NEP_NLEIGS_MATSHELL *ctxY,*ctxX;
530: PetscInt i,j;
531: PetscBool found,has;
533: PetscFunctionBeginUser;
534: PetscCall(MatShellGetContext(Y,&ctxY));
535: PetscCall(MatShellGetContext(X,&ctxX));
536: for (i=0;i<ctxX->nmat;i++) {
537: found = PETSC_FALSE;
538: for (j=0;!found&&j<ctxY->nmat;j++) {
539: if (ctxX->A[i]==ctxY->A[j]) {
540: found = PETSC_TRUE;
541: ctxY->coeff[j] += a*ctxX->coeff[i];
542: }
543: }
544: if (!found) {
545: ctxY->coeff[ctxY->nmat] = a*ctxX->coeff[i];
546: ctxY->A[ctxY->nmat] = ctxX->A[i];
547: PetscCall(MatHasOperation(ctxX->A[i],MATOP_MULT_TRANSPOSE,&has));
548: if (!has) PetscCall(MatShellSetOperation(ctxY->A[ctxY->nmat],MATOP_MULT_TRANSPOSE,NULL));
549: PetscCall(MatHasOperation(ctxX->A[i],MATOP_GET_DIAGONAL,&has));
550: if (!has) PetscCall(MatShellSetOperation(ctxY->A[ctxY->nmat],MATOP_GET_DIAGONAL,NULL));
551: ctxY->nmat++;
552: PetscCall(PetscObjectReference((PetscObject)ctxX->A[i]));
553: }
554: }
555: PetscFunctionReturn(PETSC_SUCCESS);
556: }
558: static PetscErrorCode MatScale_Fun(Mat M,PetscScalar a)
559: {
560: NEP_NLEIGS_MATSHELL *ctx;
561: PetscInt i;
563: PetscFunctionBeginUser;
564: PetscCall(MatShellGetContext(M,&ctx));
565: for (i=0;i<ctx->nmat;i++) ctx->coeff[i] *= a;
566: PetscFunctionReturn(PETSC_SUCCESS);
567: }
569: static PetscErrorCode NLEIGSMatToMatShellArray(Mat A,Mat *Ms,PetscInt maxnmat)
570: {
571: NEP_NLEIGS_MATSHELL *ctx;
572: PetscInt m,n,M,N;
573: PetscBool has;
575: PetscFunctionBegin;
576: PetscCall(MatHasOperation(A,MATOP_DUPLICATE,&has));
577: PetscCheck(has,PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"MatDuplicate operation required");
578: PetscCall(PetscNew(&ctx));
579: ctx->maxnmat = maxnmat;
580: PetscCall(PetscMalloc2(ctx->maxnmat,&ctx->A,ctx->maxnmat,&ctx->coeff));
581: PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&ctx->A[0]));
582: ctx->nmat = 1;
583: ctx->coeff[0] = 1.0;
584: PetscCall(MatCreateVecs(A,&ctx->t,NULL));
585: PetscCall(MatGetSize(A,&M,&N));
586: PetscCall(MatGetLocalSize(A,&m,&n));
587: PetscCall(MatCreateShell(PetscObjectComm((PetscObject)A),m,n,M,N,(void*)ctx,Ms));
588: PetscCall(MatShellSetManageScalingShifts(*Ms));
589: PetscCall(MatShellSetOperation(*Ms,MATOP_MULT,(PetscErrorCodeFn*)MatMult_Fun));
590: PetscCall(MatHasOperation(A,MATOP_MULT_TRANSPOSE,&has));
591: if (has) PetscCall(MatShellSetOperation(*Ms,MATOP_MULT_TRANSPOSE,(PetscErrorCodeFn*)MatMultTranspose_Fun));
592: PetscCall(MatHasOperation(A,MATOP_GET_DIAGONAL,&has));
593: if (has) PetscCall(MatShellSetOperation(*Ms,MATOP_GET_DIAGONAL,(PetscErrorCodeFn*)MatGetDiagonal_Fun));
594: PetscCall(MatShellSetOperation(*Ms,MATOP_DUPLICATE,(PetscErrorCodeFn*)MatDuplicate_Fun));
595: PetscCall(MatShellSetOperation(*Ms,MATOP_DESTROY,(PetscErrorCodeFn*)MatDestroy_Fun));
596: PetscCall(MatShellSetOperation(*Ms,MATOP_AXPY,(PetscErrorCodeFn*)MatAXPY_Fun));
597: PetscCall(MatShellSetOperation(*Ms,MATOP_SCALE,(PetscErrorCodeFn*)MatScale_Fun));
598: PetscFunctionReturn(PETSC_SUCCESS);
599: }
601: /*
602: MatIsShellAny - returns true if any of the n matrices is a shell matrix
603: */
604: static PetscErrorCode MatIsShellAny(Mat *A,PetscInt n,PetscBool *shell)
605: {
606: PetscInt i;
607: PetscBool flg;
609: PetscFunctionBegin;
610: *shell = PETSC_FALSE;
611: for (i=0;i<n;i++) {
612: PetscCall(MatIsShell(A[i],&flg));
613: if (flg) { *shell = PETSC_TRUE; break; }
614: }
615: PetscFunctionReturn(PETSC_SUCCESS);
616: }
618: static PetscErrorCode NEPNLEIGSDividedDifferences_split(NEP nep)
619: {
620: PetscErrorCode ierr;
621: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
622: PetscInt k,j,i,maxnmat,nmax;
623: PetscReal norm0,norm,*matnorm;
624: PetscScalar *s=ctx->s,*beta=ctx->beta,*xi=ctx->xi,*b,alpha,*coeffs,*pK,*pH,sone=1.0;
625: Mat T,P,Ts,K,H;
626: PetscBool shell,hasmnorm=PETSC_FALSE,matrix=PETSC_TRUE;
627: PetscBLASInt n_;
629: PetscFunctionBegin;
630: nmax = ctx->ddmaxit;
631: PetscCall(PetscMalloc1(nep->nt*nmax,&ctx->coeffD));
632: PetscCall(PetscMalloc3(nmax+1,&b,nmax+1,&coeffs,nep->nt,&matnorm));
633: for (j=0;j<nep->nt;j++) {
634: PetscCall(MatHasOperation(nep->A[j],MATOP_NORM,&hasmnorm));
635: if (!hasmnorm) break;
636: PetscCall(MatNorm(nep->A[j],NORM_INFINITY,matnorm+j));
637: }
638: /* Try matrix functions scheme */
639: PetscCall(PetscCalloc2(nmax*nmax,&pK,nmax*nmax,&pH));
640: for (i=0;i<nmax-1;i++) {
641: pK[(nmax+1)*i] = 1.0;
642: pK[(nmax+1)*i+1] = beta[i+1]/xi[i];
643: pH[(nmax+1)*i] = s[i];
644: pH[(nmax+1)*i+1] = beta[i+1];
645: }
646: pH[nmax*nmax-1] = s[nmax-1];
647: pK[nmax*nmax-1] = 1.0;
648: PetscCall(PetscBLASIntCast(nmax,&n_));
649: PetscCallBLAS("BLAStrsm",BLAStrsm_("R","L","N","U",&n_,&n_,&sone,pK,&n_,pH,&n_));
650: /* The matrix to be used is in H. K will be a work-space matrix */
651: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,nmax,nmax,pH,&H));
652: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,nmax,nmax,pK,&K));
653: for (j=0;matrix&&j<nep->nt;j++) {
654: PetscCall(PetscPushErrorHandler(PetscReturnErrorHandler,NULL));
655: ierr = FNEvaluateFunctionMat(nep->f[j],H,K);
656: PetscCall(PetscPopErrorHandler());
657: if (!ierr) {
658: for (i=0;i<nmax;i++) ctx->coeffD[j+i*nep->nt] = pK[i]*beta[0];
659: } else {
660: matrix = PETSC_FALSE;
661: PetscCall(PetscFPTrapPop());
662: }
663: }
664: PetscCall(MatDestroy(&H));
665: PetscCall(MatDestroy(&K));
666: if (!matrix) {
667: for (j=0;j<nep->nt;j++) {
668: PetscCall(FNEvaluateFunction(nep->f[j],s[0],ctx->coeffD+j));
669: ctx->coeffD[j] *= beta[0];
670: }
671: }
672: if (hasmnorm) {
673: norm0 = 0.0;
674: for (j=0;j<nep->nt;j++) norm0 += matnorm[j]*PetscAbsScalar(ctx->coeffD[j]);
675: } else {
676: norm0 = 0.0;
677: for (j=0;j<nep->nt;j++) norm0 = PetscMax(PetscAbsScalar(ctx->coeffD[j]),norm0);
678: }
679: ctx->nmat = ctx->ddmaxit;
680: for (k=1;k<ctx->ddmaxit;k++) {
681: if (!matrix) {
682: PetscCall(NEPNLEIGSEvalNRTFunct(nep,k,s[k],b));
683: for (i=0;i<nep->nt;i++) {
684: PetscCall(FNEvaluateFunction(nep->f[i],s[k],ctx->coeffD+k*nep->nt+i));
685: for (j=0;j<k;j++) {
686: ctx->coeffD[k*nep->nt+i] -= b[j]*ctx->coeffD[i+nep->nt*j];
687: }
688: ctx->coeffD[k*nep->nt+i] /= b[k];
689: }
690: }
691: if (hasmnorm) {
692: norm = 0.0;
693: for (j=0;j<nep->nt;j++) norm += matnorm[j]*PetscAbsScalar(ctx->coeffD[k*nep->nt+j]);
694: } else {
695: norm = 0.0;
696: for (j=0;j<nep->nt;j++) norm = PetscMax(PetscAbsScalar(ctx->coeffD[k*nep->nt+j]),norm);
697: }
698: if (k>1 && norm/norm0 < ctx->ddtol) {
699: ctx->nmat = k+1;
700: break;
701: }
702: }
703: if (!ctx->ksp) PetscCall(NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp));
704: PetscCall(MatIsShellAny(nep->A,nep->nt,&shell));
705: maxnmat = PetscMax(ctx->ddmaxit,nep->nt);
706: for (i=0;i<ctx->nshiftsw;i++) {
707: PetscCall(NEPNLEIGSEvalNRTFunct(nep,ctx->nmat-1,ctx->shifts[i],coeffs));
708: if (!shell) PetscCall(MatDuplicate(nep->A[0],MAT_COPY_VALUES,&T));
709: else PetscCall(NLEIGSMatToMatShellArray(nep->A[0],&T,maxnmat));
710: if (nep->P) { /* user-defined preconditioner */
711: PetscCall(MatDuplicate(nep->P[0],MAT_COPY_VALUES,&P));
712: } else P=T;
713: alpha = 0.0;
714: for (j=0;j<ctx->nmat;j++) alpha += coeffs[j]*ctx->coeffD[j*nep->nt];
715: PetscCall(MatScale(T,alpha));
716: if (nep->P) PetscCall(MatScale(P,alpha));
717: for (k=1;k<nep->nt;k++) {
718: alpha = 0.0;
719: for (j=0;j<ctx->nmat;j++) alpha += coeffs[j]*ctx->coeffD[j*nep->nt+k];
720: if (shell) PetscCall(NLEIGSMatToMatShellArray(nep->A[k],&Ts,maxnmat));
721: PetscCall(MatAXPY(T,alpha,shell?Ts:nep->A[k],nep->mstr));
722: if (nep->P) PetscCall(MatAXPY(P,alpha,nep->P[k],nep->mstrp));
723: if (shell) PetscCall(MatDestroy(&Ts));
724: }
725: PetscCall(NEP_KSPSetOperators(ctx->ksp[i],T,P));
726: PetscCall(KSPSetUp(ctx->ksp[i]));
727: PetscCall(MatDestroy(&T));
728: if (nep->P) PetscCall(MatDestroy(&P));
729: }
730: PetscCall(PetscFree3(b,coeffs,matnorm));
731: PetscCall(PetscFree2(pK,pH));
732: PetscFunctionReturn(PETSC_SUCCESS);
733: }
735: static PetscErrorCode NEPNLEIGSDividedDifferences_callback(NEP nep)
736: {
737: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
738: PetscInt k,j,i,maxnmat;
739: PetscReal norm0,norm;
740: PetscScalar *s=ctx->s,*beta=ctx->beta,*b,*coeffs;
741: Mat *D=ctx->D,*DP,T,P;
742: PetscBool shell,has,precond=(nep->function_pre!=nep->function)?PETSC_TRUE:PETSC_FALSE;
743: PetscRandom rand=NULL;
745: PetscFunctionBegin;
746: PetscCall(PetscMalloc2(ctx->ddmaxit+1,&b,ctx->ddmaxit+1,&coeffs));
747: if (nep->V) PetscCall(BVGetRandomContext(nep->V,&rand));
748: T = nep->function;
749: P = nep->function_pre;
750: PetscCall(NEPComputeFunction(nep,s[0],T,P));
751: PetscCall(MatIsShell(T,&shell));
752: maxnmat = PetscMax(ctx->ddmaxit,nep->nt);
753: if (!shell) PetscCall(MatDuplicate(T,MAT_COPY_VALUES,&D[0]));
754: else PetscCall(NLEIGSMatToMatShellArray(T,&D[0],maxnmat));
755: if (beta[0]!=1.0) PetscCall(MatScale(D[0],1.0/beta[0]));
756: PetscCall(MatHasOperation(D[0],MATOP_NORM,&has));
757: if (has) PetscCall(MatNorm(D[0],NORM_FROBENIUS,&norm0));
758: else PetscCall(MatNormApproximate(D[0],NORM_2,1,&norm0));
759: if (precond) {
760: PetscCall(PetscMalloc1(ctx->ddmaxit,&DP));
761: PetscCall(MatDuplicate(P,MAT_COPY_VALUES,&DP[0]));
762: }
763: ctx->nmat = ctx->ddmaxit;
764: for (k=1;k<ctx->ddmaxit;k++) {
765: PetscCall(NEPNLEIGSEvalNRTFunct(nep,k,s[k],b));
766: PetscCall(NEPComputeFunction(nep,s[k],T,P));
767: if (!shell) PetscCall(MatDuplicate(T,MAT_COPY_VALUES,&D[k]));
768: else PetscCall(NLEIGSMatToMatShellArray(T,&D[k],maxnmat));
769: for (j=0;j<k;j++) PetscCall(MatAXPY(D[k],-b[j],D[j],nep->mstr));
770: PetscCall(MatScale(D[k],1.0/b[k]));
771: PetscCall(MatHasOperation(D[k],MATOP_NORM,&has));
772: if (has) PetscCall(MatNorm(D[k],NORM_FROBENIUS,&norm));
773: else PetscCall(MatNormApproximate(D[k],NORM_2,1,&norm));
774: if (precond) {
775: PetscCall(MatDuplicate(P,MAT_COPY_VALUES,&DP[k]));
776: for (j=0;j<k;j++) PetscCall(MatAXPY(DP[k],-b[j],DP[j],nep->mstrp));
777: PetscCall(MatScale(DP[k],1.0/b[k]));
778: }
779: if (k>1 && norm/norm0 < ctx->ddtol && k>1) {
780: ctx->nmat = k+1;
781: break;
782: }
783: }
784: if (!ctx->ksp) PetscCall(NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp));
785: for (i=0;i<ctx->nshiftsw;i++) {
786: PetscCall(NEPNLEIGSEvalNRTFunct(nep,ctx->nmat-1,ctx->shifts[i],coeffs));
787: PetscCall(MatDuplicate(D[0],MAT_COPY_VALUES,&T));
788: if (coeffs[0]!=1.0) PetscCall(MatScale(T,coeffs[0]));
789: for (j=1;j<ctx->nmat;j++) PetscCall(MatAXPY(T,coeffs[j],D[j],nep->mstr));
790: if (precond) {
791: PetscCall(MatDuplicate(DP[0],MAT_COPY_VALUES,&P));
792: if (coeffs[0]!=1.0) PetscCall(MatScale(P,coeffs[0]));
793: for (j=1;j<ctx->nmat;j++) PetscCall(MatAXPY(P,coeffs[j],DP[j],nep->mstrp));
794: } else P=T;
795: PetscCall(NEP_KSPSetOperators(ctx->ksp[i],T,P));
796: PetscCall(KSPSetUp(ctx->ksp[i]));
797: PetscCall(MatDestroy(&T));
798: }
799: PetscCall(PetscFree2(b,coeffs));
800: if (precond) {
801: PetscCall(MatDestroy(&P));
802: PetscCall(MatDestroyMatrices(ctx->nmat,&DP));
803: }
804: PetscFunctionReturn(PETSC_SUCCESS);
805: }
807: /*
808: NEPKrylovConvergence - This is the analogue to EPSKrylovConvergence.
809: */
810: static PetscErrorCode NEPNLEIGSKrylovConvergence(NEP nep,PetscBool getall,PetscInt kini,PetscInt nits,PetscReal betah,PetscScalar betak,PetscInt *kout,Vec *w)
811: {
812: PetscInt k,newk,marker,inside;
813: PetscScalar re,im;
814: PetscReal resnorm,tt;
815: PetscBool istrivial;
816: NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
818: PetscFunctionBegin;
819: PetscCall(RGIsTrivial(nep->rg,&istrivial));
820: marker = -1;
821: if (nep->trackall) getall = PETSC_TRUE;
822: for (k=kini;k<kini+nits;k++) {
823: /* eigenvalue */
824: re = nep->eigr[k];
825: im = nep->eigi[k];
826: if (!istrivial) {
827: if (!ctx->nshifts) PetscCall(NEPNLEIGSBackTransform((PetscObject)nep,1,&re,&im));
828: PetscCall(RGCheckInside(nep->rg,1,&re,&im,&inside));
829: if (marker==-1 && inside<0) marker = k;
830: }
831: newk = k;
832: PetscCall(DSVectors(nep->ds,DS_MAT_X,&newk,&resnorm));
833: tt = ctx->nshifts?SlepcAbsEigenvalue(betak-nep->eigr[k]*betah,nep->eigi[k]*betah):betah;
834: resnorm *= PetscAbsReal(tt);
835: /* error estimate */
836: PetscCall((*nep->converged)(nep,nep->eigr[k],nep->eigi[k],resnorm,&nep->errest[k],nep->convergedctx));
837: if (marker==-1 && nep->errest[k] >= nep->tol) marker = k;
838: if (newk==k+1) {
839: nep->errest[k+1] = nep->errest[k];
840: k++;
841: }
842: if (marker!=-1 && !getall) break;
843: }
844: if (marker!=-1) k = marker;
845: *kout = k;
846: PetscFunctionReturn(PETSC_SUCCESS);
847: }
849: static PetscErrorCode NEPSetUp_NLEIGS(NEP nep)
850: {
851: PetscInt k,in;
852: PetscScalar zero=0.0;
853: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
854: SlepcSC sc;
855: PetscBool istrivial;
857: PetscFunctionBegin;
858: PetscCall(NEPSetDimensions_Default(nep,nep->nev,&nep->ncv,&nep->mpd));
859: 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");
860: if (nep->max_it==PETSC_DETERMINE) nep->max_it = PetscMax(5000,2*nep->n/nep->ncv);
861: if (!ctx->ddmaxit) ctx->ddmaxit = LBPOINTS;
862: PetscCall(RGIsTrivial(nep->rg,&istrivial));
863: PetscCheck(!istrivial,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"NEPNLEIGS requires a nontrivial region defining the target set");
864: if (!nep->which) nep->which = NEP_TARGET_MAGNITUDE;
865: 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");
867: /* Initialize the NLEIGS context structure */
868: k = ctx->ddmaxit;
869: PetscCall(PetscMalloc4(k,&ctx->s,k,&ctx->xi,k,&ctx->beta,k,&ctx->D));
870: nep->data = ctx;
871: if (nep->tol==(PetscReal)PETSC_DETERMINE) nep->tol = SLEPC_DEFAULT_TOL;
872: if (ctx->ddtol==(PetscReal)PETSC_DETERMINE) ctx->ddtol = nep->tol/10.0;
873: if (!ctx->keep) ctx->keep = 0.5;
875: /* Compute Leja-Bagby points and scaling values */
876: PetscCall(NEPNLEIGSLejaBagbyPoints(nep));
877: if (nep->problem_type!=NEP_RATIONAL) {
878: PetscCall(RGCheckInside(nep->rg,1,&nep->target,&zero,&in));
879: PetscCheck(in>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"The target is not inside the target set");
880: }
882: /* Compute the divided difference matrices */
883: if (nep->fui==NEP_USER_INTERFACE_SPLIT) PetscCall(NEPNLEIGSDividedDifferences_split(nep));
884: else PetscCall(NEPNLEIGSDividedDifferences_callback(nep));
885: PetscCall(NEPAllocateSolution(nep,ctx->nmat-1));
886: PetscCall(NEPSetWorkVecs(nep,4));
887: if (!ctx->fullbasis) {
888: PetscCheck(!nep->twosided,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Two-sided variant requires the full-basis option, rerun with -nep_nleigs_full_basis");
889: /* set-up DS and transfer split operator functions */
890: PetscCall(DSSetType(nep->ds,ctx->nshifts?DSGNHEP:DSNHEP));
891: PetscCall(DSAllocate(nep->ds,nep->ncv+1));
892: PetscCall(DSGetSlepcSC(nep->ds,&sc));
893: if (!ctx->nshifts) sc->map = NEPNLEIGSBackTransform;
894: PetscCall(DSSetExtraRow(nep->ds,PETSC_TRUE));
895: sc->mapobj = (PetscObject)nep;
896: sc->rg = nep->rg;
897: sc->comparison = nep->sc->comparison;
898: sc->comparisonctx = nep->sc->comparisonctx;
899: PetscCall(BVDestroy(&ctx->V));
900: PetscCall(BVCreateTensor(nep->V,ctx->nmat-1,&ctx->V));
901: nep->ops->solve = NEPSolve_NLEIGS;
902: nep->ops->computevectors = NEPComputeVectors_Schur;
903: } else {
904: PetscCall(NEPSetUp_NLEIGS_FullBasis(nep));
905: nep->ops->solve = NEPSolve_NLEIGS_FullBasis;
906: nep->ops->computevectors = NULL;
907: }
908: PetscFunctionReturn(PETSC_SUCCESS);
909: }
911: /*
912: Extend the TOAR basis by applying the matrix operator
913: over a vector which is decomposed on the TOAR way
914: Input:
915: - S,V: define the latest Arnoldi vector (nv vectors in V)
916: Output:
917: - t: new vector extending the TOAR basis
918: - r: temporally coefficients to compute the TOAR coefficients
919: for the new Arnoldi vector
920: Workspace: t_ (two vectors)
921: */
922: 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_)
923: {
924: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
925: PetscInt deg=ctx->nmat-1,k,j;
926: Vec v=t_[0],q=t_[1],w;
927: PetscScalar *beta=ctx->beta,*s=ctx->s,*xi=ctx->xi,*coeffs,sigma;
929: PetscFunctionBegin;
930: if (!ctx->ksp) PetscCall(NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp));
931: sigma = ctx->shifts[idxrktg];
932: PetscCall(BVSetActiveColumns(nep->V,0,nv));
933: PetscCall(PetscMalloc1(ctx->nmat,&coeffs));
934: PetscCheck(PetscAbsScalar(s[deg-2]-sigma)>100*PETSC_MACHINE_EPSILON,PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"Breakdown in NLEIGS");
935: /* i-part stored in (i-1) position */
936: for (j=0;j<nv;j++) {
937: 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);
938: }
939: PetscCall(BVSetActiveColumns(W,0,deg));
940: PetscCall(BVGetColumn(W,deg-1,&w));
941: PetscCall(BVMultVec(V,1.0/beta[deg],0,w,S+(deg-1)*ls));
942: PetscCall(BVRestoreColumn(W,deg-1,&w));
943: PetscCall(BVGetColumn(W,deg-2,&w));
944: PetscCall(BVMultVec(V,1.0,0.0,w,r+(deg-2)*lr));
945: PetscCall(BVRestoreColumn(W,deg-2,&w));
946: for (k=deg-2;k>0;k--) {
947: PetscCheck(PetscAbsScalar(s[k-1]-sigma)>100*PETSC_MACHINE_EPSILON,PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"Breakdown in NLEIGS");
948: 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);
949: PetscCall(BVGetColumn(W,k-1,&w));
950: PetscCall(BVMultVec(V,1.0,0.0,w,r+(k-1)*lr));
951: PetscCall(BVRestoreColumn(W,k-1,&w));
952: }
953: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
954: for (j=0;j<ctx->nmat-2;j++) coeffs[j] = ctx->coeffD[nep->nt*j];
955: coeffs[ctx->nmat-2] = ctx->coeffD[nep->nt*(ctx->nmat-1)];
956: PetscCall(BVMultVec(W,1.0,0.0,v,coeffs));
957: PetscCall(MatMult(nep->A[0],v,q));
958: for (k=1;k<nep->nt;k++) {
959: for (j=0;j<ctx->nmat-2;j++) coeffs[j] = ctx->coeffD[nep->nt*j+k];
960: coeffs[ctx->nmat-2] = ctx->coeffD[nep->nt*(ctx->nmat-1)+k];
961: PetscCall(BVMultVec(W,1.0,0,v,coeffs));
962: PetscCall(MatMult(nep->A[k],v,t));
963: PetscCall(VecAXPY(q,1.0,t));
964: }
965: PetscCall(KSPSolve(ctx->ksp[idxrktg],q,t));
966: PetscCall(VecScale(t,-1.0));
967: } else {
968: for (k=0;k<deg-1;k++) {
969: PetscCall(BVGetColumn(W,k,&w));
970: PetscCall(MatMult(ctx->D[k],w,q));
971: PetscCall(BVRestoreColumn(W,k,&w));
972: PetscCall(BVInsertVec(W,k,q));
973: }
974: PetscCall(BVGetColumn(W,deg-1,&w));
975: PetscCall(MatMult(ctx->D[deg],w,q));
976: PetscCall(BVRestoreColumn(W,k,&w));
977: PetscCall(BVInsertVec(W,k,q));
978: for (j=0;j<ctx->nmat-1;j++) coeffs[j] = 1.0;
979: PetscCall(BVMultVec(W,1.0,0.0,q,coeffs));
980: PetscCall(KSPSolve(ctx->ksp[idxrktg],q,t));
981: PetscCall(VecScale(t,-1.0));
982: }
983: PetscCall(PetscFree(coeffs));
984: PetscFunctionReturn(PETSC_SUCCESS);
985: }
987: /*
988: Compute TOAR coefficients of the blocks of the new Arnoldi vector computed
989: */
990: static PetscErrorCode NEPTOARCoefficients(NEP nep,PetscScalar sigma,PetscInt nv,PetscScalar *S,PetscInt ls,PetscScalar *r,PetscInt lr,PetscScalar *x,PetscScalar *work)
991: {
992: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
993: PetscInt k,j,d=ctx->nmat-1;
994: PetscScalar *t=work;
996: PetscFunctionBegin;
997: PetscCall(NEPNLEIGSEvalNRTFunct(nep,d-1,sigma,t));
998: for (k=0;k<d-1;k++) {
999: for (j=0;j<=nv;j++) r[k*lr+j] += t[k]*x[j];
1000: }
1001: for (j=0;j<=nv;j++) r[(d-1)*lr+j] = t[d-1]*x[j];
1002: PetscFunctionReturn(PETSC_SUCCESS);
1003: }
1005: /*
1006: Compute continuation vector coefficients for the Rational-Krylov run.
1007: dim(work) >= (end-ini)*(end-ini+1) + end+1 + 2*(end-ini+1), dim(t) = end.
1008: */
1009: 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)
1010: {
1011: PetscScalar *x,*W,*tau,sone=1.0,szero=0.0;
1012: PetscInt i,j,n1,n,nwu=0;
1013: PetscBLASInt n_,n1_,one=1,dim,lds_;
1014: NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
1016: PetscFunctionBegin;
1017: if (!ctx->nshifts || !end) {
1018: t[0] = 1;
1019: PetscCall(PetscArraycpy(cont,S+end*lds,lds));
1020: } else {
1021: n = end-ini;
1022: n1 = n+1;
1023: x = work+nwu;
1024: nwu += end+1;
1025: tau = work+nwu;
1026: nwu += n;
1027: W = work+nwu;
1028: nwu += n1*n;
1029: for (j=ini;j<end;j++) {
1030: for (i=ini;i<=end;i++) W[(j-ini)*n1+i-ini] = K[j*ld+i] -H[j*ld+i]*sigma;
1031: }
1032: PetscCall(PetscBLASIntCast(n,&n_));
1033: PetscCall(PetscBLASIntCast(n1,&n1_));
1034: PetscCall(PetscBLASIntCast(end+1,&dim));
1035: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
1036: PetscCallLAPACKInfo("LAPACKgeqrf",LAPACKgeqrf_(&n1_,&n_,W,&n1_,tau,work+nwu,&n1_,&info));
1037: for (i=0;i<end;i++) t[i] = 0.0;
1038: t[end] = 1.0;
1039: for (j=n-1;j>=0;j--) {
1040: for (i=0;i<ini+j;i++) x[i] = 0.0;
1041: x[ini+j] = 1.0;
1042: for (i=j+1;i<n1;i++) x[i+ini] = W[i+n1*j];
1043: tau[j] = PetscConj(tau[j]);
1044: PetscCallBLAS("LAPACKlarf",LAPACKlarf_("L",&dim,&one,x,&one,tau+j,t,&dim,work+nwu));
1045: }
1046: PetscCall(PetscBLASIntCast(lds,&lds_));
1047: PetscCallBLAS("BLASgemv",BLASgemv_("N",&lds_,&n1_,&sone,S,&lds_,t,&one,&szero,cont,&one));
1048: PetscCall(PetscFPTrapPop());
1049: }
1050: PetscFunctionReturn(PETSC_SUCCESS);
1051: }
1053: /*
1054: Compute a run of Arnoldi iterations
1055: */
1056: static PetscErrorCode NEPNLEIGSTOARrun(NEP nep,Mat MK,Mat MH,BV W,PetscInt k,PetscInt *M,PetscReal *betah,PetscScalar *betak,PetscBool *breakdown,Vec *t_)
1057: {
1058: NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
1059: PetscInt i,j,m=*M,lwa,deg=ctx->nmat-1,lds,nqt,ld,l,ldh;
1060: Vec t;
1061: PetscReal norm=0.0;
1062: PetscScalar *x,*work,*tt,sigma=1.0,*cont,*S,*K=NULL,*H;
1063: PetscBool lindep;
1064: Mat MS;
1066: PetscFunctionBegin;
1067: *betah = 0.0; *betak = 0.0;
1068: PetscCall(MatDenseGetArray(MH,&H));
1069: if (MK) PetscCall(MatDenseGetArray(MK,&K));
1070: PetscCall(MatDenseGetLDA(MH,&ldh));
1071: PetscCall(BVTensorGetFactors(ctx->V,NULL,&MS));
1072: PetscCall(MatDenseGetArray(MS,&S));
1073: PetscCall(BVGetSizes(nep->V,NULL,NULL,&ld));
1074: lds = ld*deg;
1075: PetscCall(BVGetActiveColumns(nep->V,&l,&nqt));
1076: lwa = PetscMax(ld,deg)+(m+1)*(m+1)+4*(m+1);
1077: PetscCall(PetscMalloc4(ld,&x,lwa,&work,m+1,&tt,lds,&cont));
1078: PetscCall(BVSetActiveColumns(ctx->V,0,m));
1079: for (j=k;j<m;j++) {
1080: sigma = ctx->shifts[(++ctx->idxrk)%ctx->nshiftsw];
1082: /* Continuation vector */
1083: PetscCall(NEPNLEIGS_RKcontinuation(nep,0,j,K,H,ldh,sigma,S,lds,cont,tt,work));
1085: /* apply operator */
1086: PetscCall(BVGetColumn(nep->V,nqt,&t));
1087: PetscCall(NEPTOARExtendBasis(nep,(ctx->idxrk)%ctx->nshiftsw,cont,ld,nqt,W,nep->V,t,S+(j+1)*lds,ld,t_));
1088: PetscCall(BVRestoreColumn(nep->V,nqt,&t));
1090: /* orthogonalize */
1091: PetscCall(BVOrthogonalizeColumn(nep->V,nqt,x,&norm,&lindep));
1092: if (!lindep) {
1093: x[nqt] = norm;
1094: PetscCall(BVScaleColumn(nep->V,nqt,1.0/norm));
1095: nqt++;
1096: } else x[nqt] = 0.0;
1098: PetscCall(NEPTOARCoefficients(nep,sigma,nqt-1,cont,ld,S+(j+1)*lds,ld,x,work));
1100: /* Level-2 orthogonalization */
1101: PetscCall(BVOrthogonalizeColumn(ctx->V,j+1,H+j*ldh,&norm,breakdown));
1102: H[j+1+ldh*j] = norm;
1103: if (ctx->nshifts && MK) {
1104: for (i=0;i<=j;i++) K[i+ldh*j] = sigma*H[i+ldh*j] + tt[i];
1105: K[j+1+ldh*j] = sigma*H[j+1+ldh*j];
1106: }
1107: if (*breakdown) {
1108: *M = j+1;
1109: break;
1110: }
1111: PetscCall(BVScaleColumn(ctx->V,j+1,1.0/norm));
1112: PetscCall(BVSetActiveColumns(nep->V,l,nqt));
1113: }
1114: *betah = norm;
1115: if (ctx->nshifts) *betak = norm*sigma;
1116: PetscCall(PetscFree4(x,work,tt,cont));
1117: PetscCall(MatDenseRestoreArray(MS,&S));
1118: PetscCall(MatDenseRestoreArray(MH,&H));
1119: if (MK) PetscCall(MatDenseRestoreArray(MK,&K));
1120: PetscCall(BVTensorRestoreFactors(ctx->V,NULL,&MS));
1121: PetscFunctionReturn(PETSC_SUCCESS);
1122: }
1124: PetscErrorCode NEPSolve_NLEIGS(NEP nep)
1125: {
1126: NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
1127: PetscInt i,k=0,l,nv=0,ld,lds,nq;
1128: PetscInt deg=ctx->nmat-1,nconv=0,dsn,dsk;
1129: PetscScalar *pU,betak=0,*eigr,*eigi;
1130: const PetscScalar *S;
1131: PetscReal betah;
1132: PetscBool falselock=PETSC_FALSE,breakdown=PETSC_FALSE;
1133: BV W;
1134: Mat H,K=NULL,MS,MQ,U;
1136: PetscFunctionBegin;
1137: if (ctx->lock) {
1138: /* undocumented option to use a cheaper locking instead of the true locking */
1139: PetscCall(PetscOptionsGetBool(NULL,NULL,"-nep_nleigs_falselocking",&falselock,NULL));
1140: }
1142: PetscCall(BVGetSizes(nep->V,NULL,NULL,&ld));
1143: lds = deg*ld;
1144: if (!ctx->nshifts) PetscCall(PetscMalloc2(nep->ncv,&eigr,nep->ncv,&eigi));
1145: else { eigr = nep->eigr; eigi = nep->eigi; }
1146: PetscCall(BVDuplicateResize(nep->V,PetscMax(nep->nt-1,ctx->nmat-1),&W));
1148: /* clean projected matrix (including the extra-arrow) */
1149: PetscCall(DSSetDimensions(nep->ds,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_DETERMINE));
1150: PetscCall(DSGetMat(nep->ds,DS_MAT_A,&H));
1151: PetscCall(MatZeroEntries(H));
1152: PetscCall(DSRestoreMat(nep->ds,DS_MAT_A,&H));
1153: if (ctx->nshifts) {
1154: PetscCall(DSGetMat(nep->ds,DS_MAT_B,&H));
1155: PetscCall(MatZeroEntries(H));
1156: PetscCall(DSRestoreMat(nep->ds,DS_MAT_B,&H));
1157: }
1159: /* Get the starting Arnoldi vector */
1160: PetscCall(BVTensorBuildFirstColumn(ctx->V,nep->nini));
1162: /* Restart loop */
1163: l = 0;
1164: while (nep->reason == NEP_CONVERGED_ITERATING) {
1165: nep->its++;
1167: /* Compute an nv-step Krylov relation */
1168: nv = PetscMin(nep->nconv+nep->mpd,nep->ncv);
1169: if (ctx->nshifts) PetscCall(DSGetMat(nep->ds,DS_MAT_A,&K));
1170: PetscCall(DSGetMat(nep->ds,ctx->nshifts?DS_MAT_B:DS_MAT_A,&H));
1171: PetscCall(NEPNLEIGSTOARrun(nep,K,H,W,nep->nconv+l,&nv,&betah,&betak,&breakdown,nep->work));
1172: PetscCall(DSRestoreMat(nep->ds,ctx->nshifts?DS_MAT_B:DS_MAT_A,&H));
1173: if (ctx->nshifts) PetscCall(DSRestoreMat(nep->ds,DS_MAT_A,&K));
1174: PetscCall(DSSetDimensions(nep->ds,nv,nep->nconv,nep->nconv+l));
1175: if (l==0) PetscCall(DSSetState(nep->ds,DS_STATE_INTERMEDIATE));
1176: else PetscCall(DSSetState(nep->ds,DS_STATE_RAW));
1178: /* Solve projected problem */
1179: PetscCall(DSSolve(nep->ds,nep->eigr,nep->eigi));
1180: PetscCall(DSSort(nep->ds,nep->eigr,nep->eigi,NULL,NULL,NULL));
1181: PetscCall(DSUpdateExtraRow(nep->ds));
1182: PetscCall(DSSynchronize(nep->ds,nep->eigr,nep->eigi));
1184: /* Check convergence */
1185: PetscCall(NEPNLEIGSKrylovConvergence(nep,PETSC_FALSE,nep->nconv,nv-nep->nconv,betah,betak,&k,nep->work));
1186: PetscCall((*nep->stopping)(nep,nep->its,nep->max_it,k,nep->nev,&nep->reason,nep->stoppingctx));
1188: /* Update l */
1189: if (nep->reason != NEP_CONVERGED_ITERATING || breakdown) l = 0;
1190: else {
1191: l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
1192: PetscCall(DSGetTruncateSize(nep->ds,k,nv,&l));
1193: if (!breakdown) {
1194: /* Prepare the Rayleigh quotient for restart */
1195: PetscCall(DSGetDimensions(nep->ds,&dsn,NULL,&dsk,NULL));
1196: PetscCall(DSSetDimensions(nep->ds,dsn,k,dsk));
1197: PetscCall(DSTruncate(nep->ds,k+l,PETSC_FALSE));
1198: }
1199: }
1200: nconv = k;
1201: if (!ctx->lock && nep->reason == NEP_CONVERGED_ITERATING && !breakdown) { l += k; k = 0; }
1202: if (l) PetscCall(PetscInfo(nep,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
1204: /* Update S */
1205: PetscCall(DSGetMat(nep->ds,ctx->nshifts?DS_MAT_Z:DS_MAT_Q,&MQ));
1206: PetscCall(BVMultInPlace(ctx->V,MQ,nep->nconv,k+l));
1207: PetscCall(DSRestoreMat(nep->ds,ctx->nshifts?DS_MAT_Z:DS_MAT_Q,&MQ));
1209: /* Copy last column of S */
1210: PetscCall(BVCopyColumn(ctx->V,nv,k+l));
1212: if (breakdown && nep->reason == NEP_CONVERGED_ITERATING) {
1213: /* Stop if breakdown */
1214: PetscCall(PetscInfo(nep,"Breakdown (it=%" PetscInt_FMT " norm=%g)\n",nep->its,(double)betah));
1215: nep->reason = NEP_DIVERGED_BREAKDOWN;
1216: }
1217: if (nep->reason != NEP_CONVERGED_ITERATING) l--;
1218: /* truncate S */
1219: PetscCall(BVGetActiveColumns(nep->V,NULL,&nq));
1220: if (k+l+deg<=nq) {
1221: PetscCall(BVSetActiveColumns(ctx->V,nep->nconv,k+l+1));
1222: if (!falselock && ctx->lock) PetscCall(BVTensorCompress(ctx->V,k-nep->nconv));
1223: else PetscCall(BVTensorCompress(ctx->V,0));
1224: }
1225: nep->nconv = k;
1226: if (!ctx->nshifts) {
1227: for (i=0;i<nv;i++) { eigr[i] = nep->eigr[i]; eigi[i] = nep->eigi[i]; }
1228: PetscCall(NEPNLEIGSBackTransform((PetscObject)nep,nv,eigr,eigi));
1229: }
1230: PetscCall(NEPMonitor(nep,nep->its,nconv,eigr,eigi,nep->errest,nv));
1231: }
1232: nep->nconv = nconv;
1233: if (nep->nconv>0) {
1234: PetscCall(BVSetActiveColumns(ctx->V,0,nep->nconv));
1235: PetscCall(BVGetActiveColumns(nep->V,NULL,&nq));
1236: PetscCall(BVSetActiveColumns(nep->V,0,nq));
1237: if (nq>nep->nconv) {
1238: PetscCall(BVTensorCompress(ctx->V,nep->nconv));
1239: PetscCall(BVSetActiveColumns(nep->V,0,nep->nconv));
1240: nq = nep->nconv;
1241: }
1242: if (ctx->nshifts) {
1243: PetscCall(DSGetMat(nep->ds,DS_MAT_B,&MQ));
1244: PetscCall(BVMultInPlace(ctx->V,MQ,0,nep->nconv));
1245: PetscCall(DSRestoreMat(nep->ds,DS_MAT_B,&MQ));
1246: }
1247: PetscCall(BVTensorGetFactors(ctx->V,NULL,&MS));
1248: PetscCall(MatDenseGetArrayRead(MS,&S));
1249: PetscCall(PetscMalloc1(nq*nep->nconv,&pU));
1250: for (i=0;i<nep->nconv;i++) PetscCall(PetscArraycpy(pU+i*nq,S+i*lds,nq));
1251: PetscCall(MatDenseRestoreArrayRead(MS,&S));
1252: PetscCall(BVTensorRestoreFactors(ctx->V,NULL,&MS));
1253: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,nq,nep->nconv,pU,&U));
1254: PetscCall(BVSetActiveColumns(nep->V,0,nq));
1255: PetscCall(BVMultInPlace(nep->V,U,0,nep->nconv));
1256: PetscCall(BVSetActiveColumns(nep->V,0,nep->nconv));
1257: PetscCall(MatDestroy(&U));
1258: PetscCall(PetscFree(pU));
1259: PetscCall(DSTruncate(nep->ds,nep->nconv,PETSC_TRUE));
1260: }
1262: /* Map eigenvalues back to the original problem */
1263: if (!ctx->nshifts) {
1264: PetscCall(NEPNLEIGSBackTransform((PetscObject)nep,nep->nconv,nep->eigr,nep->eigi));
1265: PetscCall(PetscFree2(eigr,eigi));
1266: }
1267: PetscCall(BVDestroy(&W));
1268: PetscFunctionReturn(PETSC_SUCCESS);
1269: }
1271: static PetscErrorCode NEPNLEIGSSetSingularitiesFunction_NLEIGS(NEP nep,NEPNLEIGSSingularitiesFn *fun,PetscCtx ctx)
1272: {
1273: NEP_NLEIGS *nepctx=(NEP_NLEIGS*)nep->data;
1275: PetscFunctionBegin;
1276: if (fun) nepctx->computesingularities = fun;
1277: if (ctx) nepctx->singularitiesctx = ctx;
1278: nep->state = NEP_STATE_INITIAL;
1279: PetscFunctionReturn(PETSC_SUCCESS);
1280: }
1282: /*@
1283: NEPNLEIGSSetSingularitiesFunction - Sets a user-defined callback function
1284: to compute a discretization of the singularity set (the values where
1285: $T(\cdot)$ is not analytic).
1287: Logically Collective
1289: Input Parameters:
1290: + nep - the nonlinear eigensolver context
1291: . fun - user function (if `NULL` then `NEP` retains any previously set value)
1292: - ctx - [optional] user-defined context for private data for the function
1293: (may be `NULL`, in which case `NEP` retains any previously set value)
1295: Notes:
1296: If the problem type has been set to `NEP_RATIONAL` with `NEPSetProblemType()`,
1297: then it is not necessary to set the singularities explicitly since the
1298: solver will try to determine them automatically.
1300: If the problem is `NEP_GENERAL`, it is also possible to omit the
1301: singularities callback. In that case, a discretization of the singularity
1302: set is approximated via the AAA algorithm {cite:p}`Nak18,Els19`.
1304: Level: intermediate
1306: .seealso: [](ch:nep), `NEPNLEIGS`, `NEPNLEIGSGetSingularitiesFunction()`, `NEPSetProblemType()`
1307: @*/
1308: PetscErrorCode NEPNLEIGSSetSingularitiesFunction(NEP nep,NEPNLEIGSSingularitiesFn *fun,PetscCtx ctx)
1309: {
1310: PetscFunctionBegin;
1312: PetscTryMethod(nep,"NEPNLEIGSSetSingularitiesFunction_C",(NEP,NEPNLEIGSSingularitiesFn*,PetscCtx),(nep,fun,ctx));
1313: PetscFunctionReturn(PETSC_SUCCESS);
1314: }
1316: static PetscErrorCode NEPNLEIGSGetSingularitiesFunction_NLEIGS(NEP nep,NEPNLEIGSSingularitiesFn **fun,PetscCtxRt ctx)
1317: {
1318: NEP_NLEIGS *nepctx=(NEP_NLEIGS*)nep->data;
1320: PetscFunctionBegin;
1321: if (fun) *fun = nepctx->computesingularities;
1322: if (ctx) *(void**)ctx = nepctx->singularitiesctx;
1323: PetscFunctionReturn(PETSC_SUCCESS);
1324: }
1326: /*@
1327: NEPNLEIGSGetSingularitiesFunction - Returns the callback function and optionally the user
1328: provided context for computing a discretization of the singularity set.
1330: Not Collective
1332: Input Parameter:
1333: . nep - the nonlinear eigensolver context
1335: Output Parameters:
1336: + fun - location to put the function (or `NULL`)
1337: - ctx - location to stash the function context (or `NULL`)
1339: Level: intermediate
1341: .seealso: [](ch:nep), `NEPNLEIGS`, `NEPNLEIGSSetSingularitiesFunction()`
1342: @*/
1343: PetscErrorCode NEPNLEIGSGetSingularitiesFunction(NEP nep,NEPNLEIGSSingularitiesFn **fun,PetscCtxRt ctx)
1344: {
1345: PetscFunctionBegin;
1347: PetscUseMethod(nep,"NEPNLEIGSGetSingularitiesFunction_C",(NEP,NEPNLEIGSSingularitiesFn**,PetscCtxRt),(nep,fun,ctx));
1348: PetscFunctionReturn(PETSC_SUCCESS);
1349: }
1351: static PetscErrorCode NEPNLEIGSSetRestart_NLEIGS(NEP nep,PetscReal keep)
1352: {
1353: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1355: PetscFunctionBegin;
1356: if (keep==(PetscReal)PETSC_DEFAULT || keep==(PetscReal)PETSC_DECIDE) ctx->keep = 0.5;
1357: else {
1358: 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]");
1359: ctx->keep = keep;
1360: }
1361: PetscFunctionReturn(PETSC_SUCCESS);
1362: }
1364: /*@
1365: NEPNLEIGSSetRestart - Sets the restart parameter for the NLEIGS
1366: method, in particular the proportion of basis vectors that must be kept
1367: after restart.
1369: Logically Collective
1371: Input Parameters:
1372: + nep - the nonlinear eigensolver context
1373: - keep - the number of vectors to be kept at restart
1375: Options Database Key:
1376: . -nep_nleigs_restart keep - sets the restart parameter
1378: Notes:
1379: Allowed values are in the range [0.1,0.9]. The default is 0.5.
1381: Level: advanced
1383: .seealso: [](ch:nep), `NEPNLEIGS`, `NEPNLEIGSGetRestart()`
1384: @*/
1385: PetscErrorCode NEPNLEIGSSetRestart(NEP nep,PetscReal keep)
1386: {
1387: PetscFunctionBegin;
1390: PetscTryMethod(nep,"NEPNLEIGSSetRestart_C",(NEP,PetscReal),(nep,keep));
1391: PetscFunctionReturn(PETSC_SUCCESS);
1392: }
1394: static PetscErrorCode NEPNLEIGSGetRestart_NLEIGS(NEP nep,PetscReal *keep)
1395: {
1396: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1398: PetscFunctionBegin;
1399: *keep = ctx->keep;
1400: PetscFunctionReturn(PETSC_SUCCESS);
1401: }
1403: /*@
1404: NEPNLEIGSGetRestart - Gets the restart parameter used in the NLEIGS method.
1406: Not Collective
1408: Input Parameter:
1409: . nep - the nonlinear eigensolver context
1411: Output Parameter:
1412: . keep - the restart parameter
1414: Level: advanced
1416: .seealso: [](ch:nep), `NEPNLEIGS`, `NEPNLEIGSSetRestart()`
1417: @*/
1418: PetscErrorCode NEPNLEIGSGetRestart(NEP nep,PetscReal *keep)
1419: {
1420: PetscFunctionBegin;
1422: PetscAssertPointer(keep,2);
1423: PetscUseMethod(nep,"NEPNLEIGSGetRestart_C",(NEP,PetscReal*),(nep,keep));
1424: PetscFunctionReturn(PETSC_SUCCESS);
1425: }
1427: static PetscErrorCode NEPNLEIGSSetLocking_NLEIGS(NEP nep,PetscBool lock)
1428: {
1429: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1431: PetscFunctionBegin;
1432: ctx->lock = lock;
1433: PetscFunctionReturn(PETSC_SUCCESS);
1434: }
1436: /*@
1437: NEPNLEIGSSetLocking - Choose between locking and non-locking variants of
1438: the NLEIGS method.
1440: Logically Collective
1442: Input Parameters:
1443: + nep - the nonlinear eigensolver context
1444: - lock - true if the locking variant must be selected
1446: Options Database Key:
1447: . -nep_nleigs_locking (true|false) - sets the locking flag
1449: Notes:
1450: The default is to lock converged eigenpairs when the method restarts.
1451: This behavior can be changed so that all directions are kept in the
1452: working subspace even if already converged to working accuracy (the
1453: non-locking variant).
1455: Level: advanced
1457: .seealso: [](ch:nep), `NEPNLEIGS`, `NEPNLEIGSGetLocking()`
1458: @*/
1459: PetscErrorCode NEPNLEIGSSetLocking(NEP nep,PetscBool lock)
1460: {
1461: PetscFunctionBegin;
1464: PetscTryMethod(nep,"NEPNLEIGSSetLocking_C",(NEP,PetscBool),(nep,lock));
1465: PetscFunctionReturn(PETSC_SUCCESS);
1466: }
1468: static PetscErrorCode NEPNLEIGSGetLocking_NLEIGS(NEP nep,PetscBool *lock)
1469: {
1470: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1472: PetscFunctionBegin;
1473: *lock = ctx->lock;
1474: PetscFunctionReturn(PETSC_SUCCESS);
1475: }
1477: /*@
1478: NEPNLEIGSGetLocking - Gets the locking flag used in the NLEIGS method.
1480: Not Collective
1482: Input Parameter:
1483: . nep - the nonlinear eigensolver context
1485: Output Parameter:
1486: . lock - the locking flag
1488: Level: advanced
1490: .seealso: [](ch:nep), `NEPNLEIGS`, `NEPNLEIGSSetLocking()`
1491: @*/
1492: PetscErrorCode NEPNLEIGSGetLocking(NEP nep,PetscBool *lock)
1493: {
1494: PetscFunctionBegin;
1496: PetscAssertPointer(lock,2);
1497: PetscUseMethod(nep,"NEPNLEIGSGetLocking_C",(NEP,PetscBool*),(nep,lock));
1498: PetscFunctionReturn(PETSC_SUCCESS);
1499: }
1501: static PetscErrorCode NEPNLEIGSSetInterpolation_NLEIGS(NEP nep,PetscReal tol,PetscInt degree)
1502: {
1503: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1505: PetscFunctionBegin;
1506: if (tol == (PetscReal)PETSC_DETERMINE) {
1507: ctx->ddtol = PETSC_DETERMINE;
1508: nep->state = NEP_STATE_INITIAL;
1509: } else if (tol != (PetscReal)PETSC_CURRENT) {
1510: PetscCheck(tol>0.0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of tol. Must be > 0");
1511: ctx->ddtol = tol;
1512: }
1513: if (degree == PETSC_DETERMINE) {
1514: ctx->ddmaxit = 0;
1515: if (nep->state) PetscCall(NEPReset(nep));
1516: nep->state = NEP_STATE_INITIAL;
1517: } else if (degree != PETSC_CURRENT) {
1518: PetscCheck(degree>0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of degree. Must be > 0");
1519: if (ctx->ddmaxit != degree) {
1520: ctx->ddmaxit = degree;
1521: if (nep->state) PetscCall(NEPReset(nep));
1522: nep->state = NEP_STATE_INITIAL;
1523: }
1524: }
1525: PetscFunctionReturn(PETSC_SUCCESS);
1526: }
1528: /*@
1529: NEPNLEIGSSetInterpolation - Sets the tolerance and maximum degree
1530: when building the interpolation via divided differences.
1532: Collective
1534: Input Parameters:
1535: + nep - the nonlinear eigensolver context
1536: . tol - tolerance to stop computing divided differences
1537: - degree - maximum degree of interpolation
1539: Options Database Keys:
1540: + -nep_nleigs_interpolation_tol tol - sets the tolerance to stop computing divided differences
1541: - -nep_nleigs_interpolation_degree degree - sets the maximum degree of interpolation
1543: Note:
1544: `PETSC_CURRENT` can be used to preserve the current value of any of the
1545: arguments, and `PETSC_DETERMINE` to set them to a default value.
1547: Level: advanced
1549: .seealso: [](ch:nep), `NEPNLEIGS`, `NEPNLEIGSGetInterpolation()`
1550: @*/
1551: PetscErrorCode NEPNLEIGSSetInterpolation(NEP nep,PetscReal tol,PetscInt degree)
1552: {
1553: PetscFunctionBegin;
1557: PetscTryMethod(nep,"NEPNLEIGSSetInterpolation_C",(NEP,PetscReal,PetscInt),(nep,tol,degree));
1558: PetscFunctionReturn(PETSC_SUCCESS);
1559: }
1561: static PetscErrorCode NEPNLEIGSGetInterpolation_NLEIGS(NEP nep,PetscReal *tol,PetscInt *degree)
1562: {
1563: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1565: PetscFunctionBegin;
1566: if (tol) *tol = ctx->ddtol;
1567: if (degree) *degree = ctx->ddmaxit;
1568: PetscFunctionReturn(PETSC_SUCCESS);
1569: }
1571: /*@
1572: NEPNLEIGSGetInterpolation - Gets the tolerance and maximum degree
1573: when building the interpolation via divided differences.
1575: Not Collective
1577: Input Parameter:
1578: . nep - the nonlinear eigensolver context
1580: Output Parameters:
1581: + tol - tolerance to stop computing divided differences
1582: - degree - maximum degree of interpolation
1584: Level: advanced
1586: .seealso: [](ch:nep), `NEPNLEIGS`, `NEPNLEIGSSetInterpolation()`
1587: @*/
1588: PetscErrorCode NEPNLEIGSGetInterpolation(NEP nep,PetscReal *tol,PetscInt *degree)
1589: {
1590: PetscFunctionBegin;
1592: PetscTryMethod(nep,"NEPNLEIGSGetInterpolation_C",(NEP,PetscReal*,PetscInt*),(nep,tol,degree));
1593: PetscFunctionReturn(PETSC_SUCCESS);
1594: }
1596: static PetscErrorCode NEPNLEIGSSetRKShifts_NLEIGS(NEP nep,PetscInt ns,PetscScalar *shifts)
1597: {
1598: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1599: PetscInt i;
1601: PetscFunctionBegin;
1602: PetscCheck(ns>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"Number of shifts must be non-negative");
1603: if (ctx->nshifts) PetscCall(PetscFree(ctx->shifts));
1604: for (i=0;i<ctx->nshiftsw;i++) PetscCall(KSPDestroy(&ctx->ksp[i]));
1605: PetscCall(PetscFree(ctx->ksp));
1606: ctx->ksp = NULL;
1607: if (ns) {
1608: PetscCall(PetscMalloc1(ns,&ctx->shifts));
1609: for (i=0;i<ns;i++) ctx->shifts[i] = shifts[i];
1610: }
1611: ctx->nshifts = ns;
1612: nep->state = NEP_STATE_INITIAL;
1613: PetscFunctionReturn(PETSC_SUCCESS);
1614: }
1616: /*@
1617: NEPNLEIGSSetRKShifts - Sets a list of shifts to be used in the Rational
1618: Krylov method.
1620: Collective
1622: Input Parameters:
1623: + nep - the nonlinear eigensolver context
1624: . ns - number of shifts
1625: - shifts - array of scalar values specifying the shifts
1627: Options Database Key:
1628: . -nep_nleigs_rk_shifts s0,s1,... - sets the list of shifts
1630: Notes:
1631: If only one shift is provided, the built subspace is equivalent to
1632: shift-and-invert Krylov-Schur (provided that the absolute convergence
1633: criterion is used). Otherwise, the rational Krylov variant is run.
1635: In the case of real scalars, complex shifts are not allowed. In the
1636: command line, a comma-separated list of complex values can be provided with
1637: the format `[+/-][realnumber][+/-]realnumberi` with no spaces, e.g.
1638: `-nep_nleigs_rk_shifts 1.0+2.0i,1.5+2.0i,1.0+1.5i`.
1640: Use `ns=0` to remove previously set shifts.
1642: Level: advanced
1644: .seealso: [](ch:nep), `NEPNLEIGS`, `NEPNLEIGSGetRKShifts()`
1645: @*/
1646: PetscErrorCode NEPNLEIGSSetRKShifts(NEP nep,PetscInt ns,PetscScalar shifts[])
1647: {
1648: PetscFunctionBegin;
1651: if (ns) PetscAssertPointer(shifts,3);
1652: PetscTryMethod(nep,"NEPNLEIGSSetRKShifts_C",(NEP,PetscInt,PetscScalar*),(nep,ns,shifts));
1653: PetscFunctionReturn(PETSC_SUCCESS);
1654: }
1656: static PetscErrorCode NEPNLEIGSGetRKShifts_NLEIGS(NEP nep,PetscInt *ns,PetscScalar **shifts)
1657: {
1658: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1659: PetscInt i;
1661: PetscFunctionBegin;
1662: *ns = ctx->nshifts;
1663: if (ctx->nshifts) {
1664: PetscCall(PetscMalloc1(ctx->nshifts,shifts));
1665: for (i=0;i<ctx->nshifts;i++) (*shifts)[i] = ctx->shifts[i];
1666: }
1667: PetscFunctionReturn(PETSC_SUCCESS);
1668: }
1670: /*@
1671: NEPNLEIGSGetRKShifts - Gets the list of shifts used in the Rational
1672: Krylov method.
1674: Not Collective
1676: Input Parameter:
1677: . nep - the nonlinear eigensolver context
1679: Output Parameters:
1680: + ns - number of shifts
1681: - shifts - array of shifts
1683: Note:
1684: The user is responsible for deallocating the returned array.
1686: Level: advanced
1688: .seealso: [](ch:nep), `NEPNLEIGS`, `NEPNLEIGSSetRKShifts()`
1689: @*/
1690: PetscErrorCode NEPNLEIGSGetRKShifts(NEP nep,PetscInt *ns,PetscScalar *shifts[]) PeNS
1691: {
1692: PetscFunctionBegin;
1694: PetscAssertPointer(ns,2);
1695: PetscAssertPointer(shifts,3);
1696: PetscTryMethod(nep,"NEPNLEIGSGetRKShifts_C",(NEP,PetscInt*,PetscScalar**),(nep,ns,shifts));
1697: PetscFunctionReturn(PETSC_SUCCESS);
1698: }
1700: static PetscErrorCode NEPNLEIGSGetKSPs_NLEIGS(NEP nep,PetscInt *nsolve,KSP **ksp)
1701: {
1702: NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
1703: PetscInt i;
1704: PC pc;
1706: PetscFunctionBegin;
1707: if (!ctx->ksp) {
1708: PetscCall(NEPNLEIGSSetShifts(nep,&ctx->nshiftsw));
1709: PetscCall(PetscMalloc1(ctx->nshiftsw,&ctx->ksp));
1710: for (i=0;i<ctx->nshiftsw;i++) {
1711: PetscCall(KSPCreate(PetscObjectComm((PetscObject)nep),&ctx->ksp[i]));
1712: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ctx->ksp[i],(PetscObject)nep,1));
1713: PetscCall(KSPSetOptionsPrefix(ctx->ksp[i],((PetscObject)nep)->prefix));
1714: PetscCall(KSPAppendOptionsPrefix(ctx->ksp[i],"nep_nleigs_"));
1715: PetscCall(PetscObjectSetOptions((PetscObject)ctx->ksp[i],((PetscObject)nep)->options));
1716: PetscCall(KSPSetErrorIfNotConverged(ctx->ksp[i],PETSC_TRUE));
1717: PetscCall(KSPSetTolerances(ctx->ksp[i],1e-3*SlepcDefaultTol(nep->tol),PETSC_CURRENT,PETSC_CURRENT,PETSC_CURRENT));
1718: PetscCall(KSPGetPC(ctx->ksp[i],&pc));
1719: if ((nep->fui==NEP_USER_INTERFACE_SPLIT && nep->P) || (nep->fui==NEP_USER_INTERFACE_CALLBACK && nep->function_pre!=nep->function)) {
1720: PetscCall(KSPSetType(ctx->ksp[i],KSPBCGS));
1721: PetscCall(PCSetType(pc,PCBJACOBI));
1722: } else {
1723: PetscCall(KSPSetType(ctx->ksp[i],KSPPREONLY));
1724: PetscCall(PCSetType(pc,PCLU));
1725: }
1726: }
1727: }
1728: if (nsolve) *nsolve = ctx->nshiftsw;
1729: if (ksp) *ksp = ctx->ksp;
1730: PetscFunctionReturn(PETSC_SUCCESS);
1731: }
1733: /*@
1734: NEPNLEIGSGetKSPs - Retrieve the array of linear solver objects associated with
1735: the nonlinear eigenvalue solver.
1737: Collective
1739: Input Parameter:
1740: . nep - the nonlinear eigensolver context
1742: Output Parameters:
1743: + nsolve - number of returned `KSP` objects
1744: - ksp - array of linear solver object
1746: Note:
1747: The number of `KSP` objects is equal to the number of shifts provided by the user,
1748: or 1 if the user did not provide shifts.
1750: Level: advanced
1752: .seealso: [](ch:nep), `NEPNLEIGS`, `NEPNLEIGSSetRKShifts()`
1753: @*/
1754: PetscErrorCode NEPNLEIGSGetKSPs(NEP nep,PetscInt *nsolve,KSP **ksp)
1755: {
1756: PetscFunctionBegin;
1758: PetscUseMethod(nep,"NEPNLEIGSGetKSPs_C",(NEP,PetscInt*,KSP**),(nep,nsolve,ksp));
1759: PetscFunctionReturn(PETSC_SUCCESS);
1760: }
1762: static PetscErrorCode NEPNLEIGSSetFullBasis_NLEIGS(NEP nep,PetscBool fullbasis)
1763: {
1764: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1766: PetscFunctionBegin;
1767: if (fullbasis!=ctx->fullbasis) {
1768: ctx->fullbasis = fullbasis;
1769: nep->state = NEP_STATE_INITIAL;
1770: nep->useds = PetscNot(fullbasis);
1771: }
1772: PetscFunctionReturn(PETSC_SUCCESS);
1773: }
1775: /*@
1776: NEPNLEIGSSetFullBasis - Choose between TOAR-basis (default) and full-basis
1777: variants of the NLEIGS method.
1779: Logically Collective
1781: Input Parameters:
1782: + nep - the nonlinear eigensolver context
1783: - fullbasis - true if the full-basis variant must be selected
1785: Options Database Key:
1786: . -nep_nleigs_full_basis (true|false) - sets the full-basis flag
1788: Notes:
1789: The default is to use a compact representation of the Krylov basis, that is,
1790: $V = (I \otimes U) S$, with a `BVTENSOR`. This behavior can be changed so that
1791: the full basis $V$ is explicitly stored and operated with. This variant is more
1792: expensive in terms of memory and computation, but is necessary in some cases,
1793: particularly for two-sided computations, see `NEPSetTwoSided()`.
1795: In the full-basis variant, the NLEIGS solver uses an `EPS` object to explicitly
1796: solve the linearized eigenproblem, see `NEPNLEIGSGetEPS()`.
1798: Level: advanced
1800: .seealso: [](ch:nep), `NEPNLEIGS`, `NEPNLEIGSGetFullBasis()`, `NEPNLEIGSGetEPS()`, `NEPSetTwoSided()`, `BVCreateTensor()`
1801: @*/
1802: PetscErrorCode NEPNLEIGSSetFullBasis(NEP nep,PetscBool fullbasis)
1803: {
1804: PetscFunctionBegin;
1807: PetscTryMethod(nep,"NEPNLEIGSSetFullBasis_C",(NEP,PetscBool),(nep,fullbasis));
1808: PetscFunctionReturn(PETSC_SUCCESS);
1809: }
1811: static PetscErrorCode NEPNLEIGSGetFullBasis_NLEIGS(NEP nep,PetscBool *fullbasis)
1812: {
1813: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1815: PetscFunctionBegin;
1816: *fullbasis = ctx->fullbasis;
1817: PetscFunctionReturn(PETSC_SUCCESS);
1818: }
1820: /*@
1821: NEPNLEIGSGetFullBasis - Gets the flag that indicates if NLEIGS is using the
1822: full-basis variant.
1824: Not Collective
1826: Input Parameter:
1827: . nep - the nonlinear eigensolver context
1829: Output Parameter:
1830: . fullbasis - the flag
1832: Level: advanced
1834: .seealso: [](ch:nep), `NEPNLEIGS`, `NEPNLEIGSSetFullBasis()`
1835: @*/
1836: PetscErrorCode NEPNLEIGSGetFullBasis(NEP nep,PetscBool *fullbasis)
1837: {
1838: PetscFunctionBegin;
1840: PetscAssertPointer(fullbasis,2);
1841: PetscUseMethod(nep,"NEPNLEIGSGetFullBasis_C",(NEP,PetscBool*),(nep,fullbasis));
1842: PetscFunctionReturn(PETSC_SUCCESS);
1843: }
1845: #define SHIFTMAX 30
1847: static PetscErrorCode NEPSetFromOptions_NLEIGS(NEP nep,PetscOptionItems PetscOptionsObject)
1848: {
1849: NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
1850: PetscInt i=0,k;
1851: PetscBool flg1,flg2,b;
1852: PetscReal r;
1853: PetscScalar array[SHIFTMAX];
1855: PetscFunctionBegin;
1856: PetscOptionsHeadBegin(PetscOptionsObject,"NEP NLEIGS Options");
1858: PetscCall(PetscOptionsReal("-nep_nleigs_restart","Proportion of vectors kept after restart","NEPNLEIGSSetRestart",0.5,&r,&flg1));
1859: if (flg1) PetscCall(NEPNLEIGSSetRestart(nep,r));
1861: PetscCall(PetscOptionsBool("-nep_nleigs_locking","Choose between locking and non-locking variants","NEPNLEIGSSetLocking",PETSC_FALSE,&b,&flg1));
1862: if (flg1) PetscCall(NEPNLEIGSSetLocking(nep,b));
1864: PetscCall(PetscOptionsBool("-nep_nleigs_full_basis","Choose between TOAR and full-basis variants","NEPNLEIGSSetFullBasis",PETSC_FALSE,&b,&flg1));
1865: if (flg1) PetscCall(NEPNLEIGSSetFullBasis(nep,b));
1867: PetscCall(NEPNLEIGSGetInterpolation(nep,&r,&i));
1868: if (!i) i = PETSC_DETERMINE;
1869: PetscCall(PetscOptionsInt("-nep_nleigs_interpolation_degree","Maximum number of terms for interpolation via divided differences","NEPNLEIGSSetInterpolation",i,&i,&flg1));
1870: PetscCall(PetscOptionsReal("-nep_nleigs_interpolation_tol","Tolerance for interpolation via divided differences","NEPNLEIGSSetInterpolation",r,&r,&flg2));
1871: if (flg1 || flg2) PetscCall(NEPNLEIGSSetInterpolation(nep,r,i));
1873: k = SHIFTMAX;
1874: for (i=0;i<k;i++) array[i] = 0;
1875: PetscCall(PetscOptionsScalarArray("-nep_nleigs_rk_shifts","Shifts for Rational Krylov","NEPNLEIGSSetRKShifts",array,&k,&flg1));
1876: if (flg1) PetscCall(NEPNLEIGSSetRKShifts(nep,k,array));
1878: PetscOptionsHeadEnd();
1880: if (!ctx->ksp) PetscCall(NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp));
1881: for (i=0;i<ctx->nshiftsw;i++) PetscCall(KSPSetFromOptions(ctx->ksp[i]));
1883: if (ctx->fullbasis) {
1884: if (!ctx->eps) PetscCall(NEPNLEIGSGetEPS(nep,&ctx->eps));
1885: PetscCall(EPSSetFromOptions(ctx->eps));
1886: }
1887: PetscFunctionReturn(PETSC_SUCCESS);
1888: }
1890: static PetscErrorCode NEPView_NLEIGS(NEP nep,PetscViewer viewer)
1891: {
1892: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1893: PetscBool isascii;
1894: PetscInt i;
1895: char str[50];
1897: PetscFunctionBegin;
1898: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
1899: if (isascii) {
1900: PetscCall(PetscViewerASCIIPrintf(viewer," %d%% of basis vectors kept after restart\n",(int)(100*ctx->keep)));
1901: if (ctx->fullbasis) PetscCall(PetscViewerASCIIPrintf(viewer," using the full-basis variant\n"));
1902: else PetscCall(PetscViewerASCIIPrintf(viewer," using the %slocking variant\n",ctx->lock?"":"non-"));
1903: PetscCall(PetscViewerASCIIPrintf(viewer," divided difference terms: used=%" PetscInt_FMT ", max=%" PetscInt_FMT "\n",ctx->nmat,ctx->ddmaxit));
1904: PetscCall(PetscViewerASCIIPrintf(viewer," tolerance for divided difference convergence: %g\n",(double)ctx->ddtol));
1905: if (ctx->nshifts) {
1906: PetscCall(PetscViewerASCIIPrintf(viewer," RK shifts: "));
1907: PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE));
1908: for (i=0;i<ctx->nshifts;i++) {
1909: PetscCall(SlepcSNPrintfScalar(str,sizeof(str),ctx->shifts[i],PETSC_FALSE));
1910: PetscCall(PetscViewerASCIIPrintf(viewer,"%s%s",str,(i<ctx->nshifts-1)?",":""));
1911: }
1912: PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
1913: PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE));
1914: }
1915: if (!ctx->ksp) PetscCall(NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp));
1916: PetscCall(PetscViewerASCIIPushTab(viewer));
1917: PetscCall(KSPView(ctx->ksp[0],viewer));
1918: PetscCall(PetscViewerASCIIPopTab(viewer));
1919: if (ctx->fullbasis) {
1920: if (!ctx->eps) PetscCall(NEPNLEIGSGetEPS(nep,&ctx->eps));
1921: PetscCall(PetscViewerASCIIPushTab(viewer));
1922: PetscCall(EPSView(ctx->eps,viewer));
1923: PetscCall(PetscViewerASCIIPopTab(viewer));
1924: }
1925: }
1926: PetscFunctionReturn(PETSC_SUCCESS);
1927: }
1929: static PetscErrorCode NEPReset_NLEIGS(NEP nep)
1930: {
1931: PetscInt k;
1932: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1934: PetscFunctionBegin;
1935: if (nep->fui==NEP_USER_INTERFACE_SPLIT) PetscCall(PetscFree(ctx->coeffD));
1936: else {
1937: for (k=0;k<ctx->nmat;k++) PetscCall(MatDestroy(&ctx->D[k]));
1938: }
1939: PetscCall(PetscFree4(ctx->s,ctx->xi,ctx->beta,ctx->D));
1940: for (k=0;k<ctx->nshiftsw;k++) PetscCall(KSPReset(ctx->ksp[k]));
1941: PetscCall(VecDestroy(&ctx->vrn));
1942: if (ctx->fullbasis) {
1943: PetscCall(MatDestroy(&ctx->A));
1944: PetscCall(EPSReset(ctx->eps));
1945: for (k=0;k<4;k++) PetscCall(VecDestroy(&ctx->w[k]));
1946: }
1947: PetscFunctionReturn(PETSC_SUCCESS);
1948: }
1950: static PetscErrorCode NEPDestroy_NLEIGS(NEP nep)
1951: {
1952: PetscInt k;
1953: NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
1955: PetscFunctionBegin;
1956: PetscCall(BVDestroy(&ctx->V));
1957: for (k=0;k<ctx->nshiftsw;k++) PetscCall(KSPDestroy(&ctx->ksp[k]));
1958: PetscCall(PetscFree(ctx->ksp));
1959: if (ctx->nshifts) PetscCall(PetscFree(ctx->shifts));
1960: if (ctx->fullbasis) PetscCall(EPSDestroy(&ctx->eps));
1961: PetscCall(PetscFree(nep->data));
1962: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetSingularitiesFunction_C",NULL));
1963: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetSingularitiesFunction_C",NULL));
1964: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRestart_C",NULL));
1965: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRestart_C",NULL));
1966: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetLocking_C",NULL));
1967: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetLocking_C",NULL));
1968: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetInterpolation_C",NULL));
1969: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetInterpolation_C",NULL));
1970: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRKShifts_C",NULL));
1971: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRKShifts_C",NULL));
1972: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetKSPs_C",NULL));
1973: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetFullBasis_C",NULL));
1974: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetFullBasis_C",NULL));
1975: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetEPS_C",NULL));
1976: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetEPS_C",NULL));
1977: PetscFunctionReturn(PETSC_SUCCESS);
1978: }
1980: /*MC
1981: NEPNLEIGS - NEPNLEIGS = "nleigs" - The NLEIGS method.
1983: Notes:
1984: This solver implements the NLEIGS method {cite:p}`Gut14`, which
1985: is based on rational interpolation followed by linearization.
1986: In our implementation, the linear eigensolver for the linearization
1987: operates with a compressed Krylov basis, as in the TOAR polynomial
1988: eigensolver, see the detailed description in {cite:p}`Cam21`.
1990: This method is particularly appropriate for nonlinear problems with
1991: singularities. The solver will try to determine the singularities
1992: automatically, but the user can also provide them with
1993: `NEPNLEIGSSetSingularitiesFunction()`.
1995: By default, the solver performs the static NLEIGS variant, with
1996: constant shift given by `NEPSetTarget()`. But the dynamic variant
1997: (rational Krylov) is also available if a list of shifts is given
1998: in `NEPNLEIGSSetRKShifts()`.
2000: `NEPNLEIGS` also implements a two-sided variant for computing left
2001: eigenvectors when `NEPSetTwoSided()` has been set.
2003: Apart from working with the compressed basis, it is also possible
2004: to enable the operation with an explicit basis for the linear
2005: eigensolver, see `NEPNLEIGSSetFullBasis()`. This allows using
2006: other eigensolvers via an `EPS` object obtained with `NEPNLEIGSGetEPS()`.
2007: Also, the explicit basis is activated in the two-sided variant.
2009: Level: beginner
2011: .seealso: [](ch:nep), `NEP`, `NEPType`, `NEPSetType()`, `NEPNLEIGSSetSingularitiesFunction()`, `NEPSetTarget()`, `NEPNLEIGSSetRKShifts()`, `NEPNLEIGSSetFullBasis()`, `NEPNLEIGSGetEPS()`, `NEPSetTwoSided()`
2012: M*/
2013: SLEPC_EXTERN PetscErrorCode NEPCreate_NLEIGS(NEP nep)
2014: {
2015: NEP_NLEIGS *ctx;
2017: PetscFunctionBegin;
2018: PetscCall(PetscNew(&ctx));
2019: nep->data = (void*)ctx;
2020: ctx->lock = PETSC_TRUE;
2021: ctx->ddtol = PETSC_DETERMINE;
2023: nep->useds = PETSC_TRUE;
2025: nep->ops->setup = NEPSetUp_NLEIGS;
2026: nep->ops->setfromoptions = NEPSetFromOptions_NLEIGS;
2027: nep->ops->view = NEPView_NLEIGS;
2028: nep->ops->destroy = NEPDestroy_NLEIGS;
2029: nep->ops->reset = NEPReset_NLEIGS;
2031: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetSingularitiesFunction_C",NEPNLEIGSSetSingularitiesFunction_NLEIGS));
2032: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetSingularitiesFunction_C",NEPNLEIGSGetSingularitiesFunction_NLEIGS));
2033: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRestart_C",NEPNLEIGSSetRestart_NLEIGS));
2034: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRestart_C",NEPNLEIGSGetRestart_NLEIGS));
2035: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetLocking_C",NEPNLEIGSSetLocking_NLEIGS));
2036: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetLocking_C",NEPNLEIGSGetLocking_NLEIGS));
2037: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetInterpolation_C",NEPNLEIGSSetInterpolation_NLEIGS));
2038: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetInterpolation_C",NEPNLEIGSGetInterpolation_NLEIGS));
2039: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRKShifts_C",NEPNLEIGSSetRKShifts_NLEIGS));
2040: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRKShifts_C",NEPNLEIGSGetRKShifts_NLEIGS));
2041: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetKSPs_C",NEPNLEIGSGetKSPs_NLEIGS));
2042: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetFullBasis_C",NEPNLEIGSSetFullBasis_NLEIGS));
2043: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetFullBasis_C",NEPNLEIGSGetFullBasis_NLEIGS));
2044: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetEPS_C",NEPNLEIGSSetEPS_NLEIGS));
2045: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetEPS_C",NEPNLEIGSGetEPS_NLEIGS));
2046: PetscFunctionReturn(PETSC_SUCCESS);
2047: }