Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
6 : This file is part of SLEPc.
7 : SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9 : */
10 :
11 : #include <slepc/private/slepccontour.h>
12 : #include <slepcblaslapack.h>
13 :
14 : /*
15 : SlepcContourDataCreate - Create a contour data structure.
16 :
17 : Input Parameters:
18 : n - the number of integration points
19 : npart - number of partitions for the subcommunicator
20 : parent - parent object
21 : */
22 26 : PetscErrorCode SlepcContourDataCreate(PetscInt n,PetscInt npart,PetscObject parent,SlepcContourData *contour)
23 : {
24 26 : PetscFunctionBegin;
25 26 : PetscCall(PetscNew(contour));
26 26 : (*contour)->parent = parent;
27 26 : PetscCall(PetscSubcommCreate(PetscObjectComm(parent),&(*contour)->subcomm));
28 26 : PetscCall(PetscSubcommSetNumber((*contour)->subcomm,npart));
29 26 : PetscCall(PetscSubcommSetType((*contour)->subcomm,PETSC_SUBCOMM_INTERLACED));
30 26 : (*contour)->npoints = n / npart;
31 26 : if (n%npart > (*contour)->subcomm->color) (*contour)->npoints++;
32 26 : PetscFunctionReturn(PETSC_SUCCESS);
33 : }
34 :
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 17 : PetscErrorCode SlepcContourDataReset(SlepcContourData contour)
40 : {
41 17 : PetscInt i;
42 :
43 17 : PetscFunctionBegin;
44 17 : if (contour->ksp) {
45 389 : for (i=0;i<contour->npoints;i++) PetscCall(KSPReset(contour->ksp[i]));
46 : }
47 17 : if (contour->pA) {
48 8 : PetscCall(MatDestroyMatrices(contour->nmat,&contour->pA));
49 8 : PetscCall(MatDestroyMatrices(contour->nmat,&contour->pP));
50 8 : contour->nmat = 0;
51 : }
52 17 : PetscCall(VecScatterDestroy(&contour->scatterin));
53 17 : PetscCall(VecDestroy(&contour->xsub));
54 17 : PetscCall(VecDestroy(&contour->xdup));
55 17 : PetscFunctionReturn(PETSC_SUCCESS);
56 : }
57 :
58 : /*
59 : SlepcContourDataDestroy - Destroys the contour data structure.
60 : */
61 38 : PetscErrorCode SlepcContourDataDestroy(SlepcContourData *contour)
62 : {
63 38 : PetscInt i;
64 :
65 38 : PetscFunctionBegin;
66 38 : if (!(*contour)) PetscFunctionReturn(PETSC_SUCCESS);
67 26 : if ((*contour)->ksp) {
68 678 : for (i=0;i<(*contour)->npoints;i++) PetscCall(KSPDestroy(&(*contour)->ksp[i]));
69 26 : PetscCall(PetscFree((*contour)->ksp));
70 : }
71 26 : PetscCall(PetscSubcommDestroy(&(*contour)->subcomm));
72 26 : PetscCall(PetscFree((*contour)));
73 26 : *contour = NULL;
74 26 : PetscFunctionReturn(PETSC_SUCCESS);
75 : }
76 :
77 : /*
78 : SlepcContourRedundantMat - Creates redundant copies of the passed matrices in the subcomm.
79 :
80 : Input Parameters:
81 : nmat - the number of matrices
82 : A - array of matrices
83 : P - array of matrices (preconditioner)
84 : */
85 26 : PetscErrorCode SlepcContourRedundantMat(SlepcContourData contour,PetscInt nmat,Mat *A,Mat *P)
86 : {
87 26 : PetscInt i;
88 26 : MPI_Comm child;
89 :
90 26 : PetscFunctionBegin;
91 26 : if (contour->pA) {
92 0 : PetscCall(MatDestroyMatrices(contour->nmat,&contour->pA));
93 0 : PetscCall(MatDestroyMatrices(contour->nmat,&contour->pP));
94 0 : contour->nmat = 0;
95 : }
96 26 : if (contour->subcomm && contour->subcomm->n != 1) {
97 8 : PetscCall(PetscSubcommGetChild(contour->subcomm,&child));
98 8 : PetscCall(PetscCalloc1(nmat,&contour->pA));
99 20 : for (i=0;i<nmat;i++) PetscCall(MatCreateRedundantMatrix(A[i],contour->subcomm->n,child,MAT_INITIAL_MATRIX,&contour->pA[i]));
100 8 : if (P) {
101 2 : PetscCall(PetscCalloc1(nmat,&contour->pP));
102 6 : for (i=0;i<nmat;i++) PetscCall(MatCreateRedundantMatrix(P[i],contour->subcomm->n,child,MAT_INITIAL_MATRIX,&contour->pP[i]));
103 : }
104 8 : contour->nmat = nmat;
105 : }
106 26 : PetscFunctionReturn(PETSC_SUCCESS);
107 : }
108 :
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).
114 :
115 : Input Parameters:
116 : v - the regular vector from which dimensions are taken
117 : */
118 8 : PetscErrorCode SlepcContourScatterCreate(SlepcContourData contour,Vec v)
119 : {
120 8 : IS is1,is2;
121 8 : PetscInt i,j,k,m,mstart,mend,mlocal;
122 8 : PetscInt *idx1,*idx2,mloc_sub;
123 8 : MPI_Comm contpar,parent;
124 :
125 8 : PetscFunctionBegin;
126 8 : PetscCall(VecDestroy(&contour->xsub));
127 8 : PetscCall(MatCreateVecsEmpty(contour->pA[0],&contour->xsub,NULL));
128 :
129 8 : PetscCall(VecDestroy(&contour->xdup));
130 8 : PetscCall(MatGetLocalSize(contour->pA[0],&mloc_sub,NULL));
131 8 : PetscCall(PetscSubcommGetContiguousParent(contour->subcomm,&contpar));
132 8 : PetscCall(VecCreate(contpar,&contour->xdup));
133 8 : PetscCall(VecSetSizes(contour->xdup,mloc_sub,PETSC_DECIDE));
134 8 : PetscCall(VecSetType(contour->xdup,((PetscObject)v)->type_name));
135 :
136 8 : PetscCall(VecScatterDestroy(&contour->scatterin));
137 8 : PetscCall(VecGetSize(v,&m));
138 8 : PetscCall(VecGetOwnershipRange(v,&mstart,&mend));
139 8 : mlocal = mend - mstart;
140 8 : PetscCall(PetscMalloc2(contour->subcomm->n*mlocal,&idx1,contour->subcomm->n*mlocal,&idx2));
141 : j = 0;
142 24 : for (k=0;k<contour->subcomm->n;k++) {
143 4968 : for (i=mstart;i<mend;i++) {
144 4952 : idx1[j] = i;
145 4952 : idx2[j++] = i + m*k;
146 : }
147 : }
148 8 : PetscCall(PetscSubcommGetParent(contour->subcomm,&parent));
149 8 : PetscCall(ISCreateGeneral(parent,contour->subcomm->n*mlocal,idx1,PETSC_COPY_VALUES,&is1));
150 8 : PetscCall(ISCreateGeneral(parent,contour->subcomm->n*mlocal,idx2,PETSC_COPY_VALUES,&is2));
151 8 : PetscCall(VecScatterCreate(v,is1,contour->xdup,is2,&contour->scatterin));
152 8 : PetscCall(ISDestroy(&is1));
153 8 : PetscCall(ISDestroy(&is2));
154 8 : PetscCall(PetscFree2(idx1,idx2));
155 8 : PetscFunctionReturn(PETSC_SUCCESS);
156 : }
157 :
158 : /*
159 : SlepcCISS_isGhost - Determine if any of the computed eigenpairs are spurious.
160 :
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
166 :
167 : Output Parameter:
168 : fl - array of n booleans
169 : */
170 26 : PetscErrorCode SlepcCISS_isGhost(Mat X,PetscInt n,PetscReal *sigma,PetscReal thresh,PetscBool *fl)
171 : {
172 26 : const PetscScalar *pX;
173 26 : PetscInt i,j,m,ld;
174 26 : PetscReal *tau,s1,s2,tau_max=0.0;
175 :
176 26 : PetscFunctionBegin;
177 26 : PetscCall(MatGetSize(X,&m,NULL));
178 26 : PetscCall(MatDenseGetLDA(X,&ld));
179 26 : PetscCall(PetscMalloc1(n,&tau));
180 26 : PetscCall(MatDenseGetArrayRead(X,&pX));
181 508 : for (j=0;j<n;j++) {
182 : s1 = 0.0;
183 : s2 = 0.0;
184 22864 : for (i=0;i<m;i++) {
185 22382 : s1 += PetscAbsScalar(PetscPowScalarInt(pX[i+j*ld],2));
186 44764 : s2 += PetscPowRealInt(PetscAbsScalar(pX[i+j*ld]),2)/sigma[i];
187 : }
188 482 : tau[j] = s1/s2;
189 596 : tau_max = PetscMax(tau_max,tau[j]);
190 : }
191 26 : PetscCall(MatDenseRestoreArrayRead(X,&pX));
192 508 : for (j=0;j<n;j++) fl[j] = (tau[j]>=thresh*tau_max)? PETSC_TRUE: PETSC_FALSE;
193 26 : PetscCall(PetscFree(tau));
194 26 : PetscFunctionReturn(PETSC_SUCCESS);
195 : }
196 :
197 : /*
198 : SlepcCISS_BH_SVD - Compute SVD of block Hankel matrix and its rank.
199 :
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
204 :
205 : Output Parameters:
206 : sigma - computed singular values
207 : rank - the rank of H
208 : */
209 6 : PetscErrorCode SlepcCISS_BH_SVD(PetscScalar *H,PetscInt ml,PetscReal delta,PetscReal *sigma,PetscInt *rank)
210 : {
211 6 : PetscInt i;
212 6 : PetscBLASInt m,n,lda,ldu,ldvt,lwork,info;
213 6 : PetscScalar *work;
214 : #if defined(PETSC_USE_COMPLEX)
215 : PetscReal *rwork;
216 : #endif
217 :
218 6 : PetscFunctionBegin;
219 6 : PetscCall(PetscMalloc1(5*ml,&work));
220 : #if defined(PETSC_USE_COMPLEX)
221 : PetscCall(PetscMalloc1(5*ml,&rwork));
222 : #endif
223 6 : PetscCall(PetscBLASIntCast(ml,&m));
224 6 : n = m; lda = m; ldu = m; ldvt = m; lwork = 5*m;
225 6 : 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 6 : PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","N",&m,&n,H,&lda,sigma,NULL,&ldu,NULL,&ldvt,work,&lwork,&info));
230 : #endif
231 6 : SlepcCheckLapackInfo("gesvd",info);
232 6 : PetscCall(PetscFPTrapPop());
233 6 : (*rank) = 0;
234 466 : for (i=0;i<ml;i++) {
235 864 : if (sigma[i]/PetscMax(sigma[0],1.0)>delta) (*rank)++;
236 : }
237 6 : PetscCall(PetscFree(work));
238 : #if defined(PETSC_USE_COMPLEX)
239 : PetscCall(PetscFree(rwork));
240 : #endif
241 6 : PetscFunctionReturn(PETSC_SUCCESS);
242 : }
|
