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