Actual source code: bvcontour.c
slepc-3.18.0 2022-10-01
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 on Vin
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;
47: const PetscScalar *array;
53: for (i=Vin->l;i<Vin->k;i++) {
54: BVGetColumn(Vin,i,&v);
55: VecScatterBegin(scat,v,xdup,INSERT_VALUES,SCATTER_FORWARD);
56: VecScatterEnd(scat,v,xdup,INSERT_VALUES,SCATTER_FORWARD);
57: BVRestoreColumn(Vin,i,&v);
58: VecGetArrayRead(xdup,&array);
59: VecPlaceArray(Vout->t,array);
60: BVInsertVec(Vout,i,Vout->t);
61: VecResetArray(Vout->t);
62: VecRestoreArrayRead(xdup,&array);
63: }
64: return 0;
65: }
67: /*@
68: BVSumQuadrature - Computes the sum of terms required in the quadrature
69: rule to approximate the contour integral.
71: Collective on S
73: Input Parameters:
74: + Y - input basis vectors
75: . M - number of moments
76: . L - block size
77: . L_max - maximum block size
78: . w - quadrature weights
79: . zn - normalized quadrature points
80: . scat - (optional) VecScatter object to communicate between subcommunicators
81: . subcomm - subcommunicator layout
82: . npoints - number of points to process by the subcommunicator
83: - useconj - whether conjugate points can be used or not
85: Output Parameter:
86: . S - output basis vectors
88: Notes:
89: This is a generalization of BVMult(). The resulting matrix S consists of M
90: panels of L columns, and the following formula is computed for each panel
91: S_k = sum_j w_j*zn_j^k*Y_j, where Y_j is the j-th panel of Y containing
92: the result of solving T(z_j)^{-1}*X for each integration point j. L_max is
93: the width of the panels in Y.
95: When using subcommunicators, Y is stored in the subcommunicators for a subset
96: of integration points. In that case, the computation is done in the subcomm
97: and then scattered to the whole communicator in S using the VecScatter scat.
98: The value npoints is the number of points to be processed in this subcomm
99: and the flag useconj indicates whether symmetric points can be reused.
101: Level: developer
103: .seealso: BVMult(), BVScatter(), BVDotQuadrature(), RGComputeQuadrature(), RGCanUseConjugates()
104: @*/
105: 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)
106: {
107: PetscInt i,j,k,nloc;
108: Vec v,sj;
109: PetscScalar *ppk,*pv,one=1.0;
115: BVGetSizes(Y,&nloc,NULL,NULL);
116: PetscMalloc1(npoints,&ppk);
117: for (i=0;i<npoints;i++) ppk[i] = 1.0;
118: BVCreateVec(Y,&v);
119: for (k=0;k<M;k++) {
120: for (j=0;j<L;j++) {
121: VecSet(v,0.0);
122: for (i=0;i<npoints;i++) {
123: BVSetActiveColumns(Y,i*L_max+j,i*L_max+j+1);
124: BVMultVec(Y,ppk[i]*w[p_id(i)],1.0,v,&one);
125: }
126: if (PetscUnlikely(useconj)) {
127: VecGetArray(v,&pv);
128: for (i=0;i<nloc;i++) pv[i] = 2.0*PetscRealPart(pv[i]);
129: VecRestoreArray(v,&pv);
130: }
131: BVGetColumn(S,k*L+j,&sj);
132: if (PetscUnlikely(scat)) {
133: VecScatterBegin(scat,v,sj,ADD_VALUES,SCATTER_REVERSE);
134: VecScatterEnd(scat,v,sj,ADD_VALUES,SCATTER_REVERSE);
135: } else VecCopy(v,sj);
136: BVRestoreColumn(S,k*L+j,&sj);
137: }
138: for (i=0;i<npoints;i++) ppk[i] *= zn[p_id(i)];
139: }
140: PetscFree(ppk);
141: VecDestroy(&v);
142: return 0;
143: }
145: /*@
146: BVDotQuadrature - Computes the projection terms required in the quadrature
147: rule to approximate the contour integral.
149: Collective on Y
151: Input Parameters:
152: + Y - first basis vectors
153: . V - second basis vectors
154: . M - number of moments
155: . L - block size
156: . L_max - maximum block size
157: . w - quadrature weights
158: . zn - normalized quadrature points
159: . subcomm - subcommunicator layout
160: . npoints - number of points to process by the subcommunicator
161: - useconj - whether conjugate points can be used or not
163: Output Parameter:
164: . Mu - computed result
166: Notes:
167: This is a generalization of BVDot(). The resulting matrix Mu consists of M
168: blocks of size LxL (placed horizontally), each of them computed as
169: Mu_k = sum_j w_j*zn_j^k*V'*Y_j, where Y_j is the j-th panel of Y containing
170: the result of solving T(z_j)^{-1}*X for each integration point j. L_max is
171: the width of the panels in Y.
173: When using subcommunicators, Y is stored in the subcommunicators for a subset
174: of integration points. In that case, the computation is done in the subcomm
175: and then the final result is combined via reduction.
176: The value npoints is the number of points to be processed in this subcomm
177: and the flag useconj indicates whether symmetric points can be reused.
179: Level: developer
181: .seealso: BVDot(), BVScatter(), BVSumQuadrature(), RGComputeQuadrature(), RGCanUseConjugates()
182: @*/
183: 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)
184: {
185: PetscMPIInt sub_size,count;
186: PetscInt i,j,k,s;
187: PetscScalar *temp,*temp2,*ppk,alp;
188: Mat H;
189: const PetscScalar *pH;
190: MPI_Comm child,parent;
195: PetscSubcommGetChild(subcomm,&child);
196: MPI_Comm_size(child,&sub_size);
197: PetscMalloc3(npoints*L*(L+1),&temp,2*M*L*L,&temp2,npoints,&ppk);
198: MatCreateSeqDense(PETSC_COMM_SELF,L,L_max*npoints,NULL,&H);
199: PetscArrayzero(temp2,2*M*L*L);
200: BVSetActiveColumns(Y,0,L_max*npoints);
201: BVSetActiveColumns(V,0,L);
202: BVDot(Y,V,H);
203: MatDenseGetArrayRead(H,&pH);
204: for (i=0;i<npoints;i++) {
205: for (j=0;j<L;j++) {
206: for (k=0;k<L;k++) {
207: temp[k+j*L+i*L*L] = pH[k+j*L+i*L*L_max];
208: }
209: }
210: }
211: MatDenseRestoreArrayRead(H,&pH);
212: for (i=0;i<npoints;i++) ppk[i] = 1;
213: for (k=0;k<2*M;k++) {
214: for (j=0;j<L;j++) {
215: for (i=0;i<npoints;i++) {
216: alp = ppk[i]*w[p_id(i)];
217: for (s=0;s<L;s++) {
218: if (!useconj) temp2[s+(j+k*L)*L] += alp*temp[s+(j+i*L)*L];
219: else temp2[s+(j+k*L)*L] += 2.0*PetscRealPart(alp*temp[s+(j+i*L)*L]);
220: }
221: }
222: }
223: for (i=0;i<npoints;i++) ppk[i] *= zn[p_id(i)];
224: }
225: for (i=0;i<2*M*L*L;i++) temp2[i] /= sub_size;
226: PetscMPIIntCast(2*M*L*L,&count);
227: PetscSubcommGetParent(subcomm,&parent);
228: MPIU_Allreduce(temp2,Mu,count,MPIU_SCALAR,MPIU_SUM,parent);
229: PetscFree3(temp,temp2,ppk);
230: MatDestroy(&H);
231: return 0;
232: }
234: /*@
235: BVTraceQuadrature - Computes an estimate of the number of eigenvalues
236: inside a region via quantities computed in the quadrature rule of
237: contour integral methods.
239: Collective on Y
241: Input Parameters:
242: + Y - first basis vectors
243: . V - second basis vectors
244: . L - block size
245: . L_max - maximum block size
246: . w - quadrature weights
247: . scat - (optional) VecScatter object to communicate between subcommunicators
248: . subcomm - subcommunicator layout
249: . npoints - number of points to process by the subcommunicator
250: - useconj - whether conjugate points can be used or not
252: Output Parameter:
253: . est_eig - estimated eigenvalue count
255: Notes:
256: This function returns an estimation of the number of eigenvalues in the
257: region, computed as trace(V'*S_0), where S_0 is the first panel of S
258: computed by BVSumQuadrature().
260: When using subcommunicators, Y is stored in the subcommunicators for a subset
261: of integration points. In that case, the computation is done in the subcomm
262: and then scattered to the whole communicator in S using the VecScatter scat.
263: The value npoints is the number of points to be processed in this subcomm
264: and the flag useconj indicates whether symmetric points can be reused.
266: Level: developer
268: .seealso: BVScatter(), BVDotQuadrature(), BVSumQuadrature(), RGComputeQuadrature(), RGCanUseConjugates()
269: @*/
270: PetscErrorCode BVTraceQuadrature(BV Y,BV V,PetscInt L,PetscInt L_max,PetscScalar *w,VecScatter scat,PetscSubcomm subcomm,PetscInt npoints,PetscBool useconj,PetscReal *est_eig)
271: {
272: PetscInt i,j;
273: Vec y,yall,vj;
274: PetscScalar dot,sum=0.0,one=1.0;
280: BVCreateVec(Y,&y);
281: BVCreateVec(V,&yall);
282: for (j=0;j<L;j++) {
283: VecSet(y,0.0);
284: for (i=0;i<npoints;i++) {
285: BVSetActiveColumns(Y,i*L_max+j,i*L_max+j+1);
286: BVMultVec(Y,w[p_id(i)],1.0,y,&one);
287: }
288: BVGetColumn(V,j,&vj);
289: if (scat) {
290: VecScatterBegin(scat,y,yall,ADD_VALUES,SCATTER_REVERSE);
291: VecScatterEnd(scat,y,yall,ADD_VALUES,SCATTER_REVERSE);
292: VecDot(vj,yall,&dot);
293: } else VecDot(vj,y,&dot);
294: BVRestoreColumn(V,j,&vj);
295: if (useconj) sum += 2.0*PetscRealPart(dot);
296: else sum += dot;
297: }
298: *est_eig = PetscAbsScalar(sum)/(PetscReal)L;
299: VecDestroy(&y);
300: VecDestroy(&yall);
301: return 0;
302: }
304: PetscErrorCode BVSVDAndRank_Refine(BV S,PetscReal delta,PetscScalar *pA,PetscReal *sigma,PetscInt *rank)
305: {
306: PetscInt i,j,k,ml=S->k;
307: PetscMPIInt len;
308: PetscScalar *work,*B,*tempB,*sarray,*Q1,*Q2,*temp2,alpha=1.0,beta=0.0;
309: PetscBLASInt l,m,n,lda,ldu,ldvt,lwork,info,ldb,ldc;
310: #if defined(PETSC_USE_COMPLEX)
311: PetscReal *rwork;
312: #endif
314: BVGetArray(S,&sarray);
315: PetscMalloc6(ml*ml,&temp2,S->n*ml,&Q1,S->n*ml,&Q2,ml*ml,&B,ml*ml,&tempB,5*ml,&work);
316: #if defined(PETSC_USE_COMPLEX)
317: PetscMalloc1(5*ml,&rwork);
318: #endif
319: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
321: PetscArrayzero(B,ml*ml);
322: for (i=0;i<ml;i++) B[i*ml+i]=1;
324: for (k=0;k<2;k++) {
325: PetscBLASIntCast(S->n,&m);
326: PetscBLASIntCast(ml,&l);
327: n = l; lda = m; ldb = m; ldc = l;
328: if (!k) PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&l,&n,&m,&alpha,sarray,&lda,sarray,&ldb,&beta,pA,&ldc));
329: else PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&l,&n,&m,&alpha,Q1,&lda,Q1,&ldb,&beta,pA,&ldc));
330: PetscArrayzero(temp2,ml*ml);
331: PetscMPIIntCast(ml*ml,&len);
332: MPIU_Allreduce(pA,temp2,len,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)S));
334: PetscBLASIntCast(ml,&m);
335: n = m; lda = m; lwork = 5*m, ldu = 1; ldvt = 1;
336: #if defined(PETSC_USE_COMPLEX)
337: PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("O","N",&m,&n,temp2,&lda,sigma,NULL,&ldu,NULL,&ldvt,work,&lwork,rwork,&info));
338: #else
339: PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("O","N",&m,&n,temp2,&lda,sigma,NULL,&ldu,NULL,&ldvt,work,&lwork,&info));
340: #endif
341: SlepcCheckLapackInfo("gesvd",info);
343: PetscBLASIntCast(S->n,&l);
344: PetscBLASIntCast(ml,&n);
345: m = n; lda = l; ldb = m; ldc = l;
346: if (!k) PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&l,&n,&m,&alpha,sarray,&lda,temp2,&ldb,&beta,Q1,&ldc));
347: else PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&l,&n,&m,&alpha,Q1,&lda,temp2,&ldb,&beta,Q2,&ldc));
349: PetscBLASIntCast(ml,&l);
350: m = l; n = l; lda = l; ldb = m; ldc = l;
351: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&l,&n,&m,&alpha,B,&lda,temp2,&ldb,&beta,tempB,&ldc));
352: for (i=0;i<ml;i++) {
353: sigma[i] = PetscSqrtReal(sigma[i]);
354: for (j=0;j<S->n;j++) {
355: if (k%2) Q2[j+i*S->n] /= sigma[i];
356: else Q1[j+i*S->n] /= sigma[i];
357: }
358: for (j=0;j<ml;j++) B[j+i*ml] = tempB[j+i*ml]*sigma[i];
359: }
360: }
362: PetscBLASIntCast(ml,&m);
363: n = m; lda = m; ldu=1; ldvt=1;
364: #if defined(PETSC_USE_COMPLEX)
365: PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","O",&m,&n,B,&lda,sigma,NULL,&ldu,NULL,&ldvt,work,&lwork,rwork,&info));
366: #else
367: PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","O",&m,&n,B,&lda,sigma,NULL,&ldu,NULL,&ldvt,work,&lwork,&info));
368: #endif
369: SlepcCheckLapackInfo("gesvd",info);
371: PetscBLASIntCast(S->n,&l);
372: PetscBLASIntCast(ml,&n);
373: m = n; lda = l; ldb = m; ldc = l;
374: if (k%2) PetscCallBLAS("BLASgemm",BLASgemm_("N","T",&l,&n,&m,&alpha,Q1,&lda,B,&ldb,&beta,sarray,&ldc));
375: else PetscCallBLAS("BLASgemm",BLASgemm_("N","T",&l,&n,&m,&alpha,Q2,&lda,B,&ldb,&beta,sarray,&ldc));
377: PetscFPTrapPop();
378: BVRestoreArray(S,&sarray);
380: if (rank) {
381: (*rank) = 0;
382: for (i=0;i<ml;i++) {
383: if (sigma[i]/PetscMax(sigma[0],1.0)>delta) (*rank)++;
384: }
385: }
386: PetscFree6(temp2,Q1,Q2,B,tempB,work);
387: #if defined(PETSC_USE_COMPLEX)
388: PetscFree(rwork);
389: #endif
390: return 0;
391: }
393: PetscErrorCode BVSVDAndRank_QR(BV S,PetscReal delta,PetscScalar *pA,PetscReal *sigma,PetscInt *rank)
394: {
395: PetscInt i,n,ml=S->k;
396: PetscBLASInt m,lda,lwork,info;
397: PetscScalar *work;
398: PetscReal *rwork;
399: Mat A;
400: Vec v;
402: /* Compute QR factorizaton of S */
403: BVGetSizes(S,NULL,&n,NULL);
404: n = PetscMin(n,ml);
405: BVSetActiveColumns(S,0,n);
406: PetscArrayzero(pA,ml*n);
407: MatCreateDense(PETSC_COMM_SELF,n,n,PETSC_DECIDE,PETSC_DECIDE,pA,&A);
408: BVOrthogonalize(S,A);
409: if (n<ml) {
410: /* the rest of the factorization */
411: for (i=n;i<ml;i++) {
412: BVGetColumn(S,i,&v);
413: BVOrthogonalizeVec(S,v,pA+i*n,NULL,NULL);
414: BVRestoreColumn(S,i,&v);
415: }
416: }
417: PetscBLASIntCast(n,&lda);
418: PetscBLASIntCast(ml,&m);
419: PetscMalloc2(5*ml,&work,5*ml,&rwork);
420: lwork = 5*m;
421: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
422: #if !defined (PETSC_USE_COMPLEX)
423: PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("O","N",&lda,&m,pA,&lda,sigma,NULL,&lda,NULL,&lda,work,&lwork,&info));
424: #else
425: PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("O","N",&lda,&m,pA,&lda,sigma,NULL,&lda,NULL,&lda,work,&lwork,rwork,&info));
426: #endif
427: SlepcCheckLapackInfo("gesvd",info);
428: PetscFPTrapPop();
429: *rank = 0;
430: for (i=0;i<n;i++) {
431: if (sigma[i]/PetscMax(sigma[0],1)>delta) (*rank)++;
432: }
433: /* n first columns of A have the left singular vectors */
434: BVMultInPlace(S,A,0,*rank);
435: PetscFree2(work,rwork);
436: MatDestroy(&A);
437: return 0;
438: }
440: PetscErrorCode BVSVDAndRank_QR_CAA(BV S,PetscInt M,PetscInt L,PetscReal delta,PetscScalar *pA,PetscReal *sigma,PetscInt *rank)
441: {
442: PetscInt i,j,n,ml=S->k;
443: PetscBLASInt m,k_,lda,lwork,info;
444: PetscScalar *work,*T,*U,*R,sone=1.0,zero=0.0;
445: PetscReal *rwork;
446: Mat A;
448: /* Compute QR factorizaton of S */
449: BVGetSizes(S,NULL,&n,NULL);
451: BVSetActiveColumns(S,0,ml);
452: PetscArrayzero(pA,ml*ml);
453: MatCreateDense(PETSC_COMM_SELF,ml,ml,PETSC_DECIDE,PETSC_DECIDE,pA,&A);
454: BVOrthogonalize(S,A);
455: MatDestroy(&A);
457: /* SVD of first (M-1)*L diagonal block */
458: PetscBLASIntCast((M-1)*L,&m);
459: PetscMalloc5(m*m,&T,m*m,&R,m*m,&U,5*ml,&work,5*ml,&rwork);
460: for (j=0;j<m;j++) PetscArraycpy(R+j*m,pA+j*ml,m);
461: lwork = 5*m;
462: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
463: #if !defined (PETSC_USE_COMPLEX)
464: PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("S","O",&m,&m,R,&m,sigma,U,&m,NULL,&m,work,&lwork,&info));
465: #else
466: PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("S","O",&m,&m,R,&m,sigma,U,&m,NULL,&m,work,&lwork,rwork,&info));
467: #endif
468: SlepcCheckLapackInfo("gesvd",info);
469: PetscFPTrapPop();
470: *rank = 0;
471: for (i=0;i<m;i++) {
472: if (sigma[i]/PetscMax(sigma[0],1)>delta) (*rank)++;
473: }
474: MatCreateDense(PETSC_COMM_SELF,m,m,PETSC_DECIDE,PETSC_DECIDE,U,&A);
475: BVSetActiveColumns(S,0,m);
476: BVMultInPlace(S,A,0,*rank);
477: MatDestroy(&A);
478: /* Projected linear system */
479: /* m first columns of A have the right singular vectors */
480: PetscBLASIntCast(*rank,&k_);
481: PetscBLASIntCast(ml,&lda);
482: PetscCallBLAS("BLASgemm",BLASgemm_("N","C",&m,&k_,&m,&sone,pA+L*lda,&lda,R,&m,&zero,T,&m));
483: PetscArrayzero(pA,ml*ml);
484: PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&k_,&k_,&m,&sone,U,&m,T,&m,&zero,pA,&k_));
485: for (j=0;j<k_;j++) for (i=0;i<k_;i++) pA[j*k_+i] /= sigma[j];
486: PetscFree5(T,R,U,work,rwork);
487: return 0;
488: }
490: /*@
491: BVSVDAndRank - Compute the SVD (left singular vectors only, and singular
492: values) and determine the numerical rank according to a tolerance.
494: Collective on S
496: Input Parameters:
497: + S - the basis vectors
498: . m - the moment degree
499: . l - the block size
500: . delta - the tolerance used to determine the rank
501: - meth - the method to be used
503: Output Parameters:
504: + A - workspace, on output contains relevant values in the CAA method
505: . sigma - computed singular values
506: - rank - estimated rank (optional)
508: Notes:
509: This function computes [U,Sigma,V] = svd(S) and replaces S with U.
510: The current implementation computes this via S'*S, and it may include
511: some kind of iterative refinement to improve accuracy in some cases.
513: The parameters m and l refer to the moment and block size of contour
514: integral methods. All columns up to m*l are modified, and the active
515: columns are set to 0..m*l.
517: The method is one of BV_SVD_METHOD_REFINE, BV_SVD_METHOD_QR, BV_SVD_METHOD_QR_CAA.
519: The A workspace should be m*l*m*l in size.
521: Once the decomposition is computed, the numerical rank is estimated
522: by counting the number of singular values that are larger than the
523: tolerance delta, relative to the first singular value.
525: Level: developer
527: .seealso: BVSetActiveColumns()
528: @*/
529: PetscErrorCode BVSVDAndRank(BV S,PetscInt m,PetscInt l,PetscReal delta,BVSVDMethod meth,PetscScalar *A,PetscReal *sigma,PetscInt *rank)
530: {
540: PetscLogEventBegin(BV_SVDAndRank,S,0,0,0);
541: BVSetActiveColumns(S,0,m*l);
542: switch (meth) {
543: case BV_SVD_METHOD_REFINE:
544: BVSVDAndRank_Refine(S,delta,A,sigma,rank);
545: break;
546: case BV_SVD_METHOD_QR:
547: BVSVDAndRank_QR(S,delta,A,sigma,rank);
548: break;
549: case BV_SVD_METHOD_QR_CAA:
550: BVSVDAndRank_QR_CAA(S,m,l,delta,A,sigma,rank);
551: break;
552: }
553: PetscLogEventEnd(BV_SVDAndRank,S,0,0,0);
554: return 0;
555: }
557: /*@
558: BVCISSResizeBases - Resize the bases involved in CISS solvers when the L grows.
560: Collective on S
562: Input Parameters:
563: + S - basis of L*M columns
564: . V - basis of L columns (may be associated to subcommunicators)
565: . Y - basis of npoints*L columns
566: . Lold - old value of L
567: . Lnew - new value of L
568: . M - the moment size
569: - npoints - number of integration points
571: Level: developer
573: .seealso: BVResize()
574: @*/
575: PetscErrorCode BVCISSResizeBases(BV S,BV V,BV Y,PetscInt Lold,PetscInt Lnew,PetscInt M,PetscInt npoints)
576: {
577: PetscInt i,j;
587: BVResize(S,Lnew*M,PETSC_FALSE);
588: BVResize(V,Lnew,PETSC_TRUE);
589: BVResize(Y,Lnew*npoints,PETSC_TRUE);
590: /* columns of Y are interleaved */
591: for (i=npoints-1;i>=0;i--) {
592: for (j=Lold-1;j>=0;j--) BVCopyColumn(Y,i*Lold+j,i*Lnew+j);
593: }
594: return 0;
595: }