Actual source code: nleigs.c
slepc-3.13.3 2020-06-14
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2020, 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: norm = 0.0;
236: for (i=0;i<ndpt;i++) {
237: mean += F[i];
238: norm = PetscMax(PetscAbsScalar(F[i]),norm);
239: }
240: mean /= ndpt;
241: PetscBLASIntCast(ndpt,&lda_);
242: for (i=0;i<ndpt;i++) R[i] = PetscAbsScalar(F[i]-mean);
243: /* next support point */
244: err = 0.0;
245: for (i=0;i<ndpt;i++) if (R[i]>=err) {idx = i; err = R[i];}
246: for (k=0;k<ndpt-1;k++) {
247: z[k] = ds[idx]; f[k] = F[idx]; R[idx] = -1.0;
248: /* next column of Cauchy matrix */
249: for (i=0;i<ndpt;i++) {
250: C[i+k*ndpt] = 1.0/(ds[i]-ds[idx]);
251: }
253: PetscArrayzero(A,ndpt*ndpt);
254: cont = 0;
255: for (i=0;i<ndpt;i++) {
256: if (R[i]!=-1.0) {
257: for (j=0;j<=k;j++)A[cont+j*ndpt] = C[i+j*ndpt]*F[i]-C[i+j*ndpt]*f[j];
258: cont++;
259: }
260: }
261: PetscBLASIntCast(cont,&m_);
262: PetscBLASIntCast(k+1,&n_);
263: #if defined(PETSC_USE_COMPLEX)
264: PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","A",&m_,&n_,A,&lda_,S,NULL,&lda_,VT,&lda_,work,&lwork,rwork,&info));
265: #else
266: PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","A",&m_,&n_,A,&lda_,S,NULL,&lda_,VT,&lda_,work,&lwork,&info));
267: #endif
268: SlepcCheckLapackInfo("gesvd",info);
269: for (i=0;i<=k;i++) {
270: ww[i] = PetscConj(VT[i*ndpt+k]);
271: D[i] = ww[i]*f[i];
272: }
273: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&lda_,&n_,&sone,C,&lda_,D,&one,&szero,N,&one));
274: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&lda_,&n_,&sone,C,&lda_,ww,&one,&szero,D,&one));
275: for (i=0;i<ndpt;i++) if (R[i]>=0) R[i] = PetscAbsScalar(F[i]-N[i]/D[i]);
276: /* next support point */
277: err = 0.0;
278: for (i=0;i<ndpt;i++) if (R[i]>=err) {idx = i; err = R[i];}
279: if (err <= ctx->ddtol*norm) break;
280: }
282: if (k==ndpt-1) SETERRQ(PetscObjectComm((PetscObject)nep),1,"Failed to determine singularities automatically in general problem");
283: /* poles */
284: PetscArrayzero(C,ndpt*ndpt);
285: PetscArrayzero(A,ndpt*ndpt);
286: for (i=0;i<=k;i++) {
287: C[i+ndpt*i] = 1.0;
288: A[(i+1)*ndpt] = ww[i];
289: A[i+1] = 1.0;
290: A[i+1+(i+1)*ndpt] = z[i];
291: }
292: C[0] = 0.0; C[k+1+(k+1)*ndpt] = 1.0;
293: n_++;
294: #if defined(PETSC_USE_COMPLEX)
295: PetscStackCallBLAS("LAPACKggev",LAPACKggev_("N","N",&n_,A,&lda_,C,&lda_,D,N,NULL,&lda_,NULL,&lda_,work,&lwork,rwork,&info));
296: #else
297: PetscStackCallBLAS("LAPACKggev",LAPACKggev_("N","N",&n_,A,&lda_,C,&lda_,D,VT,N,NULL,&lda_,NULL,&lda_,work,&lwork,&info));
298: #endif
299: SlepcCheckLapackInfo("ggev",info);
300: cont = 0.0;
301: for (i=0;i<n_;i++) if (N[i]!=0.0) {
302: dxi[cont++] = D[i]/N[i];
303: }
304: *ndptx = cont;
305: PetscFree5(R,z,f,C,ww);
306: PetscFree6(A,S,VT,work,D,N);
307: #if defined(PETSC_USE_COMPLEX)
308: PetscFree(rwork);
309: #endif
310: return(0);
311: }
313: /* Singularities using Adaptive Anderson-Antoulas algorithm */
314: static PetscErrorCode NEPNLEIGSAAASingularities(NEP nep,PetscInt ndpt,PetscScalar *ds,PetscInt *ndptx,PetscScalar *dxi)
315: {
317: Vec u,v,w;
318: PetscRandom rand;
319: PetscScalar *F,*isol;
320: PetscInt i,k,nisol,nt;
321: Mat T;
322: FN f;
325: PetscMalloc1(ndpt,&F);
326: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
327: PetscMalloc1(ndpt,&isol);
328: *ndptx = 0;
329: NEPGetSplitOperatorInfo(nep,&nt,NULL);
330: nisol = *ndptx;
331: for (k=0;k<nt;k++) {
332: NEPGetSplitOperatorTerm(nep,k,NULL,&f);
333: for (i=0;i<ndpt;i++) {
334: FNEvaluateFunction(f,ds[i],&F[i]);
335: }
336: NEPNLEIGSAAAComputation(nep,ndpt,ds,F,&nisol,isol);
337: if (nisol) {
338: NEPNLEIGSAuxiliarRmDuplicates(nisol,isol,ndptx,dxi,ndpt);
339: }
340: }
341: PetscFree(isol);
342: } else {
343: MatCreateVecs(nep->function,&u,NULL);
344: VecDuplicate(u,&v);
345: VecDuplicate(u,&w);
346: BVGetRandomContext(nep->V,&rand);
347: VecSetRandom(u,rand);
348: VecNormalize(u,NULL);
349: VecSetRandom(v,rand);
350: VecNormalize(v,NULL);
351: T = nep->function;
352: for (i=0;i<ndpt;i++) {
353: NEPComputeFunction(nep,ds[i],T,T);
354: MatMult(T,v,w);
355: VecDot(w,u,&F[i]);
356: }
357: NEPNLEIGSAAAComputation(nep,ndpt,ds,F,ndptx,dxi);
358: VecDestroy(&u);
359: VecDestroy(&v);
360: VecDestroy(&w);
361: }
362: PetscFree(F);
363: return(0);
364: }
366: static PetscErrorCode NEPNLEIGSLejaBagbyPoints(NEP nep)
367: {
369: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
370: PetscInt i,k,ndpt=NDPOINTS,ndptx=NDPOINTS;
371: PetscScalar *ds,*dsi,*dxi,*nrs,*nrxi,*s=ctx->s,*xi=ctx->xi,*beta=ctx->beta;
372: PetscReal maxnrs,minnrxi;
373: PetscBool rational;
374: #if !defined(PETSC_USE_COMPLEX)
375: PetscReal a,b,h;
376: #endif
379: if (!ctx->computesingularities && nep->problem_type!=NEP_RATIONAL) ndpt = ndptx = LBPOINTS;
380: PetscMalloc5(ndpt+1,&ds,ndpt+1,&dsi,ndpt,&dxi,ndpt+1,&nrs,ndpt,&nrxi);
382: /* Discretize the target region boundary */
383: RGComputeContour(nep->rg,ndpt,ds,dsi);
384: #if !defined(PETSC_USE_COMPLEX)
385: for (i=0;i<ndpt;i++) if (dsi[i]!=0.0) break;
386: if (i<ndpt) {
387: if (nep->problem_type==NEP_RATIONAL) {
388: /* Select a segment in the real axis */
389: RGComputeBoundingBox(nep->rg,&a,&b,NULL,NULL);
390: if (a<=-PETSC_MAX_REAL || b>=PETSC_MAX_REAL) SETERRQ(PetscObjectComm((PetscObject)nep),1,"NLEIGS requires a bounded target set");
391: h = (b-a)/ndpt;
392: for (i=0;i<ndpt;i++) {ds[i] = a+h*i; dsi[i] = 0.0;}
393: } else SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"NLEIGS with real arithmetic requires the target set to be included in the real axis");
394: }
395: #endif
396: /* Discretize the singularity region */
397: if (ctx->computesingularities) {
398: (ctx->computesingularities)(nep,&ndptx,dxi,ctx->singularitiesctx);
399: } else {
400: if (nep->problem_type==NEP_RATIONAL) {
401: NEPNLEIGSRationalSingularities(nep,&ndptx,dxi,&rational);
402: if (!rational) SETERRQ(PetscObjectComm((PetscObject)nep),1,"Failed to determine singularities automatically in rational problem; consider solving the problem as general");
403: } else {
404: /* AAA algorithm */
405: NEPNLEIGSAAASingularities(nep,ndpt,ds,&ndptx,dxi);
406: }
407: }
408: /* Look for Leja-Bagby points in the discretization sets */
409: s[0] = ds[0];
410: xi[0] = (ndptx>0)?dxi[0]:PETSC_INFINITY;
411: 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]));
412: beta[0] = 1.0; /* scaling factors are also computed here */
413: for (i=0;i<ndpt;i++) {
414: nrs[i] = 1.0;
415: nrxi[i] = 1.0;
416: }
417: for (k=1;k<ctx->ddmaxit;k++) {
418: maxnrs = 0.0;
419: minnrxi = PETSC_MAX_REAL;
420: for (i=0;i<ndpt;i++) {
421: nrs[i] *= ((ds[i]-s[k-1])/(1.0-ds[i]/xi[k-1]))/beta[k-1];
422: if (PetscAbsScalar(nrs[i])>maxnrs) {maxnrs = PetscAbsScalar(nrs[i]); s[k] = ds[i];}
423: }
424: if (ndptx>k) {
425: for (i=1;i<ndptx;i++) {
426: nrxi[i] *= ((dxi[i]-s[k-1])/(1.0-dxi[i]/xi[k-1]))/beta[k-1];
427: if (PetscAbsScalar(nrxi[i])<minnrxi) {minnrxi = PetscAbsScalar(nrxi[i]); xi[k] = dxi[i];}
428: }
429: 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]));
430: } else xi[k] = PETSC_INFINITY;
431: beta[k] = maxnrs;
432: }
433: PetscFree5(ds,dsi,dxi,nrs,nrxi);
434: return(0);
435: }
437: PetscErrorCode NEPNLEIGSEvalNRTFunct(NEP nep,PetscInt k,PetscScalar sigma,PetscScalar *b)
438: {
439: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
440: PetscInt i;
441: PetscScalar *beta=ctx->beta,*s=ctx->s,*xi=ctx->xi;
444: b[0] = 1.0/beta[0];
445: for (i=0;i<k;i++) {
446: b[i+1] = ((sigma-s[i])*b[i])/(beta[i+1]*(1.0-sigma/xi[i]));
447: }
448: return(0);
449: }
451: static PetscErrorCode MatMult_Fun(Mat A,Vec x,Vec y)
452: {
453: PetscErrorCode ierr;
454: NEP_NLEIGS_MATSHELL *ctx;
455: PetscInt i;
458: MatShellGetContext(A,(void**)&ctx);
459: MatMult(ctx->A[0],x,y);
460: if (ctx->coeff[0]!=1.0) { VecScale(y,ctx->coeff[0]); }
461: for (i=1;i<ctx->nmat;i++) {
462: MatMult(ctx->A[i],x,ctx->t);
463: VecAXPY(y,ctx->coeff[i],ctx->t);
464: }
465: return(0);
466: }
468: static PetscErrorCode MatMultTranspose_Fun(Mat A,Vec x,Vec y)
469: {
470: PetscErrorCode ierr;
471: NEP_NLEIGS_MATSHELL *ctx;
472: PetscInt i;
475: MatShellGetContext(A,(void**)&ctx);
476: MatMultTranspose(ctx->A[0],x,y);
477: if (ctx->coeff[0]!=1.0) { VecScale(y,ctx->coeff[0]); }
478: for (i=1;i<ctx->nmat;i++) {
479: MatMultTranspose(ctx->A[i],x,ctx->t);
480: VecAXPY(y,ctx->coeff[i],ctx->t);
481: }
482: return(0);
483: }
485: static PetscErrorCode MatGetDiagonal_Fun(Mat A,Vec diag)
486: {
487: PetscErrorCode ierr;
488: NEP_NLEIGS_MATSHELL *ctx;
489: PetscInt i;
492: MatShellGetContext(A,(void**)&ctx);
493: MatGetDiagonal(ctx->A[0],diag);
494: if (ctx->coeff[0]!=1.0) { VecScale(diag,ctx->coeff[0]); }
495: for (i=1;i<ctx->nmat;i++) {
496: MatGetDiagonal(ctx->A[i],ctx->t);
497: VecAXPY(diag,ctx->coeff[i],ctx->t);
498: }
499: return(0);
500: }
502: static PetscErrorCode MatDuplicate_Fun(Mat A,MatDuplicateOption op,Mat *B)
503: {
504: PetscInt n,i;
505: NEP_NLEIGS_MATSHELL *ctxnew,*ctx;
506: void (*fun)(void);
507: PetscErrorCode ierr;
510: MatShellGetContext(A,(void**)&ctx);
511: PetscNew(&ctxnew);
512: ctxnew->nmat = ctx->nmat;
513: ctxnew->maxnmat = ctx->maxnmat;
514: PetscMalloc2(ctxnew->maxnmat,&ctxnew->A,ctxnew->maxnmat,&ctxnew->coeff);
515: for (i=0;i<ctx->nmat;i++) {
516: PetscObjectReference((PetscObject)ctx->A[i]);
517: ctxnew->A[i] = ctx->A[i];
518: ctxnew->coeff[i] = ctx->coeff[i];
519: }
520: MatGetSize(ctx->A[0],&n,NULL);
521: VecDuplicate(ctx->t,&ctxnew->t);
522: MatCreateShell(PetscObjectComm((PetscObject)A),n,n,n,n,(void*)ctxnew,B);
523: MatShellSetManageScalingShifts(*B);
524: MatShellGetOperation(A,MATOP_MULT,&fun);
525: MatShellSetOperation(*B,MATOP_MULT,fun);
526: MatShellGetOperation(A,MATOP_MULT_TRANSPOSE,&fun);
527: MatShellSetOperation(*B,MATOP_MULT_TRANSPOSE,fun);
528: MatShellGetOperation(A,MATOP_GET_DIAGONAL,&fun);
529: MatShellSetOperation(*B,MATOP_GET_DIAGONAL,fun);
530: MatShellGetOperation(A,MATOP_DUPLICATE,&fun);
531: MatShellSetOperation(*B,MATOP_DUPLICATE,fun);
532: MatShellGetOperation(A,MATOP_DESTROY,&fun);
533: MatShellSetOperation(*B,MATOP_DESTROY,fun);
534: MatShellGetOperation(A,MATOP_AXPY,&fun);
535: MatShellSetOperation(*B,MATOP_AXPY,fun);
536: return(0);
537: }
539: static PetscErrorCode MatDestroy_Fun(Mat A)
540: {
541: NEP_NLEIGS_MATSHELL *ctx;
542: PetscErrorCode ierr;
543: PetscInt i;
546: if (A) {
547: MatShellGetContext(A,(void**)&ctx);
548: for (i=0;i<ctx->nmat;i++) {
549: MatDestroy(&ctx->A[i]);
550: }
551: VecDestroy(&ctx->t);
552: PetscFree2(ctx->A,ctx->coeff);
553: PetscFree(ctx);
554: }
555: return(0);
556: }
558: static PetscErrorCode MatAXPY_Fun(Mat Y,PetscScalar a,Mat X,MatStructure str)
559: {
560: NEP_NLEIGS_MATSHELL *ctxY,*ctxX;
561: PetscErrorCode ierr;
562: PetscInt i,j;
563: PetscBool found;
566: MatShellGetContext(Y,(void**)&ctxY);
567: MatShellGetContext(X,(void**)&ctxX);
568: for (i=0;i<ctxX->nmat;i++) {
569: found = PETSC_FALSE;
570: for (j=0;!found&&j<ctxY->nmat;j++) {
571: if (ctxX->A[i]==ctxY->A[j]) {
572: found = PETSC_TRUE;
573: ctxY->coeff[j] += a*ctxX->coeff[i];
574: }
575: }
576: if (!found) {
577: ctxY->coeff[ctxY->nmat] = a*ctxX->coeff[i];
578: ctxY->A[ctxY->nmat++] = ctxX->A[i];
579: PetscObjectReference((PetscObject)ctxX->A[i]);
580: }
581: }
582: return(0);
583: }
585: static PetscErrorCode MatScale_Fun(Mat M,PetscScalar a)
586: {
587: NEP_NLEIGS_MATSHELL *ctx;
588: PetscErrorCode ierr;
589: PetscInt i;
592: MatShellGetContext(M,(void**)&ctx);
593: for (i=0;i<ctx->nmat;i++) ctx->coeff[i] *= a;
594: return(0);
595: }
597: static PetscErrorCode NLEIGSMatToMatShellArray(Mat M,Mat *Ms,PetscInt maxnmat)
598: {
599: PetscErrorCode ierr;
600: NEP_NLEIGS_MATSHELL *ctx;
601: PetscInt n;
602: PetscBool has;
605: MatHasOperation(M,MATOP_DUPLICATE,&has);
606: if (!has) SETERRQ(PetscObjectComm((PetscObject)M),1,"MatDuplicate operation required");
607: PetscNew(&ctx);
608: ctx->maxnmat = maxnmat;
609: PetscMalloc2(ctx->maxnmat,&ctx->A,ctx->maxnmat,&ctx->coeff);
610: MatDuplicate(M,MAT_COPY_VALUES,&ctx->A[0]);
611: ctx->nmat = 1;
612: ctx->coeff[0] = 1.0;
613: MatCreateVecs(M,&ctx->t,NULL);
614: MatGetSize(M,&n,NULL);
615: MatCreateShell(PetscObjectComm((PetscObject)M),n,n,n,n,(void*)ctx,Ms);
616: MatShellSetManageScalingShifts(*Ms);
617: MatShellSetOperation(*Ms,MATOP_MULT,(void(*)(void))MatMult_Fun);
618: MatShellSetOperation(*Ms,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_Fun);
619: MatShellSetOperation(*Ms,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Fun);
620: MatShellSetOperation(*Ms,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_Fun);
621: MatShellSetOperation(*Ms,MATOP_DESTROY,(void(*)(void))MatDestroy_Fun);
622: MatShellSetOperation(*Ms,MATOP_AXPY,(void(*)(void))MatAXPY_Fun);
623: MatShellSetOperation(*Ms,MATOP_SCALE,(void(*)(void))MatScale_Fun);
624: return(0);
625: }
627: static PetscErrorCode NEPNLEIGSNormEstimation(NEP nep,Mat M,PetscReal *norm,Vec *w)
628: {
629: PetscScalar *z,*x,*y;
630: PetscReal tr;
631: Vec X=w[0],Y=w[1];
632: PetscInt n,i;
633: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
634: PetscRandom rand;
638: if (!ctx->vrn) {
639: /* generate a random vector with normally distributed entries with the Box-Muller transform */
640: BVGetRandomContext(nep->V,&rand);
641: MatCreateVecs(M,&ctx->vrn,NULL);
642: VecSetRandom(X,rand);
643: VecSetRandom(Y,rand);
644: VecGetLocalSize(ctx->vrn,&n);
645: VecGetArray(ctx->vrn,&z);
646: VecGetArray(X,&x);
647: VecGetArray(Y,&y);
648: for (i=0;i<n;i++) {
649: #if defined(PETSC_USE_COMPLEX)
650: z[i] = PetscCMPLX(PetscSqrtReal(-2.0*PetscLogReal(PetscRealPart(x[i])))*PetscCosReal(2.0*PETSC_PI*PetscRealPart(y[i])),PetscSqrtReal(-2.0*PetscLogReal(PetscImaginaryPart(x[i])))*PetscCosReal(2.0*PETSC_PI*PetscImaginaryPart(y[i])));
651: #else
652: z[i] = PetscSqrtReal(-2.0*PetscLogReal(x[i]))*PetscCosReal(2.0*PETSC_PI*y[i]);
653: #endif
654: }
655: VecRestoreArray(ctx->vrn,&z);
656: VecRestoreArray(X,&x);
657: VecRestoreArray(Y,&y);
658: VecNorm(ctx->vrn,NORM_2,&tr);
659: VecScale(ctx->vrn,1/tr);
660: }
661: /* matrix-free norm estimator of Ipsen https://ipsen.math.ncsu.edu/ps/slides_ima.pdf */
662: MatGetSize(M,&n,NULL);
663: MatMult(M,ctx->vrn,X);
664: VecNorm(X,NORM_2,norm);
665: *norm *= PetscSqrtReal((PetscReal)n);
666: return(0);
667: }
669: static PetscErrorCode NEPNLEIGSDividedDifferences_split(NEP nep)
670: {
672: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
673: PetscInt k,j,i,maxnmat,nmax;
674: PetscReal norm0,norm,*matnorm;
675: PetscScalar *s=ctx->s,*beta=ctx->beta,*xi=ctx->xi,*b,alpha,*coeffs,*pK,*pH,sone=1.0;
676: Mat T,Ts,K,H;
677: PetscBool shell,hasmnorm=PETSC_FALSE,matrix=PETSC_TRUE;
678: PetscBLASInt n_;
681: nmax = ctx->ddmaxit;
682: PetscMalloc1(nep->nt*nmax,&ctx->coeffD);
683: PetscMalloc3(nmax+1,&b,nmax+1,&coeffs,nep->nt,&matnorm);
684: for (j=0;j<nep->nt;j++) {
685: MatHasOperation(nep->A[j],MATOP_NORM,&hasmnorm);
686: if (!hasmnorm) break;
687: MatNorm(nep->A[j],NORM_INFINITY,matnorm+j);
688: }
689: /* Try matrix functions scheme */
690: PetscCalloc2(nmax*nmax,&pK,nmax*nmax,&pH);
691: for (i=0;i<nmax-1;i++) {
692: pK[(nmax+1)*i] = 1.0;
693: pK[(nmax+1)*i+1] = beta[i+1]/xi[i];
694: pH[(nmax+1)*i] = s[i];
695: pH[(nmax+1)*i+1] = beta[i+1];
696: }
697: pH[nmax*nmax-1] = s[nmax-1];
698: pK[nmax*nmax-1] = 1.0;
699: PetscBLASIntCast(nmax,&n_);
700: PetscStackCallBLAS("BLAStrsm",BLAStrsm_("R","L","N","U",&n_,&n_,&sone,pK,&n_,pH,&n_));
701: /* The matrix to be used is in H. K will be a work-space matrix */
702: MatCreateSeqDense(PETSC_COMM_SELF,nmax,nmax,pH,&H);
703: MatCreateSeqDense(PETSC_COMM_SELF,nmax,nmax,pK,&K);
704: for (j=0;matrix&&j<nep->nt;j++) {
705: PetscPushErrorHandler(PetscIgnoreErrorHandler,NULL);
706: FNEvaluateFunctionMat(nep->f[j],H,K);
707: PetscPopErrorHandler();
708: if (!ierr) {
709: for (i=0;i<nmax;i++) { ctx->coeffD[j+i*nep->nt] = pK[i]*beta[0]; }
710: } else {
711: matrix = PETSC_FALSE;
712: PetscFPTrapPop();
713: }
714: }
715: MatDestroy(&H);
716: MatDestroy(&K);
717: if (!matrix) {
718: for (j=0;j<nep->nt;j++) {
719: FNEvaluateFunction(nep->f[j],s[0],ctx->coeffD+j);
720: ctx->coeffD[j] *= beta[0];
721: }
722: }
723: if (hasmnorm) {
724: norm0 = 0.0;
725: for (j=0;j<nep->nt;j++) norm0 += matnorm[j]*PetscAbsScalar(ctx->coeffD[j]);
726: } else {
727: norm0 = 0.0;
728: for (j=0;j<nep->nt;j++) norm0 = PetscMax(PetscAbsScalar(ctx->coeffD[j]),norm0);
729: }
730: ctx->nmat = ctx->ddmaxit;
731: for (k=1;k<ctx->ddmaxit;k++) {
732: if (!matrix) {
733: NEPNLEIGSEvalNRTFunct(nep,k,s[k],b);
734: for (i=0;i<nep->nt;i++) {
735: FNEvaluateFunction(nep->f[i],s[k],ctx->coeffD+k*nep->nt+i);
736: for (j=0;j<k;j++) {
737: ctx->coeffD[k*nep->nt+i] -= b[j]*ctx->coeffD[i+nep->nt*j];
738: }
739: ctx->coeffD[k*nep->nt+i] /= b[k];
740: }
741: }
742: if (hasmnorm) {
743: norm = 0.0;
744: for (j=0;j<nep->nt;j++) norm += matnorm[j]*PetscAbsScalar(ctx->coeffD[k*nep->nt+j]);
745: } else {
746: norm = 0.0;
747: for (j=0;j<nep->nt;j++) norm = PetscMax(PetscAbsScalar(ctx->coeffD[k*nep->nt+j]),norm);
748: }
749: if (k>1 && norm/norm0 < ctx->ddtol) {
750: ctx->nmat = k+1;
751: break;
752: }
753: }
754: if (!ctx->ksp) { NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp); }
755: PetscObjectTypeCompare((PetscObject)nep->A[0],MATSHELL,&shell);
756: maxnmat = PetscMax(ctx->ddmaxit,nep->nt);
757: for (i=0;i<ctx->nshiftsw;i++) {
758: NEPNLEIGSEvalNRTFunct(nep,ctx->nmat-1,ctx->shifts[i],coeffs);
759: if (!shell) {
760: MatDuplicate(nep->A[0],MAT_COPY_VALUES,&T);
761: } else {
762: NLEIGSMatToMatShellArray(nep->A[0],&T,maxnmat);
763: }
764: alpha = 0.0;
765: for (j=0;j<ctx->nmat;j++) alpha += coeffs[j]*ctx->coeffD[j*nep->nt];
766: MatScale(T,alpha);
767: for (k=1;k<nep->nt;k++) {
768: alpha = 0.0;
769: for (j=0;j<ctx->nmat;j++) alpha += coeffs[j]*ctx->coeffD[j*nep->nt+k];
770: if (shell) {
771: NLEIGSMatToMatShellArray(nep->A[k],&Ts,maxnmat);
772: }
773: MatAXPY(T,alpha,shell?Ts:nep->A[k],nep->mstr);
774: if (shell) {
775: MatDestroy(&Ts);
776: }
777: }
778: KSPSetOperators(ctx->ksp[i],T,T);
779: KSPSetUp(ctx->ksp[i]);
780: MatDestroy(&T);
781: }
782: PetscFree3(b,coeffs,matnorm);
783: PetscFree2(pK,pH);
784: return(0);
785: }
787: static PetscErrorCode NEPNLEIGSDividedDifferences_callback(NEP nep)
788: {
790: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
791: PetscInt k,j,i,maxnmat;
792: PetscReal norm0,norm;
793: PetscScalar *s=ctx->s,*beta=ctx->beta,*b,*coeffs;
794: Mat *D=ctx->D,T;
795: PetscBool shell,has,vec=PETSC_FALSE;
796: Vec w[2];
799: PetscMalloc2(ctx->ddmaxit+1,&b,ctx->ddmaxit+1,&coeffs);
800: T = nep->function;
801: NEPComputeFunction(nep,s[0],T,T);
802: PetscObjectTypeCompare((PetscObject)T,MATSHELL,&shell);
803: maxnmat = PetscMax(ctx->ddmaxit,nep->nt);
804: if (!shell) {
805: MatDuplicate(T,MAT_COPY_VALUES,&D[0]);
806: } else {
807: NLEIGSMatToMatShellArray(T,&D[0],maxnmat);
808: }
809: if (beta[0]!=1.0) {
810: MatScale(D[0],1.0/beta[0]);
811: }
812: MatHasOperation(D[0],MATOP_NORM,&has);
813: if (has) {
814: MatNorm(D[0],NORM_FROBENIUS,&norm0);
815: } else {
816: MatCreateVecs(D[0],NULL,&w[0]);
817: VecDuplicate(w[0],&w[1]);
818: vec = PETSC_TRUE;
819: NEPNLEIGSNormEstimation(nep,D[0],&norm0,w);
820: }
821: ctx->nmat = ctx->ddmaxit;
822: for (k=1;k<ctx->ddmaxit;k++) {
823: NEPNLEIGSEvalNRTFunct(nep,k,s[k],b);
824: NEPComputeFunction(nep,s[k],T,T);
825: if (!shell) {
826: MatDuplicate(T,MAT_COPY_VALUES,&D[k]);
827: } else {
828: NLEIGSMatToMatShellArray(T,&D[k],maxnmat);
829: }
830: for (j=0;j<k;j++) {
831: MatAXPY(D[k],-b[j],D[j],nep->mstr);
832: }
833: MatScale(D[k],1.0/b[k]);
834: MatHasOperation(D[k],MATOP_NORM,&has);
835: if (has) {
836: MatNorm(D[k],NORM_FROBENIUS,&norm);
837: } else {
838: if (!vec) {
839: MatCreateVecs(D[k],NULL,&w[0]);
840: VecDuplicate(w[0],&w[1]);
841: vec = PETSC_TRUE;
842: }
843: NEPNLEIGSNormEstimation(nep,D[k],&norm,w);
844: }
845: if (k>1 && norm/norm0 < ctx->ddtol && k>1) {
846: ctx->nmat = k+1;
847: break;
848: }
849: }
850: if (!ctx->ksp) { NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp); }
851: for (i=0;i<ctx->nshiftsw;i++) {
852: NEPNLEIGSEvalNRTFunct(nep,ctx->nmat-1,ctx->shifts[i],coeffs);
853: MatDuplicate(ctx->D[0],MAT_COPY_VALUES,&T);
854: if (coeffs[0]!=1.0) { MatScale(T,coeffs[0]); }
855: for (j=1;j<ctx->nmat;j++) {
856: MatAXPY(T,coeffs[j],ctx->D[j],nep->mstr);
857: }
858: KSPSetOperators(ctx->ksp[i],T,T);
859: KSPSetUp(ctx->ksp[i]);
860: MatDestroy(&T);
861: }
862: PetscFree2(b,coeffs);
863: if (vec) {
864: VecDestroy(&w[0]);
865: VecDestroy(&w[1]);
866: }
867: return(0);
868: }
870: /*
871: NEPKrylovConvergence - This is the analogue to EPSKrylovConvergence.
872: */
873: PetscErrorCode NEPNLEIGSKrylovConvergence(NEP nep,PetscBool getall,PetscInt kini,PetscInt nits,PetscReal betah,PetscScalar betak,PetscInt *kout,Vec *w)
874: {
876: PetscInt k,newk,marker,inside;
877: PetscScalar re,im;
878: PetscReal resnorm,tt;
879: PetscBool istrivial;
880: NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
883: RGIsTrivial(nep->rg,&istrivial);
884: marker = -1;
885: if (nep->trackall) getall = PETSC_TRUE;
886: for (k=kini;k<kini+nits;k++) {
887: /* eigenvalue */
888: re = nep->eigr[k];
889: im = nep->eigi[k];
890: if (!istrivial) {
891: if (!ctx->nshifts) {
892: NEPNLEIGSBackTransform((PetscObject)nep,1,&re,&im);
893: }
894: RGCheckInside(nep->rg,1,&re,&im,&inside);
895: if (marker==-1 && inside<0) marker = k;
896: }
897: newk = k;
898: DSVectors(nep->ds,DS_MAT_X,&newk,&resnorm);
899: tt = (ctx->nshifts)?SlepcAbsEigenvalue(betak-nep->eigr[k]*betah,nep->eigi[k]*betah):betah;
900: resnorm *= PetscAbsReal(tt);
901: /* error estimate */
902: (*nep->converged)(nep,nep->eigr[k],nep->eigi[k],resnorm,&nep->errest[k],nep->convergedctx);
903: if (marker==-1 && nep->errest[k] >= nep->tol) marker = k;
904: if (newk==k+1) {
905: nep->errest[k+1] = nep->errest[k];
906: k++;
907: }
908: if (marker!=-1 && !getall) break;
909: }
910: if (marker!=-1) k = marker;
911: *kout = k;
912: return(0);
913: }
915: PetscErrorCode NEPSetUp_NLEIGS(NEP nep)
916: {
918: PetscInt k,in;
919: PetscScalar zero=0.0;
920: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
921: SlepcSC sc;
922: PetscBool istrivial;
925: NEPSetDimensions_Default(nep,nep->nev,&nep->ncv,&nep->mpd);
926: if (nep->ncv>nep->nev+nep->mpd) SETERRQ(PetscObjectComm((PetscObject)nep),1,"The value of ncv must not be larger than nev+mpd");
927: if (nep->max_it==PETSC_DEFAULT) nep->max_it = PetscMax(5000,2*nep->n/nep->ncv);
928: if (!ctx->ddmaxit) ctx->ddmaxit = LBPOINTS;
929: RGIsTrivial(nep->rg,&istrivial);
930: if (istrivial) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"NEPNLEIGS requires a nontrivial region defining the target set");
931: if (!nep->which) nep->which = NEP_TARGET_MAGNITUDE;
932: 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,"Should set a target selection in NEPSetWhichEigenvalues()");
934: /* Initialize the NLEIGS context structure */
935: k = ctx->ddmaxit;
936: PetscMalloc4(k,&ctx->s,k,&ctx->xi,k,&ctx->beta,k,&ctx->D);
937: nep->data = ctx;
938: if (nep->tol==PETSC_DEFAULT) nep->tol = SLEPC_DEFAULT_TOL;
939: if (ctx->ddtol==PETSC_DEFAULT) ctx->ddtol = nep->tol/10.0;
940: if (!ctx->keep) ctx->keep = 0.5;
942: /* Compute Leja-Bagby points and scaling values */
943: NEPNLEIGSLejaBagbyPoints(nep);
944: if (nep->problem_type!=NEP_RATIONAL) {
945: RGCheckInside(nep->rg,1,&nep->target,&zero,&in);
946: if (in<0) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"The target is not inside the target set");
947: }
949: /* Compute the divided difference matrices */
950: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
951: NEPNLEIGSDividedDifferences_split(nep);
952: } else {
953: NEPNLEIGSDividedDifferences_callback(nep);
954: }
955: NEPAllocateSolution(nep,ctx->nmat-1);
956: NEPSetWorkVecs(nep,4);
957: if (!ctx->fullbasis) {
958: if (nep->twosided) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Two-sided variant requires the full-basis option, rerun with -nep_nleigs_full_basis");
959: /* set-up DS and transfer split operator functions */
960: DSSetType(nep->ds,ctx->nshifts?DSGNHEP:DSNHEP);
961: DSAllocate(nep->ds,nep->ncv+1);
962: DSGetSlepcSC(nep->ds,&sc);
963: if (!ctx->nshifts) {
964: sc->map = NEPNLEIGSBackTransform;
965: DSSetExtraRow(nep->ds,PETSC_TRUE);
966: }
967: sc->mapobj = (PetscObject)nep;
968: sc->rg = nep->rg;
969: sc->comparison = nep->sc->comparison;
970: sc->comparisonctx = nep->sc->comparisonctx;
971: BVDestroy(&ctx->V);
972: BVCreateTensor(nep->V,ctx->nmat-1,&ctx->V);
973: nep->ops->solve = NEPSolve_NLEIGS;
974: nep->ops->computevectors = NEPComputeVectors_Schur;
975: } else {
976: NEPSetUp_NLEIGS_FullBasis(nep);
977: nep->ops->solve = NEPSolve_NLEIGS_FullBasis;
978: nep->ops->computevectors = NULL;
979: }
980: return(0);
981: }
983: /*
984: Extend the TOAR basis by applying the the matrix operator
985: over a vector which is decomposed on the TOAR way
986: Input:
987: - S,V: define the latest Arnoldi vector (nv vectors in V)
988: Output:
989: - t: new vector extending the TOAR basis
990: - r: temporally coefficients to compute the TOAR coefficients
991: for the new Arnoldi vector
992: Workspace: t_ (two vectors)
993: */
994: 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_)
995: {
997: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
998: PetscInt deg=ctx->nmat-1,k,j;
999: Vec v=t_[0],q=t_[1],w;
1000: PetscScalar *beta=ctx->beta,*s=ctx->s,*xi=ctx->xi,*coeffs,sigma;
1003: if (!ctx->ksp) { NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp); }
1004: sigma = ctx->shifts[idxrktg];
1005: BVSetActiveColumns(nep->V,0,nv);
1006: PetscMalloc1(ctx->nmat,&coeffs);
1007: if (PetscAbsScalar(s[deg-2]-sigma)<100*PETSC_MACHINE_EPSILON) SETERRQ(PETSC_COMM_SELF,1,"Breakdown in NLEIGS");
1008: /* i-part stored in (i-1) position */
1009: for (j=0;j<nv;j++) {
1010: 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);
1011: }
1012: BVSetActiveColumns(W,0,deg);
1013: BVGetColumn(W,deg-1,&w);
1014: BVMultVec(V,1.0/beta[deg],0,w,S+(deg-1)*ls);
1015: BVRestoreColumn(W,deg-1,&w);
1016: BVGetColumn(W,deg-2,&w);
1017: BVMultVec(V,1.0,0.0,w,r+(deg-2)*lr);
1018: BVRestoreColumn(W,deg-2,&w);
1019: for (k=deg-2;k>0;k--) {
1020: if (PetscAbsScalar(s[k-1]-sigma)<100*PETSC_MACHINE_EPSILON) SETERRQ(PETSC_COMM_SELF,1,"Breakdown in NLEIGS");
1021: 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);
1022: BVGetColumn(W,k-1,&w);
1023: BVMultVec(V,1.0,0.0,w,r+(k-1)*lr);
1024: BVRestoreColumn(W,k-1,&w);
1025: }
1026: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
1027: for (j=0;j<ctx->nmat-2;j++) coeffs[j] = ctx->coeffD[nep->nt*j];
1028: coeffs[ctx->nmat-2] = ctx->coeffD[nep->nt*(ctx->nmat-1)];
1029: BVMultVec(W,1.0,0.0,v,coeffs);
1030: MatMult(nep->A[0],v,q);
1031: for (k=1;k<nep->nt;k++) {
1032: for (j=0;j<ctx->nmat-2;j++) coeffs[j] = ctx->coeffD[nep->nt*j+k];
1033: coeffs[ctx->nmat-2] = ctx->coeffD[nep->nt*(ctx->nmat-1)+k];
1034: BVMultVec(W,1.0,0,v,coeffs);
1035: MatMult(nep->A[k],v,t);
1036: VecAXPY(q,1.0,t);
1037: }
1038: KSPSolve(ctx->ksp[idxrktg],q,t);
1039: VecScale(t,-1.0);
1040: } else {
1041: for (k=0;k<deg-1;k++) {
1042: BVGetColumn(W,k,&w);
1043: MatMult(ctx->D[k],w,q);
1044: BVRestoreColumn(W,k,&w);
1045: BVInsertVec(W,k,q);
1046: }
1047: BVGetColumn(W,deg-1,&w);
1048: MatMult(ctx->D[deg],w,q);
1049: BVRestoreColumn(W,k,&w);
1050: BVInsertVec(W,k,q);
1051: for (j=0;j<ctx->nmat-1;j++) coeffs[j] = 1.0;
1052: BVMultVec(W,1.0,0.0,q,coeffs);
1053: KSPSolve(ctx->ksp[idxrktg],q,t);
1054: VecScale(t,-1.0);
1055: }
1056: PetscFree(coeffs);
1057: return(0);
1058: }
1060: /*
1061: Compute TOAR coefficients of the blocks of the new Arnoldi vector computed
1062: */
1063: static PetscErrorCode NEPTOARCoefficients(NEP nep,PetscScalar sigma,PetscInt nv,PetscScalar *S,PetscInt ls,PetscScalar *r,PetscInt lr,PetscScalar *x,PetscScalar *work)
1064: {
1066: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1067: PetscInt k,j,d=ctx->nmat-1;
1068: PetscScalar *t=work;
1071: NEPNLEIGSEvalNRTFunct(nep,d-1,sigma,t);
1072: for (k=0;k<d-1;k++) {
1073: for (j=0;j<=nv;j++) r[k*lr+j] += t[k]*x[j];
1074: }
1075: for (j=0;j<=nv;j++) r[(d-1)*lr+j] = t[d-1]*x[j];
1076: return(0);
1077: }
1079: /*
1080: Compute continuation vector coefficients for the Rational-Krylov run.
1081: dim(work) >= (end-ini)*(end-ini+1) + end+1 + 2*(end-ini+1), dim(t) = end.
1082: */
1083: 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)
1084: {
1086: PetscScalar *x,*W,*tau,sone=1.0,szero=0.0;
1087: PetscInt i,j,n1,n,nwu=0;
1088: PetscBLASInt info,n_,n1_,one=1,dim,lds_;
1089: NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
1092: if (!ctx->nshifts || !end) {
1093: t[0] = 1;
1094: PetscArraycpy(cont,S+end*lds,lds);
1095: } else {
1096: n = end-ini;
1097: n1 = n+1;
1098: x = work+nwu;
1099: nwu += end+1;
1100: tau = work+nwu;
1101: nwu += n;
1102: W = work+nwu;
1103: nwu += n1*n;
1104: for (j=ini;j<end;j++) {
1105: for (i=ini;i<=end;i++) W[(j-ini)*n1+i-ini] = K[j*ld+i] -H[j*ld+i]*sigma;
1106: }
1107: PetscBLASIntCast(n,&n_);
1108: PetscBLASIntCast(n1,&n1_);
1109: PetscBLASIntCast(end+1,&dim);
1110: PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&n1_,&n_,W,&n1_,tau,work+nwu,&n1_,&info));
1111: SlepcCheckLapackInfo("geqrf",info);
1112: for (i=0;i<end;i++) t[i] = 0.0;
1113: t[end] = 1.0;
1114: for (j=n-1;j>=0;j--) {
1115: for (i=0;i<ini+j;i++) x[i] = 0.0;
1116: x[ini+j] = 1.0;
1117: for (i=j+1;i<n1;i++) x[i+ini] = W[i+n1*j];
1118: tau[j] = PetscConj(tau[j]);
1119: PetscStackCallBLAS("LAPACKlarf",LAPACKlarf_("L",&dim,&one,x,&one,tau+j,t,&dim,work+nwu));
1120: }
1121: PetscBLASIntCast(lds,&lds_);
1122: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&lds_,&n1_,&sone,S,&lds_,t,&one,&szero,cont,&one));
1123: }
1124: return(0);
1125: }
1127: /*
1128: Compute a run of Arnoldi iterations
1129: */
1130: PetscErrorCode NEPNLEIGSTOARrun(NEP nep,PetscScalar *K,PetscScalar *H,PetscInt ldh,BV W,PetscInt k,PetscInt *M,PetscBool *breakdown,Vec *t_)
1131: {
1133: NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
1134: PetscInt i,j,m=*M,lwa,deg=ctx->nmat-1,lds,nqt,ld,l;
1135: Vec t;
1136: PetscReal norm;
1137: PetscScalar *x,*work,*tt,sigma,*cont,*S;
1138: PetscBool lindep;
1139: Mat MS;
1142: BVTensorGetFactors(ctx->V,NULL,&MS);
1143: MatDenseGetArray(MS,&S);
1144: BVGetSizes(nep->V,NULL,NULL,&ld);
1145: lds = ld*deg;
1146: BVGetActiveColumns(nep->V,&l,&nqt);
1147: lwa = PetscMax(ld,deg)+(m+1)*(m+1)+4*(m+1);
1148: PetscMalloc4(ld,&x,lwa,&work,m+1,&tt,lds,&cont);
1149: BVSetActiveColumns(ctx->V,0,m);
1150: for (j=k;j<m;j++) {
1151: sigma = ctx->shifts[(++(ctx->idxrk))%ctx->nshiftsw];
1153: /* Continuation vector */
1154: NEPNLEIGS_RKcontinuation(nep,0,j,K,H,ldh,sigma,S,lds,cont,tt,work);
1156: /* apply operator */
1157: BVGetColumn(nep->V,nqt,&t);
1158: NEPTOARExtendBasis(nep,(ctx->idxrk)%ctx->nshiftsw,cont,ld,nqt,W,nep->V,t,S+(j+1)*lds,ld,t_);
1159: BVRestoreColumn(nep->V,nqt,&t);
1161: /* orthogonalize */
1162: BVOrthogonalizeColumn(nep->V,nqt,x,&norm,&lindep);
1163: if (!lindep) {
1164: x[nqt] = norm;
1165: BVScaleColumn(nep->V,nqt,1.0/norm);
1166: nqt++;
1167: } else x[nqt] = 0.0;
1169: NEPTOARCoefficients(nep,sigma,nqt-1,cont,ld,S+(j+1)*lds,ld,x,work);
1171: /* Level-2 orthogonalization */
1172: BVOrthogonalizeColumn(ctx->V,j+1,H+j*ldh,&norm,breakdown);
1173: H[j+1+ldh*j] = norm;
1174: if (ctx->nshifts) {
1175: for (i=0;i<=j;i++) K[i+ldh*j] = sigma*H[i+ldh*j] + tt[i];
1176: K[j+1+ldh*j] = sigma*H[j+1+ldh*j];
1177: }
1178: if (*breakdown) {
1179: *M = j+1;
1180: break;
1181: }
1182: BVScaleColumn(ctx->V,j+1,1.0/norm);
1183: BVSetActiveColumns(nep->V,l,nqt);
1184: }
1185: PetscFree4(x,work,tt,cont);
1186: MatDenseRestoreArray(MS,&S);
1187: BVTensorRestoreFactors(ctx->V,NULL,&MS);
1188: return(0);
1189: }
1191: PetscErrorCode NEPSolve_NLEIGS(NEP nep)
1192: {
1194: NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
1195: PetscInt i,k=0,l,nv=0,ld,lds,ldds,nq,newn;
1196: PetscInt deg=ctx->nmat-1,nconv=0;
1197: PetscScalar *S,*Q,*H,*pU,*K,betak=0,*eigr,*eigi;
1198: PetscReal betah;
1199: PetscBool falselock=PETSC_FALSE,breakdown=PETSC_FALSE;
1200: BV W;
1201: Mat MS,MQ,U;
1204: if (ctx->lock) {
1205: /* undocumented option to use a cheaper locking instead of the true locking */
1206: PetscOptionsGetBool(NULL,NULL,"-nep_nleigs_falselocking",&falselock,NULL);
1207: }
1209: BVGetSizes(nep->V,NULL,NULL,&ld);
1210: lds = deg*ld;
1211: DSGetLeadingDimension(nep->ds,&ldds);
1212: if (!ctx->nshifts) {
1213: PetscMalloc2(nep->ncv,&eigr,nep->ncv,&eigi);
1214: } else { eigr = nep->eigr; eigi = nep->eigi; }
1215: BVDuplicateResize(nep->V,PetscMax(nep->nt-1,ctx->nmat-1),&W);
1218: /* clean projected matrix (including the extra-arrow) */
1219: DSGetArray(nep->ds,DS_MAT_A,&H);
1220: PetscArrayzero(H,ldds*ldds);
1221: DSRestoreArray(nep->ds,DS_MAT_A,&H);
1222: if (ctx->nshifts) {
1223: DSGetArray(nep->ds,DS_MAT_B,&H);
1224: PetscArrayzero(H,ldds*ldds);
1225: DSRestoreArray(nep->ds,DS_MAT_B,&H);
1226: }
1228: /* Get the starting Arnoldi vector */
1229: BVTensorBuildFirstColumn(ctx->V,nep->nini);
1231: /* Restart loop */
1232: l = 0;
1233: while (nep->reason == NEP_CONVERGED_ITERATING) {
1234: nep->its++;
1236: /* Compute an nv-step Krylov relation */
1237: nv = PetscMin(nep->nconv+nep->mpd,nep->ncv);
1238: if (ctx->nshifts) { DSGetArray(nep->ds,DS_MAT_A,&K); }
1239: DSGetArray(nep->ds,ctx->nshifts?DS_MAT_B:DS_MAT_A,&H);
1240: NEPNLEIGSTOARrun(nep,K,H,ldds,W,nep->nconv+l,&nv,&breakdown,nep->work);
1241: betah = PetscAbsScalar(H[(nv-1)*ldds+nv]);
1242: DSRestoreArray(nep->ds,ctx->nshifts?DS_MAT_B:DS_MAT_A,&H);
1243: if (ctx->nshifts) {
1244: betak = K[(nv-1)*ldds+nv];
1245: DSRestoreArray(nep->ds,DS_MAT_A,&K);
1246: }
1247: DSSetDimensions(nep->ds,nv,0,nep->nconv,nep->nconv+l);
1248: if (l==0) {
1249: DSSetState(nep->ds,DS_STATE_INTERMEDIATE);
1250: } else {
1251: DSSetState(nep->ds,DS_STATE_RAW);
1252: }
1254: /* Solve projected problem */
1255: DSSolve(nep->ds,nep->eigr,nep->eigi);
1256: DSSort(nep->ds,nep->eigr,nep->eigi,NULL,NULL,NULL);
1257: if (!ctx->nshifts) {
1258: DSUpdateExtraRow(nep->ds);
1259: }
1260: DSSynchronize(nep->ds,nep->eigr,nep->eigi);
1262: /* Check convergence */
1263: NEPNLEIGSKrylovConvergence(nep,PETSC_FALSE,nep->nconv,nv-nep->nconv,betah,betak,&k,nep->work);
1264: (*nep->stopping)(nep,nep->its,nep->max_it,k,nep->nev,&nep->reason,nep->stoppingctx);
1266: /* Update l */
1267: if (nep->reason != NEP_CONVERGED_ITERATING || breakdown) l = 0;
1268: else {
1269: l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
1270: if (!breakdown) {
1271: /* Prepare the Rayleigh quotient for restart */
1272: if (!ctx->nshifts) {
1273: DSTruncate(nep->ds,k+l);
1274: DSGetDimensions(nep->ds,&newn,NULL,NULL,NULL,NULL);
1275: l = newn-k;
1276: } else {
1277: DSGetArray(nep->ds,DS_MAT_Q,&Q);
1278: DSGetArray(nep->ds,DS_MAT_B,&H);
1279: DSGetArray(nep->ds,DS_MAT_A,&K);
1280: for (i=ctx->lock?k:0;i<k+l;i++) {
1281: H[k+l+i*ldds] = betah*Q[nv-1+i*ldds];
1282: K[k+l+i*ldds] = betak*Q[nv-1+i*ldds];
1283: }
1284: DSRestoreArray(nep->ds,DS_MAT_B,&H);
1285: DSRestoreArray(nep->ds,DS_MAT_A,&K);
1286: DSRestoreArray(nep->ds,DS_MAT_Q,&Q);
1287: }
1288: }
1289: }
1290: nconv = k;
1291: if (!ctx->lock && nep->reason == NEP_CONVERGED_ITERATING && !breakdown) { l += k; k = 0; }
1293: /* Update S */
1294: DSGetMat(nep->ds,ctx->nshifts?DS_MAT_Z:DS_MAT_Q,&MQ);
1295: BVMultInPlace(ctx->V,MQ,nep->nconv,k+l);
1296: MatDestroy(&MQ);
1298: /* Copy last column of S */
1299: BVCopyColumn(ctx->V,nv,k+l);
1301: if (breakdown && nep->reason == NEP_CONVERGED_ITERATING) {
1302: /* Stop if breakdown */
1303: PetscInfo2(nep,"Breakdown (it=%D norm=%g)\n",nep->its,(double)betah);
1304: nep->reason = NEP_DIVERGED_BREAKDOWN;
1305: }
1306: if (nep->reason != NEP_CONVERGED_ITERATING) l--;
1307: /* truncate S */
1308: BVGetActiveColumns(nep->V,NULL,&nq);
1309: if (k+l+deg<=nq) {
1310: BVSetActiveColumns(ctx->V,nep->nconv,k+l+1);
1311: if (!falselock && ctx->lock) {
1312: BVTensorCompress(ctx->V,k-nep->nconv);
1313: } else {
1314: BVTensorCompress(ctx->V,0);
1315: }
1316: }
1317: nep->nconv = k;
1318: if (!ctx->nshifts) {
1319: for (i=0;i<nv;i++) { eigr[i] = nep->eigr[i]; eigi[i] = nep->eigi[i]; }
1320: NEPNLEIGSBackTransform((PetscObject)nep,nv,eigr,eigi);
1321: }
1322: NEPMonitor(nep,nep->its,nconv,eigr,eigi,nep->errest,nv);
1323: }
1324: nep->nconv = nconv;
1325: if (nep->nconv>0) {
1326: BVSetActiveColumns(ctx->V,0,nep->nconv);
1327: BVGetActiveColumns(nep->V,NULL,&nq);
1328: BVSetActiveColumns(nep->V,0,nq);
1329: if (nq>nep->nconv) {
1330: BVTensorCompress(ctx->V,nep->nconv);
1331: BVSetActiveColumns(nep->V,0,nep->nconv);
1332: nq = nep->nconv;
1333: }
1334: if (ctx->nshifts) {
1335: DSGetMat(nep->ds,DS_MAT_B,&MQ);
1336: BVMultInPlace(ctx->V,MQ,0,nep->nconv);
1337: MatDestroy(&MQ);
1338: }
1339: BVTensorGetFactors(ctx->V,NULL,&MS);
1340: MatDenseGetArray(MS,&S);
1341: PetscMalloc1(nq*nep->nconv,&pU);
1342: for (i=0;i<nep->nconv;i++) {
1343: PetscArraycpy(pU+i*nq,S+i*lds,nq);
1344: }
1345: MatDenseRestoreArray(MS,&S);
1346: BVTensorRestoreFactors(ctx->V,NULL,&MS);
1347: MatCreateSeqDense(PETSC_COMM_SELF,nq,nep->nconv,pU,&U);
1348: BVSetActiveColumns(nep->V,0,nq);
1349: BVMultInPlace(nep->V,U,0,nep->nconv);
1350: BVSetActiveColumns(nep->V,0,nep->nconv);
1351: MatDestroy(&U);
1352: PetscFree(pU);
1353: /* truncate Schur decomposition and change the state to raw so that
1354: DSVectors() computes eigenvectors from scratch */
1355: DSSetDimensions(nep->ds,nep->nconv,0,0,0);
1356: DSSetState(nep->ds,DS_STATE_RAW);
1357: }
1359: /* Map eigenvalues back to the original problem */
1360: if (!ctx->nshifts) {
1361: NEPNLEIGSBackTransform((PetscObject)nep,nep->nconv,nep->eigr,nep->eigi);
1362: PetscFree2(eigr,eigi);
1363: }
1364: BVDestroy(&W);
1365: return(0);
1366: }
1368: static PetscErrorCode NEPNLEIGSSetSingularitiesFunction_NLEIGS(NEP nep,PetscErrorCode (*fun)(NEP,PetscInt*,PetscScalar*,void*),void *ctx)
1369: {
1370: NEP_NLEIGS *nepctx=(NEP_NLEIGS*)nep->data;
1373: if (fun) nepctx->computesingularities = fun;
1374: if (ctx) nepctx->singularitiesctx = ctx;
1375: nep->state = NEP_STATE_INITIAL;
1376: return(0);
1377: }
1379: /*@C
1380: NEPNLEIGSSetSingularitiesFunction - Sets a user function to compute a discretization
1381: of the singularity set (where T(.) is not analytic).
1383: Logically Collective on nep
1385: Input Parameters:
1386: + nep - the NEP context
1387: . fun - user function (if NULL then NEP retains any previously set value)
1388: - ctx - [optional] user-defined context for private data for the function
1389: (may be NULL, in which case NEP retains any previously set value)
1391: Calling Sequence of fun:
1392: $ fun(NEP nep,PetscInt *maxnp,PetscScalar *xi,void *ctx)
1394: + nep - the NEP context
1395: . maxnp - on input number of requested points in the discretization (can be set)
1396: . xi - computed values of the discretization
1397: - ctx - optional context, as set by NEPNLEIGSSetSingularitiesFunction()
1399: Notes:
1400: The user-defined function can set a smaller value of maxnp if necessary.
1401: It is wrong to return a larger value.
1403: If the problem type has been set to rational with NEPSetProblemType(),
1404: then it is not necessary to set the singularities explicitly since the
1405: solver will try to determine them automatically.
1407: Level: intermediate
1409: .seealso: NEPNLEIGSGetSingularitiesFunction(), NEPSetProblemType()
1410: @*/
1411: PetscErrorCode NEPNLEIGSSetSingularitiesFunction(NEP nep,PetscErrorCode (*fun)(NEP,PetscInt*,PetscScalar*,void*),void *ctx)
1412: {
1417: PetscTryMethod(nep,"NEPNLEIGSSetSingularitiesFunction_C",(NEP,PetscErrorCode(*)(NEP,PetscInt*,PetscScalar*,void*),void*),(nep,fun,ctx));
1418: return(0);
1419: }
1421: static PetscErrorCode NEPNLEIGSGetSingularitiesFunction_NLEIGS(NEP nep,PetscErrorCode (**fun)(NEP,PetscInt*,PetscScalar*,void*),void **ctx)
1422: {
1423: NEP_NLEIGS *nepctx=(NEP_NLEIGS*)nep->data;
1426: if (fun) *fun = nepctx->computesingularities;
1427: if (ctx) *ctx = nepctx->singularitiesctx;
1428: return(0);
1429: }
1431: /*@C
1432: NEPNLEIGSGetSingularitiesFunction - Returns the Function and optionally the user
1433: provided context for computing a discretization of the singularity set.
1435: Not Collective
1437: Input Parameter:
1438: . nep - the nonlinear eigensolver context
1440: Output Parameters:
1441: + fun - location to put the function (or NULL)
1442: - ctx - location to stash the function context (or NULL)
1444: Level: advanced
1446: .seealso: NEPNLEIGSSetSingularitiesFunction()
1447: @*/
1448: PetscErrorCode NEPNLEIGSGetSingularitiesFunction(NEP nep,PetscErrorCode (**fun)(NEP,PetscInt*,PetscScalar*,void*),void **ctx)
1449: {
1454: PetscUseMethod(nep,"NEPNLEIGSGetSingularitiesFunction_C",(NEP,PetscErrorCode(**)(NEP,PetscInt*,PetscScalar*,void*),void**),(nep,fun,ctx));
1455: return(0);
1456: }
1458: static PetscErrorCode NEPNLEIGSSetRestart_NLEIGS(NEP nep,PetscReal keep)
1459: {
1460: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1463: if (keep==PETSC_DEFAULT) ctx->keep = 0.5;
1464: else {
1465: 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]");
1466: ctx->keep = keep;
1467: }
1468: return(0);
1469: }
1471: /*@
1472: NEPNLEIGSSetRestart - Sets the restart parameter for the NLEIGS
1473: method, in particular the proportion of basis vectors that must be kept
1474: after restart.
1476: Logically Collective on nep
1478: Input Parameters:
1479: + nep - the nonlinear eigensolver context
1480: - keep - the number of vectors to be kept at restart
1482: Options Database Key:
1483: . -nep_nleigs_restart - Sets the restart parameter
1485: Notes:
1486: Allowed values are in the range [0.1,0.9]. The default is 0.5.
1488: Level: advanced
1490: .seealso: NEPNLEIGSGetRestart()
1491: @*/
1492: PetscErrorCode NEPNLEIGSSetRestart(NEP nep,PetscReal keep)
1493: {
1499: PetscTryMethod(nep,"NEPNLEIGSSetRestart_C",(NEP,PetscReal),(nep,keep));
1500: return(0);
1501: }
1503: static PetscErrorCode NEPNLEIGSGetRestart_NLEIGS(NEP nep,PetscReal *keep)
1504: {
1505: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1508: *keep = ctx->keep;
1509: return(0);
1510: }
1512: /*@
1513: NEPNLEIGSGetRestart - Gets the restart parameter used in the NLEIGS method.
1515: Not Collective
1517: Input Parameter:
1518: . nep - the nonlinear eigensolver context
1520: Output Parameter:
1521: . keep - the restart parameter
1523: Level: advanced
1525: .seealso: NEPNLEIGSSetRestart()
1526: @*/
1527: PetscErrorCode NEPNLEIGSGetRestart(NEP nep,PetscReal *keep)
1528: {
1534: PetscUseMethod(nep,"NEPNLEIGSGetRestart_C",(NEP,PetscReal*),(nep,keep));
1535: return(0);
1536: }
1538: static PetscErrorCode NEPNLEIGSSetLocking_NLEIGS(NEP nep,PetscBool lock)
1539: {
1540: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1543: ctx->lock = lock;
1544: return(0);
1545: }
1547: /*@
1548: NEPNLEIGSSetLocking - Choose between locking and non-locking variants of
1549: the NLEIGS method.
1551: Logically Collective on nep
1553: Input Parameters:
1554: + nep - the nonlinear eigensolver context
1555: - lock - true if the locking variant must be selected
1557: Options Database Key:
1558: . -nep_nleigs_locking - Sets the locking flag
1560: Notes:
1561: The default is to lock converged eigenpairs when the method restarts.
1562: This behaviour can be changed so that all directions are kept in the
1563: working subspace even if already converged to working accuracy (the
1564: non-locking variant).
1566: Level: advanced
1568: .seealso: NEPNLEIGSGetLocking()
1569: @*/
1570: PetscErrorCode NEPNLEIGSSetLocking(NEP nep,PetscBool lock)
1571: {
1577: PetscTryMethod(nep,"NEPNLEIGSSetLocking_C",(NEP,PetscBool),(nep,lock));
1578: return(0);
1579: }
1581: static PetscErrorCode NEPNLEIGSGetLocking_NLEIGS(NEP nep,PetscBool *lock)
1582: {
1583: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1586: *lock = ctx->lock;
1587: return(0);
1588: }
1590: /*@
1591: NEPNLEIGSGetLocking - Gets the locking flag used in the NLEIGS method.
1593: Not Collective
1595: Input Parameter:
1596: . nep - the nonlinear eigensolver context
1598: Output Parameter:
1599: . lock - the locking flag
1601: Level: advanced
1603: .seealso: NEPNLEIGSSetLocking()
1604: @*/
1605: PetscErrorCode NEPNLEIGSGetLocking(NEP nep,PetscBool *lock)
1606: {
1612: PetscUseMethod(nep,"NEPNLEIGSGetLocking_C",(NEP,PetscBool*),(nep,lock));
1613: return(0);
1614: }
1616: static PetscErrorCode NEPNLEIGSSetInterpolation_NLEIGS(NEP nep,PetscReal tol,PetscInt degree)
1617: {
1619: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1622: if (tol == PETSC_DEFAULT) {
1623: ctx->ddtol = PETSC_DEFAULT;
1624: nep->state = NEP_STATE_INITIAL;
1625: } else {
1626: if (tol <= 0.0) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of tol. Must be > 0");
1627: ctx->ddtol = tol;
1628: }
1629: if (degree == PETSC_DEFAULT || degree == PETSC_DECIDE) {
1630: ctx->ddmaxit = 0;
1631: if (nep->state) { NEPReset(nep); }
1632: nep->state = NEP_STATE_INITIAL;
1633: } else {
1634: if (degree <= 0) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of degree. Must be > 0");
1635: if (ctx->ddmaxit != degree) {
1636: ctx->ddmaxit = degree;
1637: if (nep->state) { NEPReset(nep); }
1638: nep->state = NEP_STATE_INITIAL;
1639: }
1640: }
1641: return(0);
1642: }
1644: /*@
1645: NEPNLEIGSSetInterpolation - Sets the tolerance and maximum degree
1646: when building the interpolation via divided differences.
1648: Logically Collective on nep
1650: Input Parameters:
1651: + nep - the nonlinear eigensolver context
1652: . tol - tolerance to stop computing divided differences
1653: - degree - maximum degree of interpolation
1655: Options Database Key:
1656: + -nep_nleigs_interpolation_tol <tol> - Sets the tolerance to stop computing divided differences
1657: - -nep_nleigs_interpolation_degree <degree> - Sets the maximum degree of interpolation
1659: Notes:
1660: Use PETSC_DEFAULT for either argument to assign a reasonably good value.
1662: Level: advanced
1664: .seealso: NEPNLEIGSGetInterpolation()
1665: @*/
1666: PetscErrorCode NEPNLEIGSSetInterpolation(NEP nep,PetscReal tol,PetscInt degree)
1667: {
1674: PetscTryMethod(nep,"NEPNLEIGSSetInterpolation_C",(NEP,PetscReal,PetscInt),(nep,tol,degree));
1675: return(0);
1676: }
1678: static PetscErrorCode NEPNLEIGSGetInterpolation_NLEIGS(NEP nep,PetscReal *tol,PetscInt *degree)
1679: {
1680: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1683: if (tol) *tol = ctx->ddtol;
1684: if (degree) *degree = ctx->ddmaxit;
1685: return(0);
1686: }
1688: /*@
1689: NEPNLEIGSGetInterpolation - Gets the tolerance and maximum degree
1690: when building the interpolation via divided differences.
1692: Not Collective
1694: Input Parameter:
1695: . nep - the nonlinear eigensolver context
1697: Output Parameter:
1698: + tol - tolerance to stop computing divided differences
1699: - degree - maximum degree of interpolation
1701: Level: advanced
1703: .seealso: NEPNLEIGSSetInterpolation()
1704: @*/
1705: PetscErrorCode NEPNLEIGSGetInterpolation(NEP nep,PetscReal *tol,PetscInt *degree)
1706: {
1711: PetscTryMethod(nep,"NEPNLEIGSGetInterpolation_C",(NEP,PetscReal*,PetscInt*),(nep,tol,degree));
1712: return(0);
1713: }
1715: static PetscErrorCode NEPNLEIGSSetRKShifts_NLEIGS(NEP nep,PetscInt ns,PetscScalar *shifts)
1716: {
1718: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1719: PetscInt i;
1722: if (ns<0) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"Number of shifts must be non-negative");
1723: if (ctx->nshifts) { PetscFree(ctx->shifts); }
1724: for (i=0;i<ctx->nshiftsw;i++) { KSPDestroy(&ctx->ksp[i]); }
1725: PetscFree(ctx->ksp);
1726: ctx->ksp = NULL;
1727: if (ns) {
1728: PetscMalloc1(ns,&ctx->shifts);
1729: for (i=0;i<ns;i++) ctx->shifts[i] = shifts[i];
1730: }
1731: ctx->nshifts = ns;
1732: nep->state = NEP_STATE_INITIAL;
1733: return(0);
1734: }
1736: /*@
1737: NEPNLEIGSSetRKShifts - Sets a list of shifts to be used in the Rational
1738: Krylov method.
1740: Logically Collective on nep
1742: Input Parameters:
1743: + nep - the nonlinear eigensolver context
1744: . ns - number of shifts
1745: - shifts - array of scalar values specifying the shifts
1747: Options Database Key:
1748: . -nep_nleigs_rk_shifts - Sets the list of shifts
1750: Notes:
1751: If only one shift is provided, the built subspace built is equivalent to
1752: shift-and-invert Krylov-Schur (provided that the absolute convergence
1753: criterion is used).
1755: In the case of real scalars, complex shifts are not allowed. In the
1756: command line, a comma-separated list of complex values can be provided with
1757: the format [+/-][realnumber][+/-]realnumberi with no spaces, e.g.
1758: -nep_nleigs_rk_shifts 1.0+2.0i,1.5+2.0i,1.0+1.5i
1760: Use ns=0 to remove previously set shifts.
1762: Level: advanced
1764: .seealso: NEPNLEIGSGetRKShifts()
1765: @*/
1766: PetscErrorCode NEPNLEIGSSetRKShifts(NEP nep,PetscInt ns,PetscScalar shifts[])
1767: {
1774: PetscTryMethod(nep,"NEPNLEIGSSetRKShifts_C",(NEP,PetscInt,PetscScalar*),(nep,ns,shifts));
1775: return(0);
1776: }
1778: static PetscErrorCode NEPNLEIGSGetRKShifts_NLEIGS(NEP nep,PetscInt *ns,PetscScalar **shifts)
1779: {
1781: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1782: PetscInt i;
1785: *ns = ctx->nshifts;
1786: if (ctx->nshifts) {
1787: PetscMalloc1(ctx->nshifts,shifts);
1788: for (i=0;i<ctx->nshifts;i++) (*shifts)[i] = ctx->shifts[i];
1789: }
1790: return(0);
1791: }
1793: /*@C
1794: NEPNLEIGSGetRKShifts - Gets the list of shifts used in the Rational
1795: Krylov method.
1797: Not Collective
1799: Input Parameter:
1800: . nep - the nonlinear eigensolver context
1802: Output Parameter:
1803: + ns - number of shifts
1804: - shifts - array of shifts
1806: Note:
1807: The user is responsible for deallocating the returned array.
1809: Level: advanced
1811: .seealso: NEPNLEIGSSetRKShifts()
1812: @*/
1813: PetscErrorCode NEPNLEIGSGetRKShifts(NEP nep,PetscInt *ns,PetscScalar *shifts[])
1814: {
1821: PetscTryMethod(nep,"NEPNLEIGSGetRKShifts_C",(NEP,PetscInt*,PetscScalar**),(nep,ns,shifts));
1822: return(0);
1823: }
1825: static PetscErrorCode NEPNLEIGSGetKSPs_NLEIGS(NEP nep,PetscInt *nsolve,KSP **ksp)
1826: {
1828: NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
1829: PetscInt i;
1830: PC pc;
1833: if (!ctx->ksp) {
1834: NEPNLEIGSSetShifts(nep,&ctx->nshiftsw);
1835: PetscMalloc1(ctx->nshiftsw,&ctx->ksp);
1836: for (i=0;i<ctx->nshiftsw;i++) {
1837: KSPCreate(PetscObjectComm((PetscObject)nep),&ctx->ksp[i]);
1838: PetscObjectIncrementTabLevel((PetscObject)ctx->ksp[i],(PetscObject)nep,1);
1839: KSPSetOptionsPrefix(ctx->ksp[i],((PetscObject)nep)->prefix);
1840: KSPAppendOptionsPrefix(ctx->ksp[i],"nep_nleigs_");
1841: PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->ksp[i]);
1842: PetscObjectSetOptions((PetscObject)ctx->ksp[i],((PetscObject)nep)->options);
1843: KSPSetErrorIfNotConverged(ctx->ksp[i],PETSC_TRUE);
1844: KSPSetTolerances(ctx->ksp[i],SLEPC_DEFAULT_TOL,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
1845: KSPGetPC(ctx->ksp[i],&pc);
1846: KSPSetType(ctx->ksp[i],KSPPREONLY);
1847: PCSetType(pc,PCLU);
1848: }
1849: }
1850: if (nsolve) *nsolve = ctx->nshiftsw;
1851: if (ksp) *ksp = ctx->ksp;
1852: return(0);
1853: }
1855: /*@C
1856: NEPNLEIGSGetKSPs - Retrieve the array of linear solver objects associated with
1857: the nonlinear eigenvalue solver.
1859: Not Collective
1861: Input Parameter:
1862: . nep - nonlinear eigenvalue solver
1864: Output Parameters:
1865: + nsolve - number of returned KSP objects
1866: - ksp - array of linear solver object
1868: Notes:
1869: The number of KSP objects is equal to the number of shifts provided by the user,
1870: or 1 if the user did not provide shifts.
1872: Level: advanced
1874: .seealso: NEPNLEIGSSetRKShifts()
1875: @*/
1876: PetscErrorCode NEPNLEIGSGetKSPs(NEP nep,PetscInt *nsolve,KSP **ksp)
1877: {
1882: PetscUseMethod(nep,"NEPNLEIGSGetKSPs_C",(NEP,PetscInt*,KSP**),(nep,nsolve,ksp));
1883: return(0);
1884: }
1886: static PetscErrorCode NEPNLEIGSSetFullBasis_NLEIGS(NEP nep,PetscBool fullbasis)
1887: {
1888: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1891: if (fullbasis!=ctx->fullbasis) {
1892: ctx->fullbasis = fullbasis;
1893: nep->state = NEP_STATE_INITIAL;
1894: nep->useds = PetscNot(fullbasis);
1895: }
1896: return(0);
1897: }
1899: /*@
1900: NEPNLEIGSSetFullBasis - Choose between TOAR-basis (default) and full-basis
1901: variants of the NLEIGS method.
1903: Logically Collective on nep
1905: Input Parameters:
1906: + nep - the nonlinear eigensolver context
1907: - fullbasis - true if the full-basis variant must be selected
1909: Options Database Key:
1910: . -nep_nleigs_full_basis - Sets the full-basis flag
1912: Notes:
1913: The default is to use a compact representation of the Krylov basis, that is,
1914: V = (I otimes U) S, with a tensor BV. This behaviour can be changed so that
1915: the full basis V is explicitly stored and operated with. This variant is more
1916: expensive in terms of memory and computation, but is necessary in some cases,
1917: particularly for two-sided computations, see NEPSetTwoSided().
1919: In the full-basis variant, the NLEIGS solver uses an EPS object to explicitly
1920: solve the linearized eigenproblem, see NEPNLEIGSGetEPS().
1922: Level: advanced
1924: .seealso: NEPNLEIGSGetFullBasis(), NEPNLEIGSGetEPS(), NEPSetTwoSided(), BVCreateTensor()
1925: @*/
1926: PetscErrorCode NEPNLEIGSSetFullBasis(NEP nep,PetscBool fullbasis)
1927: {
1933: PetscTryMethod(nep,"NEPNLEIGSSetFullBasis_C",(NEP,PetscBool),(nep,fullbasis));
1934: return(0);
1935: }
1937: static PetscErrorCode NEPNLEIGSGetFullBasis_NLEIGS(NEP nep,PetscBool *fullbasis)
1938: {
1939: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1942: *fullbasis = ctx->fullbasis;
1943: return(0);
1944: }
1946: /*@
1947: NEPNLEIGSGetFullBasis - Gets the flag that indicates if NLEIGS is using the
1948: full-basis variant.
1950: Not Collective
1952: Input Parameter:
1953: . nep - the nonlinear eigensolver context
1955: Output Parameter:
1956: . fullbasis - the flag
1958: Level: advanced
1960: .seealso: NEPNLEIGSSetFullBasis()
1961: @*/
1962: PetscErrorCode NEPNLEIGSGetFullBasis(NEP nep,PetscBool *fullbasis)
1963: {
1969: PetscUseMethod(nep,"NEPNLEIGSGetFullBasis_C",(NEP,PetscBool*),(nep,fullbasis));
1970: return(0);
1971: }
1973: #define SHIFTMAX 30
1975: PetscErrorCode NEPSetFromOptions_NLEIGS(PetscOptionItems *PetscOptionsObject,NEP nep)
1976: {
1978: NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
1979: PetscInt i=0,k;
1980: PetscBool flg1,flg2,b;
1981: PetscReal r;
1982: PetscScalar array[SHIFTMAX];
1985: PetscOptionsHead(PetscOptionsObject,"NEP NLEIGS Options");
1987: PetscOptionsReal("-nep_nleigs_restart","Proportion of vectors kept after restart","NEPNLEIGSSetRestart",0.5,&r,&flg1);
1988: if (flg1) { NEPNLEIGSSetRestart(nep,r); }
1990: PetscOptionsBool("-nep_nleigs_locking","Choose between locking and non-locking variants","NEPNLEIGSSetLocking",PETSC_FALSE,&b,&flg1);
1991: if (flg1) { NEPNLEIGSSetLocking(nep,b); }
1993: PetscOptionsBool("-nep_nleigs_full_basis","Choose between TOAR and full-basis variants","NEPNLEIGSSetFullBasis",PETSC_FALSE,&b,&flg1);
1994: if (flg1) { NEPNLEIGSSetFullBasis(nep,b); }
1996: NEPNLEIGSGetInterpolation(nep,&r,&i);
1997: if (!i) i = PETSC_DEFAULT;
1998: PetscOptionsInt("-nep_nleigs_interpolation_degree","Maximum number of terms for interpolation via divided differences","NEPNLEIGSSetInterpolation",i,&i,&flg1);
1999: PetscOptionsReal("-nep_nleigs_interpolation_tol","Tolerance for interpolation via divided differences","NEPNLEIGSSetInterpolation",r,&r,&flg2);
2000: if (flg1 || flg2) { NEPNLEIGSSetInterpolation(nep,r,i); }
2002: k = SHIFTMAX;
2003: for (i=0;i<k;i++) array[i] = 0;
2004: PetscOptionsScalarArray("-nep_nleigs_rk_shifts","Shifts for Rational Krylov","NEPNLEIGSSetRKShifts",array,&k,&flg1);
2005: if (flg1) { NEPNLEIGSSetRKShifts(nep,k,array); }
2007: PetscOptionsTail();
2009: if (!ctx->ksp) { NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp); }
2010: for (i=0;i<ctx->nshiftsw;i++) {
2011: KSPSetFromOptions(ctx->ksp[i]);
2012: }
2014: if (ctx->fullbasis) {
2015: if (!ctx->eps) { NEPNLEIGSGetEPS(nep,&ctx->eps); }
2016: EPSSetFromOptions(ctx->eps);
2017: }
2018: return(0);
2019: }
2021: PetscErrorCode NEPView_NLEIGS(NEP nep,PetscViewer viewer)
2022: {
2024: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
2025: PetscBool isascii;
2026: PetscInt i;
2027: char str[50];
2030: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
2031: if (isascii) {
2032: PetscViewerASCIIPrintf(viewer," %d%% of basis vectors kept after restart\n",(int)(100*ctx->keep));
2033: if (ctx->fullbasis) {
2034: PetscViewerASCIIPrintf(viewer," using the full-basis variant\n");
2035: } else {
2036: PetscViewerASCIIPrintf(viewer," using the %slocking variant\n",ctx->lock?"":"non-");
2037: }
2038: PetscViewerASCIIPrintf(viewer," divided difference terms: used=%D, max=%D\n",ctx->nmat,ctx->ddmaxit);
2039: PetscViewerASCIIPrintf(viewer," tolerance for divided difference convergence: %g\n",(double)ctx->ddtol);
2040: if (ctx->nshifts) {
2041: PetscViewerASCIIPrintf(viewer," RK shifts: ");
2042: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
2043: for (i=0;i<ctx->nshifts;i++) {
2044: SlepcSNPrintfScalar(str,50,ctx->shifts[i],PETSC_FALSE);
2045: PetscViewerASCIIPrintf(viewer,"%s%s",str,(i<ctx->nshifts-1)?",":"");
2046: }
2047: PetscViewerASCIIPrintf(viewer,"\n");
2048: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
2049: }
2050: if (!ctx->ksp) { NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp); }
2051: PetscViewerASCIIPushTab(viewer);
2052: KSPView(ctx->ksp[0],viewer);
2053: PetscViewerASCIIPopTab(viewer);
2054: if (ctx->fullbasis) {
2055: if (!ctx->eps) { NEPNLEIGSGetEPS(nep,&ctx->eps); }
2056: PetscViewerASCIIPushTab(viewer);
2057: EPSView(ctx->eps,viewer);
2058: PetscViewerASCIIPopTab(viewer);
2059: }
2060: }
2061: return(0);
2062: }
2064: PetscErrorCode NEPReset_NLEIGS(NEP nep)
2065: {
2067: PetscInt k;
2068: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
2071: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
2072: PetscFree(ctx->coeffD);
2073: } else {
2074: for (k=0;k<ctx->nmat;k++) { MatDestroy(&ctx->D[k]); }
2075: }
2076: PetscFree4(ctx->s,ctx->xi,ctx->beta,ctx->D);
2077: for (k=0;k<ctx->nshiftsw;k++) { KSPReset(ctx->ksp[k]); }
2078: if (ctx->vrn) {
2079: VecDestroy(&ctx->vrn);
2080: }
2081: if (ctx->fullbasis) {
2082: MatDestroy(&ctx->A);
2083: EPSReset(ctx->eps);
2084: for (k=0;k<4;k++) { VecDestroy(&ctx->w[k]); }
2085: }
2086: return(0);
2087: }
2089: PetscErrorCode NEPDestroy_NLEIGS(NEP nep)
2090: {
2092: PetscInt k;
2093: NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
2096: BVDestroy(&ctx->V);
2097: for (k=0;k<ctx->nshiftsw;k++) { KSPDestroy(&ctx->ksp[k]); }
2098: PetscFree(ctx->ksp);
2099: if (ctx->nshifts) { PetscFree(ctx->shifts); }
2100: if (ctx->fullbasis) { EPSDestroy(&ctx->eps); }
2101: PetscFree(nep->data);
2102: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetSingularitiesFunction_C",NULL);
2103: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetSingularitiesFunction_C",NULL);
2104: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRestart_C",NULL);
2105: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRestart_C",NULL);
2106: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetLocking_C",NULL);
2107: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetLocking_C",NULL);
2108: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetInterpolation_C",NULL);
2109: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetInterpolation_C",NULL);
2110: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRKShifts_C",NULL);
2111: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRKShifts_C",NULL);
2112: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetKSPs_C",NULL);
2113: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetFullBasis_C",NULL);
2114: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetFullBasis_C",NULL);
2115: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetEPS_C",NULL);
2116: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetEPS_C",NULL);
2117: return(0);
2118: }
2120: SLEPC_EXTERN PetscErrorCode NEPCreate_NLEIGS(NEP nep)
2121: {
2123: NEP_NLEIGS *ctx;
2126: PetscNewLog(nep,&ctx);
2127: nep->data = (void*)ctx;
2128: ctx->lock = PETSC_TRUE;
2129: ctx->ddtol = PETSC_DEFAULT;
2131: nep->useds = PETSC_TRUE;
2132: nep->hasts = PETSC_TRUE;
2134: nep->ops->setup = NEPSetUp_NLEIGS;
2135: nep->ops->setfromoptions = NEPSetFromOptions_NLEIGS;
2136: nep->ops->view = NEPView_NLEIGS;
2137: nep->ops->destroy = NEPDestroy_NLEIGS;
2138: nep->ops->reset = NEPReset_NLEIGS;
2140: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetSingularitiesFunction_C",NEPNLEIGSSetSingularitiesFunction_NLEIGS);
2141: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetSingularitiesFunction_C",NEPNLEIGSGetSingularitiesFunction_NLEIGS);
2142: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRestart_C",NEPNLEIGSSetRestart_NLEIGS);
2143: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRestart_C",NEPNLEIGSGetRestart_NLEIGS);
2144: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetLocking_C",NEPNLEIGSSetLocking_NLEIGS);
2145: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetLocking_C",NEPNLEIGSGetLocking_NLEIGS);
2146: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetInterpolation_C",NEPNLEIGSSetInterpolation_NLEIGS);
2147: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetInterpolation_C",NEPNLEIGSGetInterpolation_NLEIGS);
2148: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRKShifts_C",NEPNLEIGSSetRKShifts_NLEIGS);
2149: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRKShifts_C",NEPNLEIGSGetRKShifts_NLEIGS);
2150: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetKSPs_C",NEPNLEIGSGetKSPs_NLEIGS);
2151: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetFullBasis_C",NEPNLEIGSSetFullBasis_NLEIGS);
2152: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetFullBasis_C",NEPNLEIGSGetFullBasis_NLEIGS);
2153: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetEPS_C",NEPNLEIGSSetEPS_NLEIGS);
2154: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetEPS_C",NEPNLEIGSGetEPS_NLEIGS);
2155: return(0);
2156: }