Actual source code: nleigs.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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> /*I "slepcnep.h" I*/
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,*work;
62: PetscBLASInt n_,info,lwork;
63: PetscInt i;
64: #if defined(PETSC_USE_COMPLEX)
65: PetscReal *rwork;
66: #endif
67: #if defined(PETSC_HAVE_ESSL)
68: PetscScalar sdummy;
69: PetscBLASInt idummy,io=0;
70: PetscScalar *wri;
71: #endif
74: #if defined(PETSC_MISSING_LAPACK_GEEV)
75: *avail = PETSC_FALSE;
76: #else
77: *avail = PETSC_TRUE;
78: if (deg>0) {
79: PetscCalloc1(deg*deg,&C);
80: PetscBLASIntCast(deg,&n_);
81: for (i=0;i<deg-1;i++) {
82: C[(deg+1)*i+1] = 1.0;
83: C[(deg-1)*deg+i] = -polcoeffs[deg-i]/polcoeffs[0];
84: }
85: C[deg*deg+-1] = -polcoeffs[1]/polcoeffs[0];
86: PetscBLASIntCast(3*deg,&lwork);
88: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
89: #if !defined(PETSC_HAVE_ESSL)
90: #if !defined(PETSC_USE_COMPLEX)
91: PetscMalloc1(lwork,&work);
92: PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_("N","N",&n_,C,&n_,wr,wi,NULL,&n_,NULL,&n_,work,&lwork,&info));
93: if (info) *avail = PETSC_FALSE;
94: PetscFree(work);
95: #else
96: PetscMalloc2(2*deg,&rwork,lwork,&work);
97: PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_("N","N",&n_,C,&n_,wr,NULL,&n_,NULL,&n_,work,&lwork,rwork,&info));
98: if (info) *avail = PETSC_FALSE;
99: PetscFree2(rwork,work);
100: #endif
101: #else
102: PetscMalloc2(lwork,&work,2*deg,&wri);
103: PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_(&io,C,&n_,wri,&sdummy,&idummy,&idummy,&n_,work,&lwork));
104: #if !defined(PETSC_USE_COMPLEX)
105: for (i=0;i<deg;i++) {
106: wr[i] = wri[2*i];
107: wi[i] = wri[2*i+1];
108: }
109: #else
110: for (i=0;i<deg;i++) wr[i] = wri[i];
111: #endif
112: PetscFree2(work,wri);
113: #endif
114: PetscFPTrapPop();
115: PetscFree(C);
116: }
117: #endif
118: return(0);
119: }
121: static PetscErrorCode NEPNLEIGSAuxiliarRmDuplicates(PetscInt nin,PetscScalar *pin,PetscInt *nout,PetscScalar *pout,PetscInt max)
122: {
123: PetscInt i,j;
126: for (i=0;i<nin;i++) {
127: if (max && *nout>max) break;
128: pout[(*nout)++] = pin[i];
129: for (j=0;j<*nout-1;j++)
130: if (PetscAbsScalar(pin[i]-pout[j])<PETSC_MACHINE_EPSILON*100) {
131: (*nout)--;
132: break;
133: }
134: }
135: return(0);
136: }
138: static PetscErrorCode NEPNLEIGSFNSingularities(FN f,PetscInt *nisol,PetscScalar **isol,PetscBool *rational)
139: {
141: FNCombineType ctype;
142: FN f1,f2;
143: PetscInt i,nq,nisol1,nisol2;
144: PetscScalar *qcoeff,*wr,*wi,*isol1,*isol2;
145: PetscBool flg,avail,rat1,rat2;
148: *rational = PETSC_FALSE;
149: PetscObjectTypeCompare((PetscObject)f,FNRATIONAL,&flg);
150: if (flg) {
151: *rational = PETSC_TRUE;
152: FNRationalGetDenominator(f,&nq,&qcoeff);
153: if (nq>1) {
154: PetscMalloc2(nq-1,&wr,nq-1,&wi);
155: NEPNLEIGSAuxiliarPRootFinder(nq-1,qcoeff,wr,wi,&avail);
156: if (avail) {
157: PetscCalloc1(nq-1,isol);
158: *nisol = 0;
159: for (i=0;i<nq-1;i++)
160: #if !defined(PETSC_USE_COMPLEX)
161: if (wi[i]==0)
162: #endif
163: (*isol)[(*nisol)++] = wr[i];
164: nq = *nisol; *nisol = 0;
165: for (i=0;i<nq;i++) wr[i] = (*isol)[i];
166: NEPNLEIGSAuxiliarRmDuplicates(nq,wr,nisol,*isol,0);
167: PetscFree2(wr,wi);
168: } else { *nisol=0; *isol = NULL; }
169: } else { *nisol = 0; *isol = NULL; }
170: PetscFree(qcoeff);
171: }
172: PetscObjectTypeCompare((PetscObject)f,FNCOMBINE,&flg);
173: if (flg) {
174: FNCombineGetChildren(f,&ctype,&f1,&f2);
175: if (ctype != FN_COMBINE_COMPOSE && ctype != FN_COMBINE_DIVIDE) {
176: NEPNLEIGSFNSingularities(f1,&nisol1,&isol1,&rat1);
177: NEPNLEIGSFNSingularities(f2,&nisol2,&isol2,&rat2);
178: if (nisol1+nisol2>0) {
179: PetscCalloc1(nisol1+nisol2,isol);
180: *nisol = 0;
181: NEPNLEIGSAuxiliarRmDuplicates(nisol1,isol1,nisol,*isol,0);
182: NEPNLEIGSAuxiliarRmDuplicates(nisol2,isol2,nisol,*isol,0);
183: }
184: *rational = (rat1&&rat2)?PETSC_TRUE:PETSC_FALSE;
185: PetscFree(isol1);
186: PetscFree(isol2);
187: }
188: }
189: return(0);
190: }
192: static PetscErrorCode NEPNLEIGSRationalSingularities(NEP nep,PetscInt *ndptx,PetscScalar *dxi,PetscBool *rational)
193: {
195: PetscInt nt,i,nisol;
196: FN f;
197: PetscScalar *isol;
198: PetscBool rat;
201: *rational = PETSC_TRUE;
202: *ndptx = 0;
203: NEPGetSplitOperatorInfo(nep,&nt,NULL);
204: for (i=0;i<nt;i++) {
205: NEPGetSplitOperatorTerm(nep,i,NULL,&f);
206: NEPNLEIGSFNSingularities(f,&nisol,&isol,&rat);
207: if (nisol) {
208: NEPNLEIGSAuxiliarRmDuplicates(nisol,isol,ndptx,dxi,0);
209: PetscFree(isol);
210: }
211: *rational = ((*rational)&&rat)?PETSC_TRUE:PETSC_FALSE;
212: }
213: return(0);
214: }
216: /* Adaptive Anderson-Antoulas algorithm */
217: static PetscErrorCode NEPNLEIGSAAAComputation(NEP nep,PetscInt ndpt,PetscScalar *ds,PetscScalar *F,PetscInt *ndptx,PetscScalar *dxi)
218: {
219: #if defined(SLEPC_MISSING_LAPACK_GGEV) || defined(PETSC_MISSING_LAPACK_GESVD)
221: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GEEV/GESVD - Lapack routines are unavailable");
222: #else
224: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
225: PetscScalar mean=0.0,*z,*f,*C,*A,*VT,*work,*ww,szero=0.0,sone=1.0;
226: PetscScalar *N,*D;
227: PetscReal *S,norm,err,*R;
228: PetscInt i,k,j,idx=0,cont;
229: PetscBLASInt n_,m_,lda_,lwork,info,one=1;
230: #if defined(PETSC_USE_COMPLEX)
231: PetscReal *rwork;
232: #endif
235: PetscBLASIntCast(8*ndpt,&lwork);
236: PetscMalloc5(ndpt,&R,ndpt,&z,ndpt,&f,ndpt*ndpt,&C,ndpt,&ww);
237: PetscMalloc6(ndpt*ndpt,&A,ndpt,&S,ndpt*ndpt,&VT,lwork,&work,ndpt,&D,ndpt,&N);
238: #if defined(PETSC_USE_COMPLEX)
239: PetscMalloc1(8*ndpt,&rwork);
240: #endif
241: norm = 0.0;
242: for (i=0;i<ndpt;i++) {
243: mean += F[i];
244: norm = PetscMax(PetscAbsScalar(F[i]),norm);
245: }
246: mean /= ndpt;
247: PetscBLASIntCast(ndpt,&lda_);
248: for (i=0;i<ndpt;i++) R[i] = PetscAbsScalar(F[i]-mean);
249: /* next support point */
250: err = 0.0;
251: for (i=0;i<ndpt;i++) if (R[i]>=err) {idx = i; err = R[i];}
252: for (k=0;k<ndpt-1;k++) {
253: z[k] = ds[idx]; f[k] = F[idx]; R[idx] = -1.0;
254: /* next column of Cauchy matrix */
255: for (i=0;i<ndpt;i++) {
256: C[i+k*ndpt] = 1.0/(ds[i]-ds[idx]);
257: }
259: PetscMemzero(A,ndpt*ndpt*sizeof(PetscScalar));
260: cont = 0;
261: for (i=0;i<ndpt;i++) {
262: if (R[i]!=-1.0) {
263: for(j=0;j<=k;j++)A[cont+j*ndpt] = C[i+j*ndpt]*F[i]-C[i+j*ndpt]*f[j];
264: cont++;
265: }
266: }
267: PetscBLASIntCast(cont,&m_);
268: PetscBLASIntCast(k+1,&n_);
269: #if defined(PETSC_USE_COMPLEX)
270: PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","A",&m_,&n_,A,&lda_,S,NULL,&lda_,VT,&lda_,work,&lwork,rwork,&info));
271: #else
272: PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","A",&m_,&n_,A,&lda_,S,NULL,&lda_,VT,&lda_,work,&lwork,&info));
273: #endif
274: SlepcCheckLapackInfo("gesvd",info);
275: for (i=0;i<=k;i++) {
276: ww[i] = PetscConj(VT[i*ndpt+k]);
277: D[i] = ww[i]*f[i];
278: }
279: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&lda_,&n_,&sone,C,&lda_,D,&one,&szero,N,&one));
280: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&lda_,&n_,&sone,C,&lda_,ww,&one,&szero,D,&one));
281: for (i=0;i<ndpt;i++) if (R[i]>=0) R[i] = PetscAbsScalar(F[i]-N[i]/D[i]);
282: /* next support point */
283: err = 0.0;
284: for (i=0;i<ndpt;i++) if (R[i]>=err) {idx = i; err = R[i];}
285: if (err <= ctx->ddtol*norm) break;
286: }
288: if (k==ndpt-1) SETERRQ(PetscObjectComm((PetscObject)nep),1,"Failed to determine singularities automatically in general problem");
289: /* poles */
290: PetscMemzero(C,ndpt*ndpt*sizeof(PetscScalar));
291: PetscMemzero(A,ndpt*ndpt*sizeof(PetscScalar));
292: for (i=0;i<=k;i++) {
293: C[i+ndpt*i] = 1.0;
294: A[(i+1)*ndpt] = ww[i];
295: A[i+1] = 1.0;
296: A[i+1+(i+1)*ndpt] = z[i];
297: }
298: C[0] = 0.0; C[k+1+(k+1)*ndpt] = 1.0;
299: n_++;
300: #if defined(PETSC_USE_COMPLEX)
301: PetscStackCallBLAS("LAPACKggev",LAPACKggev_("N","N",&n_,A,&lda_,C,&lda_,D,N,NULL,&lda_,NULL,&lda_,work,&lwork,rwork,&info));
302: #else
303: PetscStackCallBLAS("LAPACKggev",LAPACKggev_("N","N",&n_,A,&lda_,C,&lda_,D,VT,N,NULL,&lda_,NULL,&lda_,work,&lwork,&info));
304: #endif
305: SlepcCheckLapackInfo("ggev",info);
306: cont = 0.0;
307: for (i=0;i<n_;i++) if (N[i]!=0.0) {
308: dxi[cont++] = D[i]/N[i];
309: }
310: *ndptx = cont;
311: PetscFree5(R,z,f,C,ww);
312: PetscFree6(A,S,VT,work,D,N);
313: #if defined(PETSC_USE_COMPLEX)
314: PetscFree(rwork);
315: #endif
316: #endif
317: return(0);
318: }
320: /* Singularities using Adaptive Anderson-Antoulas algorithm */
321: static PetscErrorCode NEPNLEIGSAAASingularities(NEP nep,PetscInt ndpt,PetscScalar *ds,PetscInt *ndptx,PetscScalar *dxi)
322: {
324: Vec u,v,w;
325: PetscRandom rand;
326: PetscScalar *F,*isol;
327: PetscInt i,k,nisol,nt;
328: Mat T;
329: FN f;
332: PetscMalloc1(ndpt,&F);
333: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
334: PetscMalloc1(ndpt,&isol);
335: *ndptx = 0;
336: NEPGetSplitOperatorInfo(nep,&nt,NULL);
337: nisol = *ndptx;
338: for (k=0;k<nt;k++) {
339: NEPGetSplitOperatorTerm(nep,k,NULL,&f);
340: for (i=0;i<ndpt;i++) {
341: FNEvaluateFunction(f,ds[i],&F[i]);
342: }
343: NEPNLEIGSAAAComputation(nep,ndpt,ds,F,&nisol,isol);
344: if (nisol) {
345: NEPNLEIGSAuxiliarRmDuplicates(nisol,isol,ndptx,dxi,ndpt);
346: }
347: }
348: PetscFree(isol);
349: } else {
350: MatCreateVecs(nep->function,&u,NULL);
351: VecDuplicate(u,&v);
352: VecDuplicate(u,&w);
353: BVGetRandomContext(nep->V,&rand);
354: VecSetRandom(u,rand);
355: VecNormalize(u,NULL);
356: VecSetRandom(v,rand);
357: VecNormalize(v,NULL);
358: T = nep->function;
359: for (i=0;i<ndpt;i++) {
360: NEPComputeFunction(nep,ds[i],T,T);
361: MatMult(T,v,w);
362: VecDot(w,u,&F[i]);
363: }
364: NEPNLEIGSAAAComputation(nep,ndpt,ds,F,ndptx,dxi);
365: VecDestroy(&u);
366: VecDestroy(&v);
367: VecDestroy(&w);
368: }
369: PetscFree(F);
370: return(0);
371: }
373: static PetscErrorCode NEPNLEIGSLejaBagbyPoints(NEP nep)
374: {
376: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
377: PetscInt i,k,ndpt=NDPOINTS,ndptx=NDPOINTS;
378: PetscScalar *ds,*dsi,*dxi,*nrs,*nrxi,*s=ctx->s,*xi=ctx->xi,*beta=ctx->beta;
379: PetscReal maxnrs,minnrxi;
380: PetscBool rational;
381: #if !defined(PETSC_USE_COMPLEX)
382: PetscReal a,b,h;
383: #endif
386: if (!ctx->computesingularities && nep->problem_type!=NEP_RATIONAL) ndpt = ndptx = LBPOINTS;
387: PetscMalloc5(ndpt+1,&ds,ndpt+1,&dsi,ndpt,&dxi,ndpt+1,&nrs,ndpt,&nrxi);
389: /* Discretize the target region boundary */
390: RGComputeContour(nep->rg,ndpt,ds,dsi);
391: #if !defined(PETSC_USE_COMPLEX)
392: for (i=0;i<ndpt;i++) if (dsi[i]!=0.0) break;
393: if (i<ndpt) {
394: if (nep->problem_type==NEP_RATIONAL) {
395: /* Select a segment in the real axis */
396: RGComputeBoundingBox(nep->rg,&a,&b,NULL,NULL);
397: if (a<=-PETSC_MAX_REAL || b>=PETSC_MAX_REAL) SETERRQ(PetscObjectComm((PetscObject)nep),1,"NLEIGS requires a bounded target set");
398: h = (b-a)/ndpt;
399: for (i=0;i<ndpt;i++) {ds[i] = a+h*i; dsi[i] = 0.0;}
400: } else SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"NLEIGS with real arithmetic requires the target set to be included in the real axis");
401: }
402: #endif
403: /* Discretize the singularity region */
404: if (ctx->computesingularities) {
405: (ctx->computesingularities)(nep,&ndptx,dxi,ctx->singularitiesctx);
406: } else {
407: if (nep->problem_type==NEP_RATIONAL) {
408: NEPNLEIGSRationalSingularities(nep,&ndptx,dxi,&rational);
409: if (!rational) SETERRQ(PetscObjectComm((PetscObject)nep),1,"Failed to determine singularities automatically in rational problem; consider solving the problem as general");
410: } else {
411: /* AAA algorithm */
412: NEPNLEIGSAAASingularities(nep,ndpt,ds,&ndptx,dxi);
413: }
414: }
415: /* Look for Leja-Bagby points in the discretization sets */
416: s[0] = ds[0];
417: xi[0] = (ndptx>0)?dxi[0]:PETSC_INFINITY;
418: 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]));
419: beta[0] = 1.0; /* scaling factors are also computed here */
420: for (i=0;i<ndpt;i++) {
421: nrs[i] = 1.0;
422: nrxi[i] = 1.0;
423: }
424: for (k=1;k<ctx->ddmaxit;k++) {
425: maxnrs = 0.0;
426: minnrxi = PETSC_MAX_REAL;
427: for (i=0;i<ndpt;i++) {
428: nrs[i] *= ((ds[i]-s[k-1])/(1.0-ds[i]/xi[k-1]))/beta[k-1];
429: if (PetscAbsScalar(nrs[i])>maxnrs) {maxnrs = PetscAbsScalar(nrs[i]); s[k] = ds[i];}
430: }
431: if (ndptx>k) {
432: for (i=1;i<ndptx;i++) {
433: nrxi[i] *= ((dxi[i]-s[k-1])/(1.0-dxi[i]/xi[k-1]))/beta[k-1];
434: if (PetscAbsScalar(nrxi[i])<minnrxi) {minnrxi = PetscAbsScalar(nrxi[i]); xi[k] = dxi[i];}
435: }
436: 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]));
437: } else xi[k] = PETSC_INFINITY;
438: beta[k] = maxnrs;
439: }
440: PetscFree5(ds,dsi,dxi,nrs,nrxi);
441: return(0);
442: }
444: PetscErrorCode NEPNLEIGSEvalNRTFunct(NEP nep,PetscInt k,PetscScalar sigma,PetscScalar *b)
445: {
446: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
447: PetscInt i;
448: PetscScalar *beta=ctx->beta,*s=ctx->s,*xi=ctx->xi;
451: b[0] = 1.0/beta[0];
452: for (i=0;i<k;i++) {
453: b[i+1] = ((sigma-s[i])*b[i])/(beta[i+1]*(1.0-sigma/xi[i]));
454: }
455: return(0);
456: }
458: static PetscErrorCode MatMult_Fun(Mat A,Vec x,Vec y)
459: {
461: ShellMatCtx *ctx;
462: PetscInt i;
465: MatShellGetContext(A,(void**)&ctx);
466: MatMult(ctx->A[0],x,y);
467: if (ctx->coeff[0]!=1.0) { VecScale(y,ctx->coeff[0]); }
468: for (i=1;i<ctx->nmat;i++) {
469: MatMult(ctx->A[i],x,ctx->t);
470: VecAXPY(y,ctx->coeff[i],ctx->t);
471: }
472: return(0);
473: }
475: static PetscErrorCode MatMultTranspose_Fun(Mat A,Vec x,Vec y)
476: {
478: ShellMatCtx *ctx;
479: PetscInt i;
482: MatShellGetContext(A,(void**)&ctx);
483: MatMultTranspose(ctx->A[0],x,y);
484: if (ctx->coeff[0]!=1.0) { VecScale(y,ctx->coeff[0]); }
485: for (i=1;i<ctx->nmat;i++) {
486: MatMultTranspose(ctx->A[i],x,ctx->t);
487: VecAXPY(y,ctx->coeff[i],ctx->t);
488: }
489: return(0);
490: }
492: static PetscErrorCode MatGetDiagonal_Fun(Mat A,Vec diag)
493: {
495: ShellMatCtx *ctx;
496: PetscInt i;
499: MatShellGetContext(A,(void**)&ctx);
500: MatGetDiagonal(ctx->A[0],diag);
501: if (ctx->coeff[0]!=1.0) { VecScale(diag,ctx->coeff[0]); }
502: for (i=1;i<ctx->nmat;i++) {
503: MatGetDiagonal(ctx->A[i],ctx->t);
504: VecAXPY(diag,ctx->coeff[i],ctx->t);
505: }
506: return(0);
507: }
509: static PetscErrorCode MatDuplicate_Fun(Mat A,MatDuplicateOption op,Mat *B)
510: {
511: PetscInt n,i;
512: ShellMatCtx *ctxnew,*ctx;
513: void (*fun)(void);
517: MatShellGetContext(A,(void**)&ctx);
518: PetscNew(&ctxnew);
519: ctxnew->nmat = ctx->nmat;
520: ctxnew->maxnmat = ctx->maxnmat;
521: PetscMalloc2(ctxnew->maxnmat,&ctxnew->A,ctxnew->maxnmat,&ctxnew->coeff);
522: for (i=0;i<ctx->nmat;i++) {
523: PetscObjectReference((PetscObject)ctx->A[i]);
524: ctxnew->A[i] = ctx->A[i];
525: ctxnew->coeff[i] = ctx->coeff[i];
526: }
527: MatGetSize(ctx->A[0],&n,NULL);
528: VecDuplicate(ctx->t,&ctxnew->t);
529: MatCreateShell(PETSC_COMM_WORLD,n,n,n,n,(void*)ctxnew,B);
530: MatShellSetManageScalingShifts(*B);
531: MatShellGetOperation(A,MATOP_MULT,&fun);
532: MatShellSetOperation(*B,MATOP_MULT,fun);
533: MatShellGetOperation(A,MATOP_MULT_TRANSPOSE,&fun);
534: MatShellSetOperation(*B,MATOP_MULT_TRANSPOSE,fun);
535: MatShellGetOperation(A,MATOP_GET_DIAGONAL,&fun);
536: MatShellSetOperation(*B,MATOP_GET_DIAGONAL,fun);
537: MatShellGetOperation(A,MATOP_DUPLICATE,&fun);
538: MatShellSetOperation(*B,MATOP_DUPLICATE,fun);
539: MatShellGetOperation(A,MATOP_DESTROY,&fun);
540: MatShellSetOperation(*B,MATOP_DESTROY,fun);
541: MatShellGetOperation(A,MATOP_AXPY,&fun);
542: MatShellSetOperation(*B,MATOP_AXPY,fun);
543: return(0);
544: }
546: static PetscErrorCode MatDestroy_Fun(Mat A)
547: {
548: ShellMatCtx *ctx;
550: PetscInt i;
553: if (A) {
554: MatShellGetContext(A,(void**)&ctx);
555: for (i=0;i<ctx->nmat;i++) {
556: MatDestroy(&ctx->A[i]);
557: }
558: VecDestroy(&ctx->t);
559: PetscFree2(ctx->A,ctx->coeff);
560: PetscFree(ctx);
561: }
562: return(0);
563: }
565: static PetscErrorCode MatAXPY_Fun(Mat Y,PetscScalar a,Mat X,MatStructure str)
566: {
567: ShellMatCtx *ctxY,*ctxX;
569: PetscInt i,j;
570: PetscBool found;
573: MatShellGetContext(Y,(void**)&ctxY);
574: MatShellGetContext(X,(void**)&ctxX);
575: for (i=0;i<ctxX->nmat;i++) {
576: found = PETSC_FALSE;
577: for (j=0;!found&&j<ctxY->nmat;j++) {
578: if (ctxX->A[i]==ctxY->A[j]) {
579: found = PETSC_TRUE;
580: ctxY->coeff[j] += a*ctxX->coeff[i];
581: }
582: }
583: if (!found) {
584: ctxY->coeff[ctxY->nmat] = a*ctxX->coeff[i];
585: ctxY->A[ctxY->nmat++] = ctxX->A[i];
586: PetscObjectReference((PetscObject)ctxX->A[i]);
587: }
588: }
589: return(0);
590: }
592: static PetscErrorCode MatScale_Fun(Mat M,PetscScalar a)
593: {
594: ShellMatCtx *ctx;
596: PetscInt i;
599: MatShellGetContext(M,(void**)&ctx);
600: for (i=0;i<ctx->nmat;i++) ctx->coeff[i] *= a;
601: return(0);
602: }
604: static PetscErrorCode NLEIGSMatToMatShellArray(Mat M,Mat *Ms,PetscInt maxnmat)
605: {
607: ShellMatCtx *ctx;
608: PetscInt n;
609: PetscBool has;
612: MatHasOperation(M,MATOP_DUPLICATE,&has);
613: if (!has) SETERRQ(PetscObjectComm((PetscObject)M),1,"MatDuplicate operation required");
614: PetscNew(&ctx);
615: ctx->maxnmat = maxnmat;
616: PetscMalloc2(ctx->maxnmat,&ctx->A,ctx->maxnmat,&ctx->coeff);
617: MatDuplicate(M,MAT_COPY_VALUES,&ctx->A[0]);
618: ctx->nmat = 1;
619: ctx->coeff[0] = 1.0;
620: MatCreateVecs(M,&ctx->t,NULL);
621: MatGetSize(M,&n,NULL);
622: MatCreateShell(PetscObjectComm((PetscObject)M),n,n,n,n,(void*)ctx,Ms);
623: MatShellSetManageScalingShifts(*Ms);
624: MatShellSetOperation(*Ms,MATOP_MULT,(void(*)(void))MatMult_Fun);
625: MatShellSetOperation(*Ms,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_Fun);
626: MatShellSetOperation(*Ms,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Fun);
627: MatShellSetOperation(*Ms,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_Fun);
628: MatShellSetOperation(*Ms,MATOP_DESTROY,(void(*)(void))MatDestroy_Fun);
629: MatShellSetOperation(*Ms,MATOP_AXPY,(void(*)(void))MatAXPY_Fun);
630: MatShellSetOperation(*Ms,MATOP_SCALE,(void(*)(void))MatScale_Fun);
631: return(0);
632: }
634: static PetscErrorCode NEPNLEIGSNormEstimation(NEP nep,Mat M,PetscReal *norm,Vec *w)
635: {
636: PetscScalar *z,*x,*y;
637: PetscReal tr;
638: Vec X=w[0],Y=w[1];
639: PetscInt n,i;
640: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
641: PetscRandom rand;
645: if (!ctx->vrn) {
646: /* generate a random vector with normally distributed entries with the Box-Muller transform */
647: BVGetRandomContext(nep->V,&rand);
648: MatCreateVecs(M,&ctx->vrn,NULL);
649: VecSetRandom(X,rand);
650: VecSetRandom(Y,rand);
651: VecGetLocalSize(ctx->vrn,&n);
652: VecGetArray(ctx->vrn,&z);
653: VecGetArray(X,&x);
654: VecGetArray(Y,&y);
655: for (i=0;i<n;i++) {
656: #if defined(PETSC_USE_COMPLEX)
657: z[i] = PetscSqrtReal(-2.0*PetscLogReal(PetscRealPart(x[i])))*PetscCosReal(2.0*PETSC_PI*PetscRealPart(y[i]));
658: z[i] += PETSC_i*(PetscSqrtReal(-2.0*PetscLogReal(PetscImaginaryPart(x[i])))*PetscCosReal(2.0*PETSC_PI*PetscImaginaryPart(y[i])));
659: #else
660: z[i] = PetscSqrtReal(-2.0*PetscLogReal(x[i]))*PetscCosReal(2.0*PETSC_PI*y[i]);
661: #endif
662: }
663: VecRestoreArray(ctx->vrn,&z);
664: VecRestoreArray(X,&x);
665: VecRestoreArray(Y,&y);
666: VecNorm(ctx->vrn,NORM_2,&tr);
667: VecScale(ctx->vrn,1/tr);
668: }
669: /* matrix-free norm estimator of Ipsen http://www4.ncsu.edu/~ipsen/ps/slides_ima.pdf */
670: MatGetSize(M,&n,NULL);
671: MatMult(M,ctx->vrn,X);
672: VecNorm(X,NORM_2,norm);
673: *norm *= PetscSqrtReal((PetscReal)n);
674: return(0);
675: }
677: static PetscErrorCode NEPNLEIGSDividedDifferences_split(NEP nep)
678: {
680: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
681: PetscInt k,j,i,maxnmat,nmax;
682: PetscReal norm0,norm,*matnorm;
683: PetscScalar *s=ctx->s,*beta=ctx->beta,*xi=ctx->xi,*b,alpha,*coeffs,*pK,*pH,sone=1.0;
684: Mat T,Ts,K,H;
685: PetscBool shell,hasmnorm=PETSC_FALSE,matrix=PETSC_TRUE;
686: PetscBLASInt n_;
689: nmax = ctx->ddmaxit;
690: PetscMalloc1(nep->nt*nmax,&ctx->coeffD);
691: PetscMalloc3(nmax+1,&b,nmax+1,&coeffs,nep->nt,&matnorm);
692: for (j=0;j<nep->nt;j++) {
693: MatHasOperation(nep->A[j],MATOP_NORM,&hasmnorm);
694: if (!hasmnorm) break;
695: MatNorm(nep->A[j],NORM_INFINITY,matnorm+j);
696: }
697: /* Try matrix functions scheme */
698: PetscCalloc2(nmax*nmax,&pK,nmax*nmax,&pH);
699: for (i=0;i<nmax-1;i++) {
700: pK[(nmax+1)*i] = 1.0;
701: pK[(nmax+1)*i+1] = beta[i+1]/xi[i];
702: pH[(nmax+1)*i] = s[i];
703: pH[(nmax+1)*i+1] = beta[i+1];
704: }
705: pH[nmax*nmax-1] = s[nmax-1];
706: pK[nmax*nmax-1] = 1.0;
707: PetscBLASIntCast(nmax,&n_);
708: PetscStackCallBLAS("BLAStrsm",BLAStrsm_("R","L","N","U",&n_,&n_,&sone,pK,&n_,pH,&n_));
709: /* The matrix to be used is in H. K will be a work-space matrix */
710: MatCreateSeqDense(PETSC_COMM_SELF,nmax,nmax,pH,&H);
711: MatCreateSeqDense(PETSC_COMM_SELF,nmax,nmax,pK,&K);
712: for (j=0;matrix&&j<nep->nt;j++) {
713: PetscPushErrorHandler(PetscIgnoreErrorHandler,NULL);
714: FNEvaluateFunctionMat(nep->f[j],H,K);
715: PetscPopErrorHandler();
716: if (!ierr) {
717: for (i=0;i<nmax;i++) { ctx->coeffD[j+i*nep->nt] = pK[i]*beta[0]; }
718: } else {
719: matrix = PETSC_FALSE;
720: PetscFPTrapPop();
721: }
722: }
723: MatDestroy(&H);
724: MatDestroy(&K);
725: if (!matrix) {
726: for (j=0;j<nep->nt;j++) {
727: FNEvaluateFunction(nep->f[j],s[0],ctx->coeffD+j);
728: ctx->coeffD[j] *= beta[0];
729: }
730: }
731: if (hasmnorm) {
732: norm0 = 0.0;
733: for (j=0;j<nep->nt;j++) norm0 += matnorm[j]*PetscAbsScalar(ctx->coeffD[j]);
734: } else {
735: norm0 = 0.0;
736: for (j=0;j<nep->nt;j++) norm0 = PetscMax(PetscAbsScalar(ctx->coeffD[j]),norm0);
737: }
738: ctx->nmat = ctx->ddmaxit;
739: for (k=1;k<ctx->ddmaxit;k++) {
740: if (!matrix) {
741: NEPNLEIGSEvalNRTFunct(nep,k,s[k],b);
742: for (i=0;i<nep->nt;i++) {
743: FNEvaluateFunction(nep->f[i],s[k],ctx->coeffD+k*nep->nt+i);
744: for (j=0;j<k;j++) {
745: ctx->coeffD[k*nep->nt+i] -= b[j]*ctx->coeffD[i+nep->nt*j];
746: }
747: ctx->coeffD[k*nep->nt+i] /= b[k];
748: }
749: }
750: if (hasmnorm) {
751: norm = 0.0;
752: for (j=0;j<nep->nt;j++) norm += matnorm[j]*PetscAbsScalar(ctx->coeffD[k*nep->nt+j]);
753: } else {
754: norm = 0.0;
755: for (j=0;j<nep->nt;j++) norm = PetscMax(PetscAbsScalar(ctx->coeffD[k*nep->nt+j]),norm);
756: }
757: if (norm/norm0 < ctx->ddtol) {
758: ctx->nmat = k+1;
759: break;
760: }
761: }
762: if (!ctx->ksp) { NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp); }
763: PetscObjectTypeCompare((PetscObject)nep->A[0],MATSHELL,&shell);
764: maxnmat = PetscMax(ctx->ddmaxit,nep->nt);
765: for (i=0;i<ctx->nshiftsw;i++) {
766: NEPNLEIGSEvalNRTFunct(nep,ctx->nmat-1,ctx->shifts[i],coeffs);
767: if (!shell) {
768: MatDuplicate(nep->A[0],MAT_COPY_VALUES,&T);
769: } else {
770: NLEIGSMatToMatShellArray(nep->A[0],&T,maxnmat);
771: }
772: alpha = 0.0;
773: for (j=0;j<ctx->nmat;j++) alpha += coeffs[j]*ctx->coeffD[j*nep->nt];
774: MatScale(T,alpha);
775: for (k=1;k<nep->nt;k++) {
776: alpha = 0.0;
777: for (j=0;j<ctx->nmat;j++) alpha += coeffs[j]*ctx->coeffD[j*nep->nt+k];
778: if (shell) {
779: NLEIGSMatToMatShellArray(nep->A[k],&Ts,maxnmat);
780: }
781: MatAXPY(T,alpha,shell?Ts:nep->A[k],nep->mstr);
782: if (shell) {
783: MatDestroy(&Ts);
784: }
785: }
786: KSPSetOperators(ctx->ksp[i],T,T);
787: KSPSetUp(ctx->ksp[i]);
788: MatDestroy(&T);
789: }
790: PetscFree3(b,coeffs,matnorm);
791: PetscFree2(pK,pH);
792: return(0);
793: }
795: static PetscErrorCode NEPNLEIGSDividedDifferences_callback(NEP nep)
796: {
798: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
799: PetscInt k,j,i,maxnmat;
800: PetscReal norm0,norm;
801: PetscScalar *s=ctx->s,*beta=ctx->beta,*b,*coeffs;
802: Mat *D=ctx->D,T;
803: PetscBool shell,has,vec=PETSC_FALSE;
804: Vec w[2];
807: PetscMalloc2(ctx->ddmaxit+1,&b,ctx->ddmaxit+1,&coeffs);
808: T = nep->function;
809: NEPComputeFunction(nep,s[0],T,T);
810: PetscObjectTypeCompare((PetscObject)T,MATSHELL,&shell);
811: maxnmat = PetscMax(ctx->ddmaxit,nep->nt);
812: if (!shell) {
813: MatDuplicate(T,MAT_COPY_VALUES,&D[0]);
814: } else {
815: NLEIGSMatToMatShellArray(T,&D[0],maxnmat);
816: }
817: if (beta[0]!=1.0) {
818: MatScale(D[0],1.0/beta[0]);
819: }
820: MatHasOperation(D[0],MATOP_NORM,&has);
821: if (has) {
822: MatNorm(D[0],NORM_FROBENIUS,&norm0);
823: } else {
824: MatCreateVecs(D[0],NULL,&w[0]);
825: VecDuplicate(w[0],&w[1]);
826: vec = PETSC_TRUE;
827: NEPNLEIGSNormEstimation(nep,D[0],&norm0,w);
828: }
829: ctx->nmat = ctx->ddmaxit;
830: for (k=1;k<ctx->ddmaxit;k++) {
831: NEPNLEIGSEvalNRTFunct(nep,k,s[k],b);
832: NEPComputeFunction(nep,s[k],T,T);
833: if (!shell) {
834: MatDuplicate(T,MAT_COPY_VALUES,&D[k]);
835: } else {
836: NLEIGSMatToMatShellArray(T,&D[k],maxnmat);
837: }
838: for (j=0;j<k;j++) {
839: MatAXPY(D[k],-b[j],D[j],nep->mstr);
840: }
841: MatScale(D[k],1.0/b[k]);
842: MatHasOperation(D[k],MATOP_NORM,&has);
843: if (has) {
844: MatNorm(D[k],NORM_FROBENIUS,&norm);
845: } else {
846: if(!vec) {
847: MatCreateVecs(D[k],NULL,&w[0]);
848: VecDuplicate(w[0],&w[1]);
849: vec = PETSC_TRUE;
850: }
851: NEPNLEIGSNormEstimation(nep,D[k],&norm,w);
852: }
853: if (norm/norm0 < ctx->ddtol && k>1) {
854: ctx->nmat = k+1;
855: break;
856: }
857: }
858: if (!ctx->ksp) { NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp); }
859: for (i=0;i<ctx->nshiftsw;i++) {
860: NEPNLEIGSEvalNRTFunct(nep,ctx->nmat-1,ctx->shifts[i],coeffs);
861: MatDuplicate(ctx->D[0],MAT_COPY_VALUES,&T);
862: if (coeffs[0]!=1.0) { MatScale(T,coeffs[0]); }
863: for (j=1;j<ctx->nmat;j++) {
864: MatAXPY(T,coeffs[j],ctx->D[j],nep->mstr);
865: }
866: KSPSetOperators(ctx->ksp[i],T,T);
867: KSPSetUp(ctx->ksp[i]);
868: MatDestroy(&T);
869: }
870: PetscFree2(b,coeffs);
871: if (vec) {
872: VecDestroy(&w[0]);
873: VecDestroy(&w[1]);
874: }
875: return(0);
876: }
878: /*
879: NEPKrylovConvergence - This is the analogue to EPSKrylovConvergence.
880: */
881: PetscErrorCode NEPNLEIGSKrylovConvergence(NEP nep,PetscBool getall,PetscInt kini,PetscInt nits,PetscReal betah,PetscScalar betak,PetscInt *kout,Vec *w)
882: {
884: PetscInt k,newk,marker,inside;
885: PetscScalar re,im;
886: PetscReal resnorm,tt;
887: PetscBool istrivial;
888: NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
891: RGIsTrivial(nep->rg,&istrivial);
892: marker = -1;
893: if (nep->trackall) getall = PETSC_TRUE;
894: for (k=kini;k<kini+nits;k++) {
895: /* eigenvalue */
896: re = nep->eigr[k];
897: im = nep->eigi[k];
898: if (!istrivial) {
899: if (!ctx->nshifts) {
900: NEPNLEIGSBackTransform((PetscObject)nep,1,&re,&im);
901: }
902: RGCheckInside(nep->rg,1,&re,&im,&inside);
903: if (marker==-1 && inside<0) marker = k;
904: }
905: newk = k;
906: DSVectors(nep->ds,DS_MAT_X,&newk,&resnorm);
907: tt = (ctx->nshifts)?SlepcAbsEigenvalue(betak-nep->eigr[k]*betah,nep->eigi[k]*betah):betah;
908: resnorm *= PetscAbsReal(tt);
909: /* error estimate */
910: (*nep->converged)(nep,nep->eigr[k],nep->eigi[k],resnorm,&nep->errest[k],nep->convergedctx);
911: if (marker==-1 && nep->errest[k] >= nep->tol) marker = k;
912: if (newk==k+1) {
913: nep->errest[k+1] = nep->errest[k];
914: k++;
915: }
916: if (marker!=-1 && !getall) break;
917: }
918: if (marker!=-1) k = marker;
919: *kout = k;
920: return(0);
921: }
923: PetscErrorCode NEPSetUp_NLEIGS(NEP nep)
924: {
926: PetscInt k,in;
927: PetscScalar zero=0.0;
928: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
929: SlepcSC sc;
930: PetscBool istrivial;
933: NEPSetDimensions_Default(nep,nep->nev,&nep->ncv,&nep->mpd);
934: if (nep->ncv>nep->nev+nep->mpd) SETERRQ(PetscObjectComm((PetscObject)nep),1,"The value of ncv must not be larger than nev+mpd");
935: if (!nep->max_it) nep->max_it = PetscMax(5000,2*nep->n/nep->ncv);
936: if (!ctx->ddmaxit) ctx->ddmaxit = LBPOINTS;
937: RGIsTrivial(nep->rg,&istrivial);
938: if (istrivial) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"NEPNLEIGS requires a nontrivial region defining the target set");
939: if (!nep->which) nep->which = NEP_TARGET_MAGNITUDE;
940: 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()");
942: /* Initialize the NLEIGS context structure */
943: k = ctx->ddmaxit;
944: PetscMalloc4(k,&ctx->s,k,&ctx->xi,k,&ctx->beta,k,&ctx->D);
945: nep->data = ctx;
946: if (nep->tol==PETSC_DEFAULT) nep->tol = SLEPC_DEFAULT_TOL;
947: if (ctx->ddtol==PETSC_DEFAULT) ctx->ddtol = nep->tol/10.0;
948: if (!ctx->keep) ctx->keep = 0.5;
950: /* Compute Leja-Bagby points and scaling values */
951: NEPNLEIGSLejaBagbyPoints(nep);
952: if (nep->problem_type!=NEP_RATIONAL) {
953: RGCheckInside(nep->rg,1,&nep->target,&zero,&in);
954: if (in<0) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"The target is not inside the target set");
955: }
957: /* Compute the divided difference matrices */
958: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
959: NEPNLEIGSDividedDifferences_split(nep);
960: } else {
961: NEPNLEIGSDividedDifferences_callback(nep);
962: }
963: NEPAllocateSolution(nep,ctx->nmat-1);
964: NEPSetWorkVecs(nep,4);
965: if (!ctx->fullbasis) {
966: if (nep->twosided) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Two-sided variant requires the full-basis option, rerun with -nep_nleigs_full_basis");
967: /* set-up DS and transfer split operator functions */
968: DSSetType(nep->ds,ctx->nshifts?DSGNHEP:DSNHEP);
969: DSAllocate(nep->ds,nep->ncv+1);
970: DSGetSlepcSC(nep->ds,&sc);
971: if (!ctx->nshifts) {
972: sc->map = NEPNLEIGSBackTransform;
973: DSSetExtraRow(nep->ds,PETSC_TRUE);
974: }
975: sc->mapobj = (PetscObject)nep;
976: sc->rg = nep->rg;
977: sc->comparison = nep->sc->comparison;
978: sc->comparisonctx = nep->sc->comparisonctx;
979: BVDestroy(&ctx->V);
980: BVCreateTensor(nep->V,ctx->nmat-1,&ctx->V);
981: nep->ops->solve = NEPSolve_NLEIGS;
982: nep->ops->computevectors = NEPComputeVectors_Schur;
983: } else {
984: NEPSetUp_NLEIGS_FullBasis(nep);
985: nep->ops->solve = NEPSolve_NLEIGS_FullBasis;
986: nep->ops->computevectors = NULL;
987: }
988: return(0);
989: }
991: /*
992: Extend the TOAR basis by applying the the matrix operator
993: over a vector which is decomposed on the TOAR way
994: Input:
995: - S,V: define the latest Arnoldi vector (nv vectors in V)
996: Output:
997: - t: new vector extending the TOAR basis
998: - r: temporally coefficients to compute the TOAR coefficients
999: for the new Arnoldi vector
1000: Workspace: t_ (two vectors)
1001: */
1002: 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_)
1003: {
1005: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1006: PetscInt deg=ctx->nmat-1,k,j;
1007: Vec v=t_[0],q=t_[1],w;
1008: PetscScalar *beta=ctx->beta,*s=ctx->s,*xi=ctx->xi,*coeffs,sigma;
1011: if (!ctx->ksp) { NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp); }
1012: sigma = ctx->shifts[idxrktg];
1013: BVSetActiveColumns(nep->V,0,nv);
1014: PetscMalloc1(ctx->nmat,&coeffs);
1015: if (PetscAbsScalar(s[deg-2]-sigma)<100*PETSC_MACHINE_EPSILON) SETERRQ(PETSC_COMM_SELF,1,"Breakdown in NLEIGS");
1016: /* i-part stored in (i-1) position */
1017: for (j=0;j<nv;j++) {
1018: 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);
1019: }
1020: BVSetActiveColumns(W,0,deg);
1021: BVGetColumn(W,deg-1,&w);
1022: BVMultVec(V,1.0/beta[deg],0,w,S+(deg-1)*ls);
1023: BVRestoreColumn(W,deg-1,&w);
1024: BVGetColumn(W,deg-2,&w);
1025: BVMultVec(V,1.0,0.0,w,r+(deg-2)*lr);
1026: BVRestoreColumn(W,deg-2,&w);
1027: for (k=deg-2;k>0;k--) {
1028: if (PetscAbsScalar(s[k-1]-sigma)<100*PETSC_MACHINE_EPSILON) SETERRQ(PETSC_COMM_SELF,1,"Breakdown in NLEIGS");
1029: 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);
1030: BVGetColumn(W,k-1,&w);
1031: BVMultVec(V,1.0,0.0,w,r+(k-1)*lr);
1032: BVRestoreColumn(W,k-1,&w);
1033: }
1034: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
1035: for (j=0;j<ctx->nmat-2;j++) coeffs[j] = ctx->coeffD[nep->nt*j];
1036: coeffs[ctx->nmat-2] = ctx->coeffD[nep->nt*(ctx->nmat-1)];
1037: BVMultVec(W,1.0,0.0,v,coeffs);
1038: MatMult(nep->A[0],v,q);
1039: for (k=1;k<nep->nt;k++) {
1040: for (j=0;j<ctx->nmat-2;j++) coeffs[j] = ctx->coeffD[nep->nt*j+k];
1041: coeffs[ctx->nmat-2] = ctx->coeffD[nep->nt*(ctx->nmat-1)+k];
1042: BVMultVec(W,1.0,0,v,coeffs);
1043: MatMult(nep->A[k],v,t);
1044: VecAXPY(q,1.0,t);
1045: }
1046: KSPSolve(ctx->ksp[idxrktg],q,t);
1047: VecScale(t,-1.0);
1048: } else {
1049: for (k=0;k<deg-1;k++) {
1050: BVGetColumn(W,k,&w);
1051: MatMult(ctx->D[k],w,q);
1052: BVRestoreColumn(W,k,&w);
1053: BVInsertVec(W,k,q);
1054: }
1055: BVGetColumn(W,deg-1,&w);
1056: MatMult(ctx->D[deg],w,q);
1057: BVRestoreColumn(W,k,&w);
1058: BVInsertVec(W,k,q);
1059: for (j=0;j<ctx->nmat-1;j++) coeffs[j] = 1.0;
1060: BVMultVec(W,1.0,0.0,q,coeffs);
1061: KSPSolve(ctx->ksp[idxrktg],q,t);
1062: VecScale(t,-1.0);
1063: }
1064: PetscFree(coeffs);
1065: return(0);
1066: }
1068: /*
1069: Compute TOAR coefficients of the blocks of the new Arnoldi vector computed
1070: */
1071: static PetscErrorCode NEPTOARCoefficients(NEP nep,PetscScalar sigma,PetscInt nv,PetscScalar *S,PetscInt ls,PetscScalar *r,PetscInt lr,PetscScalar *x,PetscScalar *work)
1072: {
1074: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1075: PetscInt k,j,d=ctx->nmat-1;
1076: PetscScalar *t=work;
1079: NEPNLEIGSEvalNRTFunct(nep,d-1,sigma,t);
1080: for (k=0;k<d-1;k++) {
1081: for (j=0;j<=nv;j++) r[k*lr+j] += t[k]*x[j];
1082: }
1083: for (j=0;j<=nv;j++) r[(d-1)*lr+j] = t[d-1]*x[j];
1084: return(0);
1085: }
1087: /*
1088: Compute continuation vector coefficients for the Rational-Krylov run.
1089: dim(work) >= (end-ini)*(end-ini+1) + end+1 + 2*(end-ini+1), dim(t) = end.
1090: */
1091: 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)
1092: {
1093: #if defined(PETSC_MISSING_LAPACK_GEQRF) || defined(SLEPC_MISSING_LAPACK_LARF)
1095: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GEQRF/LARF - Lapack routines are unavailable");
1096: #else
1098: PetscScalar *x,*W,*tau,sone=1.0,szero=0.0;
1099: PetscInt i,j,n1,n,nwu=0;
1100: PetscBLASInt info,n_,n1_,one=1,dim,lds_;
1101: NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
1104: if (!ctx->nshifts || !end) {
1105: t[0] = 1;
1106: PetscMemcpy(cont,S+end*lds,lds*sizeof(PetscScalar));
1107: } else {
1108: n = end-ini;
1109: n1 = n+1;
1110: x = work+nwu;
1111: nwu += end+1;
1112: tau = work+nwu;
1113: nwu += n;
1114: W = work+nwu;
1115: nwu += n1*n;
1116: for (j=ini;j<end;j++) {
1117: for (i=ini;i<=end;i++) W[(j-ini)*n1+i-ini] = K[j*ld+i] -H[j*ld+i]*sigma;
1118: }
1119: PetscBLASIntCast(n,&n_);
1120: PetscBLASIntCast(n1,&n1_);
1121: PetscBLASIntCast(end+1,&dim);
1122: PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&n1_,&n_,W,&n1_,tau,work+nwu,&n1_,&info));
1123: SlepcCheckLapackInfo("geqrf",info);
1124: for (i=0;i<end;i++) t[i] = 0.0;
1125: t[end] = 1.0;
1126: for (j=n-1;j>=0;j--) {
1127: for (i=0;i<ini+j;i++) x[i] = 0.0;
1128: x[ini+j] = 1.0;
1129: for (i=j+1;i<n1;i++) x[i+ini] = W[i+n1*j];
1130: tau[j] = PetscConj(tau[j]);
1131: PetscStackCallBLAS("LAPACKlarf",LAPACKlarf_("L",&dim,&one,x,&one,tau+j,t,&dim,work+nwu));
1132: }
1133: PetscBLASIntCast(lds,&lds_);
1134: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&lds_,&n1_,&sone,S,&lds_,t,&one,&szero,cont,&one));
1135: }
1136: return(0);
1137: #endif
1138: }
1140: /*
1141: Compute a run of Arnoldi iterations
1142: */
1143: PetscErrorCode NEPNLEIGSTOARrun(NEP nep,PetscScalar *K,PetscScalar *H,PetscInt ldh,BV W,PetscInt k,PetscInt *M,PetscBool *breakdown,Vec *t_)
1144: {
1146: NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
1147: PetscInt i,j,m=*M,lwa,deg=ctx->nmat-1,lds,nqt,ld,l;
1148: Vec t;
1149: PetscReal norm;
1150: PetscScalar *x,*work,*tt,sigma,*cont,*S;
1151: PetscBool lindep;
1152: Mat MS;
1155: BVTensorGetFactors(ctx->V,NULL,&MS);
1156: MatDenseGetArray(MS,&S);
1157: BVGetSizes(nep->V,NULL,NULL,&ld);
1158: lds = ld*deg;
1159: BVGetActiveColumns(nep->V,&l,&nqt);
1160: lwa = PetscMax(ld,deg)+(m+1)*(m+1)+4*(m+1);
1161: PetscMalloc4(ld,&x,lwa,&work,m+1,&tt,lds,&cont);
1162: BVSetActiveColumns(ctx->V,0,m);
1163: for (j=k;j<m;j++) {
1164: sigma = ctx->shifts[(++(ctx->idxrk))%ctx->nshiftsw];
1166: /* Continuation vector */
1167: NEPNLEIGS_RKcontinuation(nep,0,j,K,H,ldh,sigma,S,lds,cont,tt,work);
1169: /* apply operator */
1170: BVGetColumn(nep->V,nqt,&t);
1171: NEPTOARExtendBasis(nep,(ctx->idxrk)%ctx->nshiftsw,cont,ld,nqt,W,nep->V,t,S+(j+1)*lds,ld,t_);
1172: BVRestoreColumn(nep->V,nqt,&t);
1174: /* orthogonalize */
1175: BVOrthogonalizeColumn(nep->V,nqt,x,&norm,&lindep);
1176: if (!lindep) {
1177: x[nqt] = norm;
1178: BVScaleColumn(nep->V,nqt,1.0/norm);
1179: nqt++;
1180: } else x[nqt] = 0.0;
1182: NEPTOARCoefficients(nep,sigma,nqt-1,cont,ld,S+(j+1)*lds,ld,x,work);
1184: /* Level-2 orthogonalization */
1185: BVOrthogonalizeColumn(ctx->V,j+1,H+j*ldh,&norm,breakdown);
1186: H[j+1+ldh*j] = norm;
1187: if (ctx->nshifts) {
1188: for (i=0;i<=j;i++) K[i+ldh*j] = sigma*H[i+ldh*j] + tt[i];
1189: K[j+1+ldh*j] = sigma*H[j+1+ldh*j];
1190: }
1191: if (*breakdown) {
1192: *M = j+1;
1193: break;
1194: }
1195: BVScaleColumn(ctx->V,j+1,1.0/norm);
1196: BVSetActiveColumns(nep->V,l,nqt);
1197: }
1198: PetscFree4(x,work,tt,cont);
1199: MatDenseRestoreArray(MS,&S);
1200: BVTensorRestoreFactors(ctx->V,NULL,&MS);
1201: return(0);
1202: }
1204: PetscErrorCode NEPSolve_NLEIGS(NEP nep)
1205: {
1207: NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
1208: PetscInt i,k=0,l,nv=0,ld,lds,ldds,nq,newn;
1209: PetscInt deg=ctx->nmat-1,nconv=0;
1210: PetscScalar *S,*Q,*H,*pU,*K,betak=0,*eigr,*eigi;
1211: PetscReal betah;
1212: PetscBool falselock=PETSC_FALSE,breakdown=PETSC_FALSE;
1213: BV W;
1214: Mat MS,MQ,U;
1217: if (ctx->lock) {
1218: PetscOptionsGetBool(NULL,NULL,"-nep_nleigs_falselocking",&falselock,NULL);
1219: }
1221: BVGetSizes(nep->V,NULL,NULL,&ld);
1222: lds = deg*ld;
1223: DSGetLeadingDimension(nep->ds,&ldds);
1224: if (!ctx->nshifts) {
1225: PetscMalloc2(nep->ncv,&eigr,nep->ncv,&eigi);
1226: } else { eigr = nep->eigr; eigi = nep->eigi; }
1227: BVDuplicateResize(nep->V,PetscMax(nep->nt-1,ctx->nmat-1),&W);
1230: /* clean projected matrix (including the extra-arrow) */
1231: DSGetArray(nep->ds,DS_MAT_A,&H);
1232: PetscMemzero(H,ldds*ldds*sizeof(PetscScalar));
1233: DSRestoreArray(nep->ds,DS_MAT_A,&H);
1234: if (ctx->nshifts) {
1235: DSGetArray(nep->ds,DS_MAT_B,&H);
1236: PetscMemzero(H,ldds*ldds*sizeof(PetscScalar));
1237: DSRestoreArray(nep->ds,DS_MAT_B,&H);
1238: }
1239:
1240: /* Get the starting Arnoldi vector */
1241: BVTensorBuildFirstColumn(ctx->V,nep->nini);
1242:
1243: /* Restart loop */
1244: l = 0;
1245: while (nep->reason == NEP_CONVERGED_ITERATING) {
1246: nep->its++;
1248: /* Compute an nv-step Krylov relation */
1249: nv = PetscMin(nep->nconv+nep->mpd,nep->ncv);
1250: if (ctx->nshifts) { DSGetArray(nep->ds,DS_MAT_A,&K); }
1251: DSGetArray(nep->ds,ctx->nshifts?DS_MAT_B:DS_MAT_A,&H);
1252: NEPNLEIGSTOARrun(nep,K,H,ldds,W,nep->nconv+l,&nv,&breakdown,nep->work);
1253: betah = PetscAbsScalar(H[(nv-1)*ldds+nv]);
1254: DSRestoreArray(nep->ds,ctx->nshifts?DS_MAT_B:DS_MAT_A,&H);
1255: if (ctx->nshifts) {
1256: betak = K[(nv-1)*ldds+nv];
1257: DSRestoreArray(nep->ds,DS_MAT_A,&K);
1258: }
1259: DSSetDimensions(nep->ds,nv,0,nep->nconv,nep->nconv+l);
1260: if (l==0) {
1261: DSSetState(nep->ds,DS_STATE_INTERMEDIATE);
1262: } else {
1263: DSSetState(nep->ds,DS_STATE_RAW);
1264: }
1266: /* Solve projected problem */
1267: DSSolve(nep->ds,nep->eigr,nep->eigi);
1268: DSSort(nep->ds,nep->eigr,nep->eigi,NULL,NULL,NULL);
1269: if (!ctx->nshifts) {
1270: DSUpdateExtraRow(nep->ds);
1271: }
1272: DSSynchronize(nep->ds,nep->eigr,nep->eigi);
1274: /* Check convergence */
1275: NEPNLEIGSKrylovConvergence(nep,PETSC_FALSE,nep->nconv,nv-nep->nconv,betah,betak,&k,nep->work);
1276: (*nep->stopping)(nep,nep->its,nep->max_it,k,nep->nev,&nep->reason,nep->stoppingctx);
1278: /* Update l */
1279: if (nep->reason != NEP_CONVERGED_ITERATING || breakdown) l = 0;
1280: else {
1281: l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
1282: if (!breakdown) {
1283: /* Prepare the Rayleigh quotient for restart */
1284: if (!ctx->nshifts) {
1285: DSTruncate(nep->ds,k+l);
1286: DSGetDimensions(nep->ds,&newn,NULL,NULL,NULL,NULL);
1287: l = newn-k;
1288: } else {
1289: DSGetArray(nep->ds,DS_MAT_Q,&Q);
1290: DSGetArray(nep->ds,DS_MAT_B,&H);
1291: DSGetArray(nep->ds,DS_MAT_A,&K);
1292: for (i=ctx->lock?k:0;i<k+l;i++) {
1293: H[k+l+i*ldds] = betah*Q[nv-1+i*ldds];
1294: K[k+l+i*ldds] = betak*Q[nv-1+i*ldds];
1295: }
1296: DSRestoreArray(nep->ds,DS_MAT_B,&H);
1297: DSRestoreArray(nep->ds,DS_MAT_A,&K);
1298: DSRestoreArray(nep->ds,DS_MAT_Q,&Q);
1299: }
1300: }
1301: }
1302: nconv = k;
1303: if (!ctx->lock && nep->reason == NEP_CONVERGED_ITERATING && !breakdown) { l += k; k = 0; }
1305: /* Update S */
1306: DSGetMat(nep->ds,ctx->nshifts?DS_MAT_Z:DS_MAT_Q,&MQ);
1307: BVMultInPlace(ctx->V,MQ,nep->nconv,k+l);
1308: MatDestroy(&MQ);
1310: /* Copy last column of S */
1311: BVCopyColumn(ctx->V,nv,k+l);
1313: if (breakdown && nep->reason == NEP_CONVERGED_ITERATING) {
1314: /* Stop if breakdown */
1315: PetscInfo2(nep,"Breakdown (it=%D norm=%g)\n",nep->its,(double)betah);
1316: nep->reason = NEP_DIVERGED_BREAKDOWN;
1317: }
1318: if (nep->reason != NEP_CONVERGED_ITERATING) l--;
1319: /* truncate S */
1320: BVGetActiveColumns(nep->V,NULL,&nq);
1321: if (k+l+deg<=nq) {
1322: BVSetActiveColumns(ctx->V,nep->nconv,k+l+1);
1323: if (!falselock && ctx->lock) {
1324: BVTensorCompress(ctx->V,k-nep->nconv);
1325: } else {
1326: BVTensorCompress(ctx->V,0);
1327: }
1328: }
1329: nep->nconv = k;
1330: if (!ctx->nshifts) {
1331: for (i=0;i<nv;i++) { eigr[i] = nep->eigr[i]; eigi[i] = nep->eigi[i]; }
1332: NEPNLEIGSBackTransform((PetscObject)nep,nv,eigr,eigi);
1333: }
1334: NEPMonitor(nep,nep->its,nconv,eigr,eigi,nep->errest,nv);
1335: }
1336: nep->nconv = nconv;
1337: if (nep->nconv>0) {
1338: BVSetActiveColumns(ctx->V,0,nep->nconv);
1339: BVGetActiveColumns(nep->V,NULL,&nq);
1340: BVSetActiveColumns(nep->V,0,nq);
1341: if (nq>nep->nconv) {
1342: BVTensorCompress(ctx->V,nep->nconv);
1343: BVSetActiveColumns(nep->V,0,nep->nconv);
1344: nq = nep->nconv;
1345: }
1346: if (ctx->nshifts) {
1347: DSGetMat(nep->ds,DS_MAT_B,&MQ);
1348: BVMultInPlace(ctx->V,MQ,0,nep->nconv);
1349: MatDestroy(&MQ);
1350: }
1351: BVTensorGetFactors(ctx->V,NULL,&MS);
1352: MatDenseGetArray(MS,&S);
1353: PetscMalloc1(nq*nep->nconv,&pU);
1354: for (i=0;i<nep->nconv;i++) {
1355: PetscMemcpy(pU+i*nq,S+i*lds,nq*sizeof(PetscScalar));
1356: }
1357: MatDenseRestoreArray(MS,&S);
1358: BVTensorRestoreFactors(ctx->V,NULL,&MS);
1359: MatCreateSeqDense(PETSC_COMM_SELF,nq,nep->nconv,pU,&U);
1360: BVSetActiveColumns(nep->V,0,nq);
1361: BVMultInPlace(nep->V,U,0,nep->nconv);
1362: BVSetActiveColumns(nep->V,0,nep->nconv);
1363: MatDestroy(&U);
1364: PetscFree(pU);
1365: /* truncate Schur decomposition and change the state to raw so that
1366: DSVectors() computes eigenvectors from scratch */
1367: DSSetDimensions(nep->ds,nep->nconv,0,0,0);
1368: DSSetState(nep->ds,DS_STATE_RAW);
1369: }
1371: /* Map eigenvalues back to the original problem */
1372: if (!ctx->nshifts) {
1373: NEPNLEIGSBackTransform((PetscObject)nep,nep->nconv,nep->eigr,nep->eigi);
1374: PetscFree2(eigr,eigi);
1375: }
1376: BVDestroy(&W);
1377: return(0);
1378: }
1380: static PetscErrorCode NEPNLEIGSSetSingularitiesFunction_NLEIGS(NEP nep,PetscErrorCode (*fun)(NEP,PetscInt*,PetscScalar*,void*),void *ctx)
1381: {
1382: NEP_NLEIGS *nepctx=(NEP_NLEIGS*)nep->data;
1385: if (fun) nepctx->computesingularities = fun;
1386: if (ctx) nepctx->singularitiesctx = ctx;
1387: nep->state = NEP_STATE_INITIAL;
1388: return(0);
1389: }
1391: /*@C
1392: NEPNLEIGSSetSingularitiesFunction - Sets a user function to compute a discretization
1393: of the singularity set (where T(.) is not analytic).
1395: Logically Collective on NEP
1397: Input Parameters:
1398: + nep - the NEP context
1399: . fun - user function (if NULL then NEP retains any previously set value)
1400: - ctx - [optional] user-defined context for private data for the function
1401: (may be NULL, in which case NEP retains any previously set value)
1403: Calling Sequence of fun:
1404: $ fun(NEP nep,PetscInt *maxnp,PetscScalar *xi,void *ctx)
1406: + nep - the NEP context
1407: . maxnp - on input number of requested points in the discretization (can be set)
1408: . xi - computed values of the discretization
1409: - ctx - optional context, as set by NEPNLEIGSSetSingularitiesFunction()
1411: Notes:
1412: The user-defined function can set a smaller value of maxnp if necessary.
1413: It is wrong to return a larger value.
1415: If the problem type has been set to rational with NEPSetProblemType(),
1416: then it is not necessary to set the singularities explicitly since the
1417: solver will try to determine them automatically.
1419: Level: intermediate
1421: .seealso: NEPNLEIGSGetSingularitiesFunction(), NEPSetProblemType()
1422: @*/
1423: PetscErrorCode NEPNLEIGSSetSingularitiesFunction(NEP nep,PetscErrorCode (*fun)(NEP,PetscInt*,PetscScalar*,void*),void *ctx)
1424: {
1429: PetscTryMethod(nep,"NEPNLEIGSSetSingularitiesFunction_C",(NEP,PetscErrorCode(*)(NEP,PetscInt*,PetscScalar*,void*),void*),(nep,fun,ctx));
1430: return(0);
1431: }
1433: static PetscErrorCode NEPNLEIGSGetSingularitiesFunction_NLEIGS(NEP nep,PetscErrorCode (**fun)(NEP,PetscInt*,PetscScalar*,void*),void **ctx)
1434: {
1435: NEP_NLEIGS *nepctx=(NEP_NLEIGS*)nep->data;
1438: if (fun) *fun = nepctx->computesingularities;
1439: if (ctx) *ctx = nepctx->singularitiesctx;
1440: return(0);
1441: }
1443: /*@C
1444: NEPNLEIGSGetSingularitiesFunction - Returns the Function and optionally the user
1445: provided context for computing a discretization of the singularity set.
1447: Not Collective
1449: Input Parameter:
1450: . nep - the nonlinear eigensolver context
1452: Output Parameters:
1453: + fun - location to put the function (or NULL)
1454: - ctx - location to stash the function context (or NULL)
1456: Level: advanced
1458: .seealso: NEPNLEIGSSetSingularitiesFunction()
1459: @*/
1460: PetscErrorCode NEPNLEIGSGetSingularitiesFunction(NEP nep,PetscErrorCode (**fun)(NEP,PetscInt*,PetscScalar*,void*),void **ctx)
1461: {
1466: PetscUseMethod(nep,"NEPNLEIGSGetSingularitiesFunction_C",(NEP,PetscErrorCode(**)(NEP,PetscInt*,PetscScalar*,void*),void**),(nep,fun,ctx));
1467: return(0);
1468: }
1470: static PetscErrorCode NEPNLEIGSSetRestart_NLEIGS(NEP nep,PetscReal keep)
1471: {
1472: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1475: if (keep==PETSC_DEFAULT) ctx->keep = 0.5;
1476: else {
1477: 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]");
1478: ctx->keep = keep;
1479: }
1480: return(0);
1481: }
1483: /*@
1484: NEPNLEIGSSetRestart - Sets the restart parameter for the NLEIGS
1485: method, in particular the proportion of basis vectors that must be kept
1486: after restart.
1488: Logically Collective on NEP
1490: Input Parameters:
1491: + nep - the nonlinear eigensolver context
1492: - keep - the number of vectors to be kept at restart
1494: Options Database Key:
1495: . -nep_nleigs_restart - Sets the restart parameter
1497: Notes:
1498: Allowed values are in the range [0.1,0.9]. The default is 0.5.
1500: Level: advanced
1502: .seealso: NEPNLEIGSGetRestart()
1503: @*/
1504: PetscErrorCode NEPNLEIGSSetRestart(NEP nep,PetscReal keep)
1505: {
1511: PetscTryMethod(nep,"NEPNLEIGSSetRestart_C",(NEP,PetscReal),(nep,keep));
1512: return(0);
1513: }
1515: static PetscErrorCode NEPNLEIGSGetRestart_NLEIGS(NEP nep,PetscReal *keep)
1516: {
1517: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1520: *keep = ctx->keep;
1521: return(0);
1522: }
1524: /*@
1525: NEPNLEIGSGetRestart - Gets the restart parameter used in the NLEIGS method.
1527: Not Collective
1529: Input Parameter:
1530: . nep - the nonlinear eigensolver context
1532: Output Parameter:
1533: . keep - the restart parameter
1535: Level: advanced
1537: .seealso: NEPNLEIGSSetRestart()
1538: @*/
1539: PetscErrorCode NEPNLEIGSGetRestart(NEP nep,PetscReal *keep)
1540: {
1546: PetscUseMethod(nep,"NEPNLEIGSGetRestart_C",(NEP,PetscReal*),(nep,keep));
1547: return(0);
1548: }
1550: static PetscErrorCode NEPNLEIGSSetLocking_NLEIGS(NEP nep,PetscBool lock)
1551: {
1552: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1555: ctx->lock = lock;
1556: return(0);
1557: }
1559: /*@
1560: NEPNLEIGSSetLocking - Choose between locking and non-locking variants of
1561: the NLEIGS method.
1563: Logically Collective on NEP
1565: Input Parameters:
1566: + nep - the nonlinear eigensolver context
1567: - lock - true if the locking variant must be selected
1569: Options Database Key:
1570: . -nep_nleigs_locking - Sets the locking flag
1572: Notes:
1573: The default is to lock converged eigenpairs when the method restarts.
1574: This behaviour can be changed so that all directions are kept in the
1575: working subspace even if already converged to working accuracy (the
1576: non-locking variant).
1578: Level: advanced
1580: .seealso: NEPNLEIGSGetLocking()
1581: @*/
1582: PetscErrorCode NEPNLEIGSSetLocking(NEP nep,PetscBool lock)
1583: {
1589: PetscTryMethod(nep,"NEPNLEIGSSetLocking_C",(NEP,PetscBool),(nep,lock));
1590: return(0);
1591: }
1593: static PetscErrorCode NEPNLEIGSGetLocking_NLEIGS(NEP nep,PetscBool *lock)
1594: {
1595: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1598: *lock = ctx->lock;
1599: return(0);
1600: }
1602: /*@
1603: NEPNLEIGSGetLocking - Gets the locking flag used in the NLEIGS method.
1605: Not Collective
1607: Input Parameter:
1608: . nep - the nonlinear eigensolver context
1610: Output Parameter:
1611: . lock - the locking flag
1613: Level: advanced
1615: .seealso: NEPNLEIGSSetLocking()
1616: @*/
1617: PetscErrorCode NEPNLEIGSGetLocking(NEP nep,PetscBool *lock)
1618: {
1624: PetscUseMethod(nep,"NEPNLEIGSGetLocking_C",(NEP,PetscBool*),(nep,lock));
1625: return(0);
1626: }
1628: static PetscErrorCode NEPNLEIGSSetInterpolation_NLEIGS(NEP nep,PetscReal tol,PetscInt degree)
1629: {
1631: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1634: if (tol == PETSC_DEFAULT) {
1635: ctx->ddtol = PETSC_DEFAULT;
1636: nep->state = NEP_STATE_INITIAL;
1637: } else {
1638: if (tol <= 0.0) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of tol. Must be > 0");
1639: ctx->ddtol = tol;
1640: }
1641: if (degree == PETSC_DEFAULT || degree == PETSC_DECIDE) {
1642: ctx->ddmaxit = 0;
1643: if (nep->state) { NEPReset(nep); }
1644: nep->state = NEP_STATE_INITIAL;
1645: } else {
1646: if (degree <= 0) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of degree. Must be > 0");
1647: if (ctx->ddmaxit != degree) {
1648: ctx->ddmaxit = degree;
1649: if (nep->state) { NEPReset(nep); }
1650: nep->state = NEP_STATE_INITIAL;
1651: }
1652: }
1653: return(0);
1654: }
1656: /*@
1657: NEPNLEIGSSetInterpolation - Sets the tolerance and maximum degree
1658: when building the interpolation via divided differences.
1660: Logically Collective on NEP
1662: Input Parameters:
1663: + nep - the nonlinear eigensolver context
1664: . tol - tolerance to stop computing divided differences
1665: - degree - maximum degree of interpolation
1667: Options Database Key:
1668: + -nep_nleigs_interpolation_tol <tol> - Sets the tolerance to stop computing divided differences
1669: - -nep_nleigs_interpolation_degree <degree> - Sets the maximum degree of interpolation
1671: Notes:
1672: Use PETSC_DEFAULT for either argument to assign a reasonably good value.
1674: Level: advanced
1676: .seealso: NEPNLEIGSGetInterpolation()
1677: @*/
1678: PetscErrorCode NEPNLEIGSSetInterpolation(NEP nep,PetscReal tol,PetscInt degree)
1679: {
1686: PetscTryMethod(nep,"NEPNLEIGSSetInterpolation_C",(NEP,PetscReal,PetscInt),(nep,tol,degree));
1687: return(0);
1688: }
1690: static PetscErrorCode NEPNLEIGSGetInterpolation_NLEIGS(NEP nep,PetscReal *tol,PetscInt *degree)
1691: {
1692: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1695: if (tol) *tol = ctx->ddtol;
1696: if (degree) *degree = ctx->ddmaxit;
1697: return(0);
1698: }
1700: /*@
1701: NEPNLEIGSGetInterpolation - Gets the tolerance and maximum degree
1702: when building the interpolation via divided differences.
1704: Not Collective
1706: Input Parameter:
1707: . nep - the nonlinear eigensolver context
1709: Output Parameter:
1710: + tol - tolerance to stop computing divided differences
1711: - degree - maximum degree of interpolation
1713: Level: advanced
1715: .seealso: NEPNLEIGSSetInterpolation()
1716: @*/
1717: PetscErrorCode NEPNLEIGSGetInterpolation(NEP nep,PetscReal *tol,PetscInt *degree)
1718: {
1723: PetscTryMethod(nep,"NEPNLEIGSGetInterpolation_C",(NEP,PetscReal*,PetscInt*),(nep,tol,degree));
1724: return(0);
1725: }
1727: static PetscErrorCode NEPNLEIGSSetRKShifts_NLEIGS(NEP nep,PetscInt ns,PetscScalar *shifts)
1728: {
1730: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1731: PetscInt i;
1734: if (ns<=0) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"Number of shifts must be positive");
1735: if (ctx->nshifts) { PetscFree(ctx->shifts); }
1736: for (i=0;i<ctx->nshiftsw;i++) { KSPDestroy(&ctx->ksp[i]); }
1737: PetscFree(ctx->ksp);
1738: ctx->ksp = NULL;
1739: PetscMalloc1(ns,&ctx->shifts);
1740: for (i=0;i<ns;i++) ctx->shifts[i] = shifts[i];
1741: ctx->nshifts = ns;
1742: nep->state = NEP_STATE_INITIAL;
1743: return(0);
1744: }
1746: /*@C
1747: NEPNLEIGSSetRKShifts - Sets a list of shifts to be used in the Rational
1748: Krylov method.
1750: Logically Collective on NEP
1752: Input Parameters:
1753: + nep - the nonlinear eigensolver context
1754: . ns - number of shifts
1755: - shifts - array of scalar values specifying the shifts
1757: Options Database Key:
1758: . -nep_nleigs_rk_shifts - Sets the list of shifts
1760: Notes:
1761: If only one shift is provided, the built subspace built is equivalent to
1762: shift-and-invert Krylov-Schur (provided that the absolute convergence
1763: criterion is used).
1765: In the case of real scalars, complex shifts are not allowed. In the
1766: command line, a comma-separated list of complex values can be provided with
1767: the format [+/-][realnumber][+/-]realnumberi with no spaces, e.g.
1768: -nep_nleigs_rk_shifts 1.0+2.0i,1.5+2.0i,1.0+1.5i
1770: Level: advanced
1772: .seealso: NEPNLEIGSGetRKShifts()
1773: @*/
1774: PetscErrorCode NEPNLEIGSSetRKShifts(NEP nep,PetscInt ns,PetscScalar *shifts)
1775: {
1782: PetscTryMethod(nep,"NEPNLEIGSSetRKShifts_C",(NEP,PetscInt,PetscScalar*),(nep,ns,shifts));
1783: return(0);
1784: }
1786: static PetscErrorCode NEPNLEIGSGetRKShifts_NLEIGS(NEP nep,PetscInt *ns,PetscScalar **shifts)
1787: {
1789: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1790: PetscInt i;
1793: *ns = ctx->nshifts;
1794: if (ctx->nshifts) {
1795: PetscMalloc1(ctx->nshifts,shifts);
1796: for (i=0;i<ctx->nshifts;i++) (*shifts)[i] = ctx->shifts[i];
1797: }
1798: return(0);
1799: }
1801: /*@C
1802: NEPNLEIGSGetRKShifts - Gets the list of shifts used in the Rational
1803: Krylov method.
1805: Not Collective
1807: Input Parameter:
1808: . nep - the nonlinear eigensolver context
1810: Output Parameter:
1811: + ns - number of shifts
1812: - shifts - array of shifts
1814: Level: advanced
1816: .seealso: NEPNLEIGSSetRKShifts()
1817: @*/
1818: PetscErrorCode NEPNLEIGSGetRKShifts(NEP nep,PetscInt *ns,PetscScalar **shifts)
1819: {
1826: PetscTryMethod(nep,"NEPNLEIGSGetRKShifts_C",(NEP,PetscInt*,PetscScalar**),(nep,ns,shifts));
1827: return(0);
1828: }
1830: static PetscErrorCode NEPNLEIGSGetKSPs_NLEIGS(NEP nep,PetscInt *nsolve,KSP **ksp)
1831: {
1833: NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
1834: PetscInt i;
1835: PC pc;
1838: if (!ctx->ksp) {
1839: NEPNLEIGSSetShifts(nep,&ctx->nshiftsw);
1840: PetscMalloc1(ctx->nshiftsw,&ctx->ksp);
1841: for (i=0;i<ctx->nshiftsw;i++) {
1842: KSPCreate(PetscObjectComm((PetscObject)nep),&ctx->ksp[i]);
1843: PetscObjectIncrementTabLevel((PetscObject)ctx->ksp[i],(PetscObject)nep,1);
1844: KSPSetOptionsPrefix(ctx->ksp[i],((PetscObject)nep)->prefix);
1845: KSPAppendOptionsPrefix(ctx->ksp[i],"nep_nleigs_");
1846: PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->ksp[i]);
1847: PetscObjectSetOptions((PetscObject)ctx->ksp[i],((PetscObject)nep)->options);
1848: KSPSetErrorIfNotConverged(ctx->ksp[i],PETSC_TRUE);
1849: KSPSetTolerances(ctx->ksp[i],SLEPC_DEFAULT_TOL,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
1850: KSPGetPC(ctx->ksp[i],&pc);
1851: KSPSetType(ctx->ksp[i],KSPPREONLY);
1852: PCSetType(pc,PCLU);
1853: }
1854: }
1855: if (nsolve) *nsolve = ctx->nshiftsw;
1856: if (ksp) *ksp = ctx->ksp;
1857: return(0);
1858: }
1860: /*@C
1861: NEPNLEIGSGetKSPs - Retrieve the array of linear solver objects associated with
1862: the nonlinear eigenvalue solver.
1864: Not Collective
1866: Input Parameter:
1867: . nep - nonlinear eigenvalue solver
1869: Output Parameters:
1870: + nsolve - number of returned KSP objects
1871: - ksp - array of linear solver object
1873: Notes:
1874: The number of KSP objects is equal to the number of shifts provided by the user,
1875: or 1 if the user did not provide shifts.
1877: Level: advanced
1879: .seealso: NEPNLEIGSSetRKShifts()
1880: @*/
1881: PetscErrorCode NEPNLEIGSGetKSPs(NEP nep,PetscInt *nsolve,KSP **ksp)
1882: {
1887: PetscUseMethod(nep,"NEPNLEIGSGetKSPs_C",(NEP,PetscInt*,KSP**),(nep,nsolve,ksp));
1888: return(0);
1889: }
1891: static PetscErrorCode NEPNLEIGSSetFullBasis_NLEIGS(NEP nep,PetscBool fullbasis)
1892: {
1893: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1896: if (fullbasis!=ctx->fullbasis) {
1897: ctx->fullbasis = fullbasis;
1898: nep->state = NEP_STATE_INITIAL;
1899: nep->useds = PetscNot(fullbasis);
1900: }
1901: return(0);
1902: }
1904: /*@
1905: NEPNLEIGSSetFullBasis - Choose between TOAR-basis (default) and full-basis
1906: variants of the NLEIGS method.
1908: Logically Collective on NEP
1910: Input Parameters:
1911: + nep - the nonlinear eigensolver context
1912: - fullbasis - true if the full-basis variant must be selected
1914: Options Database Key:
1915: . -nep_nleigs_full_basis - Sets the full-basis flag
1917: Notes:
1918: The default is to use a compact representation of the Krylov basis, that is,
1919: V = (I otimes U) S, with a tensor BV. This behaviour can be changed so that
1920: the full basis V is explicitly stored and operated with. This variant is more
1921: expensive in terms of memory and computation, but is necessary in some cases,
1922: particularly for two-sided computations, see NEPSetTwoSided().
1924: In the full-basis variant, the NLEIGS solver uses an EPS object to explicitly
1925: solve the linearized eigenproblem, see NEPNLEIGSGetEPS().
1927: Level: advanced
1929: .seealso: NEPNLEIGSGetFullBasis(), NEPNLEIGSGetEPS(), NEPSetTwoSided(), BVCreateTensor()
1930: @*/
1931: PetscErrorCode NEPNLEIGSSetFullBasis(NEP nep,PetscBool fullbasis)
1932: {
1938: PetscTryMethod(nep,"NEPNLEIGSSetFullBasis_C",(NEP,PetscBool),(nep,fullbasis));
1939: return(0);
1940: }
1942: static PetscErrorCode NEPNLEIGSGetFullBasis_NLEIGS(NEP nep,PetscBool *fullbasis)
1943: {
1944: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1947: *fullbasis = ctx->fullbasis;
1948: return(0);
1949: }
1951: /*@
1952: NEPNLEIGSGetFullBasis - Gets the flag that indicates if NLEIGS is using the
1953: full-basis variant.
1955: Not Collective
1957: Input Parameter:
1958: . nep - the nonlinear eigensolver context
1960: Output Parameter:
1961: . fullbasis - the flag
1963: Level: advanced
1965: .seealso: NEPNLEIGSSetFullBasis()
1966: @*/
1967: PetscErrorCode NEPNLEIGSGetFullBasis(NEP nep,PetscBool *fullbasis)
1968: {
1974: PetscUseMethod(nep,"NEPNLEIGSGetFullBasis_C",(NEP,PetscBool*),(nep,fullbasis));
1975: return(0);
1976: }
1978: #define SHIFTMAX 30
1980: PetscErrorCode NEPSetFromOptions_NLEIGS(PetscOptionItems *PetscOptionsObject,NEP nep)
1981: {
1983: NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
1984: PetscInt i=0,k;
1985: PetscBool flg1,flg2,b;
1986: PetscReal r;
1987: PetscScalar array[SHIFTMAX];
1990: PetscOptionsHead(PetscOptionsObject,"NEP NLEIGS Options");
1992: PetscOptionsReal("-nep_nleigs_restart","Proportion of vectors kept after restart","NEPNLEIGSSetRestart",0.5,&r,&flg1);
1993: if (flg1) { NEPNLEIGSSetRestart(nep,r); }
1995: PetscOptionsBool("-nep_nleigs_locking","Choose between locking and non-locking variants","NEPNLEIGSSetLocking",PETSC_FALSE,&b,&flg1);
1996: if (flg1) { NEPNLEIGSSetLocking(nep,b); }
1998: PetscOptionsBool("-nep_nleigs_full_basis","Choose between TOAR and full-basis variants","NEPNLEIGSSetFullBasis",PETSC_FALSE,&b,&flg1);
1999: if (flg1) { NEPNLEIGSSetFullBasis(nep,b); }
2001: NEPNLEIGSGetInterpolation(nep,&r,&i);
2002: if (!i) i = PETSC_DEFAULT;
2003: PetscOptionsInt("-nep_nleigs_interpolation_degree","Maximum number of terms for interpolation via divided differences","NEPNLEIGSSetInterpolation",i,&i,&flg1);
2004: PetscOptionsReal("-nep_nleigs_interpolation_tol","Tolerance for interpolation via divided differences","NEPNLEIGSSetInterpolation",r,&r,&flg2);
2005: if (flg1 || flg2) { NEPNLEIGSSetInterpolation(nep,r,i); }
2007: k = SHIFTMAX;
2008: for (i=0;i<k;i++) array[i] = 0;
2009: PetscOptionsScalarArray("-nep_nleigs_rk_shifts","Shifts for Rational Krylov","NEPNLEIGSSetRKShifts",array,&k,&flg1);
2010: if (flg1) { NEPNLEIGSSetRKShifts(nep,k,array); }
2012: PetscOptionsTail();
2014: if (!ctx->ksp) { NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp); }
2015: for (i=0;i<ctx->nshiftsw;i++) {
2016: KSPSetFromOptions(ctx->ksp[i]);
2017: }
2019: if (ctx->fullbasis) {
2020: if (!ctx->eps) { NEPNLEIGSGetEPS(nep,&ctx->eps); }
2021: EPSSetFromOptions(ctx->eps);
2022: }
2023: return(0);
2024: }
2026: PetscErrorCode NEPView_NLEIGS(NEP nep,PetscViewer viewer)
2027: {
2029: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
2030: PetscBool isascii;
2031: PetscInt i;
2032: char str[50];
2035: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
2036: if (isascii) {
2037: PetscViewerASCIIPrintf(viewer," %d%% of basis vectors kept after restart\n",(int)(100*ctx->keep));
2038: if (ctx->fullbasis) {
2039: PetscViewerASCIIPrintf(viewer," using the full-basis variant\n");
2040: } else {
2041: PetscViewerASCIIPrintf(viewer," using the %slocking variant\n",ctx->lock?"":"non-");
2042: }
2043: PetscViewerASCIIPrintf(viewer," divided difference terms: used=%D, max=%D\n",ctx->nmat,ctx->ddmaxit);
2044: PetscViewerASCIIPrintf(viewer," tolerance for divided difference convergence: %g\n",(double)ctx->ddtol);
2045: if (ctx->nshifts) {
2046: PetscViewerASCIIPrintf(viewer," RK shifts: ");
2047: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
2048: for (i=0;i<ctx->nshifts;i++) {
2049: SlepcSNPrintfScalar(str,50,ctx->shifts[i],PETSC_FALSE);
2050: PetscViewerASCIIPrintf(viewer,"%s%s",str,(i<ctx->nshifts-1)?",":"");
2051: }
2052: PetscViewerASCIIPrintf(viewer,"\n");
2053: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
2054: }
2055: if (!ctx->ksp) { NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp); }
2056: PetscViewerASCIIPushTab(viewer);
2057: KSPView(ctx->ksp[0],viewer);
2058: PetscViewerASCIIPopTab(viewer);
2059: if (ctx->fullbasis) {
2060: if (!ctx->eps) { NEPNLEIGSGetEPS(nep,&ctx->eps); }
2061: PetscViewerASCIIPushTab(viewer);
2062: EPSView(ctx->eps,viewer);
2063: PetscViewerASCIIPopTab(viewer);
2064: }
2065: }
2066: return(0);
2067: }
2069: PetscErrorCode NEPReset_NLEIGS(NEP nep)
2070: {
2072: PetscInt k;
2073: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
2076: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
2077: PetscFree(ctx->coeffD);
2078: } else {
2079: for (k=0;k<ctx->nmat;k++) { MatDestroy(&ctx->D[k]); }
2080: }
2081: PetscFree4(ctx->s,ctx->xi,ctx->beta,ctx->D);
2082: for (k=0;k<ctx->nshiftsw;k++) { KSPReset(ctx->ksp[k]); }
2083: if (ctx->vrn) {
2084: VecDestroy(&ctx->vrn);
2085: }
2086: if (ctx->fullbasis) {
2087: MatDestroy(&ctx->A);
2088: EPSReset(ctx->eps);
2089: for (k=0;k<4;k++) { VecDestroy(&ctx->w[k]); }
2090: }
2091: return(0);
2092: }
2094: PetscErrorCode NEPDestroy_NLEIGS(NEP nep)
2095: {
2097: PetscInt k;
2098: NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
2101: BVDestroy(&ctx->V);
2102: for (k=0;k<ctx->nshiftsw;k++) { KSPDestroy(&ctx->ksp[k]); }
2103: PetscFree(ctx->ksp);
2104: if (ctx->nshifts) { PetscFree(ctx->shifts); }
2105: if (ctx->fullbasis) { EPSDestroy(&ctx->eps); }
2106: PetscFree(nep->data);
2107: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetSingularitiesFunction_C",NULL);
2108: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetSingularitiesFunction_C",NULL);
2109: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRestart_C",NULL);
2110: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRestart_C",NULL);
2111: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetLocking_C",NULL);
2112: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetLocking_C",NULL);
2113: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetInterpolation_C",NULL);
2114: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetInterpolation_C",NULL);
2115: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRKShifts_C",NULL);
2116: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRKShifts_C",NULL);
2117: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetKSPs_C",NULL);
2118: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetFullBasis_C",NULL);
2119: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetFullBasis_C",NULL);
2120: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetEPS_C",NULL);
2121: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetEPS_C",NULL);
2122: return(0);
2123: }
2125: SLEPC_EXTERN PetscErrorCode NEPCreate_NLEIGS(NEP nep)
2126: {
2128: NEP_NLEIGS *ctx;
2131: PetscNewLog(nep,&ctx);
2132: nep->data = (void*)ctx;
2133: ctx->lock = PETSC_TRUE;
2134: ctx->ddtol = PETSC_DEFAULT;
2136: nep->useds = PETSC_TRUE;
2137: nep->hasts = PETSC_TRUE;
2139: nep->ops->setup = NEPSetUp_NLEIGS;
2140: nep->ops->setfromoptions = NEPSetFromOptions_NLEIGS;
2141: nep->ops->view = NEPView_NLEIGS;
2142: nep->ops->destroy = NEPDestroy_NLEIGS;
2143: nep->ops->reset = NEPReset_NLEIGS;
2145: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetSingularitiesFunction_C",NEPNLEIGSSetSingularitiesFunction_NLEIGS);
2146: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetSingularitiesFunction_C",NEPNLEIGSGetSingularitiesFunction_NLEIGS);
2147: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRestart_C",NEPNLEIGSSetRestart_NLEIGS);
2148: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRestart_C",NEPNLEIGSGetRestart_NLEIGS);
2149: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetLocking_C",NEPNLEIGSSetLocking_NLEIGS);
2150: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetLocking_C",NEPNLEIGSGetLocking_NLEIGS);
2151: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetInterpolation_C",NEPNLEIGSSetInterpolation_NLEIGS);
2152: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetInterpolation_C",NEPNLEIGSGetInterpolation_NLEIGS);
2153: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRKShifts_C",NEPNLEIGSSetRKShifts_NLEIGS);
2154: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRKShifts_C",NEPNLEIGSGetRKShifts_NLEIGS);
2155: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetKSPs_C",NEPNLEIGSGetKSPs_NLEIGS);
2156: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetFullBasis_C",NEPNLEIGSSetFullBasis_NLEIGS);
2157: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetFullBasis_C",NEPNLEIGSGetFullBasis_NLEIGS);
2158: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetEPS_C",NEPNLEIGSSetEPS_NLEIGS);
2159: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetEPS_C",NEPNLEIGSGetEPS_NLEIGS);
2160: return(0);
2161: }