Actual source code: bvops.c

slepc-3.16.2 2022-02-01
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2021, 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 operations, except those involving global communication
 12: */

 14: #include <slepc/private/bvimpl.h>
 15: #include <slepcds.h>

 17: /*@
 18:    BVMult - Computes Y = beta*Y + alpha*X*Q.

 20:    Logically Collective on Y

 22:    Input Parameters:
 23: +  Y,X        - basis vectors
 24: .  alpha,beta - scalars
 25: -  Q          - (optional) sequential dense matrix

 27:    Output Parameter:
 28: .  Y          - the modified basis vectors

 30:    Notes:
 31:    X and Y must be different objects. The case X=Y can be addressed with
 32:    BVMultInPlace().

 34:    If matrix Q is NULL, then an AXPY operation Y = beta*Y + alpha*X is done
 35:    (i.e. results as if Q = identity). If provided,
 36:    the matrix Q must be a sequential dense Mat, with all entries equal on
 37:    all processes (otherwise each process will compute a different update).
 38:    The dimensions of Q must be at least m,n where m is the number of active
 39:    columns of X and n is the number of active columns of Y.

 41:    The leading columns of Y are not modified. Also, if X has leading
 42:    columns specified, then these columns do not participate in the computation.
 43:    Hence, only rows (resp. columns) of Q starting from lx (resp. ly) are used,
 44:    where lx (resp. ly) is the number of leading columns of X (resp. Y).

 46:    Level: intermediate

 48: .seealso: BVMultVec(), BVMultColumn(), BVMultInPlace(), BVSetActiveColumns()
 49: @*/
 50: PetscErrorCode BVMult(BV Y,PetscScalar alpha,PetscScalar beta,BV X,Mat Q)
 51: {
 53:   PetscInt       m,n;

 62:   BVCheckSizes(Y,1);
 63:   BVCheckOp(Y,1,mult);
 65:   BVCheckSizes(X,4);
 68:   if (X==Y) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_WRONG,"X and Y arguments must be different");
 69:   if (Q) {
 71:     MatGetSize(Q,&m,&n);
 72:     if (m<X->k) SETERRQ2(PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_SIZ,"Mat argument has %D rows, should have at least %D",m,X->k);
 73:     if (n<Y->k) SETERRQ2(PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_SIZ,"Mat argument has %D columns, should have at least %D",n,Y->k);
 74:   }
 75:   if (X->n!=Y->n) SETERRQ2(PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension X %D, Y %D",X->n,Y->n);

 77:   PetscLogEventBegin(BV_Mult,X,Y,0,0);
 78:   (*Y->ops->mult)(Y,alpha,beta,X,Q);
 79:   PetscLogEventEnd(BV_Mult,X,Y,0,0);
 80:   PetscObjectStateIncrease((PetscObject)Y);
 81:   return(0);
 82: }

 84: /*@
 85:    BVMultVec - Computes y = beta*y + alpha*X*q.

 87:    Logically Collective on X

 89:    Input Parameters:
 90: +  X          - a basis vectors object
 91: .  alpha,beta - scalars
 92: .  y          - a vector
 93: -  q          - an array of scalars

 95:    Output Parameter:
 96: .  y          - the modified vector

 98:    Notes:
 99:    This operation is the analogue of BVMult() but with a BV and a Vec,
100:    instead of two BV. Note that arguments are listed in different order
101:    with respect to BVMult().

103:    If X has leading columns specified, then these columns do not participate
104:    in the computation.

106:    The length of array q must be equal to the number of active columns of X
107:    minus the number of leading columns, i.e. the first entry of q multiplies
108:    the first non-leading column.

110:    Level: intermediate

112: .seealso: BVMult(), BVMultColumn(), BVMultInPlace(), BVSetActiveColumns()
113: @*/
114: PetscErrorCode BVMultVec(BV X,PetscScalar alpha,PetscScalar beta,Vec y,PetscScalar q[])
115: {
117:   PetscInt       n,N;

126:   BVCheckSizes(X,1);
127:   BVCheckOp(X,1,multvec);

131:   VecGetSize(y,&N);
132:   VecGetLocalSize(y,&n);
133:   if (N!=X->N || n!=X->n) SETERRQ4(PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_INCOMP,"Vec sizes (global %D, local %D) do not match BV sizes (global %D, local %D)",N,n,X->N,X->n);

135:   PetscLogEventBegin(BV_MultVec,X,y,0,0);
136:   (*X->ops->multvec)(X,alpha,beta,y,q);
137:   PetscLogEventEnd(BV_MultVec,X,y,0,0);
138:   return(0);
139: }

141: /*@
142:    BVMultColumn - Computes y = beta*y + alpha*X*q, where y is the j-th column
143:    of X.

145:    Logically Collective on X

147:    Input Parameters:
148: +  X          - a basis vectors object
149: .  alpha,beta - scalars
150: .  j          - the column index
151: -  q          - an array of scalars

153:    Notes:
154:    This operation is equivalent to BVMultVec() but it uses column j of X
155:    rather than taking a Vec as an argument. The number of active columns of
156:    X is set to j before the computation, and restored afterwards.
157:    If X has leading columns specified, then these columns do not participate
158:    in the computation. Therefore, the length of array q must be equal to j
159:    minus the number of leading columns.

161:    Developer Notes:
162:    If q is NULL, then the coefficients are taken from position nc+l of the
163:    internal buffer vector, see BVGetBufferVec().

165:    Level: advanced

167: .seealso: BVMult(), BVMultVec(), BVMultInPlace(), BVSetActiveColumns()
168: @*/
169: PetscErrorCode BVMultColumn(BV X,PetscScalar alpha,PetscScalar beta,PetscInt j,PetscScalar *q)
170: {
172:   PetscInt       ksave;
173:   Vec            y;

181:   BVCheckSizes(X,1);

183:   if (j<0) SETERRQ(PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
184:   if (j>=X->m) SETERRQ2(PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_OUTOFRANGE,"Index j=%D but BV only has %D columns",j,X->m);

186:   PetscLogEventBegin(BV_MultVec,X,0,0,0);
187:   ksave = X->k;
188:   X->k = j;
189:   if (!q && !X->buffer) { BVGetBufferVec(X,&X->buffer); }
190:   BVGetColumn(X,j,&y);
191:   (*X->ops->multvec)(X,alpha,beta,y,q);
192:   BVRestoreColumn(X,j,&y);
193:   X->k = ksave;
194:   PetscLogEventEnd(BV_MultVec,X,0,0,0);
195:   PetscObjectStateIncrease((PetscObject)X);
196:   return(0);
197: }

199: /*@
200:    BVMultInPlace - Update a set of vectors as V(:,s:e-1) = V*Q(:,s:e-1).

202:    Logically Collective on V

204:    Input Parameters:
205: +  Q - a sequential dense matrix
206: .  s - first column of V to be overwritten
207: -  e - first column of V not to be overwritten

209:    Input/Output Parameter:
210: .  V - basis vectors

212:    Notes:
213:    The matrix Q must be a sequential dense Mat, with all entries equal on
214:    all processes (otherwise each process will compute a different update).

216:    This function computes V(:,s:e-1) = V*Q(:,s:e-1), that is, given a set of
217:    vectors V, columns from s to e-1 are overwritten with columns from s to
218:    e-1 of the matrix-matrix product V*Q. Only columns s to e-1 of Q are
219:    referenced.

221:    Level: intermediate

223: .seealso: BVMult(), BVMultVec(), BVMultInPlaceHermitianTranspose(), BVSetActiveColumns()
224: @*/
225: PetscErrorCode BVMultInPlace(BV V,Mat Q,PetscInt s,PetscInt e)
226: {
228:   PetscInt       m,n;

236:   BVCheckSizes(V,1);

240:   if (s<V->l || s>V->m) SETERRQ3(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Argument s has wrong value %D, should be between %D and %D",s,V->l,V->m);
241:   if (e<V->l || e>V->m) SETERRQ3(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Argument e has wrong value %D, should be between %D and %D",e,V->l,V->m);
242:   MatGetSize(Q,&m,&n);
243:   if (m<V->k) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"Mat argument has %D rows, should have at least %D",m,V->k);
244:   if (e>n) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"Mat argument only has %D columns, the requested value of e is larger: %D",n,e);
245:   if (s>=e) return(0);

247:   PetscLogEventBegin(BV_MultInPlace,V,Q,0,0);
248:   (*V->ops->multinplace)(V,Q,s,e);
249:   PetscLogEventEnd(BV_MultInPlace,V,Q,0,0);
250:   PetscObjectStateIncrease((PetscObject)V);
251:   return(0);
252: }

254: /*@
255:    BVMultInPlaceHermitianTranspose - Update a set of vectors as V(:,s:e-1) = V*Q'(:,s:e-1).

257:    Logically Collective on V

259:    Input Parameters:
260: +  Q - a sequential dense matrix
261: .  s - first column of V to be overwritten
262: -  e - first column of V not to be overwritten

264:    Input/Output Parameter:
265: .  V - basis vectors

267:    Notes:
268:    This is a variant of BVMultInPlace() where the conjugate transpose
269:    of Q is used.

271:    Level: intermediate

273: .seealso: BVMultInPlace()
274: @*/
275: PetscErrorCode BVMultInPlaceHermitianTranspose(BV V,Mat Q,PetscInt s,PetscInt e)
276: {
278:   PetscInt       m,n;

286:   BVCheckSizes(V,1);

290:   if (s<V->l || s>V->m) SETERRQ3(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Argument s has wrong value %D, should be between %D and %D",s,V->l,V->m);
291:   if (e<V->l || e>V->m) SETERRQ3(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Argument e has wrong value %D, should be between %D and %D",e,V->l,V->m);
292:   MatGetSize(Q,&m,&n);
293:   if (n<V->k) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"Mat argument has %D columns, should have at least %D",n,V->k);
294:   if (e>m) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"Mat argument only has %D rows, the requested value of e is larger: %D",m,e);
295:   if (s>=e || !V->n) return(0);

297:   PetscLogEventBegin(BV_MultInPlace,V,Q,0,0);
298:   (*V->ops->multinplacetrans)(V,Q,s,e);
299:   PetscLogEventEnd(BV_MultInPlace,V,Q,0,0);
300:   PetscObjectStateIncrease((PetscObject)V);
301:   return(0);
302: }

304: /*@
305:    BVScale - Multiply the BV entries by a scalar value.

307:    Logically Collective on bv

309:    Input Parameters:
310: +  bv    - basis vectors
311: -  alpha - scaling factor

313:    Note:
314:    All active columns (except the leading ones) are scaled.

316:    Level: intermediate

318: .seealso: BVScaleColumn(), BVSetActiveColumns()
319: @*/
320: PetscErrorCode BVScale(BV bv,PetscScalar alpha)
321: {

328:   BVCheckSizes(bv,1);
329:   if (alpha == (PetscScalar)1.0) return(0);

331:   PetscLogEventBegin(BV_Scale,bv,0,0,0);
332:   if (bv->n) {
333:     (*bv->ops->scale)(bv,-1,alpha);
334:   }
335:   PetscLogEventEnd(BV_Scale,bv,0,0,0);
336:   PetscObjectStateIncrease((PetscObject)bv);
337:   return(0);
338: }

340: /*@
341:    BVScaleColumn - Scale one column of a BV.

343:    Logically Collective on bv

345:    Input Parameters:
346: +  bv    - basis vectors
347: .  j     - column number to be scaled
348: -  alpha - scaling factor

350:    Level: intermediate

352: .seealso: BVScale(), BVSetActiveColumns()
353: @*/
354: PetscErrorCode BVScaleColumn(BV bv,PetscInt j,PetscScalar alpha)
355: {

363:   BVCheckSizes(bv,1);

365:   if (j<0 || j>=bv->m) SETERRQ2(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Argument j has wrong value %D, the number of columns is %D",j,bv->m);
366:   if (alpha == (PetscScalar)1.0) return(0);

368:   PetscLogEventBegin(BV_Scale,bv,0,0,0);
369:   if (bv->n) {
370:     (*bv->ops->scale)(bv,j,alpha);
371:   }
372:   PetscLogEventEnd(BV_Scale,bv,0,0,0);
373:   PetscObjectStateIncrease((PetscObject)bv);
374:   return(0);
375: }

377: PETSC_STATIC_INLINE PetscErrorCode BVSetRandomColumn_Private(BV bv,PetscInt k)
378: {
380:   PetscInt       i,low,high;
381:   PetscScalar    *px,t;
382:   Vec            x;

385:   BVGetColumn(bv,k,&x);
386:   if (bv->rrandom) {  /* generate the same vector irrespective of number of processes */
387:     VecGetOwnershipRange(x,&low,&high);
388:     VecGetArray(x,&px);
389:     for (i=0;i<bv->N;i++) {
390:       PetscRandomGetValue(bv->rand,&t);
391:       if (i>=low && i<high) px[i-low] = t;
392:     }
393:     VecRestoreArray(x,&px);
394:   } else {
395:     VecSetRandom(x,bv->rand);
396:   }
397:   BVRestoreColumn(bv,k,&x);
398:   return(0);
399: }

401: PETSC_STATIC_INLINE PetscErrorCode BVSetRandomNormalColumn_Private(BV bv,PetscInt k,Vec w1,Vec w2)
402: {
404:   PetscInt       i,low,high;
405:   PetscScalar    *px,s,t;
406:   Vec            x;

409:   BVGetColumn(bv,k,&x);
410:   if (bv->rrandom) {  /* generate the same vector irrespective of number of processes */
411:     VecGetOwnershipRange(x,&low,&high);
412:     VecGetArray(x,&px);
413:     for (i=0;i<bv->N;i++) {
414:       PetscRandomGetValue(bv->rand,&s);
415:       PetscRandomGetValue(bv->rand,&t);
416:       if (i>=low && i<high) {
417: #if defined(PETSC_USE_COMPLEX)
418:         px[i-low] = PetscCMPLX(PetscSqrtReal(-2.0*PetscLogReal(PetscRealPart(s)))*PetscCosReal(2.0*PETSC_PI*PetscRealPart(t)),PetscSqrtReal(-2.0*PetscLogReal(PetscImaginaryPart(s)))*PetscCosReal(2.0*PETSC_PI*PetscImaginaryPart(t)));
419: #else
420:         px[i-low] = PetscSqrtReal(-2.0*PetscLogReal(s))*PetscCosReal(2.0*PETSC_PI*t);
421: #endif
422:       }
423:     }
424:     VecRestoreArray(x,&px);
425:   } else {
426:     VecSetRandomNormal(x,bv->rand,w1,w2);
427:   }
428:   BVRestoreColumn(bv,k,&x);
429:   return(0);
430: }

432: PETSC_STATIC_INLINE PetscErrorCode BVSetRandomSignColumn_Private(BV bv,PetscInt k)
433: {
435:   PetscInt       i,low,high;
436:   PetscScalar    *px,t;
437:   Vec            x;

440:   BVGetColumn(bv,k,&x);
441:   VecGetOwnershipRange(x,&low,&high);
442:   if (bv->rrandom) {  /* generate the same vector irrespective of number of processes */
443:     VecGetArray(x,&px);
444:     for (i=0;i<bv->N;i++) {
445:       PetscRandomGetValue(bv->rand,&t);
446:       if (i>=low && i<high) px[i-low] = (PetscRealPart(t)<0.5)? -1.0: 1.0;
447:     }
448:     VecRestoreArray(x,&px);
449:   } else {
450:     VecSetRandom(x,bv->rand);
451:     VecGetArray(x,&px);
452:     for (i=low;i<high;i++) {
453:       px[i-low] = (PetscRealPart(px[i-low])<0.5)? -1.0: 1.0;
454:     }
455:     VecRestoreArray(x,&px);
456:   }
457:   BVRestoreColumn(bv,k,&x);
458:   return(0);
459: }

461: /*@
462:    BVSetRandom - Set the columns of a BV to random numbers.

464:    Logically Collective on bv

466:    Input Parameters:
467: .  bv - basis vectors

469:    Note:
470:    All active columns (except the leading ones) are modified.

472:    Level: advanced

474: .seealso: BVSetRandomContext(), BVSetRandomColumn(), BVSetRandomNormal(), BVSetRandomSign(), BVSetRandomCond(), BVSetActiveColumns()
475: @*/
476: PetscErrorCode BVSetRandom(BV bv)
477: {
479:   PetscInt       k;

484:   BVCheckSizes(bv,1);

486:   BVGetRandomContext(bv,&bv->rand);
487:   PetscLogEventBegin(BV_SetRandom,bv,0,0,0);
488:   for (k=bv->l;k<bv->k;k++) {
489:     BVSetRandomColumn_Private(bv,k);
490:   }
491:   PetscLogEventEnd(BV_SetRandom,bv,0,0,0);
492:   PetscObjectStateIncrease((PetscObject)bv);
493:   return(0);
494: }

496: /*@
497:    BVSetRandomColumn - Set one column of a BV to random numbers.

499:    Logically Collective on bv

501:    Input Parameters:
502: +  bv - basis vectors
503: -  j  - column number to be set

505:    Level: advanced

507: .seealso: BVSetRandomContext(), BVSetRandom(), BVSetRandomNormal(), BVSetRandomCond()
508: @*/
509: PetscErrorCode BVSetRandomColumn(BV bv,PetscInt j)
510: {

517:   BVCheckSizes(bv,1);
518:   if (j<0 || j>=bv->m) SETERRQ2(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Argument j has wrong value %D, the number of columns is %D",j,bv->m);

520:   BVGetRandomContext(bv,&bv->rand);
521:   PetscLogEventBegin(BV_SetRandom,bv,0,0,0);
522:   BVSetRandomColumn_Private(bv,j);
523:   PetscLogEventEnd(BV_SetRandom,bv,0,0,0);
524:   PetscObjectStateIncrease((PetscObject)bv);
525:   return(0);
526: }

528: /*@
529:    BVSetRandomNormal - Set the columns of a BV to random numbers with a normal
530:    distribution.

532:    Logically Collective on bv

534:    Input Parameter:
535: .  bv - basis vectors

537:    Notes:
538:    All active columns (except the leading ones) are modified.

540:    Other functions such as BVSetRandom(), BVSetRandomColumn(), and BVSetRandomCond()
541:    produce random numbers with a uniform distribution. This function returns values
542:    that fit a normal distribution (Gaussian).

544:    Developer Notes:
545:    The current implementation obtains each of the columns by applying the Box-Muller
546:    transform on two random vectors with uniformly distributed entries.

548:    Level: advanced

550: .seealso: BVSetRandomContext(), BVSetRandom(), BVSetRandomColumn(), BVSetRandomCond(), BVSetActiveColumns()
551: @*/
552: PetscErrorCode BVSetRandomNormal(BV bv)
553: {
555:   PetscInt       k;
556:   Vec            w1=NULL,w2=NULL;

561:   BVCheckSizes(bv,1);

563:   BVGetRandomContext(bv,&bv->rand);
564:   if (!bv->rrandom) {
565:     BVCreateVec(bv,&w1);
566:     BVCreateVec(bv,&w2);
567:   }
568:   PetscLogEventBegin(BV_SetRandom,bv,0,0,0);
569:   for (k=bv->l;k<bv->k;k++) {
570:     BVSetRandomNormalColumn_Private(bv,k,w1,w2);
571:   }
572:   PetscLogEventEnd(BV_SetRandom,bv,0,0,0);
573:   if (!bv->rrandom) {
574:     VecDestroy(&w1);
575:     VecDestroy(&w2);
576:   }
577:   PetscObjectStateIncrease((PetscObject)bv);
578:   return(0);
579: }

581: /*@
582:    BVSetRandomSign - Set the entries of a BV to values 1 or -1 with equal
583:    probability.

585:    Logically Collective on bv

587:    Input Parameter:
588: .  bv - basis vectors

590:    Notes:
591:    All active columns (except the leading ones) are modified.

593:    This function is used, e.g., in contour integral methods when estimating
594:    the number of eigenvalues enclosed by the contour via an unbiased
595:    estimator of tr(f(A)) [Bai et al., JCAM 1996].

597:    Developer Notes:
598:    The current implementation obtains random numbers and then replaces them
599:    with -1 or 1 depending on the value being less than 0.5 or not.

601:    Level: advanced

603: .seealso: BVSetRandomContext(), BVSetRandom(), BVSetRandomColumn(), BVSetActiveColumns()
604: @*/
605: PetscErrorCode BVSetRandomSign(BV bv)
606: {
608:   PetscScalar    low,high;
609:   PetscInt       k;

614:   BVCheckSizes(bv,1);

616:   BVGetRandomContext(bv,&bv->rand);
617:   PetscRandomGetInterval(bv->rand,&low,&high);
618:   if (PetscRealPart(low)!=0.0 || PetscRealPart(high)!=1.0) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"The PetscRandom object in the BV must have interval [0,1]");
619:   PetscLogEventBegin(BV_SetRandom,bv,0,0,0);
620:   for (k=bv->l;k<bv->k;k++) {
621:     BVSetRandomSignColumn_Private(bv,k);
622:   }
623:   PetscLogEventEnd(BV_SetRandom,bv,0,0,0);
624:   PetscObjectStateIncrease((PetscObject)bv);
625:   return(0);
626: }

628: /*@
629:    BVSetRandomCond - Set the columns of a BV to random numbers, in a way that
630:    the generated matrix has a given condition number.

632:    Logically Collective on bv

634:    Input Parameters:
635: +  bv    - basis vectors
636: -  condn - condition number

638:    Note:
639:    All active columns (except the leading ones) are modified.

641:    Level: advanced

643: .seealso: BVSetRandomContext(), BVSetRandom(), BVSetRandomColumn(), BVSetRandomNormal(), BVSetActiveColumns()
644: @*/
645: PetscErrorCode BVSetRandomCond(BV bv,PetscReal condn)
646: {
648:   PetscInt       k,i;
649:   PetscScalar    *eig,*d;
650:   DS             ds;
651:   Mat            A,X,Xt,M,G;

656:   BVCheckSizes(bv,1);

658:   BVGetRandomContext(bv,&bv->rand);
659:   PetscLogEventBegin(BV_SetRandom,bv,0,0,0);
660:   /* B = rand(n,k) */
661:   for (k=bv->l;k<bv->k;k++) {
662:     BVSetRandomColumn_Private(bv,k);
663:   }
664:   DSCreate(PetscObjectComm((PetscObject)bv),&ds);
665:   DSSetType(ds,DSHEP);
666:   DSAllocate(ds,bv->m);
667:   DSSetDimensions(ds,bv->k,bv->l,bv->k);
668:   /* [V,S] = eig(B'*B) */
669:   DSGetMat(ds,DS_MAT_A,&A);
670:   BVDot(bv,bv,A);
671:   DSRestoreMat(ds,DS_MAT_A,&A);
672:   PetscMalloc1(bv->m,&eig);
673:   DSSolve(ds,eig,NULL);
674:   DSSynchronize(ds,eig,NULL);
675:   DSVectors(ds,DS_MAT_X,NULL,NULL);
676:   /* M = diag(linspace(1/condn,1,n)./sqrt(diag(S)))' */
677:   MatCreateSeqDense(PETSC_COMM_SELF,bv->k,bv->k,NULL,&M);
678:   MatZeroEntries(M);
679:   MatDenseGetArray(M,&d);
680:   for (i=0;i<bv->k;i++) d[i+i*bv->m] = (1.0/condn+(1.0-1.0/condn)/(bv->k-1)*i)/PetscSqrtScalar(eig[i]);
681:   MatDenseRestoreArray(M,&d);
682:   /* G = X*M*X' */
683:   MatCreateSeqDense(PETSC_COMM_SELF,bv->k,bv->k,NULL,&Xt);
684:   DSGetMat(ds,DS_MAT_X,&X);
685:   MatTranspose(X,MAT_REUSE_MATRIX,&Xt);
686:   MatProductCreate(Xt,M,NULL,&G);
687:   MatProductSetType(G,MATPRODUCT_PtAP);
688:   MatProductSetFromOptions(G);
689:   MatProductSymbolic(G);
690:   MatProductNumeric(G);
691:   MatProductClear(G);
692:   MatDestroy(&X);
693:   MatDestroy(&Xt);
694:   MatDestroy(&M);
695:   /* B = B*G */
696:   BVMultInPlace(bv,G,bv->l,bv->k);
697:   MatDestroy(&G);
698:   PetscFree(eig);
699:   DSDestroy(&ds);
700:   PetscLogEventEnd(BV_SetRandom,bv,0,0,0);
701:   PetscObjectStateIncrease((PetscObject)bv);
702:   return(0);
703: }

705: /*@
706:    BVMatMult - Computes the matrix-vector product for each column, Y=A*V.

708:    Neighbor-wise Collective on A

710:    Input Parameters:
711: +  V - basis vectors context
712: -  A - the matrix

714:    Output Parameter:
715: .  Y - the result

717:    Notes:
718:    Both V and Y must be distributed in the same manner. Only active columns
719:    (excluding the leading ones) are processed.
720:    In the result Y, columns are overwritten starting from the leading ones.
721:    The number of active columns in V and Y should match, although they need
722:    not be the same columns.

724:    It is possible to choose whether the computation is done column by column
725:    or as a Mat-Mat product, see BVSetMatMultMethod().

727:    Level: beginner

729: .seealso: BVCopy(), BVSetActiveColumns(), BVMatMultColumn(), BVMatMultTranspose(), BVMatMultHermitianTranspose(), BVSetMatMultMethod()
730: @*/
731: PetscErrorCode BVMatMult(BV V,Mat A,BV Y)
732: {
734:   PetscInt       M,N,m,n;

739:   BVCheckSizes(V,1);
740:   BVCheckOp(V,1,matmult);
745:   BVCheckSizes(Y,3);

749:   MatGetSize(A,&M,&N);
750:   MatGetLocalSize(A,&m,&n);
751:   if (M!=Y->N) SETERRQ2(PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Mismatching row dimension A %D, Y %D",M,Y->N);
752:   if (m!=Y->n) SETERRQ2(PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Mismatching local row dimension A %D, Y %D",m,Y->n);
753:   if (N!=V->N) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching column dimension A %D, V %D",N,V->N);
754:   if (n!=V->n) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching local column dimension A %D, V %D",n,V->n);
755:   if (V->k-V->l!=Y->k-Y->l) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"Y has %D active columns, should match %D active columns in V",Y->k-Y->l,V->k-V->l);

757:   PetscLogEventBegin(BV_MatMult,V,A,Y,0);
758:   (*V->ops->matmult)(V,A,Y);
759:   PetscLogEventEnd(BV_MatMult,V,A,Y,0);
760:   PetscObjectStateIncrease((PetscObject)Y);
761:   return(0);
762: }

764: /*@
765:    BVMatMultTranspose - Computes the matrix-vector product with the transpose
766:    of a matrix for each column, Y=A^T*V.

768:    Neighbor-wise Collective on A

770:    Input Parameters:
771: +  V - basis vectors context
772: -  A - the matrix

774:    Output Parameter:
775: .  Y - the result

777:    Notes:
778:    Both V and Y must be distributed in the same manner. Only active columns
779:    (excluding the leading ones) are processed.
780:    In the result Y, columns are overwritten starting from the leading ones.
781:    The number of active columns in V and Y should match, although they need
782:    not be the same columns.

784:    Currently implemented via MatCreateTranspose().

786:    Level: beginner

788: .seealso: BVMatMult(), BVMatMultHermitianTranspose()
789: @*/
790: PetscErrorCode BVMatMultTranspose(BV V,Mat A,BV Y)
791: {
793:   PetscInt       M,N,m,n;
794:   Mat            AT;

799:   BVCheckSizes(V,1);
804:   BVCheckSizes(Y,3);

808:   MatGetSize(A,&M,&N);
809:   MatGetLocalSize(A,&m,&n);
810:   if (M!=V->N) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching row dimension A %D, V %D",M,V->N);
811:   if (m!=V->n) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching local row dimension A %D, V %D",m,V->n);
812:   if (N!=Y->N) SETERRQ2(PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Mismatching column dimension A %D, Y %D",N,Y->N);
813:   if (n!=Y->n) SETERRQ2(PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Mismatching local column dimension A %D, Y %D",n,Y->n);
814:   if (V->k-V->l!=Y->k-Y->l) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"Y has %D active columns, should match %D active columns in V",Y->k-Y->l,V->k-V->l);

816:   MatCreateTranspose(A,&AT);
817:   BVMatMult(V,AT,Y);
818:   MatDestroy(&AT);
819:   return(0);
820: }

822: /*@
823:    BVMatMultHermitianTranspose - Computes the matrix-vector product with the
824:    conjugate transpose of a matrix for each column, Y=A^H*V.

826:    Neighbor-wise Collective on A

828:    Input Parameters:
829: +  V - basis vectors context
830: -  A - the matrix

832:    Output Parameter:
833: .  Y - the result

835:    Note:
836:    Both V and Y must be distributed in the same manner. Only active columns
837:    (excluding the leading ones) are processed.
838:    In the result Y, columns are overwritten starting from the leading ones.
839:    The number of active columns in V and Y should match, although they need
840:    not be the same columns.

842:    Currently implemented via MatCreateHermitianTranspose().

844:    Level: beginner

846: .seealso: BVMatMult(), BVMatMultTranspose()
847: @*/
848: PetscErrorCode BVMatMultHermitianTranspose(BV V,Mat A,BV Y)
849: {
851:   PetscInt       M,N,m,n;
852:   Mat            AH;

857:   BVCheckSizes(V,1);
862:   BVCheckSizes(Y,3);

866:   MatGetSize(A,&M,&N);
867:   MatGetLocalSize(A,&m,&n);
868:   if (M!=V->N) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching row dimension A %D, V %D",M,V->N);
869:   if (m!=V->n) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching local row dimension A %D, V %D",m,V->n);
870:   if (N!=Y->N) SETERRQ2(PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Mismatching column dimension A %D, Y %D",N,Y->N);
871:   if (n!=Y->n) SETERRQ2(PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Mismatching local column dimension A %D, Y %D",n,Y->n);
872:   if (V->k-V->l!=Y->k-Y->l) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"Y has %D active columns, should match %D active columns in V",Y->k-Y->l,V->k-V->l);

874:   MatCreateHermitianTranspose(A,&AH);
875:   BVMatMult(V,AH,Y);
876:   MatDestroy(&AH);
877:   return(0);
878: }

880: /*@
881:    BVMatMultColumn - Computes the matrix-vector product for a specified
882:    column, storing the result in the next column: v_{j+1}=A*v_j.

884:    Neighbor-wise Collective on A

886:    Input Parameters:
887: +  V - basis vectors context
888: .  A - the matrix
889: -  j - the column

891:    Level: beginner

893: .seealso: BVMatMult(), BVMatMultTransposeColumn(), BVMatMultHermitianTransposeColumn()
894: @*/
895: PetscErrorCode BVMatMultColumn(BV V,Mat A,PetscInt j)
896: {
898:   Vec            vj,vj1;

903:   BVCheckSizes(V,1);
907:   if (j<0) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
908:   if (j+1>=V->m) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Result should go in index j+1=%D but BV only has %D columns",j+1,V->m);

910:   PetscLogEventBegin(BV_MatMultVec,V,A,0,0);
911:   BVGetColumn(V,j,&vj);
912:   BVGetColumn(V,j+1,&vj1);
913:   MatMult(A,vj,vj1);
914:   BVRestoreColumn(V,j,&vj);
915:   BVRestoreColumn(V,j+1,&vj1);
916:   PetscLogEventEnd(BV_MatMultVec,V,A,0,0);
917:   PetscObjectStateIncrease((PetscObject)V);
918:   return(0);
919: }

921: /*@
922:    BVMatMultTransposeColumn - Computes the transpose matrix-vector product for a
923:    specified column, storing the result in the next column: v_{j+1}=A^T*v_j.

925:    Neighbor-wise Collective on A

927:    Input Parameters:
928: +  V - basis vectors context
929: .  A - the matrix
930: -  j - the column

932:    Level: beginner

934: .seealso: BVMatMult(), BVMatMultColumn(), BVMatMultHermitianTransposeColumn()
935: @*/
936: PetscErrorCode BVMatMultTransposeColumn(BV V,Mat A,PetscInt j)
937: {
939:   Vec            vj,vj1;

944:   BVCheckSizes(V,1);
948:   if (j<0) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
949:   if (j+1>=V->m) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Result should go in index j+1=%D but BV only has %D columns",j+1,V->m);

951:   PetscLogEventBegin(BV_MatMultVec,V,A,0,0);
952:   BVGetColumn(V,j,&vj);
953:   BVGetColumn(V,j+1,&vj1);
954:   MatMultTranspose(A,vj,vj1);
955:   BVRestoreColumn(V,j,&vj);
956:   BVRestoreColumn(V,j+1,&vj1);
957:   PetscLogEventEnd(BV_MatMultVec,V,A,0,0);
958:   PetscObjectStateIncrease((PetscObject)V);
959:   return(0);
960: }

962: /*@
963:    BVMatMultHermitianTransposeColumn - Computes the conjugate-transpose matrix-vector
964:    product for a specified column, storing the result in the next column: v_{j+1}=A^H*v_j.

966:    Neighbor-wise Collective on A

968:    Input Parameters:
969: +  V - basis vectors context
970: .  A - the matrix
971: -  j - the column

973:    Level: beginner

975: .seealso: BVMatMult(), BVMatMultColumn(), BVMatMultTransposeColumn()
976: @*/
977: PetscErrorCode BVMatMultHermitianTransposeColumn(BV V,Mat A,PetscInt j)
978: {
980:   Vec            vj,vj1;

985:   BVCheckSizes(V,1);
989:   if (j<0) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
990:   if (j+1>=V->m) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Result should go in index j+1=%D but BV only has %D columns",j+1,V->m);

992:   PetscLogEventBegin(BV_MatMultVec,V,A,0,0);
993:   BVGetColumn(V,j,&vj);
994:   BVGetColumn(V,j+1,&vj1);
995:   MatMultHermitianTranspose(A,vj,vj1);
996:   BVRestoreColumn(V,j,&vj);
997:   BVRestoreColumn(V,j+1,&vj1);
998:   PetscLogEventEnd(BV_MatMultVec,V,A,0,0);
999:   PetscObjectStateIncrease((PetscObject)V);
1000:   return(0);
1001: }