Actual source code: nleigs.c
slepc-3.9.2 2018-07-02
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2018, 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>
29: #define LBPOINTS 100 /* default value of the maximum number of Leja-Bagby points */
30: #define NDPOINTS 1e4 /* number of discretization points */
32: typedef struct {
33: BV V; /* tensor vector basis for the linearization */
34: PetscInt nmat; /* number of interpolation points */
35: PetscScalar *s,*xi; /* Leja-Bagby points */
36: PetscScalar *beta; /* scaling factors */
37: Mat *D; /* divided difference matrices */
38: PetscScalar *coeffD; /* coefficients for divided differences in split form */
39: PetscInt nshifts; /* provided number of shifts */
40: PetscScalar *shifts; /* user-provided shifts for the Rational Krylov variant */
41: PetscInt nshiftsw; /* actual number of shifts (1 if Krylov-Schur) */
42: PetscReal ddtol; /* tolerance for divided difference convergence */
43: PetscInt ddmaxit; /* maximum number of divided difference terms */
44: PetscReal keep; /* restart parameter */
45: PetscBool lock; /* locking/non-locking variant */
46: PetscInt idxrk; /* index of next shift to use */
47: KSP *ksp; /* ksp array for storing shift factorizations */
48: Vec vrn; /* random vector with normally distributed value */
49: void *singularitiesctx;
50: PetscErrorCode (*computesingularities)(NEP,PetscInt*,PetscScalar*,void*);
51: } NEP_NLEIGS;
53: typedef struct {
54: PetscInt nmat,maxnmat;
55: PetscScalar *coeff;
56: Mat *A;
57: Vec t;
58: } ShellMatCtx;
60: PETSC_STATIC_INLINE PetscErrorCode NEPNLEIGSSetShifts(NEP nep,PetscInt *nshiftsw)
61: {
62: NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
65: if (!ctx->nshifts) {
66: ctx->shifts = &nep->target;
67: *nshiftsw = 1;
68: } else *nshiftsw = ctx->nshifts;
69: return(0);
70: }
72: static PetscErrorCode NEPNLEIGSBackTransform(PetscObject ob,PetscInt n,PetscScalar *valr,PetscScalar *vali)
73: {
74: NEP nep;
75: PetscInt j;
76: #if !defined(PETSC_USE_COMPLEX)
77: PetscScalar t;
78: #endif
81: nep = (NEP)ob;
82: #if !defined(PETSC_USE_COMPLEX)
83: for (j=0;j<n;j++) {
84: if (vali[j] == 0) valr[j] = 1.0 / valr[j] + nep->target;
85: else {
86: t = valr[j] * valr[j] + vali[j] * vali[j];
87: valr[j] = valr[j] / t + nep->target;
88: vali[j] = - vali[j] / t;
89: }
90: }
91: #else
92: for (j=0;j<n;j++) {
93: valr[j] = 1.0 / valr[j] + nep->target;
94: }
95: #endif
96: return(0);
97: }
99: /* Computes the roots of a polynomial */
100: static PetscErrorCode NEPNLEIGSAuxiliarPRootFinder(PetscInt deg,PetscScalar *polcoeffs,PetscScalar *wr,PetscScalar *wi,PetscBool *avail)
101: {
103: PetscScalar *C,*work;
104: PetscBLASInt n_,info,lwork;
105: PetscInt i;
106: #if defined(PETSC_USE_COMPLEX)
107: PetscReal *rwork;
108: #endif
109: #if defined(PETSC_HAVE_ESSL)
110: PetscScalar sdummy;
111: PetscBLASInt idummy,io=0;
112: PetscScalar *wri;
113: #endif
116: #if defined(PETSC_MISSING_LAPACK_GEEV)
117: *avail = PETSC_FALSE;
118: #else
119: *avail = PETSC_TRUE;
120: if (deg>0) {
121: PetscCalloc1(deg*deg,&C);
122: PetscBLASIntCast(deg,&n_);
123: for (i=0;i<deg-1;i++) {
124: C[(deg+1)*i+1] = 1.0;
125: C[(deg-1)*deg+i] = -polcoeffs[deg-i]/polcoeffs[0];
126: }
127: C[deg*deg+-1] = -polcoeffs[1]/polcoeffs[0];
128: PetscBLASIntCast(3*deg,&lwork);
130: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
131: #if !defined(PETSC_HAVE_ESSL)
132: #if !defined(PETSC_USE_COMPLEX)
133: PetscMalloc1(lwork,&work);
134: PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_("N","N",&n_,C,&n_,wr,wi,NULL,&n_,NULL,&n_,work,&lwork,&info));
135: if (info) *avail = PETSC_FALSE;
136: PetscFree(work);
137: #else
138: PetscMalloc2(2*deg,&rwork,lwork,&work);
139: PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_("N","N",&n_,C,&n_,wr,NULL,&n_,NULL,&n_,work,&lwork,rwork,&info));
140: if (info) *avail = PETSC_FALSE;
141: PetscFree2(rwork,work);
142: #endif
143: #else
144: PetscMalloc2(lwork,&work,2*deg,&wri);
145: PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_(&io,C,&n_,wri,&sdummy,&idummy,&idummy,&n_,work,&lwork));
146: #if !defined(PETSC_USE_COMPLEX)
147: for (i=0;i<deg;i++) {
148: wr[i] = wri[2*i];
149: wi[i] = wri[2*i+1];
150: }
151: #else
152: for (i=0;i<deg;i++) wr[i] = wri[i];
153: #endif
154: PetscFree2(work,wri);
155: #endif
156: PetscFPTrapPop();
157: PetscFree(C);
158: }
159: #endif
160: return(0);
161: }
163: static PetscErrorCode NEPNLEIGSAuxiliarRmDuplicates(PetscInt nin,PetscScalar *pin,PetscInt *nout,PetscScalar *pout)
164: {
165: PetscInt i,j;
168: for (i=0;i<nin;i++) {
169: pout[(*nout)++] = pin[i];
170: for (j=0;j<*nout-1;j++)
171: if (PetscAbsScalar(pin[i]-pout[j])<PETSC_MACHINE_EPSILON*100) {
172: (*nout)--;
173: break;
174: }
175: }
176: return(0);
177: }
179: static PetscErrorCode NEPNLEIGSFNSingularities(FN f,PetscInt *nisol,PetscScalar **isol,PetscBool *rational)
180: {
182: FNCombineType ctype;
183: FN f1,f2;
184: PetscInt i,nq,nisol1,nisol2;
185: PetscScalar *qcoeff,*wr,*wi,*isol1,*isol2;
186: PetscBool flg,avail,rat1,rat2;
189: *rational = PETSC_FALSE;
190: PetscObjectTypeCompare((PetscObject)f,FNRATIONAL,&flg);
191: if (flg) {
192: *rational = PETSC_TRUE;
193: FNRationalGetDenominator(f,&nq,&qcoeff);
194: if (nq>1) {
195: PetscMalloc2(nq-1,&wr,nq-1,&wi);
196: NEPNLEIGSAuxiliarPRootFinder(nq-1,qcoeff,wr,wi,&avail);
197: if (avail) {
198: PetscCalloc1(nq-1,isol);
199: *nisol = 0;
200: for (i=0;i<nq-1;i++)
201: #if !defined(PETSC_USE_COMPLEX)
202: if (wi[i]==0)
203: #endif
204: (*isol)[(*nisol)++] = wr[i];
205: nq = *nisol; *nisol = 0;
206: for (i=0;i<nq;i++) wr[i] = (*isol)[i];
207: NEPNLEIGSAuxiliarRmDuplicates(nq,wr,nisol,*isol);
208: PetscFree2(wr,wi);
209: } else { *nisol=0; *isol = NULL; }
210: } else { *nisol = 0; *isol = NULL; }
211: PetscFree(qcoeff);
212: }
213: PetscObjectTypeCompare((PetscObject)f,FNCOMBINE,&flg);
214: if (flg) {
215: FNCombineGetChildren(f,&ctype,&f1,&f2);
216: if (ctype != FN_COMBINE_COMPOSE && ctype != FN_COMBINE_DIVIDE) {
217: NEPNLEIGSFNSingularities(f1,&nisol1,&isol1,&rat1);
218: NEPNLEIGSFNSingularities(f2,&nisol2,&isol2,&rat2);
219: if (nisol1+nisol2>0) {
220: PetscCalloc1(nisol1+nisol2,isol);
221: *nisol = 0;
222: NEPNLEIGSAuxiliarRmDuplicates(nisol1,isol1,nisol,*isol);
223: NEPNLEIGSAuxiliarRmDuplicates(nisol2,isol2,nisol,*isol);
224: }
225: *rational = (rat1&&rat2)?PETSC_TRUE:PETSC_FALSE;
226: PetscFree(isol1);
227: PetscFree(isol2);
228: }
229: }
230: return(0);
231: }
233: static PetscErrorCode NEPNLEIGSRationalSingularities(NEP nep,PetscInt *ndptx,PetscScalar *dxi,PetscBool *rational)
234: {
236: PetscInt nt,i,nisol;
237: FN f;
238: PetscScalar *isol;
239: PetscBool rat;
242: *rational = PETSC_TRUE;
243: *ndptx = 0;
244: NEPGetSplitOperatorInfo(nep,&nt,NULL);
245: for (i=0;i<nt;i++) {
246: NEPGetSplitOperatorTerm(nep,i,NULL,&f);
247: NEPNLEIGSFNSingularities(f,&nisol,&isol,&rat);
248: if (nisol) {
249: NEPNLEIGSAuxiliarRmDuplicates(nisol,isol,ndptx,dxi);
250: PetscFree(isol);
251: }
252: *rational = ((*rational)&&rat)?PETSC_TRUE:PETSC_FALSE;
253: }
254: return(0);
255: }
257: static PetscErrorCode NEPNLEIGSLejaBagbyPoints(NEP nep)
258: {
260: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
261: PetscInt i,k,ndpt=NDPOINTS,ndptx=NDPOINTS;
262: PetscScalar *ds,*dsi,*dxi,*nrs,*nrxi,*s=ctx->s,*xi=ctx->xi,*beta=ctx->beta;
263: PetscReal maxnrs,minnrxi;
264: PetscBool rational;
265: #if !defined(PETSC_USE_COMPLEX)
266: PetscReal a,b,h;
267: #endif
270: PetscMalloc5(ndpt+1,&ds,ndpt+1,&dsi,ndpt,&dxi,ndpt+1,&nrs,ndpt,&nrxi);
272: /* Discretize the target region boundary */
273: RGComputeContour(nep->rg,ndpt,ds,dsi);
274: #if !defined(PETSC_USE_COMPLEX)
275: for (i=0;i<ndpt;i++) if (dsi[i]!=0.0) break;
276: if (i<ndpt) {
277: if (nep->problem_type==NEP_RATIONAL) {
278: /* Select a segment in the real axis */
279: RGComputeBoundingBox(nep->rg,&a,&b,NULL,NULL);
280: if (a<=-PETSC_MAX_REAL || b>=PETSC_MAX_REAL) SETERRQ(PetscObjectComm((PetscObject)nep),1,"NLEIGS requires a bounded target set");
281: h = (b-a)/ndpt;
282: for (i=0;i<ndpt;i++) {ds[i] = a+h*i; dsi[i] = 0.0;}
283: } else SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"NLEIGS with real arithmetic requires the target set to be included in the real axis");
284: }
285: #endif
286: /* Discretize the singularity region */
287: if (ctx->computesingularities) {
288: (ctx->computesingularities)(nep,&ndptx,dxi,ctx->singularitiesctx);
289: } else {
290: if (nep->problem_type==NEP_RATIONAL) {
291: NEPNLEIGSRationalSingularities(nep,&ndptx,dxi,&rational);
292: if (!rational) SETERRQ(PetscObjectComm((PetscObject)nep),1,"Failed to determine singularities automatically in rational problem; consider solving the problem as general");
293: } else ndptx = 0;
294: }
296: /* Look for Leja-Bagby points in the discretization sets */
297: s[0] = ds[0];
298: xi[0] = (ndptx>0)?dxi[0]:PETSC_INFINITY;
299: 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]));
300: beta[0] = 1.0; /* scaling factors are also computed here */
301: for (i=0;i<ndpt;i++) {
302: nrs[i] = 1.0;
303: nrxi[i] = 1.0;
304: }
305: for (k=1;k<ctx->ddmaxit;k++) {
306: maxnrs = 0.0;
307: minnrxi = PETSC_MAX_REAL;
308: for (i=0;i<ndpt;i++) {
309: nrs[i] *= ((ds[i]-s[k-1])/(1.0-ds[i]/xi[k-1]))/beta[k-1];
310: if (PetscAbsScalar(nrs[i])>maxnrs) {maxnrs = PetscAbsScalar(nrs[i]); s[k] = ds[i];}
311: }
312: if (ndptx>k) {
313: for (i=1;i<ndptx;i++) {
314: nrxi[i] *= ((dxi[i]-s[k-1])/(1.0-dxi[i]/xi[k-1]))/beta[k-1];
315: if (PetscAbsScalar(nrxi[i])<minnrxi) {minnrxi = PetscAbsScalar(nrxi[i]); xi[k] = dxi[i];}
316: }
317: 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]));
318: } else xi[k] = PETSC_INFINITY;
319: beta[k] = maxnrs;
320: }
321: PetscFree5(ds,dsi,dxi,nrs,nrxi);
322: return(0);
323: }
325: static PetscErrorCode NEPNLEIGSEvalNRTFunct(NEP nep,PetscInt k,PetscScalar sigma,PetscScalar *b)
326: {
327: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
328: PetscInt i;
329: PetscScalar *beta=ctx->beta,*s=ctx->s,*xi=ctx->xi;
332: b[0] = 1.0/beta[0];
333: for (i=0;i<k;i++) {
334: b[i+1] = ((sigma-s[i])*b[i])/(beta[i+1]*(1.0-sigma/xi[i]));
335: }
336: return(0);
337: }
339: static PetscErrorCode MatMult_Fun(Mat A,Vec x,Vec y)
340: {
342: ShellMatCtx *ctx;
343: PetscInt i;
346: MatShellGetContext(A,(void**)&ctx);
347: MatMult(ctx->A[0],x,y);
348: if (ctx->coeff[0]!=1.0) { VecScale(y,ctx->coeff[0]); }
349: for (i=1;i<ctx->nmat;i++) {
350: MatMult(ctx->A[i],x,ctx->t);
351: VecAXPY(y,ctx->coeff[i],ctx->t);
352: }
353: return(0);
354: }
356: static PetscErrorCode MatMultTranspose_Fun(Mat A,Vec x,Vec y)
357: {
359: ShellMatCtx *ctx;
360: PetscInt i;
363: MatShellGetContext(A,(void**)&ctx);
364: MatMultTranspose(ctx->A[0],x,y);
365: if (ctx->coeff[0]!=1.0) { VecScale(y,ctx->coeff[0]); }
366: for (i=1;i<ctx->nmat;i++) {
367: MatMultTranspose(ctx->A[i],x,ctx->t);
368: VecAXPY(y,ctx->coeff[i],ctx->t);
369: }
370: return(0);
371: }
373: static PetscErrorCode MatGetDiagonal_Fun(Mat A,Vec diag)
374: {
376: ShellMatCtx *ctx;
377: PetscInt i;
380: MatShellGetContext(A,(void**)&ctx);
381: MatGetDiagonal(ctx->A[0],diag);
382: if (ctx->coeff[0]!=1.0) { VecScale(diag,ctx->coeff[0]); }
383: for (i=1;i<ctx->nmat;i++) {
384: MatGetDiagonal(ctx->A[i],ctx->t);
385: VecAXPY(diag,ctx->coeff[i],ctx->t);
386: }
387: return(0);
388: }
390: static PetscErrorCode MatDuplicate_Fun(Mat A,MatDuplicateOption op,Mat *B)
391: {
392: PetscInt n,i;
393: ShellMatCtx *ctxnew,*ctx;
394: void (*fun)();
398: MatShellGetContext(A,(void**)&ctx);
399: PetscNew(&ctxnew);
400: ctxnew->nmat = ctx->nmat;
401: ctxnew->maxnmat = ctx->maxnmat;
402: PetscMalloc2(ctxnew->maxnmat,&ctxnew->A,ctxnew->maxnmat,&ctxnew->coeff);
403: for (i=0;i<ctx->nmat;i++) {
404: PetscObjectReference((PetscObject)ctx->A[i]);
405: ctxnew->A[i] = ctx->A[i];
406: ctxnew->coeff[i] = ctx->coeff[i];
407: }
408: MatGetSize(ctx->A[0],&n,NULL);
409: VecDuplicate(ctx->t,&ctxnew->t);
410: MatCreateShell(PETSC_COMM_WORLD,n,n,n,n,(void*)ctxnew,B);
411: MatShellSetManageScalingShifts(*B);
412: MatShellGetOperation(A,MATOP_MULT,&fun);
413: MatShellSetOperation(*B,MATOP_MULT,fun);
414: MatShellGetOperation(A,MATOP_MULT_TRANSPOSE,&fun);
415: MatShellSetOperation(*B,MATOP_MULT_TRANSPOSE,fun);
416: MatShellGetOperation(A,MATOP_GET_DIAGONAL,&fun);
417: MatShellSetOperation(*B,MATOP_GET_DIAGONAL,fun);
418: MatShellGetOperation(A,MATOP_DUPLICATE,&fun);
419: MatShellSetOperation(*B,MATOP_DUPLICATE,fun);
420: MatShellGetOperation(A,MATOP_DESTROY,&fun);
421: MatShellSetOperation(*B,MATOP_DESTROY,fun);
422: MatShellGetOperation(A,MATOP_AXPY,&fun);
423: MatShellSetOperation(*B,MATOP_AXPY,fun);
424: return(0);
425: }
427: static PetscErrorCode MatDestroy_Fun(Mat A)
428: {
429: ShellMatCtx *ctx;
431: PetscInt i;
434: if (A) {
435: MatShellGetContext(A,(void**)&ctx);
436: for (i=0;i<ctx->nmat;i++) {
437: MatDestroy(&ctx->A[i]);
438: }
439: VecDestroy(&ctx->t);
440: PetscFree2(ctx->A,ctx->coeff);
441: PetscFree(ctx);
442: }
443: return(0);
444: }
446: static PetscErrorCode MatAXPY_Fun(Mat Y,PetscScalar a,Mat X,MatStructure str)
447: {
448: ShellMatCtx *ctxY,*ctxX;
450: PetscInt i,j;
451: PetscBool found;
454: MatShellGetContext(Y,(void**)&ctxY);
455: MatShellGetContext(X,(void**)&ctxX);
456: for (i=0;i<ctxX->nmat;i++) {
457: found = PETSC_FALSE;
458: for (j=0;!found&&j<ctxY->nmat;j++) {
459: if (ctxX->A[i]==ctxY->A[j]) {
460: found = PETSC_TRUE;
461: ctxY->coeff[j] += a*ctxX->coeff[i];
462: }
463: }
464: if (!found) {
465: ctxY->coeff[ctxY->nmat] = a*ctxX->coeff[i];
466: ctxY->A[ctxY->nmat++] = ctxX->A[i];
467: PetscObjectReference((PetscObject)ctxX->A[i]);
468: }
469: }
470: return(0);
471: }
473: static PetscErrorCode MatScale_Fun(Mat M,PetscScalar a)
474: {
475: ShellMatCtx *ctx;
477: PetscInt i;
480: MatShellGetContext(M,(void**)&ctx);
481: for (i=0;i<ctx->nmat;i++) ctx->coeff[i] *= a;
482: return(0);
483: }
485: static PetscErrorCode NLEIGSMatToMatShellArray(Mat M,Mat *Ms,PetscInt maxnmat)
486: {
488: ShellMatCtx *ctx;
489: PetscInt n;
490: PetscBool has;
493: MatHasOperation(M,MATOP_DUPLICATE,&has);
494: if (!has) SETERRQ(PetscObjectComm((PetscObject)M),1,"MatDuplicate operation required");
495: PetscNew(&ctx);
496: ctx->maxnmat = maxnmat;
497: PetscMalloc2(ctx->maxnmat,&ctx->A,ctx->maxnmat,&ctx->coeff);
498: MatDuplicate(M,MAT_COPY_VALUES,&ctx->A[0]);
499: ctx->nmat = 1;
500: ctx->coeff[0] = 1.0;
501: MatCreateVecs(M,&ctx->t,NULL);
502: MatGetSize(M,&n,NULL);
503: MatCreateShell(PetscObjectComm((PetscObject)M),n,n,n,n,(void*)ctx,Ms);
504: MatShellSetManageScalingShifts(*Ms);
505: MatShellSetOperation(*Ms,MATOP_MULT,(void(*)())MatMult_Fun);
506: MatShellSetOperation(*Ms,MATOP_MULT_TRANSPOSE,(void(*)())MatMultTranspose_Fun);
507: MatShellSetOperation(*Ms,MATOP_GET_DIAGONAL,(void(*)())MatGetDiagonal_Fun);
508: MatShellSetOperation(*Ms,MATOP_DUPLICATE,(void(*)())MatDuplicate_Fun);
509: MatShellSetOperation(*Ms,MATOP_DESTROY,(void(*)())MatDestroy_Fun);
510: MatShellSetOperation(*Ms,MATOP_AXPY,(void(*)())MatAXPY_Fun);
511: MatShellSetOperation(*Ms,MATOP_SCALE,(void(*)())MatScale_Fun);
512: return(0);
513: }
515: static PetscErrorCode NEPNLEIGSNormEstimation(NEP nep,Mat M,PetscReal *norm,Vec *w)
516: {
517: PetscScalar *z,*x,*y;
518: PetscReal tr;
519: Vec X=w[0],Y=w[1];
520: PetscInt n,i;
521: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
522: PetscRandom rand;
526: if (!ctx->vrn) {
527: /* generate a random vector with normally distributed entries with the Box-Muller transform */
528: BVGetRandomContext(nep->V,&rand);
529: MatCreateVecs(M,&ctx->vrn,NULL);
530: VecSetRandom(X,rand);
531: VecSetRandom(Y,rand);
532: VecGetLocalSize(ctx->vrn,&n);
533: VecGetArray(ctx->vrn,&z);
534: VecGetArray(X,&x);
535: VecGetArray(Y,&y);
536: for (i=0;i<n;i++) {
537: #if defined(PETSC_USE_COMPLEX)
538: z[i] = PetscSqrtReal(-2.0*PetscLogReal(PetscRealPart(x[i])))*PetscCosReal(2.0*PETSC_PI*PetscRealPart(y[i]));
539: z[i] += PETSC_i*(PetscSqrtReal(-2.0*PetscLogReal(PetscImaginaryPart(x[i])))*PetscCosReal(2.0*PETSC_PI*PetscImaginaryPart(y[i])));
540: #else
541: z[i] = PetscSqrtReal(-2.0*PetscLogReal(x[i]))*PetscCosReal(2.0*PETSC_PI*y[i]);
542: #endif
543: }
544: VecRestoreArray(ctx->vrn,&z);
545: VecRestoreArray(X,&x);
546: VecRestoreArray(Y,&y);
547: VecNorm(ctx->vrn,NORM_2,&tr);
548: VecScale(ctx->vrn,1/tr);
549: }
550: /* matrix-free norm estimator of Ipsen http://www4.ncsu.edu/~ipsen/ps/slides_ima.pdf */
551: MatGetSize(M,&n,NULL);
552: MatMult(M,ctx->vrn,X);
553: VecNorm(X,NORM_2,norm);
554: *norm *= PetscSqrtReal((PetscReal)n);
555: return(0);
556: }
558: static PetscErrorCode NEPNLEIGSDividedDifferences_split(NEP nep)
559: {
561: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
562: PetscInt k,j,i,maxnmat,nmax;
563: PetscReal norm0,norm,*matnorm;
564: PetscScalar *s=ctx->s,*beta=ctx->beta,*xi=ctx->xi,*b,alpha,*coeffs,*pK,*pH,sone=1.0;
565: Mat T,Ts,K,H;
566: PetscBool shell,hasmnorm=PETSC_FALSE,matrix=PETSC_TRUE;
567: PetscBLASInt n_;
570: nmax = ctx->ddmaxit;
571: PetscMalloc1(nep->nt*nmax,&ctx->coeffD);
572: PetscMalloc3(nmax+1,&b,nmax+1,&coeffs,nep->nt,&matnorm);
573: for (j=0;j<nep->nt;j++) {
574: MatHasOperation(nep->A[j],MATOP_NORM,&hasmnorm);
575: if (!hasmnorm) break;
576: MatNorm(nep->A[j],NORM_INFINITY,matnorm+j);
577: }
578: /* Try matrix functions scheme */
579: PetscCalloc2(nmax*nmax,&pK,nmax*nmax,&pH);
580: for (i=0;i<nmax-1;i++) {
581: pK[(nmax+1)*i] = 1.0;
582: pK[(nmax+1)*i+1] = beta[i+1]/xi[i];
583: pH[(nmax+1)*i] = s[i];
584: pH[(nmax+1)*i+1] = beta[i+1];
585: }
586: pH[nmax*nmax-1] = s[nmax-1];
587: pK[nmax*nmax-1] = 1.0;
588: PetscBLASIntCast(nmax,&n_);
589: PetscStackCallBLAS("BLAStrsm",BLAStrsm_("R","L","N","U",&n_,&n_,&sone,pK,&n_,pH,&n_));
590: /* The matrix to be used is in H. K will be a work-space matrix */
591: MatCreateSeqDense(PETSC_COMM_SELF,nmax,nmax,pH,&H);
592: MatCreateSeqDense(PETSC_COMM_SELF,nmax,nmax,pK,&K);
593: for (j=0;matrix&&j<nep->nt;j++) {
594: PetscPushErrorHandler(PetscIgnoreErrorHandler,NULL);
595: FNEvaluateFunctionMat(nep->f[j],H,K);
596: PetscPopErrorHandler();
597: if (!ierr) {
598: for (i=0;i<nmax;i++) { ctx->coeffD[j+i*nep->nt] = pK[i]*beta[0]; }
599: } else {
600: matrix = PETSC_FALSE;
601: PetscFPTrapPop();
602: }
603: }
604: MatDestroy(&H);
605: MatDestroy(&K);
606: if (!matrix) {
607: for (j=0;j<nep->nt;j++) {
608: FNEvaluateFunction(nep->f[j],s[0],ctx->coeffD+j);
609: ctx->coeffD[j] *= beta[0];
610: }
611: }
612: if (hasmnorm) {
613: norm0 = 0.0;
614: for (j=0;j<nep->nt;j++) norm0 += matnorm[j]*PetscAbsScalar(ctx->coeffD[j]);
615: } else {
616: norm0 = 0.0;
617: for (j=0;j<nep->nt;j++) norm0 = PetscMax(PetscAbsScalar(ctx->coeffD[j]),norm0);
618: }
619: ctx->nmat = ctx->ddmaxit;
620: for (k=1;k<ctx->ddmaxit;k++) {
621: if (!matrix) {
622: NEPNLEIGSEvalNRTFunct(nep,k,s[k],b);
623: for (i=0;i<nep->nt;i++) {
624: FNEvaluateFunction(nep->f[i],s[k],ctx->coeffD+k*nep->nt+i);
625: for (j=0;j<k;j++) {
626: ctx->coeffD[k*nep->nt+i] -= b[j]*ctx->coeffD[i+nep->nt*j];
627: }
628: ctx->coeffD[k*nep->nt+i] /= b[k];
629: }
630: }
631: if (hasmnorm) {
632: norm = 0.0;
633: for (j=0;j<nep->nt;j++) norm += matnorm[j]*PetscAbsScalar(ctx->coeffD[k*nep->nt+j]);
634: } else {
635: norm = 0.0;
636: for (j=0;j<nep->nt;j++) norm = PetscMax(PetscAbsScalar(ctx->coeffD[k*nep->nt+j]),norm);
637: }
638: if (norm/norm0 < ctx->ddtol) {
639: ctx->nmat = k+1;
640: break;
641: }
642: }
643: if (!ctx->ksp) { NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp); }
644: PetscObjectTypeCompare((PetscObject)nep->A[0],MATSHELL,&shell);
645: maxnmat = PetscMax(ctx->ddmaxit,nep->nt);
646: for (i=0;i<ctx->nshiftsw;i++) {
647: NEPNLEIGSEvalNRTFunct(nep,ctx->nmat-1,ctx->shifts[i],coeffs);
648: if (!shell) {
649: MatDuplicate(nep->A[0],MAT_COPY_VALUES,&T);
650: } else {
651: NLEIGSMatToMatShellArray(nep->A[0],&T,maxnmat);
652: }
653: alpha = 0.0;
654: for (j=0;j<ctx->nmat;j++) alpha += coeffs[j]*ctx->coeffD[j*nep->nt];
655: MatScale(T,alpha);
656: for (k=1;k<nep->nt;k++) {
657: alpha = 0.0;
658: for (j=0;j<ctx->nmat;j++) alpha += coeffs[j]*ctx->coeffD[j*nep->nt+k];
659: if (shell) {
660: NLEIGSMatToMatShellArray(nep->A[k],&Ts,maxnmat);
661: }
662: MatAXPY(T,alpha,shell?Ts:nep->A[k],nep->mstr);
663: if (shell) {
664: MatDestroy(&Ts);
665: }
666: }
667: KSPSetOperators(ctx->ksp[i],T,T);
668: KSPSetUp(ctx->ksp[i]);
669: MatDestroy(&T);
670: }
671: PetscFree3(b,coeffs,matnorm);
672: PetscFree2(pK,pH);
673: return(0);
674: }
676: static PetscErrorCode NEPNLEIGSDividedDifferences_callback(NEP nep)
677: {
679: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
680: PetscInt k,j,i,maxnmat;
681: PetscReal norm0,norm;
682: PetscScalar *s=ctx->s,*beta=ctx->beta,*b,*coeffs;
683: Mat *D=ctx->D,T;
684: PetscBool shell,has,vec=PETSC_FALSE;
685: Vec w[2];
688: PetscMalloc2(ctx->ddmaxit+1,&b,ctx->ddmaxit+1,&coeffs);
689: T = nep->function;
690: NEPComputeFunction(nep,s[0],T,T);
691: PetscObjectTypeCompare((PetscObject)T,MATSHELL,&shell);
692: maxnmat = PetscMax(ctx->ddmaxit,nep->nt);
693: if (!shell) {
694: MatDuplicate(T,MAT_COPY_VALUES,&D[0]);
695: } else {
696: NLEIGSMatToMatShellArray(T,&D[0],maxnmat);
697: }
698: if (beta[0]!=1.0) {
699: MatScale(D[0],1.0/beta[0]);
700: }
701: MatHasOperation(D[0],MATOP_NORM,&has);
702: if (has) {
703: MatNorm(D[0],NORM_FROBENIUS,&norm0);
704: } else {
705: MatCreateVecs(D[0],NULL,&w[0]);
706: VecDuplicate(w[0],&w[1]);
707: vec = PETSC_TRUE;
708: NEPNLEIGSNormEstimation(nep,D[0],&norm0,w);
709: }
710: ctx->nmat = ctx->ddmaxit;
711: for (k=1;k<ctx->ddmaxit;k++) {
712: NEPNLEIGSEvalNRTFunct(nep,k,s[k],b);
713: NEPComputeFunction(nep,s[k],T,T);
714: if (!shell) {
715: MatDuplicate(T,MAT_COPY_VALUES,&D[k]);
716: } else {
717: NLEIGSMatToMatShellArray(T,&D[k],maxnmat);
718: }
719: for (j=0;j<k;j++) {
720: MatAXPY(D[k],-b[j],D[j],nep->mstr);
721: }
722: MatScale(D[k],1.0/b[k]);
723: MatHasOperation(D[k],MATOP_NORM,&has);
724: if (has) {
725: MatNorm(D[k],NORM_FROBENIUS,&norm);
726: } else {
727: if(!vec) {
728: MatCreateVecs(D[k],NULL,&w[0]);
729: VecDuplicate(w[0],&w[1]);
730: vec = PETSC_TRUE;
731: }
732: NEPNLEIGSNormEstimation(nep,D[k],&norm,w);
733: }
734: if (norm/norm0 < ctx->ddtol && k>1) {
735: ctx->nmat = k+1;
736: break;
737: }
738: }
739: if (!ctx->ksp) { NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp); }
740: for (i=0;i<ctx->nshiftsw;i++) {
741: NEPNLEIGSEvalNRTFunct(nep,ctx->nmat-1,ctx->shifts[i],coeffs);
742: MatDuplicate(ctx->D[0],MAT_COPY_VALUES,&T);
743: if (coeffs[0]!=1.0) { MatScale(T,coeffs[0]); }
744: for (j=1;j<ctx->nmat;j++) {
745: MatAXPY(T,coeffs[j],ctx->D[j],nep->mstr);
746: }
747: KSPSetOperators(ctx->ksp[i],T,T);
748: KSPSetUp(ctx->ksp[i]);
749: MatDestroy(&T);
750: }
751: PetscFree2(b,coeffs);
752: if (vec) {
753: VecDestroy(&w[0]);
754: VecDestroy(&w[1]);
755: }
756: return(0);
757: }
759: /*
760: NEPKrylovConvergence - This is the analogue to EPSKrylovConvergence.
761: */
762: static PetscErrorCode NEPNLEIGSKrylovConvergence(NEP nep,PetscBool getall,PetscInt kini,PetscInt nits,PetscScalar betak,PetscReal betah,PetscInt *kout,Vec *w)
763: {
765: PetscInt k,newk,marker,inside;
766: PetscScalar re,im;
767: PetscReal resnorm,tt;
768: PetscBool istrivial;
769: NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
772: RGIsTrivial(nep->rg,&istrivial);
773: marker = -1;
774: if (nep->trackall) getall = PETSC_TRUE;
775: for (k=kini;k<kini+nits;k++) {
776: /* eigenvalue */
777: re = nep->eigr[k];
778: im = nep->eigi[k];
779: if (!istrivial) {
780: if (!ctx->nshifts) {
781: NEPNLEIGSBackTransform((PetscObject)nep,1,&re,&im);
782: }
783: RGCheckInside(nep->rg,1,&re,&im,&inside);
784: if (marker==-1 && inside<0) marker = k;
785: }
786: newk = k;
787: DSVectors(nep->ds,DS_MAT_X,&newk,&resnorm);
788: tt = (ctx->nshifts)?SlepcAbsEigenvalue(betak-nep->eigr[k]*betah,nep->eigi[k]*betah):betah;
789: resnorm *= PetscAbsReal(tt);
790: /* error estimate */
791: (*nep->converged)(nep,nep->eigr[k],nep->eigi[k],resnorm,&nep->errest[k],nep->convergedctx);
792: if (marker==-1 && nep->errest[k] >= nep->tol) marker = k;
793: if (newk==k+1) {
794: nep->errest[k+1] = nep->errest[k];
795: k++;
796: }
797: if (marker!=-1 && !getall) break;
798: }
799: if (marker!=-1) k = marker;
800: *kout = k;
801: return(0);
802: }
804: PetscErrorCode NEPSetUp_NLEIGS(NEP nep)
805: {
807: PetscInt k,in;
808: PetscScalar zero=0.0;
809: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
810: SlepcSC sc;
811: PetscBool istrivial;
814: NEPSetDimensions_Default(nep,nep->nev,&nep->ncv,&nep->mpd);
815: if (nep->ncv>nep->nev+nep->mpd) SETERRQ(PetscObjectComm((PetscObject)nep),1,"The value of ncv must not be larger than nev+mpd");
816: if (!nep->max_it) nep->max_it = PetscMax(5000,2*nep->n/nep->ncv);
817: if (!ctx->ddmaxit) ctx->ddmaxit = LBPOINTS;
818: RGIsTrivial(nep->rg,&istrivial);
819: if (istrivial) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"NEPNLEIGS requires a nontrivial region defining the target set");
820: if (!nep->which) nep->which = NEP_TARGET_MAGNITUDE;
822: /* Initialize the NLEIGS context structure */
823: k = ctx->ddmaxit;
824: PetscMalloc4(k,&ctx->s,k,&ctx->xi,k,&ctx->beta,k,&ctx->D);
825: nep->data = ctx;
826: if (nep->tol==PETSC_DEFAULT) nep->tol = SLEPC_DEFAULT_TOL;
827: if (ctx->ddtol==PETSC_DEFAULT) ctx->ddtol = nep->tol/10.0;
828: if (!ctx->keep) ctx->keep = 0.5;
830: /* Compute Leja-Bagby points and scaling values */
831: NEPNLEIGSLejaBagbyPoints(nep);
832: if (nep->problem_type!=NEP_RATIONAL) {
833: RGCheckInside(nep->rg,1,&nep->target,&zero,&in);
834: if (in<0) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"The target is not inside the target set");
835: }
837: /* Compute the divided difference matrices */
838: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
839: NEPNLEIGSDividedDifferences_split(nep);
840: } else {
841: NEPNLEIGSDividedDifferences_callback(nep);
842: }
843: NEPAllocateSolution(nep,ctx->nmat-1);
844: NEPSetWorkVecs(nep,4);
846: /* set-up DS and transfer split operator functions */
847: DSSetType(nep->ds,ctx->nshifts?DSGNHEP:DSNHEP);
848: DSAllocate(nep->ds,nep->ncv+1);
849: DSGetSlepcSC(nep->ds,&sc);
850: if (!ctx->nshifts) {
851: sc->map = NEPNLEIGSBackTransform;
852: DSSetExtraRow(nep->ds,PETSC_TRUE);
853: }
854: sc->mapobj = (PetscObject)nep;
855: sc->rg = nep->rg;
856: sc->comparison = nep->sc->comparison;
857: sc->comparisonctx = nep->sc->comparisonctx;
858: BVDestroy(&ctx->V);
859: BVCreateTensor(nep->V,ctx->nmat-1,&ctx->V);
860: return(0);
861: }
863: /*
864: Extend the TOAR basis by applying the the matrix operator
865: over a vector which is decomposed on the TOAR way
866: Input:
867: - S,V: define the latest Arnoldi vector (nv vectors in V)
868: Output:
869: - t: new vector extending the TOAR basis
870: - r: temporally coefficients to compute the TOAR coefficients
871: for the new Arnoldi vector
872: Workspace: t_ (two vectors)
873: */
874: 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_)
875: {
877: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
878: PetscInt deg=ctx->nmat-1,k,j;
879: Vec v=t_[0],q=t_[1],w;
880: PetscScalar *beta=ctx->beta,*s=ctx->s,*xi=ctx->xi,*coeffs,sigma;
883: if (!ctx->ksp) { NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp); }
884: sigma = ctx->shifts[idxrktg];
885: BVSetActiveColumns(nep->V,0,nv);
886: PetscMalloc1(ctx->nmat,&coeffs);
887: if (PetscAbsScalar(s[deg-2]-sigma)<100*PETSC_MACHINE_EPSILON) SETERRQ(PETSC_COMM_SELF,1,"Breakdown in NLEIGS");
888: /* i-part stored in (i-1) position */
889: for (j=0;j<nv;j++) {
890: 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);
891: }
892: BVSetActiveColumns(W,0,deg);
893: BVGetColumn(W,deg-1,&w);
894: BVMultVec(V,1.0/beta[deg],0,w,S+(deg-1)*ls);
895: BVRestoreColumn(W,deg-1,&w);
896: BVGetColumn(W,deg-2,&w);
897: BVMultVec(V,1.0,0.0,w,r+(deg-2)*lr);
898: BVRestoreColumn(W,deg-2,&w);
899: for (k=deg-2;k>0;k--) {
900: if (PetscAbsScalar(s[k-1]-sigma)<100*PETSC_MACHINE_EPSILON) SETERRQ(PETSC_COMM_SELF,1,"Breakdown in NLEIGS");
901: 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);
902: BVGetColumn(W,k-1,&w);
903: BVMultVec(V,1.0,0.0,w,r+(k-1)*lr);
904: BVRestoreColumn(W,k-1,&w);
905: }
906: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
907: for (j=0;j<ctx->nmat-2;j++) coeffs[j] = ctx->coeffD[nep->nt*j];
908: coeffs[ctx->nmat-2] = ctx->coeffD[nep->nt*(ctx->nmat-1)];
909: BVMultVec(W,1.0,0.0,v,coeffs);
910: MatMult(nep->A[0],v,q);
911: for (k=1;k<nep->nt;k++) {
912: for (j=0;j<ctx->nmat-2;j++) coeffs[j] = ctx->coeffD[nep->nt*j+k];
913: coeffs[ctx->nmat-2] = ctx->coeffD[nep->nt*(ctx->nmat-1)+k];
914: BVMultVec(W,1.0,0,v,coeffs);
915: MatMult(nep->A[k],v,t);
916: VecAXPY(q,1.0,t);
917: }
918: KSPSolve(ctx->ksp[idxrktg],q,t);
919: VecScale(t,-1.0);
920: } else {
921: for (k=0;k<deg-1;k++) {
922: BVGetColumn(W,k,&w);
923: MatMult(ctx->D[k],w,q);
924: BVRestoreColumn(W,k,&w);
925: BVInsertVec(W,k,q);
926: }
927: BVGetColumn(W,deg-1,&w);
928: MatMult(ctx->D[deg],w,q);
929: BVRestoreColumn(W,k,&w);
930: BVInsertVec(W,k,q);
931: for (j=0;j<ctx->nmat-1;j++) coeffs[j] = 1.0;
932: BVMultVec(W,1.0,0.0,q,coeffs);
933: KSPSolve(ctx->ksp[idxrktg],q,t);
934: VecScale(t,-1.0);
935: }
936: PetscFree(coeffs);
937: return(0);
938: }
940: /*
941: Compute TOAR coefficients of the blocks of the new Arnoldi vector computed
942: */
943: static PetscErrorCode NEPTOARCoefficients(NEP nep,PetscScalar sigma,PetscInt nv,PetscScalar *S,PetscInt ls,PetscScalar *r,PetscInt lr,PetscScalar *x,PetscScalar *work)
944: {
946: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
947: PetscInt k,j,d=ctx->nmat-1;
948: PetscScalar *t=work;
951: NEPNLEIGSEvalNRTFunct(nep,d-1,sigma,t);
952: for (k=0;k<d-1;k++) {
953: for (j=0;j<=nv;j++) r[k*lr+j] += t[k]*x[j];
954: }
955: for (j=0;j<=nv;j++) r[(d-1)*lr+j] = t[d-1]*x[j];
956: return(0);
957: }
959: /*
960: Compute continuation vector coefficients for the Rational-Krylov run.
961: dim(work) >= (end-ini)*(end-ini+1) + end+1 + 2*(end-ini+1), dim(t) = end.
962: */
963: 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)
964: {
965: #if defined(PETSC_MISSING_LAPACK_GEQRF) || defined(SLEPC_MISSING_LAPACK_LARF)
967: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GEQRF/LARF - Lapack routines are unavailable");
968: #else
970: PetscScalar *x,*W,*tau,sone=1.0,szero=0.0;
971: PetscInt i,j,n1,n,nwu=0;
972: PetscBLASInt info,n_,n1_,one=1,dim,lds_;
973: NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
976: if (!ctx->nshifts || !end) {
977: t[0] = 1;
978: PetscMemcpy(cont,S+end*lds,lds*sizeof(PetscScalar));
979: } else {
980: n = end-ini;
981: n1 = n+1;
982: x = work+nwu;
983: nwu += end+1;
984: tau = work+nwu;
985: nwu += n;
986: W = work+nwu;
987: nwu += n1*n;
988: for (j=ini;j<end;j++) {
989: for (i=ini;i<=end;i++) W[(j-ini)*n1+i-ini] = K[j*ld+i] -H[j*ld+i]*sigma;
990: }
991: PetscBLASIntCast(n,&n_);
992: PetscBLASIntCast(n1,&n1_);
993: PetscBLASIntCast(end+1,&dim);
994: PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&n1_,&n_,W,&n1_,tau,work+nwu,&n1_,&info));
995: SlepcCheckLapackInfo("geqrf",info);
996: for (i=0;i<end;i++) t[i] = 0.0;
997: t[end] = 1.0;
998: for (j=n-1;j>=0;j--) {
999: for (i=0;i<ini+j;i++) x[i] = 0.0;
1000: x[ini+j] = 1.0;
1001: for (i=j+1;i<n1;i++) x[i+ini] = W[i+n1*j];
1002: tau[j] = PetscConj(tau[j]);
1003: PetscStackCallBLAS("LAPACKlarf",LAPACKlarf_("L",&dim,&one,x,&one,tau+j,t,&dim,work+nwu));
1004: }
1005: PetscBLASIntCast(lds,&lds_);
1006: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&lds_,&n1_,&sone,S,&lds_,t,&one,&szero,cont,&one));
1007: }
1008: return(0);
1009: #endif
1010: }
1012: /*
1013: Compute a run of Arnoldi iterations
1014: */
1015: static PetscErrorCode NEPNLEIGSTOARrun(NEP nep,PetscScalar *K,PetscScalar *H,PetscInt ldh,BV W,PetscInt k,PetscInt *M,PetscBool *breakdown,Vec *t_)
1016: {
1018: NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
1019: PetscInt i,j,m=*M,lwa,deg=ctx->nmat-1,lds,nqt,ld,l;
1020: Vec t;
1021: PetscReal norm;
1022: PetscScalar *x,*work,*tt,sigma,*cont,*S;
1023: PetscBool lindep;
1024: Mat MS;
1027: BVTensorGetFactors(ctx->V,NULL,&MS);
1028: MatDenseGetArray(MS,&S);
1029: BVGetSizes(nep->V,NULL,NULL,&ld);
1030: lds = ld*deg;
1031: BVGetActiveColumns(nep->V,&l,&nqt);
1032: lwa = PetscMax(ld,deg)+(m+1)*(m+1)+4*(m+1);
1033: PetscMalloc4(ld,&x,lwa,&work,m+1,&tt,lds,&cont);
1034: BVSetActiveColumns(ctx->V,0,m);
1035: for (j=k;j<m;j++) {
1036: sigma = ctx->shifts[(++(ctx->idxrk))%ctx->nshiftsw];
1038: /* Continuation vector */
1039: NEPNLEIGS_RKcontinuation(nep,0,j,K,H,ldh,sigma,S,lds,cont,tt,work);
1041: /* apply operator */
1042: BVGetColumn(nep->V,nqt,&t);
1043: NEPTOARExtendBasis(nep,(ctx->idxrk)%ctx->nshiftsw,cont,ld,nqt,W,nep->V,t,S+(j+1)*lds,ld,t_);
1044: BVRestoreColumn(nep->V,nqt,&t);
1046: /* orthogonalize */
1047: BVOrthogonalizeColumn(nep->V,nqt,x,&norm,&lindep);
1048: if (!lindep) {
1049: x[nqt] = norm;
1050: BVScaleColumn(nep->V,nqt,1.0/norm);
1051: nqt++;
1052: } else x[nqt] = 0.0;
1054: NEPTOARCoefficients(nep,sigma,nqt-1,cont,ld,S+(j+1)*lds,ld,x,work);
1056: /* Level-2 orthogonalization */
1057: BVOrthogonalizeColumn(ctx->V,j+1,H+j*ldh,&norm,breakdown);
1058: H[j+1+ldh*j] = norm;
1059: if (ctx->nshifts) {
1060: for (i=0;i<=j;i++) K[i+ldh*j] = sigma*H[i+ldh*j] + tt[i];
1061: K[j+1+ldh*j] = sigma*H[j+1+ldh*j];
1062: }
1063: if (*breakdown) {
1064: *M = j+1;
1065: break;
1066: }
1067: BVScaleColumn(ctx->V,j+1,1.0/norm);
1068: BVSetActiveColumns(nep->V,l,nqt);
1069: }
1070: PetscFree4(x,work,tt,cont);
1071: MatDenseRestoreArray(MS,&S);
1072: BVTensorRestoreFactors(ctx->V,NULL,&MS);
1073: return(0);
1074: }
1076: PetscErrorCode NEPSolve_NLEIGS(NEP nep)
1077: {
1079: NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
1080: PetscInt i,k=0,l,nv=0,ld,lds,ldds,nq,newn;
1081: PetscInt deg=ctx->nmat-1,nconv=0;
1082: PetscScalar *S,*Q,*H,*pU,*K,betak=0,*eigr,*eigi;
1083: PetscReal betah;
1084: PetscBool falselock=PETSC_FALSE,breakdown=PETSC_FALSE;
1085: BV W;
1086: Mat MS,MQ,U;
1089: if (ctx->lock) {
1090: PetscOptionsGetBool(NULL,NULL,"-nep_nleigs_falselocking",&falselock,NULL);
1091: }
1093: BVGetSizes(nep->V,NULL,NULL,&ld);
1094: lds = deg*ld;
1095: DSGetLeadingDimension(nep->ds,&ldds);
1096: if (!ctx->nshifts) {
1097: PetscMalloc2(nep->ncv,&eigr,nep->ncv,&eigi);
1098: } else { eigr = nep->eigr; eigi = nep->eigi; }
1099: BVDuplicateResize(nep->V,PetscMax(nep->nt-1,ctx->nmat-1),&W);
1102: /* clean projected matrix (including the extra-arrow) */
1103: DSGetArray(nep->ds,DS_MAT_A,&H);
1104: PetscMemzero(H,ldds*ldds*sizeof(PetscScalar));
1105: DSRestoreArray(nep->ds,DS_MAT_A,&H);
1106: if (ctx->nshifts) {
1107: DSGetArray(nep->ds,DS_MAT_B,&H);
1108: PetscMemzero(H,ldds*ldds*sizeof(PetscScalar));
1109: DSRestoreArray(nep->ds,DS_MAT_B,&H);
1110: }
1111:
1112: /* Get the starting Arnoldi vector */
1113: BVTensorBuildFirstColumn(ctx->V,nep->nini);
1114:
1115: /* Restart loop */
1116: l = 0;
1117: while (nep->reason == NEP_CONVERGED_ITERATING) {
1118: nep->its++;
1120: /* Compute an nv-step Krylov relation */
1121: nv = PetscMin(nep->nconv+nep->mpd,nep->ncv);
1122: if (ctx->nshifts) { DSGetArray(nep->ds,DS_MAT_A,&K); }
1123: DSGetArray(nep->ds,ctx->nshifts?DS_MAT_B:DS_MAT_A,&H);
1124: NEPNLEIGSTOARrun(nep,K,H,ldds,W,nep->nconv+l,&nv,&breakdown,nep->work);
1125: betah = PetscAbsScalar(H[(nv-1)*ldds+nv]);
1126: DSRestoreArray(nep->ds,ctx->nshifts?DS_MAT_B:DS_MAT_A,&H);
1127: if (ctx->nshifts) {
1128: betak = K[(nv-1)*ldds+nv];
1129: DSRestoreArray(nep->ds,DS_MAT_A,&K);
1130: }
1131: DSSetDimensions(nep->ds,nv,0,nep->nconv,nep->nconv+l);
1132: if (l==0) {
1133: DSSetState(nep->ds,DS_STATE_INTERMEDIATE);
1134: } else {
1135: DSSetState(nep->ds,DS_STATE_RAW);
1136: }
1138: /* Solve projected problem */
1139: DSSolve(nep->ds,nep->eigr,nep->eigi);
1140: DSSort(nep->ds,nep->eigr,nep->eigi,NULL,NULL,NULL);
1141: if (!ctx->nshifts) {
1142: DSUpdateExtraRow(nep->ds);
1143: }
1144: DSSynchronize(nep->ds,nep->eigr,nep->eigi);
1146: /* Check convergence */
1147: NEPNLEIGSKrylovConvergence(nep,PETSC_FALSE,nep->nconv,nv-nep->nconv,betak,betah,&k,nep->work);
1148: (*nep->stopping)(nep,nep->its,nep->max_it,k,nep->nev,&nep->reason,nep->stoppingctx);
1150: /* Update l */
1151: if (nep->reason != NEP_CONVERGED_ITERATING || breakdown) l = 0;
1152: else {
1153: l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
1154: if (!breakdown) {
1155: /* Prepare the Rayleigh quotient for restart */
1156: if (!ctx->nshifts) {
1157: DSTruncate(nep->ds,k+l);
1158: DSGetDimensions(nep->ds,&newn,NULL,NULL,NULL,NULL);
1159: l = newn-k;
1160: } else {
1161: DSGetArray(nep->ds,DS_MAT_Q,&Q);
1162: DSGetArray(nep->ds,DS_MAT_B,&H);
1163: DSGetArray(nep->ds,DS_MAT_A,&K);
1164: for (i=ctx->lock?k:0;i<k+l;i++) {
1165: H[k+l+i*ldds] = betah*Q[nv-1+i*ldds];
1166: K[k+l+i*ldds] = betak*Q[nv-1+i*ldds];
1167: }
1168: DSRestoreArray(nep->ds,DS_MAT_B,&H);
1169: DSRestoreArray(nep->ds,DS_MAT_A,&K);
1170: DSRestoreArray(nep->ds,DS_MAT_Q,&Q);
1171: }
1172: }
1173: }
1174: nconv = k;
1175: if (!ctx->lock && nep->reason == NEP_CONVERGED_ITERATING && !breakdown) { l += k; k = 0; }
1177: /* Update S */
1178: DSGetMat(nep->ds,ctx->nshifts?DS_MAT_Z:DS_MAT_Q,&MQ);
1179: BVMultInPlace(ctx->V,MQ,nep->nconv,k+l);
1180: MatDestroy(&MQ);
1182: /* Copy last column of S */
1183: BVCopyColumn(ctx->V,nv,k+l);
1185: if (breakdown && nep->reason == NEP_CONVERGED_ITERATING) {
1186: /* Stop if breakdown */
1187: PetscInfo2(nep,"Breakdown (it=%D norm=%g)\n",nep->its,(double)betah);
1188: nep->reason = NEP_DIVERGED_BREAKDOWN;
1189: }
1190: if (nep->reason != NEP_CONVERGED_ITERATING) l--;
1191: /* truncate S */
1192: BVGetActiveColumns(nep->V,NULL,&nq);
1193: if (k+l+deg<=nq) {
1194: BVSetActiveColumns(ctx->V,nep->nconv,k+l+1);
1195: if (!falselock && ctx->lock) {
1196: BVTensorCompress(ctx->V,k-nep->nconv);
1197: } else {
1198: BVTensorCompress(ctx->V,0);
1199: }
1200: }
1201: nep->nconv = k;
1202: if (!ctx->nshifts) {
1203: for (i=0;i<nv;i++) { eigr[i] = nep->eigr[i]; eigi[i] = nep->eigi[i]; }
1204: NEPNLEIGSBackTransform((PetscObject)nep,nv,eigr,eigi);
1205: }
1206: NEPMonitor(nep,nep->its,nconv,eigr,eigi,nep->errest,nv);
1207: }
1208: nep->nconv = nconv;
1209: if (nep->nconv>0) {
1210: BVSetActiveColumns(ctx->V,0,nep->nconv);
1211: BVGetActiveColumns(nep->V,NULL,&nq);
1212: BVSetActiveColumns(nep->V,0,nq);
1213: if (nq>nep->nconv) {
1214: BVTensorCompress(ctx->V,nep->nconv);
1215: BVSetActiveColumns(nep->V,0,nep->nconv);
1216: nq = nep->nconv;
1217: }
1218: if (ctx->nshifts) {
1219: DSGetMat(nep->ds,DS_MAT_B,&MQ);
1220: BVMultInPlace(ctx->V,MQ,0,nep->nconv);
1221: MatDestroy(&MQ);
1222: }
1223: BVTensorGetFactors(ctx->V,NULL,&MS);
1224: MatDenseGetArray(MS,&S);
1225: PetscMalloc1(nq*nep->nconv,&pU);
1226: for (i=0;i<nep->nconv;i++) {
1227: PetscMemcpy(pU+i*nq,S+i*lds,nq*sizeof(PetscScalar));
1228: }
1229: MatDenseRestoreArray(MS,&S);
1230: BVTensorRestoreFactors(ctx->V,NULL,&MS);
1231: MatCreateSeqDense(PETSC_COMM_SELF,nq,nep->nconv,pU,&U);
1232: BVSetActiveColumns(nep->V,0,nq);
1233: BVMultInPlace(nep->V,U,0,nep->nconv);
1234: BVSetActiveColumns(nep->V,0,nep->nconv);
1235: MatDestroy(&U);
1236: PetscFree(pU);
1237: /* truncate Schur decomposition and change the state to raw so that
1238: DSVectors() computes eigenvectors from scratch */
1239: DSSetDimensions(nep->ds,nep->nconv,0,0,0);
1240: DSSetState(nep->ds,DS_STATE_RAW);
1241: }
1243: /* Map eigenvalues back to the original problem */
1244: if (!ctx->nshifts) {
1245: NEPNLEIGSBackTransform((PetscObject)nep,nep->nconv,nep->eigr,nep->eigi);
1246: PetscFree2(eigr,eigi);
1247: }
1248: BVDestroy(&W);
1249: return(0);
1250: }
1252: static PetscErrorCode NEPNLEIGSSetSingularitiesFunction_NLEIGS(NEP nep,PetscErrorCode (*fun)(NEP,PetscInt*,PetscScalar*,void*),void *ctx)
1253: {
1254: NEP_NLEIGS *nepctx=(NEP_NLEIGS*)nep->data;
1257: if (fun) nepctx->computesingularities = fun;
1258: if (ctx) nepctx->singularitiesctx = ctx;
1259: nep->state = NEP_STATE_INITIAL;
1260: return(0);
1261: }
1263: /*@C
1264: NEPNLEIGSSetSingularitiesFunction - Sets a user function to compute a discretization
1265: of the singularity set (where T(.) is not analytic).
1267: Logically Collective on NEP
1269: Input Parameters:
1270: + nep - the NEP context
1271: . fun - user function (if NULL then NEP retains any previously set value)
1272: - ctx - [optional] user-defined context for private data for the function
1273: (may be NULL, in which case NEP retains any previously set value)
1275: Calling Sequence of fun:
1276: $ fun(NEP nep,PetscInt *maxnp,PetscScalar *xi,void *ctx)
1278: + nep - the NEP context
1279: . maxnp - on input number of requested points in the discretization (can be set)
1280: . xi - computed values of the discretization
1281: - ctx - optional context, as set by NEPNLEIGSSetSingularitiesFunction()
1283: Notes:
1284: The user-defined function can set a smaller value of maxnp if necessary.
1285: It is wrong to return a larger value.
1287: If the problem type has been set to rational with NEPSetProblemType(),
1288: then it is not necessary to set the singularities explicitly since the
1289: solver will try to determine them automatically.
1291: Level: intermediate
1293: .seealso: NEPNLEIGSGetSingularitiesFunction(), NEPSetProblemType()
1294: @*/
1295: PetscErrorCode NEPNLEIGSSetSingularitiesFunction(NEP nep,PetscErrorCode (*fun)(NEP,PetscInt*,PetscScalar*,void*),void *ctx)
1296: {
1301: PetscTryMethod(nep,"NEPNLEIGSSetSingularitiesFunction_C",(NEP,PetscErrorCode(*)(NEP,PetscInt*,PetscScalar*,void*),void*),(nep,fun,ctx));
1302: return(0);
1303: }
1305: static PetscErrorCode NEPNLEIGSGetSingularitiesFunction_NLEIGS(NEP nep,PetscErrorCode (**fun)(NEP,PetscInt*,PetscScalar*,void*),void **ctx)
1306: {
1307: NEP_NLEIGS *nepctx=(NEP_NLEIGS*)nep->data;
1310: if (fun) *fun = nepctx->computesingularities;
1311: if (ctx) *ctx = nepctx->singularitiesctx;
1312: return(0);
1313: }
1315: /*@C
1316: NEPNLEIGSGetSingularitiesFunction - Returns the Function and optionally the user
1317: provided context for computing a discretization of the singularity set.
1319: Not Collective
1321: Input Parameter:
1322: . nep - the nonlinear eigensolver context
1324: Output Parameters:
1325: + fun - location to put the function (or NULL)
1326: - ctx - location to stash the function context (or NULL)
1328: Level: advanced
1330: .seealso: NEPNLEIGSSetSingularitiesFunction()
1331: @*/
1332: PetscErrorCode NEPNLEIGSGetSingularitiesFunction(NEP nep,PetscErrorCode (**fun)(NEP,PetscInt*,PetscScalar*,void*),void **ctx)
1333: {
1338: PetscUseMethod(nep,"NEPNLEIGSGetSingularitiesFunction_C",(NEP,PetscErrorCode(**)(NEP,PetscInt*,PetscScalar*,void*),void**),(nep,fun,ctx));
1339: return(0);
1340: }
1342: static PetscErrorCode NEPNLEIGSSetRestart_NLEIGS(NEP nep,PetscReal keep)
1343: {
1344: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1347: if (keep==PETSC_DEFAULT) ctx->keep = 0.5;
1348: else {
1349: 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]");
1350: ctx->keep = keep;
1351: }
1352: return(0);
1353: }
1355: /*@
1356: NEPNLEIGSSetRestart - Sets the restart parameter for the NLEIGS
1357: method, in particular the proportion of basis vectors that must be kept
1358: after restart.
1360: Logically Collective on NEP
1362: Input Parameters:
1363: + nep - the nonlinear eigensolver context
1364: - keep - the number of vectors to be kept at restart
1366: Options Database Key:
1367: . -nep_nleigs_restart - Sets the restart parameter
1369: Notes:
1370: Allowed values are in the range [0.1,0.9]. The default is 0.5.
1372: Level: advanced
1374: .seealso: NEPNLEIGSGetRestart()
1375: @*/
1376: PetscErrorCode NEPNLEIGSSetRestart(NEP nep,PetscReal keep)
1377: {
1383: PetscTryMethod(nep,"NEPNLEIGSSetRestart_C",(NEP,PetscReal),(nep,keep));
1384: return(0);
1385: }
1387: static PetscErrorCode NEPNLEIGSGetRestart_NLEIGS(NEP nep,PetscReal *keep)
1388: {
1389: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1392: *keep = ctx->keep;
1393: return(0);
1394: }
1396: /*@
1397: NEPNLEIGSGetRestart - Gets the restart parameter used in the NLEIGS method.
1399: Not Collective
1401: Input Parameter:
1402: . nep - the nonlinear eigensolver context
1404: Output Parameter:
1405: . keep - the restart parameter
1407: Level: advanced
1409: .seealso: NEPNLEIGSSetRestart()
1410: @*/
1411: PetscErrorCode NEPNLEIGSGetRestart(NEP nep,PetscReal *keep)
1412: {
1418: PetscUseMethod(nep,"NEPNLEIGSGetRestart_C",(NEP,PetscReal*),(nep,keep));
1419: return(0);
1420: }
1422: static PetscErrorCode NEPNLEIGSSetLocking_NLEIGS(NEP nep,PetscBool lock)
1423: {
1424: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1427: ctx->lock = lock;
1428: return(0);
1429: }
1431: /*@
1432: NEPNLEIGSSetLocking - Choose between locking and non-locking variants of
1433: the NLEIGS method.
1435: Logically Collective on NEP
1437: Input Parameters:
1438: + nep - the nonlinear eigensolver context
1439: - lock - true if the locking variant must be selected
1441: Options Database Key:
1442: . -nep_nleigs_locking - Sets the locking flag
1444: Notes:
1445: The default is to lock converged eigenpairs when the method restarts.
1446: This behaviour can be changed so that all directions are kept in the
1447: working subspace even if already converged to working accuracy (the
1448: non-locking variant).
1450: Level: advanced
1452: .seealso: NEPNLEIGSGetLocking()
1453: @*/
1454: PetscErrorCode NEPNLEIGSSetLocking(NEP nep,PetscBool lock)
1455: {
1461: PetscTryMethod(nep,"NEPNLEIGSSetLocking_C",(NEP,PetscBool),(nep,lock));
1462: return(0);
1463: }
1465: static PetscErrorCode NEPNLEIGSGetLocking_NLEIGS(NEP nep,PetscBool *lock)
1466: {
1467: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1470: *lock = ctx->lock;
1471: return(0);
1472: }
1474: /*@
1475: NEPNLEIGSGetLocking - Gets the locking flag used in the NLEIGS method.
1477: Not Collective
1479: Input Parameter:
1480: . nep - the nonlinear eigensolver context
1482: Output Parameter:
1483: . lock - the locking flag
1485: Level: advanced
1487: .seealso: NEPNLEIGSSetLocking()
1488: @*/
1489: PetscErrorCode NEPNLEIGSGetLocking(NEP nep,PetscBool *lock)
1490: {
1496: PetscUseMethod(nep,"NEPNLEIGSGetLocking_C",(NEP,PetscBool*),(nep,lock));
1497: return(0);
1498: }
1500: static PetscErrorCode NEPNLEIGSSetInterpolation_NLEIGS(NEP nep,PetscReal tol,PetscInt degree)
1501: {
1503: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1506: if (tol == PETSC_DEFAULT) {
1507: ctx->ddtol = PETSC_DEFAULT;
1508: nep->state = NEP_STATE_INITIAL;
1509: } else {
1510: if (tol <= 0.0) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of tol. Must be > 0");
1511: ctx->ddtol = tol;
1512: }
1513: if (degree == PETSC_DEFAULT || degree == PETSC_DECIDE) {
1514: ctx->ddmaxit = 0;
1515: if (nep->state) { NEPReset(nep); }
1516: nep->state = NEP_STATE_INITIAL;
1517: } else {
1518: if (degree <= 0) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of degree. Must be > 0");
1519: if (ctx->ddmaxit != degree) {
1520: ctx->ddmaxit = degree;
1521: if (nep->state) { NEPReset(nep); }
1522: nep->state = NEP_STATE_INITIAL;
1523: }
1524: }
1525: return(0);
1526: }
1528: /*@
1529: NEPNLEIGSSetInterpolation - Sets the tolerance and maximum degree
1530: when building the interpolation via divided differences.
1532: Logically Collective on NEP
1534: Input Parameters:
1535: + nep - the nonlinear eigensolver context
1536: . tol - tolerance to stop computing divided differences
1537: - degree - maximum degree of interpolation
1539: Options Database Key:
1540: + -nep_nleigs_interpolation_tol <tol> - Sets the tolerance to stop computing divided differences
1541: - -nep_nleigs_interpolation_degree <degree> - Sets the maximum degree of interpolation
1543: Notes:
1544: Use PETSC_DEFAULT for either argument to assign a reasonably good value.
1546: Level: advanced
1548: .seealso: NEPNLEIGSGetInterpolation()
1549: @*/
1550: PetscErrorCode NEPNLEIGSSetInterpolation(NEP nep,PetscReal tol,PetscInt degree)
1551: {
1558: PetscTryMethod(nep,"NEPNLEIGSSetInterpolation_C",(NEP,PetscReal,PetscInt),(nep,tol,degree));
1559: return(0);
1560: }
1562: static PetscErrorCode NEPNLEIGSGetInterpolation_NLEIGS(NEP nep,PetscReal *tol,PetscInt *degree)
1563: {
1564: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1567: if (tol) *tol = ctx->ddtol;
1568: if (degree) *degree = ctx->ddmaxit;
1569: return(0);
1570: }
1572: /*@
1573: NEPNLEIGSGetInterpolation - Gets the tolerance and maximum degree
1574: when building the interpolation via divided differences.
1576: Not Collective
1578: Input Parameter:
1579: . nep - the nonlinear eigensolver context
1581: Output Parameter:
1582: + tol - tolerance to stop computing divided differences
1583: - degree - maximum degree of interpolation
1585: Level: advanced
1587: .seealso: NEPNLEIGSSetInterpolation()
1588: @*/
1589: PetscErrorCode NEPNLEIGSGetInterpolation(NEP nep,PetscReal *tol,PetscInt *degree)
1590: {
1595: PetscTryMethod(nep,"NEPNLEIGSGetInterpolation_C",(NEP,PetscReal*,PetscInt*),(nep,tol,degree));
1596: return(0);
1597: }
1599: static PetscErrorCode NEPNLEIGSSetRKShifts_NLEIGS(NEP nep,PetscInt ns,PetscScalar *shifts)
1600: {
1602: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1603: PetscInt i;
1606: if (ns<=0) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"Number of shifts must be positive");
1607: if (ctx->nshifts) { PetscFree(ctx->shifts); }
1608: for (i=0;i<ctx->nshiftsw;i++) { KSPDestroy(&ctx->ksp[i]); }
1609: PetscFree(ctx->ksp);
1610: ctx->ksp = NULL;
1611: PetscMalloc1(ns,&ctx->shifts);
1612: for (i=0;i<ns;i++) ctx->shifts[i] = shifts[i];
1613: ctx->nshifts = ns;
1614: nep->state = NEP_STATE_INITIAL;
1615: return(0);
1616: }
1618: /*@C
1619: NEPNLEIGSSetRKShifts - Sets a list of shifts to be used in the Rational
1620: Krylov method.
1622: Logically Collective on NEP
1624: Input Parameters:
1625: + nep - the nonlinear eigensolver context
1626: . ns - number of shifts
1627: - shifts - array of scalar values specifying the shifts
1629: Options Database Key:
1630: . -nep_nleigs_rk_shifts - Sets the list of shifts
1632: Notes:
1633: If only one shift is provided, the built subspace built is equivalent to
1634: shift-and-invert Krylov-Schur (provided that the absolute convergence
1635: criterion is used).
1637: In the case of real scalars, complex shifts are not allowed. In the
1638: command line, a comma-separated list of complex values can be provided with
1639: the format [+/-][realnumber][+/-]realnumberi with no spaces, e.g.
1640: -nep_nleigs_rk_shifts 1.0+2.0i,1.5+2.0i,1.0+1.5i
1642: Level: advanced
1644: .seealso: NEPNLEIGSGetRKShifts()
1645: @*/
1646: PetscErrorCode NEPNLEIGSSetRKShifts(NEP nep,PetscInt ns,PetscScalar *shifts)
1647: {
1654: PetscTryMethod(nep,"NEPNLEIGSSetRKShifts_C",(NEP,PetscInt,PetscScalar*),(nep,ns,shifts));
1655: return(0);
1656: }
1658: static PetscErrorCode NEPNLEIGSGetRKShifts_NLEIGS(NEP nep,PetscInt *ns,PetscScalar **shifts)
1659: {
1661: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1662: PetscInt i;
1665: *ns = ctx->nshifts;
1666: if (ctx->nshifts) {
1667: PetscMalloc1(ctx->nshifts,shifts);
1668: for (i=0;i<ctx->nshifts;i++) (*shifts)[i] = ctx->shifts[i];
1669: }
1670: return(0);
1671: }
1673: /*@C
1674: NEPNLEIGSGetRKShifts - Gets the list of shifts used in the Rational
1675: Krylov method.
1677: Not Collective
1679: Input Parameter:
1680: . nep - the nonlinear eigensolver context
1682: Output Parameter:
1683: + ns - number of shifts
1684: - shifts - array of shifts
1686: Level: advanced
1688: .seealso: NEPNLEIGSSetRKShifts()
1689: @*/
1690: PetscErrorCode NEPNLEIGSGetRKShifts(NEP nep,PetscInt *ns,PetscScalar **shifts)
1691: {
1698: PetscTryMethod(nep,"NEPNLEIGSGetRKShifts_C",(NEP,PetscInt*,PetscScalar**),(nep,ns,shifts));
1699: return(0);
1700: }
1702: static PetscErrorCode NEPNLEIGSGetKSPs_NLEIGS(NEP nep,PetscInt *nsolve,KSP **ksp)
1703: {
1705: NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
1706: PetscInt i;
1707: PC pc;
1710: if (!ctx->ksp) {
1711: NEPNLEIGSSetShifts(nep,&ctx->nshiftsw);
1712: PetscMalloc1(ctx->nshiftsw,&ctx->ksp);
1713: for (i=0;i<ctx->nshiftsw;i++) {
1714: KSPCreate(PetscObjectComm((PetscObject)nep),&ctx->ksp[i]);
1715: PetscObjectIncrementTabLevel((PetscObject)ctx->ksp[i],(PetscObject)nep,1);
1716: KSPSetOptionsPrefix(ctx->ksp[i],((PetscObject)nep)->prefix);
1717: KSPAppendOptionsPrefix(ctx->ksp[i],"nep_nleigs_");
1718: PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->ksp[i]);
1719: KSPSetErrorIfNotConverged(ctx->ksp[i],PETSC_TRUE);
1720: KSPSetTolerances(ctx->ksp[i],SLEPC_DEFAULT_TOL,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
1721: KSPGetPC(ctx->ksp[i],&pc);
1722: KSPSetType(ctx->ksp[i],KSPPREONLY);
1723: PCSetType(pc,PCLU);
1724: }
1725: }
1726: if (nsolve) *nsolve = ctx->nshiftsw;
1727: if (ksp) *ksp = ctx->ksp;
1728: return(0);
1729: }
1731: /*@C
1732: NEPNLEIGSGetKSPs - Retrieve the array of linear solver objects associated with
1733: the nonlinear eigenvalue solver.
1735: Not Collective
1737: Input Parameter:
1738: . nep - nonlinear eigenvalue solver
1740: Output Parameters:
1741: + nsolve - number of returned KSP objects
1742: - ksp - array of linear solver object
1744: Notes:
1745: The number of KSP objects is equal to the number of shifts provided by the user,
1746: or 1 if the user did not provide shifts.
1748: Level: advanced
1750: .seealso: NEPNLEIGSSetRKShifts()
1751: @*/
1752: PetscErrorCode NEPNLEIGSGetKSPs(NEP nep,PetscInt *nsolve,KSP **ksp)
1753: {
1758: PetscUseMethod(nep,"NEPNLEIGSGetKSPs_C",(NEP,PetscInt*,KSP**),(nep,nsolve,ksp));
1759: return(0);
1760: }
1762: #define SHIFTMAX 30
1764: PetscErrorCode NEPSetFromOptions_NLEIGS(PetscOptionItems *PetscOptionsObject,NEP nep)
1765: {
1767: NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
1768: PetscInt i=0,k;
1769: PetscBool flg1,flg2,b;
1770: PetscReal r;
1771: PetscScalar array[SHIFTMAX];
1774: PetscOptionsHead(PetscOptionsObject,"NEP NLEIGS Options");
1776: PetscOptionsReal("-nep_nleigs_restart","Proportion of vectors kept after restart","NEPNLEIGSSetRestart",0.5,&r,&flg1);
1777: if (flg1) { NEPNLEIGSSetRestart(nep,r); }
1779: PetscOptionsBool("-nep_nleigs_locking","Choose between locking and non-locking variants","NEPNLEIGSSetLocking",PETSC_FALSE,&b,&flg1);
1780: if (flg1) { NEPNLEIGSSetLocking(nep,b); }
1782: NEPNLEIGSGetInterpolation(nep,&r,&i);
1783: if (!i) i = PETSC_DEFAULT;
1784: PetscOptionsInt("-nep_nleigs_interpolation_degree","Maximum number of terms for interpolation via divided differences","NEPNLEIGSSetInterpolation",i,&i,&flg1);
1785: PetscOptionsReal("-nep_nleigs_interpolation_tol","Tolerance for interpolation via divided differences","NEPNLEIGSSetInterpolation",r,&r,&flg2);
1786: if (flg1 || flg2) { NEPNLEIGSSetInterpolation(nep,r,i); }
1788: k = SHIFTMAX;
1789: for (i=0;i<k;i++) array[i] = 0;
1790: PetscOptionsScalarArray("-nep_nleigs_rk_shifts","Shifts for Rational Krylov","NEPNLEIGSSetRKShifts",array,&k,&flg1);
1791: if (flg1) { NEPNLEIGSSetRKShifts(nep,k,array); }
1793: PetscOptionsTail();
1795: if (!ctx->ksp) { NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp); }
1796: for (i=0;i<ctx->nshiftsw;i++) {
1797: KSPSetFromOptions(ctx->ksp[i]);
1798: }
1799: return(0);
1800: }
1802: PetscErrorCode NEPView_NLEIGS(NEP nep,PetscViewer viewer)
1803: {
1805: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1806: PetscBool isascii;
1807: PetscInt i;
1808: char str[50];
1811: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
1812: if (isascii) {
1813: PetscViewerASCIIPrintf(viewer," %d%% of basis vectors kept after restart\n",(int)(100*ctx->keep));
1814: PetscViewerASCIIPrintf(viewer," using the %slocking variant\n",ctx->lock?"":"non-");
1815: PetscViewerASCIIPrintf(viewer," divided difference terms: used=%D, max=%D\n",ctx->nmat,ctx->ddmaxit);
1816: PetscViewerASCIIPrintf(viewer," tolerance for divided difference convergence: %g\n",(double)ctx->ddtol);
1817: if (ctx->nshifts) {
1818: PetscViewerASCIIPrintf(viewer," RK shifts: ");
1819: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1820: for (i=0;i<ctx->nshifts;i++) {
1821: SlepcSNPrintfScalar(str,50,ctx->shifts[i],PETSC_FALSE);
1822: PetscViewerASCIIPrintf(viewer,"%s%s",str,(i<ctx->nshifts-1)?",":"");
1823: }
1824: PetscViewerASCIIPrintf(viewer,"\n");
1825: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1826: }
1827: if (!ctx->ksp) { NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp); }
1828: KSPView(ctx->ksp[0],viewer);
1829: }
1830: return(0);
1831: }
1833: PetscErrorCode NEPReset_NLEIGS(NEP nep)
1834: {
1836: PetscInt k;
1837: NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;
1840: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
1841: PetscFree(ctx->coeffD);
1842: } else {
1843: for (k=0;k<ctx->nmat;k++) { MatDestroy(&ctx->D[k]); }
1844: }
1845: PetscFree4(ctx->s,ctx->xi,ctx->beta,ctx->D);
1846: for (k=0;k<ctx->nshiftsw;k++) { KSPReset(ctx->ksp[k]); }
1847: if (ctx->vrn) {
1848: VecDestroy(&ctx->vrn);
1849: }
1850: return(0);
1851: }
1853: PetscErrorCode NEPDestroy_NLEIGS(NEP nep)
1854: {
1856: PetscInt k;
1857: NEP_NLEIGS *ctx = (NEP_NLEIGS*)nep->data;
1860: BVDestroy(&ctx->V);
1861: for (k=0;k<ctx->nshiftsw;k++) { KSPDestroy(&ctx->ksp[k]); }
1862: PetscFree(ctx->ksp);
1863: if (ctx->nshifts) { PetscFree(ctx->shifts); }
1864: PetscFree(nep->data);
1865: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetSingularitiesFunction_C",NULL);
1866: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetSingularitiesFunction_C",NULL);
1867: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRestart_C",NULL);
1868: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRestart_C",NULL);
1869: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetLocking_C",NULL);
1870: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetLocking_C",NULL);
1871: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetInterpolation_C",NULL);
1872: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetInterpolation_C",NULL);
1873: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRKShifts_C",NULL);
1874: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRKShifts_C",NULL);
1875: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetKSPs_C",NULL);
1876: return(0);
1877: }
1879: PETSC_EXTERN PetscErrorCode NEPCreate_NLEIGS(NEP nep)
1880: {
1882: NEP_NLEIGS *ctx;
1885: PetscNewLog(nep,&ctx);
1886: nep->data = (void*)ctx;
1887: ctx->lock = PETSC_TRUE;
1888: ctx->ddtol = PETSC_DEFAULT;
1890: nep->ops->solve = NEPSolve_NLEIGS;
1891: nep->ops->setup = NEPSetUp_NLEIGS;
1892: nep->ops->setfromoptions = NEPSetFromOptions_NLEIGS;
1893: nep->ops->view = NEPView_NLEIGS;
1894: nep->ops->destroy = NEPDestroy_NLEIGS;
1895: nep->ops->reset = NEPReset_NLEIGS;
1896: nep->ops->computevectors = NEPComputeVectors_Schur;
1898: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetSingularitiesFunction_C",NEPNLEIGSSetSingularitiesFunction_NLEIGS);
1899: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetSingularitiesFunction_C",NEPNLEIGSGetSingularitiesFunction_NLEIGS);
1900: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRestart_C",NEPNLEIGSSetRestart_NLEIGS);
1901: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRestart_C",NEPNLEIGSGetRestart_NLEIGS);
1902: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetLocking_C",NEPNLEIGSSetLocking_NLEIGS);
1903: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetLocking_C",NEPNLEIGSGetLocking_NLEIGS);
1904: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetInterpolation_C",NEPNLEIGSSetInterpolation_NLEIGS);
1905: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetInterpolation_C",NEPNLEIGSGetInterpolation_NLEIGS);
1906: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRKShifts_C",NEPNLEIGSSetRKShifts_NLEIGS);
1907: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRKShifts_C",NEPNLEIGSGetRKShifts_NLEIGS);
1908: PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetKSPs_C",NEPNLEIGSGetKSPs_NLEIGS);
1909: return(0);
1910: }