Actual source code: epsdefault.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: This file contains some simple default routines for common operations
12: */
14: #include <slepc/private/epsimpl.h> /*I "slepceps.h" I*/
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: PetscInt i;
34: PetscScalar dot;
35: PetscReal norm;
36: PetscBool iscayley;
37: Mat B;
38: Vec w,x;
41: if (eps->purify) {
42: EPS_Purify(eps,eps->nconv);
43: for (i=0;i<eps->nconv;i++) {
44: BVNormColumn(eps->V,i,NORM_2,&norm);
45: BVScaleColumn(eps->V,i,1.0/norm);
46: }
47: } else {
48: /* In the case of Cayley transform, eigenvectors need to be B-normalized */
49: PetscObjectTypeCompare((PetscObject)eps->st,STCAYLEY,&iscayley);
50: if (iscayley && eps->isgeneralized) {
51: STGetMatrix(eps->st,1,&B);
52: MatCreateVecs(B,NULL,&w);
53: for (i=0;i<eps->nconv;i++) {
54: BVGetColumn(eps->V,i,&x);
55: MatMult(B,x,w);
56: VecDot(w,x,&dot);
57: VecScale(x,1.0/PetscSqrtScalar(dot));
58: BVRestoreColumn(eps->V,i,&x);
59: }
60: VecDestroy(&w);
61: }
62: }
63: return(0);
64: }
66: /*
67: EPSComputeVectors_Indefinite - similar to the Schur version but
68: for indefinite problems
69: */
70: PetscErrorCode EPSComputeVectors_Indefinite(EPS eps)
71: {
73: PetscInt n,i;
74: Mat X;
75: Vec v;
76: #if !defined(PETSC_USE_COMPLEX)
77: Vec v1;
78: #endif
81: DSGetDimensions(eps->ds,&n,NULL,NULL,NULL,NULL);
82: DSVectors(eps->ds,DS_MAT_X,NULL,NULL);
83: DSGetMat(eps->ds,DS_MAT_X,&X);
84: BVSetActiveColumns(eps->V,0,n);
85: BVMultInPlace(eps->V,X,0,n);
86: MatDestroy(&X);
88: /* purification */
89: if (eps->purify) { EPS_Purify(eps,eps->nconv); }
91: /* normalization */
92: for (i=0;i<n;i++) {
93: #if !defined(PETSC_USE_COMPLEX)
94: if (eps->eigi[i] != 0.0) {
95: BVGetColumn(eps->V,i,&v);
96: BVGetColumn(eps->V,i+1,&v1);
97: VecNormalizeComplex(v,v1,PETSC_TRUE,NULL);
98: BVRestoreColumn(eps->V,i,&v);
99: BVRestoreColumn(eps->V,i+1,&v1);
100: i++;
101: } else
102: #endif
103: {
104: BVGetColumn(eps->V,i,&v);
105: VecNormalize(v,NULL);
106: BVRestoreColumn(eps->V,i,&v);
107: }
108: }
109: return(0);
110: }
112: /*
113: EPSComputeVectors_Twosided - Adjust left eigenvectors in generalized problems: y = B^-* y.
114: */
115: PetscErrorCode EPSComputeVectors_Twosided(EPS eps)
116: {
118: PetscInt i;
119: Vec w,y;
122: if (!eps->twosided || !eps->isgeneralized) return(0);
123: EPSSetWorkVecs(eps,1);
124: w = eps->work[0];
125: for (i=0;i<eps->nconv;i++) {
126: BVCopyVec(eps->W,i,w);
127: VecConjugate(w);
128: BVGetColumn(eps->W,i,&y);
129: STMatSolveTranspose(eps->st,w,y);
130: VecConjugate(y);
131: BVRestoreColumn(eps->W,i,&y);
132: }
133: return(0);
134: }
136: /*
137: EPSComputeVectors_Schur - Compute eigenvectors from the vectors
138: provided by the eigensolver. This version is intended for solvers
139: that provide Schur vectors. Given the partial Schur decomposition
140: OP*V=V*T, the following steps are performed:
141: 1) compute eigenvectors of T: T*Z=Z*D
142: 2) compute eigenvectors of OP: X=V*Z
143: */
144: PetscErrorCode EPSComputeVectors_Schur(EPS eps)
145: {
147: PetscInt n,i;
148: Mat Z;
149: Vec z,v;
150: #if !defined(PETSC_USE_COMPLEX)
151: Vec v1;
152: #endif
155: if (eps->ishermitian) {
156: if (eps->isgeneralized && !eps->ispositive) {
157: EPSComputeVectors_Indefinite(eps);
158: } else {
159: EPSComputeVectors_Hermitian(eps);
160: }
161: return(0);
162: }
163: DSGetDimensions(eps->ds,&n,NULL,NULL,NULL,NULL);
165: /* right eigenvectors */
166: DSVectors(eps->ds,DS_MAT_X,NULL,NULL);
168: /* V = V * Z */
169: DSGetMat(eps->ds,DS_MAT_X,&Z);
170: BVSetActiveColumns(eps->V,0,n);
171: BVMultInPlace(eps->V,Z,0,n);
172: MatDestroy(&Z);
174: /* Purify eigenvectors */
175: if (eps->purify) { EPS_Purify(eps,eps->nconv); }
177: /* Fix eigenvectors if balancing was used */
178: if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
179: for (i=0;i<n;i++) {
180: BVGetColumn(eps->V,i,&z);
181: VecPointwiseDivide(z,z,eps->D);
182: BVRestoreColumn(eps->V,i,&z);
183: }
184: }
186: /* normalize eigenvectors (when using purification or balancing) */
187: if (eps->purify || (eps->balance!=EPS_BALANCE_NONE && eps->D)) {
188: for (i=0;i<n;i++) {
189: #if !defined(PETSC_USE_COMPLEX)
190: if (eps->eigi[i] != 0.0) {
191: BVGetColumn(eps->V,i,&v);
192: BVGetColumn(eps->V,i+1,&v1);
193: VecNormalizeComplex(v,v1,PETSC_TRUE,NULL);
194: BVRestoreColumn(eps->V,i,&v);
195: BVRestoreColumn(eps->V,i+1,&v1);
196: i++;
197: } else
198: #endif
199: {
200: BVGetColumn(eps->V,i,&v);
201: VecNormalize(v,NULL);
202: BVRestoreColumn(eps->V,i,&v);
203: }
204: }
205: }
207: /* left eigenvectors from eps->dsts */
208: if (eps->twosided) {
209: DSVectors(eps->dsts,DS_MAT_X,NULL,NULL);
210: /* W = W * Z */
211: DSGetMat(eps->dsts,DS_MAT_X,&Z);
212: BVSetActiveColumns(eps->W,0,n);
213: BVMultInPlace(eps->W,Z,0,n);
214: MatDestroy(&Z);
215: /* Fix left eigenvectors if balancing was used */
216: if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
217: for (i=0;i<n;i++) {
218: BVGetColumn(eps->W,i,&z);
219: VecPointwiseMult(z,z,eps->D);
220: BVRestoreColumn(eps->W,i,&z);
221: }
222: }
223: EPSComputeVectors_Twosided(eps);
224: /* normalize */
225: for (i=0;i<eps->nconv;i++) {
226: #if !defined(PETSC_USE_COMPLEX)
227: if (eps->eigi[i]!=0.0) { /* first eigenvalue of a complex conjugate pair */
228: BVGetColumn(eps->W,i,&v);
229: BVGetColumn(eps->W,i+1,&v1);
230: VecNormalizeComplex(v,v1,PETSC_TRUE,NULL);
231: BVRestoreColumn(eps->W,i,&v);
232: BVRestoreColumn(eps->W,i+1,&v1);
233: i++;
234: } else /* real eigenvalue */
235: #endif
236: {
237: BVGetColumn(eps->W,i,&v);
238: VecNormalize(v,NULL);
239: BVRestoreColumn(eps->W,i,&v);
240: }
241: }
242: #if !defined(PETSC_USE_COMPLEX)
243: for (i=0;i<eps->nconv-1;i++) {
244: if (eps->eigi[i] != 0.0) {
245: if (eps->eigi[i] > 0.0) { BVScaleColumn(eps->W,i+1,-1.0); }
246: i++;
247: }
248: }
249: #endif
250: }
251: return(0);
252: }
254: /*@
255: EPSSetWorkVecs - Sets a number of work vectors into an EPS object.
257: Collective on EPS
259: Input Parameters:
260: + eps - eigensolver context
261: - nw - number of work vectors to allocate
263: Developers Note:
264: This is SLEPC_EXTERN because it may be required by user plugin EPS
265: implementations.
267: Level: developer
268: @*/
269: PetscErrorCode EPSSetWorkVecs(EPS eps,PetscInt nw)
270: {
272: Vec t;
277: if (nw <= 0) SETERRQ1(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"nw must be > 0: nw = %D",nw);
278: if (eps->nwork < nw) {
279: VecDestroyVecs(eps->nwork,&eps->work);
280: eps->nwork = nw;
281: BVGetColumn(eps->V,0,&t);
282: VecDuplicateVecs(t,nw,&eps->work);
283: BVRestoreColumn(eps->V,0,&t);
284: PetscLogObjectParents(eps,nw,eps->work);
285: }
286: return(0);
287: }
289: /*
290: EPSSetWhichEigenpairs_Default - Sets the default value for which,
291: depending on the ST.
292: */
293: PetscErrorCode EPSSetWhichEigenpairs_Default(EPS eps)
294: {
295: PetscBool target;
299: PetscObjectTypeCompareAny((PetscObject)eps->st,&target,STSINVERT,STCAYLEY,"");
300: if (target) eps->which = EPS_TARGET_MAGNITUDE;
301: else eps->which = EPS_LARGEST_MAGNITUDE;
302: return(0);
303: }
305: /*
306: EPSConvergedRelative - Checks convergence relative to the eigenvalue.
307: */
308: PetscErrorCode EPSConvergedRelative(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
309: {
310: PetscReal w;
313: w = SlepcAbsEigenvalue(eigr,eigi);
314: *errest = res/w;
315: return(0);
316: }
318: /*
319: EPSConvergedAbsolute - Checks convergence absolutely.
320: */
321: PetscErrorCode EPSConvergedAbsolute(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
322: {
324: *errest = res;
325: return(0);
326: }
328: /*
329: EPSConvergedNorm - Checks convergence relative to the eigenvalue and
330: the matrix norms.
331: */
332: PetscErrorCode EPSConvergedNorm(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
333: {
334: PetscReal w;
337: w = SlepcAbsEigenvalue(eigr,eigi);
338: *errest = res / (eps->nrma + w*eps->nrmb);
339: return(0);
340: }
342: /*@C
343: EPSStoppingBasic - Default routine to determine whether the outer eigensolver
344: iteration must be stopped.
346: Collective on EPS
348: Input Parameters:
349: + eps - eigensolver context obtained from EPSCreate()
350: . its - current number of iterations
351: . max_it - maximum number of iterations
352: . nconv - number of currently converged eigenpairs
353: . nev - number of requested eigenpairs
354: - ctx - context (not used here)
356: Output Parameter:
357: . reason - result of the stopping test
359: Notes:
360: A positive value of reason indicates that the iteration has finished successfully
361: (converged), and a negative value indicates an error condition (diverged). If
362: the iteration needs to be continued, reason must be set to EPS_CONVERGED_ITERATING
363: (zero).
365: EPSStoppingBasic() will stop if all requested eigenvalues are converged, or if
366: the maximum number of iterations has been reached.
368: Use EPSSetStoppingTest() to provide your own test instead of using this one.
370: Level: advanced
372: .seealso: EPSSetStoppingTest(), EPSConvergedReason, EPSGetConvergedReason()
373: @*/
374: PetscErrorCode EPSStoppingBasic(EPS eps,PetscInt its,PetscInt max_it,PetscInt nconv,PetscInt nev,EPSConvergedReason *reason,void *ctx)
375: {
379: *reason = EPS_CONVERGED_ITERATING;
380: if (nconv >= nev) {
381: PetscInfo2(eps,"Linear eigensolver finished successfully: %D eigenpairs converged at iteration %D\n",nconv,its);
382: *reason = EPS_CONVERGED_TOL;
383: } else if (its >= max_it) {
384: *reason = EPS_DIVERGED_ITS;
385: PetscInfo1(eps,"Linear eigensolver iteration reached maximum number of iterations (%D)\n",its);
386: }
387: return(0);
388: }
390: /*
391: EPSComputeRitzVector - Computes the current Ritz vector.
393: Simple case (complex scalars or real scalars with Zi=NULL):
394: x = V*Zr (V is a basis of nv vectors, Zr has length nv)
396: Split case:
397: x = V*Zr y = V*Zi (Zr and Zi have length nv)
398: */
399: PetscErrorCode EPSComputeRitzVector(EPS eps,PetscScalar *Zr,PetscScalar *Zi,BV V,Vec x,Vec y)
400: {
402: PetscInt l,k;
403: PetscReal norm;
404: #if !defined(PETSC_USE_COMPLEX)
405: Vec z;
406: #endif
409: /* compute eigenvector */
410: BVGetActiveColumns(V,&l,&k);
411: BVSetActiveColumns(V,0,k);
412: BVMultVec(V,1.0,0.0,x,Zr);
414: /* purify eigenvector if necessary */
415: if (eps->purify) {
416: STApply(eps->st,x,y);
417: if (eps->ishermitian) {
418: BVNormVec(eps->V,y,NORM_2,&norm);
419: } else {
420: VecNorm(y,NORM_2,&norm);
421: }
422: VecScale(y,1.0/norm);
423: VecCopy(y,x);
424: }
425: /* fix eigenvector if balancing is used */
426: if (!eps->ishermitian && eps->balance!=EPS_BALANCE_NONE && eps->D) {
427: VecPointwiseDivide(x,x,eps->D);
428: }
429: #if !defined(PETSC_USE_COMPLEX)
430: /* compute imaginary part of eigenvector */
431: if (Zi) {
432: BVMultVec(V,1.0,0.0,y,Zi);
433: if (eps->ispositive) {
434: BVCreateVec(V,&z);
435: STApply(eps->st,y,z);
436: VecNorm(z,NORM_2,&norm);
437: VecScale(z,1.0/norm);
438: VecCopy(z,y);
439: VecDestroy(&z);
440: }
441: if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
442: VecPointwiseDivide(y,y,eps->D);
443: }
444: } else
445: #endif
446: { VecSet(y,0.0); }
448: /* normalize eigenvectors (when using balancing) */
449: if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
450: #if !defined(PETSC_USE_COMPLEX)
451: if (Zi) {
452: VecNormalizeComplex(x,y,PETSC_TRUE,NULL);
453: } else
454: #endif
455: {
456: VecNormalize(x,NULL);
457: }
458: }
459: BVSetActiveColumns(V,l,k);
460: return(0);
461: }
463: /*
464: EPSBuildBalance_Krylov - uses a Krylov subspace method to compute the
465: diagonal matrix to be applied for balancing in non-Hermitian problems.
466: */
467: PetscErrorCode EPSBuildBalance_Krylov(EPS eps)
468: {
469: Vec z,p,r;
470: PetscInt i,j;
471: PetscReal norma;
472: PetscScalar *pz,*pD;
473: const PetscScalar *pr,*pp;
474: PetscRandom rand;
475: PetscErrorCode ierr;
478: EPSSetWorkVecs(eps,3);
479: BVGetRandomContext(eps->V,&rand);
480: r = eps->work[0];
481: p = eps->work[1];
482: z = eps->work[2];
483: VecSet(eps->D,1.0);
485: for (j=0;j<eps->balance_its;j++) {
487: /* Build a random vector of +-1's */
488: VecSetRandom(z,rand);
489: VecGetArray(z,&pz);
490: for (i=0;i<eps->nloc;i++) {
491: if (PetscRealPart(pz[i])<0.5) pz[i]=-1.0;
492: else pz[i]=1.0;
493: }
494: VecRestoreArray(z,&pz);
496: /* Compute p=DA(D\z) */
497: VecPointwiseDivide(r,z,eps->D);
498: STApply(eps->st,r,p);
499: VecPointwiseMult(p,p,eps->D);
500: if (eps->balance == EPS_BALANCE_TWOSIDE) {
501: if (j==0) {
502: /* Estimate the matrix inf-norm */
503: VecAbs(p);
504: VecMax(p,NULL,&norma);
505: }
506: /* Compute r=D\(A'Dz) */
507: VecPointwiseMult(z,z,eps->D);
508: STApplyTranspose(eps->st,z,r);
509: VecPointwiseDivide(r,r,eps->D);
510: }
512: /* Adjust values of D */
513: VecGetArrayRead(r,&pr);
514: VecGetArrayRead(p,&pp);
515: VecGetArray(eps->D,&pD);
516: for (i=0;i<eps->nloc;i++) {
517: if (eps->balance == EPS_BALANCE_TWOSIDE) {
518: if (PetscAbsScalar(pp[i])>eps->balance_cutoff*norma && pr[i]!=0.0)
519: pD[i] *= PetscSqrtReal(PetscAbsScalar(pr[i]/pp[i]));
520: } else {
521: if (pp[i]!=0.0) pD[i] /= PetscAbsScalar(pp[i]);
522: }
523: }
524: VecRestoreArrayRead(r,&pr);
525: VecRestoreArrayRead(p,&pp);
526: VecRestoreArray(eps->D,&pD);
527: }
528: return(0);
529: }