Actual source code: nleigs.c

slepc-3.16.0 2021-09-30
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)
 65:   PetscReal      *rwork=NULL;
 66: #endif
 67:   PetscScalar    *work;
 68:   PetscBLASInt   info;

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

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

100: static PetscErrorCode NEPNLEIGSAuxiliarRmDuplicates(PetscInt nin,PetscScalar *pin,PetscInt *nout,PetscScalar *pout,PetscInt max)
101: {
102:   PetscInt i,j;

105:   for (i=0;i<nin;i++) {
106:     if (max && *nout>=max) break;
107:     pout[(*nout)++] = pin[i];
108:     for (j=0;j<*nout-1;j++)
109:       if (PetscAbsScalar(pin[i]-pout[j])<PETSC_MACHINE_EPSILON*100) {
110:         (*nout)--;
111:         break;
112:       }
113:   }
114:   return(0);
115: }

117: static PetscErrorCode NEPNLEIGSFNSingularities(FN f,PetscInt *nisol,PetscScalar **isol,PetscBool *rational)
118: {
120:   FNCombineType  ctype;
121:   FN             f1,f2;
122:   PetscInt       i,nq,nisol1,nisol2;
123:   PetscScalar    *qcoeff,*wr,*wi,*isol1,*isol2;
124:   PetscBool      flg,avail,rat1,rat2;

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

171: static PetscErrorCode NEPNLEIGSRationalSingularities(NEP nep,PetscInt *ndptx,PetscScalar *dxi,PetscBool *rational)
172: {
174:   PetscInt       nt,i,nisol;
175:   FN             f;
176:   PetscScalar    *isol;
177:   PetscBool      rat;

180:   *rational = PETSC_TRUE;
181:   *ndptx = 0;
182:   NEPGetSplitOperatorInfo(nep,&nt,NULL);
183:   for (i=0;i<nt;i++) {
184:     NEPGetSplitOperatorTerm(nep,i,NULL,&f);
185:     NEPNLEIGSFNSingularities(f,&nisol,&isol,&rat);
186:     if (nisol) {
187:       NEPNLEIGSAuxiliarRmDuplicates(nisol,isol,ndptx,dxi,0);
188:       PetscFree(isol);
189:     }
190:     *rational = ((*rational)&&rat)?PETSC_TRUE:PETSC_FALSE;
191:   }
192:   return(0);
193: }

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

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

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

264:   if (k==ndpt-1) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_CONV_FAILED,"Failed to determine singularities automatically in general problem");
265:   /* poles */
266:   PetscArrayzero(C,ndpt*ndpt);
267:   PetscArrayzero(A,ndpt*ndpt);
268:   for (i=0;i<=k;i++) {
269:     C[i+ndpt*i] = 1.0;
270:     A[(i+1)*ndpt] = ww[i];
271:     A[i+1] = 1.0;
272:     A[i+1+(i+1)*ndpt] = z[i];
273:   }
274:   C[0] = 0.0; C[k+1+(k+1)*ndpt] = 1.0;
275:   n_++;
276: #if defined(PETSC_USE_COMPLEX)
277:   PetscStackCallBLAS("LAPACKggev",LAPACKggev_("N","N",&n_,A,&lda_,C,&lda_,D,N,NULL,&lda_,NULL,&lda_,work,&lwork,rwork,&info));
278: #else
279:   PetscStackCallBLAS("LAPACKggev",LAPACKggev_("N","N",&n_,A,&lda_,C,&lda_,D,VT,N,NULL,&lda_,NULL,&lda_,work,&lwork,&info));
280: #endif
281:   SlepcCheckLapackInfo("ggev",info);
282:   cont = 0.0;
283:   for (i=0;i<n_;i++) if (N[i]!=0.0) {
284:     dxi[cont++] = D[i]/N[i];
285:   }
286:   *ndptx = cont;
287:   PetscFPTrapPop();
288:   PetscFree5(R,z,f,C,ww);
289:   PetscFree6(A,S,VT,work,D,N);
290: #if defined(PETSC_USE_COMPLEX)
291:   PetscFree(rwork);
292: #endif
293:   return(0);
294: }

296: /*  Singularities using Adaptive Anderson-Antoulas algorithm */
297: static PetscErrorCode NEPNLEIGSAAASingularities(NEP nep,PetscInt ndpt,PetscScalar *ds,PetscInt *ndptx,PetscScalar *dxi)
298: {
300:   Vec            u,v,w;
301:   PetscRandom    rand;
302:   PetscScalar    *F,*isol;
303:   PetscInt       i,k,nisol,nt;
304:   Mat            T;
305:   FN             f;

308:   PetscMalloc1(ndpt,&F);
309:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
310:     PetscMalloc1(ndpt,&isol);
311:     *ndptx = 0;
312:     NEPGetSplitOperatorInfo(nep,&nt,NULL);
313:     nisol = *ndptx;
314:     for (k=0;k<nt;k++) {
315:       NEPGetSplitOperatorTerm(nep,k,NULL,&f);
316:       for (i=0;i<ndpt;i++) {
317:         FNEvaluateFunction(f,ds[i],&F[i]);
318:       }
319:       NEPNLEIGSAAAComputation(nep,ndpt,ds,F,&nisol,isol);
320:       if (nisol) {
321:         NEPNLEIGSAuxiliarRmDuplicates(nisol,isol,ndptx,dxi,ndpt);
322:       }
323:     }
324:     PetscFree(isol);
325:   } else {
326:     MatCreateVecs(nep->function,&u,NULL);
327:     VecDuplicate(u,&v);
328:     VecDuplicate(u,&w);
329:     BVGetRandomContext(nep->V,&rand);
330:     VecSetRandom(u,rand);
331:     VecNormalize(u,NULL);
332:     VecSetRandom(v,rand);
333:     VecNormalize(v,NULL);
334:     T = nep->function;
335:     for (i=0;i<ndpt;i++) {
336:       NEPComputeFunction(nep,ds[i],T,T);
337:       MatMult(T,v,w);
338:       VecDot(w,u,&F[i]);
339:     }
340:     NEPNLEIGSAAAComputation(nep,ndpt,ds,F,ndptx,dxi);
341:     VecDestroy(&u);
342:     VecDestroy(&v);
343:     VecDestroy(&w);
344:   }
345:   PetscFree(F);
346:   return(0);
347: }

349: static PetscErrorCode NEPNLEIGSLejaBagbyPoints(NEP nep)
350: {
352:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
353:   PetscInt       i,k,ndpt=NDPOINTS,ndptx=NDPOINTS;
354:   PetscScalar    *ds,*dsi,*dxi,*nrs,*nrxi,*s=ctx->s,*xi=ctx->xi,*beta=ctx->beta;
355:   PetscReal      maxnrs,minnrxi;
356:   PetscBool      rational;
357: #if !defined(PETSC_USE_COMPLEX)
358:   PetscReal      a,b,h;
359: #endif

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

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

420: PetscErrorCode NEPNLEIGSEvalNRTFunct(NEP nep,PetscInt k,PetscScalar sigma,PetscScalar *b)
421: {
422:   NEP_NLEIGS  *ctx=(NEP_NLEIGS*)nep->data;
423:   PetscInt    i;
424:   PetscScalar *beta=ctx->beta,*s=ctx->s,*xi=ctx->xi;

427:   b[0] = 1.0/beta[0];
428:   for (i=0;i<k;i++) {
429:     b[i+1] = ((sigma-s[i])*b[i])/(beta[i+1]*(1.0-sigma/xi[i]));
430:   }
431:   return(0);
432: }

434: static PetscErrorCode MatMult_Fun(Mat A,Vec x,Vec y)
435: {
436:   PetscErrorCode      ierr;
437:   NEP_NLEIGS_MATSHELL *ctx;
438:   PetscInt            i;

441:   MatShellGetContext(A,&ctx);
442:   MatMult(ctx->A[0],x,y);
443:   if (ctx->coeff[0]!=1.0) { VecScale(y,ctx->coeff[0]); }
444:   for (i=1;i<ctx->nmat;i++) {
445:     MatMult(ctx->A[i],x,ctx->t);
446:     VecAXPY(y,ctx->coeff[i],ctx->t);
447:   }
448:   return(0);
449: }

451: static PetscErrorCode MatMultTranspose_Fun(Mat A,Vec x,Vec y)
452: {
453:   PetscErrorCode      ierr;
454:   NEP_NLEIGS_MATSHELL *ctx;
455:   PetscInt            i;

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

468: static PetscErrorCode MatGetDiagonal_Fun(Mat A,Vec diag)
469: {
470:   PetscErrorCode      ierr;
471:   NEP_NLEIGS_MATSHELL *ctx;
472:   PetscInt            i;

475:   MatShellGetContext(A,&ctx);
476:   MatGetDiagonal(ctx->A[0],diag);
477:   if (ctx->coeff[0]!=1.0) { VecScale(diag,ctx->coeff[0]); }
478:   for (i=1;i<ctx->nmat;i++) {
479:     MatGetDiagonal(ctx->A[i],ctx->t);
480:     VecAXPY(diag,ctx->coeff[i],ctx->t);
481:   }
482:   return(0);
483: }

485: static PetscErrorCode MatDuplicate_Fun(Mat A,MatDuplicateOption op,Mat *B)
486: {
487:   PetscInt            n,i;
488:   NEP_NLEIGS_MATSHELL *ctxnew,*ctx;
489:   void                (*fun)(void);
490:   PetscErrorCode      ierr;

493:   MatShellGetContext(A,&ctx);
494:   PetscNew(&ctxnew);
495:   ctxnew->nmat = ctx->nmat;
496:   ctxnew->maxnmat = ctx->maxnmat;
497:   PetscMalloc2(ctxnew->maxnmat,&ctxnew->A,ctxnew->maxnmat,&ctxnew->coeff);
498:   for (i=0;i<ctx->nmat;i++) {
499:     PetscObjectReference((PetscObject)ctx->A[i]);
500:     ctxnew->A[i] = ctx->A[i];
501:     ctxnew->coeff[i] = ctx->coeff[i];
502:   }
503:   MatGetSize(ctx->A[0],&n,NULL);
504:   VecDuplicate(ctx->t,&ctxnew->t);
505:   MatCreateShell(PetscObjectComm((PetscObject)A),n,n,n,n,(void*)ctxnew,B);
506:   MatShellSetManageScalingShifts(*B);
507:   MatShellGetOperation(A,MATOP_MULT,&fun);
508:   MatShellSetOperation(*B,MATOP_MULT,fun);
509:   MatShellGetOperation(A,MATOP_MULT_TRANSPOSE,&fun);
510:   MatShellSetOperation(*B,MATOP_MULT_TRANSPOSE,fun);
511:   MatShellGetOperation(A,MATOP_GET_DIAGONAL,&fun);
512:   MatShellSetOperation(*B,MATOP_GET_DIAGONAL,fun);
513:   MatShellGetOperation(A,MATOP_DUPLICATE,&fun);
514:   MatShellSetOperation(*B,MATOP_DUPLICATE,fun);
515:   MatShellGetOperation(A,MATOP_DESTROY,&fun);
516:   MatShellSetOperation(*B,MATOP_DESTROY,fun);
517:   MatShellGetOperation(A,MATOP_AXPY,&fun);
518:   MatShellSetOperation(*B,MATOP_AXPY,fun);
519:   return(0);
520: }

522: static PetscErrorCode MatDestroy_Fun(Mat A)
523: {
524:   NEP_NLEIGS_MATSHELL *ctx;
525:   PetscErrorCode      ierr;
526:   PetscInt            i;

529:   if (A) {
530:     MatShellGetContext(A,&ctx);
531:     for (i=0;i<ctx->nmat;i++) {
532:       MatDestroy(&ctx->A[i]);
533:     }
534:     VecDestroy(&ctx->t);
535:     PetscFree2(ctx->A,ctx->coeff);
536:     PetscFree(ctx);
537:   }
538:   return(0);
539: }

541: static PetscErrorCode MatAXPY_Fun(Mat Y,PetscScalar a,Mat X,MatStructure str)
542: {
543:   NEP_NLEIGS_MATSHELL *ctxY,*ctxX;
544:   PetscErrorCode      ierr;
545:   PetscInt            i,j;
546:   PetscBool           found;

549:   MatShellGetContext(Y,&ctxY);
550:   MatShellGetContext(X,&ctxX);
551:   for (i=0;i<ctxX->nmat;i++) {
552:     found = PETSC_FALSE;
553:     for (j=0;!found&&j<ctxY->nmat;j++) {
554:       if (ctxX->A[i]==ctxY->A[j]) {
555:         found = PETSC_TRUE;
556:         ctxY->coeff[j] += a*ctxX->coeff[i];
557:       }
558:     }
559:     if (!found) {
560:       ctxY->coeff[ctxY->nmat] = a*ctxX->coeff[i];
561:       ctxY->A[ctxY->nmat++] = ctxX->A[i];
562:       PetscObjectReference((PetscObject)ctxX->A[i]);
563:     }
564:   }
565:   return(0);
566: }

568: static PetscErrorCode MatScale_Fun(Mat M,PetscScalar a)
569: {
570:   NEP_NLEIGS_MATSHELL *ctx;
571:   PetscErrorCode      ierr;
572:   PetscInt            i;

575:   MatShellGetContext(M,&ctx);
576:   for (i=0;i<ctx->nmat;i++) ctx->coeff[i] *= a;
577:   return(0);
578: }

580: static PetscErrorCode NLEIGSMatToMatShellArray(Mat M,Mat *Ms,PetscInt maxnmat)
581: {
582:   PetscErrorCode      ierr;
583:   NEP_NLEIGS_MATSHELL *ctx;
584:   PetscInt            n;
585:   PetscBool           has;

588:   MatHasOperation(M,MATOP_DUPLICATE,&has);
589:   if (!has) SETERRQ(PetscObjectComm((PetscObject)M),PETSC_ERR_USER,"MatDuplicate operation required");
590:   PetscNew(&ctx);
591:   ctx->maxnmat = maxnmat;
592:   PetscMalloc2(ctx->maxnmat,&ctx->A,ctx->maxnmat,&ctx->coeff);
593:   MatDuplicate(M,MAT_COPY_VALUES,&ctx->A[0]);
594:   ctx->nmat = 1;
595:   ctx->coeff[0] = 1.0;
596:   MatCreateVecs(M,&ctx->t,NULL);
597:   MatGetSize(M,&n,NULL);
598:   MatCreateShell(PetscObjectComm((PetscObject)M),n,n,n,n,(void*)ctx,Ms);
599:   MatShellSetManageScalingShifts(*Ms);
600:   MatShellSetOperation(*Ms,MATOP_MULT,(void(*)(void))MatMult_Fun);
601:   MatShellSetOperation(*Ms,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_Fun);
602:   MatShellSetOperation(*Ms,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Fun);
603:   MatShellSetOperation(*Ms,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_Fun);
604:   MatShellSetOperation(*Ms,MATOP_DESTROY,(void(*)(void))MatDestroy_Fun);
605:   MatShellSetOperation(*Ms,MATOP_AXPY,(void(*)(void))MatAXPY_Fun);
606:   MatShellSetOperation(*Ms,MATOP_SCALE,(void(*)(void))MatScale_Fun);
607:   return(0);
608: }

610: static PetscErrorCode NEPNLEIGSDividedDifferences_split(NEP nep)
611: {
613:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
614:   PetscInt       k,j,i,maxnmat,nmax;
615:   PetscReal      norm0,norm,*matnorm;
616:   PetscScalar    *s=ctx->s,*beta=ctx->beta,*xi=ctx->xi,*b,alpha,*coeffs,*pK,*pH,sone=1.0;
617:   Mat            T,Ts,K,H;
618:   PetscBool      shell,hasmnorm=PETSC_FALSE,matrix=PETSC_TRUE;
619:   PetscBLASInt   n_;

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

728: static PetscErrorCode NEPNLEIGSDividedDifferences_callback(NEP nep)
729: {
731:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
732:   PetscInt       k,j,i,maxnmat;
733:   PetscReal      norm0,norm;
734:   PetscScalar    *s=ctx->s,*beta=ctx->beta,*b,*coeffs;
735:   Mat            *D=ctx->D,T;
736:   PetscBool      shell,has,vec=PETSC_FALSE;
737:   PetscRandom    rand;
738:   Vec            w[2];

741:   PetscMalloc2(ctx->ddmaxit+1,&b,ctx->ddmaxit+1,&coeffs);
742:   BVGetRandomContext(nep->V,&rand);
743:   T = nep->function;
744:   NEPComputeFunction(nep,s[0],T,T);
745:   MatIsShell(T,&shell);
746:   maxnmat = PetscMax(ctx->ddmaxit,nep->nt);
747:   if (!shell) {
748:     MatDuplicate(T,MAT_COPY_VALUES,&D[0]);
749:   } else {
750:     NLEIGSMatToMatShellArray(T,&D[0],maxnmat);
751:   }
752:   if (beta[0]!=1.0) {
753:     MatScale(D[0],1.0/beta[0]);
754:   }
755:   MatHasOperation(D[0],MATOP_NORM,&has);
756:   if (has) {
757:     MatNorm(D[0],NORM_FROBENIUS,&norm0);
758:   } else {
759:     MatCreateVecs(D[0],NULL,&w[0]);
760:     VecDuplicate(w[0],&w[1]);
761:     VecDuplicate(w[0],&ctx->vrn);
762:     VecSetRandomNormal(ctx->vrn,rand,w[0],w[1]);
763:     VecNormalize(ctx->vrn,NULL);
764:     vec = PETSC_TRUE;
765:     MatNormEstimate(D[0],ctx->vrn,w[0],&norm0);
766:   }
767:   ctx->nmat = ctx->ddmaxit;
768:   for (k=1;k<ctx->ddmaxit;k++) {
769:     NEPNLEIGSEvalNRTFunct(nep,k,s[k],b);
770:     NEPComputeFunction(nep,s[k],T,T);
771:     if (!shell) {
772:       MatDuplicate(T,MAT_COPY_VALUES,&D[k]);
773:     } else {
774:       NLEIGSMatToMatShellArray(T,&D[k],maxnmat);
775:     }
776:     for (j=0;j<k;j++) {
777:       MatAXPY(D[k],-b[j],D[j],nep->mstr);
778:     }
779:     MatScale(D[k],1.0/b[k]);
780:     MatHasOperation(D[k],MATOP_NORM,&has);
781:     if (has) {
782:       MatNorm(D[k],NORM_FROBENIUS,&norm);
783:     } else {
784:       if (!vec) {
785:         MatCreateVecs(D[k],NULL,&w[0]);
786:         VecDuplicate(w[0],&w[1]);
787:         VecDuplicate(w[0],&ctx->vrn);
788:         VecSetRandomNormal(ctx->vrn,rand,w[0],w[1]);
789:         VecNormalize(ctx->vrn,NULL);
790:         vec = PETSC_TRUE;
791:       }
792:       MatNormEstimate(D[k],ctx->vrn,w[0],&norm);
793:     }
794:     if (k>1 && norm/norm0 < ctx->ddtol && k>1) {
795:       ctx->nmat = k+1;
796:       break;
797:     }
798:   }
799:   if (!ctx->ksp) { NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp); }
800:   for (i=0;i<ctx->nshiftsw;i++) {
801:     NEPNLEIGSEvalNRTFunct(nep,ctx->nmat-1,ctx->shifts[i],coeffs);
802:     MatDuplicate(ctx->D[0],MAT_COPY_VALUES,&T);
803:     if (coeffs[0]!=1.0) { MatScale(T,coeffs[0]); }
804:     for (j=1;j<ctx->nmat;j++) {
805:       MatAXPY(T,coeffs[j],ctx->D[j],nep->mstr);
806:     }
807:     KSPSetOperators(ctx->ksp[i],T,T);
808:     KSPSetUp(ctx->ksp[i]);
809:     MatDestroy(&T);
810:   }
811:   PetscFree2(b,coeffs);
812:   if (vec) {
813:     VecDestroy(&w[0]);
814:     VecDestroy(&w[1]);
815:   }
816:   return(0);
817: }

819: /*
820:    NEPKrylovConvergence - This is the analogue to EPSKrylovConvergence.
821: */
822: PetscErrorCode NEPNLEIGSKrylovConvergence(NEP nep,PetscBool getall,PetscInt kini,PetscInt nits,PetscReal betah,PetscScalar betak,PetscInt *kout,Vec *w)
823: {
825:   PetscInt       k,newk,marker,inside;
826:   PetscScalar    re,im;
827:   PetscReal      resnorm,tt;
828:   PetscBool      istrivial;
829:   NEP_NLEIGS     *ctx = (NEP_NLEIGS*)nep->data;

832:   RGIsTrivial(nep->rg,&istrivial);
833:   marker = -1;
834:   if (nep->trackall) getall = PETSC_TRUE;
835:   for (k=kini;k<kini+nits;k++) {
836:     /* eigenvalue */
837:     re = nep->eigr[k];
838:     im = nep->eigi[k];
839:     if (!istrivial) {
840:       if (!ctx->nshifts) {
841:         NEPNLEIGSBackTransform((PetscObject)nep,1,&re,&im);
842:       }
843:       RGCheckInside(nep->rg,1,&re,&im,&inside);
844:       if (marker==-1 && inside<0) marker = k;
845:     }
846:     newk = k;
847:     DSVectors(nep->ds,DS_MAT_X,&newk,&resnorm);
848:     tt = (ctx->nshifts)?SlepcAbsEigenvalue(betak-nep->eigr[k]*betah,nep->eigi[k]*betah):betah;
849:     resnorm *=  PetscAbsReal(tt);
850:     /* error estimate */
851:     (*nep->converged)(nep,nep->eigr[k],nep->eigi[k],resnorm,&nep->errest[k],nep->convergedctx);
852:     if (marker==-1 && nep->errest[k] >= nep->tol) marker = k;
853:     if (newk==k+1) {
854:       nep->errest[k+1] = nep->errest[k];
855:       k++;
856:     }
857:     if (marker!=-1 && !getall) break;
858:   }
859:   if (marker!=-1) k = marker;
860:   *kout = k;
861:   return(0);
862: }

864: PetscErrorCode NEPSetUp_NLEIGS(NEP nep)
865: {
867:   PetscInt       k,in;
868:   PetscScalar    zero=0.0;
869:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
870:   SlepcSC        sc;
871:   PetscBool      istrivial;

874:   NEPSetDimensions_Default(nep,nep->nev,&nep->ncv,&nep->mpd);
875:   if (nep->ncv>nep->nev+nep->mpd) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_USER_INPUT,"The value of ncv must not be larger than nev+mpd");
876:   if (nep->max_it==PETSC_DEFAULT) nep->max_it = PetscMax(5000,2*nep->n/nep->ncv);
877:   if (!ctx->ddmaxit) ctx->ddmaxit = LBPOINTS;
878:   RGIsTrivial(nep->rg,&istrivial);
879:   if (istrivial) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"NEPNLEIGS requires a nontrivial region defining the target set");
880:   if (!nep->which) nep->which = NEP_TARGET_MAGNITUDE;
881:   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");

883:   /* Initialize the NLEIGS context structure */
884:   k = ctx->ddmaxit;
885:   PetscMalloc4(k,&ctx->s,k,&ctx->xi,k,&ctx->beta,k,&ctx->D);
886:   nep->data = ctx;
887:   if (nep->tol==PETSC_DEFAULT) nep->tol = SLEPC_DEFAULT_TOL;
888:   if (ctx->ddtol==PETSC_DEFAULT) ctx->ddtol = nep->tol/10.0;
889:   if (!ctx->keep) ctx->keep = 0.5;

891:   /* Compute Leja-Bagby points and scaling values */
892:   NEPNLEIGSLejaBagbyPoints(nep);
893:   if (nep->problem_type!=NEP_RATIONAL) {
894:     RGCheckInside(nep->rg,1,&nep->target,&zero,&in);
895:     if (in<0) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"The target is not inside the target set");
896:   }

898:   /* Compute the divided difference matrices */
899:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
900:     NEPNLEIGSDividedDifferences_split(nep);
901:   } else {
902:     NEPNLEIGSDividedDifferences_callback(nep);
903:   }
904:   NEPAllocateSolution(nep,ctx->nmat-1);
905:   NEPSetWorkVecs(nep,4);
906:   if (!ctx->fullbasis) {
907:     if (nep->twosided) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Two-sided variant requires the full-basis option, rerun with -nep_nleigs_full_basis");
908:     /* set-up DS and transfer split operator functions */
909:     DSSetType(nep->ds,ctx->nshifts?DSGNHEP:DSNHEP);
910:     DSAllocate(nep->ds,nep->ncv+1);
911:     DSGetSlepcSC(nep->ds,&sc);
912:     if (!ctx->nshifts) sc->map = NEPNLEIGSBackTransform;
913:     DSSetExtraRow(nep->ds,PETSC_TRUE);
914:     sc->mapobj        = (PetscObject)nep;
915:     sc->rg            = nep->rg;
916:     sc->comparison    = nep->sc->comparison;
917:     sc->comparisonctx = nep->sc->comparisonctx;
918:     BVDestroy(&ctx->V);
919:     BVCreateTensor(nep->V,ctx->nmat-1,&ctx->V);
920:     nep->ops->solve          = NEPSolve_NLEIGS;
921:     nep->ops->computevectors = NEPComputeVectors_Schur;
922:   } else {
923:     NEPSetUp_NLEIGS_FullBasis(nep);
924:     nep->ops->solve          = NEPSolve_NLEIGS_FullBasis;
925:     nep->ops->computevectors = NULL;
926:   }
927:   return(0);
928: }

930: /*
931:   Extend the TOAR basis by applying the the matrix operator
932:   over a vector which is decomposed on the TOAR way
933:   Input:
934:     - S,V: define the latest Arnoldi vector (nv vectors in V)
935:   Output:
936:     - t: new vector extending the TOAR basis
937:     - r: temporally coefficients to compute the TOAR coefficients
938:          for the new Arnoldi vector
939:   Workspace: t_ (two vectors)
940: */
941: 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_)
942: {
944:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
945:   PetscInt       deg=ctx->nmat-1,k,j;
946:   Vec            v=t_[0],q=t_[1],w;
947:   PetscScalar    *beta=ctx->beta,*s=ctx->s,*xi=ctx->xi,*coeffs,sigma;

950:   if (!ctx->ksp) { NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp); }
951:   sigma = ctx->shifts[idxrktg];
952:   BVSetActiveColumns(nep->V,0,nv);
953:   PetscMalloc1(ctx->nmat,&coeffs);
954:   if (PetscAbsScalar(s[deg-2]-sigma)<100*PETSC_MACHINE_EPSILON) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"Breakdown in NLEIGS");
955:   /* i-part stored in (i-1) position */
956:   for (j=0;j<nv;j++) {
957:     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);
958:   }
959:   BVSetActiveColumns(W,0,deg);
960:   BVGetColumn(W,deg-1,&w);
961:   BVMultVec(V,1.0/beta[deg],0,w,S+(deg-1)*ls);
962:   BVRestoreColumn(W,deg-1,&w);
963:   BVGetColumn(W,deg-2,&w);
964:   BVMultVec(V,1.0,0.0,w,r+(deg-2)*lr);
965:   BVRestoreColumn(W,deg-2,&w);
966:   for (k=deg-2;k>0;k--) {
967:     if (PetscAbsScalar(s[k-1]-sigma)<100*PETSC_MACHINE_EPSILON) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"Breakdown in NLEIGS");
968:     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);
969:     BVGetColumn(W,k-1,&w);
970:     BVMultVec(V,1.0,0.0,w,r+(k-1)*lr);
971:     BVRestoreColumn(W,k-1,&w);
972:   }
973:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
974:     for (j=0;j<ctx->nmat-2;j++) coeffs[j] = ctx->coeffD[nep->nt*j];
975:     coeffs[ctx->nmat-2] = ctx->coeffD[nep->nt*(ctx->nmat-1)];
976:     BVMultVec(W,1.0,0.0,v,coeffs);
977:     MatMult(nep->A[0],v,q);
978:     for (k=1;k<nep->nt;k++) {
979:       for (j=0;j<ctx->nmat-2;j++) coeffs[j] = ctx->coeffD[nep->nt*j+k];
980:       coeffs[ctx->nmat-2] = ctx->coeffD[nep->nt*(ctx->nmat-1)+k];
981:       BVMultVec(W,1.0,0,v,coeffs);
982:       MatMult(nep->A[k],v,t);
983:       VecAXPY(q,1.0,t);
984:     }
985:     KSPSolve(ctx->ksp[idxrktg],q,t);
986:     VecScale(t,-1.0);
987:   } else {
988:     for (k=0;k<deg-1;k++) {
989:       BVGetColumn(W,k,&w);
990:       MatMult(ctx->D[k],w,q);
991:       BVRestoreColumn(W,k,&w);
992:       BVInsertVec(W,k,q);
993:     }
994:     BVGetColumn(W,deg-1,&w);
995:     MatMult(ctx->D[deg],w,q);
996:     BVRestoreColumn(W,k,&w);
997:     BVInsertVec(W,k,q);
998:     for (j=0;j<ctx->nmat-1;j++) coeffs[j] = 1.0;
999:     BVMultVec(W,1.0,0.0,q,coeffs);
1000:     KSPSolve(ctx->ksp[idxrktg],q,t);
1001:     VecScale(t,-1.0);
1002:   }
1003:   PetscFree(coeffs);
1004:   return(0);
1005: }

1007: /*
1008:   Compute TOAR coefficients of the blocks of the new Arnoldi vector computed
1009: */
1010: static PetscErrorCode NEPTOARCoefficients(NEP nep,PetscScalar sigma,PetscInt nv,PetscScalar *S,PetscInt ls,PetscScalar *r,PetscInt lr,PetscScalar *x,PetscScalar *work)
1011: {
1013:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
1014:   PetscInt       k,j,d=ctx->nmat-1;
1015:   PetscScalar    *t=work;

1018:   NEPNLEIGSEvalNRTFunct(nep,d-1,sigma,t);
1019:   for (k=0;k<d-1;k++) {
1020:     for (j=0;j<=nv;j++) r[k*lr+j] += t[k]*x[j];
1021:   }
1022:   for (j=0;j<=nv;j++) r[(d-1)*lr+j] = t[d-1]*x[j];
1023:   return(0);
1024: }

1026: /*
1027:   Compute continuation vector coefficients for the Rational-Krylov run.
1028:   dim(work) >= (end-ini)*(end-ini+1) + end+1 + 2*(end-ini+1), dim(t) = end.
1029: */
1030: 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)
1031: {
1033:   PetscScalar    *x,*W,*tau,sone=1.0,szero=0.0;
1034:   PetscInt       i,j,n1,n,nwu=0;
1035:   PetscBLASInt   info,n_,n1_,one=1,dim,lds_;
1036:   NEP_NLEIGS     *ctx = (NEP_NLEIGS*)nep->data;

1039:   if (!ctx->nshifts || !end) {
1040:     t[0] = 1;
1041:     PetscArraycpy(cont,S+end*lds,lds);
1042:   } else {
1043:     n   = end-ini;
1044:     n1  = n+1;
1045:     x   = work+nwu;
1046:     nwu += end+1;
1047:     tau = work+nwu;
1048:     nwu += n;
1049:     W   = work+nwu;
1050:     nwu += n1*n;
1051:     for (j=ini;j<end;j++) {
1052:       for (i=ini;i<=end;i++) W[(j-ini)*n1+i-ini] = K[j*ld+i] -H[j*ld+i]*sigma;
1053:     }
1054:     PetscBLASIntCast(n,&n_);
1055:     PetscBLASIntCast(n1,&n1_);
1056:     PetscBLASIntCast(end+1,&dim);
1057:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1058:     PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&n1_,&n_,W,&n1_,tau,work+nwu,&n1_,&info));
1059:     SlepcCheckLapackInfo("geqrf",info);
1060:     for (i=0;i<end;i++) t[i] = 0.0;
1061:     t[end] = 1.0;
1062:     for (j=n-1;j>=0;j--) {
1063:       for (i=0;i<ini+j;i++) x[i] = 0.0;
1064:       x[ini+j] = 1.0;
1065:       for (i=j+1;i<n1;i++) x[i+ini] = W[i+n1*j];
1066:       tau[j] = PetscConj(tau[j]);
1067:       PetscStackCallBLAS("LAPACKlarf",LAPACKlarf_("L",&dim,&one,x,&one,tau+j,t,&dim,work+nwu));
1068:     }
1069:     PetscBLASIntCast(lds,&lds_);
1070:     PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&lds_,&n1_,&sone,S,&lds_,t,&one,&szero,cont,&one));
1071:     PetscFPTrapPop();
1072:   }
1073:   return(0);
1074: }

1076: /*
1077:   Compute a run of Arnoldi iterations
1078: */
1079: PetscErrorCode NEPNLEIGSTOARrun(NEP nep,PetscScalar *K,PetscScalar *H,PetscInt ldh,BV W,PetscInt k,PetscInt *M,PetscBool *breakdown,Vec *t_)
1080: {
1082:   NEP_NLEIGS     *ctx = (NEP_NLEIGS*)nep->data;
1083:   PetscInt       i,j,m=*M,lwa,deg=ctx->nmat-1,lds,nqt,ld,l;
1084:   Vec            t;
1085:   PetscReal      norm;
1086:   PetscScalar    *x,*work,*tt,sigma,*cont,*S;
1087:   PetscBool      lindep;
1088:   Mat            MS;

1091:   BVTensorGetFactors(ctx->V,NULL,&MS);
1092:   MatDenseGetArray(MS,&S);
1093:   BVGetSizes(nep->V,NULL,NULL,&ld);
1094:   lds = ld*deg;
1095:   BVGetActiveColumns(nep->V,&l,&nqt);
1096:   lwa = PetscMax(ld,deg)+(m+1)*(m+1)+4*(m+1);
1097:   PetscMalloc4(ld,&x,lwa,&work,m+1,&tt,lds,&cont);
1098:   BVSetActiveColumns(ctx->V,0,m);
1099:   for (j=k;j<m;j++) {
1100:     sigma = ctx->shifts[(++(ctx->idxrk))%ctx->nshiftsw];

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

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

1110:     /* orthogonalize */
1111:     BVOrthogonalizeColumn(nep->V,nqt,x,&norm,&lindep);
1112:     if (!lindep) {
1113:       x[nqt] = norm;
1114:       BVScaleColumn(nep->V,nqt,1.0/norm);
1115:       nqt++;
1116:     } else x[nqt] = 0.0;

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

1120:     /* Level-2 orthogonalization */
1121:     BVOrthogonalizeColumn(ctx->V,j+1,H+j*ldh,&norm,breakdown);
1122:     H[j+1+ldh*j] = norm;
1123:     if (ctx->nshifts) {
1124:       for (i=0;i<=j;i++) K[i+ldh*j] = sigma*H[i+ldh*j] + tt[i];
1125:       K[j+1+ldh*j] = sigma*H[j+1+ldh*j];
1126:     }
1127:     if (*breakdown) {
1128:       *M = j+1;
1129:       break;
1130:     }
1131:     BVScaleColumn(ctx->V,j+1,1.0/norm);
1132:     BVSetActiveColumns(nep->V,l,nqt);
1133:   }
1134:   PetscFree4(x,work,tt,cont);
1135:   MatDenseRestoreArray(MS,&S);
1136:   BVTensorRestoreFactors(ctx->V,NULL,&MS);
1137:   return(0);
1138: }

1140: PetscErrorCode NEPSolve_NLEIGS(NEP nep)
1141: {
1142:   PetscErrorCode    ierr;
1143:   NEP_NLEIGS        *ctx = (NEP_NLEIGS*)nep->data;
1144:   PetscInt          i,k=0,l,nv=0,ld,lds,ldds,nq;
1145:   PetscInt          deg=ctx->nmat-1,nconv=0,dsn,dsk;
1146:   PetscScalar       *H,*pU,*K,betak=0,*eigr,*eigi;
1147:   const PetscScalar *S;
1148:   PetscReal         betah;
1149:   PetscBool         falselock=PETSC_FALSE,breakdown=PETSC_FALSE;
1150:   BV                W;
1151:   Mat               MS,MQ,U;

1154:   if (ctx->lock) {
1155:     /* undocumented option to use a cheaper locking instead of the true locking */
1156:     PetscOptionsGetBool(NULL,NULL,"-nep_nleigs_falselocking",&falselock,NULL);
1157:   }

1159:   BVGetSizes(nep->V,NULL,NULL,&ld);
1160:   lds = deg*ld;
1161:   DSGetLeadingDimension(nep->ds,&ldds);
1162:   if (!ctx->nshifts) {
1163:     PetscMalloc2(nep->ncv,&eigr,nep->ncv,&eigi);
1164:   } else { eigr = nep->eigr; eigi = nep->eigi; }
1165:   BVDuplicateResize(nep->V,PetscMax(nep->nt-1,ctx->nmat-1),&W);

1167:   /* clean projected matrix (including the extra-arrow) */
1168:   DSGetArray(nep->ds,DS_MAT_A,&H);
1169:   PetscArrayzero(H,ldds*ldds);
1170:   DSRestoreArray(nep->ds,DS_MAT_A,&H);
1171:   if (ctx->nshifts) {
1172:     DSGetArray(nep->ds,DS_MAT_B,&H);
1173:     PetscArrayzero(H,ldds*ldds);
1174:     DSRestoreArray(nep->ds,DS_MAT_B,&H);
1175:   }

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

1180:   /* Restart loop */
1181:   l = 0;
1182:   while (nep->reason == NEP_CONVERGED_ITERATING) {
1183:     nep->its++;

1185:     /* Compute an nv-step Krylov relation */
1186:     nv = PetscMin(nep->nconv+nep->mpd,nep->ncv);
1187:     if (ctx->nshifts) { DSGetArray(nep->ds,DS_MAT_A,&K); }
1188:     DSGetArray(nep->ds,ctx->nshifts?DS_MAT_B:DS_MAT_A,&H);
1189:     NEPNLEIGSTOARrun(nep,K,H,ldds,W,nep->nconv+l,&nv,&breakdown,nep->work);
1190:     betah = PetscAbsScalar(H[(nv-1)*ldds+nv]);
1191:     DSRestoreArray(nep->ds,ctx->nshifts?DS_MAT_B:DS_MAT_A,&H);
1192:     if (ctx->nshifts) {
1193:       betak = K[(nv-1)*ldds+nv];
1194:       DSRestoreArray(nep->ds,DS_MAT_A,&K);
1195:     }
1196:     DSSetDimensions(nep->ds,nv,nep->nconv,nep->nconv+l);
1197:     if (l==0) {
1198:       DSSetState(nep->ds,DS_STATE_INTERMEDIATE);
1199:     } else {
1200:       DSSetState(nep->ds,DS_STATE_RAW);
1201:     }

1203:     /* Solve projected problem */
1204:     DSSolve(nep->ds,nep->eigr,nep->eigi);
1205:     DSSort(nep->ds,nep->eigr,nep->eigi,NULL,NULL,NULL);
1206:     DSUpdateExtraRow(nep->ds);
1207:     DSSynchronize(nep->ds,nep->eigr,nep->eigi);

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

1213:     /* Update l */
1214:     if (nep->reason != NEP_CONVERGED_ITERATING || breakdown) l = 0;
1215:     else {
1216:       l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
1217:       DSGetTruncateSize(nep->ds,k,nv,&l);
1218:       if (!breakdown) {
1219:         /* Prepare the Rayleigh quotient for restart */
1220:         DSGetDimensions(nep->ds,&dsn,NULL,&dsk,NULL);
1221:         DSSetDimensions(nep->ds,dsn,k,dsk);
1222:         DSTruncate(nep->ds,k+l,PETSC_FALSE);
1223:       }
1224:     }
1225:     nconv = k;
1226:     if (!ctx->lock && nep->reason == NEP_CONVERGED_ITERATING && !breakdown) { l += k; k = 0; }
1227:     if (l) { PetscInfo1(nep,"Preparing to restart keeping l=%D vectors\n",l); }

1229:     /* Update S */
1230:     DSGetMat(nep->ds,ctx->nshifts?DS_MAT_Z:DS_MAT_Q,&MQ);
1231:     BVMultInPlace(ctx->V,MQ,nep->nconv,k+l);
1232:     MatDestroy(&MQ);

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

1237:     if (breakdown && nep->reason == NEP_CONVERGED_ITERATING) {
1238:       /* Stop if breakdown */
1239:       PetscInfo2(nep,"Breakdown (it=%D norm=%g)\n",nep->its,(double)betah);
1240:       nep->reason = NEP_DIVERGED_BREAKDOWN;
1241:     }
1242:     if (nep->reason != NEP_CONVERGED_ITERATING) l--;
1243:     /* truncate S */
1244:     BVGetActiveColumns(nep->V,NULL,&nq);
1245:     if (k+l+deg<=nq) {
1246:       BVSetActiveColumns(ctx->V,nep->nconv,k+l+1);
1247:       if (!falselock && ctx->lock) {
1248:         BVTensorCompress(ctx->V,k-nep->nconv);
1249:       } else {
1250:         BVTensorCompress(ctx->V,0);
1251:       }
1252:     }
1253:     nep->nconv = k;
1254:     if (!ctx->nshifts) {
1255:       for (i=0;i<nv;i++) { eigr[i] = nep->eigr[i]; eigi[i] = nep->eigi[i]; }
1256:       NEPNLEIGSBackTransform((PetscObject)nep,nv,eigr,eigi);
1257:     }
1258:     NEPMonitor(nep,nep->its,nconv,eigr,eigi,nep->errest,nv);
1259:   }
1260:   nep->nconv = nconv;
1261:   if (nep->nconv>0) {
1262:     BVSetActiveColumns(ctx->V,0,nep->nconv);
1263:     BVGetActiveColumns(nep->V,NULL,&nq);
1264:     BVSetActiveColumns(nep->V,0,nq);
1265:     if (nq>nep->nconv) {
1266:       BVTensorCompress(ctx->V,nep->nconv);
1267:       BVSetActiveColumns(nep->V,0,nep->nconv);
1268:       nq = nep->nconv;
1269:     }
1270:     if (ctx->nshifts) {
1271:       DSGetMat(nep->ds,DS_MAT_B,&MQ);
1272:       BVMultInPlace(ctx->V,MQ,0,nep->nconv);
1273:       MatDestroy(&MQ);
1274:     }
1275:     BVTensorGetFactors(ctx->V,NULL,&MS);
1276:     MatDenseGetArrayRead(MS,&S);
1277:     PetscMalloc1(nq*nep->nconv,&pU);
1278:     for (i=0;i<nep->nconv;i++) {
1279:       PetscArraycpy(pU+i*nq,S+i*lds,nq);
1280:     }
1281:     MatDenseRestoreArrayRead(MS,&S);
1282:     BVTensorRestoreFactors(ctx->V,NULL,&MS);
1283:     MatCreateSeqDense(PETSC_COMM_SELF,nq,nep->nconv,pU,&U);
1284:     BVSetActiveColumns(nep->V,0,nq);
1285:     BVMultInPlace(nep->V,U,0,nep->nconv);
1286:     BVSetActiveColumns(nep->V,0,nep->nconv);
1287:     MatDestroy(&U);
1288:     PetscFree(pU);
1289:     DSTruncate(nep->ds,nep->nconv,PETSC_TRUE);
1290:   }

1292:   /* Map eigenvalues back to the original problem */
1293:   if (!ctx->nshifts) {
1294:     NEPNLEIGSBackTransform((PetscObject)nep,nep->nconv,nep->eigr,nep->eigi);
1295:     PetscFree2(eigr,eigi);
1296:   }
1297:   BVDestroy(&W);
1298:   return(0);
1299: }

1301: static PetscErrorCode NEPNLEIGSSetSingularitiesFunction_NLEIGS(NEP nep,PetscErrorCode (*fun)(NEP,PetscInt*,PetscScalar*,void*),void *ctx)
1302: {
1303:   NEP_NLEIGS *nepctx=(NEP_NLEIGS*)nep->data;

1306:   if (fun) nepctx->computesingularities = fun;
1307:   if (ctx) nepctx->singularitiesctx     = ctx;
1308:   nep->state = NEP_STATE_INITIAL;
1309:   return(0);
1310: }

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

1316:    Logically Collective on nep

1318:    Input Parameters:
1319: +  nep - the NEP context
1320: .  fun - user function (if NULL then NEP retains any previously set value)
1321: -  ctx - [optional] user-defined context for private data for the function
1322:          (may be NULL, in which case NEP retains any previously set value)

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

1327: +   nep   - the NEP context
1328: .   maxnp - on input number of requested points in the discretization (can be set)
1329: .   xi    - computed values of the discretization
1330: -   ctx   - optional context, as set by NEPNLEIGSSetSingularitiesFunction()

1332:    Notes:
1333:    The user-defined function can set a smaller value of maxnp if necessary.
1334:    It is wrong to return a larger value.

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

1340:    Level: intermediate

1342: .seealso: NEPNLEIGSGetSingularitiesFunction(), NEPSetProblemType()
1343: @*/
1344: PetscErrorCode NEPNLEIGSSetSingularitiesFunction(NEP nep,PetscErrorCode (*fun)(NEP,PetscInt*,PetscScalar*,void*),void *ctx)
1345: {

1350:   PetscTryMethod(nep,"NEPNLEIGSSetSingularitiesFunction_C",(NEP,PetscErrorCode(*)(NEP,PetscInt*,PetscScalar*,void*),void*),(nep,fun,ctx));
1351:   return(0);
1352: }

1354: static PetscErrorCode NEPNLEIGSGetSingularitiesFunction_NLEIGS(NEP nep,PetscErrorCode (**fun)(NEP,PetscInt*,PetscScalar*,void*),void **ctx)
1355: {
1356:   NEP_NLEIGS *nepctx=(NEP_NLEIGS*)nep->data;

1359:   if (fun) *fun = nepctx->computesingularities;
1360:   if (ctx) *ctx = nepctx->singularitiesctx;
1361:   return(0);
1362: }

1364: /*@C
1365:    NEPNLEIGSGetSingularitiesFunction - Returns the Function and optionally the user
1366:    provided context for computing a discretization of the singularity set.

1368:    Not Collective

1370:    Input Parameter:
1371: .  nep - the nonlinear eigensolver context

1373:    Output Parameters:
1374: +  fun - location to put the function (or NULL)
1375: -  ctx - location to stash the function context (or NULL)

1377:    Level: advanced

1379: .seealso: NEPNLEIGSSetSingularitiesFunction()
1380: @*/
1381: PetscErrorCode NEPNLEIGSGetSingularitiesFunction(NEP nep,PetscErrorCode (**fun)(NEP,PetscInt*,PetscScalar*,void*),void **ctx)
1382: {

1387:   PetscUseMethod(nep,"NEPNLEIGSGetSingularitiesFunction_C",(NEP,PetscErrorCode(**)(NEP,PetscInt*,PetscScalar*,void*),void**),(nep,fun,ctx));
1388:   return(0);
1389: }

1391: static PetscErrorCode NEPNLEIGSSetRestart_NLEIGS(NEP nep,PetscReal keep)
1392: {
1393:   NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;

1396:   if (keep==PETSC_DEFAULT) ctx->keep = 0.5;
1397:   else {
1398:     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]");
1399:     ctx->keep = keep;
1400:   }
1401:   return(0);
1402: }

1404: /*@
1405:    NEPNLEIGSSetRestart - Sets the restart parameter for the NLEIGS
1406:    method, in particular the proportion of basis vectors that must be kept
1407:    after restart.

1409:    Logically Collective on nep

1411:    Input Parameters:
1412: +  nep  - the nonlinear eigensolver context
1413: -  keep - the number of vectors to be kept at restart

1415:    Options Database Key:
1416: .  -nep_nleigs_restart - Sets the restart parameter

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

1421:    Level: advanced

1423: .seealso: NEPNLEIGSGetRestart()
1424: @*/
1425: PetscErrorCode NEPNLEIGSSetRestart(NEP nep,PetscReal keep)
1426: {

1432:   PetscTryMethod(nep,"NEPNLEIGSSetRestart_C",(NEP,PetscReal),(nep,keep));
1433:   return(0);
1434: }

1436: static PetscErrorCode NEPNLEIGSGetRestart_NLEIGS(NEP nep,PetscReal *keep)
1437: {
1438:   NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;

1441:   *keep = ctx->keep;
1442:   return(0);
1443: }

1445: /*@
1446:    NEPNLEIGSGetRestart - Gets the restart parameter used in the NLEIGS method.

1448:    Not Collective

1450:    Input Parameter:
1451: .  nep - the nonlinear eigensolver context

1453:    Output Parameter:
1454: .  keep - the restart parameter

1456:    Level: advanced

1458: .seealso: NEPNLEIGSSetRestart()
1459: @*/
1460: PetscErrorCode NEPNLEIGSGetRestart(NEP nep,PetscReal *keep)
1461: {

1467:   PetscUseMethod(nep,"NEPNLEIGSGetRestart_C",(NEP,PetscReal*),(nep,keep));
1468:   return(0);
1469: }

1471: static PetscErrorCode NEPNLEIGSSetLocking_NLEIGS(NEP nep,PetscBool lock)
1472: {
1473:   NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;

1476:   ctx->lock = lock;
1477:   return(0);
1478: }

1480: /*@
1481:    NEPNLEIGSSetLocking - Choose between locking and non-locking variants of
1482:    the NLEIGS method.

1484:    Logically Collective on nep

1486:    Input Parameters:
1487: +  nep  - the nonlinear eigensolver context
1488: -  lock - true if the locking variant must be selected

1490:    Options Database Key:
1491: .  -nep_nleigs_locking - Sets the locking flag

1493:    Notes:
1494:    The default is to lock converged eigenpairs when the method restarts.
1495:    This behaviour can be changed so that all directions are kept in the
1496:    working subspace even if already converged to working accuracy (the
1497:    non-locking variant).

1499:    Level: advanced

1501: .seealso: NEPNLEIGSGetLocking()
1502: @*/
1503: PetscErrorCode NEPNLEIGSSetLocking(NEP nep,PetscBool lock)
1504: {

1510:   PetscTryMethod(nep,"NEPNLEIGSSetLocking_C",(NEP,PetscBool),(nep,lock));
1511:   return(0);
1512: }

1514: static PetscErrorCode NEPNLEIGSGetLocking_NLEIGS(NEP nep,PetscBool *lock)
1515: {
1516:   NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;

1519:   *lock = ctx->lock;
1520:   return(0);
1521: }

1523: /*@
1524:    NEPNLEIGSGetLocking - Gets the locking flag used in the NLEIGS method.

1526:    Not Collective

1528:    Input Parameter:
1529: .  nep - the nonlinear eigensolver context

1531:    Output Parameter:
1532: .  lock - the locking flag

1534:    Level: advanced

1536: .seealso: NEPNLEIGSSetLocking()
1537: @*/
1538: PetscErrorCode NEPNLEIGSGetLocking(NEP nep,PetscBool *lock)
1539: {

1545:   PetscUseMethod(nep,"NEPNLEIGSGetLocking_C",(NEP,PetscBool*),(nep,lock));
1546:   return(0);
1547: }

1549: static PetscErrorCode NEPNLEIGSSetInterpolation_NLEIGS(NEP nep,PetscReal tol,PetscInt degree)
1550: {
1552:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;

1555:   if (tol == PETSC_DEFAULT) {
1556:     ctx->ddtol = PETSC_DEFAULT;
1557:     nep->state = NEP_STATE_INITIAL;
1558:   } else {
1559:     if (tol <= 0.0) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of tol. Must be > 0");
1560:     ctx->ddtol = tol;
1561:   }
1562:   if (degree == PETSC_DEFAULT || degree == PETSC_DECIDE) {
1563:     ctx->ddmaxit = 0;
1564:     if (nep->state) { NEPReset(nep); }
1565:     nep->state = NEP_STATE_INITIAL;
1566:   } else {
1567:     if (degree <= 0) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of degree. Must be > 0");
1568:     if (ctx->ddmaxit != degree) {
1569:       ctx->ddmaxit = degree;
1570:       if (nep->state) { NEPReset(nep); }
1571:       nep->state = NEP_STATE_INITIAL;
1572:     }
1573:   }
1574:   return(0);
1575: }

1577: /*@
1578:    NEPNLEIGSSetInterpolation - Sets the tolerance and maximum degree
1579:    when building the interpolation via divided differences.

1581:    Logically Collective on nep

1583:    Input Parameters:
1584: +  nep    - the nonlinear eigensolver context
1585: .  tol    - tolerance to stop computing divided differences
1586: -  degree - maximum degree of interpolation

1588:    Options Database Key:
1589: +  -nep_nleigs_interpolation_tol <tol> - Sets the tolerance to stop computing divided differences
1590: -  -nep_nleigs_interpolation_degree <degree> - Sets the maximum degree of interpolation

1592:    Notes:
1593:    Use PETSC_DEFAULT for either argument to assign a reasonably good value.

1595:    Level: advanced

1597: .seealso: NEPNLEIGSGetInterpolation()
1598: @*/
1599: PetscErrorCode NEPNLEIGSSetInterpolation(NEP nep,PetscReal tol,PetscInt degree)
1600: {

1607:   PetscTryMethod(nep,"NEPNLEIGSSetInterpolation_C",(NEP,PetscReal,PetscInt),(nep,tol,degree));
1608:   return(0);
1609: }

1611: static PetscErrorCode NEPNLEIGSGetInterpolation_NLEIGS(NEP nep,PetscReal *tol,PetscInt *degree)
1612: {
1613:   NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;

1616:   if (tol)    *tol    = ctx->ddtol;
1617:   if (degree) *degree = ctx->ddmaxit;
1618:   return(0);
1619: }

1621: /*@
1622:    NEPNLEIGSGetInterpolation - Gets the tolerance and maximum degree
1623:    when building the interpolation via divided differences.

1625:    Not Collective

1627:    Input Parameter:
1628: .  nep - the nonlinear eigensolver context

1630:    Output Parameters:
1631: +  tol    - tolerance to stop computing divided differences
1632: -  degree - maximum degree of interpolation

1634:    Level: advanced

1636: .seealso: NEPNLEIGSSetInterpolation()
1637: @*/
1638: PetscErrorCode NEPNLEIGSGetInterpolation(NEP nep,PetscReal *tol,PetscInt *degree)
1639: {

1644:   PetscTryMethod(nep,"NEPNLEIGSGetInterpolation_C",(NEP,PetscReal*,PetscInt*),(nep,tol,degree));
1645:   return(0);
1646: }

1648: static PetscErrorCode NEPNLEIGSSetRKShifts_NLEIGS(NEP nep,PetscInt ns,PetscScalar *shifts)
1649: {
1651:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
1652:   PetscInt       i;

1655:   if (ns<0) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"Number of shifts must be non-negative");
1656:   if (ctx->nshifts) { PetscFree(ctx->shifts); }
1657:   for (i=0;i<ctx->nshiftsw;i++) { KSPDestroy(&ctx->ksp[i]); }
1658:   PetscFree(ctx->ksp);
1659:   ctx->ksp = NULL;
1660:   if (ns) {
1661:     PetscMalloc1(ns,&ctx->shifts);
1662:     for (i=0;i<ns;i++) ctx->shifts[i] = shifts[i];
1663:   }
1664:   ctx->nshifts = ns;
1665:   nep->state   = NEP_STATE_INITIAL;
1666:   return(0);
1667: }

1669: /*@
1670:    NEPNLEIGSSetRKShifts - Sets a list of shifts to be used in the Rational
1671:    Krylov method.

1673:    Logically Collective on nep

1675:    Input Parameters:
1676: +  nep    - the nonlinear eigensolver context
1677: .  ns     - number of shifts
1678: -  shifts - array of scalar values specifying the shifts

1680:    Options Database Key:
1681: .  -nep_nleigs_rk_shifts - Sets the list of shifts

1683:    Notes:
1684:    If only one shift is provided, the built subspace built is equivalent to
1685:    shift-and-invert Krylov-Schur (provided that the absolute convergence
1686:    criterion is used).

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

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

1695:    Level: advanced

1697: .seealso: NEPNLEIGSGetRKShifts()
1698: @*/
1699: PetscErrorCode NEPNLEIGSSetRKShifts(NEP nep,PetscInt ns,PetscScalar shifts[])
1700: {

1707:   PetscTryMethod(nep,"NEPNLEIGSSetRKShifts_C",(NEP,PetscInt,PetscScalar*),(nep,ns,shifts));
1708:   return(0);
1709: }

1711: static PetscErrorCode NEPNLEIGSGetRKShifts_NLEIGS(NEP nep,PetscInt *ns,PetscScalar **shifts)
1712: {
1714:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
1715:   PetscInt       i;

1718:   *ns = ctx->nshifts;
1719:   if (ctx->nshifts) {
1720:     PetscMalloc1(ctx->nshifts,shifts);
1721:     for (i=0;i<ctx->nshifts;i++) (*shifts)[i] = ctx->shifts[i];
1722:   }
1723:   return(0);
1724: }

1726: /*@C
1727:    NEPNLEIGSGetRKShifts - Gets the list of shifts used in the Rational
1728:    Krylov method.

1730:    Not Collective

1732:    Input Parameter:
1733: .  nep - the nonlinear eigensolver context

1735:    Output Parameters:
1736: +  ns     - number of shifts
1737: -  shifts - array of shifts

1739:    Note:
1740:    The user is responsible for deallocating the returned array.

1742:    Level: advanced

1744: .seealso: NEPNLEIGSSetRKShifts()
1745: @*/
1746: PetscErrorCode NEPNLEIGSGetRKShifts(NEP nep,PetscInt *ns,PetscScalar *shifts[])
1747: {

1754:   PetscTryMethod(nep,"NEPNLEIGSGetRKShifts_C",(NEP,PetscInt*,PetscScalar**),(nep,ns,shifts));
1755:   return(0);
1756: }

1758: static PetscErrorCode NEPNLEIGSGetKSPs_NLEIGS(NEP nep,PetscInt *nsolve,KSP **ksp)
1759: {
1761:   NEP_NLEIGS     *ctx = (NEP_NLEIGS*)nep->data;
1762:   PetscInt       i;
1763:   PC             pc;

1766:   if (!ctx->ksp) {
1767:     NEPNLEIGSSetShifts(nep,&ctx->nshiftsw);
1768:     PetscMalloc1(ctx->nshiftsw,&ctx->ksp);
1769:     for (i=0;i<ctx->nshiftsw;i++) {
1770:       KSPCreate(PetscObjectComm((PetscObject)nep),&ctx->ksp[i]);
1771:       PetscObjectIncrementTabLevel((PetscObject)ctx->ksp[i],(PetscObject)nep,1);
1772:       KSPSetOptionsPrefix(ctx->ksp[i],((PetscObject)nep)->prefix);
1773:       KSPAppendOptionsPrefix(ctx->ksp[i],"nep_nleigs_");
1774:       PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->ksp[i]);
1775:       PetscObjectSetOptions((PetscObject)ctx->ksp[i],((PetscObject)nep)->options);
1776:       KSPSetErrorIfNotConverged(ctx->ksp[i],PETSC_TRUE);
1777:       KSPSetTolerances(ctx->ksp[i],SlepcDefaultTol(nep->tol),PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
1778:       KSPGetPC(ctx->ksp[i],&pc);
1779:       KSPSetType(ctx->ksp[i],KSPPREONLY);
1780:       PCSetType(pc,PCLU);
1781:     }
1782:   }
1783:   if (nsolve) *nsolve = ctx->nshiftsw;
1784:   if (ksp)    *ksp    = ctx->ksp;
1785:   return(0);
1786: }

1788: /*@C
1789:    NEPNLEIGSGetKSPs - Retrieve the array of linear solver objects associated with
1790:    the nonlinear eigenvalue solver.

1792:    Not Collective

1794:    Input Parameter:
1795: .  nep - nonlinear eigenvalue solver

1797:    Output Parameters:
1798: +  nsolve - number of returned KSP objects
1799: -  ksp - array of linear solver object

1801:    Notes:
1802:    The number of KSP objects is equal to the number of shifts provided by the user,
1803:    or 1 if the user did not provide shifts.

1805:    Level: advanced

1807: .seealso: NEPNLEIGSSetRKShifts()
1808: @*/
1809: PetscErrorCode NEPNLEIGSGetKSPs(NEP nep,PetscInt *nsolve,KSP **ksp)
1810: {

1815:   PetscUseMethod(nep,"NEPNLEIGSGetKSPs_C",(NEP,PetscInt*,KSP**),(nep,nsolve,ksp));
1816:   return(0);
1817: }

1819: static PetscErrorCode NEPNLEIGSSetFullBasis_NLEIGS(NEP nep,PetscBool fullbasis)
1820: {
1821:   NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;

1824:   if (fullbasis!=ctx->fullbasis) {
1825:     ctx->fullbasis = fullbasis;
1826:     nep->state     = NEP_STATE_INITIAL;
1827:     nep->useds     = PetscNot(fullbasis);
1828:   }
1829:   return(0);
1830: }

1832: /*@
1833:    NEPNLEIGSSetFullBasis - Choose between TOAR-basis (default) and full-basis
1834:    variants of the NLEIGS method.

1836:    Logically Collective on nep

1838:    Input Parameters:
1839: +  nep  - the nonlinear eigensolver context
1840: -  fullbasis - true if the full-basis variant must be selected

1842:    Options Database Key:
1843: .  -nep_nleigs_full_basis - Sets the full-basis flag

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

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

1855:    Level: advanced

1857: .seealso: NEPNLEIGSGetFullBasis(), NEPNLEIGSGetEPS(), NEPSetTwoSided(), BVCreateTensor()
1858: @*/
1859: PetscErrorCode NEPNLEIGSSetFullBasis(NEP nep,PetscBool fullbasis)
1860: {

1866:   PetscTryMethod(nep,"NEPNLEIGSSetFullBasis_C",(NEP,PetscBool),(nep,fullbasis));
1867:   return(0);
1868: }

1870: static PetscErrorCode NEPNLEIGSGetFullBasis_NLEIGS(NEP nep,PetscBool *fullbasis)
1871: {
1872:   NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;

1875:   *fullbasis = ctx->fullbasis;
1876:   return(0);
1877: }

1879: /*@
1880:    NEPNLEIGSGetFullBasis - Gets the flag that indicates if NLEIGS is using the
1881:    full-basis variant.

1883:    Not Collective

1885:    Input Parameter:
1886: .  nep - the nonlinear eigensolver context

1888:    Output Parameter:
1889: .  fullbasis - the flag

1891:    Level: advanced

1893: .seealso: NEPNLEIGSSetFullBasis()
1894: @*/
1895: PetscErrorCode NEPNLEIGSGetFullBasis(NEP nep,PetscBool *fullbasis)
1896: {

1902:   PetscUseMethod(nep,"NEPNLEIGSGetFullBasis_C",(NEP,PetscBool*),(nep,fullbasis));
1903:   return(0);
1904: }

1906: #define SHIFTMAX 30

1908: PetscErrorCode NEPSetFromOptions_NLEIGS(PetscOptionItems *PetscOptionsObject,NEP nep)
1909: {
1911:   NEP_NLEIGS     *ctx = (NEP_NLEIGS*)nep->data;
1912:   PetscInt       i=0,k;
1913:   PetscBool      flg1,flg2,b;
1914:   PetscReal      r;
1915:   PetscScalar    array[SHIFTMAX];

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

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

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

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

1929:     NEPNLEIGSGetInterpolation(nep,&r,&i);
1930:     if (!i) i = PETSC_DEFAULT;
1931:     PetscOptionsInt("-nep_nleigs_interpolation_degree","Maximum number of terms for interpolation via divided differences","NEPNLEIGSSetInterpolation",i,&i,&flg1);
1932:     PetscOptionsReal("-nep_nleigs_interpolation_tol","Tolerance for interpolation via divided differences","NEPNLEIGSSetInterpolation",r,&r,&flg2);
1933:     if (flg1 || flg2) { NEPNLEIGSSetInterpolation(nep,r,i); }

1935:     k = SHIFTMAX;
1936:     for (i=0;i<k;i++) array[i] = 0;
1937:     PetscOptionsScalarArray("-nep_nleigs_rk_shifts","Shifts for Rational Krylov","NEPNLEIGSSetRKShifts",array,&k,&flg1);
1938:     if (flg1) { NEPNLEIGSSetRKShifts(nep,k,array); }

1940:   PetscOptionsTail();

1942:   if (!ctx->ksp) { NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp); }
1943:   for (i=0;i<ctx->nshiftsw;i++) {
1944:     KSPSetFromOptions(ctx->ksp[i]);
1945:   }

1947:   if (ctx->fullbasis) {
1948:     if (!ctx->eps) { NEPNLEIGSGetEPS(nep,&ctx->eps); }
1949:     EPSSetFromOptions(ctx->eps);
1950:   }
1951:   return(0);
1952: }

1954: PetscErrorCode NEPView_NLEIGS(NEP nep,PetscViewer viewer)
1955: {
1957:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
1958:   PetscBool      isascii;
1959:   PetscInt       i;
1960:   char           str[50];

1963:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
1964:   if (isascii) {
1965:     PetscViewerASCIIPrintf(viewer,"  %d%% of basis vectors kept after restart\n",(int)(100*ctx->keep));
1966:     if (ctx->fullbasis) {
1967:       PetscViewerASCIIPrintf(viewer,"  using the full-basis variant\n");
1968:     } else {
1969:       PetscViewerASCIIPrintf(viewer,"  using the %slocking variant\n",ctx->lock?"":"non-");
1970:     }
1971:     PetscViewerASCIIPrintf(viewer,"  divided difference terms: used=%D, max=%D\n",ctx->nmat,ctx->ddmaxit);
1972:     PetscViewerASCIIPrintf(viewer,"  tolerance for divided difference convergence: %g\n",(double)ctx->ddtol);
1973:     if (ctx->nshifts) {
1974:       PetscViewerASCIIPrintf(viewer,"  RK shifts: ");
1975:       PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1976:       for (i=0;i<ctx->nshifts;i++) {
1977:         SlepcSNPrintfScalar(str,sizeof(str),ctx->shifts[i],PETSC_FALSE);
1978:         PetscViewerASCIIPrintf(viewer,"%s%s",str,(i<ctx->nshifts-1)?",":"");
1979:       }
1980:       PetscViewerASCIIPrintf(viewer,"\n");
1981:       PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1982:     }
1983:     if (!ctx->ksp) { NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp); }
1984:     PetscViewerASCIIPushTab(viewer);
1985:     KSPView(ctx->ksp[0],viewer);
1986:     PetscViewerASCIIPopTab(viewer);
1987:     if (ctx->fullbasis) {
1988:       if (!ctx->eps) { NEPNLEIGSGetEPS(nep,&ctx->eps); }
1989:       PetscViewerASCIIPushTab(viewer);
1990:       EPSView(ctx->eps,viewer);
1991:       PetscViewerASCIIPopTab(viewer);
1992:     }
1993:   }
1994:   return(0);
1995: }

1997: PetscErrorCode NEPReset_NLEIGS(NEP nep)
1998: {
2000:   PetscInt       k;
2001:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;

2004:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
2005:     PetscFree(ctx->coeffD);
2006:   } else {
2007:     for (k=0;k<ctx->nmat;k++) { MatDestroy(&ctx->D[k]); }
2008:   }
2009:   PetscFree4(ctx->s,ctx->xi,ctx->beta,ctx->D);
2010:   for (k=0;k<ctx->nshiftsw;k++) { KSPReset(ctx->ksp[k]); }
2011:   if (ctx->vrn) {
2012:     VecDestroy(&ctx->vrn);
2013:   }
2014:   if (ctx->fullbasis) {
2015:     MatDestroy(&ctx->A);
2016:     EPSReset(ctx->eps);
2017:     for (k=0;k<4;k++) { VecDestroy(&ctx->w[k]); }
2018:   }
2019:   return(0);
2020: }

2022: PetscErrorCode NEPDestroy_NLEIGS(NEP nep)
2023: {
2025:   PetscInt       k;
2026:   NEP_NLEIGS     *ctx = (NEP_NLEIGS*)nep->data;

2029:   BVDestroy(&ctx->V);
2030:   for (k=0;k<ctx->nshiftsw;k++) { KSPDestroy(&ctx->ksp[k]); }
2031:   PetscFree(ctx->ksp);
2032:   if (ctx->nshifts) { PetscFree(ctx->shifts); }
2033:   if (ctx->fullbasis) { EPSDestroy(&ctx->eps); }
2034:   PetscFree(nep->data);
2035:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetSingularitiesFunction_C",NULL);
2036:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetSingularitiesFunction_C",NULL);
2037:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRestart_C",NULL);
2038:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRestart_C",NULL);
2039:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetLocking_C",NULL);
2040:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetLocking_C",NULL);
2041:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetInterpolation_C",NULL);
2042:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetInterpolation_C",NULL);
2043:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRKShifts_C",NULL);
2044:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRKShifts_C",NULL);
2045:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetKSPs_C",NULL);
2046:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetFullBasis_C",NULL);
2047:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetFullBasis_C",NULL);
2048:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetEPS_C",NULL);
2049:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetEPS_C",NULL);
2050:   return(0);
2051: }

2053: SLEPC_EXTERN PetscErrorCode NEPCreate_NLEIGS(NEP nep)
2054: {
2056:   NEP_NLEIGS     *ctx;

2059:   PetscNewLog(nep,&ctx);
2060:   nep->data  = (void*)ctx;
2061:   ctx->lock  = PETSC_TRUE;
2062:   ctx->ddtol = PETSC_DEFAULT;

2064:   nep->useds = PETSC_TRUE;

2066:   nep->ops->setup          = NEPSetUp_NLEIGS;
2067:   nep->ops->setfromoptions = NEPSetFromOptions_NLEIGS;
2068:   nep->ops->view           = NEPView_NLEIGS;
2069:   nep->ops->destroy        = NEPDestroy_NLEIGS;
2070:   nep->ops->reset          = NEPReset_NLEIGS;

2072:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetSingularitiesFunction_C",NEPNLEIGSSetSingularitiesFunction_NLEIGS);
2073:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetSingularitiesFunction_C",NEPNLEIGSGetSingularitiesFunction_NLEIGS);
2074:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRestart_C",NEPNLEIGSSetRestart_NLEIGS);
2075:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRestart_C",NEPNLEIGSGetRestart_NLEIGS);
2076:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetLocking_C",NEPNLEIGSSetLocking_NLEIGS);
2077:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetLocking_C",NEPNLEIGSGetLocking_NLEIGS);
2078:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetInterpolation_C",NEPNLEIGSSetInterpolation_NLEIGS);
2079:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetInterpolation_C",NEPNLEIGSGetInterpolation_NLEIGS);
2080:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRKShifts_C",NEPNLEIGSSetRKShifts_NLEIGS);
2081:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRKShifts_C",NEPNLEIGSGetRKShifts_NLEIGS);
2082:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetKSPs_C",NEPNLEIGSGetKSPs_NLEIGS);
2083:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetFullBasis_C",NEPNLEIGSSetFullBasis_NLEIGS);
2084:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetFullBasis_C",NEPNLEIGSGetFullBasis_NLEIGS);
2085:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetEPS_C",NEPNLEIGSSetEPS_NLEIGS);
2086:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetEPS_C",NEPNLEIGSGetEPS_NLEIGS);
2087:   return(0);
2088: }