Actual source code: nleigs.c

slepc-3.15.2 2021-09-20
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2021, Universitat Politecnica de Valencia, Spain

  6:    This file is part of SLEPc.
  7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
  8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  9: */
 10: /*
 11:    SLEPc nonlinear eigensolver: "nleigs"

 13:    Method: NLEIGS

 15:    Algorithm:

 17:        Fully rational Krylov method for nonlinear eigenvalue problems.

 19:    References:

 21:        [1] S. Guttel et al., "NLEIGS: A class of robust fully rational Krylov
 22:            method for nonlinear eigenvalue problems", SIAM J. Sci. Comput.
 23:            36(6):A2842-A2864, 2014.
 24: */

 26: #include <slepc/private/nepimpl.h>
 27: #include <slepcblaslapack.h>
 28: #include "nleigs.h"

 30: PetscErrorCode NEPNLEIGSBackTransform(PetscObject ob,PetscInt n,PetscScalar *valr,PetscScalar *vali)
 31: {
 32:   NEP         nep;
 33:   PetscInt    j;
 34: #if !defined(PETSC_USE_COMPLEX)
 35:   PetscScalar t;
 36: #endif

 39:   nep = (NEP)ob;
 40: #if !defined(PETSC_USE_COMPLEX)
 41:   for (j=0;j<n;j++) {
 42:     if (vali[j] == 0) valr[j] = 1.0 / valr[j] + nep->target;
 43:     else {
 44:       t = valr[j] * valr[j] + vali[j] * vali[j];
 45:       valr[j] = valr[j] / t + nep->target;
 46:       vali[j] = - vali[j] / t;
 47:     }
 48:   }
 49: #else
 50:   for (j=0;j<n;j++) {
 51:     valr[j] = 1.0 / valr[j] + nep->target;
 52:   }
 53: #endif
 54:   return(0);
 55: }

 57: /* Computes the roots of a polynomial */
 58: static PetscErrorCode NEPNLEIGSAuxiliarPRootFinder(PetscInt deg,PetscScalar *polcoeffs,PetscScalar *wr,PetscScalar *wi,PetscBool *avail)
 59: {
 61:   PetscScalar    *C;
 62:   PetscBLASInt   n_,lwork;
 63:   PetscInt       i;
 64: #if defined(PETSC_USE_COMPLEX) || defined(PETSC_HAVE_ESSL)
 65:   PetscReal      *rwork=NULL;
 66: #endif
 67: #if defined(PETSC_HAVE_ESSL)
 68:   PetscScalar    sdummy,*wri;
 69:   PetscBLASInt   idummy,io=0;
 70: #else
 71:   PetscScalar    *work;
 72:   PetscBLASInt   info;
 73: #endif

 76:   *avail = PETSC_TRUE;
 77:   if (deg>0) {
 78:     PetscCalloc1(deg*deg,&C);
 79:     PetscBLASIntCast(deg,&n_);
 80:     for (i=0;i<deg-1;i++) {
 81:       C[(deg+1)*i+1]   = 1.0;
 82:       C[(deg-1)*deg+i] = -polcoeffs[deg-i]/polcoeffs[0];
 83:     }
 84:     C[deg*deg+-1] = -polcoeffs[1]/polcoeffs[0];
 85:     PetscBLASIntCast(3*deg,&lwork);

 87:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
 88: #if !defined(PETSC_HAVE_ESSL)
 89: #if !defined(PETSC_USE_COMPLEX)
 90:     PetscMalloc1(lwork,&work);
 91:     PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_("N","N",&n_,C,&n_,wr,wi,NULL,&n_,NULL,&n_,work,&lwork,&info));
 92:     if (info) *avail = PETSC_FALSE;
 93:     PetscFree(work);
 94: #else
 95:     PetscMalloc2(2*deg,&rwork,lwork,&work);
 96:     PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_("N","N",&n_,C,&n_,wr,NULL,&n_,NULL,&n_,work,&lwork,rwork,&info));
 97:     if (info) *avail = PETSC_FALSE;
 98:     PetscFree2(rwork,work);
 99: #endif
100: #else /* defined(PETSC_HAVE_ESSL) */
101:     PetscMalloc2(lwork,&rwork,2*deg,&wri);
102:     PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_(&io,C,&n_,wri,&sdummy,&idummy,&idummy,&n_,rwork,&lwork));
103: #if !defined(PETSC_USE_COMPLEX)
104:     for (i=0;i<deg;i++) {
105:       wr[i] = wri[2*i];
106:       wi[i] = wri[2*i+1];
107:     }
108: #else
109:     for (i=0;i<deg;i++) wr[i] = wri[i];
110: #endif
111:     PetscFree2(rwork,wri);
112: #endif
113:     PetscFPTrapPop();
114:     PetscFree(C);
115:   }
116:   return(0);
117: }

119: static PetscErrorCode NEPNLEIGSAuxiliarRmDuplicates(PetscInt nin,PetscScalar *pin,PetscInt *nout,PetscScalar *pout,PetscInt max)
120: {
121:   PetscInt i,j;

124:   for (i=0;i<nin;i++) {
125:     if (max && *nout>=max) break;
126:     pout[(*nout)++] = pin[i];
127:     for (j=0;j<*nout-1;j++)
128:       if (PetscAbsScalar(pin[i]-pout[j])<PETSC_MACHINE_EPSILON*100) {
129:         (*nout)--;
130:         break;
131:       }
132:   }
133:   return(0);
134: }

136: static PetscErrorCode NEPNLEIGSFNSingularities(FN f,PetscInt *nisol,PetscScalar **isol,PetscBool *rational)
137: {
139:   FNCombineType  ctype;
140:   FN             f1,f2;
141:   PetscInt       i,nq,nisol1,nisol2;
142:   PetscScalar    *qcoeff,*wr,*wi,*isol1,*isol2;
143:   PetscBool      flg,avail,rat1,rat2;

146:   *rational = PETSC_FALSE;
147:   PetscObjectTypeCompare((PetscObject)f,FNRATIONAL,&flg);
148:   if (flg) {
149:     *rational = PETSC_TRUE;
150:     FNRationalGetDenominator(f,&nq,&qcoeff);
151:     if (nq>1) {
152:       PetscMalloc2(nq-1,&wr,nq-1,&wi);
153:       NEPNLEIGSAuxiliarPRootFinder(nq-1,qcoeff,wr,wi,&avail);
154:       if (avail) {
155:         PetscCalloc1(nq-1,isol);
156:         *nisol = 0;
157:         for (i=0;i<nq-1;i++)
158: #if !defined(PETSC_USE_COMPLEX)
159:           if (wi[i]==0)
160: #endif
161:             (*isol)[(*nisol)++] = wr[i];
162:         nq = *nisol; *nisol = 0;
163:         for (i=0;i<nq;i++) wr[i] = (*isol)[i];
164:         NEPNLEIGSAuxiliarRmDuplicates(nq,wr,nisol,*isol,0);
165:         PetscFree2(wr,wi);
166:       } else { *nisol=0; *isol = NULL; }
167:     } else { *nisol = 0; *isol = NULL; }
168:     PetscFree(qcoeff);
169:   }
170:   PetscObjectTypeCompare((PetscObject)f,FNCOMBINE,&flg);
171:   if (flg) {
172:     FNCombineGetChildren(f,&ctype,&f1,&f2);
173:     if (ctype != FN_COMBINE_COMPOSE && ctype != FN_COMBINE_DIVIDE) {
174:       NEPNLEIGSFNSingularities(f1,&nisol1,&isol1,&rat1);
175:       NEPNLEIGSFNSingularities(f2,&nisol2,&isol2,&rat2);
176:       if (nisol1+nisol2>0) {
177:         PetscCalloc1(nisol1+nisol2,isol);
178:         *nisol = 0;
179:         NEPNLEIGSAuxiliarRmDuplicates(nisol1,isol1,nisol,*isol,0);
180:         NEPNLEIGSAuxiliarRmDuplicates(nisol2,isol2,nisol,*isol,0);
181:       }
182:       *rational = (rat1&&rat2)?PETSC_TRUE:PETSC_FALSE;
183:       PetscFree(isol1);
184:       PetscFree(isol2);
185:     }
186:   }
187:   return(0);
188: }

190: static PetscErrorCode NEPNLEIGSRationalSingularities(NEP nep,PetscInt *ndptx,PetscScalar *dxi,PetscBool *rational)
191: {
193:   PetscInt       nt,i,nisol;
194:   FN             f;
195:   PetscScalar    *isol;
196:   PetscBool      rat;

199:   *rational = PETSC_TRUE;
200:   *ndptx = 0;
201:   NEPGetSplitOperatorInfo(nep,&nt,NULL);
202:   for (i=0;i<nt;i++) {
203:     NEPGetSplitOperatorTerm(nep,i,NULL,&f);
204:     NEPNLEIGSFNSingularities(f,&nisol,&isol,&rat);
205:     if (nisol) {
206:       NEPNLEIGSAuxiliarRmDuplicates(nisol,isol,ndptx,dxi,0);
207:       PetscFree(isol);
208:     }
209:     *rational = ((*rational)&&rat)?PETSC_TRUE:PETSC_FALSE;
210:   }
211:   return(0);
212: }

214: /*  Adaptive Anderson-Antoulas algorithm */
215: static PetscErrorCode NEPNLEIGSAAAComputation(NEP nep,PetscInt ndpt,PetscScalar *ds,PetscScalar *F,PetscInt *ndptx,PetscScalar *dxi)
216: {
218:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
219:   PetscScalar    mean=0.0,*z,*f,*C,*A,*VT,*work,*ww,szero=0.0,sone=1.0;
220:   PetscScalar    *N,*D;
221:   PetscReal      *S,norm,err,*R;
222:   PetscInt       i,k,j,idx=0,cont;
223:   PetscBLASInt   n_,m_,lda_,lwork,info,one=1;
224: #if defined(PETSC_USE_COMPLEX)
225:   PetscReal      *rwork;
226: #endif

229:   PetscBLASIntCast(8*ndpt,&lwork);
230:   PetscMalloc5(ndpt,&R,ndpt,&z,ndpt,&f,ndpt*ndpt,&C,ndpt,&ww);
231:   PetscMalloc6(ndpt*ndpt,&A,ndpt,&S,ndpt*ndpt,&VT,lwork,&work,ndpt,&D,ndpt,&N);
232: #if defined(PETSC_USE_COMPLEX)
233:   PetscMalloc1(8*ndpt,&rwork);
234: #endif
235:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
236:   norm = 0.0;
237:   for (i=0;i<ndpt;i++) {
238:     mean += F[i];
239:     norm = PetscMax(PetscAbsScalar(F[i]),norm);
240:   }
241:   mean /= ndpt;
242:   PetscBLASIntCast(ndpt,&lda_);
243:   for (i=0;i<ndpt;i++) R[i] = PetscAbsScalar(F[i]-mean);
244:   /* next support point */
245:   err = 0.0;
246:   for (i=0;i<ndpt;i++) if (R[i]>=err) {idx = i; err = R[i];}
247:   for (k=0;k<ndpt-1;k++) {
248:     z[k] = ds[idx]; f[k] = F[idx]; R[idx] = -1.0;
249:     /* next column of Cauchy matrix */
250:     for (i=0;i<ndpt;i++) {
251:       C[i+k*ndpt] = 1.0/(ds[i]-ds[idx]);
252:     }

254:     PetscArrayzero(A,ndpt*ndpt);
255:     cont = 0;
256:     for (i=0;i<ndpt;i++) {
257:       if (R[i]!=-1.0) {
258:         for (j=0;j<=k;j++)A[cont+j*ndpt] = C[i+j*ndpt]*F[i]-C[i+j*ndpt]*f[j];
259:         cont++;
260:       }
261:     }
262:     PetscBLASIntCast(cont,&m_);
263:     PetscBLASIntCast(k+1,&n_);
264: #if defined(PETSC_USE_COMPLEX)
265:     PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","A",&m_,&n_,A,&lda_,S,NULL,&lda_,VT,&lda_,work,&lwork,rwork,&info));
266: #else
267:     PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","A",&m_,&n_,A,&lda_,S,NULL,&lda_,VT,&lda_,work,&lwork,&info));
268: #endif
269:     SlepcCheckLapackInfo("gesvd",info);
270:     for (i=0;i<=k;i++) {
271:       ww[i] = PetscConj(VT[i*ndpt+k]);
272:       D[i] = ww[i]*f[i];
273:     }
274:     PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&lda_,&n_,&sone,C,&lda_,D,&one,&szero,N,&one));
275:     PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&lda_,&n_,&sone,C,&lda_,ww,&one,&szero,D,&one));
276:     for (i=0;i<ndpt;i++) if (R[i]>=0) R[i] = PetscAbsScalar(F[i]-N[i]/D[i]);
277:     /* next support point */
278:     err = 0.0;
279:     for (i=0;i<ndpt;i++) if (R[i]>=err) {idx = i; err = R[i];}
280:     if (err <= ctx->ddtol*norm) break;
281:   }

283:   if (k==ndpt-1) SETERRQ(PetscObjectComm((PetscObject)nep),1,"Failed to determine singularities automatically in general problem");
284:   /* poles */
285:   PetscArrayzero(C,ndpt*ndpt);
286:   PetscArrayzero(A,ndpt*ndpt);
287:   for (i=0;i<=k;i++) {
288:     C[i+ndpt*i] = 1.0;
289:     A[(i+1)*ndpt] = ww[i];
290:     A[i+1] = 1.0;
291:     A[i+1+(i+1)*ndpt] = z[i];
292:   }
293:   C[0] = 0.0; C[k+1+(k+1)*ndpt] = 1.0;
294:   n_++;
295: #if defined(PETSC_USE_COMPLEX)
296:   PetscStackCallBLAS("LAPACKggev",LAPACKggev_("N","N",&n_,A,&lda_,C,&lda_,D,N,NULL,&lda_,NULL,&lda_,work,&lwork,rwork,&info));
297: #else
298:   PetscStackCallBLAS("LAPACKggev",LAPACKggev_("N","N",&n_,A,&lda_,C,&lda_,D,VT,N,NULL,&lda_,NULL,&lda_,work,&lwork,&info));
299: #endif
300:   SlepcCheckLapackInfo("ggev",info);
301:   cont = 0.0;
302:   for (i=0;i<n_;i++) if (N[i]!=0.0) {
303:     dxi[cont++] = D[i]/N[i];
304:   }
305:   *ndptx = cont;
306:   PetscFPTrapPop();
307:   PetscFree5(R,z,f,C,ww);
308:   PetscFree6(A,S,VT,work,D,N);
309: #if defined(PETSC_USE_COMPLEX)
310:   PetscFree(rwork);
311: #endif
312:   return(0);
313: }

315: /*  Singularities using Adaptive Anderson-Antoulas algorithm */
316: static PetscErrorCode NEPNLEIGSAAASingularities(NEP nep,PetscInt ndpt,PetscScalar *ds,PetscInt *ndptx,PetscScalar *dxi)
317: {
319:   Vec            u,v,w;
320:   PetscRandom    rand;
321:   PetscScalar    *F,*isol;
322:   PetscInt       i,k,nisol,nt;
323:   Mat            T;
324:   FN             f;

327:   PetscMalloc1(ndpt,&F);
328:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
329:     PetscMalloc1(ndpt,&isol);
330:     *ndptx = 0;
331:     NEPGetSplitOperatorInfo(nep,&nt,NULL);
332:     nisol = *ndptx;
333:     for (k=0;k<nt;k++) {
334:       NEPGetSplitOperatorTerm(nep,k,NULL,&f);
335:       for (i=0;i<ndpt;i++) {
336:         FNEvaluateFunction(f,ds[i],&F[i]);
337:       }
338:       NEPNLEIGSAAAComputation(nep,ndpt,ds,F,&nisol,isol);
339:       if (nisol) {
340:         NEPNLEIGSAuxiliarRmDuplicates(nisol,isol,ndptx,dxi,ndpt);
341:       }
342:     }
343:     PetscFree(isol);
344:   } else {
345:     MatCreateVecs(nep->function,&u,NULL);
346:     VecDuplicate(u,&v);
347:     VecDuplicate(u,&w);
348:     BVGetRandomContext(nep->V,&rand);
349:     VecSetRandom(u,rand);
350:     VecNormalize(u,NULL);
351:     VecSetRandom(v,rand);
352:     VecNormalize(v,NULL);
353:     T = nep->function;
354:     for (i=0;i<ndpt;i++) {
355:       NEPComputeFunction(nep,ds[i],T,T);
356:       MatMult(T,v,w);
357:       VecDot(w,u,&F[i]);
358:     }
359:     NEPNLEIGSAAAComputation(nep,ndpt,ds,F,ndptx,dxi);
360:     VecDestroy(&u);
361:     VecDestroy(&v);
362:     VecDestroy(&w);
363:   }
364:   PetscFree(F);
365:   return(0);
366: }

368: static PetscErrorCode NEPNLEIGSLejaBagbyPoints(NEP nep)
369: {
371:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
372:   PetscInt       i,k,ndpt=NDPOINTS,ndptx=NDPOINTS;
373:   PetscScalar    *ds,*dsi,*dxi,*nrs,*nrxi,*s=ctx->s,*xi=ctx->xi,*beta=ctx->beta;
374:   PetscReal      maxnrs,minnrxi;
375:   PetscBool      rational;
376: #if !defined(PETSC_USE_COMPLEX)
377:   PetscReal      a,b,h;
378: #endif

381:   if (!ctx->computesingularities && nep->problem_type!=NEP_RATIONAL) ndpt = ndptx = LBPOINTS;
382:   PetscMalloc5(ndpt+1,&ds,ndpt+1,&dsi,ndpt,&dxi,ndpt+1,&nrs,ndpt,&nrxi);

384:   /* Discretize the target region boundary */
385:   RGComputeContour(nep->rg,ndpt,ds,dsi);
386: #if !defined(PETSC_USE_COMPLEX)
387:   for (i=0;i<ndpt;i++) if (dsi[i]!=0.0) break;
388:   if (i<ndpt) {
389:     if (nep->problem_type==NEP_RATIONAL) {
390:       /* Select a segment in the real axis */
391:       RGComputeBoundingBox(nep->rg,&a,&b,NULL,NULL);
392:       if (a<=-PETSC_MAX_REAL || b>=PETSC_MAX_REAL) SETERRQ(PetscObjectComm((PetscObject)nep),1,"NLEIGS requires a bounded target set");
393:       h = (b-a)/ndpt;
394:       for (i=0;i<ndpt;i++) {ds[i] = a+h*i; dsi[i] = 0.0;}
395:     } else SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"NLEIGS with real arithmetic requires the target set to be included in the real axis");
396:   }
397: #endif
398:   /* Discretize the singularity region */
399:   if (ctx->computesingularities) {
400:     (ctx->computesingularities)(nep,&ndptx,dxi,ctx->singularitiesctx);
401:   } else {
402:     if (nep->problem_type==NEP_RATIONAL) {
403:       NEPNLEIGSRationalSingularities(nep,&ndptx,dxi,&rational);
404:       if (!rational) SETERRQ(PetscObjectComm((PetscObject)nep),1,"Failed to determine singularities automatically in rational problem; consider solving the problem as general");
405:     } else {
406:       /* AAA algorithm */
407:       NEPNLEIGSAAASingularities(nep,ndpt,ds,&ndptx,dxi);
408:     }
409:   }
410:   /* Look for Leja-Bagby points in the discretization sets */
411:   s[0]    = ds[0];
412:   xi[0]   = (ndptx>0)?dxi[0]:PETSC_INFINITY;
413:   if (PetscAbsScalar(xi[0])<10*PETSC_MACHINE_EPSILON) SETERRQ2(PetscObjectComm((PetscObject)nep),1,"Singularity point %D is nearly zero: %g; consider removing the singularity or shifting the problem",0,(double)PetscAbsScalar(xi[0]));
414:   beta[0] = 1.0; /* scaling factors are also computed here */
415:   for (i=0;i<ndpt;i++) {
416:     nrs[i] = 1.0;
417:     nrxi[i] = 1.0;
418:   }
419:   for (k=1;k<ctx->ddmaxit;k++) {
420:     maxnrs = 0.0;
421:     minnrxi = PETSC_MAX_REAL;
422:     for (i=0;i<ndpt;i++) {
423:       nrs[i] *= ((ds[i]-s[k-1])/(1.0-ds[i]/xi[k-1]))/beta[k-1];
424:       if (PetscAbsScalar(nrs[i])>maxnrs) {maxnrs = PetscAbsScalar(nrs[i]); s[k] = ds[i];}
425:     }
426:     if (ndptx>k) {
427:       for (i=1;i<ndptx;i++) {
428:         nrxi[i] *= ((dxi[i]-s[k-1])/(1.0-dxi[i]/xi[k-1]))/beta[k-1];
429:         if (PetscAbsScalar(nrxi[i])<minnrxi) {minnrxi = PetscAbsScalar(nrxi[i]); xi[k] = dxi[i];}
430:       }
431:       if (PetscAbsScalar(xi[k])<10*PETSC_MACHINE_EPSILON) SETERRQ2(PetscObjectComm((PetscObject)nep),1,"Singularity point %D is nearly zero: %g; consider removing the singularity or shifting the problem",k,(double)PetscAbsScalar(xi[k]));
432:     } else xi[k] = PETSC_INFINITY;
433:     beta[k] = maxnrs;
434:   }
435:   PetscFree5(ds,dsi,dxi,nrs,nrxi);
436:   return(0);
437: }

439: PetscErrorCode NEPNLEIGSEvalNRTFunct(NEP nep,PetscInt k,PetscScalar sigma,PetscScalar *b)
440: {
441:   NEP_NLEIGS  *ctx=(NEP_NLEIGS*)nep->data;
442:   PetscInt    i;
443:   PetscScalar *beta=ctx->beta,*s=ctx->s,*xi=ctx->xi;

446:   b[0] = 1.0/beta[0];
447:   for (i=0;i<k;i++) {
448:     b[i+1] = ((sigma-s[i])*b[i])/(beta[i+1]*(1.0-sigma/xi[i]));
449:   }
450:   return(0);
451: }

453: static PetscErrorCode MatMult_Fun(Mat A,Vec x,Vec y)
454: {
455:   PetscErrorCode      ierr;
456:   NEP_NLEIGS_MATSHELL *ctx;
457:   PetscInt            i;

460:   MatShellGetContext(A,(void**)&ctx);
461:   MatMult(ctx->A[0],x,y);
462:   if (ctx->coeff[0]!=1.0) { VecScale(y,ctx->coeff[0]); }
463:   for (i=1;i<ctx->nmat;i++) {
464:     MatMult(ctx->A[i],x,ctx->t);
465:     VecAXPY(y,ctx->coeff[i],ctx->t);
466:   }
467:   return(0);
468: }

470: static PetscErrorCode MatMultTranspose_Fun(Mat A,Vec x,Vec y)
471: {
472:   PetscErrorCode      ierr;
473:   NEP_NLEIGS_MATSHELL *ctx;
474:   PetscInt            i;

477:   MatShellGetContext(A,(void**)&ctx);
478:   MatMultTranspose(ctx->A[0],x,y);
479:   if (ctx->coeff[0]!=1.0) { VecScale(y,ctx->coeff[0]); }
480:   for (i=1;i<ctx->nmat;i++) {
481:     MatMultTranspose(ctx->A[i],x,ctx->t);
482:     VecAXPY(y,ctx->coeff[i],ctx->t);
483:   }
484:   return(0);
485: }

487: static PetscErrorCode MatGetDiagonal_Fun(Mat A,Vec diag)
488: {
489:   PetscErrorCode      ierr;
490:   NEP_NLEIGS_MATSHELL *ctx;
491:   PetscInt            i;

494:   MatShellGetContext(A,(void**)&ctx);
495:   MatGetDiagonal(ctx->A[0],diag);
496:   if (ctx->coeff[0]!=1.0) { VecScale(diag,ctx->coeff[0]); }
497:   for (i=1;i<ctx->nmat;i++) {
498:     MatGetDiagonal(ctx->A[i],ctx->t);
499:     VecAXPY(diag,ctx->coeff[i],ctx->t);
500:   }
501:   return(0);
502: }

504: static PetscErrorCode MatDuplicate_Fun(Mat A,MatDuplicateOption op,Mat *B)
505: {
506:   PetscInt            n,i;
507:   NEP_NLEIGS_MATSHELL *ctxnew,*ctx;
508:   void                (*fun)(void);
509:   PetscErrorCode      ierr;

512:   MatShellGetContext(A,(void**)&ctx);
513:   PetscNew(&ctxnew);
514:   ctxnew->nmat = ctx->nmat;
515:   ctxnew->maxnmat = ctx->maxnmat;
516:   PetscMalloc2(ctxnew->maxnmat,&ctxnew->A,ctxnew->maxnmat,&ctxnew->coeff);
517:   for (i=0;i<ctx->nmat;i++) {
518:     PetscObjectReference((PetscObject)ctx->A[i]);
519:     ctxnew->A[i] = ctx->A[i];
520:     ctxnew->coeff[i] = ctx->coeff[i];
521:   }
522:   MatGetSize(ctx->A[0],&n,NULL);
523:   VecDuplicate(ctx->t,&ctxnew->t);
524:   MatCreateShell(PetscObjectComm((PetscObject)A),n,n,n,n,(void*)ctxnew,B);
525:   MatShellSetManageScalingShifts(*B);
526:   MatShellGetOperation(A,MATOP_MULT,&fun);
527:   MatShellSetOperation(*B,MATOP_MULT,fun);
528:   MatShellGetOperation(A,MATOP_MULT_TRANSPOSE,&fun);
529:   MatShellSetOperation(*B,MATOP_MULT_TRANSPOSE,fun);
530:   MatShellGetOperation(A,MATOP_GET_DIAGONAL,&fun);
531:   MatShellSetOperation(*B,MATOP_GET_DIAGONAL,fun);
532:   MatShellGetOperation(A,MATOP_DUPLICATE,&fun);
533:   MatShellSetOperation(*B,MATOP_DUPLICATE,fun);
534:   MatShellGetOperation(A,MATOP_DESTROY,&fun);
535:   MatShellSetOperation(*B,MATOP_DESTROY,fun);
536:   MatShellGetOperation(A,MATOP_AXPY,&fun);
537:   MatShellSetOperation(*B,MATOP_AXPY,fun);
538:   return(0);
539: }

541: static PetscErrorCode MatDestroy_Fun(Mat A)
542: {
543:   NEP_NLEIGS_MATSHELL *ctx;
544:   PetscErrorCode      ierr;
545:   PetscInt            i;

548:   if (A) {
549:     MatShellGetContext(A,(void**)&ctx);
550:     for (i=0;i<ctx->nmat;i++) {
551:       MatDestroy(&ctx->A[i]);
552:     }
553:     VecDestroy(&ctx->t);
554:     PetscFree2(ctx->A,ctx->coeff);
555:     PetscFree(ctx);
556:   }
557:   return(0);
558: }

560: static PetscErrorCode MatAXPY_Fun(Mat Y,PetscScalar a,Mat X,MatStructure str)
561: {
562:   NEP_NLEIGS_MATSHELL *ctxY,*ctxX;
563:   PetscErrorCode      ierr;
564:   PetscInt            i,j;
565:   PetscBool           found;

568:   MatShellGetContext(Y,(void**)&ctxY);
569:   MatShellGetContext(X,(void**)&ctxX);
570:   for (i=0;i<ctxX->nmat;i++) {
571:     found = PETSC_FALSE;
572:     for (j=0;!found&&j<ctxY->nmat;j++) {
573:       if (ctxX->A[i]==ctxY->A[j]) {
574:         found = PETSC_TRUE;
575:         ctxY->coeff[j] += a*ctxX->coeff[i];
576:       }
577:     }
578:     if (!found) {
579:       ctxY->coeff[ctxY->nmat] = a*ctxX->coeff[i];
580:       ctxY->A[ctxY->nmat++] = ctxX->A[i];
581:       PetscObjectReference((PetscObject)ctxX->A[i]);
582:     }
583:   }
584:   return(0);
585: }

587: static PetscErrorCode MatScale_Fun(Mat M,PetscScalar a)
588: {
589:   NEP_NLEIGS_MATSHELL *ctx;
590:   PetscErrorCode      ierr;
591:   PetscInt            i;

594:   MatShellGetContext(M,(void**)&ctx);
595:   for (i=0;i<ctx->nmat;i++) ctx->coeff[i] *= a;
596:   return(0);
597: }

599: static PetscErrorCode NLEIGSMatToMatShellArray(Mat M,Mat *Ms,PetscInt maxnmat)
600: {
601:   PetscErrorCode      ierr;
602:   NEP_NLEIGS_MATSHELL *ctx;
603:   PetscInt            n;
604:   PetscBool           has;

607:   MatHasOperation(M,MATOP_DUPLICATE,&has);
608:   if (!has) SETERRQ(PetscObjectComm((PetscObject)M),1,"MatDuplicate operation required");
609:   PetscNew(&ctx);
610:   ctx->maxnmat = maxnmat;
611:   PetscMalloc2(ctx->maxnmat,&ctx->A,ctx->maxnmat,&ctx->coeff);
612:   MatDuplicate(M,MAT_COPY_VALUES,&ctx->A[0]);
613:   ctx->nmat = 1;
614:   ctx->coeff[0] = 1.0;
615:   MatCreateVecs(M,&ctx->t,NULL);
616:   MatGetSize(M,&n,NULL);
617:   MatCreateShell(PetscObjectComm((PetscObject)M),n,n,n,n,(void*)ctx,Ms);
618:   MatShellSetManageScalingShifts(*Ms);
619:   MatShellSetOperation(*Ms,MATOP_MULT,(void(*)(void))MatMult_Fun);
620:   MatShellSetOperation(*Ms,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_Fun);
621:   MatShellSetOperation(*Ms,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Fun);
622:   MatShellSetOperation(*Ms,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_Fun);
623:   MatShellSetOperation(*Ms,MATOP_DESTROY,(void(*)(void))MatDestroy_Fun);
624:   MatShellSetOperation(*Ms,MATOP_AXPY,(void(*)(void))MatAXPY_Fun);
625:   MatShellSetOperation(*Ms,MATOP_SCALE,(void(*)(void))MatScale_Fun);
626:   return(0);
627: }

629: static PetscErrorCode NEPNLEIGSDividedDifferences_split(NEP nep)
630: {
632:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
633:   PetscInt       k,j,i,maxnmat,nmax;
634:   PetscReal      norm0,norm,*matnorm;
635:   PetscScalar    *s=ctx->s,*beta=ctx->beta,*xi=ctx->xi,*b,alpha,*coeffs,*pK,*pH,sone=1.0;
636:   Mat            T,Ts,K,H;
637:   PetscBool      shell,hasmnorm=PETSC_FALSE,matrix=PETSC_TRUE;
638:   PetscBLASInt   n_;

641:   nmax = ctx->ddmaxit;
642:   PetscMalloc1(nep->nt*nmax,&ctx->coeffD);
643:   PetscMalloc3(nmax+1,&b,nmax+1,&coeffs,nep->nt,&matnorm);
644:   for (j=0;j<nep->nt;j++) {
645:     MatHasOperation(nep->A[j],MATOP_NORM,&hasmnorm);
646:     if (!hasmnorm) break;
647:     MatNorm(nep->A[j],NORM_INFINITY,matnorm+j);
648:   }
649:   /* Try matrix functions scheme */
650:   PetscCalloc2(nmax*nmax,&pK,nmax*nmax,&pH);
651:   for (i=0;i<nmax-1;i++) {
652:     pK[(nmax+1)*i]   = 1.0;
653:     pK[(nmax+1)*i+1] = beta[i+1]/xi[i];
654:     pH[(nmax+1)*i]   = s[i];
655:     pH[(nmax+1)*i+1] = beta[i+1];
656:   }
657:   pH[nmax*nmax-1] = s[nmax-1];
658:   pK[nmax*nmax-1] = 1.0;
659:   PetscBLASIntCast(nmax,&n_);
660:   PetscStackCallBLAS("BLAStrsm",BLAStrsm_("R","L","N","U",&n_,&n_,&sone,pK,&n_,pH,&n_));
661:   /* The matrix to be used is in H. K will be a work-space matrix */
662:   MatCreateSeqDense(PETSC_COMM_SELF,nmax,nmax,pH,&H);
663:   MatCreateSeqDense(PETSC_COMM_SELF,nmax,nmax,pK,&K);
664:   for (j=0;matrix&&j<nep->nt;j++) {
665:     PetscPushErrorHandler(PetscReturnErrorHandler,NULL);
666:     FNEvaluateFunctionMat(nep->f[j],H,K);
667:     PetscPopErrorHandler();
668:     if (!ierr) {
669:       for (i=0;i<nmax;i++) { ctx->coeffD[j+i*nep->nt] = pK[i]*beta[0]; }
670:     } else {
671:       matrix = PETSC_FALSE;
672:       PetscFPTrapPop();
673:     }
674:   }
675:   MatDestroy(&H);
676:   MatDestroy(&K);
677:   if (!matrix) {
678:     for (j=0;j<nep->nt;j++) {
679:       FNEvaluateFunction(nep->f[j],s[0],ctx->coeffD+j);
680:       ctx->coeffD[j] *= beta[0];
681:     }
682:   }
683:   if (hasmnorm) {
684:     norm0 = 0.0;
685:     for (j=0;j<nep->nt;j++) norm0 += matnorm[j]*PetscAbsScalar(ctx->coeffD[j]);
686:   } else {
687:     norm0 = 0.0;
688:     for (j=0;j<nep->nt;j++) norm0 = PetscMax(PetscAbsScalar(ctx->coeffD[j]),norm0);
689:   }
690:   ctx->nmat = ctx->ddmaxit;
691:   for (k=1;k<ctx->ddmaxit;k++) {
692:     if (!matrix) {
693:       NEPNLEIGSEvalNRTFunct(nep,k,s[k],b);
694:       for (i=0;i<nep->nt;i++) {
695:         FNEvaluateFunction(nep->f[i],s[k],ctx->coeffD+k*nep->nt+i);
696:         for (j=0;j<k;j++) {
697:           ctx->coeffD[k*nep->nt+i] -= b[j]*ctx->coeffD[i+nep->nt*j];
698:         }
699:         ctx->coeffD[k*nep->nt+i] /= b[k];
700:       }
701:     }
702:     if (hasmnorm) {
703:       norm = 0.0;
704:       for (j=0;j<nep->nt;j++) norm += matnorm[j]*PetscAbsScalar(ctx->coeffD[k*nep->nt+j]);
705:     } else {
706:       norm = 0.0;
707:       for (j=0;j<nep->nt;j++) norm = PetscMax(PetscAbsScalar(ctx->coeffD[k*nep->nt+j]),norm);
708:     }
709:     if (k>1 && norm/norm0 < ctx->ddtol) {
710:       ctx->nmat = k+1;
711:       break;
712:     }
713:   }
714:   if (!ctx->ksp) { NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp); }
715:   MatIsShell(nep->A[0],&shell);
716:   maxnmat = PetscMax(ctx->ddmaxit,nep->nt);
717:   for (i=0;i<ctx->nshiftsw;i++) {
718:     NEPNLEIGSEvalNRTFunct(nep,ctx->nmat-1,ctx->shifts[i],coeffs);
719:     if (!shell) {
720:       MatDuplicate(nep->A[0],MAT_COPY_VALUES,&T);
721:     } else {
722:       NLEIGSMatToMatShellArray(nep->A[0],&T,maxnmat);
723:     }
724:     alpha = 0.0;
725:     for (j=0;j<ctx->nmat;j++) alpha += coeffs[j]*ctx->coeffD[j*nep->nt];
726:     MatScale(T,alpha);
727:     for (k=1;k<nep->nt;k++) {
728:       alpha = 0.0;
729:       for (j=0;j<ctx->nmat;j++) alpha += coeffs[j]*ctx->coeffD[j*nep->nt+k];
730:       if (shell) {
731:         NLEIGSMatToMatShellArray(nep->A[k],&Ts,maxnmat);
732:       }
733:       MatAXPY(T,alpha,shell?Ts:nep->A[k],nep->mstr);
734:       if (shell) {
735:         MatDestroy(&Ts);
736:       }
737:     }
738:     KSPSetOperators(ctx->ksp[i],T,T);
739:     KSPSetUp(ctx->ksp[i]);
740:     MatDestroy(&T);
741:   }
742:   PetscFree3(b,coeffs,matnorm);
743:   PetscFree2(pK,pH);
744:   return(0);
745: }

747: static PetscErrorCode NEPNLEIGSDividedDifferences_callback(NEP nep)
748: {
750:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
751:   PetscInt       k,j,i,maxnmat;
752:   PetscReal      norm0,norm;
753:   PetscScalar    *s=ctx->s,*beta=ctx->beta,*b,*coeffs;
754:   Mat            *D=ctx->D,T;
755:   PetscBool      shell,has,vec=PETSC_FALSE;
756:   PetscRandom    rand;
757:   Vec            w[2];

760:   PetscMalloc2(ctx->ddmaxit+1,&b,ctx->ddmaxit+1,&coeffs);
761:   BVGetRandomContext(nep->V,&rand);
762:   T = nep->function;
763:   NEPComputeFunction(nep,s[0],T,T);
764:   MatIsShell(T,&shell);
765:   maxnmat = PetscMax(ctx->ddmaxit,nep->nt);
766:   if (!shell) {
767:     MatDuplicate(T,MAT_COPY_VALUES,&D[0]);
768:   } else {
769:     NLEIGSMatToMatShellArray(T,&D[0],maxnmat);
770:   }
771:   if (beta[0]!=1.0) {
772:     MatScale(D[0],1.0/beta[0]);
773:   }
774:   MatHasOperation(D[0],MATOP_NORM,&has);
775:   if (has) {
776:     MatNorm(D[0],NORM_FROBENIUS,&norm0);
777:   } else {
778:     MatCreateVecs(D[0],NULL,&w[0]);
779:     VecDuplicate(w[0],&w[1]);
780:     VecDuplicate(w[0],&ctx->vrn);
781:     VecSetRandomNormal(ctx->vrn,rand,w[0],w[1]);
782:     VecNormalize(ctx->vrn,NULL);
783:     vec = PETSC_TRUE;
784:     MatNormEstimate(D[0],ctx->vrn,w[0],&norm0);
785:   }
786:   ctx->nmat = ctx->ddmaxit;
787:   for (k=1;k<ctx->ddmaxit;k++) {
788:     NEPNLEIGSEvalNRTFunct(nep,k,s[k],b);
789:     NEPComputeFunction(nep,s[k],T,T);
790:     if (!shell) {
791:       MatDuplicate(T,MAT_COPY_VALUES,&D[k]);
792:     } else {
793:       NLEIGSMatToMatShellArray(T,&D[k],maxnmat);
794:     }
795:     for (j=0;j<k;j++) {
796:       MatAXPY(D[k],-b[j],D[j],nep->mstr);
797:     }
798:     MatScale(D[k],1.0/b[k]);
799:     MatHasOperation(D[k],MATOP_NORM,&has);
800:     if (has) {
801:       MatNorm(D[k],NORM_FROBENIUS,&norm);
802:     } else {
803:       if (!vec) {
804:         MatCreateVecs(D[k],NULL,&w[0]);
805:         VecDuplicate(w[0],&w[1]);
806:         VecDuplicate(w[0],&ctx->vrn);
807:         VecSetRandomNormal(ctx->vrn,rand,w[0],w[1]);
808:         VecNormalize(ctx->vrn,NULL);
809:         vec = PETSC_TRUE;
810:       }
811:       MatNormEstimate(D[k],ctx->vrn,w[0],&norm);
812:     }
813:     if (k>1 && norm/norm0 < ctx->ddtol && k>1) {
814:       ctx->nmat = k+1;
815:       break;
816:     }
817:   }
818:   if (!ctx->ksp) { NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp); }
819:   for (i=0;i<ctx->nshiftsw;i++) {
820:     NEPNLEIGSEvalNRTFunct(nep,ctx->nmat-1,ctx->shifts[i],coeffs);
821:     MatDuplicate(ctx->D[0],MAT_COPY_VALUES,&T);
822:     if (coeffs[0]!=1.0) { MatScale(T,coeffs[0]); }
823:     for (j=1;j<ctx->nmat;j++) {
824:       MatAXPY(T,coeffs[j],ctx->D[j],nep->mstr);
825:     }
826:     KSPSetOperators(ctx->ksp[i],T,T);
827:     KSPSetUp(ctx->ksp[i]);
828:     MatDestroy(&T);
829:   }
830:   PetscFree2(b,coeffs);
831:   if (vec) {
832:     VecDestroy(&w[0]);
833:     VecDestroy(&w[1]);
834:   }
835:   return(0);
836: }

838: /*
839:    NEPKrylovConvergence - This is the analogue to EPSKrylovConvergence.
840: */
841: PetscErrorCode NEPNLEIGSKrylovConvergence(NEP nep,PetscBool getall,PetscInt kini,PetscInt nits,PetscReal betah,PetscScalar betak,PetscInt *kout,Vec *w)
842: {
844:   PetscInt       k,newk,marker,inside;
845:   PetscScalar    re,im;
846:   PetscReal      resnorm,tt;
847:   PetscBool      istrivial;
848:   NEP_NLEIGS     *ctx = (NEP_NLEIGS*)nep->data;

851:   RGIsTrivial(nep->rg,&istrivial);
852:   marker = -1;
853:   if (nep->trackall) getall = PETSC_TRUE;
854:   for (k=kini;k<kini+nits;k++) {
855:     /* eigenvalue */
856:     re = nep->eigr[k];
857:     im = nep->eigi[k];
858:     if (!istrivial) {
859:       if (!ctx->nshifts) {
860:         NEPNLEIGSBackTransform((PetscObject)nep,1,&re,&im);
861:       }
862:       RGCheckInside(nep->rg,1,&re,&im,&inside);
863:       if (marker==-1 && inside<0) marker = k;
864:     }
865:     newk = k;
866:     DSVectors(nep->ds,DS_MAT_X,&newk,&resnorm);
867:     tt = (ctx->nshifts)?SlepcAbsEigenvalue(betak-nep->eigr[k]*betah,nep->eigi[k]*betah):betah;
868:     resnorm *=  PetscAbsReal(tt);
869:     /* error estimate */
870:     (*nep->converged)(nep,nep->eigr[k],nep->eigi[k],resnorm,&nep->errest[k],nep->convergedctx);
871:     if (marker==-1 && nep->errest[k] >= nep->tol) marker = k;
872:     if (newk==k+1) {
873:       nep->errest[k+1] = nep->errest[k];
874:       k++;
875:     }
876:     if (marker!=-1 && !getall) break;
877:   }
878:   if (marker!=-1) k = marker;
879:   *kout = k;
880:   return(0);
881: }

883: PetscErrorCode NEPSetUp_NLEIGS(NEP nep)
884: {
886:   PetscInt       k,in;
887:   PetscScalar    zero=0.0;
888:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
889:   SlepcSC        sc;
890:   PetscBool      istrivial;

893:   NEPSetDimensions_Default(nep,nep->nev,&nep->ncv,&nep->mpd);
894:   if (nep->ncv>nep->nev+nep->mpd) SETERRQ(PetscObjectComm((PetscObject)nep),1,"The value of ncv must not be larger than nev+mpd");
895:   if (nep->max_it==PETSC_DEFAULT) nep->max_it = PetscMax(5000,2*nep->n/nep->ncv);
896:   if (!ctx->ddmaxit) ctx->ddmaxit = LBPOINTS;
897:   RGIsTrivial(nep->rg,&istrivial);
898:   if (istrivial) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"NEPNLEIGS requires a nontrivial region defining the target set");
899:   if (!nep->which) nep->which = NEP_TARGET_MAGNITUDE;
900:   if (nep->which!=NEP_TARGET_MAGNITUDE && nep->which!=NEP_TARGET_REAL && nep->which!=NEP_TARGET_IMAGINARY && nep->which!=NEP_WHICH_USER) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"This solver supports only target selection of eigenvalues");

902:   /* Initialize the NLEIGS context structure */
903:   k = ctx->ddmaxit;
904:   PetscMalloc4(k,&ctx->s,k,&ctx->xi,k,&ctx->beta,k,&ctx->D);
905:   nep->data = ctx;
906:   if (nep->tol==PETSC_DEFAULT) nep->tol = SLEPC_DEFAULT_TOL;
907:   if (ctx->ddtol==PETSC_DEFAULT) ctx->ddtol = nep->tol/10.0;
908:   if (!ctx->keep) ctx->keep = 0.5;

910:   /* Compute Leja-Bagby points and scaling values */
911:   NEPNLEIGSLejaBagbyPoints(nep);
912:   if (nep->problem_type!=NEP_RATIONAL) {
913:     RGCheckInside(nep->rg,1,&nep->target,&zero,&in);
914:     if (in<0) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"The target is not inside the target set");
915:   }

917:   /* Compute the divided difference matrices */
918:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
919:     NEPNLEIGSDividedDifferences_split(nep);
920:   } else {
921:     NEPNLEIGSDividedDifferences_callback(nep);
922:   }
923:   NEPAllocateSolution(nep,ctx->nmat-1);
924:   NEPSetWorkVecs(nep,4);
925:   if (!ctx->fullbasis) {
926:     if (nep->twosided) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Two-sided variant requires the full-basis option, rerun with -nep_nleigs_full_basis");
927:     /* set-up DS and transfer split operator functions */
928:     DSSetType(nep->ds,ctx->nshifts?DSGNHEP:DSNHEP);
929:     DSAllocate(nep->ds,nep->ncv+1);
930:     DSGetSlepcSC(nep->ds,&sc);
931:     if (!ctx->nshifts) sc->map = NEPNLEIGSBackTransform;
932:     DSSetExtraRow(nep->ds,PETSC_TRUE);
933:     sc->mapobj        = (PetscObject)nep;
934:     sc->rg            = nep->rg;
935:     sc->comparison    = nep->sc->comparison;
936:     sc->comparisonctx = nep->sc->comparisonctx;
937:     BVDestroy(&ctx->V);
938:     BVCreateTensor(nep->V,ctx->nmat-1,&ctx->V);
939:     nep->ops->solve          = NEPSolve_NLEIGS;
940:     nep->ops->computevectors = NEPComputeVectors_Schur;
941:   } else {
942:     NEPSetUp_NLEIGS_FullBasis(nep);
943:     nep->ops->solve          = NEPSolve_NLEIGS_FullBasis;
944:     nep->ops->computevectors = NULL;
945:   }
946:   return(0);
947: }

949: /*
950:   Extend the TOAR basis by applying the the matrix operator
951:   over a vector which is decomposed on the TOAR way
952:   Input:
953:     - S,V: define the latest Arnoldi vector (nv vectors in V)
954:   Output:
955:     - t: new vector extending the TOAR basis
956:     - r: temporally coefficients to compute the TOAR coefficients
957:          for the new Arnoldi vector
958:   Workspace: t_ (two vectors)
959: */
960: static PetscErrorCode NEPTOARExtendBasis(NEP nep,PetscInt idxrktg,PetscScalar *S,PetscInt ls,PetscInt nv,BV W,BV V,Vec t,PetscScalar *r,PetscInt lr,Vec *t_)
961: {
963:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
964:   PetscInt       deg=ctx->nmat-1,k,j;
965:   Vec            v=t_[0],q=t_[1],w;
966:   PetscScalar    *beta=ctx->beta,*s=ctx->s,*xi=ctx->xi,*coeffs,sigma;

969:   if (!ctx->ksp) { NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp); }
970:   sigma = ctx->shifts[idxrktg];
971:   BVSetActiveColumns(nep->V,0,nv);
972:   PetscMalloc1(ctx->nmat,&coeffs);
973:   if (PetscAbsScalar(s[deg-2]-sigma)<100*PETSC_MACHINE_EPSILON) SETERRQ(PETSC_COMM_SELF,1,"Breakdown in NLEIGS");
974:   /* i-part stored in (i-1) position */
975:   for (j=0;j<nv;j++) {
976:     r[(deg-2)*lr+j] = (S[(deg-2)*ls+j]+(beta[deg-1]/xi[deg-2])*S[(deg-1)*ls+j])/(s[deg-2]-sigma);
977:   }
978:   BVSetActiveColumns(W,0,deg);
979:   BVGetColumn(W,deg-1,&w);
980:   BVMultVec(V,1.0/beta[deg],0,w,S+(deg-1)*ls);
981:   BVRestoreColumn(W,deg-1,&w);
982:   BVGetColumn(W,deg-2,&w);
983:   BVMultVec(V,1.0,0.0,w,r+(deg-2)*lr);
984:   BVRestoreColumn(W,deg-2,&w);
985:   for (k=deg-2;k>0;k--) {
986:     if (PetscAbsScalar(s[k-1]-sigma)<100*PETSC_MACHINE_EPSILON) SETERRQ(PETSC_COMM_SELF,1,"Breakdown in NLEIGS");
987:     for (j=0;j<nv;j++) r[(k-1)*lr+j] = (S[(k-1)*ls+j]+(beta[k]/xi[k-1])*S[k*ls+j]-beta[k]*(1.0-sigma/xi[k-1])*r[(k)*lr+j])/(s[k-1]-sigma);
988:     BVGetColumn(W,k-1,&w);
989:     BVMultVec(V,1.0,0.0,w,r+(k-1)*lr);
990:     BVRestoreColumn(W,k-1,&w);
991:   }
992:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
993:     for (j=0;j<ctx->nmat-2;j++) coeffs[j] = ctx->coeffD[nep->nt*j];
994:     coeffs[ctx->nmat-2] = ctx->coeffD[nep->nt*(ctx->nmat-1)];
995:     BVMultVec(W,1.0,0.0,v,coeffs);
996:     MatMult(nep->A[0],v,q);
997:     for (k=1;k<nep->nt;k++) {
998:       for (j=0;j<ctx->nmat-2;j++) coeffs[j] = ctx->coeffD[nep->nt*j+k];
999:       coeffs[ctx->nmat-2] = ctx->coeffD[nep->nt*(ctx->nmat-1)+k];
1000:       BVMultVec(W,1.0,0,v,coeffs);
1001:       MatMult(nep->A[k],v,t);
1002:       VecAXPY(q,1.0,t);
1003:     }
1004:     KSPSolve(ctx->ksp[idxrktg],q,t);
1005:     VecScale(t,-1.0);
1006:   } else {
1007:     for (k=0;k<deg-1;k++) {
1008:       BVGetColumn(W,k,&w);
1009:       MatMult(ctx->D[k],w,q);
1010:       BVRestoreColumn(W,k,&w);
1011:       BVInsertVec(W,k,q);
1012:     }
1013:     BVGetColumn(W,deg-1,&w);
1014:     MatMult(ctx->D[deg],w,q);
1015:     BVRestoreColumn(W,k,&w);
1016:     BVInsertVec(W,k,q);
1017:     for (j=0;j<ctx->nmat-1;j++) coeffs[j] = 1.0;
1018:     BVMultVec(W,1.0,0.0,q,coeffs);
1019:     KSPSolve(ctx->ksp[idxrktg],q,t);
1020:     VecScale(t,-1.0);
1021:   }
1022:   PetscFree(coeffs);
1023:   return(0);
1024: }

1026: /*
1027:   Compute TOAR coefficients of the blocks of the new Arnoldi vector computed
1028: */
1029: static PetscErrorCode NEPTOARCoefficients(NEP nep,PetscScalar sigma,PetscInt nv,PetscScalar *S,PetscInt ls,PetscScalar *r,PetscInt lr,PetscScalar *x,PetscScalar *work)
1030: {
1032:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
1033:   PetscInt       k,j,d=ctx->nmat-1;
1034:   PetscScalar    *t=work;

1037:   NEPNLEIGSEvalNRTFunct(nep,d-1,sigma,t);
1038:   for (k=0;k<d-1;k++) {
1039:     for (j=0;j<=nv;j++) r[k*lr+j] += t[k]*x[j];
1040:   }
1041:   for (j=0;j<=nv;j++) r[(d-1)*lr+j] = t[d-1]*x[j];
1042:   return(0);
1043: }

1045: /*
1046:   Compute continuation vector coefficients for the Rational-Krylov run.
1047:   dim(work) >= (end-ini)*(end-ini+1) + end+1 + 2*(end-ini+1), dim(t) = end.
1048: */
1049: static PetscErrorCode NEPNLEIGS_RKcontinuation(NEP nep,PetscInt ini,PetscInt end,PetscScalar *K,PetscScalar *H,PetscInt ld,PetscScalar sigma,PetscScalar *S,PetscInt lds,PetscScalar *cont,PetscScalar *t,PetscScalar *work)
1050: {
1052:   PetscScalar    *x,*W,*tau,sone=1.0,szero=0.0;
1053:   PetscInt       i,j,n1,n,nwu=0;
1054:   PetscBLASInt   info,n_,n1_,one=1,dim,lds_;
1055:   NEP_NLEIGS     *ctx = (NEP_NLEIGS*)nep->data;

1058:   if (!ctx->nshifts || !end) {
1059:     t[0] = 1;
1060:     PetscArraycpy(cont,S+end*lds,lds);
1061:   } else {
1062:     n   = end-ini;
1063:     n1  = n+1;
1064:     x   = work+nwu;
1065:     nwu += end+1;
1066:     tau = work+nwu;
1067:     nwu += n;
1068:     W   = work+nwu;
1069:     nwu += n1*n;
1070:     for (j=ini;j<end;j++) {
1071:       for (i=ini;i<=end;i++) W[(j-ini)*n1+i-ini] = K[j*ld+i] -H[j*ld+i]*sigma;
1072:     }
1073:     PetscBLASIntCast(n,&n_);
1074:     PetscBLASIntCast(n1,&n1_);
1075:     PetscBLASIntCast(end+1,&dim);
1076:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1077:     PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&n1_,&n_,W,&n1_,tau,work+nwu,&n1_,&info));
1078:     SlepcCheckLapackInfo("geqrf",info);
1079:     for (i=0;i<end;i++) t[i] = 0.0;
1080:     t[end] = 1.0;
1081:     for (j=n-1;j>=0;j--) {
1082:       for (i=0;i<ini+j;i++) x[i] = 0.0;
1083:       x[ini+j] = 1.0;
1084:       for (i=j+1;i<n1;i++) x[i+ini] = W[i+n1*j];
1085:       tau[j] = PetscConj(tau[j]);
1086:       PetscStackCallBLAS("LAPACKlarf",LAPACKlarf_("L",&dim,&one,x,&one,tau+j,t,&dim,work+nwu));
1087:     }
1088:     PetscBLASIntCast(lds,&lds_);
1089:     PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&lds_,&n1_,&sone,S,&lds_,t,&one,&szero,cont,&one));
1090:     PetscFPTrapPop();
1091:   }
1092:   return(0);
1093: }

1095: /*
1096:   Compute a run of Arnoldi iterations
1097: */
1098: PetscErrorCode NEPNLEIGSTOARrun(NEP nep,PetscScalar *K,PetscScalar *H,PetscInt ldh,BV W,PetscInt k,PetscInt *M,PetscBool *breakdown,Vec *t_)
1099: {
1101:   NEP_NLEIGS     *ctx = (NEP_NLEIGS*)nep->data;
1102:   PetscInt       i,j,m=*M,lwa,deg=ctx->nmat-1,lds,nqt,ld,l;
1103:   Vec            t;
1104:   PetscReal      norm;
1105:   PetscScalar    *x,*work,*tt,sigma,*cont,*S;
1106:   PetscBool      lindep;
1107:   Mat            MS;

1110:   BVTensorGetFactors(ctx->V,NULL,&MS);
1111:   MatDenseGetArray(MS,&S);
1112:   BVGetSizes(nep->V,NULL,NULL,&ld);
1113:   lds = ld*deg;
1114:   BVGetActiveColumns(nep->V,&l,&nqt);
1115:   lwa = PetscMax(ld,deg)+(m+1)*(m+1)+4*(m+1);
1116:   PetscMalloc4(ld,&x,lwa,&work,m+1,&tt,lds,&cont);
1117:   BVSetActiveColumns(ctx->V,0,m);
1118:   for (j=k;j<m;j++) {
1119:     sigma = ctx->shifts[(++(ctx->idxrk))%ctx->nshiftsw];

1121:     /* Continuation vector */
1122:     NEPNLEIGS_RKcontinuation(nep,0,j,K,H,ldh,sigma,S,lds,cont,tt,work);

1124:     /* apply operator */
1125:     BVGetColumn(nep->V,nqt,&t);
1126:     NEPTOARExtendBasis(nep,(ctx->idxrk)%ctx->nshiftsw,cont,ld,nqt,W,nep->V,t,S+(j+1)*lds,ld,t_);
1127:     BVRestoreColumn(nep->V,nqt,&t);

1129:     /* orthogonalize */
1130:     BVOrthogonalizeColumn(nep->V,nqt,x,&norm,&lindep);
1131:     if (!lindep) {
1132:       x[nqt] = norm;
1133:       BVScaleColumn(nep->V,nqt,1.0/norm);
1134:       nqt++;
1135:     } else x[nqt] = 0.0;

1137:     NEPTOARCoefficients(nep,sigma,nqt-1,cont,ld,S+(j+1)*lds,ld,x,work);

1139:     /* Level-2 orthogonalization */
1140:     BVOrthogonalizeColumn(ctx->V,j+1,H+j*ldh,&norm,breakdown);
1141:     H[j+1+ldh*j] = norm;
1142:     if (ctx->nshifts) {
1143:       for (i=0;i<=j;i++) K[i+ldh*j] = sigma*H[i+ldh*j] + tt[i];
1144:       K[j+1+ldh*j] = sigma*H[j+1+ldh*j];
1145:     }
1146:     if (*breakdown) {
1147:       *M = j+1;
1148:       break;
1149:     }
1150:     BVScaleColumn(ctx->V,j+1,1.0/norm);
1151:     BVSetActiveColumns(nep->V,l,nqt);
1152:   }
1153:   PetscFree4(x,work,tt,cont);
1154:   MatDenseRestoreArray(MS,&S);
1155:   BVTensorRestoreFactors(ctx->V,NULL,&MS);
1156:   return(0);
1157: }

1159: PetscErrorCode NEPSolve_NLEIGS(NEP nep)
1160: {
1161:   PetscErrorCode    ierr;
1162:   NEP_NLEIGS        *ctx = (NEP_NLEIGS*)nep->data;
1163:   PetscInt          i,k=0,l,nv=0,ld,lds,ldds,nq;
1164:   PetscInt          deg=ctx->nmat-1,nconv=0,dsn,dsk;
1165:   PetscScalar       *H,*pU,*K,betak=0,*eigr,*eigi;
1166:   const PetscScalar *S;
1167:   PetscReal         betah;
1168:   PetscBool         falselock=PETSC_FALSE,breakdown=PETSC_FALSE;
1169:   BV                W;
1170:   Mat               MS,MQ,U;

1173:   if (ctx->lock) {
1174:     /* undocumented option to use a cheaper locking instead of the true locking */
1175:     PetscOptionsGetBool(NULL,NULL,"-nep_nleigs_falselocking",&falselock,NULL);
1176:   }

1178:   BVGetSizes(nep->V,NULL,NULL,&ld);
1179:   lds = deg*ld;
1180:   DSGetLeadingDimension(nep->ds,&ldds);
1181:   if (!ctx->nshifts) {
1182:     PetscMalloc2(nep->ncv,&eigr,nep->ncv,&eigi);
1183:   } else { eigr = nep->eigr; eigi = nep->eigi; }
1184:   BVDuplicateResize(nep->V,PetscMax(nep->nt-1,ctx->nmat-1),&W);


1187:   /* clean projected matrix (including the extra-arrow) */
1188:   DSGetArray(nep->ds,DS_MAT_A,&H);
1189:   PetscArrayzero(H,ldds*ldds);
1190:   DSRestoreArray(nep->ds,DS_MAT_A,&H);
1191:   if (ctx->nshifts) {
1192:     DSGetArray(nep->ds,DS_MAT_B,&H);
1193:     PetscArrayzero(H,ldds*ldds);
1194:     DSRestoreArray(nep->ds,DS_MAT_B,&H);
1195:   }

1197:   /* Get the starting Arnoldi vector */
1198:   BVTensorBuildFirstColumn(ctx->V,nep->nini);

1200:   /* Restart loop */
1201:   l = 0;
1202:   while (nep->reason == NEP_CONVERGED_ITERATING) {
1203:     nep->its++;

1205:     /* Compute an nv-step Krylov relation */
1206:     nv = PetscMin(nep->nconv+nep->mpd,nep->ncv);
1207:     if (ctx->nshifts) { DSGetArray(nep->ds,DS_MAT_A,&K); }
1208:     DSGetArray(nep->ds,ctx->nshifts?DS_MAT_B:DS_MAT_A,&H);
1209:     NEPNLEIGSTOARrun(nep,K,H,ldds,W,nep->nconv+l,&nv,&breakdown,nep->work);
1210:     betah = PetscAbsScalar(H[(nv-1)*ldds+nv]);
1211:     DSRestoreArray(nep->ds,ctx->nshifts?DS_MAT_B:DS_MAT_A,&H);
1212:     if (ctx->nshifts) {
1213:       betak = K[(nv-1)*ldds+nv];
1214:       DSRestoreArray(nep->ds,DS_MAT_A,&K);
1215:     }
1216:     DSSetDimensions(nep->ds,nv,0,nep->nconv,nep->nconv+l);
1217:     if (l==0) {
1218:       DSSetState(nep->ds,DS_STATE_INTERMEDIATE);
1219:     } else {
1220:       DSSetState(nep->ds,DS_STATE_RAW);
1221:     }

1223:     /* Solve projected problem */
1224:     DSSolve(nep->ds,nep->eigr,nep->eigi);
1225:     DSSort(nep->ds,nep->eigr,nep->eigi,NULL,NULL,NULL);
1226:     DSUpdateExtraRow(nep->ds);
1227:     DSSynchronize(nep->ds,nep->eigr,nep->eigi);

1229:     /* Check convergence */
1230:     NEPNLEIGSKrylovConvergence(nep,PETSC_FALSE,nep->nconv,nv-nep->nconv,betah,betak,&k,nep->work);
1231:     (*nep->stopping)(nep,nep->its,nep->max_it,k,nep->nev,&nep->reason,nep->stoppingctx);

1233:     /* Update l */
1234:     if (nep->reason != NEP_CONVERGED_ITERATING || breakdown) l = 0;
1235:     else {
1236:       l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
1237:       DSGetTruncateSize(nep->ds,k,nv,&l);
1238:       if (!breakdown) {
1239:         /* Prepare the Rayleigh quotient for restart */
1240:         DSGetDimensions(nep->ds,&dsn,NULL,NULL,&dsk,NULL);
1241:         DSSetDimensions(nep->ds,dsn,0,k,dsk);
1242:         DSTruncate(nep->ds,k+l,PETSC_FALSE);
1243:       }
1244:     }
1245:     nconv = k;
1246:     if (!ctx->lock && nep->reason == NEP_CONVERGED_ITERATING && !breakdown) { l += k; k = 0; }
1247:     if (l) { PetscInfo1(nep,"Preparing to restart keeping l=%D vectors\n",l); }

1249:     /* Update S */
1250:     DSGetMat(nep->ds,ctx->nshifts?DS_MAT_Z:DS_MAT_Q,&MQ);
1251:     BVMultInPlace(ctx->V,MQ,nep->nconv,k+l);
1252:     MatDestroy(&MQ);

1254:     /* Copy last column of S */
1255:     BVCopyColumn(ctx->V,nv,k+l);

1257:     if (breakdown && nep->reason == NEP_CONVERGED_ITERATING) {
1258:       /* Stop if breakdown */
1259:       PetscInfo2(nep,"Breakdown (it=%D norm=%g)\n",nep->its,(double)betah);
1260:       nep->reason = NEP_DIVERGED_BREAKDOWN;
1261:     }
1262:     if (nep->reason != NEP_CONVERGED_ITERATING) l--;
1263:     /* truncate S */
1264:     BVGetActiveColumns(nep->V,NULL,&nq);
1265:     if (k+l+deg<=nq) {
1266:       BVSetActiveColumns(ctx->V,nep->nconv,k+l+1);
1267:       if (!falselock && ctx->lock) {
1268:         BVTensorCompress(ctx->V,k-nep->nconv);
1269:       } else {
1270:         BVTensorCompress(ctx->V,0);
1271:       }
1272:     }
1273:     nep->nconv = k;
1274:     if (!ctx->nshifts) {
1275:       for (i=0;i<nv;i++) { eigr[i] = nep->eigr[i]; eigi[i] = nep->eigi[i]; }
1276:       NEPNLEIGSBackTransform((PetscObject)nep,nv,eigr,eigi);
1277:     }
1278:     NEPMonitor(nep,nep->its,nconv,eigr,eigi,nep->errest,nv);
1279:   }
1280:   nep->nconv = nconv;
1281:   if (nep->nconv>0) {
1282:     BVSetActiveColumns(ctx->V,0,nep->nconv);
1283:     BVGetActiveColumns(nep->V,NULL,&nq);
1284:     BVSetActiveColumns(nep->V,0,nq);
1285:     if (nq>nep->nconv) {
1286:       BVTensorCompress(ctx->V,nep->nconv);
1287:       BVSetActiveColumns(nep->V,0,nep->nconv);
1288:       nq = nep->nconv;
1289:     }
1290:     if (ctx->nshifts) {
1291:       DSGetMat(nep->ds,DS_MAT_B,&MQ);
1292:       BVMultInPlace(ctx->V,MQ,0,nep->nconv);
1293:       MatDestroy(&MQ);
1294:     }
1295:     BVTensorGetFactors(ctx->V,NULL,&MS);
1296:     MatDenseGetArrayRead(MS,&S);
1297:     PetscMalloc1(nq*nep->nconv,&pU);
1298:     for (i=0;i<nep->nconv;i++) {
1299:       PetscArraycpy(pU+i*nq,S+i*lds,nq);
1300:     }
1301:     MatDenseRestoreArrayRead(MS,&S);
1302:     BVTensorRestoreFactors(ctx->V,NULL,&MS);
1303:     MatCreateSeqDense(PETSC_COMM_SELF,nq,nep->nconv,pU,&U);
1304:     BVSetActiveColumns(nep->V,0,nq);
1305:     BVMultInPlace(nep->V,U,0,nep->nconv);
1306:     BVSetActiveColumns(nep->V,0,nep->nconv);
1307:     MatDestroy(&U);
1308:     PetscFree(pU);
1309:     DSTruncate(nep->ds,nep->nconv,PETSC_TRUE);
1310:   }

1312:   /* Map eigenvalues back to the original problem */
1313:   if (!ctx->nshifts) {
1314:     NEPNLEIGSBackTransform((PetscObject)nep,nep->nconv,nep->eigr,nep->eigi);
1315:     PetscFree2(eigr,eigi);
1316:   }
1317:   BVDestroy(&W);
1318:   return(0);
1319: }

1321: static PetscErrorCode NEPNLEIGSSetSingularitiesFunction_NLEIGS(NEP nep,PetscErrorCode (*fun)(NEP,PetscInt*,PetscScalar*,void*),void *ctx)
1322: {
1323:   NEP_NLEIGS *nepctx=(NEP_NLEIGS*)nep->data;

1326:   if (fun) nepctx->computesingularities = fun;
1327:   if (ctx) nepctx->singularitiesctx     = ctx;
1328:   nep->state = NEP_STATE_INITIAL;
1329:   return(0);
1330: }

1332: /*@C
1333:    NEPNLEIGSSetSingularitiesFunction - Sets a user function to compute a discretization
1334:    of the singularity set (where T(.) is not analytic).

1336:    Logically Collective on nep

1338:    Input Parameters:
1339: +  nep - the NEP context
1340: .  fun - user function (if NULL then NEP retains any previously set value)
1341: -  ctx - [optional] user-defined context for private data for the function
1342:          (may be NULL, in which case NEP retains any previously set value)

1344:    Calling Sequence of fun:
1345: $   fun(NEP nep,PetscInt *maxnp,PetscScalar *xi,void *ctx)

1347: +   nep   - the NEP context
1348: .   maxnp - on input number of requested points in the discretization (can be set)
1349: .   xi    - computed values of the discretization
1350: -   ctx   - optional context, as set by NEPNLEIGSSetSingularitiesFunction()

1352:    Notes:
1353:    The user-defined function can set a smaller value of maxnp if necessary.
1354:    It is wrong to return a larger value.

1356:    If the problem type has been set to rational with NEPSetProblemType(),
1357:    then it is not necessary to set the singularities explicitly since the
1358:    solver will try to determine them automatically.

1360:    Level: intermediate

1362: .seealso: NEPNLEIGSGetSingularitiesFunction(), NEPSetProblemType()
1363: @*/
1364: PetscErrorCode NEPNLEIGSSetSingularitiesFunction(NEP nep,PetscErrorCode (*fun)(NEP,PetscInt*,PetscScalar*,void*),void *ctx)
1365: {

1370:   PetscTryMethod(nep,"NEPNLEIGSSetSingularitiesFunction_C",(NEP,PetscErrorCode(*)(NEP,PetscInt*,PetscScalar*,void*),void*),(nep,fun,ctx));
1371:   return(0);
1372: }

1374: static PetscErrorCode NEPNLEIGSGetSingularitiesFunction_NLEIGS(NEP nep,PetscErrorCode (**fun)(NEP,PetscInt*,PetscScalar*,void*),void **ctx)
1375: {
1376:   NEP_NLEIGS *nepctx=(NEP_NLEIGS*)nep->data;

1379:   if (fun) *fun = nepctx->computesingularities;
1380:   if (ctx) *ctx = nepctx->singularitiesctx;
1381:   return(0);
1382: }

1384: /*@C
1385:    NEPNLEIGSGetSingularitiesFunction - Returns the Function and optionally the user
1386:    provided context for computing a discretization of the singularity set.

1388:    Not Collective

1390:    Input Parameter:
1391: .  nep - the nonlinear eigensolver context

1393:    Output Parameters:
1394: +  fun - location to put the function (or NULL)
1395: -  ctx - location to stash the function context (or NULL)

1397:    Level: advanced

1399: .seealso: NEPNLEIGSSetSingularitiesFunction()
1400: @*/
1401: PetscErrorCode NEPNLEIGSGetSingularitiesFunction(NEP nep,PetscErrorCode (**fun)(NEP,PetscInt*,PetscScalar*,void*),void **ctx)
1402: {

1407:   PetscUseMethod(nep,"NEPNLEIGSGetSingularitiesFunction_C",(NEP,PetscErrorCode(**)(NEP,PetscInt*,PetscScalar*,void*),void**),(nep,fun,ctx));
1408:   return(0);
1409: }

1411: static PetscErrorCode NEPNLEIGSSetRestart_NLEIGS(NEP nep,PetscReal keep)
1412: {
1413:   NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;

1416:   if (keep==PETSC_DEFAULT) ctx->keep = 0.5;
1417:   else {
1418:     if (keep<0.1 || keep>0.9) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The keep argument must be in the range [0.1,0.9]");
1419:     ctx->keep = keep;
1420:   }
1421:   return(0);
1422: }

1424: /*@
1425:    NEPNLEIGSSetRestart - Sets the restart parameter for the NLEIGS
1426:    method, in particular the proportion of basis vectors that must be kept
1427:    after restart.

1429:    Logically Collective on nep

1431:    Input Parameters:
1432: +  nep  - the nonlinear eigensolver context
1433: -  keep - the number of vectors to be kept at restart

1435:    Options Database Key:
1436: .  -nep_nleigs_restart - Sets the restart parameter

1438:    Notes:
1439:    Allowed values are in the range [0.1,0.9]. The default is 0.5.

1441:    Level: advanced

1443: .seealso: NEPNLEIGSGetRestart()
1444: @*/
1445: PetscErrorCode NEPNLEIGSSetRestart(NEP nep,PetscReal keep)
1446: {

1452:   PetscTryMethod(nep,"NEPNLEIGSSetRestart_C",(NEP,PetscReal),(nep,keep));
1453:   return(0);
1454: }

1456: static PetscErrorCode NEPNLEIGSGetRestart_NLEIGS(NEP nep,PetscReal *keep)
1457: {
1458:   NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;

1461:   *keep = ctx->keep;
1462:   return(0);
1463: }

1465: /*@
1466:    NEPNLEIGSGetRestart - Gets the restart parameter used in the NLEIGS method.

1468:    Not Collective

1470:    Input Parameter:
1471: .  nep - the nonlinear eigensolver context

1473:    Output Parameter:
1474: .  keep - the restart parameter

1476:    Level: advanced

1478: .seealso: NEPNLEIGSSetRestart()
1479: @*/
1480: PetscErrorCode NEPNLEIGSGetRestart(NEP nep,PetscReal *keep)
1481: {

1487:   PetscUseMethod(nep,"NEPNLEIGSGetRestart_C",(NEP,PetscReal*),(nep,keep));
1488:   return(0);
1489: }

1491: static PetscErrorCode NEPNLEIGSSetLocking_NLEIGS(NEP nep,PetscBool lock)
1492: {
1493:   NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;

1496:   ctx->lock = lock;
1497:   return(0);
1498: }

1500: /*@
1501:    NEPNLEIGSSetLocking - Choose between locking and non-locking variants of
1502:    the NLEIGS method.

1504:    Logically Collective on nep

1506:    Input Parameters:
1507: +  nep  - the nonlinear eigensolver context
1508: -  lock - true if the locking variant must be selected

1510:    Options Database Key:
1511: .  -nep_nleigs_locking - Sets the locking flag

1513:    Notes:
1514:    The default is to lock converged eigenpairs when the method restarts.
1515:    This behaviour can be changed so that all directions are kept in the
1516:    working subspace even if already converged to working accuracy (the
1517:    non-locking variant).

1519:    Level: advanced

1521: .seealso: NEPNLEIGSGetLocking()
1522: @*/
1523: PetscErrorCode NEPNLEIGSSetLocking(NEP nep,PetscBool lock)
1524: {

1530:   PetscTryMethod(nep,"NEPNLEIGSSetLocking_C",(NEP,PetscBool),(nep,lock));
1531:   return(0);
1532: }

1534: static PetscErrorCode NEPNLEIGSGetLocking_NLEIGS(NEP nep,PetscBool *lock)
1535: {
1536:   NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;

1539:   *lock = ctx->lock;
1540:   return(0);
1541: }

1543: /*@
1544:    NEPNLEIGSGetLocking - Gets the locking flag used in the NLEIGS method.

1546:    Not Collective

1548:    Input Parameter:
1549: .  nep - the nonlinear eigensolver context

1551:    Output Parameter:
1552: .  lock - the locking flag

1554:    Level: advanced

1556: .seealso: NEPNLEIGSSetLocking()
1557: @*/
1558: PetscErrorCode NEPNLEIGSGetLocking(NEP nep,PetscBool *lock)
1559: {

1565:   PetscUseMethod(nep,"NEPNLEIGSGetLocking_C",(NEP,PetscBool*),(nep,lock));
1566:   return(0);
1567: }

1569: static PetscErrorCode NEPNLEIGSSetInterpolation_NLEIGS(NEP nep,PetscReal tol,PetscInt degree)
1570: {
1572:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;

1575:   if (tol == PETSC_DEFAULT) {
1576:     ctx->ddtol = PETSC_DEFAULT;
1577:     nep->state = NEP_STATE_INITIAL;
1578:   } else {
1579:     if (tol <= 0.0) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of tol. Must be > 0");
1580:     ctx->ddtol = tol;
1581:   }
1582:   if (degree == PETSC_DEFAULT || degree == PETSC_DECIDE) {
1583:     ctx->ddmaxit = 0;
1584:     if (nep->state) { NEPReset(nep); }
1585:     nep->state = NEP_STATE_INITIAL;
1586:   } else {
1587:     if (degree <= 0) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of degree. Must be > 0");
1588:     if (ctx->ddmaxit != degree) {
1589:       ctx->ddmaxit = degree;
1590:       if (nep->state) { NEPReset(nep); }
1591:       nep->state = NEP_STATE_INITIAL;
1592:     }
1593:   }
1594:   return(0);
1595: }

1597: /*@
1598:    NEPNLEIGSSetInterpolation - Sets the tolerance and maximum degree
1599:    when building the interpolation via divided differences.

1601:    Logically Collective on nep

1603:    Input Parameters:
1604: +  nep    - the nonlinear eigensolver context
1605: .  tol    - tolerance to stop computing divided differences
1606: -  degree - maximum degree of interpolation

1608:    Options Database Key:
1609: +  -nep_nleigs_interpolation_tol <tol> - Sets the tolerance to stop computing divided differences
1610: -  -nep_nleigs_interpolation_degree <degree> - Sets the maximum degree of interpolation

1612:    Notes:
1613:    Use PETSC_DEFAULT for either argument to assign a reasonably good value.

1615:    Level: advanced

1617: .seealso: NEPNLEIGSGetInterpolation()
1618: @*/
1619: PetscErrorCode NEPNLEIGSSetInterpolation(NEP nep,PetscReal tol,PetscInt degree)
1620: {

1627:   PetscTryMethod(nep,"NEPNLEIGSSetInterpolation_C",(NEP,PetscReal,PetscInt),(nep,tol,degree));
1628:   return(0);
1629: }

1631: static PetscErrorCode NEPNLEIGSGetInterpolation_NLEIGS(NEP nep,PetscReal *tol,PetscInt *degree)
1632: {
1633:   NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;

1636:   if (tol)    *tol    = ctx->ddtol;
1637:   if (degree) *degree = ctx->ddmaxit;
1638:   return(0);
1639: }

1641: /*@
1642:    NEPNLEIGSGetInterpolation - Gets the tolerance and maximum degree
1643:    when building the interpolation via divided differences.

1645:    Not Collective

1647:    Input Parameter:
1648: .  nep - the nonlinear eigensolver context

1650:    Output Parameter:
1651: +  tol    - tolerance to stop computing divided differences
1652: -  degree - maximum degree of interpolation

1654:    Level: advanced

1656: .seealso: NEPNLEIGSSetInterpolation()
1657: @*/
1658: PetscErrorCode NEPNLEIGSGetInterpolation(NEP nep,PetscReal *tol,PetscInt *degree)
1659: {

1664:   PetscTryMethod(nep,"NEPNLEIGSGetInterpolation_C",(NEP,PetscReal*,PetscInt*),(nep,tol,degree));
1665:   return(0);
1666: }

1668: static PetscErrorCode NEPNLEIGSSetRKShifts_NLEIGS(NEP nep,PetscInt ns,PetscScalar *shifts)
1669: {
1671:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
1672:   PetscInt       i;

1675:   if (ns<0) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"Number of shifts must be non-negative");
1676:   if (ctx->nshifts) { PetscFree(ctx->shifts); }
1677:   for (i=0;i<ctx->nshiftsw;i++) { KSPDestroy(&ctx->ksp[i]); }
1678:   PetscFree(ctx->ksp);
1679:   ctx->ksp = NULL;
1680:   if (ns) {
1681:     PetscMalloc1(ns,&ctx->shifts);
1682:     for (i=0;i<ns;i++) ctx->shifts[i] = shifts[i];
1683:   }
1684:   ctx->nshifts = ns;
1685:   nep->state   = NEP_STATE_INITIAL;
1686:   return(0);
1687: }

1689: /*@
1690:    NEPNLEIGSSetRKShifts - Sets a list of shifts to be used in the Rational
1691:    Krylov method.

1693:    Logically Collective on nep

1695:    Input Parameters:
1696: +  nep    - the nonlinear eigensolver context
1697: .  ns     - number of shifts
1698: -  shifts - array of scalar values specifying the shifts

1700:    Options Database Key:
1701: .  -nep_nleigs_rk_shifts - Sets the list of shifts

1703:    Notes:
1704:    If only one shift is provided, the built subspace built is equivalent to
1705:    shift-and-invert Krylov-Schur (provided that the absolute convergence
1706:    criterion is used).

1708:    In the case of real scalars, complex shifts are not allowed. In the
1709:    command line, a comma-separated list of complex values can be provided with
1710:    the format [+/-][realnumber][+/-]realnumberi with no spaces, e.g.
1711:    -nep_nleigs_rk_shifts 1.0+2.0i,1.5+2.0i,1.0+1.5i

1713:    Use ns=0 to remove previously set shifts.

1715:    Level: advanced

1717: .seealso: NEPNLEIGSGetRKShifts()
1718: @*/
1719: PetscErrorCode NEPNLEIGSSetRKShifts(NEP nep,PetscInt ns,PetscScalar shifts[])
1720: {

1727:   PetscTryMethod(nep,"NEPNLEIGSSetRKShifts_C",(NEP,PetscInt,PetscScalar*),(nep,ns,shifts));
1728:   return(0);
1729: }

1731: static PetscErrorCode NEPNLEIGSGetRKShifts_NLEIGS(NEP nep,PetscInt *ns,PetscScalar **shifts)
1732: {
1734:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
1735:   PetscInt       i;

1738:   *ns = ctx->nshifts;
1739:   if (ctx->nshifts) {
1740:     PetscMalloc1(ctx->nshifts,shifts);
1741:     for (i=0;i<ctx->nshifts;i++) (*shifts)[i] = ctx->shifts[i];
1742:   }
1743:   return(0);
1744: }

1746: /*@C
1747:    NEPNLEIGSGetRKShifts - Gets the list of shifts used in the Rational
1748:    Krylov method.

1750:    Not Collective

1752:    Input Parameter:
1753: .  nep - the nonlinear eigensolver context

1755:    Output Parameter:
1756: +  ns     - number of shifts
1757: -  shifts - array of shifts

1759:    Note:
1760:    The user is responsible for deallocating the returned array.

1762:    Level: advanced

1764: .seealso: NEPNLEIGSSetRKShifts()
1765: @*/
1766: PetscErrorCode NEPNLEIGSGetRKShifts(NEP nep,PetscInt *ns,PetscScalar *shifts[])
1767: {

1774:   PetscTryMethod(nep,"NEPNLEIGSGetRKShifts_C",(NEP,PetscInt*,PetscScalar**),(nep,ns,shifts));
1775:   return(0);
1776: }

1778: static PetscErrorCode NEPNLEIGSGetKSPs_NLEIGS(NEP nep,PetscInt *nsolve,KSP **ksp)
1779: {
1781:   NEP_NLEIGS     *ctx = (NEP_NLEIGS*)nep->data;
1782:   PetscInt       i;
1783:   PC             pc;

1786:   if (!ctx->ksp) {
1787:     NEPNLEIGSSetShifts(nep,&ctx->nshiftsw);
1788:     PetscMalloc1(ctx->nshiftsw,&ctx->ksp);
1789:     for (i=0;i<ctx->nshiftsw;i++) {
1790:       KSPCreate(PetscObjectComm((PetscObject)nep),&ctx->ksp[i]);
1791:       PetscObjectIncrementTabLevel((PetscObject)ctx->ksp[i],(PetscObject)nep,1);
1792:       KSPSetOptionsPrefix(ctx->ksp[i],((PetscObject)nep)->prefix);
1793:       KSPAppendOptionsPrefix(ctx->ksp[i],"nep_nleigs_");
1794:       PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->ksp[i]);
1795:       PetscObjectSetOptions((PetscObject)ctx->ksp[i],((PetscObject)nep)->options);
1796:       KSPSetErrorIfNotConverged(ctx->ksp[i],PETSC_TRUE);
1797:       KSPSetTolerances(ctx->ksp[i],SLEPC_DEFAULT_TOL,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
1798:       KSPGetPC(ctx->ksp[i],&pc);
1799:       KSPSetType(ctx->ksp[i],KSPPREONLY);
1800:       PCSetType(pc,PCLU);
1801:     }
1802:   }
1803:   if (nsolve) *nsolve = ctx->nshiftsw;
1804:   if (ksp)    *ksp    = ctx->ksp;
1805:   return(0);
1806: }

1808: /*@C
1809:    NEPNLEIGSGetKSPs - Retrieve the array of linear solver objects associated with
1810:    the nonlinear eigenvalue solver.

1812:    Not Collective

1814:    Input Parameter:
1815: .  nep - nonlinear eigenvalue solver

1817:    Output Parameters:
1818: +  nsolve - number of returned KSP objects
1819: -  ksp - array of linear solver object

1821:    Notes:
1822:    The number of KSP objects is equal to the number of shifts provided by the user,
1823:    or 1 if the user did not provide shifts.

1825:    Level: advanced

1827: .seealso: NEPNLEIGSSetRKShifts()
1828: @*/
1829: PetscErrorCode NEPNLEIGSGetKSPs(NEP nep,PetscInt *nsolve,KSP **ksp)
1830: {

1835:   PetscUseMethod(nep,"NEPNLEIGSGetKSPs_C",(NEP,PetscInt*,KSP**),(nep,nsolve,ksp));
1836:   return(0);
1837: }

1839: static PetscErrorCode NEPNLEIGSSetFullBasis_NLEIGS(NEP nep,PetscBool fullbasis)
1840: {
1841:   NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;

1844:   if (fullbasis!=ctx->fullbasis) {
1845:     ctx->fullbasis = fullbasis;
1846:     nep->state     = NEP_STATE_INITIAL;
1847:     nep->useds     = PetscNot(fullbasis);
1848:   }
1849:   return(0);
1850: }

1852: /*@
1853:    NEPNLEIGSSetFullBasis - Choose between TOAR-basis (default) and full-basis
1854:    variants of the NLEIGS method.

1856:    Logically Collective on nep

1858:    Input Parameters:
1859: +  nep  - the nonlinear eigensolver context
1860: -  fullbasis - true if the full-basis variant must be selected

1862:    Options Database Key:
1863: .  -nep_nleigs_full_basis - Sets the full-basis flag

1865:    Notes:
1866:    The default is to use a compact representation of the Krylov basis, that is,
1867:    V = (I otimes U) S, with a tensor BV. This behaviour can be changed so that
1868:    the full basis V is explicitly stored and operated with. This variant is more
1869:    expensive in terms of memory and computation, but is necessary in some cases,
1870:    particularly for two-sided computations, see NEPSetTwoSided().

1872:    In the full-basis variant, the NLEIGS solver uses an EPS object to explicitly
1873:    solve the linearized eigenproblem, see NEPNLEIGSGetEPS().

1875:    Level: advanced

1877: .seealso: NEPNLEIGSGetFullBasis(), NEPNLEIGSGetEPS(), NEPSetTwoSided(), BVCreateTensor()
1878: @*/
1879: PetscErrorCode NEPNLEIGSSetFullBasis(NEP nep,PetscBool fullbasis)
1880: {

1886:   PetscTryMethod(nep,"NEPNLEIGSSetFullBasis_C",(NEP,PetscBool),(nep,fullbasis));
1887:   return(0);
1888: }

1890: static PetscErrorCode NEPNLEIGSGetFullBasis_NLEIGS(NEP nep,PetscBool *fullbasis)
1891: {
1892:   NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;

1895:   *fullbasis = ctx->fullbasis;
1896:   return(0);
1897: }

1899: /*@
1900:    NEPNLEIGSGetFullBasis - Gets the flag that indicates if NLEIGS is using the
1901:    full-basis variant.

1903:    Not Collective

1905:    Input Parameter:
1906: .  nep - the nonlinear eigensolver context

1908:    Output Parameter:
1909: .  fullbasis - the flag

1911:    Level: advanced

1913: .seealso: NEPNLEIGSSetFullBasis()
1914: @*/
1915: PetscErrorCode NEPNLEIGSGetFullBasis(NEP nep,PetscBool *fullbasis)
1916: {

1922:   PetscUseMethod(nep,"NEPNLEIGSGetFullBasis_C",(NEP,PetscBool*),(nep,fullbasis));
1923:   return(0);
1924: }

1926: #define SHIFTMAX 30

1928: PetscErrorCode NEPSetFromOptions_NLEIGS(PetscOptionItems *PetscOptionsObject,NEP nep)
1929: {
1931:   NEP_NLEIGS     *ctx = (NEP_NLEIGS*)nep->data;
1932:   PetscInt       i=0,k;
1933:   PetscBool      flg1,flg2,b;
1934:   PetscReal      r;
1935:   PetscScalar    array[SHIFTMAX];

1938:   PetscOptionsHead(PetscOptionsObject,"NEP NLEIGS Options");

1940:     PetscOptionsReal("-nep_nleigs_restart","Proportion of vectors kept after restart","NEPNLEIGSSetRestart",0.5,&r,&flg1);
1941:     if (flg1) { NEPNLEIGSSetRestart(nep,r); }

1943:     PetscOptionsBool("-nep_nleigs_locking","Choose between locking and non-locking variants","NEPNLEIGSSetLocking",PETSC_FALSE,&b,&flg1);
1944:     if (flg1) { NEPNLEIGSSetLocking(nep,b); }

1946:     PetscOptionsBool("-nep_nleigs_full_basis","Choose between TOAR and full-basis variants","NEPNLEIGSSetFullBasis",PETSC_FALSE,&b,&flg1);
1947:     if (flg1) { NEPNLEIGSSetFullBasis(nep,b); }

1949:     NEPNLEIGSGetInterpolation(nep,&r,&i);
1950:     if (!i) i = PETSC_DEFAULT;
1951:     PetscOptionsInt("-nep_nleigs_interpolation_degree","Maximum number of terms for interpolation via divided differences","NEPNLEIGSSetInterpolation",i,&i,&flg1);
1952:     PetscOptionsReal("-nep_nleigs_interpolation_tol","Tolerance for interpolation via divided differences","NEPNLEIGSSetInterpolation",r,&r,&flg2);
1953:     if (flg1 || flg2) { NEPNLEIGSSetInterpolation(nep,r,i); }

1955:     k = SHIFTMAX;
1956:     for (i=0;i<k;i++) array[i] = 0;
1957:     PetscOptionsScalarArray("-nep_nleigs_rk_shifts","Shifts for Rational Krylov","NEPNLEIGSSetRKShifts",array,&k,&flg1);
1958:     if (flg1) { NEPNLEIGSSetRKShifts(nep,k,array); }

1960:   PetscOptionsTail();

1962:   if (!ctx->ksp) { NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp); }
1963:   for (i=0;i<ctx->nshiftsw;i++) {
1964:     KSPSetFromOptions(ctx->ksp[i]);
1965:   }

1967:   if (ctx->fullbasis) {
1968:     if (!ctx->eps) { NEPNLEIGSGetEPS(nep,&ctx->eps); }
1969:     EPSSetFromOptions(ctx->eps);
1970:   }
1971:   return(0);
1972: }

1974: PetscErrorCode NEPView_NLEIGS(NEP nep,PetscViewer viewer)
1975: {
1977:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
1978:   PetscBool      isascii;
1979:   PetscInt       i;
1980:   char           str[50];

1983:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
1984:   if (isascii) {
1985:     PetscViewerASCIIPrintf(viewer,"  %d%% of basis vectors kept after restart\n",(int)(100*ctx->keep));
1986:     if (ctx->fullbasis) {
1987:       PetscViewerASCIIPrintf(viewer,"  using the full-basis variant\n");
1988:     } else {
1989:       PetscViewerASCIIPrintf(viewer,"  using the %slocking variant\n",ctx->lock?"":"non-");
1990:     }
1991:     PetscViewerASCIIPrintf(viewer,"  divided difference terms: used=%D, max=%D\n",ctx->nmat,ctx->ddmaxit);
1992:     PetscViewerASCIIPrintf(viewer,"  tolerance for divided difference convergence: %g\n",(double)ctx->ddtol);
1993:     if (ctx->nshifts) {
1994:       PetscViewerASCIIPrintf(viewer,"  RK shifts: ");
1995:       PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1996:       for (i=0;i<ctx->nshifts;i++) {
1997:         SlepcSNPrintfScalar(str,sizeof(str),ctx->shifts[i],PETSC_FALSE);
1998:         PetscViewerASCIIPrintf(viewer,"%s%s",str,(i<ctx->nshifts-1)?",":"");
1999:       }
2000:       PetscViewerASCIIPrintf(viewer,"\n");
2001:       PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
2002:     }
2003:     if (!ctx->ksp) { NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp); }
2004:     PetscViewerASCIIPushTab(viewer);
2005:     KSPView(ctx->ksp[0],viewer);
2006:     PetscViewerASCIIPopTab(viewer);
2007:     if (ctx->fullbasis) {
2008:       if (!ctx->eps) { NEPNLEIGSGetEPS(nep,&ctx->eps); }
2009:       PetscViewerASCIIPushTab(viewer);
2010:       EPSView(ctx->eps,viewer);
2011:       PetscViewerASCIIPopTab(viewer);
2012:     }
2013:   }
2014:   return(0);
2015: }

2017: PetscErrorCode NEPReset_NLEIGS(NEP nep)
2018: {
2020:   PetscInt       k;
2021:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;

2024:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
2025:     PetscFree(ctx->coeffD);
2026:   } else {
2027:     for (k=0;k<ctx->nmat;k++) { MatDestroy(&ctx->D[k]); }
2028:   }
2029:   PetscFree4(ctx->s,ctx->xi,ctx->beta,ctx->D);
2030:   for (k=0;k<ctx->nshiftsw;k++) { KSPReset(ctx->ksp[k]); }
2031:   if (ctx->vrn) {
2032:     VecDestroy(&ctx->vrn);
2033:   }
2034:   if (ctx->fullbasis) {
2035:     MatDestroy(&ctx->A);
2036:     EPSReset(ctx->eps);
2037:     for (k=0;k<4;k++) { VecDestroy(&ctx->w[k]); }
2038:   }
2039:   return(0);
2040: }

2042: PetscErrorCode NEPDestroy_NLEIGS(NEP nep)
2043: {
2045:   PetscInt       k;
2046:   NEP_NLEIGS     *ctx = (NEP_NLEIGS*)nep->data;

2049:   BVDestroy(&ctx->V);
2050:   for (k=0;k<ctx->nshiftsw;k++) { KSPDestroy(&ctx->ksp[k]); }
2051:   PetscFree(ctx->ksp);
2052:   if (ctx->nshifts) { PetscFree(ctx->shifts); }
2053:   if (ctx->fullbasis) { EPSDestroy(&ctx->eps); }
2054:   PetscFree(nep->data);
2055:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetSingularitiesFunction_C",NULL);
2056:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetSingularitiesFunction_C",NULL);
2057:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRestart_C",NULL);
2058:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRestart_C",NULL);
2059:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetLocking_C",NULL);
2060:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetLocking_C",NULL);
2061:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetInterpolation_C",NULL);
2062:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetInterpolation_C",NULL);
2063:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRKShifts_C",NULL);
2064:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRKShifts_C",NULL);
2065:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetKSPs_C",NULL);
2066:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetFullBasis_C",NULL);
2067:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetFullBasis_C",NULL);
2068:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetEPS_C",NULL);
2069:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetEPS_C",NULL);
2070:   return(0);
2071: }

2073: SLEPC_EXTERN PetscErrorCode NEPCreate_NLEIGS(NEP nep)
2074: {
2076:   NEP_NLEIGS     *ctx;

2079:   PetscNewLog(nep,&ctx);
2080:   nep->data  = (void*)ctx;
2081:   ctx->lock  = PETSC_TRUE;
2082:   ctx->ddtol = PETSC_DEFAULT;

2084:   nep->useds = PETSC_TRUE;

2086:   nep->ops->setup          = NEPSetUp_NLEIGS;
2087:   nep->ops->setfromoptions = NEPSetFromOptions_NLEIGS;
2088:   nep->ops->view           = NEPView_NLEIGS;
2089:   nep->ops->destroy        = NEPDestroy_NLEIGS;
2090:   nep->ops->reset          = NEPReset_NLEIGS;

2092:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetSingularitiesFunction_C",NEPNLEIGSSetSingularitiesFunction_NLEIGS);
2093:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetSingularitiesFunction_C",NEPNLEIGSGetSingularitiesFunction_NLEIGS);
2094:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRestart_C",NEPNLEIGSSetRestart_NLEIGS);
2095:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRestart_C",NEPNLEIGSGetRestart_NLEIGS);
2096:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetLocking_C",NEPNLEIGSSetLocking_NLEIGS);
2097:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetLocking_C",NEPNLEIGSGetLocking_NLEIGS);
2098:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetInterpolation_C",NEPNLEIGSSetInterpolation_NLEIGS);
2099:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetInterpolation_C",NEPNLEIGSGetInterpolation_NLEIGS);
2100:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRKShifts_C",NEPNLEIGSSetRKShifts_NLEIGS);
2101:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRKShifts_C",NEPNLEIGSGetRKShifts_NLEIGS);
2102:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetKSPs_C",NEPNLEIGSGetKSPs_NLEIGS);
2103:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetFullBasis_C",NEPNLEIGSSetFullBasis_NLEIGS);
2104:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetFullBasis_C",NEPNLEIGSGetFullBasis_NLEIGS);
2105:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetEPS_C",NEPNLEIGSSetEPS_NLEIGS);
2106:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetEPS_C",NEPNLEIGSGetEPS_NLEIGS);
2107:   return(0);
2108: }