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