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