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 90 : PetscErrorCode SlepcContourDataCreate(PetscInt n,PetscInt npart,PetscObject parent,SlepcContourData *contour)
23 : {
24 90 : PetscFunctionBegin;
25 90 : PetscCall(PetscNew(contour));
26 90 : (*contour)->parent = parent;
27 90 : PetscCall(PetscSubcommCreate(PetscObjectComm(parent),&(*contour)->subcomm));
28 90 : PetscCall(PetscSubcommSetNumber((*contour)->subcomm,npart));
29 90 : PetscCall(PetscSubcommSetType((*contour)->subcomm,PETSC_SUBCOMM_INTERLACED));
30 90 : (*contour)->npoints = n / npart;
31 90 : if (n%npart > (*contour)->subcomm->color) (*contour)->npoints++;
32 90 : 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 71 : PetscErrorCode SlepcContourDataReset(SlepcContourData contour)
40 : {
41 71 : PetscInt i;
42 :
43 71 : PetscFunctionBegin;
44 71 : if (contour->ksp) {
45 1611 : for (i=0;i<contour->npoints;i++) PetscCall(KSPReset(contour->ksp[i]));
46 : }
47 71 : if (contour->pA) {
48 30 : PetscCall(MatDestroyMatrices(contour->nmat,&contour->pA));
49 30 : PetscCall(MatDestroyMatrices(contour->nmat,&contour->pP));
50 30 : contour->nmat = 0;
51 : }
52 71 : PetscCall(VecScatterDestroy(&contour->scatterin));
53 71 : PetscCall(VecDestroy(&contour->xsub));
54 71 : PetscCall(VecDestroy(&contour->xdup));
55 71 : PetscFunctionReturn(PETSC_SUCCESS);
56 : }
57 :
58 : /*
59 : SlepcContourDataDestroy - Destroys the contour data structure.
60 : */
61 135 : PetscErrorCode SlepcContourDataDestroy(SlepcContourData *contour)
62 : {
63 135 : PetscInt i;
64 :
65 135 : PetscFunctionBegin;
66 135 : if (!(*contour)) PetscFunctionReturn(PETSC_SUCCESS);
67 90 : if ((*contour)->ksp) {
68 2102 : for (i=0;i<(*contour)->npoints;i++) PetscCall(KSPDestroy(&(*contour)->ksp[i]));
69 90 : PetscCall(PetscFree((*contour)->ksp));
70 : }
71 90 : PetscCall(PetscSubcommDestroy(&(*contour)->subcomm));
72 90 : PetscCall(PetscFree((*contour)));
73 90 : *contour = NULL;
74 90 : 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 86 : PetscErrorCode SlepcContourRedundantMat(SlepcContourData contour,PetscInt nmat,Mat *A,Mat *P)
86 : {
87 86 : PetscInt i;
88 86 : MPI_Comm child;
89 :
90 86 : PetscFunctionBegin;
91 86 : 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 86 : if (contour->subcomm && contour->subcomm->n != 1) {
97 30 : PetscCall(PetscSubcommGetChild(contour->subcomm,&child));
98 30 : PetscCall(PetscCalloc1(nmat,&contour->pA));
99 96 : for (i=0;i<nmat;i++) PetscCall(MatCreateRedundantMatrix(A[i],contour->subcomm->n,child,MAT_INITIAL_MATRIX,&contour->pA[i]));
100 30 : if (P) {
101 6 : PetscCall(PetscCalloc1(nmat,&contour->pP));
102 22 : for (i=0;i<nmat;i++) PetscCall(MatCreateRedundantMatrix(P[i],contour->subcomm->n,child,MAT_INITIAL_MATRIX,&contour->pP[i]));
103 : }
104 30 : contour->nmat = nmat;
105 : }
106 86 : 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 30 : PetscErrorCode SlepcContourScatterCreate(SlepcContourData contour,Vec v)
119 : {
120 30 : IS is1,is2;
121 30 : PetscInt i,j,k,m,mstart,mend,mlocal;
122 30 : PetscInt *idx1,*idx2,mloc_sub;
123 30 : MPI_Comm contpar,parent;
124 :
125 30 : PetscFunctionBegin;
126 30 : PetscCall(VecDestroy(&contour->xsub));
127 30 : PetscCall(MatCreateVecsEmpty(contour->pA[0],&contour->xsub,NULL));
128 :
129 30 : PetscCall(VecDestroy(&contour->xdup));
130 30 : PetscCall(MatGetLocalSize(contour->pA[0],&mloc_sub,NULL));
131 30 : PetscCall(PetscSubcommGetContiguousParent(contour->subcomm,&contpar));
132 30 : PetscCall(VecCreate(contpar,&contour->xdup));
133 30 : PetscCall(VecSetSizes(contour->xdup,mloc_sub,PETSC_DECIDE));
134 30 : PetscCall(VecSetType(contour->xdup,((PetscObject)v)->type_name));
135 :
136 30 : PetscCall(VecScatterDestroy(&contour->scatterin));
137 30 : PetscCall(VecGetSize(v,&m));
138 30 : PetscCall(VecGetOwnershipRange(v,&mstart,&mend));
139 30 : mlocal = mend - mstart;
140 30 : PetscCall(PetscMalloc2(contour->subcomm->n*mlocal,&idx1,contour->subcomm->n*mlocal,&idx2));
141 : j = 0;
142 98 : for (k=0;k<contour->subcomm->n;k++) {
143 7460 : for (i=mstart;i<mend;i++) {
144 7392 : idx1[j] = i;
145 7392 : idx2[j++] = i + m*k;
146 : }
147 : }
148 30 : PetscCall(PetscSubcommGetParent(contour->subcomm,&parent));
149 30 : PetscCall(ISCreateGeneral(parent,contour->subcomm->n*mlocal,idx1,PETSC_COPY_VALUES,&is1));
150 30 : PetscCall(ISCreateGeneral(parent,contour->subcomm->n*mlocal,idx2,PETSC_COPY_VALUES,&is2));
151 30 : PetscCall(VecScatterCreate(v,is1,contour->xdup,is2,&contour->scatterin));
152 30 : PetscCall(ISDestroy(&is1));
153 30 : PetscCall(ISDestroy(&is2));
154 30 : PetscCall(PetscFree2(idx1,idx2));
155 30 : 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 88 : PetscErrorCode SlepcCISS_isGhost(Mat X,PetscInt n,PetscReal *sigma,PetscReal thresh,PetscBool *fl)
171 : {
172 88 : const PetscScalar *pX;
173 88 : PetscInt i,j,m,ld;
174 88 : PetscReal *tau,s1,s2,tau_max=0.0;
175 :
176 88 : PetscFunctionBegin;
177 88 : PetscCall(MatGetSize(X,&m,NULL));
178 88 : PetscCall(MatDenseGetLDA(X,&ld));
179 88 : PetscCall(PetscMalloc1(n,&tau));
180 88 : PetscCall(MatDenseGetArrayRead(X,&pX));
181 2829 : for (j=0;j<n;j++) {
182 : s1 = 0.0;
183 : s2 = 0.0;
184 157342 : for (i=0;i<m;i++) {
185 154601 : s1 += PetscAbsScalar(PetscPowScalarInt(pX[i+j*ld],2));
186 309202 : s2 += PetscPowRealInt(PetscAbsScalar(pX[i+j*ld]),2)/sigma[i];
187 : }
188 2741 : tau[j] = s1/s2;
189 3102 : tau_max = PetscMax(tau_max,tau[j]);
190 : }
191 88 : PetscCall(MatDenseRestoreArrayRead(X,&pX));
192 2829 : for (j=0;j<n;j++) fl[j] = (tau[j]>=thresh*tau_max)? PETSC_TRUE: PETSC_FALSE;
193 88 : PetscCall(PetscFree(tau));
194 88 : 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 26 : PetscErrorCode SlepcCISS_BH_SVD(PetscScalar *H,PetscInt ml,PetscReal delta,PetscReal *sigma,PetscInt *rank)
210 : {
211 26 : PetscInt i;
212 26 : PetscBLASInt m,n,lda,ldu,ldvt,lwork,info;
213 26 : PetscScalar *work;
214 : #if defined(PETSC_USE_COMPLEX)
215 26 : PetscReal *rwork;
216 : #endif
217 :
218 26 : PetscFunctionBegin;
219 26 : PetscCall(PetscMalloc1(5*ml,&work));
220 : #if defined(PETSC_USE_COMPLEX)
221 26 : PetscCall(PetscMalloc1(5*ml,&rwork));
222 : #endif
223 26 : PetscCall(PetscBLASIntCast(ml,&m));
224 26 : n = m; lda = m; ldu = m; ldvt = m; lwork = 5*m;
225 26 : PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
226 : #if defined(PETSC_USE_COMPLEX)
227 26 : 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 26 : SlepcCheckLapackInfo("gesvd",info);
232 26 : PetscCall(PetscFPTrapPop());
233 26 : (*rank) = 0;
234 1534 : for (i=0;i<ml;i++) {
235 3016 : if (sigma[i]/PetscMax(sigma[0],1.0)>delta) (*rank)++;
236 : }
237 26 : PetscCall(PetscFree(work));
238 : #if defined(PETSC_USE_COMPLEX)
239 26 : PetscCall(PetscFree(rwork));
240 : #endif
241 26 : PetscFunctionReturn(PETSC_SUCCESS);
242 : }
|
