Actual source code: nleigs.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2019, Universitat Politecnica de Valencia, Spain

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

 13:    Method: NLEIGS

 15:    Algorithm:

 17:        Fully rational Krylov method for nonlinear eigenvalue problems.

 19:    References:

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

 26: #include <slepc/private/nepimpl.h>         /*I "slepcnep.h" I*/
 27: #include <slepcblaslapack.h>
 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,*work;
 62:   PetscBLASInt   n_,info,lwork;
 63:   PetscInt       i;
 64: #if defined(PETSC_USE_COMPLEX)
 65:   PetscReal      *rwork;
 66: #endif
 67: #if defined(PETSC_HAVE_ESSL)
 68:   PetscScalar    sdummy;
 69:   PetscBLASInt   idummy,io=0;
 70:   PetscScalar    *wri;
 71: #endif

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

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

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

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

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

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

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

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

216: /*  Adaptive Anderson-Antoulas algorithm */
217: static PetscErrorCode NEPNLEIGSAAAComputation(NEP nep,PetscInt ndpt,PetscScalar *ds,PetscScalar *F,PetscInt *ndptx,PetscScalar *dxi)
218: {
219: #if defined(SLEPC_MISSING_LAPACK_GGEV) || defined(PETSC_MISSING_LAPACK_GESVD)
221:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GEEV/GESVD - Lapack routines are unavailable");
222: #else
224:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
225:   PetscScalar    mean=0.0,*z,*f,*C,*A,*VT,*work,*ww,szero=0.0,sone=1.0;
226:   PetscScalar    *N,*D;
227:   PetscReal      *S,norm,err,*R;
228:   PetscInt       i,k,j,idx=0,cont;
229:   PetscBLASInt   n_,m_,lda_,lwork,info,one=1;
230: #if defined(PETSC_USE_COMPLEX)
231:   PetscReal      *rwork;
232: #endif

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

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

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

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

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

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

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

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

444: PetscErrorCode NEPNLEIGSEvalNRTFunct(NEP nep,PetscInt k,PetscScalar sigma,PetscScalar *b)
445: {
446:   NEP_NLEIGS  *ctx=(NEP_NLEIGS*)nep->data;
447:   PetscInt    i;
448:   PetscScalar *beta=ctx->beta,*s=ctx->s,*xi=ctx->xi;

451:   b[0] = 1.0/beta[0];
452:   for (i=0;i<k;i++) {
453:     b[i+1] = ((sigma-s[i])*b[i])/(beta[i+1]*(1.0-sigma/xi[i]));
454:   }
455:   return(0);
456: }

458: static PetscErrorCode MatMult_Fun(Mat A,Vec x,Vec y)
459: {
461:   ShellMatCtx    *ctx;
462:   PetscInt       i;

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

475: static PetscErrorCode MatMultTranspose_Fun(Mat A,Vec x,Vec y)
476: {
478:   ShellMatCtx    *ctx;
479:   PetscInt       i;

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

492: static PetscErrorCode MatGetDiagonal_Fun(Mat A,Vec diag)
493: {
495:   ShellMatCtx    *ctx;
496:   PetscInt       i;

499:   MatShellGetContext(A,(void**)&ctx);
500:   MatGetDiagonal(ctx->A[0],diag);
501:   if (ctx->coeff[0]!=1.0) { VecScale(diag,ctx->coeff[0]); }
502:   for (i=1;i<ctx->nmat;i++) {
503:     MatGetDiagonal(ctx->A[i],ctx->t);
504:     VecAXPY(diag,ctx->coeff[i],ctx->t);
505:   }
506:   return(0);
507: }

509: static PetscErrorCode MatDuplicate_Fun(Mat A,MatDuplicateOption op,Mat *B)
510: {
511:   PetscInt       n,i;
512:   ShellMatCtx    *ctxnew,*ctx;
513:   void           (*fun)(void);

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

546: static PetscErrorCode MatDestroy_Fun(Mat A)
547: {
548:   ShellMatCtx    *ctx;
550:   PetscInt       i;

553:   if (A) {
554:     MatShellGetContext(A,(void**)&ctx);
555:     for (i=0;i<ctx->nmat;i++) {
556:       MatDestroy(&ctx->A[i]);
557:     }
558:     VecDestroy(&ctx->t);
559:     PetscFree2(ctx->A,ctx->coeff);
560:     PetscFree(ctx);
561:   }
562:   return(0);
563: }

565: static PetscErrorCode MatAXPY_Fun(Mat Y,PetscScalar a,Mat X,MatStructure str)
566: {
567:   ShellMatCtx    *ctxY,*ctxX;
569:   PetscInt       i,j;
570:   PetscBool      found;

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

592: static PetscErrorCode MatScale_Fun(Mat M,PetscScalar a)
593: {
594:   ShellMatCtx    *ctx;
596:   PetscInt       i;

599:   MatShellGetContext(M,(void**)&ctx);
600:   for (i=0;i<ctx->nmat;i++) ctx->coeff[i] *= a;
601:   return(0);
602: }

604: static PetscErrorCode NLEIGSMatToMatShellArray(Mat M,Mat *Ms,PetscInt maxnmat)
605: {
607:   ShellMatCtx    *ctx;
608:   PetscInt       n;
609:   PetscBool      has;

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

634: static PetscErrorCode NEPNLEIGSNormEstimation(NEP nep,Mat M,PetscReal *norm,Vec *w)
635: {
636:   PetscScalar    *z,*x,*y;
637:   PetscReal      tr;
638:   Vec            X=w[0],Y=w[1];
639:   PetscInt       n,i;
640:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
641:   PetscRandom    rand;

645:   if (!ctx->vrn) {
646:     /* generate a random vector with normally distributed entries with the Box-Muller transform */
647:     BVGetRandomContext(nep->V,&rand);
648:     MatCreateVecs(M,&ctx->vrn,NULL);
649:     VecSetRandom(X,rand);
650:     VecSetRandom(Y,rand);
651:     VecGetLocalSize(ctx->vrn,&n);
652:     VecGetArray(ctx->vrn,&z);
653:     VecGetArray(X,&x);
654:     VecGetArray(Y,&y);
655:     for (i=0;i<n;i++) {
656: #if defined(PETSC_USE_COMPLEX)
657:       z[i] = PetscSqrtReal(-2.0*PetscLogReal(PetscRealPart(x[i])))*PetscCosReal(2.0*PETSC_PI*PetscRealPart(y[i]));
658:       z[i] += PETSC_i*(PetscSqrtReal(-2.0*PetscLogReal(PetscImaginaryPart(x[i])))*PetscCosReal(2.0*PETSC_PI*PetscImaginaryPart(y[i])));
659: #else
660:       z[i] = PetscSqrtReal(-2.0*PetscLogReal(x[i]))*PetscCosReal(2.0*PETSC_PI*y[i]);
661: #endif
662:     }
663:     VecRestoreArray(ctx->vrn,&z);
664:     VecRestoreArray(X,&x);
665:     VecRestoreArray(Y,&y);
666:     VecNorm(ctx->vrn,NORM_2,&tr);
667:     VecScale(ctx->vrn,1/tr);
668:   }
669:   /* matrix-free norm estimator of Ipsen http://www4.ncsu.edu/~ipsen/ps/slides_ima.pdf */
670:   MatGetSize(M,&n,NULL);
671:   MatMult(M,ctx->vrn,X);
672:   VecNorm(X,NORM_2,norm);
673:   *norm *= PetscSqrtReal((PetscReal)n);
674:   return(0);
675: }

677: static PetscErrorCode NEPNLEIGSDividedDifferences_split(NEP nep)
678: {
680:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
681:   PetscInt       k,j,i,maxnmat,nmax;
682:   PetscReal      norm0,norm,*matnorm;
683:   PetscScalar    *s=ctx->s,*beta=ctx->beta,*xi=ctx->xi,*b,alpha,*coeffs,*pK,*pH,sone=1.0;
684:   Mat            T,Ts,K,H;
685:   PetscBool      shell,hasmnorm=PETSC_FALSE,matrix=PETSC_TRUE;
686:   PetscBLASInt   n_;

689:   nmax = ctx->ddmaxit;
690:   PetscMalloc1(nep->nt*nmax,&ctx->coeffD);
691:   PetscMalloc3(nmax+1,&b,nmax+1,&coeffs,nep->nt,&matnorm);
692:   for (j=0;j<nep->nt;j++) {
693:     MatHasOperation(nep->A[j],MATOP_NORM,&hasmnorm);
694:     if (!hasmnorm) break;
695:     MatNorm(nep->A[j],NORM_INFINITY,matnorm+j);
696:   }
697:   /* Try matrix functions scheme */
698:   PetscCalloc2(nmax*nmax,&pK,nmax*nmax,&pH);
699:   for (i=0;i<nmax-1;i++) {
700:     pK[(nmax+1)*i]   = 1.0;
701:     pK[(nmax+1)*i+1] = beta[i+1]/xi[i];
702:     pH[(nmax+1)*i]   = s[i];
703:     pH[(nmax+1)*i+1] = beta[i+1];
704:   }
705:   pH[nmax*nmax-1] = s[nmax-1];
706:   pK[nmax*nmax-1] = 1.0;
707:   PetscBLASIntCast(nmax,&n_);
708:   PetscStackCallBLAS("BLAStrsm",BLAStrsm_("R","L","N","U",&n_,&n_,&sone,pK,&n_,pH,&n_));
709:   /* The matrix to be used is in H. K will be a work-space matrix */
710:   MatCreateSeqDense(PETSC_COMM_SELF,nmax,nmax,pH,&H);
711:   MatCreateSeqDense(PETSC_COMM_SELF,nmax,nmax,pK,&K);
712:   for (j=0;matrix&&j<nep->nt;j++) {
713:     PetscPushErrorHandler(PetscIgnoreErrorHandler,NULL);
714:     FNEvaluateFunctionMat(nep->f[j],H,K);
715:     PetscPopErrorHandler();
716:     if (!ierr) { 
717:       for (i=0;i<nmax;i++) { ctx->coeffD[j+i*nep->nt] = pK[i]*beta[0]; }
718:     } else {
719:       matrix = PETSC_FALSE;
720:       PetscFPTrapPop();
721:     }
722:   }
723:   MatDestroy(&H);
724:   MatDestroy(&K);
725:   if (!matrix) {
726:     for (j=0;j<nep->nt;j++) {
727:       FNEvaluateFunction(nep->f[j],s[0],ctx->coeffD+j);
728:       ctx->coeffD[j] *= beta[0];
729:     }
730:   }
731:   if (hasmnorm) {
732:     norm0 = 0.0;
733:     for (j=0;j<nep->nt;j++) norm0 += matnorm[j]*PetscAbsScalar(ctx->coeffD[j]);
734:   } else {
735:     norm0 = 0.0;
736:     for (j=0;j<nep->nt;j++) norm0 = PetscMax(PetscAbsScalar(ctx->coeffD[j]),norm0);
737:   }
738:   ctx->nmat = ctx->ddmaxit;
739:   for (k=1;k<ctx->ddmaxit;k++) {
740:     if (!matrix) {
741:       NEPNLEIGSEvalNRTFunct(nep,k,s[k],b);
742:       for (i=0;i<nep->nt;i++) {
743:         FNEvaluateFunction(nep->f[i],s[k],ctx->coeffD+k*nep->nt+i);
744:         for (j=0;j<k;j++) {
745:           ctx->coeffD[k*nep->nt+i] -= b[j]*ctx->coeffD[i+nep->nt*j];
746:         }
747:         ctx->coeffD[k*nep->nt+i] /= b[k];
748:       }
749:     }
750:     if (hasmnorm) {
751:       norm = 0.0;
752:       for (j=0;j<nep->nt;j++) norm += matnorm[j]*PetscAbsScalar(ctx->coeffD[k*nep->nt+j]);
753:     } else {
754:       norm = 0.0;
755:       for (j=0;j<nep->nt;j++) norm = PetscMax(PetscAbsScalar(ctx->coeffD[k*nep->nt+j]),norm);
756:     }
757:     if (norm/norm0 < ctx->ddtol) {
758:       ctx->nmat = k+1;
759:       break;
760:     }
761:   }
762:   if (!ctx->ksp) { NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp); }
763:   PetscObjectTypeCompare((PetscObject)nep->A[0],MATSHELL,&shell);
764:   maxnmat = PetscMax(ctx->ddmaxit,nep->nt);
765:   for (i=0;i<ctx->nshiftsw;i++) {
766:     NEPNLEIGSEvalNRTFunct(nep,ctx->nmat-1,ctx->shifts[i],coeffs);
767:     if (!shell) {
768:       MatDuplicate(nep->A[0],MAT_COPY_VALUES,&T);
769:     } else {
770:       NLEIGSMatToMatShellArray(nep->A[0],&T,maxnmat);
771:     }
772:     alpha = 0.0;
773:     for (j=0;j<ctx->nmat;j++) alpha += coeffs[j]*ctx->coeffD[j*nep->nt];
774:     MatScale(T,alpha);
775:     for (k=1;k<nep->nt;k++) {
776:       alpha = 0.0;
777:       for (j=0;j<ctx->nmat;j++) alpha += coeffs[j]*ctx->coeffD[j*nep->nt+k];
778:       if (shell) {
779:         NLEIGSMatToMatShellArray(nep->A[k],&Ts,maxnmat);
780:       }
781:       MatAXPY(T,alpha,shell?Ts:nep->A[k],nep->mstr);
782:       if (shell) {
783:         MatDestroy(&Ts);
784:       }
785:     }
786:     KSPSetOperators(ctx->ksp[i],T,T);
787:     KSPSetUp(ctx->ksp[i]);
788:     MatDestroy(&T);
789:   }
790:   PetscFree3(b,coeffs,matnorm);
791:   PetscFree2(pK,pH);
792:   return(0);
793: }

795: static PetscErrorCode NEPNLEIGSDividedDifferences_callback(NEP nep)
796: {
798:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
799:   PetscInt       k,j,i,maxnmat;
800:   PetscReal      norm0,norm;
801:   PetscScalar    *s=ctx->s,*beta=ctx->beta,*b,*coeffs;
802:   Mat            *D=ctx->D,T;
803:   PetscBool      shell,has,vec=PETSC_FALSE;
804:   Vec            w[2];

807:   PetscMalloc2(ctx->ddmaxit+1,&b,ctx->ddmaxit+1,&coeffs);
808:   T = nep->function;
809:   NEPComputeFunction(nep,s[0],T,T);
810:   PetscObjectTypeCompare((PetscObject)T,MATSHELL,&shell);
811:   maxnmat = PetscMax(ctx->ddmaxit,nep->nt);
812:   if (!shell) {
813:     MatDuplicate(T,MAT_COPY_VALUES,&D[0]);
814:   } else {
815:     NLEIGSMatToMatShellArray(T,&D[0],maxnmat);
816:   }
817:   if (beta[0]!=1.0) {
818:     MatScale(D[0],1.0/beta[0]);
819:   }
820:   MatHasOperation(D[0],MATOP_NORM,&has);
821:   if (has) {
822:     MatNorm(D[0],NORM_FROBENIUS,&norm0);
823:   } else {
824:     MatCreateVecs(D[0],NULL,&w[0]);
825:     VecDuplicate(w[0],&w[1]);
826:     vec = PETSC_TRUE;
827:     NEPNLEIGSNormEstimation(nep,D[0],&norm0,w);
828:   }
829:   ctx->nmat = ctx->ddmaxit;
830:   for (k=1;k<ctx->ddmaxit;k++) {
831:     NEPNLEIGSEvalNRTFunct(nep,k,s[k],b);
832:     NEPComputeFunction(nep,s[k],T,T);
833:     if (!shell) {
834:       MatDuplicate(T,MAT_COPY_VALUES,&D[k]);
835:     } else {
836:       NLEIGSMatToMatShellArray(T,&D[k],maxnmat);
837:     }
838:     for (j=0;j<k;j++) {
839:       MatAXPY(D[k],-b[j],D[j],nep->mstr);
840:     }
841:     MatScale(D[k],1.0/b[k]);
842:     MatHasOperation(D[k],MATOP_NORM,&has);
843:     if (has) {
844:       MatNorm(D[k],NORM_FROBENIUS,&norm);
845:     } else {
846:       if(!vec) {
847:         MatCreateVecs(D[k],NULL,&w[0]);
848:         VecDuplicate(w[0],&w[1]);
849:         vec = PETSC_TRUE;
850:       }
851:       NEPNLEIGSNormEstimation(nep,D[k],&norm,w);
852:     }
853:     if (norm/norm0 < ctx->ddtol && k>1) {
854:       ctx->nmat = k+1;
855:       break;
856:     }
857:   }
858:   if (!ctx->ksp) { NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp); }
859:   for (i=0;i<ctx->nshiftsw;i++) {
860:     NEPNLEIGSEvalNRTFunct(nep,ctx->nmat-1,ctx->shifts[i],coeffs);
861:     MatDuplicate(ctx->D[0],MAT_COPY_VALUES,&T);
862:     if (coeffs[0]!=1.0) { MatScale(T,coeffs[0]); }
863:     for (j=1;j<ctx->nmat;j++) {
864:       MatAXPY(T,coeffs[j],ctx->D[j],nep->mstr);
865:     }
866:     KSPSetOperators(ctx->ksp[i],T,T);
867:     KSPSetUp(ctx->ksp[i]);
868:     MatDestroy(&T);
869:   }
870:   PetscFree2(b,coeffs);
871:   if (vec) {
872:     VecDestroy(&w[0]);
873:     VecDestroy(&w[1]);
874:   }
875:   return(0);
876: }

878: /*
879:    NEPKrylovConvergence - This is the analogue to EPSKrylovConvergence.
880: */
881: PetscErrorCode NEPNLEIGSKrylovConvergence(NEP nep,PetscBool getall,PetscInt kini,PetscInt nits,PetscReal betah,PetscScalar betak,PetscInt *kout,Vec *w)
882: {
884:   PetscInt       k,newk,marker,inside;
885:   PetscScalar    re,im;
886:   PetscReal      resnorm,tt;
887:   PetscBool      istrivial;
888:   NEP_NLEIGS     *ctx = (NEP_NLEIGS*)nep->data;

891:   RGIsTrivial(nep->rg,&istrivial);
892:   marker = -1;
893:   if (nep->trackall) getall = PETSC_TRUE;
894:   for (k=kini;k<kini+nits;k++) {
895:     /* eigenvalue */
896:     re = nep->eigr[k];
897:     im = nep->eigi[k];
898:     if (!istrivial) {
899:       if (!ctx->nshifts) {
900:         NEPNLEIGSBackTransform((PetscObject)nep,1,&re,&im);
901:       }
902:       RGCheckInside(nep->rg,1,&re,&im,&inside);
903:       if (marker==-1 && inside<0) marker = k;
904:     }
905:     newk = k;
906:     DSVectors(nep->ds,DS_MAT_X,&newk,&resnorm);
907:     tt = (ctx->nshifts)?SlepcAbsEigenvalue(betak-nep->eigr[k]*betah,nep->eigi[k]*betah):betah;
908:     resnorm *=  PetscAbsReal(tt);
909:     /* error estimate */
910:     (*nep->converged)(nep,nep->eigr[k],nep->eigi[k],resnorm,&nep->errest[k],nep->convergedctx);
911:     if (marker==-1 && nep->errest[k] >= nep->tol) marker = k;
912:     if (newk==k+1) {
913:       nep->errest[k+1] = nep->errest[k];
914:       k++;
915:     }
916:     if (marker!=-1 && !getall) break;
917:   }
918:   if (marker!=-1) k = marker;
919:   *kout = k;
920:   return(0);
921: }

923: PetscErrorCode NEPSetUp_NLEIGS(NEP nep)
924: {
926:   PetscInt       k,in;
927:   PetscScalar    zero=0.0;
928:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
929:   SlepcSC        sc;
930:   PetscBool      istrivial;

933:   NEPSetDimensions_Default(nep,nep->nev,&nep->ncv,&nep->mpd);
934:   if (nep->ncv>nep->nev+nep->mpd) SETERRQ(PetscObjectComm((PetscObject)nep),1,"The value of ncv must not be larger than nev+mpd");
935:   if (!nep->max_it) nep->max_it = PetscMax(5000,2*nep->n/nep->ncv);
936:   if (!ctx->ddmaxit) ctx->ddmaxit = LBPOINTS;
937:   RGIsTrivial(nep->rg,&istrivial);
938:   if (istrivial) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"NEPNLEIGS requires a nontrivial region defining the target set");
939:   if (!nep->which) nep->which = NEP_TARGET_MAGNITUDE;
940:   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,"Should set a target selection in NEPSetWhichEigenvalues()");

942:   /* Initialize the NLEIGS context structure */
943:   k = ctx->ddmaxit;
944:   PetscMalloc4(k,&ctx->s,k,&ctx->xi,k,&ctx->beta,k,&ctx->D);
945:   nep->data = ctx;
946:   if (nep->tol==PETSC_DEFAULT) nep->tol = SLEPC_DEFAULT_TOL;
947:   if (ctx->ddtol==PETSC_DEFAULT) ctx->ddtol = nep->tol/10.0;
948:   if (!ctx->keep) ctx->keep = 0.5;

950:   /* Compute Leja-Bagby points and scaling values */
951:   NEPNLEIGSLejaBagbyPoints(nep);
952:   if (nep->problem_type!=NEP_RATIONAL) {
953:     RGCheckInside(nep->rg,1,&nep->target,&zero,&in);
954:     if (in<0) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"The target is not inside the target set");
955:   }

957:   /* Compute the divided difference matrices */
958:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
959:     NEPNLEIGSDividedDifferences_split(nep);
960:   } else {
961:     NEPNLEIGSDividedDifferences_callback(nep);
962:   }
963:   NEPAllocateSolution(nep,ctx->nmat-1);
964:   NEPSetWorkVecs(nep,4);
965:   if (!ctx->fullbasis) {
966:     if (nep->twosided) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Two-sided variant requires the full-basis option, rerun with -nep_nleigs_full_basis");
967:     /* set-up DS and transfer split operator functions */
968:     DSSetType(nep->ds,ctx->nshifts?DSGNHEP:DSNHEP);
969:     DSAllocate(nep->ds,nep->ncv+1);
970:     DSGetSlepcSC(nep->ds,&sc);
971:     if (!ctx->nshifts) {
972:       sc->map = NEPNLEIGSBackTransform;
973:       DSSetExtraRow(nep->ds,PETSC_TRUE);
974:     }
975:     sc->mapobj        = (PetscObject)nep;
976:     sc->rg            = nep->rg;
977:     sc->comparison    = nep->sc->comparison;
978:     sc->comparisonctx = nep->sc->comparisonctx;
979:     BVDestroy(&ctx->V);
980:     BVCreateTensor(nep->V,ctx->nmat-1,&ctx->V);
981:     nep->ops->solve          = NEPSolve_NLEIGS;
982:     nep->ops->computevectors = NEPComputeVectors_Schur;
983:   } else {
984:     NEPSetUp_NLEIGS_FullBasis(nep);
985:     nep->ops->solve          = NEPSolve_NLEIGS_FullBasis;
986:     nep->ops->computevectors = NULL;
987:   }
988:   return(0);
989: }

991: /*
992:   Extend the TOAR basis by applying the the matrix operator
993:   over a vector which is decomposed on the TOAR way
994:   Input:
995:     - S,V: define the latest Arnoldi vector (nv vectors in V)
996:   Output:
997:     - t: new vector extending the TOAR basis
998:     - r: temporally coefficients to compute the TOAR coefficients
999:          for the new Arnoldi vector
1000:   Workspace: t_ (two vectors)
1001: */
1002: 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_)
1003: {
1005:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
1006:   PetscInt       deg=ctx->nmat-1,k,j;
1007:   Vec            v=t_[0],q=t_[1],w;
1008:   PetscScalar    *beta=ctx->beta,*s=ctx->s,*xi=ctx->xi,*coeffs,sigma;

1011:   if (!ctx->ksp) { NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp); }
1012:   sigma = ctx->shifts[idxrktg];
1013:   BVSetActiveColumns(nep->V,0,nv);
1014:   PetscMalloc1(ctx->nmat,&coeffs);
1015:   if (PetscAbsScalar(s[deg-2]-sigma)<100*PETSC_MACHINE_EPSILON) SETERRQ(PETSC_COMM_SELF,1,"Breakdown in NLEIGS");
1016:   /* i-part stored in (i-1) position */
1017:   for (j=0;j<nv;j++) {
1018:     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);
1019:   }
1020:   BVSetActiveColumns(W,0,deg);
1021:   BVGetColumn(W,deg-1,&w);
1022:   BVMultVec(V,1.0/beta[deg],0,w,S+(deg-1)*ls);
1023:   BVRestoreColumn(W,deg-1,&w);
1024:   BVGetColumn(W,deg-2,&w);
1025:   BVMultVec(V,1.0,0.0,w,r+(deg-2)*lr);
1026:   BVRestoreColumn(W,deg-2,&w);
1027:   for (k=deg-2;k>0;k--) {
1028:     if (PetscAbsScalar(s[k-1]-sigma)<100*PETSC_MACHINE_EPSILON) SETERRQ(PETSC_COMM_SELF,1,"Breakdown in NLEIGS");
1029:     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);
1030:     BVGetColumn(W,k-1,&w);
1031:     BVMultVec(V,1.0,0.0,w,r+(k-1)*lr);
1032:     BVRestoreColumn(W,k-1,&w);
1033:   }
1034:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
1035:     for (j=0;j<ctx->nmat-2;j++) coeffs[j] = ctx->coeffD[nep->nt*j];
1036:     coeffs[ctx->nmat-2] = ctx->coeffD[nep->nt*(ctx->nmat-1)];
1037:     BVMultVec(W,1.0,0.0,v,coeffs);
1038:     MatMult(nep->A[0],v,q);
1039:     for (k=1;k<nep->nt;k++) {
1040:       for (j=0;j<ctx->nmat-2;j++) coeffs[j] = ctx->coeffD[nep->nt*j+k];
1041:       coeffs[ctx->nmat-2] = ctx->coeffD[nep->nt*(ctx->nmat-1)+k];
1042:       BVMultVec(W,1.0,0,v,coeffs);
1043:       MatMult(nep->A[k],v,t);
1044:       VecAXPY(q,1.0,t);
1045:     }
1046:     KSPSolve(ctx->ksp[idxrktg],q,t);
1047:     VecScale(t,-1.0);
1048:   } else {
1049:     for (k=0;k<deg-1;k++) {
1050:       BVGetColumn(W,k,&w);
1051:       MatMult(ctx->D[k],w,q);
1052:       BVRestoreColumn(W,k,&w);
1053:       BVInsertVec(W,k,q);
1054:     }
1055:     BVGetColumn(W,deg-1,&w);
1056:     MatMult(ctx->D[deg],w,q);
1057:     BVRestoreColumn(W,k,&w);
1058:     BVInsertVec(W,k,q);
1059:     for (j=0;j<ctx->nmat-1;j++) coeffs[j] = 1.0;
1060:     BVMultVec(W,1.0,0.0,q,coeffs);
1061:     KSPSolve(ctx->ksp[idxrktg],q,t);
1062:     VecScale(t,-1.0);
1063:   }
1064:   PetscFree(coeffs);
1065:   return(0);
1066: }

1068: /*
1069:   Compute TOAR coefficients of the blocks of the new Arnoldi vector computed
1070: */
1071: static PetscErrorCode NEPTOARCoefficients(NEP nep,PetscScalar sigma,PetscInt nv,PetscScalar *S,PetscInt ls,PetscScalar *r,PetscInt lr,PetscScalar *x,PetscScalar *work)
1072: {
1074:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
1075:   PetscInt       k,j,d=ctx->nmat-1;
1076:   PetscScalar    *t=work;

1079:   NEPNLEIGSEvalNRTFunct(nep,d-1,sigma,t);
1080:   for (k=0;k<d-1;k++) {
1081:     for (j=0;j<=nv;j++) r[k*lr+j] += t[k]*x[j];
1082:   }
1083:   for (j=0;j<=nv;j++) r[(d-1)*lr+j] = t[d-1]*x[j];
1084:   return(0);
1085: }

1087: /*
1088:   Compute continuation vector coefficients for the Rational-Krylov run.
1089:   dim(work) >= (end-ini)*(end-ini+1) + end+1 + 2*(end-ini+1), dim(t) = end.
1090: */
1091: 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)
1092: {
1093: #if defined(PETSC_MISSING_LAPACK_GEQRF) || defined(SLEPC_MISSING_LAPACK_LARF)
1095:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GEQRF/LARF - Lapack routines are unavailable");
1096: #else
1098:   PetscScalar    *x,*W,*tau,sone=1.0,szero=0.0;
1099:   PetscInt       i,j,n1,n,nwu=0;
1100:   PetscBLASInt   info,n_,n1_,one=1,dim,lds_;
1101:   NEP_NLEIGS     *ctx = (NEP_NLEIGS*)nep->data;

1104:   if (!ctx->nshifts || !end) {
1105:     t[0] = 1;
1106:     PetscMemcpy(cont,S+end*lds,lds*sizeof(PetscScalar));
1107:   } else {
1108:     n   = end-ini;
1109:     n1  = n+1;
1110:     x   = work+nwu;
1111:     nwu += end+1;
1112:     tau = work+nwu;
1113:     nwu += n;
1114:     W   = work+nwu;
1115:     nwu += n1*n;
1116:     for (j=ini;j<end;j++) {
1117:       for (i=ini;i<=end;i++) W[(j-ini)*n1+i-ini] = K[j*ld+i] -H[j*ld+i]*sigma;
1118:     }
1119:     PetscBLASIntCast(n,&n_);
1120:     PetscBLASIntCast(n1,&n1_);
1121:     PetscBLASIntCast(end+1,&dim);
1122:     PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&n1_,&n_,W,&n1_,tau,work+nwu,&n1_,&info));
1123:     SlepcCheckLapackInfo("geqrf",info);
1124:     for (i=0;i<end;i++) t[i] = 0.0;
1125:     t[end] = 1.0;
1126:     for (j=n-1;j>=0;j--) {
1127:       for (i=0;i<ini+j;i++) x[i] = 0.0;
1128:       x[ini+j] = 1.0;
1129:       for (i=j+1;i<n1;i++) x[i+ini] = W[i+n1*j];
1130:       tau[j] = PetscConj(tau[j]);
1131:       PetscStackCallBLAS("LAPACKlarf",LAPACKlarf_("L",&dim,&one,x,&one,tau+j,t,&dim,work+nwu));
1132:     }
1133:     PetscBLASIntCast(lds,&lds_);
1134:     PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&lds_,&n1_,&sone,S,&lds_,t,&one,&szero,cont,&one));
1135:   }
1136:   return(0);
1137: #endif
1138: }

1140: /*
1141:   Compute a run of Arnoldi iterations
1142: */
1143: PetscErrorCode NEPNLEIGSTOARrun(NEP nep,PetscScalar *K,PetscScalar *H,PetscInt ldh,BV W,PetscInt k,PetscInt *M,PetscBool *breakdown,Vec *t_)
1144: {
1146:   NEP_NLEIGS     *ctx = (NEP_NLEIGS*)nep->data;
1147:   PetscInt       i,j,m=*M,lwa,deg=ctx->nmat-1,lds,nqt,ld,l;
1148:   Vec            t;
1149:   PetscReal      norm;
1150:   PetscScalar    *x,*work,*tt,sigma,*cont,*S;
1151:   PetscBool      lindep;
1152:   Mat            MS;

1155:   BVTensorGetFactors(ctx->V,NULL,&MS);
1156:   MatDenseGetArray(MS,&S);
1157:   BVGetSizes(nep->V,NULL,NULL,&ld);
1158:   lds = ld*deg;
1159:   BVGetActiveColumns(nep->V,&l,&nqt);
1160:   lwa = PetscMax(ld,deg)+(m+1)*(m+1)+4*(m+1);
1161:   PetscMalloc4(ld,&x,lwa,&work,m+1,&tt,lds,&cont);
1162:   BVSetActiveColumns(ctx->V,0,m);
1163:   for (j=k;j<m;j++) {
1164:     sigma = ctx->shifts[(++(ctx->idxrk))%ctx->nshiftsw];

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

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

1174:     /* orthogonalize */
1175:     BVOrthogonalizeColumn(nep->V,nqt,x,&norm,&lindep);
1176:     if (!lindep) {
1177:       x[nqt] = norm;
1178:       BVScaleColumn(nep->V,nqt,1.0/norm);
1179:       nqt++;
1180:     } else x[nqt] = 0.0;

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

1184:     /* Level-2 orthogonalization */
1185:     BVOrthogonalizeColumn(ctx->V,j+1,H+j*ldh,&norm,breakdown);
1186:     H[j+1+ldh*j] = norm;
1187:     if (ctx->nshifts) {
1188:       for (i=0;i<=j;i++) K[i+ldh*j] = sigma*H[i+ldh*j] + tt[i];
1189:       K[j+1+ldh*j] = sigma*H[j+1+ldh*j];
1190:     }
1191:     if (*breakdown) {
1192:       *M = j+1;
1193:       break;
1194:     }
1195:     BVScaleColumn(ctx->V,j+1,1.0/norm);
1196:     BVSetActiveColumns(nep->V,l,nqt);
1197:   }
1198:   PetscFree4(x,work,tt,cont);
1199:   MatDenseRestoreArray(MS,&S);
1200:   BVTensorRestoreFactors(ctx->V,NULL,&MS);
1201:   return(0);
1202: }

1204: PetscErrorCode NEPSolve_NLEIGS(NEP nep)
1205: {
1207:   NEP_NLEIGS     *ctx = (NEP_NLEIGS*)nep->data;
1208:   PetscInt       i,k=0,l,nv=0,ld,lds,ldds,nq,newn;
1209:   PetscInt       deg=ctx->nmat-1,nconv=0;
1210:   PetscScalar    *S,*Q,*H,*pU,*K,betak=0,*eigr,*eigi;
1211:   PetscReal      betah;
1212:   PetscBool      falselock=PETSC_FALSE,breakdown=PETSC_FALSE;
1213:   BV             W;
1214:   Mat            MS,MQ,U;

1217:   if (ctx->lock) {
1218:     PetscOptionsGetBool(NULL,NULL,"-nep_nleigs_falselocking",&falselock,NULL);
1219:   }

1221:   BVGetSizes(nep->V,NULL,NULL,&ld);
1222:   lds = deg*ld;
1223:   DSGetLeadingDimension(nep->ds,&ldds);
1224:   if (!ctx->nshifts) {
1225:     PetscMalloc2(nep->ncv,&eigr,nep->ncv,&eigi);
1226:   } else { eigr = nep->eigr; eigi = nep->eigi; }
1227:   BVDuplicateResize(nep->V,PetscMax(nep->nt-1,ctx->nmat-1),&W);


1230:   /* clean projected matrix (including the extra-arrow) */
1231:   DSGetArray(nep->ds,DS_MAT_A,&H);
1232:   PetscMemzero(H,ldds*ldds*sizeof(PetscScalar));
1233:   DSRestoreArray(nep->ds,DS_MAT_A,&H);
1234:   if (ctx->nshifts) {
1235:     DSGetArray(nep->ds,DS_MAT_B,&H);
1236:     PetscMemzero(H,ldds*ldds*sizeof(PetscScalar));
1237:     DSRestoreArray(nep->ds,DS_MAT_B,&H);
1238:   }
1239:   
1240:   /* Get the starting Arnoldi vector */
1241:   BVTensorBuildFirstColumn(ctx->V,nep->nini);
1242:   
1243:   /* Restart loop */
1244:   l = 0;
1245:   while (nep->reason == NEP_CONVERGED_ITERATING) {
1246:     nep->its++;

1248:     /* Compute an nv-step Krylov relation */
1249:     nv = PetscMin(nep->nconv+nep->mpd,nep->ncv);
1250:     if (ctx->nshifts) { DSGetArray(nep->ds,DS_MAT_A,&K); }
1251:     DSGetArray(nep->ds,ctx->nshifts?DS_MAT_B:DS_MAT_A,&H);
1252:     NEPNLEIGSTOARrun(nep,K,H,ldds,W,nep->nconv+l,&nv,&breakdown,nep->work);
1253:     betah = PetscAbsScalar(H[(nv-1)*ldds+nv]);
1254:     DSRestoreArray(nep->ds,ctx->nshifts?DS_MAT_B:DS_MAT_A,&H);
1255:     if (ctx->nshifts) {
1256:       betak = K[(nv-1)*ldds+nv];
1257:       DSRestoreArray(nep->ds,DS_MAT_A,&K);
1258:     }
1259:     DSSetDimensions(nep->ds,nv,0,nep->nconv,nep->nconv+l);
1260:     if (l==0) {
1261:       DSSetState(nep->ds,DS_STATE_INTERMEDIATE);
1262:     } else {
1263:       DSSetState(nep->ds,DS_STATE_RAW);
1264:     }

1266:     /* Solve projected problem */
1267:     DSSolve(nep->ds,nep->eigr,nep->eigi);
1268:     DSSort(nep->ds,nep->eigr,nep->eigi,NULL,NULL,NULL);
1269:     if (!ctx->nshifts) {
1270:       DSUpdateExtraRow(nep->ds);
1271:     }
1272:     DSSynchronize(nep->ds,nep->eigr,nep->eigi);

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

1278:     /* Update l */
1279:     if (nep->reason != NEP_CONVERGED_ITERATING || breakdown) l = 0;
1280:     else {
1281:       l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
1282:       if (!breakdown) {
1283:         /* Prepare the Rayleigh quotient for restart */
1284:         if (!ctx->nshifts) {
1285:           DSTruncate(nep->ds,k+l);
1286:           DSGetDimensions(nep->ds,&newn,NULL,NULL,NULL,NULL);
1287:           l = newn-k;
1288:         } else {
1289:           DSGetArray(nep->ds,DS_MAT_Q,&Q);
1290:           DSGetArray(nep->ds,DS_MAT_B,&H);
1291:           DSGetArray(nep->ds,DS_MAT_A,&K);
1292:           for (i=ctx->lock?k:0;i<k+l;i++) {
1293:             H[k+l+i*ldds] = betah*Q[nv-1+i*ldds];
1294:             K[k+l+i*ldds] = betak*Q[nv-1+i*ldds];
1295:           }
1296:           DSRestoreArray(nep->ds,DS_MAT_B,&H);
1297:           DSRestoreArray(nep->ds,DS_MAT_A,&K);
1298:           DSRestoreArray(nep->ds,DS_MAT_Q,&Q);
1299:         }
1300:       }
1301:     }
1302:     nconv = k;
1303:     if (!ctx->lock && nep->reason == NEP_CONVERGED_ITERATING && !breakdown) { l += k; k = 0; }

1305:     /* Update S */
1306:     DSGetMat(nep->ds,ctx->nshifts?DS_MAT_Z:DS_MAT_Q,&MQ);
1307:     BVMultInPlace(ctx->V,MQ,nep->nconv,k+l);
1308:     MatDestroy(&MQ);

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

1313:     if (breakdown && nep->reason == NEP_CONVERGED_ITERATING) {
1314:       /* Stop if breakdown */
1315:       PetscInfo2(nep,"Breakdown (it=%D norm=%g)\n",nep->its,(double)betah);
1316:       nep->reason = NEP_DIVERGED_BREAKDOWN;
1317:     }
1318:     if (nep->reason != NEP_CONVERGED_ITERATING) l--;
1319:     /* truncate S */
1320:     BVGetActiveColumns(nep->V,NULL,&nq);
1321:     if (k+l+deg<=nq) {
1322:       BVSetActiveColumns(ctx->V,nep->nconv,k+l+1);
1323:       if (!falselock && ctx->lock) {
1324:         BVTensorCompress(ctx->V,k-nep->nconv);
1325:       } else {
1326:         BVTensorCompress(ctx->V,0);
1327:       }
1328:     }
1329:     nep->nconv = k;
1330:     if (!ctx->nshifts) {
1331:       for (i=0;i<nv;i++) { eigr[i] = nep->eigr[i]; eigi[i] = nep->eigi[i]; }
1332:       NEPNLEIGSBackTransform((PetscObject)nep,nv,eigr,eigi);
1333:     }
1334:     NEPMonitor(nep,nep->its,nconv,eigr,eigi,nep->errest,nv);
1335:   }
1336:   nep->nconv = nconv;
1337:   if (nep->nconv>0) {
1338:     BVSetActiveColumns(ctx->V,0,nep->nconv);
1339:     BVGetActiveColumns(nep->V,NULL,&nq);
1340:     BVSetActiveColumns(nep->V,0,nq);
1341:     if (nq>nep->nconv) {
1342:       BVTensorCompress(ctx->V,nep->nconv);
1343:       BVSetActiveColumns(nep->V,0,nep->nconv);
1344:       nq = nep->nconv;
1345:     }
1346:     if (ctx->nshifts) {
1347:       DSGetMat(nep->ds,DS_MAT_B,&MQ);
1348:       BVMultInPlace(ctx->V,MQ,0,nep->nconv);
1349:       MatDestroy(&MQ);
1350:     }
1351:     BVTensorGetFactors(ctx->V,NULL,&MS);
1352:     MatDenseGetArray(MS,&S);
1353:     PetscMalloc1(nq*nep->nconv,&pU);
1354:     for (i=0;i<nep->nconv;i++) {
1355:       PetscMemcpy(pU+i*nq,S+i*lds,nq*sizeof(PetscScalar));
1356:     }
1357:     MatDenseRestoreArray(MS,&S);
1358:     BVTensorRestoreFactors(ctx->V,NULL,&MS);
1359:     MatCreateSeqDense(PETSC_COMM_SELF,nq,nep->nconv,pU,&U);
1360:     BVSetActiveColumns(nep->V,0,nq);
1361:     BVMultInPlace(nep->V,U,0,nep->nconv);
1362:     BVSetActiveColumns(nep->V,0,nep->nconv);
1363:     MatDestroy(&U);
1364:     PetscFree(pU);
1365:     /* truncate Schur decomposition and change the state to raw so that
1366:        DSVectors() computes eigenvectors from scratch */
1367:     DSSetDimensions(nep->ds,nep->nconv,0,0,0);
1368:     DSSetState(nep->ds,DS_STATE_RAW);
1369:   }

1371:   /* Map eigenvalues back to the original problem */
1372:   if (!ctx->nshifts) {
1373:     NEPNLEIGSBackTransform((PetscObject)nep,nep->nconv,nep->eigr,nep->eigi);
1374:     PetscFree2(eigr,eigi);
1375:   }
1376:   BVDestroy(&W);
1377:   return(0);
1378: }

1380: static PetscErrorCode NEPNLEIGSSetSingularitiesFunction_NLEIGS(NEP nep,PetscErrorCode (*fun)(NEP,PetscInt*,PetscScalar*,void*),void *ctx)
1381: {
1382:   NEP_NLEIGS *nepctx=(NEP_NLEIGS*)nep->data;

1385:   if (fun) nepctx->computesingularities = fun;
1386:   if (ctx) nepctx->singularitiesctx     = ctx;
1387:   nep->state = NEP_STATE_INITIAL;
1388:   return(0);
1389: }

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

1395:    Logically Collective on NEP

1397:    Input Parameters:
1398: +  nep - the NEP context
1399: .  fun - user function (if NULL then NEP retains any previously set value)
1400: -  ctx - [optional] user-defined context for private data for the function
1401:          (may be NULL, in which case NEP retains any previously set value)

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

1406: +   nep   - the NEP context
1407: .   maxnp - on input number of requested points in the discretization (can be set)
1408: .   xi    - computed values of the discretization
1409: -   ctx   - optional context, as set by NEPNLEIGSSetSingularitiesFunction()

1411:    Notes:
1412:    The user-defined function can set a smaller value of maxnp if necessary.
1413:    It is wrong to return a larger value.

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

1419:    Level: intermediate

1421: .seealso: NEPNLEIGSGetSingularitiesFunction(), NEPSetProblemType()
1422: @*/
1423: PetscErrorCode NEPNLEIGSSetSingularitiesFunction(NEP nep,PetscErrorCode (*fun)(NEP,PetscInt*,PetscScalar*,void*),void *ctx)
1424: {

1429:   PetscTryMethod(nep,"NEPNLEIGSSetSingularitiesFunction_C",(NEP,PetscErrorCode(*)(NEP,PetscInt*,PetscScalar*,void*),void*),(nep,fun,ctx));
1430:   return(0);
1431: }

1433: static PetscErrorCode NEPNLEIGSGetSingularitiesFunction_NLEIGS(NEP nep,PetscErrorCode (**fun)(NEP,PetscInt*,PetscScalar*,void*),void **ctx)
1434: {
1435:   NEP_NLEIGS *nepctx=(NEP_NLEIGS*)nep->data;

1438:   if (fun) *fun = nepctx->computesingularities;
1439:   if (ctx) *ctx = nepctx->singularitiesctx;
1440:   return(0);
1441: }

1443: /*@C
1444:    NEPNLEIGSGetSingularitiesFunction - Returns the Function and optionally the user
1445:    provided context for computing a discretization of the singularity set.

1447:    Not Collective

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

1452:    Output Parameters:
1453: +  fun - location to put the function (or NULL)
1454: -  ctx - location to stash the function context (or NULL)

1456:    Level: advanced

1458: .seealso: NEPNLEIGSSetSingularitiesFunction()
1459: @*/
1460: PetscErrorCode NEPNLEIGSGetSingularitiesFunction(NEP nep,PetscErrorCode (**fun)(NEP,PetscInt*,PetscScalar*,void*),void **ctx)
1461: {

1466:   PetscUseMethod(nep,"NEPNLEIGSGetSingularitiesFunction_C",(NEP,PetscErrorCode(**)(NEP,PetscInt*,PetscScalar*,void*),void**),(nep,fun,ctx));
1467:   return(0);
1468: }

1470: static PetscErrorCode NEPNLEIGSSetRestart_NLEIGS(NEP nep,PetscReal keep)
1471: {
1472:   NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;

1475:   if (keep==PETSC_DEFAULT) ctx->keep = 0.5;
1476:   else {
1477:     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]");
1478:     ctx->keep = keep;
1479:   }
1480:   return(0);
1481: }

1483: /*@
1484:    NEPNLEIGSSetRestart - Sets the restart parameter for the NLEIGS
1485:    method, in particular the proportion of basis vectors that must be kept
1486:    after restart.

1488:    Logically Collective on NEP

1490:    Input Parameters:
1491: +  nep  - the nonlinear eigensolver context
1492: -  keep - the number of vectors to be kept at restart

1494:    Options Database Key:
1495: .  -nep_nleigs_restart - Sets the restart parameter

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

1500:    Level: advanced

1502: .seealso: NEPNLEIGSGetRestart()
1503: @*/
1504: PetscErrorCode NEPNLEIGSSetRestart(NEP nep,PetscReal keep)
1505: {

1511:   PetscTryMethod(nep,"NEPNLEIGSSetRestart_C",(NEP,PetscReal),(nep,keep));
1512:   return(0);
1513: }

1515: static PetscErrorCode NEPNLEIGSGetRestart_NLEIGS(NEP nep,PetscReal *keep)
1516: {
1517:   NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;

1520:   *keep = ctx->keep;
1521:   return(0);
1522: }

1524: /*@
1525:    NEPNLEIGSGetRestart - Gets the restart parameter used in the NLEIGS method.

1527:    Not Collective

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

1532:    Output Parameter:
1533: .  keep - the restart parameter

1535:    Level: advanced

1537: .seealso: NEPNLEIGSSetRestart()
1538: @*/
1539: PetscErrorCode NEPNLEIGSGetRestart(NEP nep,PetscReal *keep)
1540: {

1546:   PetscUseMethod(nep,"NEPNLEIGSGetRestart_C",(NEP,PetscReal*),(nep,keep));
1547:   return(0);
1548: }

1550: static PetscErrorCode NEPNLEIGSSetLocking_NLEIGS(NEP nep,PetscBool lock)
1551: {
1552:   NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;

1555:   ctx->lock = lock;
1556:   return(0);
1557: }

1559: /*@
1560:    NEPNLEIGSSetLocking - Choose between locking and non-locking variants of
1561:    the NLEIGS method.

1563:    Logically Collective on NEP

1565:    Input Parameters:
1566: +  nep  - the nonlinear eigensolver context
1567: -  lock - true if the locking variant must be selected

1569:    Options Database Key:
1570: .  -nep_nleigs_locking - Sets the locking flag

1572:    Notes:
1573:    The default is to lock converged eigenpairs when the method restarts.
1574:    This behaviour can be changed so that all directions are kept in the
1575:    working subspace even if already converged to working accuracy (the
1576:    non-locking variant).

1578:    Level: advanced

1580: .seealso: NEPNLEIGSGetLocking()
1581: @*/
1582: PetscErrorCode NEPNLEIGSSetLocking(NEP nep,PetscBool lock)
1583: {

1589:   PetscTryMethod(nep,"NEPNLEIGSSetLocking_C",(NEP,PetscBool),(nep,lock));
1590:   return(0);
1591: }

1593: static PetscErrorCode NEPNLEIGSGetLocking_NLEIGS(NEP nep,PetscBool *lock)
1594: {
1595:   NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;

1598:   *lock = ctx->lock;
1599:   return(0);
1600: }

1602: /*@
1603:    NEPNLEIGSGetLocking - Gets the locking flag used in the NLEIGS method.

1605:    Not Collective

1607:    Input Parameter:
1608: .  nep - the nonlinear eigensolver context

1610:    Output Parameter:
1611: .  lock - the locking flag

1613:    Level: advanced

1615: .seealso: NEPNLEIGSSetLocking()
1616: @*/
1617: PetscErrorCode NEPNLEIGSGetLocking(NEP nep,PetscBool *lock)
1618: {

1624:   PetscUseMethod(nep,"NEPNLEIGSGetLocking_C",(NEP,PetscBool*),(nep,lock));
1625:   return(0);
1626: }

1628: static PetscErrorCode NEPNLEIGSSetInterpolation_NLEIGS(NEP nep,PetscReal tol,PetscInt degree)
1629: {
1631:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;

1634:   if (tol == PETSC_DEFAULT) {
1635:     ctx->ddtol = PETSC_DEFAULT;
1636:     nep->state = NEP_STATE_INITIAL;
1637:   } else {
1638:     if (tol <= 0.0) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of tol. Must be > 0");
1639:     ctx->ddtol = tol;
1640:   }
1641:   if (degree == PETSC_DEFAULT || degree == PETSC_DECIDE) {
1642:     ctx->ddmaxit = 0;
1643:     if (nep->state) { NEPReset(nep); }
1644:     nep->state = NEP_STATE_INITIAL;
1645:   } else {
1646:     if (degree <= 0) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of degree. Must be > 0");
1647:     if (ctx->ddmaxit != degree) {
1648:       ctx->ddmaxit = degree;
1649:       if (nep->state) { NEPReset(nep); }
1650:       nep->state = NEP_STATE_INITIAL;
1651:     }
1652:   }
1653:   return(0);
1654: }

1656: /*@
1657:    NEPNLEIGSSetInterpolation - Sets the tolerance and maximum degree
1658:    when building the interpolation via divided differences.

1660:    Logically Collective on NEP

1662:    Input Parameters:
1663: +  nep    - the nonlinear eigensolver context
1664: .  tol    - tolerance to stop computing divided differences
1665: -  degree - maximum degree of interpolation

1667:    Options Database Key:
1668: +  -nep_nleigs_interpolation_tol <tol> - Sets the tolerance to stop computing divided differences
1669: -  -nep_nleigs_interpolation_degree <degree> - Sets the maximum degree of interpolation

1671:    Notes:
1672:    Use PETSC_DEFAULT for either argument to assign a reasonably good value.

1674:    Level: advanced

1676: .seealso: NEPNLEIGSGetInterpolation()
1677: @*/
1678: PetscErrorCode NEPNLEIGSSetInterpolation(NEP nep,PetscReal tol,PetscInt degree)
1679: {

1686:   PetscTryMethod(nep,"NEPNLEIGSSetInterpolation_C",(NEP,PetscReal,PetscInt),(nep,tol,degree));
1687:   return(0);
1688: }

1690: static PetscErrorCode NEPNLEIGSGetInterpolation_NLEIGS(NEP nep,PetscReal *tol,PetscInt *degree)
1691: {
1692:   NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;

1695:   if (tol)    *tol    = ctx->ddtol;
1696:   if (degree) *degree = ctx->ddmaxit;
1697:   return(0);
1698: }

1700: /*@
1701:    NEPNLEIGSGetInterpolation - Gets the tolerance and maximum degree
1702:    when building the interpolation via divided differences.

1704:    Not Collective

1706:    Input Parameter:
1707: .  nep - the nonlinear eigensolver context

1709:    Output Parameter:
1710: +  tol    - tolerance to stop computing divided differences
1711: -  degree - maximum degree of interpolation

1713:    Level: advanced

1715: .seealso: NEPNLEIGSSetInterpolation()
1716: @*/
1717: PetscErrorCode NEPNLEIGSGetInterpolation(NEP nep,PetscReal *tol,PetscInt *degree)
1718: {

1723:   PetscTryMethod(nep,"NEPNLEIGSGetInterpolation_C",(NEP,PetscReal*,PetscInt*),(nep,tol,degree));
1724:   return(0);
1725: }

1727: static PetscErrorCode NEPNLEIGSSetRKShifts_NLEIGS(NEP nep,PetscInt ns,PetscScalar *shifts)
1728: {
1730:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
1731:   PetscInt       i;

1734:   if (ns<=0) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"Number of shifts must be positive");
1735:   if (ctx->nshifts) { PetscFree(ctx->shifts); }
1736:   for (i=0;i<ctx->nshiftsw;i++) { KSPDestroy(&ctx->ksp[i]); }
1737:   PetscFree(ctx->ksp);
1738:   ctx->ksp = NULL;
1739:   PetscMalloc1(ns,&ctx->shifts);
1740:   for (i=0;i<ns;i++) ctx->shifts[i] = shifts[i];
1741:   ctx->nshifts = ns;
1742:   nep->state   = NEP_STATE_INITIAL;
1743:   return(0);
1744: }

1746: /*@C
1747:    NEPNLEIGSSetRKShifts - Sets a list of shifts to be used in the Rational
1748:    Krylov method.

1750:    Logically Collective on NEP

1752:    Input Parameters:
1753: +  nep    - the nonlinear eigensolver context
1754: .  ns     - number of shifts
1755: -  shifts - array of scalar values specifying the shifts

1757:    Options Database Key:
1758: .  -nep_nleigs_rk_shifts - Sets the list of shifts

1760:    Notes:
1761:    If only one shift is provided, the built subspace built is equivalent to
1762:    shift-and-invert Krylov-Schur (provided that the absolute convergence
1763:    criterion is used).

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

1770:    Level: advanced

1772: .seealso: NEPNLEIGSGetRKShifts()
1773: @*/
1774: PetscErrorCode NEPNLEIGSSetRKShifts(NEP nep,PetscInt ns,PetscScalar *shifts)
1775: {

1782:   PetscTryMethod(nep,"NEPNLEIGSSetRKShifts_C",(NEP,PetscInt,PetscScalar*),(nep,ns,shifts));
1783:   return(0);
1784: }

1786: static PetscErrorCode NEPNLEIGSGetRKShifts_NLEIGS(NEP nep,PetscInt *ns,PetscScalar **shifts)
1787: {
1789:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
1790:   PetscInt       i;

1793:   *ns = ctx->nshifts;
1794:   if (ctx->nshifts) {
1795:     PetscMalloc1(ctx->nshifts,shifts);
1796:     for (i=0;i<ctx->nshifts;i++) (*shifts)[i] = ctx->shifts[i];
1797:   }
1798:   return(0);
1799: }

1801: /*@C
1802:    NEPNLEIGSGetRKShifts - Gets the list of shifts used in the Rational
1803:    Krylov method.

1805:    Not Collective

1807:    Input Parameter:
1808: .  nep - the nonlinear eigensolver context

1810:    Output Parameter:
1811: +  ns     - number of shifts
1812: -  shifts - array of shifts

1814:    Level: advanced

1816: .seealso: NEPNLEIGSSetRKShifts()
1817: @*/
1818: PetscErrorCode NEPNLEIGSGetRKShifts(NEP nep,PetscInt *ns,PetscScalar **shifts)
1819: {

1826:   PetscTryMethod(nep,"NEPNLEIGSGetRKShifts_C",(NEP,PetscInt*,PetscScalar**),(nep,ns,shifts));
1827:   return(0);
1828: }

1830: static PetscErrorCode NEPNLEIGSGetKSPs_NLEIGS(NEP nep,PetscInt *nsolve,KSP **ksp)
1831: {
1833:   NEP_NLEIGS     *ctx = (NEP_NLEIGS*)nep->data;
1834:   PetscInt       i;
1835:   PC             pc;

1838:   if (!ctx->ksp) {
1839:     NEPNLEIGSSetShifts(nep,&ctx->nshiftsw);
1840:     PetscMalloc1(ctx->nshiftsw,&ctx->ksp);
1841:     for (i=0;i<ctx->nshiftsw;i++) {
1842:       KSPCreate(PetscObjectComm((PetscObject)nep),&ctx->ksp[i]);
1843:       PetscObjectIncrementTabLevel((PetscObject)ctx->ksp[i],(PetscObject)nep,1);
1844:       KSPSetOptionsPrefix(ctx->ksp[i],((PetscObject)nep)->prefix);
1845:       KSPAppendOptionsPrefix(ctx->ksp[i],"nep_nleigs_");
1846:       PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->ksp[i]);
1847:       PetscObjectSetOptions((PetscObject)ctx->ksp[i],((PetscObject)nep)->options);
1848:       KSPSetErrorIfNotConverged(ctx->ksp[i],PETSC_TRUE);
1849:       KSPSetTolerances(ctx->ksp[i],SLEPC_DEFAULT_TOL,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
1850:       KSPGetPC(ctx->ksp[i],&pc);
1851:       KSPSetType(ctx->ksp[i],KSPPREONLY);
1852:       PCSetType(pc,PCLU);
1853:     }
1854:   }
1855:   if (nsolve) *nsolve = ctx->nshiftsw;
1856:   if (ksp)    *ksp    = ctx->ksp;
1857:   return(0);
1858: }

1860: /*@C
1861:    NEPNLEIGSGetKSPs - Retrieve the array of linear solver objects associated with
1862:    the nonlinear eigenvalue solver.

1864:    Not Collective

1866:    Input Parameter:
1867: .  nep - nonlinear eigenvalue solver

1869:    Output Parameters:
1870: +  nsolve - number of returned KSP objects
1871: -  ksp - array of linear solver object

1873:    Notes:
1874:    The number of KSP objects is equal to the number of shifts provided by the user,
1875:    or 1 if the user did not provide shifts.

1877:    Level: advanced

1879: .seealso: NEPNLEIGSSetRKShifts()
1880: @*/
1881: PetscErrorCode NEPNLEIGSGetKSPs(NEP nep,PetscInt *nsolve,KSP **ksp)
1882: {

1887:   PetscUseMethod(nep,"NEPNLEIGSGetKSPs_C",(NEP,PetscInt*,KSP**),(nep,nsolve,ksp));
1888:   return(0);
1889: }

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

1896:   if (fullbasis!=ctx->fullbasis) {
1897:     ctx->fullbasis = fullbasis;
1898:     nep->state     = NEP_STATE_INITIAL;
1899:     nep->useds     = PetscNot(fullbasis);
1900:   }
1901:   return(0);
1902: }

1904: /*@
1905:    NEPNLEIGSSetFullBasis - Choose between TOAR-basis (default) and full-basis
1906:    variants of the NLEIGS method.

1908:    Logically Collective on NEP

1910:    Input Parameters:
1911: +  nep  - the nonlinear eigensolver context
1912: -  fullbasis - true if the full-basis variant must be selected

1914:    Options Database Key:
1915: .  -nep_nleigs_full_basis - Sets the full-basis flag

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

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

1927:    Level: advanced

1929: .seealso: NEPNLEIGSGetFullBasis(), NEPNLEIGSGetEPS(), NEPSetTwoSided(), BVCreateTensor()
1930: @*/
1931: PetscErrorCode NEPNLEIGSSetFullBasis(NEP nep,PetscBool fullbasis)
1932: {

1938:   PetscTryMethod(nep,"NEPNLEIGSSetFullBasis_C",(NEP,PetscBool),(nep,fullbasis));
1939:   return(0);
1940: }

1942: static PetscErrorCode NEPNLEIGSGetFullBasis_NLEIGS(NEP nep,PetscBool *fullbasis)
1943: {
1944:   NEP_NLEIGS *ctx=(NEP_NLEIGS*)nep->data;

1947:   *fullbasis = ctx->fullbasis;
1948:   return(0);
1949: }

1951: /*@
1952:    NEPNLEIGSGetFullBasis - Gets the flag that indicates if NLEIGS is using the
1953:    full-basis variant.

1955:    Not Collective

1957:    Input Parameter:
1958: .  nep - the nonlinear eigensolver context

1960:    Output Parameter:
1961: .  fullbasis - the flag

1963:    Level: advanced

1965: .seealso: NEPNLEIGSSetFullBasis()
1966: @*/
1967: PetscErrorCode NEPNLEIGSGetFullBasis(NEP nep,PetscBool *fullbasis)
1968: {

1974:   PetscUseMethod(nep,"NEPNLEIGSGetFullBasis_C",(NEP,PetscBool*),(nep,fullbasis));
1975:   return(0);
1976: }

1978: #define SHIFTMAX 30

1980: PetscErrorCode NEPSetFromOptions_NLEIGS(PetscOptionItems *PetscOptionsObject,NEP nep)
1981: {
1983:   NEP_NLEIGS     *ctx = (NEP_NLEIGS*)nep->data;
1984:   PetscInt       i=0,k;
1985:   PetscBool      flg1,flg2,b;
1986:   PetscReal      r;
1987:   PetscScalar    array[SHIFTMAX];

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

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

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

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

2001:     NEPNLEIGSGetInterpolation(nep,&r,&i);
2002:     if (!i) i = PETSC_DEFAULT;
2003:     PetscOptionsInt("-nep_nleigs_interpolation_degree","Maximum number of terms for interpolation via divided differences","NEPNLEIGSSetInterpolation",i,&i,&flg1);
2004:     PetscOptionsReal("-nep_nleigs_interpolation_tol","Tolerance for interpolation via divided differences","NEPNLEIGSSetInterpolation",r,&r,&flg2);
2005:     if (flg1 || flg2) { NEPNLEIGSSetInterpolation(nep,r,i); }

2007:     k = SHIFTMAX;
2008:     for (i=0;i<k;i++) array[i] = 0;
2009:     PetscOptionsScalarArray("-nep_nleigs_rk_shifts","Shifts for Rational Krylov","NEPNLEIGSSetRKShifts",array,&k,&flg1);
2010:     if (flg1) { NEPNLEIGSSetRKShifts(nep,k,array); }

2012:   PetscOptionsTail();

2014:   if (!ctx->ksp) { NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp); }
2015:   for (i=0;i<ctx->nshiftsw;i++) {
2016:     KSPSetFromOptions(ctx->ksp[i]);
2017:   }

2019:   if (ctx->fullbasis) {
2020:     if (!ctx->eps) { NEPNLEIGSGetEPS(nep,&ctx->eps); }
2021:     EPSSetFromOptions(ctx->eps);
2022:   }
2023:   return(0);
2024: }

2026: PetscErrorCode NEPView_NLEIGS(NEP nep,PetscViewer viewer)
2027: {
2029:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;
2030:   PetscBool      isascii;
2031:   PetscInt       i;
2032:   char           str[50];

2035:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
2036:   if (isascii) {
2037:     PetscViewerASCIIPrintf(viewer,"  %d%% of basis vectors kept after restart\n",(int)(100*ctx->keep));
2038:     if (ctx->fullbasis) {
2039:       PetscViewerASCIIPrintf(viewer,"  using the full-basis variant\n");
2040:     } else {
2041:       PetscViewerASCIIPrintf(viewer,"  using the %slocking variant\n",ctx->lock?"":"non-");
2042:     }
2043:     PetscViewerASCIIPrintf(viewer,"  divided difference terms: used=%D, max=%D\n",ctx->nmat,ctx->ddmaxit);
2044:     PetscViewerASCIIPrintf(viewer,"  tolerance for divided difference convergence: %g\n",(double)ctx->ddtol);
2045:     if (ctx->nshifts) {
2046:       PetscViewerASCIIPrintf(viewer,"  RK shifts: ");
2047:       PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
2048:       for (i=0;i<ctx->nshifts;i++) {
2049:         SlepcSNPrintfScalar(str,50,ctx->shifts[i],PETSC_FALSE);
2050:         PetscViewerASCIIPrintf(viewer,"%s%s",str,(i<ctx->nshifts-1)?",":"");
2051:       }
2052:       PetscViewerASCIIPrintf(viewer,"\n");
2053:       PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
2054:     }
2055:     if (!ctx->ksp) { NEPNLEIGSGetKSPs(nep,&ctx->nshiftsw,&ctx->ksp); }
2056:     PetscViewerASCIIPushTab(viewer);
2057:     KSPView(ctx->ksp[0],viewer);
2058:     PetscViewerASCIIPopTab(viewer);
2059:     if (ctx->fullbasis) {
2060:       if (!ctx->eps) { NEPNLEIGSGetEPS(nep,&ctx->eps); }
2061:       PetscViewerASCIIPushTab(viewer);
2062:       EPSView(ctx->eps,viewer);
2063:       PetscViewerASCIIPopTab(viewer);
2064:     }
2065:   }
2066:   return(0);
2067: }

2069: PetscErrorCode NEPReset_NLEIGS(NEP nep)
2070: {
2072:   PetscInt       k;
2073:   NEP_NLEIGS     *ctx=(NEP_NLEIGS*)nep->data;

2076:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
2077:     PetscFree(ctx->coeffD);
2078:   } else {
2079:     for (k=0;k<ctx->nmat;k++) { MatDestroy(&ctx->D[k]); }
2080:   }
2081:   PetscFree4(ctx->s,ctx->xi,ctx->beta,ctx->D);
2082:   for (k=0;k<ctx->nshiftsw;k++) { KSPReset(ctx->ksp[k]); }
2083:   if (ctx->vrn) {
2084:     VecDestroy(&ctx->vrn);
2085:   }
2086:   if (ctx->fullbasis) {
2087:     MatDestroy(&ctx->A);
2088:     EPSReset(ctx->eps);
2089:     for (k=0;k<4;k++) { VecDestroy(&ctx->w[k]); }
2090:   }
2091:   return(0);
2092: }

2094: PetscErrorCode NEPDestroy_NLEIGS(NEP nep)
2095: {
2097:   PetscInt       k;
2098:   NEP_NLEIGS     *ctx = (NEP_NLEIGS*)nep->data;

2101:   BVDestroy(&ctx->V);
2102:   for (k=0;k<ctx->nshiftsw;k++) { KSPDestroy(&ctx->ksp[k]); }
2103:   PetscFree(ctx->ksp);
2104:   if (ctx->nshifts) { PetscFree(ctx->shifts); }
2105:   if (ctx->fullbasis) { EPSDestroy(&ctx->eps); }
2106:   PetscFree(nep->data);
2107:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetSingularitiesFunction_C",NULL);
2108:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetSingularitiesFunction_C",NULL);
2109:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRestart_C",NULL);
2110:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRestart_C",NULL);
2111:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetLocking_C",NULL);
2112:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetLocking_C",NULL);
2113:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetInterpolation_C",NULL);
2114:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetInterpolation_C",NULL);
2115:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRKShifts_C",NULL);
2116:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRKShifts_C",NULL);
2117:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetKSPs_C",NULL);
2118:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetFullBasis_C",NULL);
2119:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetFullBasis_C",NULL);
2120:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetEPS_C",NULL);
2121:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetEPS_C",NULL);
2122:   return(0);
2123: }

2125: SLEPC_EXTERN PetscErrorCode NEPCreate_NLEIGS(NEP nep)
2126: {
2128:   NEP_NLEIGS     *ctx;

2131:   PetscNewLog(nep,&ctx);
2132:   nep->data  = (void*)ctx;
2133:   ctx->lock  = PETSC_TRUE;
2134:   ctx->ddtol = PETSC_DEFAULT;

2136:   nep->useds = PETSC_TRUE;
2137:   nep->hasts = PETSC_TRUE;

2139:   nep->ops->setup          = NEPSetUp_NLEIGS;
2140:   nep->ops->setfromoptions = NEPSetFromOptions_NLEIGS;
2141:   nep->ops->view           = NEPView_NLEIGS;
2142:   nep->ops->destroy        = NEPDestroy_NLEIGS;
2143:   nep->ops->reset          = NEPReset_NLEIGS;

2145:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetSingularitiesFunction_C",NEPNLEIGSSetSingularitiesFunction_NLEIGS);
2146:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetSingularitiesFunction_C",NEPNLEIGSGetSingularitiesFunction_NLEIGS);
2147:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRestart_C",NEPNLEIGSSetRestart_NLEIGS);
2148:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRestart_C",NEPNLEIGSGetRestart_NLEIGS);
2149:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetLocking_C",NEPNLEIGSSetLocking_NLEIGS);
2150:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetLocking_C",NEPNLEIGSGetLocking_NLEIGS);
2151:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetInterpolation_C",NEPNLEIGSSetInterpolation_NLEIGS);
2152:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetInterpolation_C",NEPNLEIGSGetInterpolation_NLEIGS);
2153:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetRKShifts_C",NEPNLEIGSSetRKShifts_NLEIGS);
2154:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetRKShifts_C",NEPNLEIGSGetRKShifts_NLEIGS);
2155:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetKSPs_C",NEPNLEIGSGetKSPs_NLEIGS);
2156:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetFullBasis_C",NEPNLEIGSSetFullBasis_NLEIGS);
2157:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetFullBasis_C",NEPNLEIGSGetFullBasis_NLEIGS);
2158:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSSetEPS_C",NEPNLEIGSSetEPS_NLEIGS);
2159:   PetscObjectComposeFunction((PetscObject)nep,"NEPNLEIGSGetEPS_C",NEPNLEIGSGetEPS_NLEIGS);
2160:   return(0);
2161: }