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