Actual source code: epsdefault.c
slepc-3.16.2 2022-02-01
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2021, 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: {
22: STBackTransform(eps->st,eps->nconv,eps->eigr,eps->eigi);
23: return(0);
24: }
26: /*
27: EPSComputeVectors_Hermitian - Copies the Lanczos vectors as eigenvectors
28: using purification for generalized eigenproblems.
29: */
30: PetscErrorCode EPSComputeVectors_Hermitian(EPS eps)
31: {
33: PetscBool iscayley,indef;
34: Mat B,C;
37: if (eps->purify) {
38: EPS_Purify(eps,eps->nconv);
39: BVNormalize(eps->V,NULL);
40: } else {
41: /* In the case of Cayley transform, eigenvectors need to be B-normalized */
42: PetscObjectTypeCompare((PetscObject)eps->st,STCAYLEY,&iscayley);
43: if (iscayley && eps->isgeneralized) {
44: STGetMatrix(eps->st,1,&B);
45: BVGetMatrix(eps->V,&C,&indef);
46: if (indef) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"The inner product should not be indefinite");
47: BVSetMatrix(eps->V,B,PETSC_FALSE);
48: BVNormalize(eps->V,NULL);
49: BVSetMatrix(eps->V,C,PETSC_FALSE); /* restore original matrix */
50: }
51: }
52: return(0);
53: }
55: /*
56: EPSComputeVectors_Indefinite - similar to the Schur version but
57: for indefinite problems
58: */
59: PetscErrorCode EPSComputeVectors_Indefinite(EPS eps)
60: {
62: PetscInt n;
63: Mat X;
66: DSGetDimensions(eps->ds,&n,NULL,NULL,NULL);
67: DSVectors(eps->ds,DS_MAT_X,NULL,NULL);
68: DSGetMat(eps->ds,DS_MAT_X,&X);
69: BVMultInPlace(eps->V,X,0,n);
70: MatDestroy(&X);
72: /* purification */
73: if (eps->purify) { EPS_Purify(eps,eps->nconv); }
75: /* normalization */
76: BVNormalize(eps->V,eps->eigi);
77: return(0);
78: }
80: /*
81: EPSComputeVectors_Twosided - Adjust left eigenvectors in generalized problems: y = B^-* y.
82: */
83: PetscErrorCode EPSComputeVectors_Twosided(EPS eps)
84: {
86: PetscInt i;
87: Vec w,y;
90: if (!eps->twosided || !eps->isgeneralized) return(0);
91: EPSSetWorkVecs(eps,1);
92: w = eps->work[0];
93: for (i=0;i<eps->nconv;i++) {
94: BVCopyVec(eps->W,i,w);
95: VecConjugate(w);
96: BVGetColumn(eps->W,i,&y);
97: STMatSolveTranspose(eps->st,w,y);
98: VecConjugate(y);
99: BVRestoreColumn(eps->W,i,&y);
100: }
101: return(0);
102: }
104: /*
105: EPSComputeVectors_Schur - Compute eigenvectors from the vectors
106: provided by the eigensolver. This version is intended for solvers
107: that provide Schur vectors. Given the partial Schur decomposition
108: OP*V=V*T, the following steps are performed:
109: 1) compute eigenvectors of T: T*Z=Z*D
110: 2) compute eigenvectors of OP: X=V*Z
111: */
112: PetscErrorCode EPSComputeVectors_Schur(EPS eps)
113: {
115: PetscInt i;
116: Mat Z;
117: Vec z;
120: if (eps->ishermitian) {
121: if (eps->isgeneralized && !eps->ispositive) {
122: EPSComputeVectors_Indefinite(eps);
123: } else {
124: EPSComputeVectors_Hermitian(eps);
125: }
126: return(0);
127: }
129: /* right eigenvectors */
130: DSVectors(eps->ds,DS_MAT_X,NULL,NULL);
132: /* V = V * Z */
133: DSGetMat(eps->ds,DS_MAT_X,&Z);
134: BVMultInPlace(eps->V,Z,0,eps->nconv);
135: MatDestroy(&Z);
137: /* Purify eigenvectors */
138: if (eps->purify) { EPS_Purify(eps,eps->nconv); }
140: /* Fix eigenvectors if balancing was used */
141: if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
142: for (i=0;i<eps->nconv;i++) {
143: BVGetColumn(eps->V,i,&z);
144: VecPointwiseDivide(z,z,eps->D);
145: BVRestoreColumn(eps->V,i,&z);
146: }
147: }
149: /* normalize eigenvectors (when using purification or balancing) */
150: if (eps->purify || (eps->balance!=EPS_BALANCE_NONE && eps->D)) {
151: BVNormalize(eps->V,eps->eigi);
152: }
154: /* left eigenvectors */
155: if (eps->twosided) {
156: DSVectors(eps->ds,DS_MAT_Y,NULL,NULL);
157: /* W = W * Z */
158: DSGetMat(eps->ds,DS_MAT_Y,&Z);
159: BVMultInPlace(eps->W,Z,0,eps->nconv);
160: MatDestroy(&Z);
161: /* Fix left eigenvectors if balancing was used */
162: if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
163: for (i=0;i<eps->nconv;i++) {
164: BVGetColumn(eps->W,i,&z);
165: VecPointwiseMult(z,z,eps->D);
166: BVRestoreColumn(eps->W,i,&z);
167: }
168: }
169: EPSComputeVectors_Twosided(eps);
170: /* normalize */
171: BVNormalize(eps->W,eps->eigi);
172: #if !defined(PETSC_USE_COMPLEX)
173: for (i=0;i<eps->nconv-1;i++) {
174: if (eps->eigi[i] != 0.0) {
175: if (eps->eigi[i] > 0.0) { BVScaleColumn(eps->W,i+1,-1.0); }
176: i++;
177: }
178: }
179: #endif
180: }
181: return(0);
182: }
184: /*@
185: EPSSetWorkVecs - Sets a number of work vectors into an EPS object.
187: Collective on eps
189: Input Parameters:
190: + eps - eigensolver context
191: - nw - number of work vectors to allocate
193: Developers Note:
194: This is SLEPC_EXTERN because it may be required by user plugin EPS
195: implementations.
197: Level: developer
198: @*/
199: PetscErrorCode EPSSetWorkVecs(EPS eps,PetscInt nw)
200: {
202: Vec t;
207: if (nw <= 0) SETERRQ1(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"nw must be > 0: nw = %D",nw);
208: if (eps->nwork < nw) {
209: VecDestroyVecs(eps->nwork,&eps->work);
210: eps->nwork = nw;
211: BVGetColumn(eps->V,0,&t);
212: VecDuplicateVecs(t,nw,&eps->work);
213: BVRestoreColumn(eps->V,0,&t);
214: PetscLogObjectParents(eps,nw,eps->work);
215: }
216: return(0);
217: }
219: /*
220: EPSSetWhichEigenpairs_Default - Sets the default value for which,
221: depending on the ST.
222: */
223: PetscErrorCode EPSSetWhichEigenpairs_Default(EPS eps)
224: {
225: PetscBool target;
229: PetscObjectTypeCompareAny((PetscObject)eps->st,&target,STSINVERT,STCAYLEY,"");
230: if (target) eps->which = EPS_TARGET_MAGNITUDE;
231: else eps->which = EPS_LARGEST_MAGNITUDE;
232: return(0);
233: }
235: /*
236: EPSConvergedRelative - Checks convergence relative to the eigenvalue.
237: */
238: PetscErrorCode EPSConvergedRelative(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
239: {
240: PetscReal w;
243: w = SlepcAbsEigenvalue(eigr,eigi);
244: *errest = res/w;
245: return(0);
246: }
248: /*
249: EPSConvergedAbsolute - Checks convergence absolutely.
250: */
251: PetscErrorCode EPSConvergedAbsolute(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
252: {
254: *errest = res;
255: return(0);
256: }
258: /*
259: EPSConvergedNorm - Checks convergence relative to the eigenvalue and
260: the matrix norms.
261: */
262: PetscErrorCode EPSConvergedNorm(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
263: {
264: PetscReal w;
267: w = SlepcAbsEigenvalue(eigr,eigi);
268: *errest = res / (eps->nrma + w*eps->nrmb);
269: return(0);
270: }
272: /*@C
273: EPSStoppingBasic - Default routine to determine whether the outer eigensolver
274: iteration must be stopped.
276: Collective on eps
278: Input Parameters:
279: + eps - eigensolver context obtained from EPSCreate()
280: . its - current number of iterations
281: . max_it - maximum number of iterations
282: . nconv - number of currently converged eigenpairs
283: . nev - number of requested eigenpairs
284: - ctx - context (not used here)
286: Output Parameter:
287: . reason - result of the stopping test
289: Notes:
290: A positive value of reason indicates that the iteration has finished successfully
291: (converged), and a negative value indicates an error condition (diverged). If
292: the iteration needs to be continued, reason must be set to EPS_CONVERGED_ITERATING
293: (zero).
295: EPSStoppingBasic() will stop if all requested eigenvalues are converged, or if
296: the maximum number of iterations has been reached.
298: Use EPSSetStoppingTest() to provide your own test instead of using this one.
300: Level: advanced
302: .seealso: EPSSetStoppingTest(), EPSConvergedReason, EPSGetConvergedReason()
303: @*/
304: PetscErrorCode EPSStoppingBasic(EPS eps,PetscInt its,PetscInt max_it,PetscInt nconv,PetscInt nev,EPSConvergedReason *reason,void *ctx)
305: {
309: *reason = EPS_CONVERGED_ITERATING;
310: if (nconv >= nev) {
311: PetscInfo2(eps,"Linear eigensolver finished successfully: %D eigenpairs converged at iteration %D\n",nconv,its);
312: *reason = EPS_CONVERGED_TOL;
313: } else if (its >= max_it) {
314: *reason = EPS_DIVERGED_ITS;
315: PetscInfo1(eps,"Linear eigensolver iteration reached maximum number of iterations (%D)\n",its);
316: }
317: return(0);
318: }
320: /*
321: EPSComputeRitzVector - Computes the current Ritz vector.
323: Simple case (complex scalars or real scalars with Zi=NULL):
324: x = V*Zr (V is a basis of nv vectors, Zr has length nv)
326: Split case:
327: x = V*Zr y = V*Zi (Zr and Zi have length nv)
328: */
329: PetscErrorCode EPSComputeRitzVector(EPS eps,PetscScalar *Zr,PetscScalar *Zi,BV V,Vec x,Vec y)
330: {
332: PetscInt l,k;
333: PetscReal norm;
334: #if !defined(PETSC_USE_COMPLEX)
335: Vec z;
336: #endif
339: /* compute eigenvector */
340: BVGetActiveColumns(V,&l,&k);
341: BVSetActiveColumns(V,0,k);
342: BVMultVec(V,1.0,0.0,x,Zr);
344: /* purify eigenvector if necessary */
345: if (eps->purify) {
346: STApply(eps->st,x,y);
347: if (eps->ishermitian) {
348: BVNormVec(eps->V,y,NORM_2,&norm);
349: } else {
350: VecNorm(y,NORM_2,&norm);
351: }
352: VecScale(y,1.0/norm);
353: VecCopy(y,x);
354: }
355: /* fix eigenvector if balancing is used */
356: if (!eps->ishermitian && eps->balance!=EPS_BALANCE_NONE && eps->D) {
357: VecPointwiseDivide(x,x,eps->D);
358: }
359: #if !defined(PETSC_USE_COMPLEX)
360: /* compute imaginary part of eigenvector */
361: if (Zi) {
362: BVMultVec(V,1.0,0.0,y,Zi);
363: if (eps->ispositive) {
364: BVCreateVec(V,&z);
365: STApply(eps->st,y,z);
366: VecNorm(z,NORM_2,&norm);
367: VecScale(z,1.0/norm);
368: VecCopy(z,y);
369: VecDestroy(&z);
370: }
371: if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
372: VecPointwiseDivide(y,y,eps->D);
373: }
374: } else
375: #endif
376: { VecSet(y,0.0); }
378: /* normalize eigenvectors (when using balancing) */
379: if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
380: #if !defined(PETSC_USE_COMPLEX)
381: if (Zi) {
382: VecNormalizeComplex(x,y,PETSC_TRUE,NULL);
383: } else
384: #endif
385: {
386: VecNormalize(x,NULL);
387: }
388: }
389: BVSetActiveColumns(V,l,k);
390: return(0);
391: }
393: /*
394: EPSBuildBalance_Krylov - uses a Krylov subspace method to compute the
395: diagonal matrix to be applied for balancing in non-Hermitian problems.
396: */
397: PetscErrorCode EPSBuildBalance_Krylov(EPS eps)
398: {
399: Vec z,p,r;
400: PetscInt i,j;
401: PetscReal norma;
402: PetscScalar *pz,*pD;
403: const PetscScalar *pr,*pp;
404: PetscRandom rand;
405: PetscErrorCode ierr;
408: EPSSetWorkVecs(eps,3);
409: BVGetRandomContext(eps->V,&rand);
410: r = eps->work[0];
411: p = eps->work[1];
412: z = eps->work[2];
413: VecSet(eps->D,1.0);
415: for (j=0;j<eps->balance_its;j++) {
417: /* Build a random vector of +-1's */
418: VecSetRandom(z,rand);
419: VecGetArray(z,&pz);
420: for (i=0;i<eps->nloc;i++) {
421: if (PetscRealPart(pz[i])<0.5) pz[i]=-1.0;
422: else pz[i]=1.0;
423: }
424: VecRestoreArray(z,&pz);
426: /* Compute p=DA(D\z) */
427: VecPointwiseDivide(r,z,eps->D);
428: STApply(eps->st,r,p);
429: VecPointwiseMult(p,p,eps->D);
430: if (eps->balance == EPS_BALANCE_TWOSIDE) {
431: if (j==0) {
432: /* Estimate the matrix inf-norm */
433: VecAbs(p);
434: VecMax(p,NULL,&norma);
435: }
436: /* Compute r=D\(A'Dz) */
437: VecPointwiseMult(z,z,eps->D);
438: STApplyHermitianTranspose(eps->st,z,r);
439: VecPointwiseDivide(r,r,eps->D);
440: }
442: /* Adjust values of D */
443: VecGetArrayRead(r,&pr);
444: VecGetArrayRead(p,&pp);
445: VecGetArray(eps->D,&pD);
446: for (i=0;i<eps->nloc;i++) {
447: if (eps->balance == EPS_BALANCE_TWOSIDE) {
448: if (PetscAbsScalar(pp[i])>eps->balance_cutoff*norma && pr[i]!=0.0)
449: pD[i] *= PetscSqrtReal(PetscAbsScalar(pr[i]/pp[i]));
450: } else {
451: if (pp[i]!=0.0) pD[i] /= PetscAbsScalar(pp[i]);
452: }
453: }
454: VecRestoreArrayRead(r,&pr);
455: VecRestoreArrayRead(p,&pp);
456: VecRestoreArray(eps->D,&pD);
457: }
458: return(0);
459: }