Actual source code: bvcontour.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: */
10: /*
11: BV developer functions needed in contour integral methods
12: */
14: #include <slepc/private/bvimpl.h>
15: #include <slepcblaslapack.h>
17: #define p_id(i) (i*subcomm->n + subcomm->color)
19: /*@
20: BVScatter - Scatters the columns of a BV to another BV created in a
21: subcommunicator.
23: Collective
25: Input Parameters:
26: + Vin - input basis vectors (defined on the whole communicator)
27: . scat - VecScatter object that contains the info for the communication
28: - xdup - an auxiliary vector
30: Output Parameter:
31: . Vout - output basis vectors (defined on the subcommunicator)
33: Notes:
34: Currently implemented as a loop for each the active column, where each
35: column is scattered independently. The vector xdup is defined on the
36: contiguous parent communicator and have enough space to store one
37: duplicate of the original vector per each subcommunicator.
39: Level: developer
41: .seealso: BVGetColumn()
42: @*/
43: PetscErrorCode BVScatter(BV Vin,BV Vout,VecScatter scat,Vec xdup)
44: {
45: PetscInt i;
46: Vec v,z;
47: const PetscScalar *array;
49: PetscFunctionBegin;
54: PetscCall(BVCreateVec(Vout,&z));
55: for (i=Vin->l;i<Vin->k;i++) {
56: PetscCall(BVGetColumn(Vin,i,&v));
57: PetscCall(VecScatterBegin(scat,v,xdup,INSERT_VALUES,SCATTER_FORWARD));
58: PetscCall(VecScatterEnd(scat,v,xdup,INSERT_VALUES,SCATTER_FORWARD));
59: PetscCall(BVRestoreColumn(Vin,i,&v));
60: PetscCall(VecGetArrayRead(xdup,&array));
61: PetscCall(VecPlaceArray(z,array));
62: PetscCall(BVInsertVec(Vout,i,z));
63: PetscCall(VecResetArray(z));
64: PetscCall(VecRestoreArrayRead(xdup,&array));
65: }
66: PetscCall(VecDestroy(&z));
67: PetscFunctionReturn(PETSC_SUCCESS);
68: }
70: /*@
71: BVSumQuadrature - Computes the sum of terms required in the quadrature
72: rule to approximate the contour integral.
74: Collective
76: Input Parameters:
77: + Y - input basis vectors
78: . M - number of moments
79: . L - block size
80: . L_max - maximum block size
81: . w - quadrature weights
82: . zn - normalized quadrature points
83: . scat - (optional) VecScatter object to communicate between subcommunicators
84: . subcomm - subcommunicator layout
85: . npoints - number of points to process by the subcommunicator
86: - useconj - whether conjugate points can be used or not
88: Output Parameter:
89: . S - output basis vectors
91: Notes:
92: This is a generalization of BVMult(). The resulting matrix S consists of M
93: panels of L columns, and the following formula is computed for each panel
94: S_k = sum_j w_j*zn_j^k*Y_j, where Y_j is the j-th panel of Y containing
95: the result of solving T(z_j)^{-1}*X for each integration point j. L_max is
96: the width of the panels in Y.
98: When using subcommunicators, Y is stored in the subcommunicators for a subset
99: of integration points. In that case, the computation is done in the subcomm
100: and then scattered to the whole communicator in S using the VecScatter scat.
101: The value npoints is the number of points to be processed in this subcomm
102: and the flag useconj indicates whether symmetric points can be reused.
104: Level: developer
106: .seealso: BVMult(), BVScatter(), BVDotQuadrature(), RGComputeQuadrature(), RGCanUseConjugates()
107: @*/
108: PetscErrorCode BVSumQuadrature(BV S,BV Y,PetscInt M,PetscInt L,PetscInt L_max,PetscScalar *w,PetscScalar *zn,VecScatter scat,PetscSubcomm subcomm,PetscInt npoints,PetscBool useconj)
109: {
110: PetscInt i,j,k,nloc;
111: Vec v,sj;
112: PetscScalar *ppk,*pv,one=1.0;
114: PetscFunctionBegin;
119: PetscCall(BVGetSizes(Y,&nloc,NULL,NULL));
120: PetscCall(PetscMalloc1(npoints,&ppk));
121: for (i=0;i<npoints;i++) ppk[i] = 1.0;
122: PetscCall(BVCreateVec(Y,&v));
123: for (k=0;k<M;k++) {
124: for (j=0;j<L;j++) {
125: PetscCall(VecSet(v,0.0));
126: for (i=0;i<npoints;i++) {
127: PetscCall(BVSetActiveColumns(Y,i*L_max+j,i*L_max+j+1));
128: PetscCall(BVMultVec(Y,ppk[i]*w[p_id(i)],1.0,v,&one));
129: }
130: if (PetscUnlikely(useconj)) {
131: PetscCall(VecGetArray(v,&pv));
132: for (i=0;i<nloc;i++) pv[i] = 2.0*PetscRealPart(pv[i]);
133: PetscCall(VecRestoreArray(v,&pv));
134: }
135: PetscCall(BVGetColumn(S,k*L+j,&sj));
136: if (PetscUnlikely(scat)) {
137: PetscCall(VecScatterBegin(scat,v,sj,ADD_VALUES,SCATTER_REVERSE));
138: PetscCall(VecScatterEnd(scat,v,sj,ADD_VALUES,SCATTER_REVERSE));
139: } else PetscCall(VecCopy(v,sj));
140: PetscCall(BVRestoreColumn(S,k*L+j,&sj));
141: }
142: for (i=0;i<npoints;i++) ppk[i] *= zn[p_id(i)];
143: }
144: PetscCall(PetscFree(ppk));
145: PetscCall(VecDestroy(&v));
146: PetscFunctionReturn(PETSC_SUCCESS);
147: }
149: /*@
150: BVDotQuadrature - Computes the projection terms required in the quadrature
151: rule to approximate the contour integral.
153: Collective
155: Input Parameters:
156: + Y - first basis vectors
157: . V - second basis vectors
158: . M - number of moments
159: . L - block size
160: . L_max - maximum block size
161: . w - quadrature weights
162: . zn - normalized quadrature points
163: . subcomm - subcommunicator layout
164: . npoints - number of points to process by the subcommunicator
165: - useconj - whether conjugate points can be used or not
167: Output Parameter:
168: . Mu - computed result
170: Notes:
171: This is a generalization of BVDot(). The resulting matrix Mu consists of M
172: blocks of size LxL (placed horizontally), each of them computed as
173: Mu_k = sum_j w_j*zn_j^k*V'*Y_j, where Y_j is the j-th panel of Y containing
174: the result of solving T(z_j)^{-1}*X for each integration point j. L_max is
175: the width of the panels in Y.
177: When using subcommunicators, Y is stored in the subcommunicators for a subset
178: of integration points. In that case, the computation is done in the subcomm
179: and then the final result is combined via reduction.
180: The value npoints is the number of points to be processed in this subcomm
181: and the flag useconj indicates whether symmetric points can be reused.
183: Level: developer
185: .seealso: BVDot(), BVScatter(), BVSumQuadrature(), RGComputeQuadrature(), RGCanUseConjugates()
186: @*/
187: PetscErrorCode BVDotQuadrature(BV Y,BV V,PetscScalar *Mu,PetscInt M,PetscInt L,PetscInt L_max,PetscScalar *w,PetscScalar *zn,PetscSubcomm subcomm,PetscInt npoints,PetscBool useconj)
188: {
189: PetscMPIInt sub_size,count;
190: PetscInt i,j,k,s;
191: PetscScalar *temp,*temp2,*ppk,alp;
192: Mat H;
193: const PetscScalar *pH;
194: MPI_Comm child,parent;
196: PetscFunctionBegin;
200: PetscCall(PetscSubcommGetChild(subcomm,&child));
201: PetscCallMPI(MPI_Comm_size(child,&sub_size));
202: PetscCall(PetscMalloc3(npoints*L*(L+1),&temp,2*M*L*L,&temp2,npoints,&ppk));
203: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,L,L_max*npoints,NULL,&H));
204: PetscCall(PetscArrayzero(temp2,2*M*L*L));
205: PetscCall(BVSetActiveColumns(Y,0,L_max*npoints));
206: PetscCall(BVSetActiveColumns(V,0,L));
207: PetscCall(BVDot(Y,V,H));
208: PetscCall(MatDenseGetArrayRead(H,&pH));
209: for (i=0;i<npoints;i++) {
210: for (j=0;j<L;j++) {
211: for (k=0;k<L;k++) {
212: temp[k+j*L+i*L*L] = pH[k+j*L+i*L*L_max];
213: }
214: }
215: }
216: PetscCall(MatDenseRestoreArrayRead(H,&pH));
217: for (i=0;i<npoints;i++) ppk[i] = 1;
218: for (k=0;k<2*M;k++) {
219: for (j=0;j<L;j++) {
220: for (i=0;i<npoints;i++) {
221: alp = ppk[i]*w[p_id(i)];
222: for (s=0;s<L;s++) {
223: if (!useconj) temp2[s+(j+k*L)*L] += alp*temp[s+(j+i*L)*L];
224: else temp2[s+(j+k*L)*L] += 2.0*PetscRealPart(alp*temp[s+(j+i*L)*L]);
225: }
226: }
227: }
228: for (i=0;i<npoints;i++) ppk[i] *= zn[p_id(i)];
229: }
230: for (i=0;i<2*M*L*L;i++) temp2[i] /= sub_size;
231: PetscCall(PetscMPIIntCast(2*M*L*L,&count));
232: PetscCall(PetscSubcommGetParent(subcomm,&parent));
233: PetscCall(MPIU_Allreduce(temp2,Mu,count,MPIU_SCALAR,MPIU_SUM,parent));
234: PetscCall(PetscFree3(temp,temp2,ppk));
235: PetscCall(MatDestroy(&H));
236: PetscFunctionReturn(PETSC_SUCCESS);
237: }
239: /*@
240: BVTraceQuadrature - Computes an estimate of the number of eigenvalues
241: inside a region via quantities computed in the quadrature rule of
242: contour integral methods.
244: Collective
246: Input Parameters:
247: + Y - first basis vectors
248: . V - second basis vectors
249: . L - block size
250: . L_max - maximum block size
251: . w - quadrature weights
252: . scat - (optional) VecScatter object to communicate between subcommunicators
253: . subcomm - subcommunicator layout
254: . npoints - number of points to process by the subcommunicator
255: - useconj - whether conjugate points can be used or not
257: Output Parameter:
258: . est_eig - estimated eigenvalue count
260: Notes:
261: This function returns an estimation of the number of eigenvalues in the
262: region, computed as trace(V'*S_0), where S_0 is the first panel of S
263: computed by BVSumQuadrature().
265: When using subcommunicators, Y is stored in the subcommunicators for a subset
266: of integration points. In that case, the computation is done in the subcomm
267: and then scattered to the whole communicator in S using the VecScatter scat.
268: The value npoints is the number of points to be processed in this subcomm
269: and the flag useconj indicates whether symmetric points can be reused.
271: Level: developer
273: .seealso: BVScatter(), BVDotQuadrature(), BVSumQuadrature(), RGComputeQuadrature(), RGCanUseConjugates()
274: @*/
275: PetscErrorCode BVTraceQuadrature(BV Y,BV V,PetscInt L,PetscInt L_max,PetscScalar *w,VecScatter scat,PetscSubcomm subcomm,PetscInt npoints,PetscBool useconj,PetscReal *est_eig)
276: {
277: PetscInt i,j;
278: Vec y,yall,vj;
279: PetscScalar dot,sum=0.0,one=1.0;
281: PetscFunctionBegin;
286: PetscCall(BVCreateVec(Y,&y));
287: PetscCall(BVCreateVec(V,&yall));
288: for (j=0;j<L;j++) {
289: PetscCall(VecSet(y,0.0));
290: for (i=0;i<npoints;i++) {
291: PetscCall(BVSetActiveColumns(Y,i*L_max+j,i*L_max+j+1));
292: PetscCall(BVMultVec(Y,w[p_id(i)],1.0,y,&one));
293: }
294: PetscCall(BVGetColumn(V,j,&vj));
295: if (scat) {
296: PetscCall(VecScatterBegin(scat,y,yall,ADD_VALUES,SCATTER_REVERSE));
297: PetscCall(VecScatterEnd(scat,y,yall,ADD_VALUES,SCATTER_REVERSE));
298: PetscCall(VecDot(vj,yall,&dot));
299: } else PetscCall(VecDot(vj,y,&dot));
300: PetscCall(BVRestoreColumn(V,j,&vj));
301: if (useconj) sum += 2.0*PetscRealPart(dot);
302: else sum += dot;
303: }
304: *est_eig = PetscAbsScalar(sum)/(PetscReal)L;
305: PetscCall(VecDestroy(&y));
306: PetscCall(VecDestroy(&yall));
307: PetscFunctionReturn(PETSC_SUCCESS);
308: }
310: static PetscErrorCode BVSVDAndRank_Refine(BV S,PetscReal delta,PetscScalar *pA,PetscReal *sigma,PetscInt *rank)
311: {
312: PetscInt i,j,k,ml=S->k;
313: PetscMPIInt len;
314: PetscScalar *work,*B,*tempB,*sarray,*Q1,*Q2,*temp2,alpha=1.0,beta=0.0;
315: PetscBLASInt l,m,n,lda,ldu,ldvt,lwork,info,ldb,ldc,lds;
316: #if defined(PETSC_USE_COMPLEX)
317: PetscReal *rwork;
318: #endif
320: PetscFunctionBegin;
321: PetscCall(PetscBLASIntCast(S->ld,&lds));
322: PetscCall(BVGetArray(S,&sarray));
323: PetscCall(PetscMalloc6(ml*ml,&temp2,S->n*ml,&Q1,S->n*ml,&Q2,ml*ml,&B,ml*ml,&tempB,5*ml,&work));
324: #if defined(PETSC_USE_COMPLEX)
325: PetscCall(PetscMalloc1(5*ml,&rwork));
326: #endif
327: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
329: PetscCall(PetscArrayzero(B,ml*ml));
330: for (i=0;i<ml;i++) B[i*ml+i]=1;
332: for (k=0;k<2;k++) {
333: PetscCall(PetscBLASIntCast(S->n,&m));
334: PetscCall(PetscBLASIntCast(ml,&l));
335: n = l; lda = m; ldb = m; ldc = l;
336: if (!k) PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&l,&n,&m,&alpha,sarray,&lds,sarray,&lds,&beta,pA,&ldc));
337: else PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&l,&n,&m,&alpha,Q1,&lda,Q1,&ldb,&beta,pA,&ldc));
338: PetscCall(PetscArrayzero(temp2,ml*ml));
339: PetscCall(PetscMPIIntCast(ml*ml,&len));
340: PetscCall(MPIU_Allreduce(pA,temp2,len,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)S)));
342: PetscCall(PetscBLASIntCast(ml,&m));
343: n = m; lda = m; lwork = 5*m, ldu = 1; ldvt = 1;
344: #if defined(PETSC_USE_COMPLEX)
345: PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("O","N",&m,&n,temp2,&lda,sigma,NULL,&ldu,NULL,&ldvt,work,&lwork,rwork,&info));
346: #else
347: PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("O","N",&m,&n,temp2,&lda,sigma,NULL,&ldu,NULL,&ldvt,work,&lwork,&info));
348: #endif
349: SlepcCheckLapackInfo("gesvd",info);
351: PetscCall(PetscBLASIntCast(S->n,&l));
352: PetscCall(PetscBLASIntCast(ml,&n));
353: m = n; lda = l; ldb = m; ldc = l;
354: if (!k) PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&l,&n,&m,&alpha,sarray,&lds,temp2,&ldb,&beta,Q1,&ldc));
355: else PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&l,&n,&m,&alpha,Q1,&lda,temp2,&ldb,&beta,Q2,&ldc));
357: PetscCall(PetscBLASIntCast(ml,&l));
358: m = l; n = l; lda = l; ldb = m; ldc = l;
359: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&l,&n,&m,&alpha,B,&lda,temp2,&ldb,&beta,tempB,&ldc));
360: for (i=0;i<ml;i++) {
361: sigma[i] = PetscSqrtReal(sigma[i]);
362: for (j=0;j<S->n;j++) {
363: if (k%2) Q2[j+i*S->n] /= sigma[i];
364: else Q1[j+i*S->n] /= sigma[i];
365: }
366: for (j=0;j<ml;j++) B[j+i*ml] = tempB[j+i*ml]*sigma[i];
367: }
368: }
370: PetscCall(PetscBLASIntCast(ml,&m));
371: n = m; lda = m; ldu=1; ldvt=1;
372: #if defined(PETSC_USE_COMPLEX)
373: PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","O",&m,&n,B,&lda,sigma,NULL,&ldu,NULL,&ldvt,work,&lwork,rwork,&info));
374: #else
375: PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","O",&m,&n,B,&lda,sigma,NULL,&ldu,NULL,&ldvt,work,&lwork,&info));
376: #endif
377: SlepcCheckLapackInfo("gesvd",info);
379: PetscCall(PetscBLASIntCast(S->n,&l));
380: PetscCall(PetscBLASIntCast(ml,&n));
381: m = n; lda = l; ldb = m;
382: if (k%2) PetscCallBLAS("BLASgemm",BLASgemm_("N","T",&l,&n,&m,&alpha,Q1,&lda,B,&ldb,&beta,sarray,&lds));
383: else PetscCallBLAS("BLASgemm",BLASgemm_("N","T",&l,&n,&m,&alpha,Q2,&lda,B,&ldb,&beta,sarray,&lds));
385: PetscCall(PetscFPTrapPop());
386: PetscCall(BVRestoreArray(S,&sarray));
388: if (rank) {
389: (*rank) = 0;
390: for (i=0;i<ml;i++) {
391: if (sigma[i]/PetscMax(sigma[0],1.0)>delta) (*rank)++;
392: }
393: }
394: PetscCall(PetscFree6(temp2,Q1,Q2,B,tempB,work));
395: #if defined(PETSC_USE_COMPLEX)
396: PetscCall(PetscFree(rwork));
397: #endif
398: PetscFunctionReturn(PETSC_SUCCESS);
399: }
401: static PetscErrorCode BVSVDAndRank_QR(BV S,PetscReal delta,PetscScalar *pA,PetscReal *sigma,PetscInt *rank)
402: {
403: PetscInt i,n,ml=S->k;
404: PetscBLASInt m,lda,lwork,info;
405: PetscScalar *work;
406: PetscReal *rwork;
407: Mat A;
408: Vec v;
410: PetscFunctionBegin;
411: /* Compute QR factorizaton of S */
412: PetscCall(BVGetSizes(S,NULL,&n,NULL));
413: n = PetscMin(n,ml);
414: PetscCall(BVSetActiveColumns(S,0,n));
415: PetscCall(PetscArrayzero(pA,ml*n));
416: PetscCall(MatCreateDense(PETSC_COMM_SELF,n,n,PETSC_DECIDE,PETSC_DECIDE,pA,&A));
417: PetscCall(BVOrthogonalize(S,A));
418: if (n<ml) {
419: /* the rest of the factorization */
420: for (i=n;i<ml;i++) {
421: PetscCall(BVGetColumn(S,i,&v));
422: PetscCall(BVOrthogonalizeVec(S,v,pA+i*n,NULL,NULL));
423: PetscCall(BVRestoreColumn(S,i,&v));
424: }
425: }
426: PetscCall(PetscBLASIntCast(n,&lda));
427: PetscCall(PetscBLASIntCast(ml,&m));
428: PetscCall(PetscMalloc2(5*ml,&work,5*ml,&rwork));
429: lwork = 5*m;
430: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
431: #if !defined (PETSC_USE_COMPLEX)
432: PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("O","N",&lda,&m,pA,&lda,sigma,NULL,&lda,NULL,&lda,work,&lwork,&info));
433: #else
434: PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("O","N",&lda,&m,pA,&lda,sigma,NULL,&lda,NULL,&lda,work,&lwork,rwork,&info));
435: #endif
436: SlepcCheckLapackInfo("gesvd",info);
437: PetscCall(PetscFPTrapPop());
438: *rank = 0;
439: for (i=0;i<n;i++) {
440: if (sigma[i]/PetscMax(sigma[0],1)>delta) (*rank)++;
441: }
442: /* n first columns of A have the left singular vectors */
443: PetscCall(BVMultInPlace(S,A,0,*rank));
444: PetscCall(PetscFree2(work,rwork));
445: PetscCall(MatDestroy(&A));
446: PetscFunctionReturn(PETSC_SUCCESS);
447: }
449: static PetscErrorCode BVSVDAndRank_QR_CAA(BV S,PetscInt M,PetscInt L,PetscReal delta,PetscScalar *pA,PetscReal *sigma,PetscInt *rank)
450: {
451: PetscInt i,j,n,ml=S->k;
452: PetscBLASInt m,k_,lda,lwork,info;
453: PetscScalar *work,*T,*U,*R,sone=1.0,zero=0.0;
454: PetscReal *rwork;
455: Mat A;
457: PetscFunctionBegin;
458: /* Compute QR factorizaton of S */
459: PetscCall(BVGetSizes(S,NULL,&n,NULL));
460: PetscCheck(n>=ml,PetscObjectComm((PetscObject)S),PETSC_ERR_SUP,"The QR_CAA method does not support problem size n < m*L");
461: PetscCall(BVSetActiveColumns(S,0,ml));
462: PetscCall(PetscArrayzero(pA,ml*ml));
463: PetscCall(MatCreateDense(PETSC_COMM_SELF,ml,ml,PETSC_DECIDE,PETSC_DECIDE,pA,&A));
464: PetscCall(BVOrthogonalize(S,A));
465: PetscCall(MatDestroy(&A));
467: /* SVD of first (M-1)*L diagonal block */
468: PetscCall(PetscBLASIntCast((M-1)*L,&m));
469: PetscCall(PetscMalloc5(m*m,&T,m*m,&R,m*m,&U,5*ml,&work,5*ml,&rwork));
470: for (j=0;j<m;j++) PetscCall(PetscArraycpy(R+j*m,pA+j*ml,m));
471: lwork = 5*m;
472: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
473: #if !defined (PETSC_USE_COMPLEX)
474: PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("S","O",&m,&m,R,&m,sigma,U,&m,NULL,&m,work,&lwork,&info));
475: #else
476: PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("S","O",&m,&m,R,&m,sigma,U,&m,NULL,&m,work,&lwork,rwork,&info));
477: #endif
478: SlepcCheckLapackInfo("gesvd",info);
479: PetscCall(PetscFPTrapPop());
480: *rank = 0;
481: for (i=0;i<m;i++) {
482: if (sigma[i]/PetscMax(sigma[0],1)>delta) (*rank)++;
483: }
484: PetscCall(MatCreateDense(PETSC_COMM_SELF,m,m,PETSC_DECIDE,PETSC_DECIDE,U,&A));
485: PetscCall(BVSetActiveColumns(S,0,m));
486: PetscCall(BVMultInPlace(S,A,0,*rank));
487: PetscCall(MatDestroy(&A));
488: /* Projected linear system */
489: /* m first columns of A have the right singular vectors */
490: PetscCall(PetscBLASIntCast(*rank,&k_));
491: PetscCall(PetscBLASIntCast(ml,&lda));
492: PetscCallBLAS("BLASgemm",BLASgemm_("N","C",&m,&k_,&m,&sone,pA+L*lda,&lda,R,&m,&zero,T,&m));
493: PetscCall(PetscArrayzero(pA,ml*ml));
494: PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&k_,&k_,&m,&sone,U,&m,T,&m,&zero,pA,&k_));
495: for (j=0;j<k_;j++) for (i=0;i<k_;i++) pA[j*k_+i] /= sigma[j];
496: PetscCall(PetscFree5(T,R,U,work,rwork));
497: PetscFunctionReturn(PETSC_SUCCESS);
498: }
500: /*@
501: BVSVDAndRank - Compute the SVD (left singular vectors only, and singular
502: values) and determine the numerical rank according to a tolerance.
504: Collective
506: Input Parameters:
507: + S - the basis vectors
508: . m - the moment degree
509: . l - the block size
510: . delta - the tolerance used to determine the rank
511: - meth - the method to be used
513: Output Parameters:
514: + A - workspace, on output contains relevant values in the CAA method
515: . sigma - computed singular values
516: - rank - estimated rank (optional)
518: Notes:
519: This function computes [U,Sigma,V] = svd(S) and replaces S with U.
520: The current implementation computes this via S'*S, and it may include
521: some kind of iterative refinement to improve accuracy in some cases.
523: The parameters m and l refer to the moment and block size of contour
524: integral methods. All columns up to m*l are modified, and the active
525: columns are set to 0..m*l.
527: The method is one of BV_SVD_METHOD_REFINE, BV_SVD_METHOD_QR, BV_SVD_METHOD_QR_CAA.
529: The A workspace should be m*l*m*l in size.
531: Once the decomposition is computed, the numerical rank is estimated
532: by counting the number of singular values that are larger than the
533: tolerance delta, relative to the first singular value.
535: Level: developer
537: .seealso: BVSetActiveColumns()
538: @*/
539: PetscErrorCode BVSVDAndRank(BV S,PetscInt m,PetscInt l,PetscReal delta,BVSVDMethod meth,PetscScalar *A,PetscReal *sigma,PetscInt *rank)
540: {
541: PetscFunctionBegin;
547: PetscAssertPointer(A,6);
548: PetscAssertPointer(sigma,7);
549: PetscAssertPointer(rank,8);
551: PetscCall(PetscLogEventBegin(BV_SVDAndRank,S,0,0,0));
552: PetscCall(BVSetActiveColumns(S,0,m*l));
553: switch (meth) {
554: case BV_SVD_METHOD_REFINE:
555: PetscCall(BVSVDAndRank_Refine(S,delta,A,sigma,rank));
556: break;
557: case BV_SVD_METHOD_QR:
558: PetscCall(BVSVDAndRank_QR(S,delta,A,sigma,rank));
559: break;
560: case BV_SVD_METHOD_QR_CAA:
561: PetscCall(BVSVDAndRank_QR_CAA(S,m,l,delta,A,sigma,rank));
562: break;
563: }
564: PetscCall(PetscLogEventEnd(BV_SVDAndRank,S,0,0,0));
565: PetscFunctionReturn(PETSC_SUCCESS);
566: }
568: /*@
569: BVCISSResizeBases - Resize the bases involved in CISS solvers when the L grows.
571: Logically Collective
573: Input Parameters:
574: + S - basis of L*M columns
575: . V - basis of L columns (may be associated to subcommunicators)
576: . Y - basis of npoints*L columns
577: . Lold - old value of L
578: . Lnew - new value of L
579: . M - the moment size
580: - npoints - number of integration points
582: Level: developer
584: .seealso: BVResize()
585: @*/
586: PetscErrorCode BVCISSResizeBases(BV S,BV V,BV Y,PetscInt Lold,PetscInt Lnew,PetscInt M,PetscInt npoints)
587: {
588: PetscInt i,j;
590: PetscFunctionBegin;
599: PetscCall(BVResize(S,Lnew*M,PETSC_FALSE));
600: PetscCall(BVResize(V,Lnew,PETSC_TRUE));
601: PetscCall(BVResize(Y,Lnew*npoints,PETSC_TRUE));
602: /* columns of Y are interleaved */
603: for (i=npoints-1;i>=0;i--) {
604: for (j=Lold-1;j>=0;j--) PetscCall(BVCopyColumn(Y,i*Lold+j,i*Lnew+j));
605: }
606: PetscFunctionReturn(PETSC_SUCCESS);
607: }