Actual source code: epsdefault.c
slepc-3.20.2 2024-03-15
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(VecConjugate(w));
91: PetscCall(BVGetColumn(eps->W,i,&y));
92: PetscCall(STMatSolveTranspose(eps->st,w,y));
93: PetscCall(VecConjugate(y));
94: PetscCall(BVRestoreColumn(eps->W,i,&y));
95: }
96: PetscFunctionReturn(PETSC_SUCCESS);
97: }
99: /*
100: EPSComputeVectors_Schur - Compute eigenvectors from the vectors
101: provided by the eigensolver. This version is intended for solvers
102: that provide Schur vectors. Given the partial Schur decomposition
103: OP*V=V*T, the following steps are performed:
104: 1) compute eigenvectors of T: T*Z=Z*D
105: 2) compute eigenvectors of OP: X=V*Z
106: */
107: PetscErrorCode EPSComputeVectors_Schur(EPS eps)
108: {
109: PetscInt i;
110: Mat Z;
111: Vec z;
113: PetscFunctionBegin;
114: if (eps->ishermitian) {
115: if (eps->isgeneralized && !eps->ispositive) PetscCall(EPSComputeVectors_Indefinite(eps));
116: else PetscCall(EPSComputeVectors_Hermitian(eps));
117: PetscFunctionReturn(PETSC_SUCCESS);
118: }
120: /* right eigenvectors */
121: PetscCall(DSVectors(eps->ds,DS_MAT_X,NULL,NULL));
123: /* V = V * Z */
124: PetscCall(DSGetMat(eps->ds,DS_MAT_X,&Z));
125: PetscCall(BVMultInPlace(eps->V,Z,0,eps->nconv));
126: PetscCall(DSRestoreMat(eps->ds,DS_MAT_X,&Z));
128: /* Purify eigenvectors */
129: if (eps->purify) PetscCall(EPS_Purify(eps,eps->nconv));
131: /* Fix eigenvectors if balancing was used */
132: if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
133: for (i=0;i<eps->nconv;i++) {
134: PetscCall(BVGetColumn(eps->V,i,&z));
135: PetscCall(VecPointwiseDivide(z,z,eps->D));
136: PetscCall(BVRestoreColumn(eps->V,i,&z));
137: }
138: }
140: /* normalize eigenvectors (when using purification or balancing) */
141: if (eps->purify || (eps->balance!=EPS_BALANCE_NONE && eps->D)) PetscCall(BVNormalize(eps->V,eps->eigi));
143: /* left eigenvectors */
144: if (eps->twosided) {
145: PetscCall(DSVectors(eps->ds,DS_MAT_Y,NULL,NULL));
146: /* W = W * Z */
147: PetscCall(DSGetMat(eps->ds,DS_MAT_Y,&Z));
148: PetscCall(BVMultInPlace(eps->W,Z,0,eps->nconv));
149: PetscCall(DSRestoreMat(eps->ds,DS_MAT_Y,&Z));
150: /* Fix left eigenvectors if balancing was used */
151: if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
152: for (i=0;i<eps->nconv;i++) {
153: PetscCall(BVGetColumn(eps->W,i,&z));
154: PetscCall(VecPointwiseMult(z,z,eps->D));
155: PetscCall(BVRestoreColumn(eps->W,i,&z));
156: }
157: }
158: PetscCall(EPSComputeVectors_Twosided(eps));
159: /* normalize */
160: PetscCall(BVNormalize(eps->W,eps->eigi));
161: #if !defined(PETSC_USE_COMPLEX)
162: for (i=0;i<eps->nconv-1;i++) {
163: if (eps->eigi[i] != 0.0) {
164: if (eps->eigi[i] > 0.0) PetscCall(BVScaleColumn(eps->W,i+1,-1.0));
165: i++;
166: }
167: }
168: #endif
169: }
170: PetscFunctionReturn(PETSC_SUCCESS);
171: }
173: /*@
174: EPSSetWorkVecs - Sets a number of work vectors into an EPS object.
176: Collective
178: Input Parameters:
179: + eps - eigensolver context
180: - nw - number of work vectors to allocate
182: Developer Notes:
183: This is SLEPC_EXTERN because it may be required by user plugin EPS
184: implementations.
186: Level: developer
188: .seealso: EPSSetUp()
189: @*/
190: PetscErrorCode EPSSetWorkVecs(EPS eps,PetscInt nw)
191: {
192: Vec t;
194: PetscFunctionBegin;
197: PetscCheck(nw>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"nw must be > 0: nw = %" PetscInt_FMT,nw);
198: if (eps->nwork < nw) {
199: PetscCall(VecDestroyVecs(eps->nwork,&eps->work));
200: eps->nwork = nw;
201: PetscCall(BVGetColumn(eps->V,0,&t));
202: PetscCall(VecDuplicateVecs(t,nw,&eps->work));
203: PetscCall(BVRestoreColumn(eps->V,0,&t));
204: }
205: PetscFunctionReturn(PETSC_SUCCESS);
206: }
208: /*
209: EPSSetWhichEigenpairs_Default - Sets the default value for which,
210: depending on the ST.
211: */
212: PetscErrorCode EPSSetWhichEigenpairs_Default(EPS eps)
213: {
214: PetscBool target;
216: PetscFunctionBegin;
217: PetscCall(PetscObjectTypeCompareAny((PetscObject)eps->st,&target,STSINVERT,STCAYLEY,""));
218: if (target) eps->which = EPS_TARGET_MAGNITUDE;
219: else eps->which = EPS_LARGEST_MAGNITUDE;
220: PetscFunctionReturn(PETSC_SUCCESS);
221: }
223: /*
224: EPSConvergedRelative - Checks convergence relative to the eigenvalue.
225: */
226: PetscErrorCode EPSConvergedRelative(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
227: {
228: PetscReal w;
230: PetscFunctionBegin;
231: w = SlepcAbsEigenvalue(eigr,eigi);
232: *errest = res/w;
233: PetscFunctionReturn(PETSC_SUCCESS);
234: }
236: /*
237: EPSConvergedAbsolute - Checks convergence absolutely.
238: */
239: PetscErrorCode EPSConvergedAbsolute(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
240: {
241: PetscFunctionBegin;
242: *errest = res;
243: PetscFunctionReturn(PETSC_SUCCESS);
244: }
246: /*
247: EPSConvergedNorm - Checks convergence relative to the eigenvalue and
248: the matrix norms.
249: */
250: PetscErrorCode EPSConvergedNorm(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
251: {
252: PetscReal w;
254: PetscFunctionBegin;
255: w = SlepcAbsEigenvalue(eigr,eigi);
256: *errest = res / (eps->nrma + w*eps->nrmb);
257: PetscFunctionReturn(PETSC_SUCCESS);
258: }
260: /*@C
261: EPSStoppingBasic - Default routine to determine whether the outer eigensolver
262: iteration must be stopped.
264: Collective
266: Input Parameters:
267: + eps - eigensolver context obtained from EPSCreate()
268: . its - current number of iterations
269: . max_it - maximum number of iterations
270: . nconv - number of currently converged eigenpairs
271: . nev - number of requested eigenpairs
272: - ctx - context (not used here)
274: Output Parameter:
275: . reason - result of the stopping test
277: Notes:
278: A positive value of reason indicates that the iteration has finished successfully
279: (converged), and a negative value indicates an error condition (diverged). If
280: the iteration needs to be continued, reason must be set to EPS_CONVERGED_ITERATING
281: (zero).
283: EPSStoppingBasic() will stop if all requested eigenvalues are converged, or if
284: the maximum number of iterations has been reached.
286: Use EPSSetStoppingTest() to provide your own test instead of using this one.
288: Level: advanced
290: .seealso: EPSSetStoppingTest(), EPSConvergedReason, EPSGetConvergedReason()
291: @*/
292: PetscErrorCode EPSStoppingBasic(EPS eps,PetscInt its,PetscInt max_it,PetscInt nconv,PetscInt nev,EPSConvergedReason *reason,void *ctx)
293: {
294: PetscFunctionBegin;
295: *reason = EPS_CONVERGED_ITERATING;
296: if (nconv >= nev) {
297: PetscCall(PetscInfo(eps,"Linear eigensolver finished successfully: %" PetscInt_FMT " eigenpairs converged at iteration %" PetscInt_FMT "\n",nconv,its));
298: *reason = EPS_CONVERGED_TOL;
299: } else if (its >= max_it) {
300: *reason = EPS_DIVERGED_ITS;
301: PetscCall(PetscInfo(eps,"Linear eigensolver iteration reached maximum number of iterations (%" PetscInt_FMT ")\n",its));
302: }
303: PetscFunctionReturn(PETSC_SUCCESS);
304: }
306: /*
307: EPSComputeRitzVector - Computes the current Ritz vector.
309: Simple case (complex scalars or real scalars with Zi=NULL):
310: x = V*Zr (V is a basis of nv vectors, Zr has length nv)
312: Split case:
313: x = V*Zr y = V*Zi (Zr and Zi have length nv)
314: */
315: PetscErrorCode EPSComputeRitzVector(EPS eps,PetscScalar *Zr,PetscScalar *Zi,BV V,Vec x,Vec y)
316: {
317: PetscInt l,k;
318: PetscReal norm;
319: #if !defined(PETSC_USE_COMPLEX)
320: Vec z;
321: #endif
323: PetscFunctionBegin;
324: /* compute eigenvector */
325: PetscCall(BVGetActiveColumns(V,&l,&k));
326: PetscCall(BVSetActiveColumns(V,0,k));
327: PetscCall(BVMultVec(V,1.0,0.0,x,Zr));
329: /* purify eigenvector if necessary */
330: if (eps->purify) {
331: PetscCall(STApply(eps->st,x,y));
332: if (eps->ishermitian) PetscCall(BVNormVec(eps->V,y,NORM_2,&norm));
333: else PetscCall(VecNorm(y,NORM_2,&norm));
334: PetscCall(VecScale(y,1.0/norm));
335: PetscCall(VecCopy(y,x));
336: }
337: /* fix eigenvector if balancing is used */
338: if (!eps->ishermitian && eps->balance!=EPS_BALANCE_NONE && eps->D) PetscCall(VecPointwiseDivide(x,x,eps->D));
339: #if !defined(PETSC_USE_COMPLEX)
340: /* compute imaginary part of eigenvector */
341: if (Zi) {
342: PetscCall(BVMultVec(V,1.0,0.0,y,Zi));
343: if (eps->ispositive) {
344: PetscCall(BVCreateVec(V,&z));
345: PetscCall(STApply(eps->st,y,z));
346: PetscCall(VecNorm(z,NORM_2,&norm));
347: PetscCall(VecScale(z,1.0/norm));
348: PetscCall(VecCopy(z,y));
349: PetscCall(VecDestroy(&z));
350: }
351: if (eps->balance!=EPS_BALANCE_NONE && eps->D) PetscCall(VecPointwiseDivide(y,y,eps->D));
352: } else
353: #endif
354: PetscCall(VecSet(y,0.0));
356: /* normalize eigenvectors (when using balancing) */
357: if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
358: #if !defined(PETSC_USE_COMPLEX)
359: if (Zi) PetscCall(VecNormalizeComplex(x,y,PETSC_TRUE,NULL));
360: else
361: #endif
362: PetscCall(VecNormalize(x,NULL));
363: }
364: PetscCall(BVSetActiveColumns(V,l,k));
365: PetscFunctionReturn(PETSC_SUCCESS);
366: }
368: /*
369: EPSBuildBalance_Krylov - uses a Krylov subspace method to compute the
370: diagonal matrix to be applied for balancing in non-Hermitian problems.
371: */
372: PetscErrorCode EPSBuildBalance_Krylov(EPS eps)
373: {
374: Vec z,p,r;
375: PetscInt i,j;
376: PetscReal norma;
377: PetscScalar *pz,*pD;
378: const PetscScalar *pr,*pp;
379: PetscRandom rand;
381: PetscFunctionBegin;
382: PetscCall(EPSSetWorkVecs(eps,3));
383: PetscCall(BVGetRandomContext(eps->V,&rand));
384: r = eps->work[0];
385: p = eps->work[1];
386: z = eps->work[2];
387: PetscCall(VecSet(eps->D,1.0));
389: for (j=0;j<eps->balance_its;j++) {
391: /* Build a random vector of +-1's */
392: PetscCall(VecSetRandom(z,rand));
393: PetscCall(VecGetArray(z,&pz));
394: for (i=0;i<eps->nloc;i++) {
395: if (PetscRealPart(pz[i])<0.5) pz[i]=-1.0;
396: else pz[i]=1.0;
397: }
398: PetscCall(VecRestoreArray(z,&pz));
400: /* Compute p=DA(D\z) */
401: PetscCall(VecPointwiseDivide(r,z,eps->D));
402: PetscCall(STApply(eps->st,r,p));
403: PetscCall(VecPointwiseMult(p,p,eps->D));
404: if (eps->balance == EPS_BALANCE_TWOSIDE) {
405: if (j==0) {
406: /* Estimate the matrix inf-norm */
407: PetscCall(VecAbs(p));
408: PetscCall(VecMax(p,NULL,&norma));
409: }
410: /* Compute r=D\(A'Dz) */
411: PetscCall(VecPointwiseMult(z,z,eps->D));
412: PetscCall(STApplyHermitianTranspose(eps->st,z,r));
413: PetscCall(VecPointwiseDivide(r,r,eps->D));
414: }
416: /* Adjust values of D */
417: PetscCall(VecGetArrayRead(r,&pr));
418: PetscCall(VecGetArrayRead(p,&pp));
419: PetscCall(VecGetArray(eps->D,&pD));
420: for (i=0;i<eps->nloc;i++) {
421: if (eps->balance == EPS_BALANCE_TWOSIDE) {
422: if (PetscAbsScalar(pp[i])>eps->balance_cutoff*norma && pr[i]!=0.0)
423: pD[i] *= PetscSqrtReal(PetscAbsScalar(pr[i]/pp[i]));
424: } else {
425: if (pp[i]!=0.0) pD[i] /= PetscAbsScalar(pp[i]);
426: }
427: }
428: PetscCall(VecRestoreArrayRead(r,&pr));
429: PetscCall(VecRestoreArrayRead(p,&pp));
430: PetscCall(VecRestoreArray(eps->D,&pD));
431: }
432: PetscFunctionReturn(PETSC_SUCCESS);
433: }