Actual source code: pepdefault.c
slepc-3.21.0 2024-03-30
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: Simple default routines for common PEP operations
12: */
14: #include <slepc/private/pepimpl.h>
16: /*@
17: PEPSetWorkVecs - Sets a number of work vectors into a PEP object.
19: Collective
21: Input Parameters:
22: + pep - polynomial eigensolver context
23: - nw - number of work vectors to allocate
25: Developer Notes:
26: This is SLEPC_EXTERN because it may be required by user plugin PEP
27: implementations.
29: Level: developer
31: .seealso: PEPSetUp()
32: @*/
33: PetscErrorCode PEPSetWorkVecs(PEP pep,PetscInt nw)
34: {
35: Vec t;
37: PetscFunctionBegin;
40: PetscCheck(nw>0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"nw must be > 0: nw = %" PetscInt_FMT,nw);
41: if (pep->nwork < nw) {
42: PetscCall(VecDestroyVecs(pep->nwork,&pep->work));
43: pep->nwork = nw;
44: PetscCall(BVGetColumn(pep->V,0,&t));
45: PetscCall(VecDuplicateVecs(t,nw,&pep->work));
46: PetscCall(BVRestoreColumn(pep->V,0,&t));
47: }
48: PetscFunctionReturn(PETSC_SUCCESS);
49: }
51: /*
52: PEPConvergedRelative - Checks convergence relative to the eigenvalue.
53: */
54: PetscErrorCode PEPConvergedRelative(PEP pep,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
55: {
56: PetscReal w;
58: PetscFunctionBegin;
59: w = SlepcAbsEigenvalue(eigr,eigi);
60: *errest = res/w;
61: PetscFunctionReturn(PETSC_SUCCESS);
62: }
64: /*
65: PEPConvergedNorm - Checks convergence relative to the matrix norms.
66: */
67: PetscErrorCode PEPConvergedNorm(PEP pep,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
68: {
69: PetscReal w=0.0,t;
70: PetscInt j;
71: PetscBool flg;
73: PetscFunctionBegin;
74: /* initialization of matrix norms */
75: if (!pep->nrma[pep->nmat-1]) {
76: for (j=0;j<pep->nmat;j++) {
77: PetscCall(MatHasOperation(pep->A[j],MATOP_NORM,&flg));
78: PetscCheck(flg,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONG,"The convergence test related to the matrix norms requires a matrix norm operation");
79: PetscCall(MatNorm(pep->A[j],NORM_INFINITY,&pep->nrma[j]));
80: }
81: }
82: t = SlepcAbsEigenvalue(eigr,eigi);
83: for (j=pep->nmat-1;j>=0;j--) {
84: w = w*t+pep->nrma[j];
85: }
86: *errest = res/w;
87: PetscFunctionReturn(PETSC_SUCCESS);
88: }
90: /*
91: PEPSetWhichEigenpairs_Default - Sets the default value for which,
92: depending on the ST.
93: */
94: PetscErrorCode PEPSetWhichEigenpairs_Default(PEP pep)
95: {
96: PetscBool target;
98: PetscFunctionBegin;
99: PetscCall(PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&target));
100: if (target) pep->which = PEP_TARGET_MAGNITUDE;
101: else pep->which = PEP_LARGEST_MAGNITUDE;
102: PetscFunctionReturn(PETSC_SUCCESS);
103: }
105: /*
106: PEPConvergedAbsolute - Checks convergence absolutely.
107: */
108: PetscErrorCode PEPConvergedAbsolute(PEP pep,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
109: {
110: PetscFunctionBegin;
111: *errest = res;
112: PetscFunctionReturn(PETSC_SUCCESS);
113: }
115: /*@C
116: PEPStoppingBasic - Default routine to determine whether the outer eigensolver
117: iteration must be stopped.
119: Collective
121: Input Parameters:
122: + pep - eigensolver context obtained from PEPCreate()
123: . its - current number of iterations
124: . max_it - maximum number of iterations
125: . nconv - number of currently converged eigenpairs
126: . nev - number of requested eigenpairs
127: - ctx - context (not used here)
129: Output Parameter:
130: . reason - result of the stopping test
132: Notes:
133: A positive value of reason indicates that the iteration has finished successfully
134: (converged), and a negative value indicates an error condition (diverged). If
135: the iteration needs to be continued, reason must be set to PEP_CONVERGED_ITERATING
136: (zero).
138: PEPStoppingBasic() will stop if all requested eigenvalues are converged, or if
139: the maximum number of iterations has been reached.
141: Use PEPSetStoppingTest() to provide your own test instead of using this one.
143: Level: advanced
145: .seealso: PEPSetStoppingTest(), PEPConvergedReason, PEPGetConvergedReason()
146: @*/
147: PetscErrorCode PEPStoppingBasic(PEP pep,PetscInt its,PetscInt max_it,PetscInt nconv,PetscInt nev,PEPConvergedReason *reason,void *ctx)
148: {
149: PetscFunctionBegin;
150: *reason = PEP_CONVERGED_ITERATING;
151: if (nconv >= nev) {
152: PetscCall(PetscInfo(pep,"Polynomial eigensolver finished successfully: %" PetscInt_FMT " eigenpairs converged at iteration %" PetscInt_FMT "\n",nconv,its));
153: *reason = PEP_CONVERGED_TOL;
154: } else if (its >= max_it) {
155: *reason = PEP_DIVERGED_ITS;
156: PetscCall(PetscInfo(pep,"Polynomial eigensolver iteration reached maximum number of iterations (%" PetscInt_FMT ")\n",its));
157: }
158: PetscFunctionReturn(PETSC_SUCCESS);
159: }
161: PetscErrorCode PEPBackTransform_Default(PEP pep)
162: {
163: PetscFunctionBegin;
164: PetscCall(STBackTransform(pep->st,pep->nconv,pep->eigr,pep->eigi));
165: PetscFunctionReturn(PETSC_SUCCESS);
166: }
168: PetscErrorCode PEPComputeVectors_Default(PEP pep)
169: {
170: PetscInt i;
171: Vec v;
173: PetscFunctionBegin;
174: PetscCall(PEPExtractVectors(pep));
176: /* Fix eigenvectors if balancing was used */
177: if ((pep->scale==PEP_SCALE_DIAGONAL || pep->scale==PEP_SCALE_BOTH) && pep->Dr && (pep->refine!=PEP_REFINE_MULTIPLE)) {
178: for (i=0;i<pep->nconv;i++) {
179: PetscCall(BVGetColumn(pep->V,i,&v));
180: PetscCall(VecPointwiseMult(v,v,pep->Dr));
181: PetscCall(BVRestoreColumn(pep->V,i,&v));
182: }
183: }
185: /* normalization */
186: PetscCall(BVNormalize(pep->V,pep->eigi));
187: PetscFunctionReturn(PETSC_SUCCESS);
188: }
190: /*
191: PEPBuildDiagonalScaling - compute two diagonal matrices to be applied for balancing
192: in polynomial eigenproblems.
193: */
194: PetscErrorCode PEPBuildDiagonalScaling(PEP pep)
195: {
196: PetscInt it,i,j,k,nmat,nr,e,nz,lst,lend,nc=0,*cols,emax,emin,emaxl,eminl;
197: const PetscInt *cidx,*ridx;
198: Mat M,*T,A;
199: PetscMPIInt n;
200: PetscBool cont=PETSC_TRUE,flg=PETSC_FALSE;
201: PetscScalar *array,*Dr,*Dl,t;
202: PetscReal l2,d,*rsum,*aux,*csum,w=1.0;
203: MatStructure str;
204: MatInfo info;
206: PetscFunctionBegin;
207: l2 = 2*PetscLogReal(2.0);
208: nmat = pep->nmat;
209: PetscCall(PetscMPIIntCast(pep->n,&n));
210: PetscCall(STGetMatStructure(pep->st,&str));
211: PetscCall(PetscMalloc1(nmat,&T));
212: for (k=0;k<nmat;k++) PetscCall(STGetMatrixTransformed(pep->st,k,&T[k]));
213: /* Form local auxiliary matrix M */
214: PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)T[0],&cont,MATMPIAIJ,MATSEQAIJ,""));
215: PetscCheck(cont,PetscObjectComm((PetscObject)T[0]),PETSC_ERR_SUP,"Only for MPIAIJ or SEQAIJ matrix types");
216: PetscCall(PetscObjectBaseTypeCompare((PetscObject)T[0],MATMPIAIJ,&cont));
217: if (cont) {
218: PetscCall(MatMPIAIJGetLocalMat(T[0],MAT_INITIAL_MATRIX,&M));
219: flg = PETSC_TRUE;
220: } else PetscCall(MatDuplicate(T[0],MAT_COPY_VALUES,&M));
221: PetscCall(MatGetInfo(M,MAT_LOCAL,&info));
222: nz = (PetscInt)info.nz_used;
223: PetscCall(MatSeqAIJGetArray(M,&array));
224: for (i=0;i<nz;i++) {
225: t = PetscAbsScalar(array[i]);
226: array[i] = t*t;
227: }
228: PetscCall(MatSeqAIJRestoreArray(M,&array));
229: for (k=1;k<nmat;k++) {
230: if (flg) PetscCall(MatMPIAIJGetLocalMat(T[k],MAT_INITIAL_MATRIX,&A));
231: else {
232: if (str==SAME_NONZERO_PATTERN) PetscCall(MatCopy(T[k],A,SAME_NONZERO_PATTERN));
233: else PetscCall(MatDuplicate(T[k],MAT_COPY_VALUES,&A));
234: }
235: PetscCall(MatGetInfo(A,MAT_LOCAL,&info));
236: nz = (PetscInt)info.nz_used;
237: PetscCall(MatSeqAIJGetArray(A,&array));
238: for (i=0;i<nz;i++) {
239: t = PetscAbsScalar(array[i]);
240: array[i] = t*t;
241: }
242: PetscCall(MatSeqAIJRestoreArray(A,&array));
243: w *= pep->slambda*pep->slambda*pep->sfactor;
244: PetscCall(MatAXPY(M,w,A,str));
245: if (flg || str!=SAME_NONZERO_PATTERN || k==nmat-2) PetscCall(MatDestroy(&A));
246: }
247: PetscCall(MatGetRowIJ(M,0,PETSC_FALSE,PETSC_FALSE,&nr,&ridx,&cidx,&cont));
248: PetscCheck(cont,PetscObjectComm((PetscObject)T[0]),PETSC_ERR_SUP,"It is not possible to compute scaling diagonals for these PEP matrices");
249: PetscCall(MatGetInfo(M,MAT_LOCAL,&info));
250: nz = (PetscInt)info.nz_used;
251: PetscCall(VecGetOwnershipRange(pep->Dl,&lst,&lend));
252: PetscCall(PetscMalloc4(nr,&rsum,pep->n,&csum,pep->n,&aux,PetscMin(pep->n-lend+lst,nz),&cols));
253: PetscCall(VecSet(pep->Dr,1.0));
254: PetscCall(VecSet(pep->Dl,1.0));
255: PetscCall(VecGetArray(pep->Dl,&Dl));
256: PetscCall(VecGetArray(pep->Dr,&Dr));
257: PetscCall(MatSeqAIJGetArray(M,&array));
258: PetscCall(PetscArrayzero(aux,pep->n));
259: for (j=0;j<nz;j++) {
260: /* Search non-zero columns outsize lst-lend */
261: if (aux[cidx[j]]==0 && (cidx[j]<lst || lend<=cidx[j])) cols[nc++] = cidx[j];
262: /* Local column sums */
263: aux[cidx[j]] += PetscAbsScalar(array[j]);
264: }
265: for (it=0;it<pep->sits && cont;it++) {
266: emaxl = 0; eminl = 0;
267: /* Column sum */
268: if (it>0) { /* it=0 has been already done*/
269: PetscCall(MatSeqAIJGetArray(M,&array));
270: PetscCall(PetscArrayzero(aux,pep->n));
271: for (j=0;j<nz;j++) aux[cidx[j]] += PetscAbsScalar(array[j]);
272: PetscCall(MatSeqAIJRestoreArray(M,&array));
273: }
274: PetscCall(MPIU_Allreduce(aux,csum,n,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)pep->Dr)));
275: /* Update Dr */
276: for (j=lst;j<lend;j++) {
277: d = PetscLogReal(csum[j])/l2;
278: e = -(PetscInt)((d < 0)?(d-0.5):(d+0.5));
279: d = PetscPowReal(2.0,e);
280: Dr[j-lst] *= d;
281: aux[j] = d*d;
282: emaxl = PetscMax(emaxl,e);
283: eminl = PetscMin(eminl,e);
284: }
285: for (j=0;j<nc;j++) {
286: d = PetscLogReal(csum[cols[j]])/l2;
287: e = -(PetscInt)((d < 0)?(d-0.5):(d+0.5));
288: d = PetscPowReal(2.0,e);
289: aux[cols[j]] = d*d;
290: emaxl = PetscMax(emaxl,e);
291: eminl = PetscMin(eminl,e);
292: }
293: /* Scale M */
294: PetscCall(MatSeqAIJGetArray(M,&array));
295: for (j=0;j<nz;j++) {
296: array[j] *= aux[cidx[j]];
297: }
298: PetscCall(MatSeqAIJRestoreArray(M,&array));
299: /* Row sum */
300: PetscCall(PetscArrayzero(rsum,nr));
301: PetscCall(MatSeqAIJGetArray(M,&array));
302: for (i=0;i<nr;i++) {
303: for (j=ridx[i];j<ridx[i+1];j++) rsum[i] += PetscAbsScalar(array[j]);
304: /* Update Dl */
305: d = PetscLogReal(rsum[i])/l2;
306: e = -(PetscInt)((d < 0)?(d-0.5):(d+0.5));
307: d = PetscPowReal(2.0,e);
308: Dl[i] *= d;
309: /* Scale M */
310: for (j=ridx[i];j<ridx[i+1];j++) array[j] *= d*d;
311: emaxl = PetscMax(emaxl,e);
312: eminl = PetscMin(eminl,e);
313: }
314: PetscCall(MatSeqAIJRestoreArray(M,&array));
315: /* Compute global max and min */
316: PetscCall(MPIU_Allreduce(&emaxl,&emax,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)pep->Dl)));
317: PetscCall(MPIU_Allreduce(&eminl,&emin,1,MPIU_INT,MPI_MIN,PetscObjectComm((PetscObject)pep->Dl)));
318: if (emax<=emin+2) cont = PETSC_FALSE;
319: }
320: PetscCall(VecRestoreArray(pep->Dr,&Dr));
321: PetscCall(VecRestoreArray(pep->Dl,&Dl));
322: /* Free memory*/
323: PetscCall(MatDestroy(&M));
324: PetscCall(PetscFree4(rsum,csum,aux,cols));
325: PetscCall(PetscFree(T));
326: PetscFunctionReturn(PETSC_SUCCESS);
327: }
329: /*
330: PEPComputeScaleFactor - compute sfactor as described in [Betcke 2008].
331: */
332: PetscErrorCode PEPComputeScaleFactor(PEP pep)
333: {
334: PetscBool has0,has1,flg;
335: PetscReal norm0,norm1;
336: Mat T[2];
337: PEPBasis basis;
338: PetscInt i;
340: PetscFunctionBegin;
341: if (pep->scale==PEP_SCALE_NONE || pep->scale==PEP_SCALE_DIAGONAL) { /* no scalar scaling */
342: pep->sfactor = 1.0;
343: pep->dsfactor = 1.0;
344: PetscFunctionReturn(PETSC_SUCCESS);
345: }
346: if (pep->sfactor_set) PetscFunctionReturn(PETSC_SUCCESS); /* user provided value */
347: pep->sfactor = 1.0;
348: pep->dsfactor = 1.0;
349: PetscCall(PEPGetBasis(pep,&basis));
350: if (basis==PEP_BASIS_MONOMIAL) {
351: PetscCall(STGetTransform(pep->st,&flg));
352: if (flg) {
353: PetscCall(STGetMatrixTransformed(pep->st,0,&T[0]));
354: PetscCall(STGetMatrixTransformed(pep->st,pep->nmat-1,&T[1]));
355: } else {
356: T[0] = pep->A[0];
357: T[1] = pep->A[pep->nmat-1];
358: }
359: if (pep->nmat>2) {
360: PetscCall(MatHasOperation(T[0],MATOP_NORM,&has0));
361: PetscCall(MatHasOperation(T[1],MATOP_NORM,&has1));
362: if (has0 && has1) {
363: PetscCall(MatNorm(T[0],NORM_INFINITY,&norm0));
364: PetscCall(MatNorm(T[1],NORM_INFINITY,&norm1));
365: pep->sfactor = PetscPowReal(norm0/norm1,1.0/(pep->nmat-1));
366: pep->dsfactor = norm1;
367: for (i=pep->nmat-2;i>0;i--) {
368: PetscCall(STGetMatrixTransformed(pep->st,i,&T[1]));
369: PetscCall(MatHasOperation(T[1],MATOP_NORM,&has1));
370: if (has1) {
371: PetscCall(MatNorm(T[1],NORM_INFINITY,&norm1));
372: pep->dsfactor = pep->dsfactor*pep->sfactor+norm1;
373: } else break;
374: }
375: if (has1) {
376: pep->dsfactor = pep->dsfactor*pep->sfactor+norm0;
377: pep->dsfactor = pep->nmat/pep->dsfactor;
378: } else pep->dsfactor = 1.0;
379: }
380: }
381: }
382: PetscFunctionReturn(PETSC_SUCCESS);
383: }
385: /*
386: PEPBasisCoefficients - compute polynomial basis coefficients
387: */
388: PetscErrorCode PEPBasisCoefficients(PEP pep,PetscReal *pbc)
389: {
390: PetscReal *ca,*cb,*cg;
391: PetscInt k,nmat=pep->nmat;
393: PetscFunctionBegin;
394: ca = pbc;
395: cb = pbc+nmat;
396: cg = pbc+2*nmat;
397: switch (pep->basis) {
398: case PEP_BASIS_MONOMIAL:
399: for (k=0;k<nmat;k++) {
400: ca[k] = 1.0; cb[k] = 0.0; cg[k] = 0.0;
401: }
402: break;
403: case PEP_BASIS_CHEBYSHEV1:
404: ca[0] = 1.0; cb[0] = 0.0; cg[0] = 0.0;
405: for (k=1;k<nmat;k++) {
406: ca[k] = .5; cb[k] = 0.0; cg[k] = .5;
407: }
408: break;
409: case PEP_BASIS_CHEBYSHEV2:
410: ca[0] = .5; cb[0] = 0.0; cg[0] = 0.0;
411: for (k=1;k<nmat;k++) {
412: ca[k] = .5; cb[k] = 0.0; cg[k] = .5;
413: }
414: break;
415: case PEP_BASIS_LEGENDRE:
416: ca[0] = 1.0; cb[0] = 0.0; cg[0] = 0.0;
417: for (k=1;k<nmat;k++) {
418: ca[k] = k+1; cb[k] = -2*k; cg[k] = k;
419: }
420: break;
421: case PEP_BASIS_LAGUERRE:
422: ca[0] = -1.0; cb[0] = 0.0; cg[0] = 0.0;
423: for (k=1;k<nmat;k++) {
424: ca[k] = -(k+1); cb[k] = 2*k+1; cg[k] = -k;
425: }
426: break;
427: case PEP_BASIS_HERMITE:
428: ca[0] = .5; cb[0] = 0.0; cg[0] = 0.0;
429: for (k=1;k<nmat;k++) {
430: ca[k] = .5; cb[k] = 0.0; cg[k] = -k;
431: }
432: break;
433: }
434: PetscFunctionReturn(PETSC_SUCCESS);
435: }