Actual source code: slp.c
slepc-3.11.2 2019-07-30
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: "slp"
13: Method: Succesive linear problems
15: Algorithm:
17: Newton-type iteration based on first order Taylor approximation.
19: References:
21: [1] A. Ruhe, "Algorithms for the nonlinear eigenvalue problem", SIAM J.
22: Numer. Anal. 10(4):674-689, 1973.
23: */
25: #include <slepc/private/nepimpl.h> /*I "slepcnep.h" I*/
26: #include <../src/nep/impls/nepdefl.h>
28: typedef struct {
29: EPS eps; /* linear eigensolver for T*z = mu*Tp*z */
30: KSP ksp;
31: } NEP_SLP;
33: typedef struct {
34: NEP_EXT_OP extop;
35: Vec w;
36: } NEP_SLP_EPS_MSHELL;
38: PetscErrorCode NEPSetUp_SLP(NEP nep)
39: {
41: NEP_SLP *ctx = (NEP_SLP*)nep->data;
42: PetscBool istrivial,flg;
43: ST st;
46: if (nep->ncv) { PetscInfo(nep,"Setting ncv = nev, ignoring user-provided value\n"); }
47: nep->ncv = nep->nev;
48: if (nep->mpd) { PetscInfo(nep,"Setting mpd = nev, ignoring user-provided value\n"); }
49: nep->mpd = nep->nev;
50: if (nep->ncv>nep->nev+nep->mpd) SETERRQ(PetscObjectComm((PetscObject)nep),1,"The value of ncv must not be larger than nev+mpd");
51: if (!nep->max_it) nep->max_it = PetscMax(5000,2*nep->n/nep->ncv);
52: if (nep->which && nep->which!=NEP_TARGET_MAGNITUDE) SETERRQ(PetscObjectComm((PetscObject)nep),1,"Wrong value of which");
54: RGIsTrivial(nep->rg,&istrivial);
55: if (!istrivial) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"This solver does not support region filtering");
57: if (!ctx->eps) { NEPSLPGetEPS(nep,&ctx->eps); }
58: EPSGetST(ctx->eps,&st);
59: PetscObjectTypeCompareAny((PetscObject)st,&flg,STSINVERT,STCAYLEY,"");
60: if (flg) SETERRQ(PetscObjectComm((PetscObject)nep),1,"SLP does not support spectral transformation");
61: EPSSetDimensions(ctx->eps,1,PETSC_DECIDE,PETSC_DECIDE);
62: EPSSetWhichEigenpairs(ctx->eps,EPS_LARGEST_MAGNITUDE);
63: EPSSetTolerances(ctx->eps,nep->tol==PETSC_DEFAULT?SLEPC_DEFAULT_TOL/10.0:nep->tol/10.0,nep->max_it?nep->max_it:PETSC_DEFAULT);
65: NEPAllocateSolution(nep,0);
66: return(0);
67: }
69: PetscErrorCode NEPSLPEPSMatShell_MatMult(Mat M,Vec x,Vec y)
70: {
71: PetscErrorCode ierr;
72: NEP_SLP_EPS_MSHELL *ctx;
75: MatShellGetContext(M,(void**)&ctx);
76: MatMult(ctx->extop->MJ,x,ctx->w);
77: NEPDeflationFunctionSolve(ctx->extop,ctx->w,y);
78: return(0);
79: }
81: PetscErrorCode NEPSLPEPSMatShell_Destroy(Mat M)
82: {
83: PetscErrorCode ierr;
84: NEP_SLP_EPS_MSHELL *ctx;
87: MatShellGetContext(M,(void**)&ctx);
88: VecDestroy(&ctx->w);
89: PetscFree(ctx);
90: return(0);
91: }
93: PetscErrorCode NEPSLPEPSMatShell_CreateVecs(Mat M,Vec *left,Vec *right)
94: {
95: PetscErrorCode ierr;
96: NEP_SLP_EPS_MSHELL *ctx;
99: MatShellGetContext(M,(void**)&ctx);
100: if (right) {
101: VecDuplicate(ctx->w,right);
102: }
103: if (left) {
104: VecDuplicate(ctx->w,left);
105: }
106: return(0);
107: }
109: static PetscErrorCode NEPSLPSetUpLinearEP(NEP nep,NEP_EXT_OP extop,PetscScalar lambda,Vec u,PetscBool ini)
110: {
111: PetscErrorCode ierr;
112: NEP_SLP *slpctx = (NEP_SLP*)nep->data;
113: Mat Mshell;
114: PetscInt nloc,mloc;
115: NEP_SLP_EPS_MSHELL *shellctx;
118: if (ini) {
119: /* Create mat shell */
120: PetscNew(&shellctx);
121: shellctx->extop = extop;
122: NEPDeflationCreateVec(extop,&shellctx->w);
123: MatGetLocalSize(nep->function,&mloc,&nloc);
124: nloc += extop->szd; mloc += extop->szd;
125: MatCreateShell(PetscObjectComm((PetscObject)nep),nloc,mloc,PETSC_DETERMINE,PETSC_DETERMINE,shellctx,&Mshell);
126: MatShellSetOperation(Mshell,MATOP_MULT,(void(*)(void))NEPSLPEPSMatShell_MatMult);
127: MatShellSetOperation(Mshell,MATOP_DESTROY,(void(*)(void))NEPSLPEPSMatShell_Destroy);
128: MatShellSetOperation(Mshell,MATOP_CREATE_VECS,(void(*)(void))NEPSLPEPSMatShell_CreateVecs);
129: EPSSetOperators(slpctx->eps,Mshell,NULL);
130: MatDestroy(&Mshell);
131: }
132: NEPDeflationSolveSetUp(extop,lambda);
133: NEPDeflationComputeJacobian(extop,lambda,NULL);
134: EPSSetInitialSpace(slpctx->eps,1,&u);
135: return(0);
136: }
138: PetscErrorCode NEPSolve_SLP(NEP nep)
139: {
141: NEP_SLP *ctx = (NEP_SLP*)nep->data;
142: Mat F,H;
143: Vec uu,u,r;
144: PetscScalar sigma,lambda,mu,im,*Hp,*Ap;
145: PetscReal resnorm;
146: PetscInt nconv,ldh,ldds,i,j;
147: PetscBool skip=PETSC_FALSE;
148: NEP_EXT_OP extop=NULL; /* Extended operator for deflation */
151: /* get initial approximation of eigenvalue and eigenvector */
152: NEPGetDefaultShift(nep,&sigma);
153: if (!nep->nini) {
154: BVSetRandomColumn(nep->V,0);
155: }
156: lambda = sigma;
157: if (!ctx->ksp) { NEPSLPGetKSP(nep,&ctx->ksp); }
158: NEPDeflationInitialize(nep,nep->V,ctx->ksp,PETSC_TRUE,nep->nev,&extop);
159: NEPDeflationCreateVec(extop,&u);
160: VecDuplicate(u,&r);
161: BVGetColumn(nep->V,0,&uu);
162: NEPDeflationCopyToExtendedVec(extop,uu,NULL,u,PETSC_FALSE);
163: BVRestoreColumn(nep->V,0,&uu);
165: /* Restart loop */
166: while (nep->reason == NEP_CONVERGED_ITERATING) {
167: nep->its++;
169: /* form residual, r = T(lambda)*u (used in convergence test only) */
170: NEPDeflationComputeFunction(extop,lambda,&F);
171: MatMult(F,u,r);
173: /* convergence test */
174: VecNorm(r,NORM_2,&resnorm);
175: (*nep->converged)(nep,lambda,0,resnorm,&nep->errest[nep->nconv],nep->convergedctx);
176: nep->eigr[nep->nconv] = lambda;
177: if (nep->errest[nep->nconv]<=nep->tol) {
178: nep->nconv = nep->nconv + 1;
179: skip = PETSC_TRUE;
180: NEPDeflationLocking(extop,u,lambda);
181: }
182: (*nep->stopping)(nep,nep->its,nep->max_it,nep->nconv,nep->nev,&nep->reason,nep->stoppingctx);
183: if (!skip || nep->reason>0) {
184: NEPMonitor(nep,nep->its,nep->nconv,nep->eigr,nep->eigi,nep->errest,(nep->reason>0)?nep->nconv:nep->nconv+1);
185: }
187: if (nep->reason == NEP_CONVERGED_ITERATING) {
188: if (!skip) {
189: /* evaluate T(lambda) and T'(lambda) */
190: NEPSLPSetUpLinearEP(nep,extop,lambda,u,nep->its==1?PETSC_TRUE:PETSC_FALSE);
191: /* compute new eigenvalue correction mu and eigenvector approximation u */
192: EPSSolve(ctx->eps);
193: EPSGetConverged(ctx->eps,&nconv);
194: if (!nconv) {
195: PetscInfo1(nep,"iter=%D, inner iteration failed, stopping solve\n",nep->its);
196: nep->reason = NEP_DIVERGED_LINEAR_SOLVE;
197: break;
198: }
199: EPSGetEigenpair(ctx->eps,0,&mu,&im,u,NULL);
200: mu = 1.0/mu;
201: if (PetscAbsScalar(im)>PETSC_MACHINE_EPSILON) SETERRQ(PetscObjectComm((PetscObject)nep),1,"Complex eigenvalue approximation - not implemented in real scalars");
202: } else {
203: nep->its--; /* do not count this as a full iteration */
204: /* use second eigenpair computed in previous iteration */
205: EPSGetConverged(ctx->eps,&nconv);
206: if (nconv>=2) {
207: EPSGetEigenpair(ctx->eps,1,&mu,&im,u,NULL);
208: mu = 1.0/mu;
209: } else {
210: NEPDeflationSetRandomVec(extop,u);
211: mu = lambda-sigma;
212: }
213: skip = PETSC_FALSE;
214: }
215: /* correct eigenvalue */
216: lambda = lambda - mu;
217: }
218: }
219: NEPDeflationGetInvariantPair(extop,NULL,&H);
220: MatGetSize(H,NULL,&ldh);
221: DSSetType(nep->ds,DSNHEP);
222: DSAllocate(nep->ds,PetscMax(nep->nconv,1));
223: DSGetLeadingDimension(nep->ds,&ldds);
224: MatDenseGetArray(H,&Hp);
225: DSGetArray(nep->ds,DS_MAT_A,&Ap);
226: for (j=0;j<nep->nconv;j++)
227: for (i=0;i<nep->nconv;i++) Ap[j*ldds+i] = Hp[j*ldh+i];
228: DSRestoreArray(nep->ds,DS_MAT_A,&Ap);
229: MatDenseRestoreArray(H,&Hp);
230: MatDestroy(&H);
231: DSSetDimensions(nep->ds,nep->nconv,0,0,nep->nconv);
232: DSSolve(nep->ds,nep->eigr,nep->eigi);
233: NEPDeflationReset(extop);
234: VecDestroy(&u);
235: VecDestroy(&r);
236: return(0);
237: }
239: PetscErrorCode NEPSetFromOptions_SLP(PetscOptionItems *PetscOptionsObject,NEP nep)
240: {
242: NEP_SLP *ctx = (NEP_SLP*)nep->data;
245: PetscOptionsHead(PetscOptionsObject,"NEP SLP Options");
246: PetscOptionsTail();
248: if (!ctx->eps) { NEPSLPGetEPS(nep,&ctx->eps); }
249: EPSSetFromOptions(ctx->eps);
250: if (!ctx->ksp) { NEPSLPGetKSP(nep,&ctx->ksp); }
251: KSPSetFromOptions(ctx->ksp);
252: return(0);
253: }
255: static PetscErrorCode NEPSLPSetEPS_SLP(NEP nep,EPS eps)
256: {
258: NEP_SLP *ctx = (NEP_SLP*)nep->data;
261: PetscObjectReference((PetscObject)eps);
262: EPSDestroy(&ctx->eps);
263: ctx->eps = eps;
264: PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->eps);
265: nep->state = NEP_STATE_INITIAL;
266: return(0);
267: }
269: /*@
270: NEPSLPSetEPS - Associate a linear eigensolver object (EPS) to the
271: nonlinear eigenvalue solver.
273: Collective on NEP
275: Input Parameters:
276: + nep - nonlinear eigenvalue solver
277: - eps - the eigensolver object
279: Level: advanced
281: .seealso: NEPSLPGetEPS()
282: @*/
283: PetscErrorCode NEPSLPSetEPS(NEP nep,EPS eps)
284: {
291: PetscTryMethod(nep,"NEPSLPSetEPS_C",(NEP,EPS),(nep,eps));
292: return(0);
293: }
295: static PetscErrorCode NEPSLPGetEPS_SLP(NEP nep,EPS *eps)
296: {
298: NEP_SLP *ctx = (NEP_SLP*)nep->data;
301: if (!ctx->eps) {
302: EPSCreate(PetscObjectComm((PetscObject)nep),&ctx->eps);
303: PetscObjectIncrementTabLevel((PetscObject)ctx->eps,(PetscObject)nep,1);
304: EPSSetOptionsPrefix(ctx->eps,((PetscObject)nep)->prefix);
305: EPSAppendOptionsPrefix(ctx->eps,"nep_slp_");
306: PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->eps);
307: PetscObjectSetOptions((PetscObject)ctx->eps,((PetscObject)nep)->options);
308: }
309: *eps = ctx->eps;
310: return(0);
311: }
313: /*@
314: NEPSLPGetEPS - Retrieve the linear eigensolver object (EPS) associated
315: to the nonlinear eigenvalue solver.
317: Not Collective
319: Input Parameter:
320: . nep - nonlinear eigenvalue solver
322: Output Parameter:
323: . eps - the eigensolver object
325: Level: advanced
327: .seealso: NEPSLPSetEPS()
328: @*/
329: PetscErrorCode NEPSLPGetEPS(NEP nep,EPS *eps)
330: {
336: PetscUseMethod(nep,"NEPSLPGetEPS_C",(NEP,EPS*),(nep,eps));
337: return(0);
338: }
340: static PetscErrorCode NEPSLPSetKSP_SLP(NEP nep,KSP ksp)
341: {
343: NEP_SLP *ctx = (NEP_SLP*)nep->data;
346: PetscObjectReference((PetscObject)ksp);
347: KSPDestroy(&ctx->ksp);
348: ctx->ksp = ksp;
349: PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->ksp);
350: nep->state = NEP_STATE_INITIAL;
351: return(0);
352: }
354: /*@
355: NEPSLPSetKSP - Associate a linear solver object (KSP) to the nonlinear
356: eigenvalue solver.
358: Collective on NEP
360: Input Parameters:
361: + nep - eigenvalue solver
362: - ksp - the linear solver object
364: Level: advanced
366: .seealso: NEPSLPGetKSP()
367: @*/
368: PetscErrorCode NEPSLPSetKSP(NEP nep,KSP ksp)
369: {
376: PetscTryMethod(nep,"NEPSLPSetKSP_C",(NEP,KSP),(nep,ksp));
377: return(0);
378: }
380: static PetscErrorCode NEPSLPGetKSP_SLP(NEP nep,KSP *ksp)
381: {
383: NEP_SLP *ctx = (NEP_SLP*)nep->data;
386: if (!ctx->ksp) {
387: KSPCreate(PetscObjectComm((PetscObject)nep),&ctx->ksp);
388: PetscObjectIncrementTabLevel((PetscObject)ctx->ksp,(PetscObject)nep,1);
389: KSPSetOptionsPrefix(ctx->ksp,((PetscObject)nep)->prefix);
390: KSPAppendOptionsPrefix(ctx->ksp,"nep_slp_");
391: PetscLogObjectParent((PetscObject)nep,(PetscObject)ctx->ksp);
392: KSPSetErrorIfNotConverged(ctx->ksp,PETSC_TRUE);
393: KSPSetTolerances(ctx->ksp,SLEPC_DEFAULT_TOL,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
394: }
395: *ksp = ctx->ksp;
396: return(0);
397: }
399: /*@
400: NEPSLPGetKSP - Retrieve the linear solver object (KSP) associated with
401: the nonlinear eigenvalue solver.
403: Not Collective
405: Input Parameter:
406: . nep - nonlinear eigenvalue solver
408: Output Parameter:
409: . ksp - the linear solver object
411: Level: advanced
413: .seealso: NEPSLPSetKSP()
414: @*/
415: PetscErrorCode NEPSLPGetKSP(NEP nep,KSP *ksp)
416: {
422: PetscUseMethod(nep,"NEPSLPGetKSP_C",(NEP,KSP*),(nep,ksp));
423: return(0);
424: }
426: PetscErrorCode NEPView_SLP(NEP nep,PetscViewer viewer)
427: {
429: NEP_SLP *ctx = (NEP_SLP*)nep->data;
430: PetscBool isascii;
433: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
434: if (isascii) {
435: if (!ctx->eps) { NEPSLPGetEPS(nep,&ctx->eps); }
436: PetscViewerASCIIPushTab(viewer);
437: EPSView(ctx->eps,viewer);
438: PetscViewerASCIIPopTab(viewer);
439: }
440: return(0);
441: }
443: PetscErrorCode NEPReset_SLP(NEP nep)
444: {
446: NEP_SLP *ctx = (NEP_SLP*)nep->data;
449: EPSReset(ctx->eps);
450: KSPReset(ctx->ksp);
451: return(0);
452: }
454: PetscErrorCode NEPDestroy_SLP(NEP nep)
455: {
457: NEP_SLP *ctx = (NEP_SLP*)nep->data;
460: KSPDestroy(&ctx->ksp);
461: EPSDestroy(&ctx->eps);
462: PetscFree(nep->data);
463: PetscObjectComposeFunction((PetscObject)nep,"NEPSLPSetEPS_C",NULL);
464: PetscObjectComposeFunction((PetscObject)nep,"NEPSLPGetEPS_C",NULL);
465: PetscObjectComposeFunction((PetscObject)nep,"NEPSLPSetKSP_C",NULL);
466: PetscObjectComposeFunction((PetscObject)nep,"NEPSLPGetKSP_C",NULL);
467: return(0);
468: }
470: SLEPC_EXTERN PetscErrorCode NEPCreate_SLP(NEP nep)
471: {
473: NEP_SLP *ctx;
476: PetscNewLog(nep,&ctx);
477: nep->data = (void*)ctx;
479: nep->useds = PETSC_TRUE;
481: nep->ops->solve = NEPSolve_SLP;
482: nep->ops->setup = NEPSetUp_SLP;
483: nep->ops->setfromoptions = NEPSetFromOptions_SLP;
484: nep->ops->reset = NEPReset_SLP;
485: nep->ops->destroy = NEPDestroy_SLP;
486: nep->ops->view = NEPView_SLP;
487: nep->ops->computevectors = NEPComputeVectors_Schur;
489: PetscObjectComposeFunction((PetscObject)nep,"NEPSLPSetEPS_C",NEPSLPSetEPS_SLP);
490: PetscObjectComposeFunction((PetscObject)nep,"NEPSLPGetEPS_C",NEPSLPGetEPS_SLP);
491: PetscObjectComposeFunction((PetscObject)nep,"NEPSLPSetKSP_C",NEPSLPSetKSP_SLP);
492: PetscObjectComposeFunction((PetscObject)nep,"NEPSLPGetKSP_C",NEPSLPGetKSP_SLP);
493: return(0);
494: }