Actual source code: nleigs.c
slepc-3.21.1 2024-04-26
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
10: /*
11: 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: void (*fun)(void);
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,(void(*)(void))MatMult_Fun));
587: PetscCall(MatShellSetOperation(*Ms,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_Fun));
588: PetscCall(MatShellSetOperation(*Ms,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Fun));
589: PetscCall(MatShellSetOperation(*Ms,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_Fun));
590: PetscCall(MatShellSetOperation(*Ms,MATOP_DESTROY,(void(*)(void))MatDestroy_Fun));
591: PetscCall(MatShellSetOperation(*Ms,MATOP_AXPY,(void(*)(void))MatAXPY_Fun));
592: PetscCall(MatShellSetOperation(*Ms,MATOP_SCALE,(void(*)(void))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_DEFAULT) 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_DEFAULT) nep->tol = SLEPC_DEFAULT_TOL;
890: if (ctx->ddtol==(PetscReal)PETSC_DEFAULT) 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_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT));
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,PetscErrorCode (*fun)(NEP,PetscInt*,PetscScalar*,void*),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 function to compute a discretization
1303: of the singularity set (where T(.) is not analytic).
1305: Logically Collective
1307: Input Parameters:
1308: + nep - the NEP context
1309: . fun - user function (if NULL then NEP retains any previously set value)
1310: - ctx - [optional] user-defined context for private data for the function
1311: (may be NULL, in which case NEP retains any previously set value)
1313: Calling sequence of fun:
1314: $ PetscErrorCode fun(NEP nep,PetscInt *maxnp,PetscScalar *xi,void *ctx)
1315: + nep - the NEP context
1316: . maxnp - on input number of requested points in the discretization (can be set)
1317: . xi - computed values of the discretization
1318: - ctx - optional context, as set by NEPNLEIGSSetSingularitiesFunction()
1320: Notes:
1321: The user-defined function can set a smaller value of maxnp if necessary.
1322: It is wrong to return a larger value.
1324: If the problem type has been set to rational with NEPSetProblemType(),
1325: then it is not necessary to set the singularities explicitly since the
1326: solver will try to determine them automatically.
1328: Level: intermediate
1330: .seealso: NEPNLEIGSGetSingularitiesFunction(), NEPSetProblemType()
1331: @*/
1332: PetscErrorCode NEPNLEIGSSetSingularitiesFunction(NEP nep,PetscErrorCode (*fun)(NEP nep,PetscInt *maxnp,PetscScalar *xi,void *ctx),void *ctx)
1333: {
1334: PetscFunctionBegin;
1336: PetscTryMethod(nep,"NEPNLEIGSSetSingularitiesFunction_C",(NEP,PetscErrorCode(*)(NEP,PetscInt*,PetscScalar*,void*),void*),(nep,fun,ctx));
1337: PetscFunctionReturn(PETSC_SUCCESS);
1338: }
1340: static PetscErrorCode NEPNLEIGSGetSingularitiesFunction_NLEIGS(NEP nep,PetscErrorCode (**fun)(NEP,PetscInt*,PetscScalar*,void*),void **ctx)
1341: {
1342: NEP_NLEIGS *nepctx=(NEP_NLEIGS*)nep->data;
1344: PetscFunctionBegin;
1345: if (fun) *fun = nepctx->computesingularities;
1346: if (ctx) *ctx = nepctx->singularitiesctx;
1347: PetscFunctionReturn(PETSC_SUCCESS);
1348: }
1350: /*@C
1351: NEPNLEIGSGetSingularitiesFunction - Returns the Function and optionally the user
1352: provided context for computing a discretization of the singularity set.
1354: Not Collective
1356: Input Parameter:
1357: . nep - the nonlinear eigensolver context
1359: Output Parameters:
1360: + fun - location to put the function (or NULL)
1361: - ctx - location to stash the function context (or NULL)
1363: Level: advanced
1365: .seealso: NEPNLEIGSSetSingularitiesFunction()
1366: @*/
1367: PetscErrorCode NEPNLEIGSGetSingularitiesFunction(NEP nep,PetscErrorCode (**fun)(NEP,PetscInt*,PetscScalar*,void*),void **ctx)
1368: {
1369: PetscFunctionBegin;
1371: PetscUseMethod(nep,"NEPNLEIGSGetSingularitiesFunction_C",(NEP,PetscErrorCode(**)(NEP,PetscInt*,PetscScalar*,void*),void**),(nep,fun,ctx));
1372: PetscFunctionReturn(PETSC_SUCCESS);
1373: }
1375: static PetscErrorCode NEPNLEIGSSetRestart_NLEIGS(NEP nep,PetscReal keep)
1376: {
1377: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1379: PetscFunctionBegin;
1380: if (keep==(PetscReal)PETSC_DEFAULT) ctx->keep = 0.5;
1381: else {
1382: PetscCheck(keep>=0.1 && keep<=0.9,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The keep argument must be in the range [0.1,0.9]");
1383: ctx->keep = keep;
1384: }
1385: PetscFunctionReturn(PETSC_SUCCESS);
1386: }
1388: /*@
1389: NEPNLEIGSSetRestart - Sets the restart parameter for the NLEIGS
1390: method, in particular the proportion of basis vectors that must be kept
1391: after restart.
1393: Logically Collective
1395: Input Parameters:
1396: + nep - the nonlinear eigensolver context
1397: - keep - the number of vectors to be kept at restart
1399: Options Database Key:
1400: . -nep_nleigs_restart - Sets the restart parameter
1402: Notes:
1403: Allowed values are in the range [0.1,0.9]. The default is 0.5.
1405: Level: advanced
1407: .seealso: NEPNLEIGSGetRestart()
1408: @*/
1409: PetscErrorCode NEPNLEIGSSetRestart(NEP nep,PetscReal keep)
1410: {
1411: PetscFunctionBegin;
1414: PetscTryMethod(nep,"NEPNLEIGSSetRestart_C",(NEP,PetscReal),(nep,keep));
1415: PetscFunctionReturn(PETSC_SUCCESS);
1416: }
1418: static PetscErrorCode NEPNLEIGSGetRestart_NLEIGS(NEP nep,PetscReal *keep)
1419: {
1420: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1422: PetscFunctionBegin;
1423: *keep = ctx->keep;
1424: PetscFunctionReturn(PETSC_SUCCESS);
1425: }
1427: /*@
1428: NEPNLEIGSGetRestart - Gets the restart parameter used in the NLEIGS method.
1430: Not Collective
1432: Input Parameter:
1433: . nep - the nonlinear eigensolver context
1435: Output Parameter:
1436: . keep - the restart parameter
1438: Level: advanced
1440: .seealso: NEPNLEIGSSetRestart()
1441: @*/
1442: PetscErrorCode NEPNLEIGSGetRestart(NEP nep,PetscReal *keep)
1443: {
1444: PetscFunctionBegin;
1446: PetscAssertPointer(keep,2);
1447: PetscUseMethod(nep,"NEPNLEIGSGetRestart_C",(NEP,PetscReal*),(nep,keep));
1448: PetscFunctionReturn(PETSC_SUCCESS);
1449: }
1451: static PetscErrorCode NEPNLEIGSSetLocking_NLEIGS(NEP nep,PetscBool lock)
1452: {
1453: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1455: PetscFunctionBegin;
1456: ctx->lock = lock;
1457: PetscFunctionReturn(PETSC_SUCCESS);
1458: }
1460: /*@
1461: NEPNLEIGSSetLocking - Choose between locking and non-locking variants of
1462: the NLEIGS method.
1464: Logically Collective
1466: Input Parameters:
1467: + nep - the nonlinear eigensolver context
1468: - lock - true if the locking variant must be selected
1470: Options Database Key:
1471: . -nep_nleigs_locking - Sets the locking flag
1473: Notes:
1474: The default is to lock converged eigenpairs when the method restarts.
1475: This behaviour can be changed so that all directions are kept in the
1476: working subspace even if already converged to working accuracy (the
1477: non-locking variant).
1479: Level: advanced
1481: .seealso: NEPNLEIGSGetLocking()
1482: @*/
1483: PetscErrorCode NEPNLEIGSSetLocking(NEP nep,PetscBool lock)
1484: {
1485: PetscFunctionBegin;
1488: PetscTryMethod(nep,"NEPNLEIGSSetLocking_C",(NEP,PetscBool),(nep,lock));
1489: PetscFunctionReturn(PETSC_SUCCESS);
1490: }
1492: static PetscErrorCode NEPNLEIGSGetLocking_NLEIGS(NEP nep,PetscBool *lock)
1493: {
1494: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1496: PetscFunctionBegin;
1497: *lock = ctx->lock;
1498: PetscFunctionReturn(PETSC_SUCCESS);
1499: }
1501: /*@
1502: NEPNLEIGSGetLocking - Gets the locking flag used in the NLEIGS method.
1504: Not Collective
1506: Input Parameter:
1507: . nep - the nonlinear eigensolver context
1509: Output Parameter:
1510: . lock - the locking flag
1512: Level: advanced
1514: .seealso: NEPNLEIGSSetLocking()
1515: @*/
1516: PetscErrorCode NEPNLEIGSGetLocking(NEP nep,PetscBool *lock)
1517: {
1518: PetscFunctionBegin;
1520: PetscAssertPointer(lock,2);
1521: PetscUseMethod(nep,"NEPNLEIGSGetLocking_C",(NEP,PetscBool*),(nep,lock));
1522: PetscFunctionReturn(PETSC_SUCCESS);
1523: }
1525: static PetscErrorCode NEPNLEIGSSetInterpolation_NLEIGS(NEP nep,PetscReal tol,PetscInt degree)
1526: {
1527: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1529: PetscFunctionBegin;
1530: if (tol == (PetscReal)PETSC_DEFAULT) {
1531: ctx->ddtol = PETSC_DEFAULT;
1532: nep->state = NEP_STATE_INITIAL;
1533: } else {
1534: PetscCheck(tol>0.0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of tol. Must be > 0");
1535: ctx->ddtol = tol;
1536: }
1537: if (degree == PETSC_DEFAULT || degree == PETSC_DECIDE) {
1538: ctx->ddmaxit = 0;
1539: if (nep->state) PetscCall(NEPReset(nep));
1540: nep->state = NEP_STATE_INITIAL;
1541: } else {
1542: PetscCheck(degree>0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of degree. Must be > 0");
1543: if (ctx->ddmaxit != degree) {
1544: ctx->ddmaxit = degree;
1545: if (nep->state) PetscCall(NEPReset(nep));
1546: nep->state = NEP_STATE_INITIAL;
1547: }
1548: }
1549: PetscFunctionReturn(PETSC_SUCCESS);
1550: }
1552: /*@
1553: NEPNLEIGSSetInterpolation - Sets the tolerance and maximum degree
1554: when building the interpolation via divided differences.
1556: Collective
1558: Input Parameters:
1559: + nep - the nonlinear eigensolver context
1560: . tol - tolerance to stop computing divided differences
1561: - degree - maximum degree of interpolation
1563: Options Database Key:
1564: + -nep_nleigs_interpolation_tol <tol> - Sets the tolerance to stop computing divided differences
1565: - -nep_nleigs_interpolation_degree <degree> - Sets the maximum degree of interpolation
1567: Notes:
1568: Use PETSC_DEFAULT for either argument to assign a reasonably good value.
1570: Level: advanced
1572: .seealso: NEPNLEIGSGetInterpolation()
1573: @*/
1574: PetscErrorCode NEPNLEIGSSetInterpolation(NEP nep,PetscReal tol,PetscInt degree)
1575: {
1576: PetscFunctionBegin;
1580: PetscTryMethod(nep,"NEPNLEIGSSetInterpolation_C",(NEP,PetscReal,PetscInt),(nep,tol,degree));
1581: PetscFunctionReturn(PETSC_SUCCESS);
1582: }
1584: static PetscErrorCode NEPNLEIGSGetInterpolation_NLEIGS(NEP nep,PetscReal *tol,PetscInt *degree)
1585: {
1586: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1588: PetscFunctionBegin;
1589: if (tol) *tol = ctx->ddtol;
1590: if (degree) *degree = ctx->ddmaxit;
1591: PetscFunctionReturn(PETSC_SUCCESS);
1592: }
1594: /*@
1595: NEPNLEIGSGetInterpolation - Gets the tolerance and maximum degree
1596: when building the interpolation via divided differences.
1598: Not Collective
1600: Input Parameter:
1601: . nep - the nonlinear eigensolver context
1603: Output Parameters:
1604: + tol - tolerance to stop computing divided differences
1605: - degree - maximum degree of interpolation
1607: Level: advanced
1609: .seealso: NEPNLEIGSSetInterpolation()
1610: @*/
1611: PetscErrorCode NEPNLEIGSGetInterpolation(NEP nep,PetscReal *tol,PetscInt *degree)
1612: {
1613: PetscFunctionBegin;
1615: PetscTryMethod(nep,"NEPNLEIGSGetInterpolation_C",(NEP,PetscReal*,PetscInt*),(nep,tol,degree));
1616: PetscFunctionReturn(PETSC_SUCCESS);
1617: }
1619: static PetscErrorCode NEPNLEIGSSetRKShifts_NLEIGS(NEP nep,PetscInt ns,PetscScalar *shifts)
1620: {
1621: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1622: PetscInt i;
1624: PetscFunctionBegin;
1625: PetscCheck(ns>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"Number of shifts must be non-negative");
1626: if (ctx->nshifts) PetscCall(PetscFree(ctx->shifts));
1627: for (i=0;i<ctx->nshiftsw;i++) PetscCall(KSPDestroy(&ctx->ksp[i]));
1628: PetscCall(PetscFree(ctx->ksp));
1629: ctx->ksp = NULL;
1630: if (ns) {
1631: PetscCall(PetscMalloc1(ns,&ctx->shifts));
1632: for (i=0;i<ns;i++) ctx->shifts[i] = shifts[i];
1633: }
1634: ctx->nshifts = ns;
1635: nep->state = NEP_STATE_INITIAL;
1636: PetscFunctionReturn(PETSC_SUCCESS);
1637: }
1639: /*@
1640: NEPNLEIGSSetRKShifts - Sets a list of shifts to be used in the Rational
1641: Krylov method.
1643: Collective
1645: Input Parameters:
1646: + nep - the nonlinear eigensolver context
1647: . ns - number of shifts
1648: - shifts - array of scalar values specifying the shifts
1650: Options Database Key:
1651: . -nep_nleigs_rk_shifts - Sets the list of shifts
1653: Notes:
1654: If only one shift is provided, the built subspace built is equivalent to
1655: shift-and-invert Krylov-Schur (provided that the absolute convergence
1656: criterion is used).
1658: In the case of real scalars, complex shifts are not allowed. In the
1659: command line, a comma-separated list of complex values can be provided with
1660: the format [+/-][realnumber][+/-]realnumberi with no spaces, e.g.
1661: -nep_nleigs_rk_shifts 1.0+2.0i,1.5+2.0i,1.0+1.5i
1663: Use ns=0 to remove previously set shifts.
1665: Level: advanced
1667: .seealso: NEPNLEIGSGetRKShifts()
1668: @*/
1669: PetscErrorCode NEPNLEIGSSetRKShifts(NEP nep,PetscInt ns,PetscScalar shifts[])
1670: {
1671: PetscFunctionBegin;
1674: if (ns) PetscAssertPointer(shifts,3);
1675: PetscTryMethod(nep,"NEPNLEIGSSetRKShifts_C",(NEP,PetscInt,PetscScalar*),(nep,ns,shifts));
1676: PetscFunctionReturn(PETSC_SUCCESS);
1677: }
1679: static PetscErrorCode NEPNLEIGSGetRKShifts_NLEIGS(NEP nep,PetscInt *ns,PetscScalar **shifts)
1680: {
1681: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1682: PetscInt i;
1684: PetscFunctionBegin;
1685: *ns = ctx->nshifts;
1686: if (ctx->nshifts) {
1687: PetscCall(PetscMalloc1(ctx->nshifts,shifts));
1688: for (i=0;i<ctx->nshifts;i++) (*shifts)[i] = ctx->shifts[i];
1689: }
1690: PetscFunctionReturn(PETSC_SUCCESS);
1691: }
1693: /*@C
1694: NEPNLEIGSGetRKShifts - Gets the list of shifts used in the Rational
1695: Krylov method.
1697: Not Collective
1699: Input Parameter:
1700: . nep - the nonlinear eigensolver context
1702: Output Parameters:
1703: + ns - number of shifts
1704: - shifts - array of shifts
1706: Note:
1707: The user is responsible for deallocating the returned array.
1709: Level: advanced
1711: .seealso: NEPNLEIGSSetRKShifts()
1712: @*/
1713: PetscErrorCode NEPNLEIGSGetRKShifts(NEP nep,PetscInt *ns,PetscScalar *shifts[])
1714: {
1715: PetscFunctionBegin;
1717: PetscAssertPointer(ns,2);
1718: PetscAssertPointer(shifts,3);
1719: PetscTryMethod(nep,"NEPNLEIGSGetRKShifts_C",(NEP,PetscInt*,PetscScalar**),(nep,ns,shifts));
1720: PetscFunctionReturn(PETSC_SUCCESS);
1721: }
1723: static PetscErrorCode NEPNLEIGSGetKSPs_NLEIGS(NEP nep,PetscInt *nsolve,KSP **ksp)
1724: {
1725: NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
1726: PetscInt i;
1727: PC pc;
1729: PetscFunctionBegin;
1730: if (!ctx->ksp) {
1731: PetscCall(NEPNLEIGSSetShifts(nep,&ctx->nshiftsw));
1732: PetscCall(PetscMalloc1(ctx->nshiftsw,&ctx->ksp));
1733: for (i=0;i<ctx->nshiftsw;i++) {
1734: PetscCall(KSPCreate(PetscObjectComm((PetscObject)nep),&ctx->ksp[i]));
1735: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ctx->ksp[i],(PetscObject)nep,1));
1736: PetscCall(KSPSetOptionsPrefix(ctx->ksp[i],((PetscObject)nep)->prefix));
1737: PetscCall(KSPAppendOptionsPrefix(ctx->ksp[i],"nep_nleigs_"));
1738: PetscCall(PetscObjectSetOptions((PetscObject)ctx->ksp[i],((PetscObject)nep)->options));
1739: PetscCall(KSPSetErrorIfNotConverged(ctx->ksp[i],PETSC_TRUE));
1740: PetscCall(KSPSetTolerances(ctx->ksp[i],1e-3*SlepcDefaultTol(nep->tol),PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT));
1741: PetscCall(KSPGetPC(ctx->ksp[i],&pc));
1742: if ((nep->fui==NEP_USER_INTERFACE_SPLIT && nep->P) || (nep->fui==NEP_USER_INTERFACE_CALLBACK && nep->function_pre!=nep->function)) {
1743: PetscCall(KSPSetType(ctx->ksp[i],KSPBCGS));
1744: PetscCall(PCSetType(pc,PCBJACOBI));
1745: } else {
1746: PetscCall(KSPSetType(ctx->ksp[i],KSPPREONLY));
1747: PetscCall(PCSetType(pc,PCLU));
1748: }
1749: }
1750: }
1751: if (nsolve) *nsolve = ctx->nshiftsw;
1752: if (ksp) *ksp = ctx->ksp;
1753: PetscFunctionReturn(PETSC_SUCCESS);
1754: }
1756: /*@C
1757: NEPNLEIGSGetKSPs - Retrieve the array of linear solver objects associated with
1758: the nonlinear eigenvalue solver.
1760: Collective
1762: Input Parameter:
1763: . nep - nonlinear eigenvalue solver
1765: Output Parameters:
1766: + nsolve - number of returned KSP objects
1767: - ksp - array of linear solver object
1769: Notes:
1770: The number of KSP objects is equal to the number of shifts provided by the user,
1771: or 1 if the user did not provide shifts.
1773: Level: advanced
1775: .seealso: NEPNLEIGSSetRKShifts()
1776: @*/
1777: PetscErrorCode NEPNLEIGSGetKSPs(NEP nep,PetscInt *nsolve,KSP **ksp)
1778: {
1779: PetscFunctionBegin;
1781: PetscUseMethod(nep,"NEPNLEIGSGetKSPs_C",(NEP,PetscInt*,KSP**),(nep,nsolve,ksp));
1782: PetscFunctionReturn(PETSC_SUCCESS);
1783: }
1785: static PetscErrorCode NEPNLEIGSSetFullBasis_NLEIGS(NEP nep,PetscBool fullbasis)
1786: {
1787: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1789: PetscFunctionBegin;
1790: if (fullbasis!=ctx->fullbasis) {
1791: ctx->fullbasis = fullbasis;
1792: nep->state = NEP_STATE_INITIAL;
1793: nep->useds = PetscNot(fullbasis);
1794: }
1795: PetscFunctionReturn(PETSC_SUCCESS);
1796: }
1798: /*@
1799: NEPNLEIGSSetFullBasis - Choose between TOAR-basis (default) and full-basis
1800: variants of the NLEIGS method.
1802: Logically Collective
1804: Input Parameters:
1805: + nep - the nonlinear eigensolver context
1806: - fullbasis - true if the full-basis variant must be selected
1808: Options Database Key:
1809: . -nep_nleigs_full_basis - Sets the full-basis flag
1811: Notes:
1812: The default is to use a compact representation of the Krylov basis, that is,
1813: V = (I otimes U) S, with a tensor BV. This behaviour can be changed so that
1814: the full basis V is explicitly stored and operated with. This variant is more
1815: expensive in terms of memory and computation, but is necessary in some cases,
1816: particularly for two-sided computations, see NEPSetTwoSided().
1818: In the full-basis variant, the NLEIGS solver uses an EPS object to explicitly
1819: solve the linearized eigenproblem, see NEPNLEIGSGetEPS().
1821: Level: advanced
1823: .seealso: NEPNLEIGSGetFullBasis(), NEPNLEIGSGetEPS(), NEPSetTwoSided(), BVCreateTensor()
1824: @*/
1825: PetscErrorCode NEPNLEIGSSetFullBasis(NEP nep,PetscBool fullbasis)
1826: {
1827: PetscFunctionBegin;
1830: PetscTryMethod(nep,"NEPNLEIGSSetFullBasis_C",(NEP,PetscBool),(nep,fullbasis));
1831: PetscFunctionReturn(PETSC_SUCCESS);
1832: }
1834: static PetscErrorCode NEPNLEIGSGetFullBasis_NLEIGS(NEP nep,PetscBool *fullbasis)
1835: {
1836: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1838: PetscFunctionBegin;
1839: *fullbasis = ctx->fullbasis;
1840: PetscFunctionReturn(PETSC_SUCCESS);
1841: }
1843: /*@
1844: NEPNLEIGSGetFullBasis - Gets the flag that indicates if NLEIGS is using the
1845: full-basis variant.
1847: Not Collective
1849: Input Parameter:
1850: . nep - the nonlinear eigensolver context
1852: Output Parameter:
1853: . fullbasis - the flag
1855: Level: advanced
1857: .seealso: NEPNLEIGSSetFullBasis()
1858: @*/
1859: PetscErrorCode NEPNLEIGSGetFullBasis(NEP nep,PetscBool *fullbasis)
1860: {
1861: PetscFunctionBegin;
1863: PetscAssertPointer(fullbasis,2);
1864: PetscUseMethod(nep,"NEPNLEIGSGetFullBasis_C",(NEP,PetscBool*),(nep,fullbasis));
1865: PetscFunctionReturn(PETSC_SUCCESS);
1866: }
1868: #define SHIFTMAX 30
1870: static PetscErrorCode NEPSetFromOptions_NLEIGS(NEP nep,PetscOptionItems *PetscOptionsObject)
1871: {
1872: NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
1873: PetscInt i=0,k;
1874: PetscBool flg1,flg2,b;
1875: PetscReal r;
1876: PetscScalar array[SHIFTMAX];
1878: PetscFunctionBegin;
1879: PetscOptionsHeadBegin(PetscOptionsObject,"NEP NLEIGS Options");
1881: PetscCall(PetscOptionsReal("-nep_nleigs_restart","Proportion of vectors kept after restart","NEPNLEIGSSetRestart",0.5,&r,&flg1));
1882: if (flg1) PetscCall(NEPNLEIGSSetRestart(nep,r));
1884: PetscCall(PetscOptionsBool("-nep_nleigs_locking","Choose between locking and non-locking variants","NEPNLEIGSSetLocking",PETSC_FALSE,&b,&flg1));
1885: if (flg1) PetscCall(NEPNLEIGSSetLocking(nep,b));
1887: PetscCall(PetscOptionsBool("-nep_nleigs_full_basis","Choose between TOAR and full-basis variants","NEPNLEIGSSetFullBasis",PETSC_FALSE,&b,&flg1));
1888: if (flg1) PetscCall(NEPNLEIGSSetFullBasis(nep,b));
1890: PetscCall(NEPNLEIGSGetInterpolation(nep,&r,&i));
1891: if (!i) i = PETSC_DEFAULT;
1892: PetscCall(PetscOptionsInt("-nep_nleigs_interpolation_degree","Maximum number of terms for interpolation via divided differences","NEPNLEIGSSetInterpolation",i,&i,&flg1));
1893: PetscCall(PetscOptionsReal("-nep_nleigs_interpolation_tol","Tolerance for interpolation via divided differences","NEPNLEIGSSetInterpolation",r,&r,&flg2));
1894: if (flg1 || flg2) PetscCall(NEPNLEIGSSetInterpolation(nep,r,i));
1896: k = SHIFTMAX;
1897: for (i=0;i<k;i++) array[i] = 0;
1898: PetscCall(PetscOptionsScalarArray("-nep_nleigs_rk_shifts","Shifts for Rational Krylov","NEPNLEIGSSetRKShifts",array,&k,&flg1));
1899: if (flg1) PetscCall(NEPNLEIGSSetRKShifts(nep,k,array));
1901: PetscOptionsHeadEnd();
1903: if (!ctx->ksp) PetscCall(NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp));
1904: for (i=0;i<ctx->nshiftsw;i++) PetscCall(KSPSetFromOptions(ctx->ksp[i]));
1906: if (ctx->fullbasis) {
1907: if (!ctx->eps) PetscCall(NEPNLEIGSGetEPS(nep,&ctx->eps));
1908: PetscCall(EPSSetFromOptions(ctx->eps));
1909: }
1910: PetscFunctionReturn(PETSC_SUCCESS);
1911: }
1913: static PetscErrorCode NEPView_NLEIGS(NEP nep,PetscViewer viewer)
1914: {
1915: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1916: PetscBool isascii;
1917: PetscInt i;
1918: char str[50];
1920: PetscFunctionBegin;
1921: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
1922: if (isascii) {
1923: PetscCall(PetscViewerASCIIPrintf(viewer," %d%% of basis vectors kept after restart\n",(int)(100*ctx->keep)));
1924: if (ctx->fullbasis) PetscCall(PetscViewerASCIIPrintf(viewer," using the full-basis variant\n"));
1925: else PetscCall(PetscViewerASCIIPrintf(viewer," using the %slocking variant\n",ctx->lock?"":"non-"));
1926: PetscCall(PetscViewerASCIIPrintf(viewer," divided difference terms: used=%" PetscInt_FMT ", max=%" PetscInt_FMT "\n",ctx->nmat,ctx->ddmaxit));
1927: PetscCall(PetscViewerASCIIPrintf(viewer," tolerance for divided difference convergence: %g\n",(double)ctx->ddtol));
1928: if (ctx->nshifts) {
1929: PetscCall(PetscViewerASCIIPrintf(viewer," RK shifts: "));
1930: PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE));
1931: for (i=0;i<ctx->nshifts;i++) {
1932: PetscCall(SlepcSNPrintfScalar(str,sizeof(str),ctx->shifts[i],PETSC_FALSE));
1933: PetscCall(PetscViewerASCIIPrintf(viewer,"%s%s",str,(i<ctx->nshifts-1)?",":""));
1934: }
1935: PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
1936: PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE));
1937: }
1938: if (!ctx->ksp) PetscCall(NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp));
1939: PetscCall(PetscViewerASCIIPushTab(viewer));
1940: PetscCall(KSPView(ctx->ksp[0],viewer));
1941: PetscCall(PetscViewerASCIIPopTab(viewer));
1942: if (ctx->fullbasis) {
1943: if (!ctx->eps) PetscCall(NEPNLEIGSGetEPS(nep,&ctx->eps));
1944: PetscCall(PetscViewerASCIIPushTab(viewer));
1945: PetscCall(EPSView(ctx->eps,viewer));
1946: PetscCall(PetscViewerASCIIPopTab(viewer));
1947: }
1948: }
1949: PetscFunctionReturn(PETSC_SUCCESS);
1950: }
1952: static PetscErrorCode NEPReset_NLEIGS(NEP nep)
1953: {
1954: PetscInt k;
1955: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1957: PetscFunctionBegin;
1958: if (nep->fui==NEP_USER_INTERFACE_SPLIT) PetscCall(PetscFree(ctx->coeffD));
1959: else {
1960: for (k=0;k<ctx->nmat;k++) PetscCall(MatDestroy(&ctx->D[k]));
1961: }
1962: PetscCall(PetscFree4(ctx->s,ctx->xi,ctx->beta,ctx->D));
1963: for (k=0;k<ctx->nshiftsw;k++) PetscCall(KSPReset(ctx->ksp[k]));
1964: PetscCall(VecDestroy(&ctx->vrn));
1965: if (ctx->fullbasis) {
1966: PetscCall(MatDestroy(&ctx->A));
1967: PetscCall(EPSReset(ctx->eps));
1968: for (k=0;k<4;k++) PetscCall(VecDestroy(&ctx->w[k]));
1969: }
1970: PetscFunctionReturn(PETSC_SUCCESS);
1971: }
1973: static PetscErrorCode NEPDestroy_NLEIGS(NEP nep)
1974: {
1975: PetscInt k;
1976: NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
1978: PetscFunctionBegin;
1979: PetscCall(BVDestroy(&ctx->V));
1980: for (k=0;k<ctx->nshiftsw;k++) PetscCall(KSPDestroy(&ctx->ksp[k]));
1981: PetscCall(PetscFree(ctx->ksp));
1982: if (ctx->nshifts) PetscCall(PetscFree(ctx->shifts));
1983: if (ctx->fullbasis) PetscCall(EPSDestroy(&ctx->eps));
1984: PetscCall(PetscFree(nep->data));
1985: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetSingularitiesFunction_C",NULL));
1986: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetSingularitiesFunction_C",NULL));
1987: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRestart_C",NULL));
1988: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRestart_C",NULL));
1989: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetLocking_C",NULL));
1990: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetLocking_C",NULL));
1991: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetInterpolation_C",NULL));
1992: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetInterpolation_C",NULL));
1993: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRKShifts_C",NULL));
1994: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRKShifts_C",NULL));
1995: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetKSPs_C",NULL));
1996: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetFullBasis_C",NULL));
1997: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetFullBasis_C",NULL));
1998: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetEPS_C",NULL));
1999: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetEPS_C",NULL));
2000: PetscFunctionReturn(PETSC_SUCCESS);
2001: }
2003: SLEPC_EXTERN PetscErrorCode NEPCreate_NLEIGS(NEP nep)
2004: {
2005: NEP_NLEIGS *ctx;
2007: PetscFunctionBegin;
2008: PetscCall(PetscNew(&ctx));
2009: nep->data = (void*)ctx;
2010: ctx->lock = PETSC_TRUE;
2011: ctx->ddtol = PETSC_DEFAULT;
2013: nep->useds = PETSC_TRUE;
2015: nep->ops->setup = NEPSetUp_NLEIGS;
2016: nep->ops->setfromoptions = NEPSetFromOptions_NLEIGS;
2017: nep->ops->view = NEPView_NLEIGS;
2018: nep->ops->destroy = NEPDestroy_NLEIGS;
2019: nep->ops->reset = NEPReset_NLEIGS;
2021: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetSingularitiesFunction_C",NEPNLEIGSSetSingularitiesFunction_NLEIGS));
2022: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetSingularitiesFunction_C",NEPNLEIGSGetSingularitiesFunction_NLEIGS));
2023: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRestart_C",NEPNLEIGSSetRestart_NLEIGS));
2024: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRestart_C",NEPNLEIGSGetRestart_NLEIGS));
2025: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetLocking_C",NEPNLEIGSSetLocking_NLEIGS));
2026: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetLocking_C",NEPNLEIGSGetLocking_NLEIGS));
2027: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetInterpolation_C",NEPNLEIGSSetInterpolation_NLEIGS));
2028: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetInterpolation_C",NEPNLEIGSGetInterpolation_NLEIGS));
2029: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRKShifts_C",NEPNLEIGSSetRKShifts_NLEIGS));
2030: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRKShifts_C",NEPNLEIGSGetRKShifts_NLEIGS));
2031: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetKSPs_C",NEPNLEIGSGetKSPs_NLEIGS));
2032: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetFullBasis_C",NEPNLEIGSSetFullBasis_NLEIGS));
2033: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetFullBasis_C",NEPNLEIGSGetFullBasis_NLEIGS));
2034: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetEPS_C",NEPNLEIGSSetEPS_NLEIGS));
2035: PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetEPS_C",NEPNLEIGSGetEPS_NLEIGS));
2036: PetscFunctionReturn(PETSC_SUCCESS);
2037: }