Actual source code: epsdefault.c
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: This file contains some simple default routines for common operations
12: */
14: #include <slepc/private/epsimpl.h>
15: #include <slepcvec.h>
17: PetscErrorCode EPSBackTransform_Default(EPS eps)
18: {
19: PetscFunctionBegin;
20: PetscCall(STBackTransform(eps->st,eps->nconv,eps->eigr,eps->eigi));
21: PetscFunctionReturn(PETSC_SUCCESS);
22: }
24: /*
25: EPSComputeVectors_Hermitian - Copies the Lanczos vectors as eigenvectors
26: using purification for generalized eigenproblems.
27: */
28: PetscErrorCode EPSComputeVectors_Hermitian(EPS eps)
29: {
30: PetscBool iscayley,indef;
31: Mat B,C;
33: PetscFunctionBegin;
34: if (eps->purify) {
35: PetscCall(EPS_Purify(eps,eps->nconv));
36: PetscCall(BVNormalize(eps->V,NULL));
37: } else {
38: /* In the case of Cayley transform, eigenvectors need to be B-normalized */
39: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STCAYLEY,&iscayley));
40: if (iscayley && eps->isgeneralized) {
41: PetscCall(STGetMatrix(eps->st,1,&B));
42: PetscCall(BVGetMatrix(eps->V,&C,&indef));
43: PetscCheck(!indef,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"The inner product should not be indefinite");
44: PetscCall(BVSetMatrix(eps->V,B,PETSC_FALSE));
45: PetscCall(BVNormalize(eps->V,NULL));
46: PetscCall(BVSetMatrix(eps->V,C,PETSC_FALSE)); /* restore original matrix */
47: }
48: }
49: PetscFunctionReturn(PETSC_SUCCESS);
50: }
52: /*
53: EPSComputeVectors_Indefinite - similar to the Schur version but
54: for indefinite problems
55: */
56: PetscErrorCode EPSComputeVectors_Indefinite(EPS eps)
57: {
58: PetscInt n;
59: Mat X;
61: PetscFunctionBegin;
62: PetscCall(DSGetDimensions(eps->ds,&n,NULL,NULL,NULL));
63: PetscCall(DSVectors(eps->ds,DS_MAT_X,NULL,NULL));
64: PetscCall(DSGetMat(eps->ds,DS_MAT_X,&X));
65: PetscCall(BVMultInPlace(eps->V,X,0,n));
66: PetscCall(DSRestoreMat(eps->ds,DS_MAT_X,&X));
68: /* purification */
69: if (eps->purify) PetscCall(EPS_Purify(eps,eps->nconv));
71: /* normalization */
72: PetscCall(BVNormalize(eps->V,eps->eigi));
73: PetscFunctionReturn(PETSC_SUCCESS);
74: }
76: /*
77: EPSComputeVectors_Twosided - Adjust left eigenvectors in generalized problems: y = B^-* y.
78: */
79: PetscErrorCode EPSComputeVectors_Twosided(EPS eps)
80: {
81: PetscInt i;
82: Vec w,y;
84: PetscFunctionBegin;
85: if (!eps->twosided || !eps->isgeneralized) PetscFunctionReturn(PETSC_SUCCESS);
86: PetscCall(EPSSetWorkVecs(eps,1));
87: w = eps->work[0];
88: for (i=0;i<eps->nconv;i++) {
89: PetscCall(BVCopyVec(eps->W,i,w));
90: PetscCall(BVGetColumn(eps->W,i,&y));
91: PetscCall(STMatSolveHermitianTranspose(eps->st,w,y));
92: PetscCall(BVRestoreColumn(eps->W,i,&y));
93: }
94: PetscFunctionReturn(PETSC_SUCCESS);
95: }
97: /*
98: EPSComputeVectors_Schur - Compute eigenvectors from the vectors
99: provided by the eigensolver. This version is intended for solvers
100: that provide Schur vectors. Given the partial Schur decomposition
101: OP*V=V*T, the following steps are performed:
102: 1) compute eigenvectors of T: T*Z=Z*D
103: 2) compute eigenvectors of OP: X=V*Z
104: */
105: PetscErrorCode EPSComputeVectors_Schur(EPS eps)
106: {
107: PetscInt i;
108: Mat Z;
109: Vec z;
111: PetscFunctionBegin;
112: if (eps->ishermitian) {
113: if (eps->isgeneralized && !eps->ispositive) PetscCall(EPSComputeVectors_Indefinite(eps));
114: else PetscCall(EPSComputeVectors_Hermitian(eps));
115: PetscFunctionReturn(PETSC_SUCCESS);
116: }
118: /* right eigenvectors */
119: PetscCall(DSVectors(eps->ds,DS_MAT_X,NULL,NULL));
121: /* V = V * Z */
122: PetscCall(DSGetMat(eps->ds,DS_MAT_X,&Z));
123: PetscCall(BVMultInPlace(eps->V,Z,0,eps->nconv));
124: PetscCall(DSRestoreMat(eps->ds,DS_MAT_X,&Z));
126: /* Purify eigenvectors */
127: if (eps->purify) PetscCall(EPS_Purify(eps,eps->nconv));
129: /* Fix eigenvectors if balancing was used */
130: if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
131: for (i=0;i<eps->nconv;i++) {
132: PetscCall(BVGetColumn(eps->V,i,&z));
133: PetscCall(VecPointwiseDivide(z,z,eps->D));
134: PetscCall(BVRestoreColumn(eps->V,i,&z));
135: }
136: }
138: /* normalize eigenvectors (when using purification or balancing) */
139: if (eps->purify || (eps->balance!=EPS_BALANCE_NONE && eps->D)) PetscCall(BVNormalize(eps->V,eps->eigi));
141: /* left eigenvectors */
142: if (eps->twosided) {
143: PetscCall(DSVectors(eps->ds,DS_MAT_Y,NULL,NULL));
144: /* W = W * Z */
145: PetscCall(DSGetMat(eps->ds,DS_MAT_Y,&Z));
146: PetscCall(BVMultInPlace(eps->W,Z,0,eps->nconv));
147: PetscCall(DSRestoreMat(eps->ds,DS_MAT_Y,&Z));
148: /* Fix left eigenvectors if balancing was used */
149: if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
150: for (i=0;i<eps->nconv;i++) {
151: PetscCall(BVGetColumn(eps->W,i,&z));
152: PetscCall(VecPointwiseMult(z,z,eps->D));
153: PetscCall(BVRestoreColumn(eps->W,i,&z));
154: }
155: }
156: PetscCall(EPSComputeVectors_Twosided(eps));
157: /* normalize */
158: PetscCall(BVNormalize(eps->W,eps->eigi));
159: #if !defined(PETSC_USE_COMPLEX)
160: for (i=0;i<eps->nconv-1;i++) {
161: if (eps->eigi[i] != 0.0) {
162: if (eps->eigi[i] > 0.0) PetscCall(BVScaleColumn(eps->W,i+1,-1.0));
163: i++;
164: }
165: }
166: #endif
167: }
168: PetscFunctionReturn(PETSC_SUCCESS);
169: }
171: /*@
172: EPSSetWorkVecs - Sets a number of work vectors into an `EPS` object.
174: Collective
176: Input Parameters:
177: + eps - the linear eigensolver context
178: - nw - number of work vectors to allocate
180: Developer Note:
181: This is `SLEPC_EXTERN` because it may be required by user plugin `EPS`
182: implementations.
184: Level: developer
186: .seealso: [](ch:eps), `EPSSetUp()`
187: @*/
188: PetscErrorCode EPSSetWorkVecs(EPS eps,PetscInt nw)
189: {
190: Vec t;
192: PetscFunctionBegin;
195: PetscCheck(nw>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"nw must be > 0: nw = %" PetscInt_FMT,nw);
196: if (eps->nwork < nw) {
197: PetscCall(VecDestroyVecs(eps->nwork,&eps->work));
198: eps->nwork = nw;
199: PetscCall(BVGetColumn(eps->V,0,&t));
200: PetscCall(VecDuplicateVecs(t,nw,&eps->work));
201: PetscCall(BVRestoreColumn(eps->V,0,&t));
202: }
203: PetscFunctionReturn(PETSC_SUCCESS);
204: }
206: /*
207: EPSSetWhichEigenpairs_Default - Sets the default value for which,
208: depending on the ST.
209: */
210: PetscErrorCode EPSSetWhichEigenpairs_Default(EPS eps)
211: {
212: PetscBool target;
214: PetscFunctionBegin;
215: PetscCall(PetscObjectTypeCompareAny((PetscObject)eps->st,&target,STSINVERT,STCAYLEY,""));
216: if (target) eps->which = EPS_TARGET_MAGNITUDE;
217: else eps->which = EPS_LARGEST_MAGNITUDE;
218: PetscFunctionReturn(PETSC_SUCCESS);
219: }
221: /*
222: EPSConvergedRelative - Checks convergence relative to the eigenvalue.
223: */
224: PetscErrorCode EPSConvergedRelative(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
225: {
226: PetscReal w;
228: PetscFunctionBegin;
229: w = SlepcAbsEigenvalue(eigr,eigi);
230: *errest = (w!=0.0)? res/w: PETSC_MAX_REAL;
231: PetscFunctionReturn(PETSC_SUCCESS);
232: }
234: /*
235: EPSConvergedAbsolute - Checks convergence absolutely.
236: */
237: PetscErrorCode EPSConvergedAbsolute(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
238: {
239: PetscFunctionBegin;
240: *errest = res;
241: PetscFunctionReturn(PETSC_SUCCESS);
242: }
244: /*
245: EPSConvergedNorm - Checks convergence relative to the eigenvalue and
246: the matrix norms.
247: */
248: PetscErrorCode EPSConvergedNorm(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
249: {
250: PetscReal w;
252: PetscFunctionBegin;
253: w = SlepcAbsEigenvalue(eigr,eigi);
254: *errest = res / (eps->nrma + w*eps->nrmb);
255: PetscFunctionReturn(PETSC_SUCCESS);
256: }
258: /*@C
259: EPSStoppingBasic - Default routine to determine whether the outer eigensolver
260: iteration must be stopped.
262: Collective
264: Input Parameters:
265: + eps - the linear eigensolver context
266: . its - current number of iterations
267: . max_it - maximum number of iterations
268: . nconv - number of currently converged eigenpairs
269: . nev - number of requested eigenpairs
270: - ctx - context (not used here)
272: Output Parameter:
273: . reason - result of the stopping test
275: Notes:
276: `EPSStoppingBasic()` will stop if all requested eigenvalues are converged, or if
277: the maximum number of iterations has been reached.
279: This is the default stopping test.
280: Use `EPSSetStoppingTest()` to provide your own test instead of using this one.
282: Level: advanced
284: .seealso: [](ch:eps), `EPSSetStoppingTest()`, `EPSStoppingThreshold()`, `EPSConvergedReason`, `EPSGetConvergedReason()`
285: @*/
286: PetscErrorCode EPSStoppingBasic(EPS eps,PetscInt its,PetscInt max_it,PetscInt nconv,PetscInt nev,EPSConvergedReason *reason,void *ctx)
287: {
288: PetscFunctionBegin;
289: *reason = EPS_CONVERGED_ITERATING;
290: if (nconv >= nev) {
291: PetscCall(PetscInfo(eps,"Linear eigensolver finished successfully: %" PetscInt_FMT " eigenpairs converged at iteration %" PetscInt_FMT "\n",nconv,its));
292: *reason = EPS_CONVERGED_TOL;
293: } else if (its >= max_it) {
294: *reason = EPS_DIVERGED_ITS;
295: PetscCall(PetscInfo(eps,"Linear eigensolver iteration reached maximum number of iterations (%" PetscInt_FMT ")\n",its));
296: }
297: PetscFunctionReturn(PETSC_SUCCESS);
298: }
300: /*@C
301: EPSStoppingThreshold - Routine to determine whether the outer eigenvalue solver
302: iteration must be stopped, according to some threshold for the computed values.
304: Collective
306: Input Parameters:
307: + eps - the linear eigensolver context
308: . its - current number of iterations
309: . max_it - maximum number of iterations
310: . nconv - number of currently converged eigenpairs
311: . nev - number of requested eigenpairs
312: - ctx - context containing additional data (`EPSStoppingCtx`)
314: Output Parameter:
315: . reason - result of the stopping test
317: Notes:
318: `EPSStoppingThreshold()` will stop when one of the computed eigenvalues is not
319: above/below the threshold given at `EPSSetThreshold()`. If a number of wanted
320: eigenvalues has been specified via `EPSSetDimensions()` then it is also taken into
321: account, and the solver will stop when one of the two conditions (threshold or
322: number of converged values) is met.
324: Use `EPSSetStoppingTest()` to provide your own test instead of using this one.
326: Level: advanced
328: .seealso: [](ch:eps), `EPSSetStoppingTest()`, `EPSStoppingBasic()`, `EPSSetThreshold()`, `EPSSetDimensions()`, `EPSConvergedReason`, `EPSGetConvergedReason()`
329: @*/
330: PetscErrorCode EPSStoppingThreshold(EPS eps,PetscInt its,PetscInt max_it,PetscInt nconv,PetscInt nev,EPSConvergedReason *reason,void *ctx)
331: {
332: PetscReal thres,firstev,lastev,firstnc,errest;
333: PetscBool magnit,rel,isfilter;
334: PetscInt napprox;
335: EPSWhich which;
337: PetscFunctionBegin;
338: *reason = EPS_CONVERGED_ITERATING;
339: firstev = ((EPSStoppingCtx)ctx)->firstev;
340: lastev = ((EPSStoppingCtx)ctx)->lastev;
341: firstnc = ((EPSStoppingCtx)ctx)->firstnc;
342: errest = ((EPSStoppingCtx)ctx)->errest;
343: thres = ((EPSStoppingCtx)ctx)->thres;
344: rel = ((EPSStoppingCtx)ctx)->threlative;
345: napprox = ((EPSStoppingCtx)ctx)->napprox;
346: which = ((EPSStoppingCtx)ctx)->which;
348: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilter));
349: if (isfilter) {
350: /*
351: ((EPSStoppingCtx)ctx)->threlative should be PETSC_FALSE
352: ((EPSStoppingCtx)ctx)->which should be EPS_LARGEST_MAGNITUDE
353: */
354: if (nconv && (nconv==napprox || firstnc+errest<thres)) {
355: if (its==((EPSStoppingCtx)ctx)->its+1) {
356: if (nconv<napprox) PetscCall(PetscInfo(eps,"Linear eigensolver finished successfully: the eigenvalue magnitude %g (plus error estimate %g) is below the threshold %g\n",(double)firstnc,(double)errest,(double)thres));
357: else PetscCall(PetscInfo(eps,"Linear eigensolver finished successfully: all available eigenvalue approximations have converged\n"));
358: *reason = EPS_CONVERGED_TOL;
359: } else ((EPSStoppingCtx)ctx)->its = its; /* wait until next iteration */
360: } else if (nev && nconv >= nev) {
361: PetscCall(PetscInfo(eps,"Linear eigensolver finished successfully: %" PetscInt_FMT " eigenpairs converged at iteration %" PetscInt_FMT "\n",nconv,its));
362: *reason = EPS_CONVERGED_TOL;
363: } else if (its >= max_it) {
364: *reason = EPS_DIVERGED_ITS;
365: PetscCall(PetscInfo(eps,"Linear eigensolver iteration reached maximum number of iterations (%" PetscInt_FMT ")\n",its));
366: }
367: } else {
368: magnit = (which==EPS_SMALLEST_MAGNITUDE || which==EPS_LARGEST_MAGNITUDE || which==EPS_TARGET_MAGNITUDE)? PETSC_TRUE: PETSC_FALSE;
369: if (nconv && magnit && which==EPS_TARGET_MAGNITUDE && ((rel && ((thres>1.0 && lastev>thres*firstev) || (thres<1.0 && lastev<thres*firstev))) || (!rel && lastev>thres))) {
370: if (!rel) PetscCall(PetscInfo(eps,"Linear eigensolver finished successfully: the eigenvalue magnitude %g is above the threshold %g\n",(double)lastev,(double)thres));
371: else if (thres>1.0) PetscCall(PetscInfo(eps,"Linear eigensolver finished successfully: the ratio %g/%g is above the threshold %g\n",(double)lastev,(double)firstev,(double)thres));
372: else PetscCall(PetscInfo(eps,"Linear eigensolver finished successfully: the ratio %g/%g is below the threshold %g\n",(double)lastev,(double)firstev,(double)thres));
373: *reason = EPS_CONVERGED_TOL;
374: } else if (nconv && magnit && ((which==EPS_LARGEST_MAGNITUDE && ((rel && lastev<thres*firstev) || (!rel && lastev<thres))) || (which==EPS_SMALLEST_MAGNITUDE && lastev>thres))) {
375: if (which==EPS_SMALLEST_MAGNITUDE) PetscCall(PetscInfo(eps,"Linear eigensolver finished successfully: the eigenvalue magnitude %g is above the threshold %g\n",(double)lastev,(double)thres));
376: else if (!rel) PetscCall(PetscInfo(eps,"Linear eigensolver finished successfully: the eigenvalue magnitude %g is below the threshold %g\n",(double)lastev,(double)thres));
377: else PetscCall(PetscInfo(eps,"Linear eigensolver finished successfully: the ratio %g/%g is below the threshold %g\n",(double)lastev,(double)firstev,(double)thres));
378: *reason = EPS_CONVERGED_TOL;
379: } else if (nconv && !magnit && ((which==EPS_LARGEST_REAL && lastev<thres) || (which==EPS_SMALLEST_REAL && lastev>thres))) {
380: if (which==EPS_LARGEST_REAL) PetscCall(PetscInfo(eps,"Linear eigensolver finished successfully: eigenvalue %g is below the threshold %g\n",(double)lastev,(double)thres));
381: else PetscCall(PetscInfo(eps,"Linear eigensolver finished successfully: eigenvalue %g is above the threshold %g\n",(double)lastev,(double)thres));
382: *reason = EPS_CONVERGED_TOL;
383: } else if (nev && nconv >= nev) {
384: PetscCall(PetscInfo(eps,"Linear eigensolver finished successfully: %" PetscInt_FMT " eigenpairs converged at iteration %" PetscInt_FMT "\n",nconv,its));
385: *reason = EPS_CONVERGED_TOL;
386: } else if (its >= max_it) {
387: *reason = EPS_DIVERGED_ITS;
388: PetscCall(PetscInfo(eps,"Linear eigensolver iteration reached maximum number of iterations (%" PetscInt_FMT ")\n",its));
389: }
390: }
391: PetscFunctionReturn(PETSC_SUCCESS);
392: }
394: /*
395: EPSComputeRitzVector - Computes the current Ritz vector.
397: Simple case (complex scalars or real scalars with Zi=NULL):
398: x = V*Zr (V is a basis of nv vectors, Zr has length nv)
400: Split case:
401: x = V*Zr y = V*Zi (Zr and Zi have length nv)
402: */
403: PetscErrorCode EPSComputeRitzVector(EPS eps,PetscScalar *Zr,PetscScalar *Zi,BV V,Vec x,Vec y)
404: {
405: PetscInt l,k;
406: PetscReal norm;
407: #if !defined(PETSC_USE_COMPLEX)
408: Vec z;
409: #endif
411: PetscFunctionBegin;
412: /* compute eigenvector */
413: PetscCall(BVGetActiveColumns(V,&l,&k));
414: PetscCall(BVSetActiveColumns(V,0,k));
415: PetscCall(BVMultVec(V,1.0,0.0,x,Zr));
417: /* purify eigenvector if necessary */
418: if (eps->purify) {
419: PetscCall(STApply(eps->st,x,y));
420: if (eps->ishermitian) PetscCall(BVNormVec(eps->V,y,NORM_2,&norm));
421: else PetscCall(VecNorm(y,NORM_2,&norm));
422: PetscCall(VecScale(y,1.0/norm));
423: PetscCall(VecCopy(y,x));
424: }
425: /* fix eigenvector if balancing is used */
426: if (!eps->ishermitian && eps->balance!=EPS_BALANCE_NONE && eps->D) PetscCall(VecPointwiseDivide(x,x,eps->D));
427: #if !defined(PETSC_USE_COMPLEX)
428: /* compute imaginary part of eigenvector */
429: if (Zi) {
430: PetscCall(BVMultVec(V,1.0,0.0,y,Zi));
431: if (eps->ispositive) {
432: PetscCall(BVCreateVec(V,&z));
433: PetscCall(STApply(eps->st,y,z));
434: PetscCall(VecNorm(z,NORM_2,&norm));
435: PetscCall(VecScale(z,1.0/norm));
436: PetscCall(VecCopy(z,y));
437: PetscCall(VecDestroy(&z));
438: }
439: if (eps->balance!=EPS_BALANCE_NONE && eps->D) PetscCall(VecPointwiseDivide(y,y,eps->D));
440: } else
441: #endif
442: PetscCall(VecSet(y,0.0));
444: /* normalize eigenvectors (when using balancing) */
445: if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
446: #if !defined(PETSC_USE_COMPLEX)
447: if (Zi) PetscCall(VecNormalizeComplex(x,y,PETSC_TRUE,NULL));
448: else
449: #endif
450: PetscCall(VecNormalize(x,NULL));
451: }
452: PetscCall(BVSetActiveColumns(V,l,k));
453: PetscFunctionReturn(PETSC_SUCCESS);
454: }
456: /*
457: EPSBuildBalance_Krylov - uses a Krylov subspace method to compute the
458: diagonal matrix to be applied for balancing in non-Hermitian problems.
459: */
460: PetscErrorCode EPSBuildBalance_Krylov(EPS eps)
461: {
462: Vec z,p,r;
463: PetscInt i,j;
464: PetscReal norma;
465: PetscScalar *pz,*pD;
466: const PetscScalar *pr,*pp;
467: PetscRandom rand;
469: PetscFunctionBegin;
470: PetscCall(EPSSetWorkVecs(eps,3));
471: PetscCall(BVGetRandomContext(eps->V,&rand));
472: r = eps->work[0];
473: p = eps->work[1];
474: z = eps->work[2];
475: PetscCall(VecSet(eps->D,1.0));
477: for (j=0;j<eps->balance_its;j++) {
479: /* Build a random vector of +-1's */
480: PetscCall(VecSetRandom(z,rand));
481: PetscCall(VecGetArray(z,&pz));
482: for (i=0;i<eps->nloc;i++) {
483: if (PetscRealPart(pz[i])<0.5) pz[i]=-1.0;
484: else pz[i]=1.0;
485: }
486: PetscCall(VecRestoreArray(z,&pz));
488: /* Compute p=DA(D\z) */
489: PetscCall(VecPointwiseDivide(r,z,eps->D));
490: PetscCall(STApply(eps->st,r,p));
491: PetscCall(VecPointwiseMult(p,p,eps->D));
492: if (eps->balance == EPS_BALANCE_TWOSIDE) {
493: if (j==0) {
494: /* Estimate the matrix inf-norm */
495: PetscCall(VecAbs(p));
496: PetscCall(VecMax(p,NULL,&norma));
497: }
498: /* Compute r=D\(A'Dz) */
499: PetscCall(VecPointwiseMult(z,z,eps->D));
500: PetscCall(STApplyHermitianTranspose(eps->st,z,r));
501: PetscCall(VecPointwiseDivide(r,r,eps->D));
502: }
504: /* Adjust values of D */
505: PetscCall(VecGetArrayRead(r,&pr));
506: PetscCall(VecGetArrayRead(p,&pp));
507: PetscCall(VecGetArray(eps->D,&pD));
508: for (i=0;i<eps->nloc;i++) {
509: if (eps->balance == EPS_BALANCE_TWOSIDE) {
510: if (PetscAbsScalar(pp[i])>eps->balance_cutoff*norma && pr[i]!=0.0)
511: pD[i] *= PetscSqrtReal(PetscAbsScalar(pr[i]/pp[i]));
512: } else {
513: if (pp[i]!=0.0) pD[i] /= PetscAbsScalar(pp[i]);
514: }
515: }
516: PetscCall(VecRestoreArrayRead(r,&pr));
517: PetscCall(VecRestoreArrayRead(p,&pp));
518: PetscCall(VecRestoreArray(eps->D,&pD));
519: }
520: PetscFunctionReturn(PETSC_SUCCESS);
521: }