Actual source code: slepccontour.c
slepc-3.21.1 2024-04-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: */
11: #include <slepc/private/slepccontour.h>
12: #include <slepcblaslapack.h>
14: /*
15: SlepcContourDataCreate - Create a contour data structure.
17: Input Parameters:
18: n - the number of integration points
19: npart - number of partitions for the subcommunicator
20: parent - parent object
21: */
22: PetscErrorCode SlepcContourDataCreate(PetscInt n,PetscInt npart,PetscObject parent,SlepcContourData *contour)
23: {
24: PetscFunctionBegin;
25: PetscCall(PetscNew(contour));
26: (*contour)->parent = parent;
27: PetscCall(PetscSubcommCreate(PetscObjectComm(parent),&(*contour)->subcomm));
28: PetscCall(PetscSubcommSetNumber((*contour)->subcomm,npart));
29: PetscCall(PetscSubcommSetType((*contour)->subcomm,PETSC_SUBCOMM_INTERLACED));
30: (*contour)->npoints = n / npart;
31: if (n%npart > (*contour)->subcomm->color) (*contour)->npoints++;
32: PetscFunctionReturn(PETSC_SUCCESS);
33: }
35: /*
36: SlepcContourDataReset - Resets the KSP objects in a contour data structure,
37: and destroys any objects whose size depends on the problem size.
38: */
39: PetscErrorCode SlepcContourDataReset(SlepcContourData contour)
40: {
41: PetscInt i;
43: PetscFunctionBegin;
44: if (contour->ksp) {
45: for (i=0;i<contour->npoints;i++) PetscCall(KSPReset(contour->ksp[i]));
46: }
47: if (contour->pA) {
48: PetscCall(MatDestroyMatrices(contour->nmat,&contour->pA));
49: PetscCall(MatDestroyMatrices(contour->nmat,&contour->pP));
50: contour->nmat = 0;
51: }
52: PetscCall(VecScatterDestroy(&contour->scatterin));
53: PetscCall(VecDestroy(&contour->xsub));
54: PetscCall(VecDestroy(&contour->xdup));
55: PetscFunctionReturn(PETSC_SUCCESS);
56: }
58: /*
59: SlepcContourDataDestroy - Destroys the contour data structure.
60: */
61: PetscErrorCode SlepcContourDataDestroy(SlepcContourData *contour)
62: {
63: PetscInt i;
65: PetscFunctionBegin;
66: if (!(*contour)) PetscFunctionReturn(PETSC_SUCCESS);
67: if ((*contour)->ksp) {
68: for (i=0;i<(*contour)->npoints;i++) PetscCall(KSPDestroy(&(*contour)->ksp[i]));
69: PetscCall(PetscFree((*contour)->ksp));
70: }
71: PetscCall(PetscSubcommDestroy(&(*contour)->subcomm));
72: PetscCall(PetscFree((*contour)));
73: *contour = NULL;
74: PetscFunctionReturn(PETSC_SUCCESS);
75: }
77: /*
78: SlepcContourRedundantMat - Creates redundant copies of the passed matrices in the subcomm.
80: Input Parameters:
81: nmat - the number of matrices
82: A - array of matrices
83: P - array of matrices (preconditioner)
84: */
85: PetscErrorCode SlepcContourRedundantMat(SlepcContourData contour,PetscInt nmat,Mat *A,Mat *P)
86: {
87: PetscInt i;
88: MPI_Comm child;
90: PetscFunctionBegin;
91: if (contour->pA) {
92: PetscCall(MatDestroyMatrices(contour->nmat,&contour->pA));
93: PetscCall(MatDestroyMatrices(contour->nmat,&contour->pP));
94: contour->nmat = 0;
95: }
96: if (contour->subcomm && contour->subcomm->n != 1) {
97: PetscCall(PetscSubcommGetChild(contour->subcomm,&child));
98: PetscCall(PetscCalloc1(nmat,&contour->pA));
99: for (i=0;i<nmat;i++) PetscCall(MatCreateRedundantMatrix(A[i],contour->subcomm->n,child,MAT_INITIAL_MATRIX,&contour->pA[i]));
100: if (P) {
101: PetscCall(PetscCalloc1(nmat,&contour->pP));
102: for (i=0;i<nmat;i++) PetscCall(MatCreateRedundantMatrix(P[i],contour->subcomm->n,child,MAT_INITIAL_MATRIX,&contour->pP[i]));
103: }
104: contour->nmat = nmat;
105: }
106: PetscFunctionReturn(PETSC_SUCCESS);
107: }
109: /*
110: SlepcContourScatterCreate - Creates a scatter context to communicate between a
111: regular vector and a vector xdup that can hold one duplicate per each subcommunicator
112: on the contiguous parent communicator. Also creates auxiliary vectors xdup and xsub
113: (the latter with the same layout as the redundant matrices in the subcommunicator).
115: Input Parameters:
116: v - the regular vector from which dimensions are taken
117: */
118: PetscErrorCode SlepcContourScatterCreate(SlepcContourData contour,Vec v)
119: {
120: IS is1,is2;
121: PetscInt i,j,k,m,mstart,mend,mlocal;
122: PetscInt *idx1,*idx2,mloc_sub;
123: MPI_Comm contpar,parent;
125: PetscFunctionBegin;
126: PetscCall(VecDestroy(&contour->xsub));
127: PetscCall(MatCreateVecsEmpty(contour->pA[0],&contour->xsub,NULL));
129: PetscCall(VecDestroy(&contour->xdup));
130: PetscCall(MatGetLocalSize(contour->pA[0],&mloc_sub,NULL));
131: PetscCall(PetscSubcommGetContiguousParent(contour->subcomm,&contpar));
132: PetscCall(VecCreate(contpar,&contour->xdup));
133: PetscCall(VecSetSizes(contour->xdup,mloc_sub,PETSC_DECIDE));
134: PetscCall(VecSetType(contour->xdup,((PetscObject)v)->type_name));
136: PetscCall(VecScatterDestroy(&contour->scatterin));
137: PetscCall(VecGetSize(v,&m));
138: PetscCall(VecGetOwnershipRange(v,&mstart,&mend));
139: mlocal = mend - mstart;
140: PetscCall(PetscMalloc2(contour->subcomm->n*mlocal,&idx1,contour->subcomm->n*mlocal,&idx2));
141: j = 0;
142: for (k=0;k<contour->subcomm->n;k++) {
143: for (i=mstart;i<mend;i++) {
144: idx1[j] = i;
145: idx2[j++] = i + m*k;
146: }
147: }
148: PetscCall(PetscSubcommGetParent(contour->subcomm,&parent));
149: PetscCall(ISCreateGeneral(parent,contour->subcomm->n*mlocal,idx1,PETSC_COPY_VALUES,&is1));
150: PetscCall(ISCreateGeneral(parent,contour->subcomm->n*mlocal,idx2,PETSC_COPY_VALUES,&is2));
151: PetscCall(VecScatterCreate(v,is1,contour->xdup,is2,&contour->scatterin));
152: PetscCall(ISDestroy(&is1));
153: PetscCall(ISDestroy(&is2));
154: PetscCall(PetscFree2(idx1,idx2));
155: PetscFunctionReturn(PETSC_SUCCESS);
156: }
158: /*
159: SlepcCISS_isGhost - Determine if any of the computed eigenpairs are spurious.
161: Input Parameters:
162: X - the matrix of eigenvectors (MATSEQDENSE)
163: n - the number of columns to consider
164: sigma - the singular values
165: thresh - threshold to decide whether a value is spurious
167: Output Parameter:
168: fl - array of n booleans
169: */
170: PetscErrorCode SlepcCISS_isGhost(Mat X,PetscInt n,PetscReal *sigma,PetscReal thresh,PetscBool *fl)
171: {
172: const PetscScalar *pX;
173: PetscInt i,j,m,ld;
174: PetscReal *tau,s1,s2,tau_max=0.0;
176: PetscFunctionBegin;
177: PetscCall(MatGetSize(X,&m,NULL));
178: PetscCall(MatDenseGetLDA(X,&ld));
179: PetscCall(PetscMalloc1(n,&tau));
180: PetscCall(MatDenseGetArrayRead(X,&pX));
181: for (j=0;j<n;j++) {
182: s1 = 0.0;
183: s2 = 0.0;
184: for (i=0;i<m;i++) {
185: s1 += PetscAbsScalar(PetscPowScalarInt(pX[i+j*ld],2));
186: s2 += PetscPowRealInt(PetscAbsScalar(pX[i+j*ld]),2)/sigma[i];
187: }
188: tau[j] = s1/s2;
189: tau_max = PetscMax(tau_max,tau[j]);
190: }
191: PetscCall(MatDenseRestoreArrayRead(X,&pX));
192: for (j=0;j<n;j++) fl[j] = (tau[j]>=thresh*tau_max)? PETSC_TRUE: PETSC_FALSE;
193: PetscCall(PetscFree(tau));
194: PetscFunctionReturn(PETSC_SUCCESS);
195: }
197: /*
198: SlepcCISS_BH_SVD - Compute SVD of block Hankel matrix and its rank.
200: Input Parameters:
201: H - block Hankel matrix obtained via CISS_BlockHankel()
202: ml - dimension of rows and columns, equal to M*L
203: delta - the tolerance used to determine the rank
205: Output Parameters:
206: sigma - computed singular values
207: rank - the rank of H
208: */
209: PetscErrorCode SlepcCISS_BH_SVD(PetscScalar *H,PetscInt ml,PetscReal delta,PetscReal *sigma,PetscInt *rank)
210: {
211: PetscInt i;
212: PetscBLASInt m,n,lda,ldu,ldvt,lwork,info;
213: PetscScalar *work;
214: #if defined(PETSC_USE_COMPLEX)
215: PetscReal *rwork;
216: #endif
218: PetscFunctionBegin;
219: PetscCall(PetscMalloc1(5*ml,&work));
220: #if defined(PETSC_USE_COMPLEX)
221: PetscCall(PetscMalloc1(5*ml,&rwork));
222: #endif
223: PetscCall(PetscBLASIntCast(ml,&m));
224: n = m; lda = m; ldu = m; ldvt = m; lwork = 5*m;
225: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
226: #if defined(PETSC_USE_COMPLEX)
227: PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","N",&m,&n,H,&lda,sigma,NULL,&ldu,NULL,&ldvt,work,&lwork,rwork,&info));
228: #else
229: PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","N",&m,&n,H,&lda,sigma,NULL,&ldu,NULL,&ldvt,work,&lwork,&info));
230: #endif
231: SlepcCheckLapackInfo("gesvd",info);
232: PetscCall(PetscFPTrapPop());
233: (*rank) = 0;
234: for (i=0;i<ml;i++) {
235: if (sigma[i]/PetscMax(sigma[0],1.0)>delta) (*rank)++;
236: }
237: PetscCall(PetscFree(work));
238: #if defined(PETSC_USE_COMPLEX)
239: PetscCall(PetscFree(rwork));
240: #endif
241: PetscFunctionReturn(PETSC_SUCCESS);
242: }