Actual source code: nleigs.c

slepc-3.17.0 2022-03-31
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain

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

 13:    Method: NLEIGS

 15:    Algorithm:

 17:        Fully rational Krylov method for nonlinear eigenvalue problems.

 19:    References:

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

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

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

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

 56: /* Computes the roots of a polynomial */
 57: static PetscErrorCode NEPNLEIGSAuxiliarPRootFinder(PetscInt deg,PetscScalar *polcoeffs,PetscScalar *wr,PetscScalar *wi,PetscBool *avail)
 58: {
 59:   PetscScalar    *C;
 60:   PetscBLASInt   n_,lwork;
 61:   PetscInt       i;
 62: #if defined(PETSC_USE_COMPLEX)
 63:   PetscReal      *rwork=NULL;
 64: #endif
 65:   PetscScalar    *work;
 66:   PetscBLASInt   info;

 68:   *avail = PETSC_TRUE;
 69:   if (deg>0) {
 70:     PetscCalloc1(deg*deg,&C);
 71:     PetscBLASIntCast(deg,&n_);
 72:     for (i=0;i<deg-1;i++) {
 73:       C[(deg+1)*i+1]   = 1.0;
 74:       C[(deg-1)*deg+i] = -polcoeffs[deg-i]/polcoeffs[0];
 75:     }
 76:     C[deg*deg+-1] = -polcoeffs[1]/polcoeffs[0];
 77:     PetscBLASIntCast(3*deg,&lwork);

 79:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
 80: #if !defined(PETSC_USE_COMPLEX)
 81:     PetscMalloc1(lwork,&work);
 82:     PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_("N","N",&n_,C,&n_,wr,wi,NULL,&n_,NULL,&n_,work,&lwork,&info));
 83:     if (info) *avail = PETSC_FALSE;
 84:     PetscFree(work);
 85: #else
 86:     PetscMalloc2(2*deg,&rwork,lwork,&work);
 87:     PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_("N","N",&n_,C,&n_,wr,NULL,&n_,NULL,&n_,work,&lwork,rwork,&info));
 88:     if (info) *avail = PETSC_FALSE;
 89:     PetscFree2(rwork,work);
 90: #endif
 91:     PetscFPTrapPop();
 92:     PetscFree(C);
 93:   }
 94:   PetscFunctionReturn(0);
 95: }

 97: static PetscErrorCode NEPNLEIGSAuxiliarRmDuplicates(PetscInt nin,PetscScalar *pin,PetscInt *nout,PetscScalar *pout,PetscInt max)
 98: {
 99:   PetscInt i,j;

101:   for (i=0;i<nin;i++) {
102:     if (max && *nout>=max) break;
103:     pout[(*nout)++] = pin[i];
104:     for (j=0;j<*nout-1;j++)
105:       if (PetscAbsScalar(pin[i]-pout[j])<PETSC_MACHINE_EPSILON*100) {
106:         (*nout)--;
107:         break;
108:       }
109:   }
110:   PetscFunctionReturn(0);
111: }

113: static PetscErrorCode NEPNLEIGSFNSingularities(FN f,PetscInt *nisol,PetscScalar **isol,PetscBool *rational)
114: {
115:   FNCombineType  ctype;
116:   FN             f1,f2;
117:   PetscInt       i,nq,nisol1,nisol2;
118:   PetscScalar    *qcoeff,*wr,*wi,*isol1,*isol2;
119:   PetscBool      flg,avail,rat1,rat2;

121:   *rational = PETSC_FALSE;
122:   PetscObjectTypeCompare((PetscObject)f,FNRATIONAL,&flg);
123:   if (flg) {
124:     *rational = PETSC_TRUE;
125:     FNRationalGetDenominator(f,&nq,&qcoeff);
126:     if (nq>1) {
127:       PetscMalloc2(nq-1,&wr,nq-1,&wi);
128:       NEPNLEIGSAuxiliarPRootFinder(nq-1,qcoeff,wr,wi,&avail);
129:       if (avail) {
130:         PetscCalloc1(nq-1,isol);
131:         *nisol = 0;
132:         for (i=0;i<nq-1;i++)
133: #if !defined(PETSC_USE_COMPLEX)
134:           if (wi[i]==0)
135: #endif
136:             (*isol)[(*nisol)++] = wr[i];
137:         nq = *nisol; *nisol = 0;
138:         for (i=0;i<nq;i++) wr[i] = (*isol)[i];
139:         NEPNLEIGSAuxiliarRmDuplicates(nq,wr,nisol,*isol,0);
140:         PetscFree2(wr,wi);
141:       } else { *nisol=0; *isol = NULL; }
142:     } else { *nisol = 0; *isol = NULL; }
143:     PetscFree(qcoeff);
144:   }
145:   PetscObjectTypeCompare((PetscObject)f,FNCOMBINE,&flg);
146:   if (flg) {
147:     FNCombineGetChildren(f,&ctype,&f1,&f2);
148:     if (ctype != FN_COMBINE_COMPOSE && ctype != FN_COMBINE_DIVIDE) {
149:       NEPNLEIGSFNSingularities(f1,&nisol1,&isol1,&rat1);
150:       NEPNLEIGSFNSingularities(f2,&nisol2,&isol2,&rat2);
151:       if (nisol1+nisol2>0) {
152:         PetscCalloc1(nisol1+nisol2,isol);
153:         *nisol = 0;
154:         NEPNLEIGSAuxiliarRmDuplicates(nisol1,isol1,nisol,*isol,0);
155:         NEPNLEIGSAuxiliarRmDuplicates(nisol2,isol2,nisol,*isol,0);
156:       }
157:       *rational = (rat1&&rat2)?PETSC_TRUE:PETSC_FALSE;
158:       PetscFree(isol1);
159:       PetscFree(isol2);
160:     }
161:   }
162:   PetscFunctionReturn(0);
163: }

165: static PetscErrorCode NEPNLEIGSRationalSingularities(NEP nep,PetscInt *ndptx,PetscScalar *dxi,PetscBool *rational)
166: {
167:   PetscInt       nt,i,nisol;
168:   FN             f;
169:   PetscScalar    *isol;
170:   PetscBool      rat;

172:   *rational = PETSC_TRUE;
173:   *ndptx = 0;
174:   NEPGetSplitOperatorInfo(nep,&nt,NULL);
175:   for (i=0;i<nt;i++) {
176:     NEPGetSplitOperatorTerm(nep,i,NULL,&f);
177:     NEPNLEIGSFNSingularities(f,&nisol,&isol,&rat);
178:     if (nisol) {
179:       NEPNLEIGSAuxiliarRmDuplicates(nisol,isol,ndptx,dxi,0);
180:       PetscFree(isol);
181:     }
182:     *rational = ((*rational)&&rat)?PETSC_TRUE:PETSC_FALSE;
183:   }
184:   PetscFunctionReturn(0);
185: }

187: /*  Adaptive Anderson-Antoulas algorithm */
188: static PetscErrorCode NEPNLEIGSAAAComputation(NEP nep,PetscInt ndpt,PetscScalar *ds,PetscScalar *F,PetscInt *ndptx,PetscScalar *dxi)
189: {
190:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
191:   PetscScalar    mean=0.0,*z,*f,*C,*A,*VT,*work,*ww,szero=0.0,sone=1.0;
192:   PetscScalar    *N,*D;
193:   PetscReal      *S,norm,err,*R;
194:   PetscInt       i,k,j,idx=0,cont;
195:   PetscBLASInt   n_,m_,lda_,lwork,info,one=1;
196: #if defined(PETSC_USE_COMPLEX)
197:   PetscReal      *rwork;
198: #endif

200:   PetscBLASIntCast(8*ndpt,&lwork);
201:   PetscMalloc5(ndpt,&R,ndpt,&z,ndpt,&f,ndpt*ndpt,&C,ndpt,&ww);
202:   PetscMalloc6(ndpt*ndpt,&A,ndpt,&S,ndpt*ndpt,&VT,lwork,&work,ndpt,&D,ndpt,&N);
203: #if defined(PETSC_USE_COMPLEX)
204:   PetscMalloc1(8*ndpt,&rwork);
205: #endif
206:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
207:   norm = 0.0;
208:   for (i=0;i<ndpt;i++) {
209:     mean += F[i];
210:     norm = PetscMax(PetscAbsScalar(F[i]),norm);
211:   }
212:   mean /= ndpt;
213:   PetscBLASIntCast(ndpt,&lda_);
214:   for (i=0;i<ndpt;i++) R[i] = PetscAbsScalar(F[i]-mean);
215:   /* next support point */
216:   err = 0.0;
217:   for (i=0;i<ndpt;i++) if (R[i]>=err) {idx = i; err = R[i];}
218:   for (k=0;k<ndpt-1;k++) {
219:     z[k] = ds[idx]; f[k] = F[idx]; R[idx] = -1.0;
220:     /* next column of Cauchy matrix */
221:     for (i=0;i<ndpt;i++) {
222:       C[i+k*ndpt] = 1.0/(ds[i]-ds[idx]);
223:     }

225:     PetscArrayzero(A,ndpt*ndpt);
226:     cont = 0;
227:     for (i=0;i<ndpt;i++) {
228:       if (R[i]!=-1.0) {
229:         for (j=0;j<=k;j++)A[cont+j*ndpt] = C[i+j*ndpt]*F[i]-C[i+j*ndpt]*f[j];
230:         cont++;
231:       }
232:     }
233:     PetscBLASIntCast(cont,&m_);
234:     PetscBLASIntCast(k+1,&n_);
235: #if defined(PETSC_USE_COMPLEX)
236:     PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","A",&m_,&n_,A,&lda_,S,NULL,&lda_,VT,&lda_,work,&lwork,rwork,&info));
237: #else
238:     PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","A",&m_,&n_,A,&lda_,S,NULL,&lda_,VT,&lda_,work,&lwork,&info));
239: #endif
240:     SlepcCheckLapackInfo("gesvd",info);
241:     for (i=0;i<=k;i++) {
242:       ww[i] = PetscConj(VT[i*ndpt+k]);
243:       D[i] = ww[i]*f[i];
244:     }
245:     PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&lda_,&n_,&sone,C,&lda_,D,&one,&szero,N,&one));
246:     PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&lda_,&n_,&sone,C,&lda_,ww,&one,&szero,D,&one));
247:     for (i=0;i<ndpt;i++) if (R[i]>=0) R[i] = PetscAbsScalar(F[i]-N[i]/D[i]);
248:     /* next support point */
249:     err = 0.0;
250:     for (i=0;i<ndpt;i++) if (R[i]>=err) {idx = i; err = R[i];}
251:     if (err <= ctx->ddtol*norm) break;
252:   }

255:   /* poles */
256:   PetscArrayzero(C,ndpt*ndpt);
257:   PetscArrayzero(A,ndpt*ndpt);
258:   for (i=0;i<=k;i++) {
259:     C[i+ndpt*i] = 1.0;
260:     A[(i+1)*ndpt] = ww[i];
261:     A[i+1] = 1.0;
262:     A[i+1+(i+1)*ndpt] = z[i];
263:   }
264:   C[0] = 0.0; C[k+1+(k+1)*ndpt] = 1.0;
265:   n_++;
266: #if defined(PETSC_USE_COMPLEX)
267:   PetscStackCallBLAS("LAPACKggev",LAPACKggev_("N","N",&n_,A,&lda_,C,&lda_,D,N,NULL,&lda_,NULL,&lda_,work,&lwork,rwork,&info));
268: #else
269:   PetscStackCallBLAS("LAPACKggev",LAPACKggev_("N","N",&n_,A,&lda_,C,&lda_,D,VT,N,NULL,&lda_,NULL,&lda_,work,&lwork,&info));
270: #endif
271:   SlepcCheckLapackInfo("ggev",info);
272:   cont = 0.0;
273:   for (i=0;i<n_;i++) if (N[i]!=0.0) {
274:     dxi[cont++] = D[i]/N[i];
275:   }
276:   *ndptx = cont;
277:   PetscFPTrapPop();
278:   PetscFree5(R,z,f,C,ww);
279:   PetscFree6(A,S,VT,work,D,N);
280: #if defined(PETSC_USE_COMPLEX)
281:   PetscFree(rwork);
282: #endif
283:   PetscFunctionReturn(0);
284: }

286: /*  Singularities using Adaptive Anderson-Antoulas algorithm */
287: static PetscErrorCode NEPNLEIGSAAASingularities(NEP nep,PetscInt ndpt,PetscScalar *ds,PetscInt *ndptx,PetscScalar *dxi)
288: {
289:   Vec            u,v,w;
290:   PetscRandom    rand=NULL;
291:   PetscScalar    *F,*isol;
292:   PetscInt       i,k,nisol,nt;
293:   Mat            T;
294:   FN             f;

296:   PetscMalloc1(ndpt,&F);
297:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
298:     PetscMalloc1(ndpt,&isol);
299:     *ndptx = 0;
300:     NEPGetSplitOperatorInfo(nep,&nt,NULL);
301:     nisol = *ndptx;
302:     for (k=0;k<nt;k++) {
303:       NEPGetSplitOperatorTerm(nep,k,NULL,&f);
304:       for (i=0;i<ndpt;i++) FNEvaluateFunction(f,ds[i],&F[i]);
305:       NEPNLEIGSAAAComputation(nep,ndpt,ds,F,&nisol,isol);
306:       if (nisol) NEPNLEIGSAuxiliarRmDuplicates(nisol,isol,ndptx,dxi,ndpt);
307:     }
308:     PetscFree(isol);
309:   } else {
310:     MatCreateVecs(nep->function,&u,NULL);
311:     VecDuplicate(u,&v);
312:     VecDuplicate(u,&w);
313:     if (nep->V) BVGetRandomContext(nep->V,&rand);
314:     VecSetRandom(u,rand);
315:     VecNormalize(u,NULL);
316:     VecSetRandom(v,rand);
317:     VecNormalize(v,NULL);
318:     T = nep->function;
319:     for (i=0;i<ndpt;i++) {
320:       NEPComputeFunction(nep,ds[i],T,T);
321:       MatMult(T,v,w);
322:       VecDot(w,u,&F[i]);
323:     }
324:     NEPNLEIGSAAAComputation(nep,ndpt,ds,F,ndptx,dxi);
325:     VecDestroy(&u);
326:     VecDestroy(&v);
327:     VecDestroy(&w);
328:   }
329:   PetscFree(F);
330:   PetscFunctionReturn(0);
331: }

333: static PetscErrorCode NEPNLEIGSLejaBagbyPoints(NEP nep)
334: {
335:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
336:   PetscInt       i,k,ndpt=NDPOINTS,ndptx=NDPOINTS;
337:   PetscScalar    *ds,*dsi,*dxi,*nrs,*nrxi,*s=ctx->s,*xi=ctx->xi,*beta=ctx->beta;
338:   PetscReal      maxnrs,minnrxi;
339:   PetscBool      rational;
340: #if !defined(PETSC_USE_COMPLEX)
341:   PetscReal      a,b,h;
342: #endif

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

347:   /* Discretize the target region boundary */
348:   RGComputeContour(nep->rg,ndpt,ds,dsi);
349: #if !defined(PETSC_USE_COMPLEX)
350:   for (i=0;i<ndpt;i++) if (dsi[i]!=0.0) break;
351:   if (i<ndpt) {
353:     /* Select a segment in the real axis */
354:     RGComputeBoundingBox(nep->rg,&a,&b,NULL,NULL);
356:     h = (b-a)/ndpt;
357:     for (i=0;i<ndpt;i++) {ds[i] = a+h*i; dsi[i] = 0.0;}
358:   }
359: #endif
360:   /* Discretize the singularity region */
361:   if (ctx->computesingularities) (ctx->computesingularities)(nep,&ndptx,dxi,ctx->singularitiesctx);
362:   else {
363:     if (nep->problem_type==NEP_RATIONAL) {
364:       NEPNLEIGSRationalSingularities(nep,&ndptx,dxi,&rational);
366:     } else {
367:       /* AAA algorithm */
368:       NEPNLEIGSAAASingularities(nep,ndpt,ds,&ndptx,dxi);
369:     }
370:   }
371:   /* Look for Leja-Bagby points in the discretization sets */
372:   s[0]  = ds[0];
373:   xi[0] = (ndptx>0)?dxi[0]:PETSC_INFINITY;
375:   beta[0] = 1.0; /* scaling factors are also computed here */
376:   for (i=0;i<ndpt;i++) {
377:     nrs[i] = 1.0;
378:     nrxi[i] = 1.0;
379:   }
380:   for (k=1;k<ctx->ddmaxit;k++) {
381:     maxnrs = 0.0;
382:     minnrxi = PETSC_MAX_REAL;
383:     for (i=0;i<ndpt;i++) {
384:       nrs[i] *= ((ds[i]-s[k-1])/(1.0-ds[i]/xi[k-1]))/beta[k-1];
385:       if (PetscAbsScalar(nrs[i])>maxnrs) {maxnrs = PetscAbsScalar(nrs[i]); s[k] = ds[i];}
386:     }
387:     if (ndptx>k) {
388:       for (i=1;i<ndptx;i++) {
389:         nrxi[i] *= ((dxi[i]-s[k-1])/(1.0-dxi[i]/xi[k-1]))/beta[k-1];
390:         if (PetscAbsScalar(nrxi[i])<minnrxi) {minnrxi = PetscAbsScalar(nrxi[i]); xi[k] = dxi[i];}
391:       }
393:     } else xi[k] = PETSC_INFINITY;
394:     beta[k] = maxnrs;
395:   }
396:   PetscFree5(ds,dsi,dxi,nrs,nrxi);
397:   PetscFunctionReturn(0);
398: }

400: PetscErrorCode NEPNLEIGSEvalNRTFunct(NEP nep,PetscInt k,PetscScalar sigma,PetscScalar *b)
401: {
402:   NEP_NLEIGS  *ctx=(NEP_NLEIGS*)nep->data;
403:   PetscInt    i;
404:   PetscScalar *beta=ctx->beta,*s=ctx->s,*xi=ctx->xi;

406:   b[0] = 1.0/beta[0];
407:   for (i=0;i<k;i++) {
408:     b[i+1] = ((sigma-s[i])*b[i])/(beta[i+1]*(1.0-sigma/xi[i]));
409:   }
410:   PetscFunctionReturn(0);
411: }

413: static PetscErrorCode MatMult_Fun(Mat A,Vec x,Vec y)
414: {
415:   NEP_NLEIGS_MATSHELL *ctx;
416:   PetscInt            i;

419:   MatShellGetContext(A,&ctx);
420:   MatMult(ctx->A[0],x,y);
421:   if (ctx->coeff[0]!=1.0) VecScale(y,ctx->coeff[0]);
422:   for (i=1;i<ctx->nmat;i++) {
423:     MatMult(ctx->A[i],x,ctx->t);
424:     VecAXPY(y,ctx->coeff[i],ctx->t);
425:   }
426:   PetscFunctionReturn(0);
427: }

429: static PetscErrorCode MatMultTranspose_Fun(Mat A,Vec x,Vec y)
430: {
431:   NEP_NLEIGS_MATSHELL *ctx;
432:   PetscInt            i;

435:   MatShellGetContext(A,&ctx);
436:   MatMultTranspose(ctx->A[0],x,y);
437:   if (ctx->coeff[0]!=1.0) VecScale(y,ctx->coeff[0]);
438:   for (i=1;i<ctx->nmat;i++) {
439:     MatMultTranspose(ctx->A[i],x,ctx->t);
440:     VecAXPY(y,ctx->coeff[i],ctx->t);
441:   }
442:   PetscFunctionReturn(0);
443: }

445: static PetscErrorCode MatGetDiagonal_Fun(Mat A,Vec diag)
446: {
447:   NEP_NLEIGS_MATSHELL *ctx;
448:   PetscInt            i;

451:   MatShellGetContext(A,&ctx);
452:   MatGetDiagonal(ctx->A[0],diag);
453:   if (ctx->coeff[0]!=1.0) VecScale(diag,ctx->coeff[0]);
454:   for (i=1;i<ctx->nmat;i++) {
455:     MatGetDiagonal(ctx->A[i],ctx->t);
456:     VecAXPY(diag,ctx->coeff[i],ctx->t);
457:   }
458:   PetscFunctionReturn(0);
459: }

461: static PetscErrorCode MatDuplicate_Fun(Mat A,MatDuplicateOption op,Mat *B)
462: {
463:   PetscInt            m,n,M,N,i;
464:   NEP_NLEIGS_MATSHELL *ctxnew,*ctx;
465:   void                (*fun)(void);

468:   MatShellGetContext(A,&ctx);
469:   PetscNew(&ctxnew);
470:   ctxnew->nmat = ctx->nmat;
471:   ctxnew->maxnmat = ctx->maxnmat;
472:   PetscMalloc2(ctxnew->maxnmat,&ctxnew->A,ctxnew->maxnmat,&ctxnew->coeff);
473:   for (i=0;i<ctx->nmat;i++) {
474:     PetscObjectReference((PetscObject)ctx->A[i]);
475:     ctxnew->A[i] = ctx->A[i];
476:     ctxnew->coeff[i] = ctx->coeff[i];
477:   }
478:   MatGetSize(ctx->A[0],&M,&N);
479:   MatGetLocalSize(ctx->A[0],&m,&n);
480:   VecDuplicate(ctx->t,&ctxnew->t);
481:   MatCreateShell(PetscObjectComm((PetscObject)A),m,n,M,N,(void*)ctxnew,B);
482:   MatShellSetManageScalingShifts(*B);
483:   MatShellGetOperation(A,MATOP_MULT,&fun);
484:   MatShellSetOperation(*B,MATOP_MULT,fun);
485:   MatShellGetOperation(A,MATOP_MULT_TRANSPOSE,&fun);
486:   MatShellSetOperation(*B,MATOP_MULT_TRANSPOSE,fun);
487:   MatShellGetOperation(A,MATOP_GET_DIAGONAL,&fun);
488:   MatShellSetOperation(*B,MATOP_GET_DIAGONAL,fun);
489:   MatShellGetOperation(A,MATOP_DUPLICATE,&fun);
490:   MatShellSetOperation(*B,MATOP_DUPLICATE,fun);
491:   MatShellGetOperation(A,MATOP_DESTROY,&fun);
492:   MatShellSetOperation(*B,MATOP_DESTROY,fun);
493:   MatShellGetOperation(A,MATOP_AXPY,&fun);
494:   MatShellSetOperation(*B,MATOP_AXPY,fun);
495:   PetscFunctionReturn(0);
496: }

498: static PetscErrorCode MatDestroy_Fun(Mat A)
499: {
500:   NEP_NLEIGS_MATSHELL *ctx;
501:   PetscInt            i;

504:   if (A) {
505:     MatShellGetContext(A,&ctx);
506:     for (i=0;i<ctx->nmat;i++) MatDestroy(&ctx->A[i]);
507:     VecDestroy(&ctx->t);
508:     PetscFree2(ctx->A,ctx->coeff);
509:     PetscFree(ctx);
510:   }
511:   PetscFunctionReturn(0);
512: }

514: static PetscErrorCode MatAXPY_Fun(Mat Y,PetscScalar a,Mat X,MatStructure str)
515: {
516:   NEP_NLEIGS_MATSHELL *ctxY,*ctxX;
517:   PetscInt            i,j;
518:   PetscBool           found;

521:   MatShellGetContext(Y,&ctxY);
522:   MatShellGetContext(X,&ctxX);
523:   for (i=0;i<ctxX->nmat;i++) {
524:     found = PETSC_FALSE;
525:     for (j=0;!found&&j<ctxY->nmat;j++) {
526:       if (ctxX->A[i]==ctxY->A[j]) {
527:         found = PETSC_TRUE;
528:         ctxY->coeff[j] += a*ctxX->coeff[i];
529:       }
530:     }
531:     if (!found) {
532:       ctxY->coeff[ctxY->nmat] = a*ctxX->coeff[i];
533:       ctxY->A[ctxY->nmat++] = ctxX->A[i];
534:       PetscObjectReference((PetscObject)ctxX->A[i]);
535:     }
536:   }
537:   PetscFunctionReturn(0);
538: }

540: static PetscErrorCode MatScale_Fun(Mat M,PetscScalar a)
541: {
542:   NEP_NLEIGS_MATSHELL *ctx;
543:   PetscInt            i;

546:   MatShellGetContext(M,&ctx);
547:   for (i=0;i<ctx->nmat;i++) ctx->coeff[i] *= a;
548:   PetscFunctionReturn(0);
549: }

551: static PetscErrorCode NLEIGSMatToMatShellArray(Mat A,Mat *Ms,PetscInt maxnmat)
552: {
553:   NEP_NLEIGS_MATSHELL *ctx;
554:   PetscInt            m,n,M,N;
555:   PetscBool           has;

557:   MatHasOperation(A,MATOP_DUPLICATE,&has);
559:   PetscNew(&ctx);
560:   ctx->maxnmat = maxnmat;
561:   PetscMalloc2(ctx->maxnmat,&ctx->A,ctx->maxnmat,&ctx->coeff);
562:   MatDuplicate(A,MAT_COPY_VALUES,&ctx->A[0]);
563:   ctx->nmat = 1;
564:   ctx->coeff[0] = 1.0;
565:   MatCreateVecs(A,&ctx->t,NULL);
566:   MatGetSize(A,&M,&N);
567:   MatGetLocalSize(A,&m,&n);
568:   MatCreateShell(PetscObjectComm((PetscObject)A),m,n,M,N,(void*)ctx,Ms);
569:   MatShellSetManageScalingShifts(*Ms);
570:   MatShellSetOperation(*Ms,MATOP_MULT,(void(*)(void))MatMult_Fun);
571:   MatShellSetOperation(*Ms,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_Fun);
572:   MatShellSetOperation(*Ms,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Fun);
573:   MatShellSetOperation(*Ms,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_Fun);
574:   MatShellSetOperation(*Ms,MATOP_DESTROY,(void(*)(void))MatDestroy_Fun);
575:   MatShellSetOperation(*Ms,MATOP_AXPY,(void(*)(void))MatAXPY_Fun);
576:   MatShellSetOperation(*Ms,MATOP_SCALE,(void(*)(void))MatScale_Fun);
577:   PetscFunctionReturn(0);
578: }

580: /*
581:    MatIsShellAny - returns true if any of the n matrices is a shell matrix
582:  */
583: static PetscErrorCode MatIsShellAny(Mat *A,PetscInt n,PetscBool *shell)
584: {
585:   PetscInt       i;
586:   PetscBool      flg;

588:   *shell = PETSC_FALSE;
589:   for (i=0;i<n;i++) {
590:     MatIsShell(A[i],&flg);
591:     if (flg) { *shell = PETSC_TRUE; break; }
592:   }
593:   PetscFunctionReturn(0);
594: }

596: static PetscErrorCode NEPNLEIGSDividedDifferences_split(NEP nep)
597: {
599:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
600:   PetscInt       k,j,i,maxnmat,nmax;
601:   PetscReal      norm0,norm,*matnorm;
602:   PetscScalar    *s=ctx->s,*beta=ctx->beta,*xi=ctx->xi,*b,alpha,*coeffs,*pK,*pH,sone=1.0;
603:   Mat            T,P,Ts,K,H;
604:   PetscBool      shell,hasmnorm=PETSC_FALSE,matrix=PETSC_TRUE;
605:   PetscBLASInt   n_;

607:   nmax = ctx->ddmaxit;
608:   PetscMalloc1(nep->nt*nmax,&ctx->coeffD);
609:   PetscMalloc3(nmax+1,&b,nmax+1,&coeffs,nep->nt,&matnorm);
610:   for (j=0;j<nep->nt;j++) {
611:     MatHasOperation(nep->A[j],MATOP_NORM,&hasmnorm);
612:     if (!hasmnorm) break;
613:     MatNorm(nep->A[j],NORM_INFINITY,matnorm+j);
614:   }
615:   /* Try matrix functions scheme */
616:   PetscCalloc2(nmax*nmax,&pK,nmax*nmax,&pH);
617:   for (i=0;i<nmax-1;i++) {
618:     pK[(nmax+1)*i]   = 1.0;
619:     pK[(nmax+1)*i+1] = beta[i+1]/xi[i];
620:     pH[(nmax+1)*i]   = s[i];
621:     pH[(nmax+1)*i+1] = beta[i+1];
622:   }
623:   pH[nmax*nmax-1] = s[nmax-1];
624:   pK[nmax*nmax-1] = 1.0;
625:   PetscBLASIntCast(nmax,&n_);
626:   PetscStackCallBLAS("BLAStrsm",BLAStrsm_("R","L","N","U",&n_,&n_,&sone,pK,&n_,pH,&n_));
627:   /* The matrix to be used is in H. K will be a work-space matrix */
628:   MatCreateSeqDense(PETSC_COMM_SELF,nmax,nmax,pH,&H);
629:   MatCreateSeqDense(PETSC_COMM_SELF,nmax,nmax,pK,&K);
630:   for (j=0;matrix&&j<nep->nt;j++) {
631:     PetscPushErrorHandler(PetscReturnErrorHandler,NULL);
632:     ierr = FNEvaluateFunctionMat(nep->f[j],H,K);
633:     PetscPopErrorHandler();
634:     if (!ierr) {
635:       for (i=0;i<nmax;i++) ctx->coeffD[j+i*nep->nt] = pK[i]*beta[0];
636:     } else {
637:       matrix = PETSC_FALSE;
638:       PetscFPTrapPop();
639:     }
640:   }
641:   MatDestroy(&H);
642:   MatDestroy(&K);
643:   if (!matrix) {
644:     for (j=0;j<nep->nt;j++) {
645:       FNEvaluateFunction(nep->f[j],s[0],ctx->coeffD+j);
646:       ctx->coeffD[j] *= beta[0];
647:     }
648:   }
649:   if (hasmnorm) {
650:     norm0 = 0.0;
651:     for (j=0;j<nep->nt;j++) norm0 += matnorm[j]*PetscAbsScalar(ctx->coeffD[j]);
652:   } else {
653:     norm0 = 0.0;
654:     for (j=0;j<nep->nt;j++) norm0 = PetscMax(PetscAbsScalar(ctx->coeffD[j]),norm0);
655:   }
656:   ctx->nmat = ctx->ddmaxit;
657:   for (k=1;k<ctx->ddmaxit;k++) {
658:     if (!matrix) {
659:       NEPNLEIGSEvalNRTFunct(nep,k,s[k],b);
660:       for (i=0;i<nep->nt;i++) {
661:         FNEvaluateFunction(nep->f[i],s[k],ctx->coeffD+k*nep->nt+i);
662:         for (j=0;j<k;j++) {
663:           ctx->coeffD[k*nep->nt+i] -= b[j]*ctx->coeffD[i+nep->nt*j];
664:         }
665:         ctx->coeffD[k*nep->nt+i] /= b[k];
666:       }
667:     }
668:     if (hasmnorm) {
669:       norm = 0.0;
670:       for (j=0;j<nep->nt;j++) norm += matnorm[j]*PetscAbsScalar(ctx->coeffD[k*nep->nt+j]);
671:     } else {
672:       norm = 0.0;
673:       for (j=0;j<nep->nt;j++) norm = PetscMax(PetscAbsScalar(ctx->coeffD[k*nep->nt+j]),norm);
674:     }
675:     if (k>1 && norm/norm0 < ctx->ddtol) {
676:       ctx->nmat = k+1;
677:       break;
678:     }
679:   }
680:   if (!ctx->ksp) NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp);
681:   MatIsShellAny(nep->A,nep->nt,&shell);
682:   maxnmat = PetscMax(ctx->ddmaxit,nep->nt);
683:   for (i=0;i<ctx->nshiftsw;i++) {
684:     NEPNLEIGSEvalNRTFunct(nep,ctx->nmat-1,ctx->shifts[i],coeffs);
685:     if (!shell) MatDuplicate(nep->A[0],MAT_COPY_VALUES,&T);
686:     else NLEIGSMatToMatShellArray(nep->A[0],&T,maxnmat);
687:     if (nep->P) { /* user-defined preconditioner */
688:       MatDuplicate(nep->P[0],MAT_COPY_VALUES,&P);
689:     } else P=T;
690:     alpha = 0.0;
691:     for (j=0;j<ctx->nmat;j++) alpha += coeffs[j]*ctx->coeffD[j*nep->nt];
692:     MatScale(T,alpha);
693:     if (nep->P) MatScale(P,alpha);
694:     for (k=1;k<nep->nt;k++) {
695:       alpha = 0.0;
696:       for (j=0;j<ctx->nmat;j++) alpha += coeffs[j]*ctx->coeffD[j*nep->nt+k];
697:       if (shell) NLEIGSMatToMatShellArray(nep->A[k],&Ts,maxnmat);
698:       MatAXPY(T,alpha,shell?Ts:nep->A[k],nep->mstr);
699:       if (nep->P) MatAXPY(P,alpha,nep->P[k],nep->mstrp);
700:       if (shell) MatDestroy(&Ts);
701:     }
702:     NEP_KSPSetOperators(ctx->ksp[i],T,P);
703:     KSPSetUp(ctx->ksp[i]);
704:     MatDestroy(&T);
705:     if (nep->P) MatDestroy(&P);
706:   }
707:   PetscFree3(b,coeffs,matnorm);
708:   PetscFree2(pK,pH);
709:   PetscFunctionReturn(0);
710: }

712: static PetscErrorCode NEPNLEIGSDividedDifferences_callback(NEP nep)
713: {
714:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
715:   PetscInt       k,j,i,maxnmat;
716:   PetscReal      norm0,norm;
717:   PetscScalar    *s=ctx->s,*beta=ctx->beta,*b,*coeffs;
718:   Mat            *D=ctx->D,*DP,T,P;
719:   PetscBool      shell,has,vec=PETSC_FALSE,precond=(nep->function_pre!=nep->function)?PETSC_TRUE:PETSC_FALSE;
720:   PetscRandom    rand=NULL;
721:   Vec            w[2];

723:   PetscMalloc2(ctx->ddmaxit+1,&b,ctx->ddmaxit+1,&coeffs);
724:   if (nep->V) BVGetRandomContext(nep->V,&rand);
725:   T = nep->function;
726:   P = nep->function_pre;
727:   NEPComputeFunction(nep,s[0],T,P);
728:   MatIsShell(T,&shell);
729:   maxnmat = PetscMax(ctx->ddmaxit,nep->nt);
730:   if (!shell) MatDuplicate(T,MAT_COPY_VALUES,&D[0]);
731:   else NLEIGSMatToMatShellArray(T,&D[0],maxnmat);
732:   if (beta[0]!=1.0) MatScale(D[0],1.0/beta[0]);
733:   MatHasOperation(D[0],MATOP_NORM,&has);
734:   if (has) MatNorm(D[0],NORM_FROBENIUS,&norm0);
735:   else {
736:     MatCreateVecs(D[0],NULL,&w[0]);
737:     VecDuplicate(w[0],&w[1]);
738:     VecDuplicate(w[0],&ctx->vrn);
739:     VecSetRandomNormal(ctx->vrn,rand,w[0],w[1]);
740:     VecNormalize(ctx->vrn,NULL);
741:     vec = PETSC_TRUE;
742:     MatNormEstimate(D[0],ctx->vrn,w[0],&norm0);
743:   }
744:   if (precond) {
745:     PetscMalloc1(ctx->ddmaxit,&DP);
746:     MatDuplicate(P,MAT_COPY_VALUES,&DP[0]);
747:   }
748:   ctx->nmat = ctx->ddmaxit;
749:   for (k=1;k<ctx->ddmaxit;k++) {
750:     NEPNLEIGSEvalNRTFunct(nep,k,s[k],b);
751:     NEPComputeFunction(nep,s[k],T,P);
752:     if (!shell) MatDuplicate(T,MAT_COPY_VALUES,&D[k]);
753:     else NLEIGSMatToMatShellArray(T,&D[k],maxnmat);
754:     for (j=0;j<k;j++) MatAXPY(D[k],-b[j],D[j],nep->mstr);
755:     MatScale(D[k],1.0/b[k]);
756:     MatHasOperation(D[k],MATOP_NORM,&has);
757:     if (has) MatNorm(D[k],NORM_FROBENIUS,&norm);
758:     else {
759:       if (!vec) {
760:         MatCreateVecs(D[k],NULL,&w[0]);
761:         VecDuplicate(w[0],&w[1]);
762:         VecDuplicate(w[0],&ctx->vrn);
763:         VecSetRandomNormal(ctx->vrn,rand,w[0],w[1]);
764:         VecNormalize(ctx->vrn,NULL);
765:         vec = PETSC_TRUE;
766:       }
767:       MatNormEstimate(D[k],ctx->vrn,w[0],&norm);
768:     }
769:     if (precond) {
770:       MatDuplicate(P,MAT_COPY_VALUES,&DP[k]);
771:       for (j=0;j<k;j++) MatAXPY(DP[k],-b[j],DP[j],nep->mstrp);
772:       MatScale(DP[k],1.0/b[k]);
773:     }
774:     if (k>1 && norm/norm0 < ctx->ddtol && k>1) {
775:       ctx->nmat = k+1;
776:       break;
777:     }
778:   }
779:   if (!ctx->ksp) NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp);
780:   for (i=0;i<ctx->nshiftsw;i++) {
781:     NEPNLEIGSEvalNRTFunct(nep,ctx->nmat-1,ctx->shifts[i],coeffs);
782:     MatDuplicate(D[0],MAT_COPY_VALUES,&T);
783:     if (coeffs[0]!=1.0) MatScale(T,coeffs[0]);
784:     for (j=1;j<ctx->nmat;j++) MatAXPY(T,coeffs[j],D[j],nep->mstr);
785:     if (precond) {
786:       MatDuplicate(DP[0],MAT_COPY_VALUES,&P);
787:       if (coeffs[0]!=1.0) MatScale(P,coeffs[0]);
788:       for (j=1;j<ctx->nmat;j++) MatAXPY(P,coeffs[j],DP[j],nep->mstrp);
789:     } else P=T;
790:     NEP_KSPSetOperators(ctx->ksp[i],T,P);
791:     KSPSetUp(ctx->ksp[i]);
792:     MatDestroy(&T);
793:   }
794:   PetscFree2(b,coeffs);
795:   if (vec) {
796:     VecDestroy(&w[0]);
797:     VecDestroy(&w[1]);
798:   }
799:   if (precond) {
800:     MatDestroy(&P);
801:     MatDestroyMatrices(ctx->nmat,&DP);
802:   }
803:   PetscFunctionReturn(0);
804: }

806: /*
807:    NEPKrylovConvergence - This is the analogue to EPSKrylovConvergence.
808: */
809: PetscErrorCode NEPNLEIGSKrylovConvergence(NEP nep,PetscBool getall,PetscInt kini,PetscInt nits,PetscReal betah,PetscScalar betak,PetscInt *kout,Vec *w)
810: {
811:   PetscInt       k,newk,marker,inside;
812:   PetscScalar    re,im;
813:   PetscReal      resnorm,tt;
814:   PetscBool      istrivial;
815:   NEP_NLEIGS     *ctx = (NEP_NLEIGS*)nep->data;

817:   RGIsTrivial(nep->rg,&istrivial);
818:   marker = -1;
819:   if (nep->trackall) getall = PETSC_TRUE;
820:   for (k=kini;k<kini+nits;k++) {
821:     /* eigenvalue */
822:     re = nep->eigr[k];
823:     im = nep->eigi[k];
824:     if (!istrivial) {
825:       if (!ctx->nshifts) NEPNLEIGSBackTransform((PetscObject)nep,1,&re,&im);
826:       RGCheckInside(nep->rg,1,&re,&im,&inside);
827:       if (marker==-1 && inside<0) marker = k;
828:     }
829:     newk = k;
830:     DSVectors(nep->ds,DS_MAT_X,&newk,&resnorm);
831:     tt = (ctx->nshifts)?SlepcAbsEigenvalue(betak-nep->eigr[k]*betah,nep->eigi[k]*betah):betah;
832:     resnorm *=  PetscAbsReal(tt);
833:     /* error estimate */
834:     (*nep->converged)(nep,nep->eigr[k],nep->eigi[k],resnorm,&nep->errest[k],nep->convergedctx);
835:     if (marker==-1 && nep->errest[k] >= nep->tol) marker = k;
836:     if (newk==k+1) {
837:       nep->errest[k+1] = nep->errest[k];
838:       k++;
839:     }
840:     if (marker!=-1 && !getall) break;
841:   }
842:   if (marker!=-1) k = marker;
843:   *kout = k;
844:   PetscFunctionReturn(0);
845: }

847: PetscErrorCode NEPSetUp_NLEIGS(NEP nep)
848: {
849:   PetscInt       k,in;
850:   PetscScalar    zero=0.0;
851:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
852:   SlepcSC        sc;
853:   PetscBool      istrivial;

855:   NEPSetDimensions_Default(nep,nep->nev,&nep->ncv,&nep->mpd);
857:   if (nep->max_it==PETSC_DEFAULT) nep->max_it = PetscMax(5000,2*nep->n/nep->ncv);
858:   if (!ctx->ddmaxit) ctx->ddmaxit = LBPOINTS;
859:   RGIsTrivial(nep->rg,&istrivial);
861:   if (!nep->which) nep->which = NEP_TARGET_MAGNITUDE;

864:   /* Initialize the NLEIGS context structure */
865:   k = ctx->ddmaxit;
866:   PetscMalloc4(k,&ctx->s,k,&ctx->xi,k,&ctx->beta,k,&ctx->D);
867:   nep->data = ctx;
868:   if (nep->tol==PETSC_DEFAULT) nep->tol = SLEPC_DEFAULT_TOL;
869:   if (ctx->ddtol==PETSC_DEFAULT) ctx->ddtol = nep->tol/10.0;
870:   if (!ctx->keep) ctx->keep = 0.5;

872:   /* Compute Leja-Bagby points and scaling values */
873:   NEPNLEIGSLejaBagbyPoints(nep);
874:   if (nep->problem_type!=NEP_RATIONAL) {
875:     RGCheckInside(nep->rg,1,&nep->target,&zero,&in);
877:   }

879:   /* Compute the divided difference matrices */
880:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) NEPNLEIGSDividedDifferences_split(nep);
881:   else NEPNLEIGSDividedDifferences_callback(nep);
882:   NEPAllocateSolution(nep,ctx->nmat-1);
883:   NEPSetWorkVecs(nep,4);
884:   if (!ctx->fullbasis) {
886:     /* set-up DS and transfer split operator functions */
887:     DSSetType(nep->ds,ctx->nshifts?DSGNHEP:DSNHEP);
888:     DSAllocate(nep->ds,nep->ncv+1);
889:     DSGetSlepcSC(nep->ds,&sc);
890:     if (!ctx->nshifts) sc->map = NEPNLEIGSBackTransform;
891:     DSSetExtraRow(nep->ds,PETSC_TRUE);
892:     sc->mapobj        = (PetscObject)nep;
893:     sc->rg            = nep->rg;
894:     sc->comparison    = nep->sc->comparison;
895:     sc->comparisonctx = nep->sc->comparisonctx;
896:     BVDestroy(&ctx->V);
897:     BVCreateTensor(nep->V,ctx->nmat-1,&ctx->V);
898:     nep->ops->solve          = NEPSolve_NLEIGS;
899:     nep->ops->computevectors = NEPComputeVectors_Schur;
900:   } else {
901:     NEPSetUp_NLEIGS_FullBasis(nep);
902:     nep->ops->solve          = NEPSolve_NLEIGS_FullBasis;
903:     nep->ops->computevectors = NULL;
904:   }
905:   PetscFunctionReturn(0);
906: }

908: /*
909:   Extend the TOAR basis by applying the the matrix operator
910:   over a vector which is decomposed on the TOAR way
911:   Input:
912:     - S,V: define the latest Arnoldi vector (nv vectors in V)
913:   Output:
914:     - t: new vector extending the TOAR basis
915:     - r: temporally coefficients to compute the TOAR coefficients
916:          for the new Arnoldi vector
917:   Workspace: t_ (two vectors)
918: */
919: static PetscErrorCode NEPTOARExtendBasis(NEP nep,PetscInt idxrktg,PetscScalar *S,PetscInt ls,PetscInt nv,BV W,BV V,Vec t,PetscScalar *r,PetscInt lr,Vec *t_)
920: {
921:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
922:   PetscInt       deg=ctx->nmat-1,k,j;
923:   Vec            v=t_[0],q=t_[1],w;
924:   PetscScalar    *beta=ctx->beta,*s=ctx->s,*xi=ctx->xi,*coeffs,sigma;

926:   if (!ctx->ksp) NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp);
927:   sigma = ctx->shifts[idxrktg];
928:   BVSetActiveColumns(nep->V,0,nv);
929:   PetscMalloc1(ctx->nmat,&coeffs);
931:   /* i-part stored in (i-1) position */
932:   for (j=0;j<nv;j++) {
933:     r[(deg-2)*lr+j] = (S[(deg-2)*ls+j]+(beta[deg-1]/xi[deg-2])*S[(deg-1)*ls+j])/(s[deg-2]-sigma);
934:   }
935:   BVSetActiveColumns(W,0,deg);
936:   BVGetColumn(W,deg-1,&w);
937:   BVMultVec(V,1.0/beta[deg],0,w,S+(deg-1)*ls);
938:   BVRestoreColumn(W,deg-1,&w);
939:   BVGetColumn(W,deg-2,&w);
940:   BVMultVec(V,1.0,0.0,w,r+(deg-2)*lr);
941:   BVRestoreColumn(W,deg-2,&w);
942:   for (k=deg-2;k>0;k--) {
944:     for (j=0;j<nv;j++) r[(k-1)*lr+j] = (S[(k-1)*ls+j]+(beta[k]/xi[k-1])*S[k*ls+j]-beta[k]*(1.0-sigma/xi[k-1])*r[(k)*lr+j])/(s[k-1]-sigma);
945:     BVGetColumn(W,k-1,&w);
946:     BVMultVec(V,1.0,0.0,w,r+(k-1)*lr);
947:     BVRestoreColumn(W,k-1,&w);
948:   }
949:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
950:     for (j=0;j<ctx->nmat-2;j++) coeffs[j] = ctx->coeffD[nep->nt*j];
951:     coeffs[ctx->nmat-2] = ctx->coeffD[nep->nt*(ctx->nmat-1)];
952:     BVMultVec(W,1.0,0.0,v,coeffs);
953:     MatMult(nep->A[0],v,q);
954:     for (k=1;k<nep->nt;k++) {
955:       for (j=0;j<ctx->nmat-2;j++) coeffs[j] = ctx->coeffD[nep->nt*j+k];
956:       coeffs[ctx->nmat-2] = ctx->coeffD[nep->nt*(ctx->nmat-1)+k];
957:       BVMultVec(W,1.0,0,v,coeffs);
958:       MatMult(nep->A[k],v,t);
959:       VecAXPY(q,1.0,t);
960:     }
961:     KSPSolve(ctx->ksp[idxrktg],q,t);
962:     VecScale(t,-1.0);
963:   } else {
964:     for (k=0;k<deg-1;k++) {
965:       BVGetColumn(W,k,&w);
966:       MatMult(ctx->D[k],w,q);
967:       BVRestoreColumn(W,k,&w);
968:       BVInsertVec(W,k,q);
969:     }
970:     BVGetColumn(W,deg-1,&w);
971:     MatMult(ctx->D[deg],w,q);
972:     BVRestoreColumn(W,k,&w);
973:     BVInsertVec(W,k,q);
974:     for (j=0;j<ctx->nmat-1;j++) coeffs[j] = 1.0;
975:     BVMultVec(W,1.0,0.0,q,coeffs);
976:     KSPSolve(ctx->ksp[idxrktg],q,t);
977:     VecScale(t,-1.0);
978:   }
979:   PetscFree(coeffs);
980:   PetscFunctionReturn(0);
981: }

983: /*
984:   Compute TOAR coefficients of the blocks of the new Arnoldi vector computed
985: */
986: static PetscErrorCode NEPTOARCoefficients(NEP nep,PetscScalar sigma,PetscInt nv,PetscScalar *S,PetscInt ls,PetscScalar *r,PetscInt lr,PetscScalar *x,PetscScalar *work)
987: {
988:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
989:   PetscInt       k,j,d=ctx->nmat-1;
990:   PetscScalar    *t=work;

992:   NEPNLEIGSEvalNRTFunct(nep,d-1,sigma,t);
993:   for (k=0;k<d-1;k++) {
994:     for (j=0;j<=nv;j++) r[k*lr+j] += t[k]*x[j];
995:   }
996:   for (j=0;j<=nv;j++) r[(d-1)*lr+j] = t[d-1]*x[j];
997:   PetscFunctionReturn(0);
998: }

1000: /*
1001:   Compute continuation vector coefficients for the Rational-Krylov run.
1002:   dim(work) >= (end-ini)*(end-ini+1) + end+1 + 2*(end-ini+1), dim(t) = end.
1003: */
1004: static PetscErrorCode NEPNLEIGS_RKcontinuation(NEP nep,PetscInt ini,PetscInt end,PetscScalar *K,PetscScalar *H,PetscInt ld,PetscScalar sigma,PetscScalar *S,PetscInt lds,PetscScalar *cont,PetscScalar *t,PetscScalar *work)
1005: {
1006:   PetscScalar    *x,*W,*tau,sone=1.0,szero=0.0;
1007:   PetscInt       i,j,n1,n,nwu=0;
1008:   PetscBLASInt   info,n_,n1_,one=1,dim,lds_;
1009:   NEP_NLEIGS     *ctx = (NEP_NLEIGS*)nep->data;

1011:   if (!ctx->nshifts || !end) {
1012:     t[0] = 1;
1013:     PetscArraycpy(cont,S+end*lds,lds);
1014:   } else {
1015:     n   = end-ini;
1016:     n1  = n+1;
1017:     x   = work+nwu;
1018:     nwu += end+1;
1019:     tau = work+nwu;
1020:     nwu += n;
1021:     W   = work+nwu;
1022:     nwu += n1*n;
1023:     for (j=ini;j<end;j++) {
1024:       for (i=ini;i<=end;i++) W[(j-ini)*n1+i-ini] = K[j*ld+i] -H[j*ld+i]*sigma;
1025:     }
1026:     PetscBLASIntCast(n,&n_);
1027:     PetscBLASIntCast(n1,&n1_);
1028:     PetscBLASIntCast(end+1,&dim);
1029:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1030:     PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&n1_,&n_,W,&n1_,tau,work+nwu,&n1_,&info));
1031:     SlepcCheckLapackInfo("geqrf",info);
1032:     for (i=0;i<end;i++) t[i] = 0.0;
1033:     t[end] = 1.0;
1034:     for (j=n-1;j>=0;j--) {
1035:       for (i=0;i<ini+j;i++) x[i] = 0.0;
1036:       x[ini+j] = 1.0;
1037:       for (i=j+1;i<n1;i++) x[i+ini] = W[i+n1*j];
1038:       tau[j] = PetscConj(tau[j]);
1039:       PetscStackCallBLAS("LAPACKlarf",LAPACKlarf_("L",&dim,&one,x,&one,tau+j,t,&dim,work+nwu));
1040:     }
1041:     PetscBLASIntCast(lds,&lds_);
1042:     PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&lds_,&n1_,&sone,S,&lds_,t,&one,&szero,cont,&one));
1043:     PetscFPTrapPop();
1044:   }
1045:   PetscFunctionReturn(0);
1046: }

1048: /*
1049:   Compute a run of Arnoldi iterations
1050: */
1051: PetscErrorCode NEPNLEIGSTOARrun(NEP nep,PetscScalar *K,PetscScalar *H,PetscInt ldh,BV W,PetscInt k,PetscInt *M,PetscBool *breakdown,Vec *t_)
1052: {
1053:   NEP_NLEIGS     *ctx = (NEP_NLEIGS*)nep->data;
1054:   PetscInt       i,j,m=*M,lwa,deg=ctx->nmat-1,lds,nqt,ld,l;
1055:   Vec            t;
1056:   PetscReal      norm;
1057:   PetscScalar    *x,*work,*tt,sigma,*cont,*S;
1058:   PetscBool      lindep;
1059:   Mat            MS;

1061:   BVTensorGetFactors(ctx->V,NULL,&MS);
1062:   MatDenseGetArray(MS,&S);
1063:   BVGetSizes(nep->V,NULL,NULL,&ld);
1064:   lds = ld*deg;
1065:   BVGetActiveColumns(nep->V,&l,&nqt);
1066:   lwa = PetscMax(ld,deg)+(m+1)*(m+1)+4*(m+1);
1067:   PetscMalloc4(ld,&x,lwa,&work,m+1,&tt,lds,&cont);
1068:   BVSetActiveColumns(ctx->V,0,m);
1069:   for (j=k;j<m;j++) {
1070:     sigma = ctx->shifts[(++(ctx->idxrk))%ctx->nshiftsw];

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

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

1080:     /* orthogonalize */
1081:     BVOrthogonalizeColumn(nep->V,nqt,x,&norm,&lindep);
1082:     if (!lindep) {
1083:       x[nqt] = norm;
1084:       BVScaleColumn(nep->V,nqt,1.0/norm);
1085:       nqt++;
1086:     } else x[nqt] = 0.0;

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

1090:     /* Level-2 orthogonalization */
1091:     BVOrthogonalizeColumn(ctx->V,j+1,H+j*ldh,&norm,breakdown);
1092:     H[j+1+ldh*j] = norm;
1093:     if (ctx->nshifts) {
1094:       for (i=0;i<=j;i++) K[i+ldh*j] = sigma*H[i+ldh*j] + tt[i];
1095:       K[j+1+ldh*j] = sigma*H[j+1+ldh*j];
1096:     }
1097:     if (*breakdown) {
1098:       *M = j+1;
1099:       break;
1100:     }
1101:     BVScaleColumn(ctx->V,j+1,1.0/norm);
1102:     BVSetActiveColumns(nep->V,l,nqt);
1103:   }
1104:   PetscFree4(x,work,tt,cont);
1105:   MatDenseRestoreArray(MS,&S);
1106:   BVTensorRestoreFactors(ctx->V,NULL,&MS);
1107:   PetscFunctionReturn(0);
1108: }

1110: PetscErrorCode NEPSolve_NLEIGS(NEP nep)
1111: {
1112:   NEP_NLEIGS        *ctx = (NEP_NLEIGS*)nep->data;
1113:   PetscInt          i,k=0,l,nv=0,ld,lds,ldds,nq;
1114:   PetscInt          deg=ctx->nmat-1,nconv=0,dsn,dsk;
1115:   PetscScalar       *H,*pU,*K,betak=0,*eigr,*eigi;
1116:   const PetscScalar *S;
1117:   PetscReal         betah;
1118:   PetscBool         falselock=PETSC_FALSE,breakdown=PETSC_FALSE;
1119:   BV                W;
1120:   Mat               MS,MQ,U;

1122:   if (ctx->lock) {
1123:     /* undocumented option to use a cheaper locking instead of the true locking */
1124:     PetscOptionsGetBool(NULL,NULL,"-nep_nleigs_falselocking",&falselock,NULL);
1125:   }

1127:   BVGetSizes(nep->V,NULL,NULL,&ld);
1128:   lds = deg*ld;
1129:   DSGetLeadingDimension(nep->ds,&ldds);
1130:   if (!ctx->nshifts) PetscMalloc2(nep->ncv,&eigr,nep->ncv,&eigi);
1131:   else { eigr = nep->eigr; eigi = nep->eigi; }
1132:   BVDuplicateResize(nep->V,PetscMax(nep->nt-1,ctx->nmat-1),&W);

1134:   /* clean projected matrix (including the extra-arrow) */
1135:   DSGetArray(nep->ds,DS_MAT_A,&H);
1136:   PetscArrayzero(H,ldds*ldds);
1137:   DSRestoreArray(nep->ds,DS_MAT_A,&H);
1138:   if (ctx->nshifts) {
1139:     DSGetArray(nep->ds,DS_MAT_B,&H);
1140:     PetscArrayzero(H,ldds*ldds);
1141:     DSRestoreArray(nep->ds,DS_MAT_B,&H);
1142:   }

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

1147:   /* Restart loop */
1148:   l = 0;
1149:   while (nep->reason == NEP_CONVERGED_ITERATING) {
1150:     nep->its++;

1152:     /* Compute an nv-step Krylov relation */
1153:     nv = PetscMin(nep->nconv+nep->mpd,nep->ncv);
1154:     if (ctx->nshifts) DSGetArray(nep->ds,DS_MAT_A,&K);
1155:     DSGetArray(nep->ds,ctx->nshifts?DS_MAT_B:DS_MAT_A,&H);
1156:     NEPNLEIGSTOARrun(nep,K,H,ldds,W,nep->nconv+l,&nv,&breakdown,nep->work);
1157:     betah = PetscAbsScalar(H[(nv-1)*ldds+nv]);
1158:     DSRestoreArray(nep->ds,ctx->nshifts?DS_MAT_B:DS_MAT_A,&H);
1159:     if (ctx->nshifts) {
1160:       betak = K[(nv-1)*ldds+nv];
1161:       DSRestoreArray(nep->ds,DS_MAT_A,&K);
1162:     }
1163:     DSSetDimensions(nep->ds,nv,nep->nconv,nep->nconv+l);
1164:     if (l==0) DSSetState(nep->ds,DS_STATE_INTERMEDIATE);
1165:     else DSSetState(nep->ds,DS_STATE_RAW);

1167:     /* Solve projected problem */
1168:     DSSolve(nep->ds,nep->eigr,nep->eigi);
1169:     DSSort(nep->ds,nep->eigr,nep->eigi,NULL,NULL,NULL);
1170:     DSUpdateExtraRow(nep->ds);
1171:     DSSynchronize(nep->ds,nep->eigr,nep->eigi);

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

1177:     /* Update l */
1178:     if (nep->reason != NEP_CONVERGED_ITERATING || breakdown) l = 0;
1179:     else {
1180:       l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
1181:       DSGetTruncateSize(nep->ds,k,nv,&l);
1182:       if (!breakdown) {
1183:         /* Prepare the Rayleigh quotient for restart */
1184:         DSGetDimensions(nep->ds,&dsn,NULL,&dsk,NULL);
1185:         DSSetDimensions(nep->ds,dsn,k,dsk);
1186:         DSTruncate(nep->ds,k+l,PETSC_FALSE);
1187:       }
1188:     }
1189:     nconv = k;
1190:     if (!ctx->lock && nep->reason == NEP_CONVERGED_ITERATING && !breakdown) { l += k; k = 0; }
1191:     if (l) PetscInfo(nep,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l);

1193:     /* Update S */
1194:     DSGetMat(nep->ds,ctx->nshifts?DS_MAT_Z:DS_MAT_Q,&MQ);
1195:     BVMultInPlace(ctx->V,MQ,nep->nconv,k+l);
1196:     MatDestroy(&MQ);

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

1201:     if (breakdown && nep->reason == NEP_CONVERGED_ITERATING) {
1202:       /* Stop if breakdown */
1203:       PetscInfo(nep,"Breakdown (it=%" PetscInt_FMT " norm=%g)\n",nep->its,(double)betah);
1204:       nep->reason = NEP_DIVERGED_BREAKDOWN;
1205:     }
1206:     if (nep->reason != NEP_CONVERGED_ITERATING) l--;
1207:     /* truncate S */
1208:     BVGetActiveColumns(nep->V,NULL,&nq);
1209:     if (k+l+deg<=nq) {
1210:       BVSetActiveColumns(ctx->V,nep->nconv,k+l+1);
1211:       if (!falselock && ctx->lock) BVTensorCompress(ctx->V,k-nep->nconv);
1212:       else BVTensorCompress(ctx->V,0);
1213:     }
1214:     nep->nconv = k;
1215:     if (!ctx->nshifts) {
1216:       for (i=0;i<nv;i++) { eigr[i] = nep->eigr[i]; eigi[i] = nep->eigi[i]; }
1217:       NEPNLEIGSBackTransform((PetscObject)nep,nv,eigr,eigi);
1218:     }
1219:     NEPMonitor(nep,nep->its,nconv,eigr,eigi,nep->errest,nv);
1220:   }
1221:   nep->nconv = nconv;
1222:   if (nep->nconv>0) {
1223:     BVSetActiveColumns(ctx->V,0,nep->nconv);
1224:     BVGetActiveColumns(nep->V,NULL,&nq);
1225:     BVSetActiveColumns(nep->V,0,nq);
1226:     if (nq>nep->nconv) {
1227:       BVTensorCompress(ctx->V,nep->nconv);
1228:       BVSetActiveColumns(nep->V,0,nep->nconv);
1229:       nq = nep->nconv;
1230:     }
1231:     if (ctx->nshifts) {
1232:       DSGetMat(nep->ds,DS_MAT_B,&MQ);
1233:       BVMultInPlace(ctx->V,MQ,0,nep->nconv);
1234:       MatDestroy(&MQ);
1235:     }
1236:     BVTensorGetFactors(ctx->V,NULL,&MS);
1237:     MatDenseGetArrayRead(MS,&S);
1238:     PetscMalloc1(nq*nep->nconv,&pU);
1239:     for (i=0;i<nep->nconv;i++) PetscArraycpy(pU+i*nq,S+i*lds,nq);
1240:     MatDenseRestoreArrayRead(MS,&S);
1241:     BVTensorRestoreFactors(ctx->V,NULL,&MS);
1242:     MatCreateSeqDense(PETSC_COMM_SELF,nq,nep->nconv,pU,&U);
1243:     BVSetActiveColumns(nep->V,0,nq);
1244:     BVMultInPlace(nep->V,U,0,nep->nconv);
1245:     BVSetActiveColumns(nep->V,0,nep->nconv);
1246:     MatDestroy(&U);
1247:     PetscFree(pU);
1248:     DSTruncate(nep->ds,nep->nconv,PETSC_TRUE);
1249:   }

1251:   /* Map eigenvalues back to the original problem */
1252:   if (!ctx->nshifts) {
1253:     NEPNLEIGSBackTransform((PetscObject)nep,nep->nconv,nep->eigr,nep->eigi);
1254:     PetscFree2(eigr,eigi);
1255:   }
1256:   BVDestroy(&W);
1257:   PetscFunctionReturn(0);
1258: }

1260: static PetscErrorCode NEPNLEIGSSetSingularitiesFunction_NLEIGS(NEP nep,PetscErrorCode (*fun)(NEP,PetscInt*,PetscScalar*,void*),void *ctx)
1261: {
1262:   NEP_NLEIGS *nepctx=(NEP_NLEIGS*)nep->data;

1264:   if (fun) nepctx->computesingularities = fun;
1265:   if (ctx) nepctx->singularitiesctx     = ctx;
1266:   nep->state = NEP_STATE_INITIAL;
1267:   PetscFunctionReturn(0);
1268: }

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

1274:    Logically Collective on nep

1276:    Input Parameters:
1277: +  nep - the NEP context
1278: .  fun - user function (if NULL then NEP retains any previously set value)
1279: -  ctx - [optional] user-defined context for private data for the function
1280:          (may be NULL, in which case NEP retains any previously set value)

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

1285: +   nep   - the NEP context
1286: .   maxnp - on input number of requested points in the discretization (can be set)
1287: .   xi    - computed values of the discretization
1288: -   ctx   - optional context, as set by NEPNLEIGSSetSingularitiesFunction()

1290:    Notes:
1291:    The user-defined function can set a smaller value of maxnp if necessary.
1292:    It is wrong to return a larger value.

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

1298:    Level: intermediate

1300: .seealso: NEPNLEIGSGetSingularitiesFunction(), NEPSetProblemType()
1301: @*/
1302: PetscErrorCode NEPNLEIGSSetSingularitiesFunction(NEP nep,PetscErrorCode (*fun)(NEP,PetscInt*,PetscScalar*,void*),void *ctx)
1303: {
1305:   PetscTryMethod(nep,"NEPNLEIGSSetSingularitiesFunction_C",(NEP,PetscErrorCode(*)(NEP,PetscInt*,PetscScalar*,void*),void*),(nep,fun,ctx));
1306:   PetscFunctionReturn(0);
1307: }

1309: static PetscErrorCode NEPNLEIGSGetSingularitiesFunction_NLEIGS(NEP nep,PetscErrorCode (**fun)(NEP,PetscInt*,PetscScalar*,void*),void **ctx)
1310: {
1311:   NEP_NLEIGS *nepctx=(NEP_NLEIGS*)nep->data;

1313:   if (fun) *fun = nepctx->computesingularities;
1314:   if (ctx) *ctx = nepctx->singularitiesctx;
1315:   PetscFunctionReturn(0);
1316: }

1318: /*@C
1319:    NEPNLEIGSGetSingularitiesFunction - Returns the Function and optionally the user
1320:    provided context for computing a discretization of the singularity set.

1322:    Not Collective

1324:    Input Parameter:
1325: .  nep - the nonlinear eigensolver context

1327:    Output Parameters:
1328: +  fun - location to put the function (or NULL)
1329: -  ctx - location to stash the function context (or NULL)

1331:    Level: advanced

1333: .seealso: NEPNLEIGSSetSingularitiesFunction()
1334: @*/
1335: PetscErrorCode NEPNLEIGSGetSingularitiesFunction(NEP nep,PetscErrorCode (**fun)(NEP,PetscInt*,PetscScalar*,void*),void **ctx)
1336: {
1338:   PetscUseMethod(nep,"NEPNLEIGSGetSingularitiesFunction_C",(NEP,PetscErrorCode(**)(NEP,PetscInt*,PetscScalar*,void*),void**),(nep,fun,ctx));
1339:   PetscFunctionReturn(0);
1340: }

1342: static PetscErrorCode NEPNLEIGSSetRestart_NLEIGS(NEP nep,PetscReal keep)
1343: {
1344:   NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;

1346:   if (keep==PETSC_DEFAULT) ctx->keep = 0.5;
1347:   else {
1349:     ctx->keep = keep;
1350:   }
1351:   PetscFunctionReturn(0);
1352: }

1354: /*@
1355:    NEPNLEIGSSetRestart - Sets the restart parameter for the NLEIGS
1356:    method, in particular the proportion of basis vectors that must be kept
1357:    after restart.

1359:    Logically Collective on nep

1361:    Input Parameters:
1362: +  nep  - the nonlinear eigensolver context
1363: -  keep - the number of vectors to be kept at restart

1365:    Options Database Key:
1366: .  -nep_nleigs_restart - Sets the restart parameter

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

1371:    Level: advanced

1373: .seealso: NEPNLEIGSGetRestart()
1374: @*/
1375: PetscErrorCode NEPNLEIGSSetRestart(NEP nep,PetscReal keep)
1376: {
1379:   PetscTryMethod(nep,"NEPNLEIGSSetRestart_C",(NEP,PetscReal),(nep,keep));
1380:   PetscFunctionReturn(0);
1381: }

1383: static PetscErrorCode NEPNLEIGSGetRestart_NLEIGS(NEP nep,PetscReal *keep)
1384: {
1385:   NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;

1387:   *keep = ctx->keep;
1388:   PetscFunctionReturn(0);
1389: }

1391: /*@
1392:    NEPNLEIGSGetRestart - Gets the restart parameter used in the NLEIGS method.

1394:    Not Collective

1396:    Input Parameter:
1397: .  nep - the nonlinear eigensolver context

1399:    Output Parameter:
1400: .  keep - the restart parameter

1402:    Level: advanced

1404: .seealso: NEPNLEIGSSetRestart()
1405: @*/
1406: PetscErrorCode NEPNLEIGSGetRestart(NEP nep,PetscReal *keep)
1407: {
1410:   PetscUseMethod(nep,"NEPNLEIGSGetRestart_C",(NEP,PetscReal*),(nep,keep));
1411:   PetscFunctionReturn(0);
1412: }

1414: static PetscErrorCode NEPNLEIGSSetLocking_NLEIGS(NEP nep,PetscBool lock)
1415: {
1416:   NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;

1418:   ctx->lock = lock;
1419:   PetscFunctionReturn(0);
1420: }

1422: /*@
1423:    NEPNLEIGSSetLocking - Choose between locking and non-locking variants of
1424:    the NLEIGS method.

1426:    Logically Collective on nep

1428:    Input Parameters:
1429: +  nep  - the nonlinear eigensolver context
1430: -  lock - true if the locking variant must be selected

1432:    Options Database Key:
1433: .  -nep_nleigs_locking - Sets the locking flag

1435:    Notes:
1436:    The default is to lock converged eigenpairs when the method restarts.
1437:    This behaviour can be changed so that all directions are kept in the
1438:    working subspace even if already converged to working accuracy (the
1439:    non-locking variant).

1441:    Level: advanced

1443: .seealso: NEPNLEIGSGetLocking()
1444: @*/
1445: PetscErrorCode NEPNLEIGSSetLocking(NEP nep,PetscBool lock)
1446: {
1449:   PetscTryMethod(nep,"NEPNLEIGSSetLocking_C",(NEP,PetscBool),(nep,lock));
1450:   PetscFunctionReturn(0);
1451: }

1453: static PetscErrorCode NEPNLEIGSGetLocking_NLEIGS(NEP nep,PetscBool *lock)
1454: {
1455:   NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;

1457:   *lock = ctx->lock;
1458:   PetscFunctionReturn(0);
1459: }

1461: /*@
1462:    NEPNLEIGSGetLocking - Gets the locking flag used in the NLEIGS method.

1464:    Not Collective

1466:    Input Parameter:
1467: .  nep - the nonlinear eigensolver context

1469:    Output Parameter:
1470: .  lock - the locking flag

1472:    Level: advanced

1474: .seealso: NEPNLEIGSSetLocking()
1475: @*/
1476: PetscErrorCode NEPNLEIGSGetLocking(NEP nep,PetscBool *lock)
1477: {
1480:   PetscUseMethod(nep,"NEPNLEIGSGetLocking_C",(NEP,PetscBool*),(nep,lock));
1481:   PetscFunctionReturn(0);
1482: }

1484: static PetscErrorCode NEPNLEIGSSetInterpolation_NLEIGS(NEP nep,PetscReal tol,PetscInt degree)
1485: {
1486:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;

1488:   if (tol == PETSC_DEFAULT) {
1489:     ctx->ddtol = PETSC_DEFAULT;
1490:     nep->state = NEP_STATE_INITIAL;
1491:   } else {
1493:     ctx->ddtol = tol;
1494:   }
1495:   if (degree == PETSC_DEFAULT || degree == PETSC_DECIDE) {
1496:     ctx->ddmaxit = 0;
1497:     if (nep->state) NEPReset(nep);
1498:     nep->state = NEP_STATE_INITIAL;
1499:   } else {
1501:     if (ctx->ddmaxit != degree) {
1502:       ctx->ddmaxit = degree;
1503:       if (nep->state) NEPReset(nep);
1504:       nep->state = NEP_STATE_INITIAL;
1505:     }
1506:   }
1507:   PetscFunctionReturn(0);
1508: }

1510: /*@
1511:    NEPNLEIGSSetInterpolation - Sets the tolerance and maximum degree
1512:    when building the interpolation via divided differences.

1514:    Logically Collective on nep

1516:    Input Parameters:
1517: +  nep    - the nonlinear eigensolver context
1518: .  tol    - tolerance to stop computing divided differences
1519: -  degree - maximum degree of interpolation

1521:    Options Database Key:
1522: +  -nep_nleigs_interpolation_tol <tol> - Sets the tolerance to stop computing divided differences
1523: -  -nep_nleigs_interpolation_degree <degree> - Sets the maximum degree of interpolation

1525:    Notes:
1526:    Use PETSC_DEFAULT for either argument to assign a reasonably good value.

1528:    Level: advanced

1530: .seealso: NEPNLEIGSGetInterpolation()
1531: @*/
1532: PetscErrorCode NEPNLEIGSSetInterpolation(NEP nep,PetscReal tol,PetscInt degree)
1533: {
1537:   PetscTryMethod(nep,"NEPNLEIGSSetInterpolation_C",(NEP,PetscReal,PetscInt),(nep,tol,degree));
1538:   PetscFunctionReturn(0);
1539: }

1541: static PetscErrorCode NEPNLEIGSGetInterpolation_NLEIGS(NEP nep,PetscReal *tol,PetscInt *degree)
1542: {
1543:   NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;

1545:   if (tol)    *tol    = ctx->ddtol;
1546:   if (degree) *degree = ctx->ddmaxit;
1547:   PetscFunctionReturn(0);
1548: }

1550: /*@
1551:    NEPNLEIGSGetInterpolation - Gets the tolerance and maximum degree
1552:    when building the interpolation via divided differences.

1554:    Not Collective

1556:    Input Parameter:
1557: .  nep - the nonlinear eigensolver context

1559:    Output Parameters:
1560: +  tol    - tolerance to stop computing divided differences
1561: -  degree - maximum degree of interpolation

1563:    Level: advanced

1565: .seealso: NEPNLEIGSSetInterpolation()
1566: @*/
1567: PetscErrorCode NEPNLEIGSGetInterpolation(NEP nep,PetscReal *tol,PetscInt *degree)
1568: {
1570:   PetscTryMethod(nep,"NEPNLEIGSGetInterpolation_C",(NEP,PetscReal*,PetscInt*),(nep,tol,degree));
1571:   PetscFunctionReturn(0);
1572: }

1574: static PetscErrorCode NEPNLEIGSSetRKShifts_NLEIGS(NEP nep,PetscInt ns,PetscScalar *shifts)
1575: {
1576:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
1577:   PetscInt       i;

1580:   if (ctx->nshifts) PetscFree(ctx->shifts);
1581:   for (i=0;i<ctx->nshiftsw;i++) KSPDestroy(&ctx->ksp[i]);
1582:   PetscFree(ctx->ksp);
1583:   ctx->ksp = NULL;
1584:   if (ns) {
1585:     PetscMalloc1(ns,&ctx->shifts);
1586:     for (i=0;i<ns;i++) ctx->shifts[i] = shifts[i];
1587:   }
1588:   ctx->nshifts = ns;
1589:   nep->state   = NEP_STATE_INITIAL;
1590:   PetscFunctionReturn(0);
1591: }

1593: /*@
1594:    NEPNLEIGSSetRKShifts - Sets a list of shifts to be used in the Rational
1595:    Krylov method.

1597:    Logically Collective on nep

1599:    Input Parameters:
1600: +  nep    - the nonlinear eigensolver context
1601: .  ns     - number of shifts
1602: -  shifts - array of scalar values specifying the shifts

1604:    Options Database Key:
1605: .  -nep_nleigs_rk_shifts - Sets the list of shifts

1607:    Notes:
1608:    If only one shift is provided, the built subspace built is equivalent to
1609:    shift-and-invert Krylov-Schur (provided that the absolute convergence
1610:    criterion is used).

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

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

1619:    Level: advanced

1621: .seealso: NEPNLEIGSGetRKShifts()
1622: @*/
1623: PetscErrorCode NEPNLEIGSSetRKShifts(NEP nep,PetscInt ns,PetscScalar shifts[])
1624: {
1628:   PetscTryMethod(nep,"NEPNLEIGSSetRKShifts_C",(NEP,PetscInt,PetscScalar*),(nep,ns,shifts));
1629:   PetscFunctionReturn(0);
1630: }

1632: static PetscErrorCode NEPNLEIGSGetRKShifts_NLEIGS(NEP nep,PetscInt *ns,PetscScalar **shifts)
1633: {
1634:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
1635:   PetscInt       i;

1637:   *ns = ctx->nshifts;
1638:   if (ctx->nshifts) {
1639:     PetscMalloc1(ctx->nshifts,shifts);
1640:     for (i=0;i<ctx->nshifts;i++) (*shifts)[i] = ctx->shifts[i];
1641:   }
1642:   PetscFunctionReturn(0);
1643: }

1645: /*@C
1646:    NEPNLEIGSGetRKShifts - Gets the list of shifts used in the Rational
1647:    Krylov method.

1649:    Not Collective

1651:    Input Parameter:
1652: .  nep - the nonlinear eigensolver context

1654:    Output Parameters:
1655: +  ns     - number of shifts
1656: -  shifts - array of shifts

1658:    Note:
1659:    The user is responsible for deallocating the returned array.

1661:    Level: advanced

1663: .seealso: NEPNLEIGSSetRKShifts()
1664: @*/
1665: PetscErrorCode NEPNLEIGSGetRKShifts(NEP nep,PetscInt *ns,PetscScalar *shifts[])
1666: {
1670:   PetscTryMethod(nep,"NEPNLEIGSGetRKShifts_C",(NEP,PetscInt*,PetscScalar**),(nep,ns,shifts));
1671:   PetscFunctionReturn(0);
1672: }

1674: static PetscErrorCode NEPNLEIGSGetKSPs_NLEIGS(NEP nep,PetscInt *nsolve,KSP **ksp)
1675: {
1676:   NEP_NLEIGS     *ctx = (NEP_NLEIGS*)nep->data;
1677:   PetscInt       i;
1678:   PC             pc;

1680:   if (!ctx->ksp) {
1681:     NEPNLEIGSSetShifts(nep,&ctx->nshiftsw);
1682:     PetscMalloc1(ctx->nshiftsw,&ctx->ksp);
1683:     for (i=0;i<ctx->nshiftsw;i++) {
1684:       KSPCreate(PetscObjectComm((PetscObject)nep),&ctx->ksp[i]);
1685:       PetscObjectIncrementTabLevel((PetscObject)ctx->ksp[i],(PetscObject)nep,1);
1686:       KSPSetOptionsPrefix(ctx->ksp[i],((PetscObject)nep)->prefix);
1687:       KSPAppendOptionsPrefix(ctx->ksp[i],"nep_nleigs_");
1688:       PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->ksp[i]);
1689:       PetscObjectSetOptions((PetscObject)ctx->ksp[i],((PetscObject)nep)->options);
1690:       KSPSetErrorIfNotConverged(ctx->ksp[i],PETSC_TRUE);
1691:       KSPSetTolerances(ctx->ksp[i],1e-3*SlepcDefaultTol(nep->tol),PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
1692:       KSPGetPC(ctx->ksp[i],&pc);
1693:       if ((nep->fui==NEP_USER_INTERFACE_SPLIT && nep->P) || (nep->fui==NEP_USER_INTERFACE_CALLBACK && nep->function_pre!=nep->function)) {
1694:         KSPSetType(ctx->ksp[i],KSPBCGS);
1695:         PCSetType(pc,PCBJACOBI);
1696:       } else {
1697:         KSPSetType(ctx->ksp[i],KSPPREONLY);
1698:         PCSetType(pc,PCLU);
1699:       }
1700:     }
1701:   }
1702:   if (nsolve) *nsolve = ctx->nshiftsw;
1703:   if (ksp)    *ksp    = ctx->ksp;
1704:   PetscFunctionReturn(0);
1705: }

1707: /*@C
1708:    NEPNLEIGSGetKSPs - Retrieve the array of linear solver objects associated with
1709:    the nonlinear eigenvalue solver.

1711:    Not Collective

1713:    Input Parameter:
1714: .  nep - nonlinear eigenvalue solver

1716:    Output Parameters:
1717: +  nsolve - number of returned KSP objects
1718: -  ksp - array of linear solver object

1720:    Notes:
1721:    The number of KSP objects is equal to the number of shifts provided by the user,
1722:    or 1 if the user did not provide shifts.

1724:    Level: advanced

1726: .seealso: NEPNLEIGSSetRKShifts()
1727: @*/
1728: PetscErrorCode NEPNLEIGSGetKSPs(NEP nep,PetscInt *nsolve,KSP **ksp)
1729: {
1731:   PetscUseMethod(nep,"NEPNLEIGSGetKSPs_C",(NEP,PetscInt*,KSP**),(nep,nsolve,ksp));
1732:   PetscFunctionReturn(0);
1733: }

1735: static PetscErrorCode NEPNLEIGSSetFullBasis_NLEIGS(NEP nep,PetscBool fullbasis)
1736: {
1737:   NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;

1739:   if (fullbasis!=ctx->fullbasis) {
1740:     ctx->fullbasis = fullbasis;
1741:     nep->state     = NEP_STATE_INITIAL;
1742:     nep->useds     = PetscNot(fullbasis);
1743:   }
1744:   PetscFunctionReturn(0);
1745: }

1747: /*@
1748:    NEPNLEIGSSetFullBasis - Choose between TOAR-basis (default) and full-basis
1749:    variants of the NLEIGS method.

1751:    Logically Collective on nep

1753:    Input Parameters:
1754: +  nep  - the nonlinear eigensolver context
1755: -  fullbasis - true if the full-basis variant must be selected

1757:    Options Database Key:
1758: .  -nep_nleigs_full_basis - Sets the full-basis flag

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

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

1770:    Level: advanced

1772: .seealso: NEPNLEIGSGetFullBasis(), NEPNLEIGSGetEPS(), NEPSetTwoSided(), BVCreateTensor()
1773: @*/
1774: PetscErrorCode NEPNLEIGSSetFullBasis(NEP nep,PetscBool fullbasis)
1775: {
1778:   PetscTryMethod(nep,"NEPNLEIGSSetFullBasis_C",(NEP,PetscBool),(nep,fullbasis));
1779:   PetscFunctionReturn(0);
1780: }

1782: static PetscErrorCode NEPNLEIGSGetFullBasis_NLEIGS(NEP nep,PetscBool *fullbasis)
1783: {
1784:   NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;

1786:   *fullbasis = ctx->fullbasis;
1787:   PetscFunctionReturn(0);
1788: }

1790: /*@
1791:    NEPNLEIGSGetFullBasis - Gets the flag that indicates if NLEIGS is using the
1792:    full-basis variant.

1794:    Not Collective

1796:    Input Parameter:
1797: .  nep - the nonlinear eigensolver context

1799:    Output Parameter:
1800: .  fullbasis - the flag

1802:    Level: advanced

1804: .seealso: NEPNLEIGSSetFullBasis()
1805: @*/
1806: PetscErrorCode NEPNLEIGSGetFullBasis(NEP nep,PetscBool *fullbasis)
1807: {
1810:   PetscUseMethod(nep,"NEPNLEIGSGetFullBasis_C",(NEP,PetscBool*),(nep,fullbasis));
1811:   PetscFunctionReturn(0);
1812: }

1814: #define SHIFTMAX 30

1816: PetscErrorCode NEPSetFromOptions_NLEIGS(PetscOptionItems *PetscOptionsObject,NEP nep)
1817: {
1818:   NEP_NLEIGS     *ctx = (NEP_NLEIGS*)nep->data;
1819:   PetscInt       i=0,k;
1820:   PetscBool      flg1,flg2,b;
1821:   PetscReal      r;
1822:   PetscScalar    array[SHIFTMAX];

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

1826:     PetscOptionsReal("-nep_nleigs_restart","Proportion of vectors kept after restart","NEPNLEIGSSetRestart",0.5,&r,&flg1);
1827:     if (flg1) NEPNLEIGSSetRestart(nep,r);

1829:     PetscOptionsBool("-nep_nleigs_locking","Choose between locking and non-locking variants","NEPNLEIGSSetLocking",PETSC_FALSE,&b,&flg1);
1830:     if (flg1) NEPNLEIGSSetLocking(nep,b);

1832:     PetscOptionsBool("-nep_nleigs_full_basis","Choose between TOAR and full-basis variants","NEPNLEIGSSetFullBasis",PETSC_FALSE,&b,&flg1);
1833:     if (flg1) NEPNLEIGSSetFullBasis(nep,b);

1835:     NEPNLEIGSGetInterpolation(nep,&r,&i);
1836:     if (!i) i = PETSC_DEFAULT;
1837:     PetscOptionsInt("-nep_nleigs_interpolation_degree","Maximum number of terms for interpolation via divided differences","NEPNLEIGSSetInterpolation",i,&i,&flg1);
1838:     PetscOptionsReal("-nep_nleigs_interpolation_tol","Tolerance for interpolation via divided differences","NEPNLEIGSSetInterpolation",r,&r,&flg2);
1839:     if (flg1 || flg2) NEPNLEIGSSetInterpolation(nep,r,i);

1841:     k = SHIFTMAX;
1842:     for (i=0;i<k;i++) array[i] = 0;
1843:     PetscOptionsScalarArray("-nep_nleigs_rk_shifts","Shifts for Rational Krylov","NEPNLEIGSSetRKShifts",array,&k,&flg1);
1844:     if (flg1) NEPNLEIGSSetRKShifts(nep,k,array);

1846:   PetscOptionsTail();

1848:   if (!ctx->ksp) NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp);
1849:   for (i=0;i<ctx->nshiftsw;i++) KSPSetFromOptions(ctx->ksp[i]);

1851:   if (ctx->fullbasis) {
1852:     if (!ctx->eps) NEPNLEIGSGetEPS(nep,&ctx->eps);
1853:     EPSSetFromOptions(ctx->eps);
1854:   }
1855:   PetscFunctionReturn(0);
1856: }

1858: PetscErrorCode NEPView_NLEIGS(NEP nep,PetscViewer viewer)
1859: {
1860:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
1861:   PetscBool      isascii;
1862:   PetscInt       i;
1863:   char           str[50];

1865:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
1866:   if (isascii) {
1867:     PetscViewerASCIIPrintf(viewer,"  %d%% of basis vectors kept after restart\n",(int)(100*ctx->keep));
1868:     if (ctx->fullbasis) PetscViewerASCIIPrintf(viewer,"  using the full-basis variant\n");
1869:     else PetscViewerASCIIPrintf(viewer,"  using the %slocking variant\n",ctx->lock?"":"non-");
1870:     PetscViewerASCIIPrintf(viewer,"  divided difference terms: used=%" PetscInt_FMT ", max=%" PetscInt_FMT "\n",ctx->nmat,ctx->ddmaxit);
1871:     PetscViewerASCIIPrintf(viewer,"  tolerance for divided difference convergence: %g\n",(double)ctx->ddtol);
1872:     if (ctx->nshifts) {
1873:       PetscViewerASCIIPrintf(viewer,"  RK shifts: ");
1874:       PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1875:       for (i=0;i<ctx->nshifts;i++) {
1876:         SlepcSNPrintfScalar(str,sizeof(str),ctx->shifts[i],PETSC_FALSE);
1877:         PetscViewerASCIIPrintf(viewer,"%s%s",str,(i<ctx->nshifts-1)?",":"");
1878:       }
1879:       PetscViewerASCIIPrintf(viewer,"\n");
1880:       PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1881:     }
1882:     if (!ctx->ksp) NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp);
1883:     PetscViewerASCIIPushTab(viewer);
1884:     KSPView(ctx->ksp[0],viewer);
1885:     PetscViewerASCIIPopTab(viewer);
1886:     if (ctx->fullbasis) {
1887:       if (!ctx->eps) NEPNLEIGSGetEPS(nep,&ctx->eps);
1888:       PetscViewerASCIIPushTab(viewer);
1889:       EPSView(ctx->eps,viewer);
1890:       PetscViewerASCIIPopTab(viewer);
1891:     }
1892:   }
1893:   PetscFunctionReturn(0);
1894: }

1896: PetscErrorCode NEPReset_NLEIGS(NEP nep)
1897: {
1898:   PetscInt       k;
1899:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;

1901:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) PetscFree(ctx->coeffD);
1902:   else {
1903:     for (k=0;k<ctx->nmat;k++) MatDestroy(&ctx->D[k]);
1904:   }
1905:   PetscFree4(ctx->s,ctx->xi,ctx->beta,ctx->D);
1906:   for (k=0;k<ctx->nshiftsw;k++) KSPReset(ctx->ksp[k]);
1907:   VecDestroy(&ctx->vrn);
1908:   if (ctx->fullbasis) {
1909:     MatDestroy(&ctx->A);
1910:     EPSReset(ctx->eps);
1911:     for (k=0;k<4;k++) VecDestroy(&ctx->w[k]);
1912:   }
1913:   PetscFunctionReturn(0);
1914: }

1916: PetscErrorCode NEPDestroy_NLEIGS(NEP nep)
1917: {
1918:   PetscInt       k;
1919:   NEP_NLEIGS     *ctx = (NEP_NLEIGS*)nep->data;

1921:   BVDestroy(&ctx->V);
1922:   for (k=0;k<ctx->nshiftsw;k++) KSPDestroy(&ctx->ksp[k]);
1923:   PetscFree(ctx->ksp);
1924:   if (ctx->nshifts) PetscFree(ctx->shifts);
1925:   if (ctx->fullbasis) EPSDestroy(&ctx->eps);
1926:   PetscFree(nep->data);
1927:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetSingularitiesFunction_C",NULL);
1928:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetSingularitiesFunction_C",NULL);
1929:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRestart_C",NULL);
1930:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRestart_C",NULL);
1931:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetLocking_C",NULL);
1932:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetLocking_C",NULL);
1933:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetInterpolation_C",NULL);
1934:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetInterpolation_C",NULL);
1935:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRKShifts_C",NULL);
1936:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRKShifts_C",NULL);
1937:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetKSPs_C",NULL);
1938:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetFullBasis_C",NULL);
1939:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetFullBasis_C",NULL);
1940:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetEPS_C",NULL);
1941:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetEPS_C",NULL);
1942:   PetscFunctionReturn(0);
1943: }

1945: SLEPC_EXTERN PetscErrorCode NEPCreate_NLEIGS(NEP nep)
1946: {
1947:   NEP_NLEIGS     *ctx;

1949:   PetscNewLog(nep,&ctx);
1950:   nep->data  = (void*)ctx;
1951:   ctx->lock  = PETSC_TRUE;
1952:   ctx->ddtol = PETSC_DEFAULT;

1954:   nep->useds = PETSC_TRUE;

1956:   nep->ops->setup          = NEPSetUp_NLEIGS;
1957:   nep->ops->setfromoptions = NEPSetFromOptions_NLEIGS;
1958:   nep->ops->view           = NEPView_NLEIGS;
1959:   nep->ops->destroy        = NEPDestroy_NLEIGS;
1960:   nep->ops->reset          = NEPReset_NLEIGS;

1962:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetSingularitiesFunction_C",NEPNLEIGSSetSingularitiesFunction_NLEIGS);
1963:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetSingularitiesFunction_C",NEPNLEIGSGetSingularitiesFunction_NLEIGS);
1964:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRestart_C",NEPNLEIGSSetRestart_NLEIGS);
1965:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRestart_C",NEPNLEIGSGetRestart_NLEIGS);
1966:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetLocking_C",NEPNLEIGSSetLocking_NLEIGS);
1967:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetLocking_C",NEPNLEIGSGetLocking_NLEIGS);
1968:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetInterpolation_C",NEPNLEIGSSetInterpolation_NLEIGS);
1969:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetInterpolation_C",NEPNLEIGSGetInterpolation_NLEIGS);
1970:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRKShifts_C",NEPNLEIGSSetRKShifts_NLEIGS);
1971:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRKShifts_C",NEPNLEIGSGetRKShifts_NLEIGS);
1972:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetKSPs_C",NEPNLEIGSGetKSPs_NLEIGS);
1973:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetFullBasis_C",NEPNLEIGSSetFullBasis_NLEIGS);
1974:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetFullBasis_C",NEPNLEIGSGetFullBasis_NLEIGS);
1975:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetEPS_C",NEPNLEIGSSetEPS_NLEIGS);
1976:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetEPS_C",NEPNLEIGSGetEPS_NLEIGS);
1977:   PetscFunctionReturn(0);
1978: }