Actual source code: qarnoldi.c

slepc-main 2024-12-17
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain

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

 13:    Method: Q-Arnoldi

 15:    Algorithm:

 17:        Quadratic Arnoldi with Krylov-Schur type restart.

 19:    References:

 21:        [1] K. Meerbergen, "The Quadratic Arnoldi method for the solution
 22:            of the quadratic eigenvalue problem", SIAM J. Matrix Anal.
 23:            Appl. 30(4):1462-1482, 2008.
 24: */

 26: #include <slepc/private/pepimpl.h>
 27: #include <petscblaslapack.h>

 29: typedef struct {
 30:   PetscReal keep;         /* restart parameter */
 31:   PetscBool lock;         /* locking/non-locking variant */
 32: } PEP_QARNOLDI;

 34: static PetscErrorCode PEPSetUp_QArnoldi(PEP pep)
 35: {
 36:   PEP_QARNOLDI   *ctx = (PEP_QARNOLDI*)pep->data;
 37:   PetscBool      flg;

 39:   PetscFunctionBegin;
 40:   PEPCheckQuadratic(pep);
 41:   PEPCheckShiftSinvert(pep);
 42:   PetscCall(PEPSetDimensions_Default(pep,pep->nev,&pep->ncv,&pep->mpd));
 43:   PetscCheck(ctx->lock || pep->mpd>=pep->ncv,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Should not use mpd parameter in non-locking variant");
 44:   if (pep->max_it==PETSC_DETERMINE) pep->max_it = PetscMax(100,4*pep->n/pep->ncv);
 45:   if (!pep->which) PetscCall(PEPSetWhichEigenpairs_Default(pep));
 46:   PetscCheck(pep->which!=PEP_ALL,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"This solver does not support computing all eigenvalues");

 48:   PetscCall(STGetTransform(pep->st,&flg));
 49:   PetscCheck(flg,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Solver requires the ST transformation flag set, see STSetTransform()");

 51:   /* set default extraction */
 52:   if (!pep->extract) pep->extract = PEP_EXTRACT_NONE;
 53:   PEPCheckUnsupported(pep,PEP_FEATURE_NONMONOMIAL | PEP_FEATURE_EXTRACT);

 55:   if (!ctx->keep) ctx->keep = 0.5;

 57:   PetscCall(PEPAllocateSolution(pep,0));
 58:   PetscCall(PEPSetWorkVecs(pep,4));

 60:   PetscCall(DSSetType(pep->ds,DSNHEP));
 61:   PetscCall(DSSetExtraRow(pep->ds,PETSC_TRUE));
 62:   PetscCall(DSAllocate(pep->ds,pep->ncv+1));
 63:   PetscFunctionReturn(PETSC_SUCCESS);
 64: }

 66: static PetscErrorCode PEPExtractVectors_QArnoldi(PEP pep)
 67: {
 68:   PetscInt       k=pep->nconv;
 69:   Mat            X,X0;

 71:   PetscFunctionBegin;
 72:   if (pep->nconv==0) PetscFunctionReturn(PETSC_SUCCESS);
 73:   PetscCall(DSVectors(pep->ds,DS_MAT_X,NULL,NULL));

 75:   /* update vectors V = V*X */
 76:   PetscCall(DSGetMat(pep->ds,DS_MAT_X,&X));
 77:   PetscCall(MatDenseGetSubMatrix(X,0,k,0,k,&X0));
 78:   PetscCall(BVMultInPlace(pep->V,X0,0,k));
 79:   PetscCall(MatDenseRestoreSubMatrix(X,&X0));
 80:   PetscCall(DSRestoreMat(pep->ds,DS_MAT_X,&X));
 81:   PetscFunctionReturn(PETSC_SUCCESS);
 82: }

 84: /*
 85:   Compute a step of Classical Gram-Schmidt orthogonalization
 86: */
 87: static PetscErrorCode PEPQArnoldiCGS(PEP pep,PetscScalar *H,PetscBLASInt ldh,PetscScalar *h,PetscBLASInt j,BV V,Vec t,Vec v,Vec w,PetscReal *onorm,PetscReal *norm,PetscScalar *work)
 88: {
 89:   PetscBLASInt   ione = 1,j_1 = j+1;
 90:   PetscReal      x,y;
 91:   PetscScalar    dot,one = 1.0,zero = 0.0;

 93:   PetscFunctionBegin;
 94:   /* compute norm of v and w */
 95:   if (onorm) {
 96:     PetscCall(VecNorm(v,NORM_2,&x));
 97:     PetscCall(VecNorm(w,NORM_2,&y));
 98:     *onorm = SlepcAbs(x,y);
 99:   }

101:   /* orthogonalize: compute h */
102:   PetscCall(BVDotVec(V,v,h));
103:   PetscCall(BVDotVec(V,w,work));
104:   if (j>0) PetscCallBLAS("BLASgemv",BLASgemv_("C",&j_1,&j,&one,H,&ldh,work,&ione,&one,h,&ione));
105:   PetscCall(VecDot(w,t,&dot));
106:   h[j] += dot;

108:   /* orthogonalize: update v and w */
109:   PetscCall(BVMultVec(V,-1.0,1.0,v,h));
110:   if (j>0) {
111:     PetscCallBLAS("BLASgemv",BLASgemv_("N",&j_1,&j,&one,H,&ldh,h,&ione,&zero,work,&ione));
112:     PetscCall(BVMultVec(V,-1.0,1.0,w,work));
113:   }
114:   PetscCall(VecAXPY(w,-h[j],t));

116:   /* compute norm of v and w */
117:   if (norm) {
118:     PetscCall(VecNorm(v,NORM_2,&x));
119:     PetscCall(VecNorm(w,NORM_2,&y));
120:     *norm = SlepcAbs(x,y);
121:   }
122:   PetscFunctionReturn(PETSC_SUCCESS);
123: }

125: /*
126:   Compute a run of Q-Arnoldi iterations
127: */
128: static PetscErrorCode PEPQArnoldi(PEP pep,Mat A,PetscInt k,PetscInt *M,Vec v,Vec w,PetscReal *beta,PetscBool *breakdown,PetscScalar *work)
129: {
130:   PetscInt           i,j,l,m = *M,ldh;
131:   PetscBLASInt       jj,ldhh;
132:   Vec                t = pep->work[2],u = pep->work[3];
133:   BVOrthogRefineType refinement;
134:   PetscReal          norm=0.0,onorm,eta;
135:   PetscScalar        *H,*c = work + m;

137:   PetscFunctionBegin;
138:   *beta = 0.0;
139:   PetscCall(MatDenseGetArray(A,&H));
140:   PetscCall(MatDenseGetLDA(A,&ldh));
141:   PetscCall(BVGetOrthogonalization(pep->V,NULL,&refinement,&eta,NULL));
142:   PetscCall(BVInsertVec(pep->V,k,v));
143:   for (j=k;j<m;j++) {
144:     /* apply operator */
145:     PetscCall(VecCopy(w,t));
146:     if (pep->Dr) PetscCall(VecPointwiseMult(v,v,pep->Dr));
147:     PetscCall(STMatMult(pep->st,0,v,u));
148:     PetscCall(VecCopy(t,v));
149:     if (pep->Dr) PetscCall(VecPointwiseMult(t,t,pep->Dr));
150:     PetscCall(STMatMult(pep->st,1,t,w));
151:     PetscCall(VecAXPY(u,pep->sfactor,w));
152:     PetscCall(STMatSolve(pep->st,u,w));
153:     PetscCall(VecScale(w,-1.0/(pep->sfactor*pep->sfactor)));
154:     if (pep->Dr) PetscCall(VecPointwiseDivide(w,w,pep->Dr));
155:     PetscCall(VecCopy(v,t));
156:     PetscCall(BVSetActiveColumns(pep->V,0,j+1));

158:     /* orthogonalize */
159:     PetscCall(PetscBLASIntCast(j,&jj));
160:     PetscCall(PetscBLASIntCast(ldh,&ldhh));
161:     switch (refinement) {
162:       case BV_ORTHOG_REFINE_NEVER:
163:         PetscCall(PEPQArnoldiCGS(pep,H,ldhh,H+ldh*j,jj,pep->V,t,v,w,NULL,&norm,work));
164:         *breakdown = PETSC_FALSE;
165:         break;
166:       case BV_ORTHOG_REFINE_ALWAYS:
167:         PetscCall(PEPQArnoldiCGS(pep,H,ldhh,H+ldh*j,jj,pep->V,t,v,w,NULL,NULL,work));
168:         PetscCall(PEPQArnoldiCGS(pep,H,ldhh,c,jj,pep->V,t,v,w,&onorm,&norm,work));
169:         for (i=0;i<=j;i++) H[ldh*j+i] += c[i];
170:         if (norm < eta * onorm) *breakdown = PETSC_TRUE;
171:         else *breakdown = PETSC_FALSE;
172:         break;
173:       case BV_ORTHOG_REFINE_IFNEEDED:
174:         PetscCall(PEPQArnoldiCGS(pep,H,ldhh,H+ldh*j,jj,pep->V,t,v,w,&onorm,&norm,work));
175:         /* ||q|| < eta ||h|| */
176:         l = 1;
177:         while (l<3 && norm < eta * onorm) {
178:           l++;
179:           onorm = norm;
180:           PetscCall(PEPQArnoldiCGS(pep,H,ldhh,c,jj,pep->V,t,v,w,NULL,&norm,work));
181:           for (i=0;i<=j;i++) H[ldh*j+i] += c[i];
182:         }
183:         if (norm < eta * onorm) *breakdown = PETSC_TRUE;
184:         else *breakdown = PETSC_FALSE;
185:         break;
186:     }
187:     PetscCall(VecScale(v,1.0/norm));
188:     PetscCall(VecScale(w,1.0/norm));

190:     H[j+1+ldh*j] = norm;
191:     if (j<m-1) PetscCall(BVInsertVec(pep->V,j+1,v));
192:   }
193:   *beta = norm;
194:   PetscCall(MatDenseRestoreArray(A,&H));
195:   PetscFunctionReturn(PETSC_SUCCESS);
196: }

198: static PetscErrorCode PEPSolve_QArnoldi(PEP pep)
199: {
200:   PEP_QARNOLDI   *ctx = (PEP_QARNOLDI*)pep->data;
201:   PetscInt       j,k,l,lwork,nv,nconv;
202:   Vec            v=pep->work[0],w=pep->work[1];
203:   Mat            Q,S;
204:   PetscScalar    *work;
205:   PetscReal      beta,norm,x,y;
206:   PetscBool      breakdown=PETSC_FALSE,sinv;

208:   PetscFunctionBegin;
209:   lwork = 7*pep->ncv;
210:   PetscCall(PetscMalloc1(lwork,&work));
211:   PetscCall(PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinv));
212:   PetscCall(RGPushScale(pep->rg,sinv?pep->sfactor:1.0/pep->sfactor));
213:   PetscCall(STScaleShift(pep->st,sinv?pep->sfactor:1.0/pep->sfactor));

215:   /* Get the starting Arnoldi vector */
216:   for (j=0;j<2;j++) {
217:     if (j>=pep->nini) PetscCall(BVSetRandomColumn(pep->V,j));
218:   }
219:   PetscCall(BVCopyVec(pep->V,0,v));
220:   PetscCall(BVCopyVec(pep->V,1,w));
221:   PetscCall(VecNorm(v,NORM_2,&x));
222:   PetscCall(VecNorm(w,NORM_2,&y));
223:   norm = SlepcAbs(x,y);
224:   PetscCall(VecScale(v,1.0/norm));
225:   PetscCall(VecScale(w,1.0/norm));

227:   /* clean projected matrix (including the extra-arrow) */
228:   PetscCall(DSSetDimensions(pep->ds,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_DETERMINE));
229:   PetscCall(DSGetMat(pep->ds,DS_MAT_A,&S));
230:   PetscCall(MatZeroEntries(S));
231:   PetscCall(DSRestoreMat(pep->ds,DS_MAT_A,&S));

233:    /* Restart loop */
234:   l = 0;
235:   while (pep->reason == PEP_CONVERGED_ITERATING) {
236:     pep->its++;

238:     /* Compute an nv-step Arnoldi factorization */
239:     nv = PetscMin(pep->nconv+pep->mpd,pep->ncv);
240:     PetscCall(DSGetMat(pep->ds,DS_MAT_A,&S));
241:     PetscCall(PEPQArnoldi(pep,S,pep->nconv+l,&nv,v,w,&beta,&breakdown,work));
242:     PetscCall(DSRestoreMat(pep->ds,DS_MAT_A,&S));
243:     PetscCall(DSSetDimensions(pep->ds,nv,pep->nconv,pep->nconv+l));
244:     PetscCall(DSSetState(pep->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
245:     PetscCall(BVSetActiveColumns(pep->V,pep->nconv,nv));

247:     /* Solve projected problem */
248:     PetscCall(DSSolve(pep->ds,pep->eigr,pep->eigi));
249:     PetscCall(DSSort(pep->ds,pep->eigr,pep->eigi,NULL,NULL,NULL));
250:     PetscCall(DSUpdateExtraRow(pep->ds));
251:     PetscCall(DSSynchronize(pep->ds,pep->eigr,pep->eigi));

253:     /* Check convergence */
254:     PetscCall(PEPKrylovConvergence(pep,PETSC_FALSE,pep->nconv,nv-pep->nconv,beta,&k));
255:     PetscCall((*pep->stopping)(pep,pep->its,pep->max_it,k,pep->nev,&pep->reason,pep->stoppingctx));
256:     nconv = k;

258:     /* Update l */
259:     if (pep->reason != PEP_CONVERGED_ITERATING || breakdown) l = 0;
260:     else {
261:       l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
262:       PetscCall(DSGetTruncateSize(pep->ds,k,nv,&l));
263:     }
264:     if (!ctx->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
265:     if (l) PetscCall(PetscInfo(pep,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));

267:     if (pep->reason == PEP_CONVERGED_ITERATING) {
268:       if (PetscUnlikely(breakdown)) {
269:         /* Stop if breakdown */
270:         PetscCall(PetscInfo(pep,"Breakdown Quadratic Arnoldi method (it=%" PetscInt_FMT " norm=%g)\n",pep->its,(double)beta));
271:         pep->reason = PEP_DIVERGED_BREAKDOWN;
272:       } else {
273:         /* Prepare the Rayleigh quotient for restart */
274:         PetscCall(DSTruncate(pep->ds,k+l,PETSC_FALSE));
275:       }
276:     }
277:     /* Update the corresponding vectors V(:,idx) = V*Q(:,idx) */
278:     PetscCall(DSGetMat(pep->ds,DS_MAT_Q,&Q));
279:     PetscCall(BVMultInPlace(pep->V,Q,pep->nconv,k+l));
280:     PetscCall(DSRestoreMat(pep->ds,DS_MAT_Q,&Q));

282:     pep->nconv = k;
283:     PetscCall(PEPMonitor(pep,pep->its,nconv,pep->eigr,pep->eigi,pep->errest,nv));
284:   }
285:   PetscCall(BVSetActiveColumns(pep->V,0,pep->nconv));
286:   for (j=0;j<pep->nconv;j++) {
287:     pep->eigr[j] *= pep->sfactor;
288:     pep->eigi[j] *= pep->sfactor;
289:   }

291:   PetscCall(STScaleShift(pep->st,sinv?1.0/pep->sfactor:pep->sfactor));
292:   PetscCall(RGPopScale(pep->rg));

294:   PetscCall(DSTruncate(pep->ds,pep->nconv,PETSC_TRUE));
295:   PetscCall(PetscFree(work));
296:   PetscFunctionReturn(PETSC_SUCCESS);
297: }

299: static PetscErrorCode PEPQArnoldiSetRestart_QArnoldi(PEP pep,PetscReal keep)
300: {
301:   PEP_QARNOLDI *ctx = (PEP_QARNOLDI*)pep->data;

303:   PetscFunctionBegin;
304:   if (keep==(PetscReal)PETSC_DEFAULT || keep==(PetscReal)PETSC_DECIDE) ctx->keep = 0.5;
305:   else {
306:     PetscCheck(keep>=0.1 && keep<=0.9,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"The keep argument must be in the range [0.1,0.9]");
307:     ctx->keep = keep;
308:   }
309:   PetscFunctionReturn(PETSC_SUCCESS);
310: }

312: /*@
313:    PEPQArnoldiSetRestart - Sets the restart parameter for the Q-Arnoldi
314:    method, in particular the proportion of basis vectors that must be kept
315:    after restart.

317:    Logically Collective

319:    Input Parameters:
320: +  pep  - the eigenproblem solver context
321: -  keep - the number of vectors to be kept at restart

323:    Options Database Key:
324: .  -pep_qarnoldi_restart - Sets the restart parameter

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

329:    Level: advanced

331: .seealso: PEPQArnoldiGetRestart()
332: @*/
333: PetscErrorCode PEPQArnoldiSetRestart(PEP pep,PetscReal keep)
334: {
335:   PetscFunctionBegin;
338:   PetscTryMethod(pep,"PEPQArnoldiSetRestart_C",(PEP,PetscReal),(pep,keep));
339:   PetscFunctionReturn(PETSC_SUCCESS);
340: }

342: static PetscErrorCode PEPQArnoldiGetRestart_QArnoldi(PEP pep,PetscReal *keep)
343: {
344:   PEP_QARNOLDI *ctx = (PEP_QARNOLDI*)pep->data;

346:   PetscFunctionBegin;
347:   *keep = ctx->keep;
348:   PetscFunctionReturn(PETSC_SUCCESS);
349: }

351: /*@
352:    PEPQArnoldiGetRestart - Gets the restart parameter used in the Q-Arnoldi method.

354:    Not Collective

356:    Input Parameter:
357: .  pep - the eigenproblem solver context

359:    Output Parameter:
360: .  keep - the restart parameter

362:    Level: advanced

364: .seealso: PEPQArnoldiSetRestart()
365: @*/
366: PetscErrorCode PEPQArnoldiGetRestart(PEP pep,PetscReal *keep)
367: {
368:   PetscFunctionBegin;
370:   PetscAssertPointer(keep,2);
371:   PetscUseMethod(pep,"PEPQArnoldiGetRestart_C",(PEP,PetscReal*),(pep,keep));
372:   PetscFunctionReturn(PETSC_SUCCESS);
373: }

375: static PetscErrorCode PEPQArnoldiSetLocking_QArnoldi(PEP pep,PetscBool lock)
376: {
377:   PEP_QARNOLDI *ctx = (PEP_QARNOLDI*)pep->data;

379:   PetscFunctionBegin;
380:   ctx->lock = lock;
381:   PetscFunctionReturn(PETSC_SUCCESS);
382: }

384: /*@
385:    PEPQArnoldiSetLocking - Choose between locking and non-locking variants of
386:    the Q-Arnoldi method.

388:    Logically Collective

390:    Input Parameters:
391: +  pep  - the eigenproblem solver context
392: -  lock - true if the locking variant must be selected

394:    Options Database Key:
395: .  -pep_qarnoldi_locking - Sets the locking flag

397:    Notes:
398:    The default is to lock converged eigenpairs when the method restarts.
399:    This behaviour can be changed so that all directions are kept in the
400:    working subspace even if already converged to working accuracy (the
401:    non-locking variant).

403:    Level: advanced

405: .seealso: PEPQArnoldiGetLocking()
406: @*/
407: PetscErrorCode PEPQArnoldiSetLocking(PEP pep,PetscBool lock)
408: {
409:   PetscFunctionBegin;
412:   PetscTryMethod(pep,"PEPQArnoldiSetLocking_C",(PEP,PetscBool),(pep,lock));
413:   PetscFunctionReturn(PETSC_SUCCESS);
414: }

416: static PetscErrorCode PEPQArnoldiGetLocking_QArnoldi(PEP pep,PetscBool *lock)
417: {
418:   PEP_QARNOLDI *ctx = (PEP_QARNOLDI*)pep->data;

420:   PetscFunctionBegin;
421:   *lock = ctx->lock;
422:   PetscFunctionReturn(PETSC_SUCCESS);
423: }

425: /*@
426:    PEPQArnoldiGetLocking - Gets the locking flag used in the Q-Arnoldi method.

428:    Not Collective

430:    Input Parameter:
431: .  pep - the eigenproblem solver context

433:    Output Parameter:
434: .  lock - the locking flag

436:    Level: advanced

438: .seealso: PEPQArnoldiSetLocking()
439: @*/
440: PetscErrorCode PEPQArnoldiGetLocking(PEP pep,PetscBool *lock)
441: {
442:   PetscFunctionBegin;
444:   PetscAssertPointer(lock,2);
445:   PetscUseMethod(pep,"PEPQArnoldiGetLocking_C",(PEP,PetscBool*),(pep,lock));
446:   PetscFunctionReturn(PETSC_SUCCESS);
447: }

449: static PetscErrorCode PEPSetFromOptions_QArnoldi(PEP pep,PetscOptionItems *PetscOptionsObject)
450: {
451:   PetscBool      flg,lock;
452:   PetscReal      keep;

454:   PetscFunctionBegin;
455:   PetscOptionsHeadBegin(PetscOptionsObject,"PEP Q-Arnoldi Options");

457:     PetscCall(PetscOptionsReal("-pep_qarnoldi_restart","Proportion of vectors kept after restart","PEPQArnoldiSetRestart",0.5,&keep,&flg));
458:     if (flg) PetscCall(PEPQArnoldiSetRestart(pep,keep));

460:     PetscCall(PetscOptionsBool("-pep_qarnoldi_locking","Choose between locking and non-locking variants","PEPQArnoldiSetLocking",PETSC_FALSE,&lock,&flg));
461:     if (flg) PetscCall(PEPQArnoldiSetLocking(pep,lock));

463:   PetscOptionsHeadEnd();
464:   PetscFunctionReturn(PETSC_SUCCESS);
465: }

467: static PetscErrorCode PEPView_QArnoldi(PEP pep,PetscViewer viewer)
468: {
469:   PEP_QARNOLDI   *ctx = (PEP_QARNOLDI*)pep->data;
470:   PetscBool      isascii;

472:   PetscFunctionBegin;
473:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
474:   if (isascii) {
475:     PetscCall(PetscViewerASCIIPrintf(viewer,"  %d%% of basis vectors kept after restart\n",(int)(100*ctx->keep)));
476:     PetscCall(PetscViewerASCIIPrintf(viewer,"  using the %slocking variant\n",ctx->lock?"":"non-"));
477:   }
478:   PetscFunctionReturn(PETSC_SUCCESS);
479: }

481: static PetscErrorCode PEPDestroy_QArnoldi(PEP pep)
482: {
483:   PetscFunctionBegin;
484:   PetscCall(PetscFree(pep->data));
485:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiSetRestart_C",NULL));
486:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiGetRestart_C",NULL));
487:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiSetLocking_C",NULL));
488:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiGetLocking_C",NULL));
489:   PetscFunctionReturn(PETSC_SUCCESS);
490: }

492: SLEPC_EXTERN PetscErrorCode PEPCreate_QArnoldi(PEP pep)
493: {
494:   PEP_QARNOLDI   *ctx;

496:   PetscFunctionBegin;
497:   PetscCall(PetscNew(&ctx));
498:   pep->data = (void*)ctx;

500:   pep->lineariz = PETSC_TRUE;
501:   ctx->lock     = PETSC_TRUE;

503:   pep->ops->solve          = PEPSolve_QArnoldi;
504:   pep->ops->setup          = PEPSetUp_QArnoldi;
505:   pep->ops->setfromoptions = PEPSetFromOptions_QArnoldi;
506:   pep->ops->destroy        = PEPDestroy_QArnoldi;
507:   pep->ops->view           = PEPView_QArnoldi;
508:   pep->ops->backtransform  = PEPBackTransform_Default;
509:   pep->ops->computevectors = PEPComputeVectors_Default;
510:   pep->ops->extractvectors = PEPExtractVectors_QArnoldi;
511:   pep->ops->setdefaultst   = PEPSetDefaultST_Transform;

513:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiSetRestart_C",PEPQArnoldiSetRestart_QArnoldi));
514:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiGetRestart_C",PEPQArnoldiGetRestart_QArnoldi));
515:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiSetLocking_C",PEPQArnoldiSetLocking_QArnoldi));
516:   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPQArnoldiGetLocking_C",PEPQArnoldiGetLocking_QArnoldi));
517:   PetscFunctionReturn(PETSC_SUCCESS);
518: }