Actual source code: bvcontour.c

slepc-3.17.0 2022-03-31
Report Typos and Errors
  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:   PetscFunctionReturn(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:   PetscFunctionReturn(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:   PetscFunctionReturn(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:   PetscFunctionReturn(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) {
329:       PetscStackCallBLAS("BLASgemm",BLASgemm_("C","N",&l,&n,&m,&alpha,sarray,&lda,sarray,&ldb,&beta,pA,&ldc));
330:     } else {
331:       PetscStackCallBLAS("BLASgemm",BLASgemm_("C","N",&l,&n,&m,&alpha,Q1,&lda,Q1,&ldb,&beta,pA,&ldc));
332:     }
333:     PetscArrayzero(temp2,ml*ml);
334:     PetscMPIIntCast(ml*ml,&len);
335:     MPIU_Allreduce(pA,temp2,len,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)S));

337:     PetscBLASIntCast(ml,&m);
338:     n = m; lda = m; lwork = 5*m, ldu = 1; ldvt = 1;
339: #if defined(PETSC_USE_COMPLEX)
340:     PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("O","N",&m,&n,temp2,&lda,sigma,NULL,&ldu,NULL,&ldvt,work,&lwork,rwork,&info));
341: #else
342:     PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("O","N",&m,&n,temp2,&lda,sigma,NULL,&ldu,NULL,&ldvt,work,&lwork,&info));
343: #endif
344:     SlepcCheckLapackInfo("gesvd",info);

346:     PetscBLASIntCast(S->n,&l);
347:     PetscBLASIntCast(ml,&n);
348:     m = n; lda = l; ldb = m; ldc = l;
349:     if (!k) {
350:       PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&l,&n,&m,&alpha,sarray,&lda,temp2,&ldb,&beta,Q1,&ldc));
351:     } else {
352:       PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&l,&n,&m,&alpha,Q1,&lda,temp2,&ldb,&beta,Q2,&ldc));
353:     }

355:     PetscBLASIntCast(ml,&l);
356:     m = l; n = l; lda = l; ldb = m; ldc = l;
357:     PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&l,&n,&m,&alpha,B,&lda,temp2,&ldb,&beta,tempB,&ldc));
358:     for (i=0;i<ml;i++) {
359:       sigma[i] = PetscSqrtReal(sigma[i]);
360:       for (j=0;j<S->n;j++) {
361:         if (k%2) Q2[j+i*S->n] /= sigma[i];
362:         else Q1[j+i*S->n] /= sigma[i];
363:       }
364:       for (j=0;j<ml;j++) B[j+i*ml] = tempB[j+i*ml]*sigma[i];
365:     }
366:   }

368:   PetscBLASIntCast(ml,&m);
369:   n = m; lda = m; ldu=1; ldvt=1;
370: #if defined(PETSC_USE_COMPLEX)
371:   PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","O",&m,&n,B,&lda,sigma,NULL,&ldu,NULL,&ldvt,work,&lwork,rwork,&info));
372: #else
373:   PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","O",&m,&n,B,&lda,sigma,NULL,&ldu,NULL,&ldvt,work,&lwork,&info));
374: #endif
375:   SlepcCheckLapackInfo("gesvd",info);

377:   PetscBLASIntCast(S->n,&l);
378:   PetscBLASIntCast(ml,&n);
379:   m = n; lda = l; ldb = m; ldc = l;
380:   if (k%2) {
381:     PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&l,&n,&m,&alpha,Q1,&lda,B,&ldb,&beta,sarray,&ldc));
382:   } else {
383:     PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&l,&n,&m,&alpha,Q2,&lda,B,&ldb,&beta,sarray,&ldc));
384:   }

386:   PetscFPTrapPop();
387:   BVRestoreArray(S,&sarray);

389:   if (rank) {
390:     (*rank) = 0;
391:     for (i=0;i<ml;i++) {
392:       if (sigma[i]/PetscMax(sigma[0],1.0)>delta) (*rank)++;
393:     }
394:   }
395:   PetscFree6(temp2,Q1,Q2,B,tempB,work);
396: #if defined(PETSC_USE_COMPLEX)
397:   PetscFree(rwork);
398: #endif
399:   PetscFunctionReturn(0);
400: }

402: PetscErrorCode BVSVDAndRank_QR(BV S,PetscReal delta,PetscScalar *pA,PetscReal *sigma,PetscInt *rank)
403: {
404:   PetscInt       i,n,ml=S->k;
405:   PetscBLASInt   m,lda,lwork,info;
406:   PetscScalar    *work;
407:   PetscReal      *rwork;
408:   Mat            A;
409:   Vec            v;

411:   /* Compute QR factorizaton of S */
412:   BVGetSizes(S,NULL,&n,NULL);
413:   n    = PetscMin(n,ml);
414:   BVSetActiveColumns(S,0,n);
415:   PetscArrayzero(pA,ml*n);
416:   MatCreateDense(PETSC_COMM_SELF,n,n,PETSC_DECIDE,PETSC_DECIDE,pA,&A);
417:   BVOrthogonalize(S,A);
418:   if (n<ml) {
419:     /* the rest of the factorization */
420:     for (i=n;i<ml;i++) {
421:       BVGetColumn(S,i,&v);
422:       BVOrthogonalizeVec(S,v,pA+i*n,NULL,NULL);
423:       BVRestoreColumn(S,i,&v);
424:     }
425:   }
426:   PetscBLASIntCast(n,&lda);
427:   PetscBLASIntCast(ml,&m);
428:   PetscMalloc2(5*ml,&work,5*ml,&rwork);
429:   lwork = 5*m;
430:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
431: #if !defined (PETSC_USE_COMPLEX)
432:   PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("O","N",&lda,&m,pA,&lda,sigma,NULL,&lda,NULL,&lda,work,&lwork,&info));
433: #else
434:   PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("O","N",&lda,&m,pA,&lda,sigma,NULL,&lda,NULL,&lda,work,&lwork,rwork,&info));
435: #endif
436:   SlepcCheckLapackInfo("gesvd",info);
437:   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:   BVMultInPlace(S,A,0,*rank);
444:   PetscFree2(work,rwork);
445:   MatDestroy(&A);
446:   PetscFunctionReturn(0);
447: }

449: 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:   /* Compute QR factorizaton of S */
458:   BVGetSizes(S,NULL,&n,NULL);
460:   BVSetActiveColumns(S,0,ml);
461:   PetscArrayzero(pA,ml*ml);
462:   MatCreateDense(PETSC_COMM_SELF,ml,ml,PETSC_DECIDE,PETSC_DECIDE,pA,&A);
463:   BVOrthogonalize(S,A);
464:   MatDestroy(&A);

466:   /* SVD of first (M-1)*L diagonal block */
467:   PetscBLASIntCast((M-1)*L,&m);
468:   PetscMalloc5(m*m,&T,m*m,&R,m*m,&U,5*ml,&work,5*ml,&rwork);
469:   for (j=0;j<m;j++) PetscArraycpy(R+j*m,pA+j*ml,m);
470:   lwork = 5*m;
471:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
472: #if !defined (PETSC_USE_COMPLEX)
473:   PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("S","O",&m,&m,R,&m,sigma,U,&m,NULL,&m,work,&lwork,&info));
474: #else
475:   PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("S","O",&m,&m,R,&m,sigma,U,&m,NULL,&m,work,&lwork,rwork,&info));
476: #endif
477:   SlepcCheckLapackInfo("gesvd",info);
478:   PetscFPTrapPop();
479:   *rank = 0;
480:   for (i=0;i<m;i++) {
481:     if (sigma[i]/PetscMax(sigma[0],1)>delta) (*rank)++;
482:   }
483:   MatCreateDense(PETSC_COMM_SELF,m,m,PETSC_DECIDE,PETSC_DECIDE,U,&A);
484:   BVSetActiveColumns(S,0,m);
485:   BVMultInPlace(S,A,0,*rank);
486:   MatDestroy(&A);
487:   /* Projected linear system */
488:   /* m first columns of A have the right singular vectors */
489:   PetscBLASIntCast(*rank,&k_);
490:   PetscBLASIntCast(ml,&lda);
491:   PetscStackCallBLAS("BLASgemm",BLASgemm_("N","C",&m,&k_,&m,&sone,pA+L*lda,&lda,R,&m,&zero,T,&m));
492:   PetscArrayzero(pA,ml*ml);
493:   PetscStackCallBLAS("BLASgemm",BLASgemm_("C","N",&k_,&k_,&m,&sone,U,&m,T,&m,&zero,pA,&k_));
494:   for (j=0;j<k_;j++) for (i=0;i<k_;i++) pA[j*k_+i] /= sigma[j];
495:   PetscFree5(T,R,U,work,rwork);
496:   PetscFunctionReturn(0);
497: }

499: /*@
500:    BVSVDAndRank - Compute the SVD (left singular vectors only, and singular
501:    values) and determine the numerical rank according to a tolerance.

503:    Collective on S

505:    Input Parameters:
506: +  S     - the basis vectors
507: .  m     - the moment degree
508: .  l     - the block size
509: .  delta - the tolerance used to determine the rank
510: -  meth  - the method to be used

512:    Output Parameters:
513: +  A     - workspace, on output contains relevant values in the CAA method
514: .  sigma - computed singular values
515: -  rank  - estimated rank (optional)

517:    Notes:
518:    This function computes [U,Sigma,V] = svd(S) and replaces S with U.
519:    The current implementation computes this via S'*S, and it may include
520:    some kind of iterative refinement to improve accuracy in some cases.

522:    The parameters m and l refer to the moment and block size of contour
523:    integral methods. All columns up to m*l are modified, and the active
524:    columns are set to 0..m*l.

526:    The method is one of BV_SVD_METHOD_REFINE, BV_SVD_METHOD_QR, BV_SVD_METHOD_QR_CAA.

528:    The A workspace should be m*l*m*l in size.

530:    Once the decomposition is computed, the numerical rank is estimated
531:    by counting the number of singular values that are larger than the
532:    tolerance delta, relative to the first singular value.

534:    Level: developer

536: .seealso: BVSetActiveColumns()
537: @*/
538: PetscErrorCode BVSVDAndRank(BV S,PetscInt m,PetscInt l,PetscReal delta,BVSVDMethod meth,PetscScalar *A,PetscReal *sigma,PetscInt *rank)
539: {

549:   PetscLogEventBegin(BV_SVDAndRank,S,0,0,0);
550:   BVSetActiveColumns(S,0,m*l);
551:   switch (meth) {
552:     case BV_SVD_METHOD_REFINE:
553:       BVSVDAndRank_Refine(S,delta,A,sigma,rank);
554:       break;
555:     case BV_SVD_METHOD_QR:
556:       BVSVDAndRank_QR(S,delta,A,sigma,rank);
557:       break;
558:     case BV_SVD_METHOD_QR_CAA:
559:       BVSVDAndRank_QR_CAA(S,m,l,delta,A,sigma,rank);
560:       break;
561:   }
562:   PetscLogEventEnd(BV_SVDAndRank,S,0,0,0);
563:   PetscFunctionReturn(0);
564: }

566: /*@
567:    BVCISSResizeBases - Resize the bases involved in CISS solvers when the L grows.

569:    Collective on S

571:    Input Parameters:
572: +  S      - basis of L*M columns
573: .  V      - basis of L columns (may be associated to subcommunicators)
574: .  Y      - basis of npoints*L columns
575: .  Lold   - old value of L
576: .  Lnew   - new value of L
577: .  M      - the moment size
578: -  npoints - number of integration points

580:    Level: developer

582: .seealso: BVResize()
583: @*/
584: PetscErrorCode BVCISSResizeBases(BV S,BV V,BV Y,PetscInt Lold,PetscInt Lnew,PetscInt M,PetscInt npoints)
585: {
586:   PetscInt       i,j;


596:   BVResize(S,Lnew*M,PETSC_FALSE);
597:   BVResize(V,Lnew,PETSC_TRUE);
598:   BVResize(Y,Lnew*npoints,PETSC_TRUE);
599:   /* columns of Y are interleaved */
600:   for (i=npoints-1;i>=0;i--) {
601:     for (j=Lold-1;j>=0;j--) BVCopyColumn(Y,i*Lold+j,i*Lnew+j);
602:   }
603:   PetscFunctionReturn(0);
604: }