Actual source code: nleigs.c
slepc-3.17.2 2022-08-09
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
10: /*
11: SLEPc nonlinear eigensolver: "nleigs"
13: Method: NLEIGS
15: Algorithm:
17: Fully rational Krylov method for nonlinear eigenvalue problems.
19: References:
21: [1] S. Guttel et al., "NLEIGS: A class of robust fully rational Krylov
22: method for nonlinear eigenvalue problems", SIAM J. Sci. Comput.
23: 36(6):A2842-A2864, 2014.
24: */
26: #include <slepc/private/nepimpl.h>
27: #include <slepcblaslapack.h>
28: #include "nleigs.h"
30: PetscErrorCode NEPNLEIGSBackTransform(PetscObject ob,PetscInt n,PetscScalar *valr,PetscScalar *vali)
31: {
32: NEP nep;
33: PetscInt j;
34: #if !defined(PETSC_USE_COMPLEX)
35: PetscScalar t;
36: #endif
38: nep = (NEP)ob;
39: #if !defined(PETSC_USE_COMPLEX)
40: for (j=0;j<n;j++) {
41: if (vali[j] == 0) valr[j] = 1.0 / valr[j] + nep->target;
42: else {
43: t = valr[j] * valr[j] + vali[j] * vali[j];
44: valr[j] = valr[j] / t + nep->target;
45: vali[j] = - vali[j] / t;
46: }
47: }
48: #else
49: for (j=0;j<n;j++) {
50: valr[j] = 1.0 / valr[j] + nep->target;
51: }
52: #endif
53: PetscFunctionReturn(0);
54: }
56: /* Computes the roots of a polynomial */
57: static PetscErrorCode NEPNLEIGSAuxiliarPRootFinder(PetscInt deg,PetscScalar *polcoeffs,PetscScalar *wr,PetscScalar *wi,PetscBool *avail)
58: {
59: PetscScalar *C;
60: PetscBLASInt n_,lwork;
61: PetscInt i;
62: #if defined(PETSC_USE_COMPLEX)
63: PetscReal *rwork=NULL;
64: #endif
65: PetscScalar *work;
66: PetscBLASInt info;
68: *avail = PETSC_TRUE;
69: if (deg>0) {
70: PetscCalloc1(deg*deg,&C);
71: PetscBLASIntCast(deg,&n_);
72: for (i=0;i<deg-1;i++) {
73: C[(deg+1)*i+1] = 1.0;
74: C[(deg-1)*deg+i] = -polcoeffs[deg-i]/polcoeffs[0];
75: }
76: C[deg*deg+-1] = -polcoeffs[1]/polcoeffs[0];
77: PetscBLASIntCast(3*deg,&lwork);
79: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
80: #if !defined(PETSC_USE_COMPLEX)
81: PetscMalloc1(lwork,&work);
82: PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_("N","N",&n_,C,&n_,wr,wi,NULL,&n_,NULL,&n_,work,&lwork,&info));
83: if (info) *avail = PETSC_FALSE;
84: PetscFree(work);
85: #else
86: PetscMalloc2(2*deg,&rwork,lwork,&work);
87: PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_("N","N",&n_,C,&n_,wr,NULL,&n_,NULL,&n_,work,&lwork,rwork,&info));
88: if (info) *avail = PETSC_FALSE;
89: PetscFree2(rwork,work);
90: #endif
91: PetscFPTrapPop();
92: PetscFree(C);
93: }
94: PetscFunctionReturn(0);
95: }
97: static PetscErrorCode NEPNLEIGSAuxiliarRmDuplicates(PetscInt nin,PetscScalar *pin,PetscInt *nout,PetscScalar *pout,PetscInt max)
98: {
99: PetscInt i,j;
101: for (i=0;i<nin;i++) {
102: if (max && *nout>=max) break;
103: pout[(*nout)++] = pin[i];
104: for (j=0;j<*nout-1;j++)
105: if (PetscAbsScalar(pin[i]-pout[j])<PETSC_MACHINE_EPSILON*100) {
106: (*nout)--;
107: break;
108: }
109: }
110: PetscFunctionReturn(0);
111: }
113: static PetscErrorCode NEPNLEIGSFNSingularities(FN f,PetscInt *nisol,PetscScalar **isol,PetscBool *rational)
114: {
115: FNCombineType ctype;
116: FN f1,f2;
117: PetscInt i,nq,nisol1,nisol2;
118: PetscScalar *qcoeff,*wr,*wi,*isol1,*isol2;
119: PetscBool flg,avail,rat1,rat2;
121: *rational = PETSC_FALSE;
122: PetscObjectTypeCompare((PetscObject)f,FNRATIONAL,&flg);
123: if (flg) {
124: *rational = PETSC_TRUE;
125: FNRationalGetDenominator(f,&nq,&qcoeff);
126: if (nq>1) {
127: PetscMalloc2(nq-1,&wr,nq-1,&wi);
128: NEPNLEIGSAuxiliarPRootFinder(nq-1,qcoeff,wr,wi,&avail);
129: if (avail) {
130: PetscCalloc1(nq-1,isol);
131: *nisol = 0;
132: for (i=0;i<nq-1;i++)
133: #if !defined(PETSC_USE_COMPLEX)
134: if (wi[i]==0)
135: #endif
136: (*isol)[(*nisol)++] = wr[i];
137: nq = *nisol; *nisol = 0;
138: for (i=0;i<nq;i++) wr[i] = (*isol)[i];
139: NEPNLEIGSAuxiliarRmDuplicates(nq,wr,nisol,*isol,0);
140: PetscFree2(wr,wi);
141: } else { *nisol=0; *isol = NULL; }
142: } else { *nisol = 0; *isol = NULL; }
143: PetscFree(qcoeff);
144: }
145: PetscObjectTypeCompare((PetscObject)f,FNCOMBINE,&flg);
146: if (flg) {
147: FNCombineGetChildren(f,&ctype,&f1,&f2);
148: if (ctype != FN_COMBINE_COMPOSE && ctype != FN_COMBINE_DIVIDE) {
149: NEPNLEIGSFNSingularities(f1,&nisol1,&isol1,&rat1);
150: NEPNLEIGSFNSingularities(f2,&nisol2,&isol2,&rat2);
151: if (nisol1+nisol2>0) {
152: PetscCalloc1(nisol1+nisol2,isol);
153: *nisol = 0;
154: NEPNLEIGSAuxiliarRmDuplicates(nisol1,isol1,nisol,*isol,0);
155: NEPNLEIGSAuxiliarRmDuplicates(nisol2,isol2,nisol,*isol,0);
156: }
157: *rational = (rat1&&rat2)?PETSC_TRUE:PETSC_FALSE;
158: PetscFree(isol1);
159: PetscFree(isol2);
160: }
161: }
162: PetscFunctionReturn(0);
163: }
165: static PetscErrorCode NEPNLEIGSRationalSingularities(NEP nep,PetscInt *ndptx,PetscScalar *dxi,PetscBool *rational)
166: {
167: PetscInt nt,i,nisol;
168: FN f;
169: PetscScalar *isol;
170: PetscBool rat;
172: *rational = PETSC_TRUE;
173: *ndptx = 0;
174: NEPGetSplitOperatorInfo(nep,&nt,NULL);
175: for (i=0;i<nt;i++) {
176: NEPGetSplitOperatorTerm(nep,i,NULL,&f);
177: NEPNLEIGSFNSingularities(f,&nisol,&isol,&rat);
178: if (nisol) {
179: NEPNLEIGSAuxiliarRmDuplicates(nisol,isol,ndptx,dxi,0);
180: PetscFree(isol);
181: }
182: *rational = ((*rational)&&rat)?PETSC_TRUE:PETSC_FALSE;
183: }
184: PetscFunctionReturn(0);
185: }
187: /* Adaptive Anderson-Antoulas algorithm */
188: static PetscErrorCode NEPNLEIGSAAAComputation(NEP nep,PetscInt ndpt,PetscScalar *ds,PetscScalar *F,PetscInt *ndptx,PetscScalar *dxi)
189: {
190: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
191: PetscScalar mean=0.0,*z,*f,*C,*A,*VT,*work,*ww,szero=0.0,sone=1.0;
192: PetscScalar *N,*D;
193: PetscReal *S,norm,err,*R;
194: PetscInt i,k,j,idx=0,cont;
195: PetscBLASInt n_,m_,lda_,lwork,info,one=1;
196: #if defined(PETSC_USE_COMPLEX)
197: PetscReal *rwork;
198: #endif
200: PetscBLASIntCast(8*ndpt,&lwork);
201: PetscMalloc5(ndpt,&R,ndpt,&z,ndpt,&f,ndpt*ndpt,&C,ndpt,&ww);
202: PetscMalloc6(ndpt*ndpt,&A,ndpt,&S,ndpt*ndpt,&VT,lwork,&work,ndpt,&D,ndpt,&N);
203: #if defined(PETSC_USE_COMPLEX)
204: PetscMalloc1(8*ndpt,&rwork);
205: #endif
206: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
207: norm = 0.0;
208: for (i=0;i<ndpt;i++) {
209: mean += F[i];
210: norm = PetscMax(PetscAbsScalar(F[i]),norm);
211: }
212: mean /= ndpt;
213: PetscBLASIntCast(ndpt,&lda_);
214: for (i=0;i<ndpt;i++) R[i] = PetscAbsScalar(F[i]-mean);
215: /* next support point */
216: err = 0.0;
217: for (i=0;i<ndpt;i++) if (R[i]>=err) {idx = i; err = R[i];}
218: for (k=0;k<ndpt-1;k++) {
219: z[k] = ds[idx]; f[k] = F[idx]; R[idx] = -1.0;
220: /* next column of Cauchy matrix */
221: for (i=0;i<ndpt;i++) {
222: C[i+k*ndpt] = 1.0/(ds[i]-ds[idx]);
223: }
225: PetscArrayzero(A,ndpt*ndpt);
226: cont = 0;
227: for (i=0;i<ndpt;i++) {
228: if (R[i]!=-1.0) {
229: for (j=0;j<=k;j++)A[cont+j*ndpt] = C[i+j*ndpt]*F[i]-C[i+j*ndpt]*f[j];
230: cont++;
231: }
232: }
233: PetscBLASIntCast(cont,&m_);
234: PetscBLASIntCast(k+1,&n_);
235: #if defined(PETSC_USE_COMPLEX)
236: PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","A",&m_,&n_,A,&lda_,S,NULL,&lda_,VT,&lda_,work,&lwork,rwork,&info));
237: #else
238: PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","A",&m_,&n_,A,&lda_,S,NULL,&lda_,VT,&lda_,work,&lwork,&info));
239: #endif
240: SlepcCheckLapackInfo("gesvd",info);
241: for (i=0;i<=k;i++) {
242: ww[i] = PetscConj(VT[i*ndpt+k]);
243: D[i] = ww[i]*f[i];
244: }
245: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&lda_,&n_,&sone,C,&lda_,D,&one,&szero,N,&one));
246: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&lda_,&n_,&sone,C,&lda_,ww,&one,&szero,D,&one));
247: for (i=0;i<ndpt;i++) if (R[i]>=0) R[i] = PetscAbsScalar(F[i]-N[i]/D[i]);
248: /* next support point */
249: err = 0.0;
250: for (i=0;i<ndpt;i++) if (R[i]>=err) {idx = i; err = R[i];}
251: if (err <= ctx->ddtol*norm) break;
252: }
255: /* poles */
256: PetscArrayzero(C,ndpt*ndpt);
257: PetscArrayzero(A,ndpt*ndpt);
258: for (i=0;i<=k;i++) {
259: C[i+ndpt*i] = 1.0;
260: A[(i+1)*ndpt] = ww[i];
261: A[i+1] = 1.0;
262: A[i+1+(i+1)*ndpt] = z[i];
263: }
264: C[0] = 0.0; C[k+1+(k+1)*ndpt] = 1.0;
265: n_++;
266: #if defined(PETSC_USE_COMPLEX)
267: PetscStackCallBLAS("LAPACKggev",LAPACKggev_("N","N",&n_,A,&lda_,C,&lda_,D,N,NULL,&lda_,NULL,&lda_,work,&lwork,rwork,&info));
268: #else
269: PetscStackCallBLAS("LAPACKggev",LAPACKggev_("N","N",&n_,A,&lda_,C,&lda_,D,VT,N,NULL,&lda_,NULL,&lda_,work,&lwork,&info));
270: #endif
271: SlepcCheckLapackInfo("ggev",info);
272: cont = 0.0;
273: for (i=0;i<n_;i++) if (N[i]!=0.0) {
274: dxi[cont++] = D[i]/N[i];
275: }
276: *ndptx = cont;
277: PetscFPTrapPop();
278: PetscFree5(R,z,f,C,ww);
279: PetscFree6(A,S,VT,work,D,N);
280: #if defined(PETSC_USE_COMPLEX)
281: PetscFree(rwork);
282: #endif
283: PetscFunctionReturn(0);
284: }
286: /* Singularities using Adaptive Anderson-Antoulas algorithm */
287: static PetscErrorCode NEPNLEIGSAAASingularities(NEP nep,PetscInt ndpt,PetscScalar *ds,PetscInt *ndptx,PetscScalar *dxi)
288: {
289: Vec u,v,w;
290: PetscRandom rand=NULL;
291: PetscScalar *F,*isol;
292: PetscInt i,k,nisol,nt;
293: Mat T;
294: FN f;
296: PetscMalloc1(ndpt,&F);
297: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
298: PetscMalloc1(ndpt,&isol);
299: *ndptx = 0;
300: NEPGetSplitOperatorInfo(nep,&nt,NULL);
301: nisol = *ndptx;
302: for (k=0;k<nt;k++) {
303: NEPGetSplitOperatorTerm(nep,k,NULL,&f);
304: for (i=0;i<ndpt;i++) FNEvaluateFunction(f,ds[i],&F[i]);
305: NEPNLEIGSAAAComputation(nep,ndpt,ds,F,&nisol,isol);
306: if (nisol) NEPNLEIGSAuxiliarRmDuplicates(nisol,isol,ndptx,dxi,ndpt);
307: }
308: PetscFree(isol);
309: } else {
310: MatCreateVecs(nep->function,&u,NULL);
311: VecDuplicate(u,&v);
312: VecDuplicate(u,&w);
313: if (nep->V) BVGetRandomContext(nep->V,&rand);
314: VecSetRandom(u,rand);
315: VecNormalize(u,NULL);
316: VecSetRandom(v,rand);
317: VecNormalize(v,NULL);
318: T = nep->function;
319: for (i=0;i<ndpt;i++) {
320: NEPComputeFunction(nep,ds[i],T,T);
321: MatMult(T,v,w);
322: VecDot(w,u,&F[i]);
323: }
324: NEPNLEIGSAAAComputation(nep,ndpt,ds,F,ndptx,dxi);
325: VecDestroy(&u);
326: VecDestroy(&v);
327: VecDestroy(&w);
328: }
329: PetscFree(F);
330: PetscFunctionReturn(0);
331: }
333: static PetscErrorCode NEPNLEIGSLejaBagbyPoints(NEP nep)
334: {
335: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
336: PetscInt i,k,ndpt=NDPOINTS,ndptx=NDPOINTS;
337: PetscScalar *ds,*dsi,*dxi,*nrs,*nrxi,*s=ctx->s,*xi=ctx->xi,*beta=ctx->beta;
338: PetscReal maxnrs,minnrxi;
339: PetscBool rational;
340: #if !defined(PETSC_USE_COMPLEX)
341: PetscReal a,b,h;
342: #endif
344: if (!ctx->computesingularities && nep->problem_type!=NEP_RATIONAL) ndpt = ndptx = LBPOINTS;
345: PetscMalloc5(ndpt+1,&ds,ndpt+1,&dsi,ndpt,&dxi,ndpt+1,&nrs,ndpt,&nrxi);
347: /* Discretize the target region boundary */
348: RGComputeContour(nep->rg,ndpt,ds,dsi);
349: #if !defined(PETSC_USE_COMPLEX)
350: for (i=0;i<ndpt;i++) if (dsi[i]!=0.0) break;
351: if (i<ndpt) {
353: /* Select a segment in the real axis */
354: RGComputeBoundingBox(nep->rg,&a,&b,NULL,NULL);
356: h = (b-a)/ndpt;
357: for (i=0;i<ndpt;i++) {ds[i] = a+h*i; dsi[i] = 0.0;}
358: }
359: #endif
360: /* Discretize the singularity region */
361: if (ctx->computesingularities) (ctx->computesingularities)(nep,&ndptx,dxi,ctx->singularitiesctx);
362: else {
363: if (nep->problem_type==NEP_RATIONAL) {
364: NEPNLEIGSRationalSingularities(nep,&ndptx,dxi,&rational);
366: } else {
367: /* AAA algorithm */
368: NEPNLEIGSAAASingularities(nep,ndpt,ds,&ndptx,dxi);
369: }
370: }
371: /* Look for Leja-Bagby points in the discretization sets */
372: s[0] = ds[0];
373: xi[0] = (ndptx>0)?dxi[0]:PETSC_INFINITY;
375: beta[0] = 1.0; /* scaling factors are also computed here */
376: for (i=0;i<ndpt;i++) {
377: nrs[i] = 1.0;
378: nrxi[i] = 1.0;
379: }
380: for (k=1;k<ctx->ddmaxit;k++) {
381: maxnrs = 0.0;
382: minnrxi = PETSC_MAX_REAL;
383: for (i=0;i<ndpt;i++) {
384: nrs[i] *= ((ds[i]-s[k-1])/(1.0-ds[i]/xi[k-1]))/beta[k-1];
385: if (PetscAbsScalar(nrs[i])>maxnrs) {maxnrs = PetscAbsScalar(nrs[i]); s[k] = ds[i];}
386: }
387: if (ndptx>k) {
388: for (i=1;i<ndptx;i++) {
389: nrxi[i] *= ((dxi[i]-s[k-1])/(1.0-dxi[i]/xi[k-1]))/beta[k-1];
390: if (PetscAbsScalar(nrxi[i])<minnrxi) {minnrxi = PetscAbsScalar(nrxi[i]); xi[k] = dxi[i];}
391: }
393: } else xi[k] = PETSC_INFINITY;
394: beta[k] = maxnrs;
395: }
396: PetscFree5(ds,dsi,dxi,nrs,nrxi);
397: PetscFunctionReturn(0);
398: }
400: PetscErrorCode NEPNLEIGSEvalNRTFunct(NEP nep,PetscInt k,PetscScalar sigma,PetscScalar *b)
401: {
402: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
403: PetscInt i;
404: PetscScalar *beta=ctx->beta,*s=ctx->s,*xi=ctx->xi;
406: b[0] = 1.0/beta[0];
407: for (i=0;i<k;i++) {
408: b[i+1] = ((sigma-s[i])*b[i])/(beta[i+1]*(1.0-sigma/xi[i]));
409: }
410: PetscFunctionReturn(0);
411: }
413: static PetscErrorCode MatMult_Fun(Mat A,Vec x,Vec y)
414: {
415: NEP_NLEIGS_MATSHELL *ctx;
416: PetscInt i;
419: MatShellGetContext(A,&ctx);
420: MatMult(ctx->A[0],x,y);
421: if (ctx->coeff[0]!=1.0) VecScale(y,ctx->coeff[0]);
422: for (i=1;i<ctx->nmat;i++) {
423: MatMult(ctx->A[i],x,ctx->t);
424: VecAXPY(y,ctx->coeff[i],ctx->t);
425: }
426: PetscFunctionReturn(0);
427: }
429: static PetscErrorCode MatMultTranspose_Fun(Mat A,Vec x,Vec y)
430: {
431: NEP_NLEIGS_MATSHELL *ctx;
432: PetscInt i;
435: MatShellGetContext(A,&ctx);
436: MatMultTranspose(ctx->A[0],x,y);
437: if (ctx->coeff[0]!=1.0) VecScale(y,ctx->coeff[0]);
438: for (i=1;i<ctx->nmat;i++) {
439: MatMultTranspose(ctx->A[i],x,ctx->t);
440: VecAXPY(y,ctx->coeff[i],ctx->t);
441: }
442: PetscFunctionReturn(0);
443: }
445: static PetscErrorCode MatGetDiagonal_Fun(Mat A,Vec diag)
446: {
447: NEP_NLEIGS_MATSHELL *ctx;
448: PetscInt i;
451: MatShellGetContext(A,&ctx);
452: MatGetDiagonal(ctx->A[0],diag);
453: if (ctx->coeff[0]!=1.0) VecScale(diag,ctx->coeff[0]);
454: for (i=1;i<ctx->nmat;i++) {
455: MatGetDiagonal(ctx->A[i],ctx->t);
456: VecAXPY(diag,ctx->coeff[i],ctx->t);
457: }
458: PetscFunctionReturn(0);
459: }
461: static PetscErrorCode MatDuplicate_Fun(Mat A,MatDuplicateOption op,Mat *B)
462: {
463: PetscInt m,n,M,N,i;
464: NEP_NLEIGS_MATSHELL *ctxnew,*ctx;
465: void (*fun)(void);
468: MatShellGetContext(A,&ctx);
469: PetscNew(&ctxnew);
470: ctxnew->nmat = ctx->nmat;
471: ctxnew->maxnmat = ctx->maxnmat;
472: PetscMalloc2(ctxnew->maxnmat,&ctxnew->A,ctxnew->maxnmat,&ctxnew->coeff);
473: for (i=0;i<ctx->nmat;i++) {
474: PetscObjectReference((PetscObject)ctx->A[i]);
475: ctxnew->A[i] = ctx->A[i];
476: ctxnew->coeff[i] = ctx->coeff[i];
477: }
478: MatGetSize(ctx->A[0],&M,&N);
479: MatGetLocalSize(ctx->A[0],&m,&n);
480: VecDuplicate(ctx->t,&ctxnew->t);
481: MatCreateShell(PetscObjectComm((PetscObject)A),m,n,M,N,(void*)ctxnew,B);
482: MatShellSetManageScalingShifts(*B);
483: MatShellGetOperation(A,MATOP_MULT,&fun);
484: MatShellSetOperation(*B,MATOP_MULT,fun);
485: MatShellGetOperation(A,MATOP_MULT_TRANSPOSE,&fun);
486: MatShellSetOperation(*B,MATOP_MULT_TRANSPOSE,fun);
487: MatShellGetOperation(A,MATOP_GET_DIAGONAL,&fun);
488: MatShellSetOperation(*B,MATOP_GET_DIAGONAL,fun);
489: MatShellGetOperation(A,MATOP_DUPLICATE,&fun);
490: MatShellSetOperation(*B,MATOP_DUPLICATE,fun);
491: MatShellGetOperation(A,MATOP_DESTROY,&fun);
492: MatShellSetOperation(*B,MATOP_DESTROY,fun);
493: MatShellGetOperation(A,MATOP_AXPY,&fun);
494: MatShellSetOperation(*B,MATOP_AXPY,fun);
495: PetscFunctionReturn(0);
496: }
498: static PetscErrorCode MatDestroy_Fun(Mat A)
499: {
500: NEP_NLEIGS_MATSHELL *ctx;
501: PetscInt i;
504: if (A) {
505: MatShellGetContext(A,&ctx);
506: for (i=0;i<ctx->nmat;i++) MatDestroy(&ctx->A[i]);
507: VecDestroy(&ctx->t);
508: PetscFree2(ctx->A,ctx->coeff);
509: PetscFree(ctx);
510: }
511: PetscFunctionReturn(0);
512: }
514: static PetscErrorCode MatAXPY_Fun(Mat Y,PetscScalar a,Mat X,MatStructure str)
515: {
516: NEP_NLEIGS_MATSHELL *ctxY,*ctxX;
517: PetscInt i,j;
518: PetscBool found;
521: MatShellGetContext(Y,&ctxY);
522: MatShellGetContext(X,&ctxX);
523: for (i=0;i<ctxX->nmat;i++) {
524: found = PETSC_FALSE;
525: for (j=0;!found&&j<ctxY->nmat;j++) {
526: if (ctxX->A[i]==ctxY->A[j]) {
527: found = PETSC_TRUE;
528: ctxY->coeff[j] += a*ctxX->coeff[i];
529: }
530: }
531: if (!found) {
532: ctxY->coeff[ctxY->nmat] = a*ctxX->coeff[i];
533: ctxY->A[ctxY->nmat++] = ctxX->A[i];
534: PetscObjectReference((PetscObject)ctxX->A[i]);
535: }
536: }
537: PetscFunctionReturn(0);
538: }
540: static PetscErrorCode MatScale_Fun(Mat M,PetscScalar a)
541: {
542: NEP_NLEIGS_MATSHELL *ctx;
543: PetscInt i;
546: MatShellGetContext(M,&ctx);
547: for (i=0;i<ctx->nmat;i++) ctx->coeff[i] *= a;
548: PetscFunctionReturn(0);
549: }
551: static PetscErrorCode NLEIGSMatToMatShellArray(Mat A,Mat *Ms,PetscInt maxnmat)
552: {
553: NEP_NLEIGS_MATSHELL *ctx;
554: PetscInt m,n,M,N;
555: PetscBool has;
557: MatHasOperation(A,MATOP_DUPLICATE,&has);
559: PetscNew(&ctx);
560: ctx->maxnmat = maxnmat;
561: PetscMalloc2(ctx->maxnmat,&ctx->A,ctx->maxnmat,&ctx->coeff);
562: MatDuplicate(A,MAT_COPY_VALUES,&ctx->A[0]);
563: ctx->nmat = 1;
564: ctx->coeff[0] = 1.0;
565: MatCreateVecs(A,&ctx->t,NULL);
566: MatGetSize(A,&M,&N);
567: MatGetLocalSize(A,&m,&n);
568: MatCreateShell(PetscObjectComm((PetscObject)A),m,n,M,N,(void*)ctx,Ms);
569: MatShellSetManageScalingShifts(*Ms);
570: MatShellSetOperation(*Ms,MATOP_MULT,(void(*)(void))MatMult_Fun);
571: MatShellSetOperation(*Ms,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_Fun);
572: MatShellSetOperation(*Ms,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Fun);
573: MatShellSetOperation(*Ms,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_Fun);
574: MatShellSetOperation(*Ms,MATOP_DESTROY,(void(*)(void))MatDestroy_Fun);
575: MatShellSetOperation(*Ms,MATOP_AXPY,(void(*)(void))MatAXPY_Fun);
576: MatShellSetOperation(*Ms,MATOP_SCALE,(void(*)(void))MatScale_Fun);
577: PetscFunctionReturn(0);
578: }
580: /*
581: MatIsShellAny - returns true if any of the n matrices is a shell matrix
582: */
583: static PetscErrorCode MatIsShellAny(Mat *A,PetscInt n,PetscBool *shell)
584: {
585: PetscInt i;
586: PetscBool flg;
588: *shell = PETSC_FALSE;
589: for (i=0;i<n;i++) {
590: MatIsShell(A[i],&flg);
591: if (flg) { *shell = PETSC_TRUE; break; }
592: }
593: PetscFunctionReturn(0);
594: }
596: static PetscErrorCode NEPNLEIGSDividedDifferences_split(NEP nep)
597: {
599: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
600: PetscInt k,j,i,maxnmat,nmax;
601: PetscReal norm0,norm,*matnorm;
602: PetscScalar *s=ctx->s,*beta=ctx->beta,*xi=ctx->xi,*b,alpha,*coeffs,*pK,*pH,sone=1.0;
603: Mat T,P,Ts,K,H;
604: PetscBool shell,hasmnorm=PETSC_FALSE,matrix=PETSC_TRUE;
605: PetscBLASInt n_;
607: nmax = ctx->ddmaxit;
608: PetscMalloc1(nep->nt*nmax,&ctx->coeffD);
609: PetscMalloc3(nmax+1,&b,nmax+1,&coeffs,nep->nt,&matnorm);
610: for (j=0;j<nep->nt;j++) {
611: MatHasOperation(nep->A[j],MATOP_NORM,&hasmnorm);
612: if (!hasmnorm) break;
613: MatNorm(nep->A[j],NORM_INFINITY,matnorm+j);
614: }
615: /* Try matrix functions scheme */
616: PetscCalloc2(nmax*nmax,&pK,nmax*nmax,&pH);
617: for (i=0;i<nmax-1;i++) {
618: pK[(nmax+1)*i] = 1.0;
619: pK[(nmax+1)*i+1] = beta[i+1]/xi[i];
620: pH[(nmax+1)*i] = s[i];
621: pH[(nmax+1)*i+1] = beta[i+1];
622: }
623: pH[nmax*nmax-1] = s[nmax-1];
624: pK[nmax*nmax-1] = 1.0;
625: PetscBLASIntCast(nmax,&n_);
626: PetscStackCallBLAS("BLAStrsm",BLAStrsm_("R","L","N","U",&n_,&n_,&sone,pK,&n_,pH,&n_));
627: /* The matrix to be used is in H. K will be a work-space matrix */
628: MatCreateSeqDense(PETSC_COMM_SELF,nmax,nmax,pH,&H);
629: MatCreateSeqDense(PETSC_COMM_SELF,nmax,nmax,pK,&K);
630: for (j=0;matrix&&j<nep->nt;j++) {
631: PetscPushErrorHandler(PetscReturnErrorHandler,NULL);
632: ierr = FNEvaluateFunctionMat(nep->f[j],H,K);
633: PetscPopErrorHandler();
634: if (!ierr) {
635: for (i=0;i<nmax;i++) ctx->coeffD[j+i*nep->nt] = pK[i]*beta[0];
636: } else {
637: matrix = PETSC_FALSE;
638: PetscFPTrapPop();
639: }
640: }
641: MatDestroy(&H);
642: MatDestroy(&K);
643: if (!matrix) {
644: for (j=0;j<nep->nt;j++) {
645: FNEvaluateFunction(nep->f[j],s[0],ctx->coeffD+j);
646: ctx->coeffD[j] *= beta[0];
647: }
648: }
649: if (hasmnorm) {
650: norm0 = 0.0;
651: for (j=0;j<nep->nt;j++) norm0 += matnorm[j]*PetscAbsScalar(ctx->coeffD[j]);
652: } else {
653: norm0 = 0.0;
654: for (j=0;j<nep->nt;j++) norm0 = PetscMax(PetscAbsScalar(ctx->coeffD[j]),norm0);
655: }
656: ctx->nmat = ctx->ddmaxit;
657: for (k=1;k<ctx->ddmaxit;k++) {
658: if (!matrix) {
659: NEPNLEIGSEvalNRTFunct(nep,k,s[k],b);
660: for (i=0;i<nep->nt;i++) {
661: FNEvaluateFunction(nep->f[i],s[k],ctx->coeffD+k*nep->nt+i);
662: for (j=0;j<k;j++) {
663: ctx->coeffD[k*nep->nt+i] -= b[j]*ctx->coeffD[i+nep->nt*j];
664: }
665: ctx->coeffD[k*nep->nt+i] /= b[k];
666: }
667: }
668: if (hasmnorm) {
669: norm = 0.0;
670: for (j=0;j<nep->nt;j++) norm += matnorm[j]*PetscAbsScalar(ctx->coeffD[k*nep->nt+j]);
671: } else {
672: norm = 0.0;
673: for (j=0;j<nep->nt;j++) norm = PetscMax(PetscAbsScalar(ctx->coeffD[k*nep->nt+j]),norm);
674: }
675: if (k>1 && norm/norm0 < ctx->ddtol) {
676: ctx->nmat = k+1;
677: break;
678: }
679: }
680: if (!ctx->ksp) NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp);
681: MatIsShellAny(nep->A,nep->nt,&shell);
682: maxnmat = PetscMax(ctx->ddmaxit,nep->nt);
683: for (i=0;i<ctx->nshiftsw;i++) {
684: NEPNLEIGSEvalNRTFunct(nep,ctx->nmat-1,ctx->shifts[i],coeffs);
685: if (!shell) MatDuplicate(nep->A[0],MAT_COPY_VALUES,&T);
686: else NLEIGSMatToMatShellArray(nep->A[0],&T,maxnmat);
687: if (nep->P) { /* user-defined preconditioner */
688: MatDuplicate(nep->P[0],MAT_COPY_VALUES,&P);
689: } else P=T;
690: alpha = 0.0;
691: for (j=0;j<ctx->nmat;j++) alpha += coeffs[j]*ctx->coeffD[j*nep->nt];
692: MatScale(T,alpha);
693: if (nep->P) MatScale(P,alpha);
694: for (k=1;k<nep->nt;k++) {
695: alpha = 0.0;
696: for (j=0;j<ctx->nmat;j++) alpha += coeffs[j]*ctx->coeffD[j*nep->nt+k];
697: if (shell) NLEIGSMatToMatShellArray(nep->A[k],&Ts,maxnmat);
698: MatAXPY(T,alpha,shell?Ts:nep->A[k],nep->mstr);
699: if (nep->P) MatAXPY(P,alpha,nep->P[k],nep->mstrp);
700: if (shell) MatDestroy(&Ts);
701: }
702: NEP_KSPSetOperators(ctx->ksp[i],T,P);
703: KSPSetUp(ctx->ksp[i]);
704: MatDestroy(&T);
705: if (nep->P) MatDestroy(&P);
706: }
707: PetscFree3(b,coeffs,matnorm);
708: PetscFree2(pK,pH);
709: PetscFunctionReturn(0);
710: }
712: static PetscErrorCode NEPNLEIGSDividedDifferences_callback(NEP nep)
713: {
714: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
715: PetscInt k,j,i,maxnmat;
716: PetscReal norm0,norm;
717: PetscScalar *s=ctx->s,*beta=ctx->beta,*b,*coeffs;
718: Mat *D=ctx->D,*DP,T,P;
719: PetscBool shell,has,vec=PETSC_FALSE,precond=(nep->function_pre!=nep->function)?PETSC_TRUE:PETSC_FALSE;
720: PetscRandom rand=NULL;
721: Vec w[2];
723: PetscMalloc2(ctx->ddmaxit+1,&b,ctx->ddmaxit+1,&coeffs);
724: if (nep->V) BVGetRandomContext(nep->V,&rand);
725: T = nep->function;
726: P = nep->function_pre;
727: NEPComputeFunction(nep,s[0],T,P);
728: MatIsShell(T,&shell);
729: maxnmat = PetscMax(ctx->ddmaxit,nep->nt);
730: if (!shell) MatDuplicate(T,MAT_COPY_VALUES,&D[0]);
731: else NLEIGSMatToMatShellArray(T,&D[0],maxnmat);
732: if (beta[0]!=1.0) MatScale(D[0],1.0/beta[0]);
733: MatHasOperation(D[0],MATOP_NORM,&has);
734: if (has) MatNorm(D[0],NORM_FROBENIUS,&norm0);
735: else {
736: MatCreateVecs(D[0],NULL,&w[0]);
737: VecDuplicate(w[0],&w[1]);
738: VecDuplicate(w[0],&ctx->vrn);
739: VecSetRandomNormal(ctx->vrn,rand,w[0],w[1]);
740: VecNormalize(ctx->vrn,NULL);
741: vec = PETSC_TRUE;
742: MatNormEstimate(D[0],ctx->vrn,w[0],&norm0);
743: }
744: if (precond) {
745: PetscMalloc1(ctx->ddmaxit,&DP);
746: MatDuplicate(P,MAT_COPY_VALUES,&DP[0]);
747: }
748: ctx->nmat = ctx->ddmaxit;
749: for (k=1;k<ctx->ddmaxit;k++) {
750: NEPNLEIGSEvalNRTFunct(nep,k,s[k],b);
751: NEPComputeFunction(nep,s[k],T,P);
752: if (!shell) MatDuplicate(T,MAT_COPY_VALUES,&D[k]);
753: else NLEIGSMatToMatShellArray(T,&D[k],maxnmat);
754: for (j=0;j<k;j++) MatAXPY(D[k],-b[j],D[j],nep->mstr);
755: MatScale(D[k],1.0/b[k]);
756: MatHasOperation(D[k],MATOP_NORM,&has);
757: if (has) MatNorm(D[k],NORM_FROBENIUS,&norm);
758: else {
759: if (!vec) {
760: MatCreateVecs(D[k],NULL,&w[0]);
761: VecDuplicate(w[0],&w[1]);
762: VecDuplicate(w[0],&ctx->vrn);
763: VecSetRandomNormal(ctx->vrn,rand,w[0],w[1]);
764: VecNormalize(ctx->vrn,NULL);
765: vec = PETSC_TRUE;
766: }
767: MatNormEstimate(D[k],ctx->vrn,w[0],&norm);
768: }
769: if (precond) {
770: MatDuplicate(P,MAT_COPY_VALUES,&DP[k]);
771: for (j=0;j<k;j++) MatAXPY(DP[k],-b[j],DP[j],nep->mstrp);
772: MatScale(DP[k],1.0/b[k]);
773: }
774: if (k>1 && norm/norm0 < ctx->ddtol && k>1) {
775: ctx->nmat = k+1;
776: break;
777: }
778: }
779: if (!ctx->ksp) NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp);
780: for (i=0;i<ctx->nshiftsw;i++) {
781: NEPNLEIGSEvalNRTFunct(nep,ctx->nmat-1,ctx->shifts[i],coeffs);
782: MatDuplicate(D[0],MAT_COPY_VALUES,&T);
783: if (coeffs[0]!=1.0) MatScale(T,coeffs[0]);
784: for (j=1;j<ctx->nmat;j++) MatAXPY(T,coeffs[j],D[j],nep->mstr);
785: if (precond) {
786: MatDuplicate(DP[0],MAT_COPY_VALUES,&P);
787: if (coeffs[0]!=1.0) MatScale(P,coeffs[0]);
788: for (j=1;j<ctx->nmat;j++) MatAXPY(P,coeffs[j],DP[j],nep->mstrp);
789: } else P=T;
790: NEP_KSPSetOperators(ctx->ksp[i],T,P);
791: KSPSetUp(ctx->ksp[i]);
792: MatDestroy(&T);
793: }
794: PetscFree2(b,coeffs);
795: if (vec) {
796: VecDestroy(&w[0]);
797: VecDestroy(&w[1]);
798: }
799: if (precond) {
800: MatDestroy(&P);
801: MatDestroyMatrices(ctx->nmat,&DP);
802: }
803: PetscFunctionReturn(0);
804: }
806: /*
807: NEPKrylovConvergence - This is the analogue to EPSKrylovConvergence.
808: */
809: PetscErrorCode NEPNLEIGSKrylovConvergence(NEP nep,PetscBool getall,PetscInt kini,PetscInt nits,PetscReal betah,PetscScalar betak,PetscInt *kout,Vec *w)
810: {
811: PetscInt k,newk,marker,inside;
812: PetscScalar re,im;
813: PetscReal resnorm,tt;
814: PetscBool istrivial;
815: NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
817: RGIsTrivial(nep->rg,&istrivial);
818: marker = -1;
819: if (nep->trackall) getall = PETSC_TRUE;
820: for (k=kini;k<kini+nits;k++) {
821: /* eigenvalue */
822: re = nep->eigr[k];
823: im = nep->eigi[k];
824: if (!istrivial) {
825: if (!ctx->nshifts) NEPNLEIGSBackTransform((PetscObject)nep,1,&re,&im);
826: RGCheckInside(nep->rg,1,&re,&im,&inside);
827: if (marker==-1 && inside<0) marker = k;
828: }
829: newk = k;
830: DSVectors(nep->ds,DS_MAT_X,&newk,&resnorm);
831: tt = (ctx->nshifts)?SlepcAbsEigenvalue(betak-nep->eigr[k]*betah,nep->eigi[k]*betah):betah;
832: resnorm *= PetscAbsReal(tt);
833: /* error estimate */
834: (*nep->converged)(nep,nep->eigr[k],nep->eigi[k],resnorm,&nep->errest[k],nep->convergedctx);
835: if (marker==-1 && nep->errest[k] >= nep->tol) marker = k;
836: if (newk==k+1) {
837: nep->errest[k+1] = nep->errest[k];
838: k++;
839: }
840: if (marker!=-1 && !getall) break;
841: }
842: if (marker!=-1) k = marker;
843: *kout = k;
844: PetscFunctionReturn(0);
845: }
847: PetscErrorCode NEPSetUp_NLEIGS(NEP nep)
848: {
849: PetscInt k,in;
850: PetscScalar zero=0.0;
851: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
852: SlepcSC sc;
853: PetscBool istrivial;
855: NEPSetDimensions_Default(nep,nep->nev,&nep->ncv,&nep->mpd);
857: if (nep->max_it==PETSC_DEFAULT) nep->max_it = PetscMax(5000,2*nep->n/nep->ncv);
858: if (!ctx->ddmaxit) ctx->ddmaxit = LBPOINTS;
859: RGIsTrivial(nep->rg,&istrivial);
861: if (!nep->which) nep->which = NEP_TARGET_MAGNITUDE;
864: /* Initialize the NLEIGS context structure */
865: k = ctx->ddmaxit;
866: PetscMalloc4(k,&ctx->s,k,&ctx->xi,k,&ctx->beta,k,&ctx->D);
867: nep->data = ctx;
868: if (nep->tol==PETSC_DEFAULT) nep->tol = SLEPC_DEFAULT_TOL;
869: if (ctx->ddtol==PETSC_DEFAULT) ctx->ddtol = nep->tol/10.0;
870: if (!ctx->keep) ctx->keep = 0.5;
872: /* Compute Leja-Bagby points and scaling values */
873: NEPNLEIGSLejaBagbyPoints(nep);
874: if (nep->problem_type!=NEP_RATIONAL) {
875: RGCheckInside(nep->rg,1,&nep->target,&zero,&in);
877: }
879: /* Compute the divided difference matrices */
880: if (nep->fui==NEP_USER_INTERFACE_SPLIT) NEPNLEIGSDividedDifferences_split(nep);
881: else NEPNLEIGSDividedDifferences_callback(nep);
882: NEPAllocateSolution(nep,ctx->nmat-1);
883: NEPSetWorkVecs(nep,4);
884: if (!ctx->fullbasis) {
886: /* set-up DS and transfer split operator functions */
887: DSSetType(nep->ds,ctx->nshifts?DSGNHEP:DSNHEP);
888: DSAllocate(nep->ds,nep->ncv+1);
889: DSGetSlepcSC(nep->ds,&sc);
890: if (!ctx->nshifts) sc->map = NEPNLEIGSBackTransform;
891: DSSetExtraRow(nep->ds,PETSC_TRUE);
892: sc->mapobj = (PetscObject)nep;
893: sc->rg = nep->rg;
894: sc->comparison = nep->sc->comparison;
895: sc->comparisonctx = nep->sc->comparisonctx;
896: BVDestroy(&ctx->V);
897: BVCreateTensor(nep->V,ctx->nmat-1,&ctx->V);
898: nep->ops->solve = NEPSolve_NLEIGS;
899: nep->ops->computevectors = NEPComputeVectors_Schur;
900: } else {
901: NEPSetUp_NLEIGS_FullBasis(nep);
902: nep->ops->solve = NEPSolve_NLEIGS_FullBasis;
903: nep->ops->computevectors = NULL;
904: }
905: PetscFunctionReturn(0);
906: }
908: /*
909: Extend the TOAR basis by applying the the matrix operator
910: over a vector which is decomposed on the TOAR way
911: Input:
912: - S,V: define the latest Arnoldi vector (nv vectors in V)
913: Output:
914: - t: new vector extending the TOAR basis
915: - r: temporally coefficients to compute the TOAR coefficients
916: for the new Arnoldi vector
917: Workspace: t_ (two vectors)
918: */
919: 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_)
920: {
921: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
922: PetscInt deg=ctx->nmat-1,k,j;
923: Vec v=t_[0],q=t_[1],w;
924: PetscScalar *beta=ctx->beta,*s=ctx->s,*xi=ctx->xi,*coeffs,sigma;
926: if (!ctx->ksp) NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp);
927: sigma = ctx->shifts[idxrktg];
928: BVSetActiveColumns(nep->V,0,nv);
929: PetscMalloc1(ctx->nmat,&coeffs);
931: /* i-part stored in (i-1) position */
932: for (j=0;j<nv;j++) {
933: 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);
934: }
935: BVSetActiveColumns(W,0,deg);
936: BVGetColumn(W,deg-1,&w);
937: BVMultVec(V,1.0/beta[deg],0,w,S+(deg-1)*ls);
938: BVRestoreColumn(W,deg-1,&w);
939: BVGetColumn(W,deg-2,&w);
940: BVMultVec(V,1.0,0.0,w,r+(deg-2)*lr);
941: BVRestoreColumn(W,deg-2,&w);
942: for (k=deg-2;k>0;k--) {
944: 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);
945: BVGetColumn(W,k-1,&w);
946: BVMultVec(V,1.0,0.0,w,r+(k-1)*lr);
947: BVRestoreColumn(W,k-1,&w);
948: }
949: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
950: for (j=0;j<ctx->nmat-2;j++) coeffs[j] = ctx->coeffD[nep->nt*j];
951: coeffs[ctx->nmat-2] = ctx->coeffD[nep->nt*(ctx->nmat-1)];
952: BVMultVec(W,1.0,0.0,v,coeffs);
953: MatMult(nep->A[0],v,q);
954: for (k=1;k<nep->nt;k++) {
955: for (j=0;j<ctx->nmat-2;j++) coeffs[j] = ctx->coeffD[nep->nt*j+k];
956: coeffs[ctx->nmat-2] = ctx->coeffD[nep->nt*(ctx->nmat-1)+k];
957: BVMultVec(W,1.0,0,v,coeffs);
958: MatMult(nep->A[k],v,t);
959: VecAXPY(q,1.0,t);
960: }
961: KSPSolve(ctx->ksp[idxrktg],q,t);
962: VecScale(t,-1.0);
963: } else {
964: for (k=0;k<deg-1;k++) {
965: BVGetColumn(W,k,&w);
966: MatMult(ctx->D[k],w,q);
967: BVRestoreColumn(W,k,&w);
968: BVInsertVec(W,k,q);
969: }
970: BVGetColumn(W,deg-1,&w);
971: MatMult(ctx->D[deg],w,q);
972: BVRestoreColumn(W,k,&w);
973: BVInsertVec(W,k,q);
974: for (j=0;j<ctx->nmat-1;j++) coeffs[j] = 1.0;
975: BVMultVec(W,1.0,0.0,q,coeffs);
976: KSPSolve(ctx->ksp[idxrktg],q,t);
977: VecScale(t,-1.0);
978: }
979: PetscFree(coeffs);
980: PetscFunctionReturn(0);
981: }
983: /*
984: Compute TOAR coefficients of the blocks of the new Arnoldi vector computed
985: */
986: static PetscErrorCode NEPTOARCoefficients(NEP nep,PetscScalar sigma,PetscInt nv,PetscScalar *S,PetscInt ls,PetscScalar *r,PetscInt lr,PetscScalar *x,PetscScalar *work)
987: {
988: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
989: PetscInt k,j,d=ctx->nmat-1;
990: PetscScalar *t=work;
992: NEPNLEIGSEvalNRTFunct(nep,d-1,sigma,t);
993: for (k=0;k<d-1;k++) {
994: for (j=0;j<=nv;j++) r[k*lr+j] += t[k]*x[j];
995: }
996: for (j=0;j<=nv;j++) r[(d-1)*lr+j] = t[d-1]*x[j];
997: PetscFunctionReturn(0);
998: }
1000: /*
1001: Compute continuation vector coefficients for the Rational-Krylov run.
1002: dim(work) >= (end-ini)*(end-ini+1) + end+1 + 2*(end-ini+1), dim(t) = end.
1003: */
1004: 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)
1005: {
1006: PetscScalar *x,*W,*tau,sone=1.0,szero=0.0;
1007: PetscInt i,j,n1,n,nwu=0;
1008: PetscBLASInt info,n_,n1_,one=1,dim,lds_;
1009: NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
1011: if (!ctx->nshifts || !end) {
1012: t[0] = 1;
1013: PetscArraycpy(cont,S+end*lds,lds);
1014: } else {
1015: n = end-ini;
1016: n1 = n+1;
1017: x = work+nwu;
1018: nwu += end+1;
1019: tau = work+nwu;
1020: nwu += n;
1021: W = work+nwu;
1022: nwu += n1*n;
1023: for (j=ini;j<end;j++) {
1024: for (i=ini;i<=end;i++) W[(j-ini)*n1+i-ini] = K[j*ld+i] -H[j*ld+i]*sigma;
1025: }
1026: PetscBLASIntCast(n,&n_);
1027: PetscBLASIntCast(n1,&n1_);
1028: PetscBLASIntCast(end+1,&dim);
1029: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1030: PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&n1_,&n_,W,&n1_,tau,work+nwu,&n1_,&info));
1031: SlepcCheckLapackInfo("geqrf",info);
1032: for (i=0;i<end;i++) t[i] = 0.0;
1033: t[end] = 1.0;
1034: for (j=n-1;j>=0;j--) {
1035: for (i=0;i<ini+j;i++) x[i] = 0.0;
1036: x[ini+j] = 1.0;
1037: for (i=j+1;i<n1;i++) x[i+ini] = W[i+n1*j];
1038: tau[j] = PetscConj(tau[j]);
1039: PetscStackCallBLAS("LAPACKlarf",LAPACKlarf_("L",&dim,&one,x,&one,tau+j,t,&dim,work+nwu));
1040: }
1041: PetscBLASIntCast(lds,&lds_);
1042: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&lds_,&n1_,&sone,S,&lds_,t,&one,&szero,cont,&one));
1043: PetscFPTrapPop();
1044: }
1045: PetscFunctionReturn(0);
1046: }
1048: /*
1049: Compute a run of Arnoldi iterations
1050: */
1051: PetscErrorCode NEPNLEIGSTOARrun(NEP nep,PetscScalar *K,PetscScalar *H,PetscInt ldh,BV W,PetscInt k,PetscInt *M,PetscBool *breakdown,Vec *t_)
1052: {
1053: NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
1054: PetscInt i,j,m=*M,lwa,deg=ctx->nmat-1,lds,nqt,ld,l;
1055: Vec t;
1056: PetscReal norm;
1057: PetscScalar *x,*work,*tt,sigma,*cont,*S;
1058: PetscBool lindep;
1059: Mat MS;
1061: BVTensorGetFactors(ctx->V,NULL,&MS);
1062: MatDenseGetArray(MS,&S);
1063: BVGetSizes(nep->V,NULL,NULL,&ld);
1064: lds = ld*deg;
1065: BVGetActiveColumns(nep->V,&l,&nqt);
1066: lwa = PetscMax(ld,deg)+(m+1)*(m+1)+4*(m+1);
1067: PetscMalloc4(ld,&x,lwa,&work,m+1,&tt,lds,&cont);
1068: BVSetActiveColumns(ctx->V,0,m);
1069: for (j=k;j<m;j++) {
1070: sigma = ctx->shifts[(++(ctx->idxrk))%ctx->nshiftsw];
1072: /* Continuation vector */
1073: NEPNLEIGS_RKcontinuation(nep,0,j,K,H,ldh,sigma,S,lds,cont,tt,work);
1075: /* apply operator */
1076: BVGetColumn(nep->V,nqt,&t);
1077: NEPTOARExtendBasis(nep,(ctx->idxrk)%ctx->nshiftsw,cont,ld,nqt,W,nep->V,t,S+(j+1)*lds,ld,t_);
1078: BVRestoreColumn(nep->V,nqt,&t);
1080: /* orthogonalize */
1081: BVOrthogonalizeColumn(nep->V,nqt,x,&norm,&lindep);
1082: if (!lindep) {
1083: x[nqt] = norm;
1084: BVScaleColumn(nep->V,nqt,1.0/norm);
1085: nqt++;
1086: } else x[nqt] = 0.0;
1088: NEPTOARCoefficients(nep,sigma,nqt-1,cont,ld,S+(j+1)*lds,ld,x,work);
1090: /* Level-2 orthogonalization */
1091: BVOrthogonalizeColumn(ctx->V,j+1,H+j*ldh,&norm,breakdown);
1092: H[j+1+ldh*j] = norm;
1093: if (ctx->nshifts) {
1094: for (i=0;i<=j;i++) K[i+ldh*j] = sigma*H[i+ldh*j] + tt[i];
1095: K[j+1+ldh*j] = sigma*H[j+1+ldh*j];
1096: }
1097: if (*breakdown) {
1098: *M = j+1;
1099: break;
1100: }
1101: BVScaleColumn(ctx->V,j+1,1.0/norm);
1102: BVSetActiveColumns(nep->V,l,nqt);
1103: }
1104: PetscFree4(x,work,tt,cont);
1105: MatDenseRestoreArray(MS,&S);
1106: BVTensorRestoreFactors(ctx->V,NULL,&MS);
1107: PetscFunctionReturn(0);
1108: }
1110: PetscErrorCode NEPSolve_NLEIGS(NEP nep)
1111: {
1112: NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
1113: PetscInt i,k=0,l,nv=0,ld,lds,ldds,nq;
1114: PetscInt deg=ctx->nmat-1,nconv=0,dsn,dsk;
1115: PetscScalar *H,*pU,*K,betak=0,*eigr,*eigi;
1116: const PetscScalar *S;
1117: PetscReal betah;
1118: PetscBool falselock=PETSC_FALSE,breakdown=PETSC_FALSE;
1119: BV W;
1120: Mat MS,MQ,U;
1122: if (ctx->lock) {
1123: /* undocumented option to use a cheaper locking instead of the true locking */
1124: PetscOptionsGetBool(NULL,NULL,"-nep_nleigs_falselocking",&falselock,NULL);
1125: }
1127: BVGetSizes(nep->V,NULL,NULL,&ld);
1128: lds = deg*ld;
1129: DSGetLeadingDimension(nep->ds,&ldds);
1130: if (!ctx->nshifts) PetscMalloc2(nep->ncv,&eigr,nep->ncv,&eigi);
1131: else { eigr = nep->eigr; eigi = nep->eigi; }
1132: BVDuplicateResize(nep->V,PetscMax(nep->nt-1,ctx->nmat-1),&W);
1134: /* clean projected matrix (including the extra-arrow) */
1135: DSGetArray(nep->ds,DS_MAT_A,&H);
1136: PetscArrayzero(H,ldds*ldds);
1137: DSRestoreArray(nep->ds,DS_MAT_A,&H);
1138: if (ctx->nshifts) {
1139: DSGetArray(nep->ds,DS_MAT_B,&H);
1140: PetscArrayzero(H,ldds*ldds);
1141: DSRestoreArray(nep->ds,DS_MAT_B,&H);
1142: }
1144: /* Get the starting Arnoldi vector */
1145: BVTensorBuildFirstColumn(ctx->V,nep->nini);
1147: /* Restart loop */
1148: l = 0;
1149: while (nep->reason == NEP_CONVERGED_ITERATING) {
1150: nep->its++;
1152: /* Compute an nv-step Krylov relation */
1153: nv = PetscMin(nep->nconv+nep->mpd,nep->ncv);
1154: if (ctx->nshifts) DSGetArray(nep->ds,DS_MAT_A,&K);
1155: DSGetArray(nep->ds,ctx->nshifts?DS_MAT_B:DS_MAT_A,&H);
1156: NEPNLEIGSTOARrun(nep,K,H,ldds,W,nep->nconv+l,&nv,&breakdown,nep->work);
1157: betah = PetscAbsScalar(H[(nv-1)*ldds+nv]);
1158: DSRestoreArray(nep->ds,ctx->nshifts?DS_MAT_B:DS_MAT_A,&H);
1159: if (ctx->nshifts) {
1160: betak = K[(nv-1)*ldds+nv];
1161: DSRestoreArray(nep->ds,DS_MAT_A,&K);
1162: }
1163: DSSetDimensions(nep->ds,nv,nep->nconv,nep->nconv+l);
1164: if (l==0) DSSetState(nep->ds,DS_STATE_INTERMEDIATE);
1165: else DSSetState(nep->ds,DS_STATE_RAW);
1167: /* Solve projected problem */
1168: DSSolve(nep->ds,nep->eigr,nep->eigi);
1169: DSSort(nep->ds,nep->eigr,nep->eigi,NULL,NULL,NULL);
1170: DSUpdateExtraRow(nep->ds);
1171: DSSynchronize(nep->ds,nep->eigr,nep->eigi);
1173: /* Check convergence */
1174: NEPNLEIGSKrylovConvergence(nep,PETSC_FALSE,nep->nconv,nv-nep->nconv,betah,betak,&k,nep->work);
1175: (*nep->stopping)(nep,nep->its,nep->max_it,k,nep->nev,&nep->reason,nep->stoppingctx);
1177: /* Update l */
1178: if (nep->reason != NEP_CONVERGED_ITERATING || breakdown) l = 0;
1179: else {
1180: l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
1181: DSGetTruncateSize(nep->ds,k,nv,&l);
1182: if (!breakdown) {
1183: /* Prepare the Rayleigh quotient for restart */
1184: DSGetDimensions(nep->ds,&dsn,NULL,&dsk,NULL);
1185: DSSetDimensions(nep->ds,dsn,k,dsk);
1186: DSTruncate(nep->ds,k+l,PETSC_FALSE);
1187: }
1188: }
1189: nconv = k;
1190: if (!ctx->lock && nep->reason == NEP_CONVERGED_ITERATING && !breakdown) { l += k; k = 0; }
1191: if (l) PetscInfo(nep,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l);
1193: /* Update S */
1194: DSGetMat(nep->ds,ctx->nshifts?DS_MAT_Z:DS_MAT_Q,&MQ);
1195: BVMultInPlace(ctx->V,MQ,nep->nconv,k+l);
1196: MatDestroy(&MQ);
1198: /* Copy last column of S */
1199: BVCopyColumn(ctx->V,nv,k+l);
1201: if (breakdown && nep->reason == NEP_CONVERGED_ITERATING) {
1202: /* Stop if breakdown */
1203: PetscInfo(nep,"Breakdown (it=%" PetscInt_FMT " norm=%g)\n",nep->its,(double)betah);
1204: nep->reason = NEP_DIVERGED_BREAKDOWN;
1205: }
1206: if (nep->reason != NEP_CONVERGED_ITERATING) l--;
1207: /* truncate S */
1208: BVGetActiveColumns(nep->V,NULL,&nq);
1209: if (k+l+deg<=nq) {
1210: BVSetActiveColumns(ctx->V,nep->nconv,k+l+1);
1211: if (!falselock && ctx->lock) BVTensorCompress(ctx->V,k-nep->nconv);
1212: else BVTensorCompress(ctx->V,0);
1213: }
1214: nep->nconv = k;
1215: if (!ctx->nshifts) {
1216: for (i=0;i<nv;i++) { eigr[i] = nep->eigr[i]; eigi[i] = nep->eigi[i]; }
1217: NEPNLEIGSBackTransform((PetscObject)nep,nv,eigr,eigi);
1218: }
1219: NEPMonitor(nep,nep->its,nconv,eigr,eigi,nep->errest,nv);
1220: }
1221: nep->nconv = nconv;
1222: if (nep->nconv>0) {
1223: BVSetActiveColumns(ctx->V,0,nep->nconv);
1224: BVGetActiveColumns(nep->V,NULL,&nq);
1225: BVSetActiveColumns(nep->V,0,nq);
1226: if (nq>nep->nconv) {
1227: BVTensorCompress(ctx->V,nep->nconv);
1228: BVSetActiveColumns(nep->V,0,nep->nconv);
1229: nq = nep->nconv;
1230: }
1231: if (ctx->nshifts) {
1232: DSGetMat(nep->ds,DS_MAT_B,&MQ);
1233: BVMultInPlace(ctx->V,MQ,0,nep->nconv);
1234: MatDestroy(&MQ);
1235: }
1236: BVTensorGetFactors(ctx->V,NULL,&MS);
1237: MatDenseGetArrayRead(MS,&S);
1238: PetscMalloc1(nq*nep->nconv,&pU);
1239: for (i=0;i<nep->nconv;i++) PetscArraycpy(pU+i*nq,S+i*lds,nq);
1240: MatDenseRestoreArrayRead(MS,&S);
1241: BVTensorRestoreFactors(ctx->V,NULL,&MS);
1242: MatCreateSeqDense(PETSC_COMM_SELF,nq,nep->nconv,pU,&U);
1243: BVSetActiveColumns(nep->V,0,nq);
1244: BVMultInPlace(nep->V,U,0,nep->nconv);
1245: BVSetActiveColumns(nep->V,0,nep->nconv);
1246: MatDestroy(&U);
1247: PetscFree(pU);
1248: DSTruncate(nep->ds,nep->nconv,PETSC_TRUE);
1249: }
1251: /* Map eigenvalues back to the original problem */
1252: if (!ctx->nshifts) {
1253: NEPNLEIGSBackTransform((PetscObject)nep,nep->nconv,nep->eigr,nep->eigi);
1254: PetscFree2(eigr,eigi);
1255: }
1256: BVDestroy(&W);
1257: PetscFunctionReturn(0);
1258: }
1260: static PetscErrorCode NEPNLEIGSSetSingularitiesFunction_NLEIGS(NEP nep,PetscErrorCode (*fun)(NEP,PetscInt*,PetscScalar*,void*),void *ctx)
1261: {
1262: NEP_NLEIGS *nepctx=(NEP_NLEIGS*)nep->data;
1264: if (fun) nepctx->computesingularities = fun;
1265: if (ctx) nepctx->singularitiesctx = ctx;
1266: nep->state = NEP_STATE_INITIAL;
1267: PetscFunctionReturn(0);
1268: }
1270: /*@C
1271: NEPNLEIGSSetSingularitiesFunction - Sets a user function to compute a discretization
1272: of the singularity set (where T(.) is not analytic).
1274: Logically Collective on nep
1276: Input Parameters:
1277: + nep - the NEP context
1278: . fun - user function (if NULL then NEP retains any previously set value)
1279: - ctx - [optional] user-defined context for private data for the function
1280: (may be NULL, in which case NEP retains any previously set value)
1282: Calling Sequence of fun:
1283: $ fun(NEP nep,PetscInt *maxnp,PetscScalar *xi,void *ctx)
1285: + nep - the NEP context
1286: . maxnp - on input number of requested points in the discretization (can be set)
1287: . xi - computed values of the discretization
1288: - ctx - optional context, as set by NEPNLEIGSSetSingularitiesFunction()
1290: Notes:
1291: The user-defined function can set a smaller value of maxnp if necessary.
1292: It is wrong to return a larger value.
1294: If the problem type has been set to rational with NEPSetProblemType(),
1295: then it is not necessary to set the singularities explicitly since the
1296: solver will try to determine them automatically.
1298: Level: intermediate
1300: .seealso: NEPNLEIGSGetSingularitiesFunction(), NEPSetProblemType()
1301: @*/
1302: PetscErrorCode NEPNLEIGSSetSingularitiesFunction(NEP nep,PetscErrorCode (*fun)(NEP,PetscInt*,PetscScalar*,void*),void *ctx)
1303: {
1305: PetscTryMethod(nep,"NEPNLEIGSSetSingularitiesFunction_C",(NEP,PetscErrorCode(*)(NEP,PetscInt*,PetscScalar*,void*),void*),(nep,fun,ctx));
1306: PetscFunctionReturn(0);
1307: }
1309: static PetscErrorCode NEPNLEIGSGetSingularitiesFunction_NLEIGS(NEP nep,PetscErrorCode (**fun)(NEP,PetscInt*,PetscScalar*,void*),void **ctx)
1310: {
1311: NEP_NLEIGS *nepctx=(NEP_NLEIGS*)nep->data;
1313: if (fun) *fun = nepctx->computesingularities;
1314: if (ctx) *ctx = nepctx->singularitiesctx;
1315: PetscFunctionReturn(0);
1316: }
1318: /*@C
1319: NEPNLEIGSGetSingularitiesFunction - Returns the Function and optionally the user
1320: provided context for computing a discretization of the singularity set.
1322: Not Collective
1324: Input Parameter:
1325: . nep - the nonlinear eigensolver context
1327: Output Parameters:
1328: + fun - location to put the function (or NULL)
1329: - ctx - location to stash the function context (or NULL)
1331: Level: advanced
1333: .seealso: NEPNLEIGSSetSingularitiesFunction()
1334: @*/
1335: PetscErrorCode NEPNLEIGSGetSingularitiesFunction(NEP nep,PetscErrorCode (**fun)(NEP,PetscInt*,PetscScalar*,void*),void **ctx)
1336: {
1338: PetscUseMethod(nep,"NEPNLEIGSGetSingularitiesFunction_C",(NEP,PetscErrorCode(**)(NEP,PetscInt*,PetscScalar*,void*),void**),(nep,fun,ctx));
1339: PetscFunctionReturn(0);
1340: }
1342: static PetscErrorCode NEPNLEIGSSetRestart_NLEIGS(NEP nep,PetscReal keep)
1343: {
1344: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1346: if (keep==PETSC_DEFAULT) ctx->keep = 0.5;
1347: else {
1349: ctx->keep = keep;
1350: }
1351: PetscFunctionReturn(0);
1352: }
1354: /*@
1355: NEPNLEIGSSetRestart - Sets the restart parameter for the NLEIGS
1356: method, in particular the proportion of basis vectors that must be kept
1357: after restart.
1359: Logically Collective on nep
1361: Input Parameters:
1362: + nep - the nonlinear eigensolver context
1363: - keep - the number of vectors to be kept at restart
1365: Options Database Key:
1366: . -nep_nleigs_restart - Sets the restart parameter
1368: Notes:
1369: Allowed values are in the range [0.1,0.9]. The default is 0.5.
1371: Level: advanced
1373: .seealso: NEPNLEIGSGetRestart()
1374: @*/
1375: PetscErrorCode NEPNLEIGSSetRestart(NEP nep,PetscReal keep)
1376: {
1379: PetscTryMethod(nep,"NEPNLEIGSSetRestart_C",(NEP,PetscReal),(nep,keep));
1380: PetscFunctionReturn(0);
1381: }
1383: static PetscErrorCode NEPNLEIGSGetRestart_NLEIGS(NEP nep,PetscReal *keep)
1384: {
1385: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1387: *keep = ctx->keep;
1388: PetscFunctionReturn(0);
1389: }
1391: /*@
1392: NEPNLEIGSGetRestart - Gets the restart parameter used in the NLEIGS method.
1394: Not Collective
1396: Input Parameter:
1397: . nep - the nonlinear eigensolver context
1399: Output Parameter:
1400: . keep - the restart parameter
1402: Level: advanced
1404: .seealso: NEPNLEIGSSetRestart()
1405: @*/
1406: PetscErrorCode NEPNLEIGSGetRestart(NEP nep,PetscReal *keep)
1407: {
1410: PetscUseMethod(nep,"NEPNLEIGSGetRestart_C",(NEP,PetscReal*),(nep,keep));
1411: PetscFunctionReturn(0);
1412: }
1414: static PetscErrorCode NEPNLEIGSSetLocking_NLEIGS(NEP nep,PetscBool lock)
1415: {
1416: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1418: ctx->lock = lock;
1419: PetscFunctionReturn(0);
1420: }
1422: /*@
1423: NEPNLEIGSSetLocking - Choose between locking and non-locking variants of
1424: the NLEIGS method.
1426: Logically Collective on nep
1428: Input Parameters:
1429: + nep - the nonlinear eigensolver context
1430: - lock - true if the locking variant must be selected
1432: Options Database Key:
1433: . -nep_nleigs_locking - Sets the locking flag
1435: Notes:
1436: The default is to lock converged eigenpairs when the method restarts.
1437: This behaviour can be changed so that all directions are kept in the
1438: working subspace even if already converged to working accuracy (the
1439: non-locking variant).
1441: Level: advanced
1443: .seealso: NEPNLEIGSGetLocking()
1444: @*/
1445: PetscErrorCode NEPNLEIGSSetLocking(NEP nep,PetscBool lock)
1446: {
1449: PetscTryMethod(nep,"NEPNLEIGSSetLocking_C",(NEP,PetscBool),(nep,lock));
1450: PetscFunctionReturn(0);
1451: }
1453: static PetscErrorCode NEPNLEIGSGetLocking_NLEIGS(NEP nep,PetscBool *lock)
1454: {
1455: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1457: *lock = ctx->lock;
1458: PetscFunctionReturn(0);
1459: }
1461: /*@
1462: NEPNLEIGSGetLocking - Gets the locking flag used in the NLEIGS method.
1464: Not Collective
1466: Input Parameter:
1467: . nep - the nonlinear eigensolver context
1469: Output Parameter:
1470: . lock - the locking flag
1472: Level: advanced
1474: .seealso: NEPNLEIGSSetLocking()
1475: @*/
1476: PetscErrorCode NEPNLEIGSGetLocking(NEP nep,PetscBool *lock)
1477: {
1480: PetscUseMethod(nep,"NEPNLEIGSGetLocking_C",(NEP,PetscBool*),(nep,lock));
1481: PetscFunctionReturn(0);
1482: }
1484: static PetscErrorCode NEPNLEIGSSetInterpolation_NLEIGS(NEP nep,PetscReal tol,PetscInt degree)
1485: {
1486: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1488: if (tol == PETSC_DEFAULT) {
1489: ctx->ddtol = PETSC_DEFAULT;
1490: nep->state = NEP_STATE_INITIAL;
1491: } else {
1493: ctx->ddtol = tol;
1494: }
1495: if (degree == PETSC_DEFAULT || degree == PETSC_DECIDE) {
1496: ctx->ddmaxit = 0;
1497: if (nep->state) NEPReset(nep);
1498: nep->state = NEP_STATE_INITIAL;
1499: } else {
1501: if (ctx->ddmaxit != degree) {
1502: ctx->ddmaxit = degree;
1503: if (nep->state) NEPReset(nep);
1504: nep->state = NEP_STATE_INITIAL;
1505: }
1506: }
1507: PetscFunctionReturn(0);
1508: }
1510: /*@
1511: NEPNLEIGSSetInterpolation - Sets the tolerance and maximum degree
1512: when building the interpolation via divided differences.
1514: Logically Collective on nep
1516: Input Parameters:
1517: + nep - the nonlinear eigensolver context
1518: . tol - tolerance to stop computing divided differences
1519: - degree - maximum degree of interpolation
1521: Options Database Key:
1522: + -nep_nleigs_interpolation_tol <tol> - Sets the tolerance to stop computing divided differences
1523: - -nep_nleigs_interpolation_degree <degree> - Sets the maximum degree of interpolation
1525: Notes:
1526: Use PETSC_DEFAULT for either argument to assign a reasonably good value.
1528: Level: advanced
1530: .seealso: NEPNLEIGSGetInterpolation()
1531: @*/
1532: PetscErrorCode NEPNLEIGSSetInterpolation(NEP nep,PetscReal tol,PetscInt degree)
1533: {
1537: PetscTryMethod(nep,"NEPNLEIGSSetInterpolation_C",(NEP,PetscReal,PetscInt),(nep,tol,degree));
1538: PetscFunctionReturn(0);
1539: }
1541: static PetscErrorCode NEPNLEIGSGetInterpolation_NLEIGS(NEP nep,PetscReal *tol,PetscInt *degree)
1542: {
1543: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1545: if (tol) *tol = ctx->ddtol;
1546: if (degree) *degree = ctx->ddmaxit;
1547: PetscFunctionReturn(0);
1548: }
1550: /*@
1551: NEPNLEIGSGetInterpolation - Gets the tolerance and maximum degree
1552: when building the interpolation via divided differences.
1554: Not Collective
1556: Input Parameter:
1557: . nep - the nonlinear eigensolver context
1559: Output Parameters:
1560: + tol - tolerance to stop computing divided differences
1561: - degree - maximum degree of interpolation
1563: Level: advanced
1565: .seealso: NEPNLEIGSSetInterpolation()
1566: @*/
1567: PetscErrorCode NEPNLEIGSGetInterpolation(NEP nep,PetscReal *tol,PetscInt *degree)
1568: {
1570: PetscTryMethod(nep,"NEPNLEIGSGetInterpolation_C",(NEP,PetscReal*,PetscInt*),(nep,tol,degree));
1571: PetscFunctionReturn(0);
1572: }
1574: static PetscErrorCode NEPNLEIGSSetRKShifts_NLEIGS(NEP nep,PetscInt ns,PetscScalar *shifts)
1575: {
1576: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1577: PetscInt i;
1580: if (ctx->nshifts) PetscFree(ctx->shifts);
1581: for (i=0;i<ctx->nshiftsw;i++) KSPDestroy(&ctx->ksp[i]);
1582: PetscFree(ctx->ksp);
1583: ctx->ksp = NULL;
1584: if (ns) {
1585: PetscMalloc1(ns,&ctx->shifts);
1586: for (i=0;i<ns;i++) ctx->shifts[i] = shifts[i];
1587: }
1588: ctx->nshifts = ns;
1589: nep->state = NEP_STATE_INITIAL;
1590: PetscFunctionReturn(0);
1591: }
1593: /*@
1594: NEPNLEIGSSetRKShifts - Sets a list of shifts to be used in the Rational
1595: Krylov method.
1597: Logically Collective on nep
1599: Input Parameters:
1600: + nep - the nonlinear eigensolver context
1601: . ns - number of shifts
1602: - shifts - array of scalar values specifying the shifts
1604: Options Database Key:
1605: . -nep_nleigs_rk_shifts - Sets the list of shifts
1607: Notes:
1608: If only one shift is provided, the built subspace built is equivalent to
1609: shift-and-invert Krylov-Schur (provided that the absolute convergence
1610: criterion is used).
1612: In the case of real scalars, complex shifts are not allowed. In the
1613: command line, a comma-separated list of complex values can be provided with
1614: the format [+/-][realnumber][+/-]realnumberi with no spaces, e.g.
1615: -nep_nleigs_rk_shifts 1.0+2.0i,1.5+2.0i,1.0+1.5i
1617: Use ns=0 to remove previously set shifts.
1619: Level: advanced
1621: .seealso: NEPNLEIGSGetRKShifts()
1622: @*/
1623: PetscErrorCode NEPNLEIGSSetRKShifts(NEP nep,PetscInt ns,PetscScalar shifts[])
1624: {
1628: PetscTryMethod(nep,"NEPNLEIGSSetRKShifts_C",(NEP,PetscInt,PetscScalar*),(nep,ns,shifts));
1629: PetscFunctionReturn(0);
1630: }
1632: static PetscErrorCode NEPNLEIGSGetRKShifts_NLEIGS(NEP nep,PetscInt *ns,PetscScalar **shifts)
1633: {
1634: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1635: PetscInt i;
1637: *ns = ctx->nshifts;
1638: if (ctx->nshifts) {
1639: PetscMalloc1(ctx->nshifts,shifts);
1640: for (i=0;i<ctx->nshifts;i++) (*shifts)[i] = ctx->shifts[i];
1641: }
1642: PetscFunctionReturn(0);
1643: }
1645: /*@C
1646: NEPNLEIGSGetRKShifts - Gets the list of shifts used in the Rational
1647: Krylov method.
1649: Not Collective
1651: Input Parameter:
1652: . nep - the nonlinear eigensolver context
1654: Output Parameters:
1655: + ns - number of shifts
1656: - shifts - array of shifts
1658: Note:
1659: The user is responsible for deallocating the returned array.
1661: Level: advanced
1663: .seealso: NEPNLEIGSSetRKShifts()
1664: @*/
1665: PetscErrorCode NEPNLEIGSGetRKShifts(NEP nep,PetscInt *ns,PetscScalar *shifts[])
1666: {
1670: PetscTryMethod(nep,"NEPNLEIGSGetRKShifts_C",(NEP,PetscInt*,PetscScalar**),(nep,ns,shifts));
1671: PetscFunctionReturn(0);
1672: }
1674: static PetscErrorCode NEPNLEIGSGetKSPs_NLEIGS(NEP nep,PetscInt *nsolve,KSP **ksp)
1675: {
1676: NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
1677: PetscInt i;
1678: PC pc;
1680: if (!ctx->ksp) {
1681: NEPNLEIGSSetShifts(nep,&ctx->nshiftsw);
1682: PetscMalloc1(ctx->nshiftsw,&ctx->ksp);
1683: for (i=0;i<ctx->nshiftsw;i++) {
1684: KSPCreate(PetscObjectComm((PetscObject)nep),&ctx->ksp[i]);
1685: PetscObjectIncrementTabLevel((PetscObject)ctx->ksp[i],(PetscObject)nep,1);
1686: KSPSetOptionsPrefix(ctx->ksp[i],((PetscObject)nep)->prefix);
1687: KSPAppendOptionsPrefix(ctx->ksp[i],"nep_nleigs_");
1688: PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->ksp[i]);
1689: PetscObjectSetOptions((PetscObject)ctx->ksp[i],((PetscObject)nep)->options);
1690: KSPSetErrorIfNotConverged(ctx->ksp[i],PETSC_TRUE);
1691: KSPSetTolerances(ctx->ksp[i],1e-3*SlepcDefaultTol(nep->tol),PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
1692: KSPGetPC(ctx->ksp[i],&pc);
1693: if ((nep->fui==NEP_USER_INTERFACE_SPLIT && nep->P) || (nep->fui==NEP_USER_INTERFACE_CALLBACK && nep->function_pre!=nep->function)) {
1694: KSPSetType(ctx->ksp[i],KSPBCGS);
1695: PCSetType(pc,PCBJACOBI);
1696: } else {
1697: KSPSetType(ctx->ksp[i],KSPPREONLY);
1698: PCSetType(pc,PCLU);
1699: }
1700: }
1701: }
1702: if (nsolve) *nsolve = ctx->nshiftsw;
1703: if (ksp) *ksp = ctx->ksp;
1704: PetscFunctionReturn(0);
1705: }
1707: /*@C
1708: NEPNLEIGSGetKSPs - Retrieve the array of linear solver objects associated with
1709: the nonlinear eigenvalue solver.
1711: Not Collective
1713: Input Parameter:
1714: . nep - nonlinear eigenvalue solver
1716: Output Parameters:
1717: + nsolve - number of returned KSP objects
1718: - ksp - array of linear solver object
1720: Notes:
1721: The number of KSP objects is equal to the number of shifts provided by the user,
1722: or 1 if the user did not provide shifts.
1724: Level: advanced
1726: .seealso: NEPNLEIGSSetRKShifts()
1727: @*/
1728: PetscErrorCode NEPNLEIGSGetKSPs(NEP nep,PetscInt *nsolve,KSP **ksp)
1729: {
1731: PetscUseMethod(nep,"NEPNLEIGSGetKSPs_C",(NEP,PetscInt*,KSP**),(nep,nsolve,ksp));
1732: PetscFunctionReturn(0);
1733: }
1735: static PetscErrorCode NEPNLEIGSSetFullBasis_NLEIGS(NEP nep,PetscBool fullbasis)
1736: {
1737: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1739: if (fullbasis!=ctx->fullbasis) {
1740: ctx->fullbasis = fullbasis;
1741: nep->state = NEP_STATE_INITIAL;
1742: nep->useds = PetscNot(fullbasis);
1743: }
1744: PetscFunctionReturn(0);
1745: }
1747: /*@
1748: NEPNLEIGSSetFullBasis - Choose between TOAR-basis (default) and full-basis
1749: variants of the NLEIGS method.
1751: Logically Collective on nep
1753: Input Parameters:
1754: + nep - the nonlinear eigensolver context
1755: - fullbasis - true if the full-basis variant must be selected
1757: Options Database Key:
1758: . -nep_nleigs_full_basis - Sets the full-basis flag
1760: Notes:
1761: The default is to use a compact representation of the Krylov basis, that is,
1762: V = (I otimes U) S, with a tensor BV. This behaviour can be changed so that
1763: the full basis V is explicitly stored and operated with. This variant is more
1764: expensive in terms of memory and computation, but is necessary in some cases,
1765: particularly for two-sided computations, see NEPSetTwoSided().
1767: In the full-basis variant, the NLEIGS solver uses an EPS object to explicitly
1768: solve the linearized eigenproblem, see NEPNLEIGSGetEPS().
1770: Level: advanced
1772: .seealso: NEPNLEIGSGetFullBasis(), NEPNLEIGSGetEPS(), NEPSetTwoSided(), BVCreateTensor()
1773: @*/
1774: PetscErrorCode NEPNLEIGSSetFullBasis(NEP nep,PetscBool fullbasis)
1775: {
1778: PetscTryMethod(nep,"NEPNLEIGSSetFullBasis_C",(NEP,PetscBool),(nep,fullbasis));
1779: PetscFunctionReturn(0);
1780: }
1782: static PetscErrorCode NEPNLEIGSGetFullBasis_NLEIGS(NEP nep,PetscBool *fullbasis)
1783: {
1784: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1786: *fullbasis = ctx->fullbasis;
1787: PetscFunctionReturn(0);
1788: }
1790: /*@
1791: NEPNLEIGSGetFullBasis - Gets the flag that indicates if NLEIGS is using the
1792: full-basis variant.
1794: Not Collective
1796: Input Parameter:
1797: . nep - the nonlinear eigensolver context
1799: Output Parameter:
1800: . fullbasis - the flag
1802: Level: advanced
1804: .seealso: NEPNLEIGSSetFullBasis()
1805: @*/
1806: PetscErrorCode NEPNLEIGSGetFullBasis(NEP nep,PetscBool *fullbasis)
1807: {
1810: PetscUseMethod(nep,"NEPNLEIGSGetFullBasis_C",(NEP,PetscBool*),(nep,fullbasis));
1811: PetscFunctionReturn(0);
1812: }
1814: #define SHIFTMAX 30
1816: PetscErrorCode NEPSetFromOptions_NLEIGS(PetscOptionItems *PetscOptionsObject,NEP nep)
1817: {
1818: NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
1819: PetscInt i=0,k;
1820: PetscBool flg1,flg2,b;
1821: PetscReal r;
1822: PetscScalar array[SHIFTMAX];
1824: PetscOptionsHead(PetscOptionsObject,"NEP NLEIGS Options");
1826: PetscOptionsReal("-nep_nleigs_restart","Proportion of vectors kept after restart","NEPNLEIGSSetRestart",0.5,&r,&flg1);
1827: if (flg1) NEPNLEIGSSetRestart(nep,r);
1829: PetscOptionsBool("-nep_nleigs_locking","Choose between locking and non-locking variants","NEPNLEIGSSetLocking",PETSC_FALSE,&b,&flg1);
1830: if (flg1) NEPNLEIGSSetLocking(nep,b);
1832: PetscOptionsBool("-nep_nleigs_full_basis","Choose between TOAR and full-basis variants","NEPNLEIGSSetFullBasis",PETSC_FALSE,&b,&flg1);
1833: if (flg1) NEPNLEIGSSetFullBasis(nep,b);
1835: NEPNLEIGSGetInterpolation(nep,&r,&i);
1836: if (!i) i = PETSC_DEFAULT;
1837: PetscOptionsInt("-nep_nleigs_interpolation_degree","Maximum number of terms for interpolation via divided differences","NEPNLEIGSSetInterpolation",i,&i,&flg1);
1838: PetscOptionsReal("-nep_nleigs_interpolation_tol","Tolerance for interpolation via divided differences","NEPNLEIGSSetInterpolation",r,&r,&flg2);
1839: if (flg1 || flg2) NEPNLEIGSSetInterpolation(nep,r,i);
1841: k = SHIFTMAX;
1842: for (i=0;i<k;i++) array[i] = 0;
1843: PetscOptionsScalarArray("-nep_nleigs_rk_shifts","Shifts for Rational Krylov","NEPNLEIGSSetRKShifts",array,&k,&flg1);
1844: if (flg1) NEPNLEIGSSetRKShifts(nep,k,array);
1846: PetscOptionsTail();
1848: if (!ctx->ksp) NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp);
1849: for (i=0;i<ctx->nshiftsw;i++) KSPSetFromOptions(ctx->ksp[i]);
1851: if (ctx->fullbasis) {
1852: if (!ctx->eps) NEPNLEIGSGetEPS(nep,&ctx->eps);
1853: EPSSetFromOptions(ctx->eps);
1854: }
1855: PetscFunctionReturn(0);
1856: }
1858: PetscErrorCode NEPView_NLEIGS(NEP nep,PetscViewer viewer)
1859: {
1860: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1861: PetscBool isascii;
1862: PetscInt i;
1863: char str[50];
1865: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
1866: if (isascii) {
1867: PetscViewerASCIIPrintf(viewer," %d%% of basis vectors kept after restart\n",(int)(100*ctx->keep));
1868: if (ctx->fullbasis) PetscViewerASCIIPrintf(viewer," using the full-basis variant\n");
1869: else PetscViewerASCIIPrintf(viewer," using the %slocking variant\n",ctx->lock?"":"non-");
1870: PetscViewerASCIIPrintf(viewer," divided difference terms: used=%" PetscInt_FMT ", max=%" PetscInt_FMT "\n",ctx->nmat,ctx->ddmaxit);
1871: PetscViewerASCIIPrintf(viewer," tolerance for divided difference convergence: %g\n",(double)ctx->ddtol);
1872: if (ctx->nshifts) {
1873: PetscViewerASCIIPrintf(viewer," RK shifts: ");
1874: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1875: for (i=0;i<ctx->nshifts;i++) {
1876: SlepcSNPrintfScalar(str,sizeof(str),ctx->shifts[i],PETSC_FALSE);
1877: PetscViewerASCIIPrintf(viewer,"%s%s",str,(i<ctx->nshifts-1)?",":"");
1878: }
1879: PetscViewerASCIIPrintf(viewer,"\n");
1880: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1881: }
1882: if (!ctx->ksp) NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp);
1883: PetscViewerASCIIPushTab(viewer);
1884: KSPView(ctx->ksp[0],viewer);
1885: PetscViewerASCIIPopTab(viewer);
1886: if (ctx->fullbasis) {
1887: if (!ctx->eps) NEPNLEIGSGetEPS(nep,&ctx->eps);
1888: PetscViewerASCIIPushTab(viewer);
1889: EPSView(ctx->eps,viewer);
1890: PetscViewerASCIIPopTab(viewer);
1891: }
1892: }
1893: PetscFunctionReturn(0);
1894: }
1896: PetscErrorCode NEPReset_NLEIGS(NEP nep)
1897: {
1898: PetscInt k;
1899: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1901: if (nep->fui==NEP_USER_INTERFACE_SPLIT) PetscFree(ctx->coeffD);
1902: else {
1903: for (k=0;k<ctx->nmat;k++) MatDestroy(&ctx->D[k]);
1904: }
1905: PetscFree4(ctx->s,ctx->xi,ctx->beta,ctx->D);
1906: for (k=0;k<ctx->nshiftsw;k++) KSPReset(ctx->ksp[k]);
1907: VecDestroy(&ctx->vrn);
1908: if (ctx->fullbasis) {
1909: MatDestroy(&ctx->A);
1910: EPSReset(ctx->eps);
1911: for (k=0;k<4;k++) VecDestroy(&ctx->w[k]);
1912: }
1913: PetscFunctionReturn(0);
1914: }
1916: PetscErrorCode NEPDestroy_NLEIGS(NEP nep)
1917: {
1918: PetscInt k;
1919: NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
1921: BVDestroy(&ctx->V);
1922: for (k=0;k<ctx->nshiftsw;k++) KSPDestroy(&ctx->ksp[k]);
1923: PetscFree(ctx->ksp);
1924: if (ctx->nshifts) PetscFree(ctx->shifts);
1925: if (ctx->fullbasis) EPSDestroy(&ctx->eps);
1926: PetscFree(nep->data);
1927: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetSingularitiesFunction_C",NULL);
1928: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetSingularitiesFunction_C",NULL);
1929: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRestart_C",NULL);
1930: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRestart_C",NULL);
1931: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetLocking_C",NULL);
1932: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetLocking_C",NULL);
1933: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetInterpolation_C",NULL);
1934: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetInterpolation_C",NULL);
1935: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRKShifts_C",NULL);
1936: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRKShifts_C",NULL);
1937: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetKSPs_C",NULL);
1938: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetFullBasis_C",NULL);
1939: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetFullBasis_C",NULL);
1940: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetEPS_C",NULL);
1941: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetEPS_C",NULL);
1942: PetscFunctionReturn(0);
1943: }
1945: SLEPC_EXTERN PetscErrorCode NEPCreate_NLEIGS(NEP nep)
1946: {
1947: NEP_NLEIGS *ctx;
1949: PetscNewLog(nep,&ctx);
1950: nep->data = (void*)ctx;
1951: ctx->lock = PETSC_TRUE;
1952: ctx->ddtol = PETSC_DEFAULT;
1954: nep->useds = PETSC_TRUE;
1956: nep->ops->setup = NEPSetUp_NLEIGS;
1957: nep->ops->setfromoptions = NEPSetFromOptions_NLEIGS;
1958: nep->ops->view = NEPView_NLEIGS;
1959: nep->ops->destroy = NEPDestroy_NLEIGS;
1960: nep->ops->reset = NEPReset_NLEIGS;
1962: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetSingularitiesFunction_C",NEPNLEIGSSetSingularitiesFunction_NLEIGS);
1963: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetSingularitiesFunction_C",NEPNLEIGSGetSingularitiesFunction_NLEIGS);
1964: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRestart_C",NEPNLEIGSSetRestart_NLEIGS);
1965: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRestart_C",NEPNLEIGSGetRestart_NLEIGS);
1966: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetLocking_C",NEPNLEIGSSetLocking_NLEIGS);
1967: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetLocking_C",NEPNLEIGSGetLocking_NLEIGS);
1968: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetInterpolation_C",NEPNLEIGSSetInterpolation_NLEIGS);
1969: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetInterpolation_C",NEPNLEIGSGetInterpolation_NLEIGS);
1970: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRKShifts_C",NEPNLEIGSSetRKShifts_NLEIGS);
1971: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRKShifts_C",NEPNLEIGSGetRKShifts_NLEIGS);
1972: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetKSPs_C",NEPNLEIGSGetKSPs_NLEIGS);
1973: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetFullBasis_C",NEPNLEIGSSetFullBasis_NLEIGS);
1974: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetFullBasis_C",NEPNLEIGSGetFullBasis_NLEIGS);
1975: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetEPS_C",NEPNLEIGSSetEPS_NLEIGS);
1976: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetEPS_C",NEPNLEIGSGetEPS_NLEIGS);
1977: PetscFunctionReturn(0);
1978: }