Actual source code: bvops.c

slepc-3.14.1 2020-12-08
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2020, 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:   PetscBool      match;
 54:   PetscInt       m,n;

 63:   BVCheckSizes(Y,1);
 64:   BVCheckOp(Y,1,mult);
 66:   BVCheckSizes(X,4);
 69:   if (X==Y) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_WRONG,"X and Y arguments must be different");
 70:   if (Q) {
 71:     PetscObjectTypeCompare((PetscObject)Q,MATSEQDENSE,&match);
 72:     if (!match) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_SUP,"Mat argument must be of type seqdense");
 73:     MatGetSize(Q,&m,&n);
 74:     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);
 75:     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);
 76:   }
 77:   if (X->n!=Y->n) SETERRQ2(PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension X %D, Y %D",X->n,Y->n);

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

 86: /*@
 87:    BVMultVec - Computes y = beta*y + alpha*X*q.

 89:    Logically Collective on X

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

 97:    Output Parameter:
 98: .  y          - the modified vector

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

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

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

112:    Level: intermediate

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

128:   BVCheckSizes(X,1);
129:   BVCheckOp(X,1,multvec);

133:   VecGetSize(y,&N);
134:   VecGetLocalSize(y,&n);
135:   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);

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

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

147:    Logically Collective on X

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

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

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

167:    Level: advanced

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

183:   BVCheckSizes(X,1);

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

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

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

204:    Logically Collective on V

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

211:    Input/Output Parameter:
212: .  V - basis vectors

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

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

223:    Level: intermediate

225: .seealso: BVMult(), BVMultVec(), BVMultInPlaceTranspose(), BVSetActiveColumns()
226: @*/
227: PetscErrorCode BVMultInPlace(BV V,Mat Q,PetscInt s,PetscInt e)
228: {
230:   PetscBool      match;
231:   PetscInt       m,n;

239:   BVCheckSizes(V,1);
241:   PetscObjectTypeCompare((PetscObject)Q,MATSEQDENSE,&match);
242:   if (!match) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Mat argument must be of type seqdense");

244:   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);
245:   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);
246:   MatGetSize(Q,&m,&n);
247:   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);
248:   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);
249:   if (s>=e) return(0);

251:   PetscLogEventBegin(BV_MultInPlace,V,Q,0,0);
252:   (*V->ops->multinplace)(V,Q,s,e);
253:   PetscLogEventEnd(BV_MultInPlace,V,Q,0,0);
254:   PetscObjectStateIncrease((PetscObject)V);
255:   return(0);
256: }

258: /*@
259:    BVMultInPlaceTranspose - Update a set of vectors as V(:,s:e-1) = V*Q'(:,s:e-1).

261:    Logically Collective on V

263:    Input Parameters:
264: +  Q - a sequential dense matrix
265: .  s - first column of V to be overwritten
266: -  e - first column of V not to be overwritten

268:    Input/Output Parameter:
269: .  V - basis vectors

271:    Notes:
272:    This is a variant of BVMultInPlace() where the conjugate transpose
273:    of Q is used.

275:    Level: intermediate

277: .seealso: BVMultInPlace()
278: @*/
279: PetscErrorCode BVMultInPlaceTranspose(BV V,Mat Q,PetscInt s,PetscInt e)
280: {
282:   PetscBool      match;
283:   PetscInt       m,n;

291:   BVCheckSizes(V,1);
293:   PetscObjectTypeCompare((PetscObject)Q,MATSEQDENSE,&match);
294:   if (!match) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Mat argument must be of type seqdense");

296:   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);
297:   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);
298:   MatGetSize(Q,&m,&n);
299:   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);
300:   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);
301:   if (s>=e || !V->n) return(0);

303:   PetscLogEventBegin(BV_MultInPlace,V,Q,0,0);
304:   (*V->ops->multinplacetrans)(V,Q,s,e);
305:   PetscLogEventEnd(BV_MultInPlace,V,Q,0,0);
306:   PetscObjectStateIncrease((PetscObject)V);
307:   return(0);
308: }

310: /*@
311:    BVScale - Multiply the BV entries by a scalar value.

313:    Logically Collective on bv

315:    Input Parameters:
316: +  bv    - basis vectors
317: -  alpha - scaling factor

319:    Note:
320:    All active columns (except the leading ones) are scaled.

322:    Level: intermediate

324: .seealso: BVScaleColumn(), BVSetActiveColumns()
325: @*/
326: PetscErrorCode BVScale(BV bv,PetscScalar alpha)
327: {

334:   BVCheckSizes(bv,1);
335:   if (alpha == (PetscScalar)1.0) return(0);

337:   PetscLogEventBegin(BV_Scale,bv,0,0,0);
338:   if (bv->n) {
339:     (*bv->ops->scale)(bv,-1,alpha);
340:   }
341:   PetscLogEventEnd(BV_Scale,bv,0,0,0);
342:   PetscObjectStateIncrease((PetscObject)bv);
343:   return(0);
344: }

346: /*@
347:    BVScaleColumn - Scale one column of a BV.

349:    Logically Collective on bv

351:    Input Parameters:
352: +  bv    - basis vectors
353: .  j     - column number to be scaled
354: -  alpha - scaling factor

356:    Level: intermediate

358: .seealso: BVScale(), BVSetActiveColumns()
359: @*/
360: PetscErrorCode BVScaleColumn(BV bv,PetscInt j,PetscScalar alpha)
361: {

369:   BVCheckSizes(bv,1);

371:   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);
372:   if (alpha == (PetscScalar)1.0) return(0);

374:   PetscLogEventBegin(BV_Scale,bv,0,0,0);
375:   if (bv->n) {
376:     (*bv->ops->scale)(bv,j,alpha);
377:   }
378:   PetscLogEventEnd(BV_Scale,bv,0,0,0);
379:   PetscObjectStateIncrease((PetscObject)bv);
380:   return(0);
381: }

383: PETSC_STATIC_INLINE PetscErrorCode BVSetRandomColumn_Private(BV bv,PetscInt k)
384: {
386:   PetscInt       i,low,high;
387:   PetscScalar    *px,t;
388:   Vec            x;

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

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

410:    Logically Collective on bv

412:    Input Parameters:
413: .  bv - basis vectors

415:    Note:
416:    All active columns (except the leading ones) are modified.

418:    Level: advanced

420: .seealso: BVSetRandomContext(), BVSetRandomColumn(), BVSetRandomCond(), BVSetActiveColumns()
421: @*/
422: PetscErrorCode BVSetRandom(BV bv)
423: {
425:   PetscInt       k;

430:   BVCheckSizes(bv,1);

432:   BVGetRandomContext(bv,&bv->rand);
433:   PetscLogEventBegin(BV_SetRandom,bv,0,0,0);
434:   for (k=bv->l;k<bv->k;k++) {
435:     BVSetRandomColumn_Private(bv,k);
436:   }
437:   PetscLogEventEnd(BV_SetRandom,bv,0,0,0);
438:   PetscObjectStateIncrease((PetscObject)bv);
439:   return(0);
440: }

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

445:    Logically Collective on bv

447:    Input Parameters:
448: +  bv - basis vectors
449: -  j  - column number to be set

451:    Level: advanced

453: .seealso: BVSetRandomContext(), BVSetRandom(), BVSetActiveColumns()
454: @*/
455: PetscErrorCode BVSetRandomColumn(BV bv,PetscInt j)
456: {

463:   BVCheckSizes(bv,1);
464:   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);

466:   BVGetRandomContext(bv,&bv->rand);
467:   PetscLogEventBegin(BV_SetRandom,bv,0,0,0);
468:   BVSetRandomColumn_Private(bv,j);
469:   PetscLogEventEnd(BV_SetRandom,bv,0,0,0);
470:   PetscObjectStateIncrease((PetscObject)bv);
471:   return(0);
472: }

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

478:    Logically Collective on bv

480:    Input Parameters:
481: +  bv    - basis vectors
482: -  condn - condition number

484:    Note:
485:    All active columns (except the leading ones) are modified.

487:    Level: advanced

489: .seealso: BVSetRandomContext(), BVSetRandomColumn(), BVSetActiveColumns()
490: @*/
491: PetscErrorCode BVSetRandomCond(BV bv,PetscReal condn)
492: {
494:   PetscInt       k,i;
495:   PetscScalar    *eig,*d;
496:   DS             ds;
497:   Mat            A,X,Xt,M,G;

502:   BVCheckSizes(bv,1);

504:   BVGetRandomContext(bv,&bv->rand);
505:   PetscLogEventBegin(BV_SetRandom,bv,0,0,0);
506:   /* B = rand(n,k) */
507:   for (k=bv->l;k<bv->k;k++) {
508:     BVSetRandomColumn_Private(bv,k);
509:   }
510:   DSCreate(PetscObjectComm((PetscObject)bv),&ds);
511:   DSSetType(ds,DSHEP);
512:   DSAllocate(ds,bv->m);
513:   DSSetDimensions(ds,bv->k,0,bv->l,bv->k);
514:   /* [V,S] = eig(B'*B) */
515:   DSGetMat(ds,DS_MAT_A,&A);
516:   BVDot(bv,bv,A);
517:   DSRestoreMat(ds,DS_MAT_A,&A);
518:   PetscMalloc1(bv->m,&eig);
519:   DSSolve(ds,eig,NULL);
520:   DSSynchronize(ds,eig,NULL);
521:   DSVectors(ds,DS_MAT_X,NULL,NULL);
522:   /* M = diag(linspace(1/condn,1,n)./sqrt(diag(S)))' */
523:   MatCreateSeqDense(PETSC_COMM_SELF,bv->k,bv->k,NULL,&M);
524:   MatZeroEntries(M);
525:   MatDenseGetArray(M,&d);
526:   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]);
527:   MatDenseRestoreArray(M,&d);
528:   /* G = X*M*X' */
529:   MatCreateSeqDense(PETSC_COMM_SELF,bv->k,bv->k,NULL,&Xt);
530:   DSGetMat(ds,DS_MAT_X,&X);
531:   MatTranspose(X,MAT_REUSE_MATRIX,&Xt);
532:   MatProductCreate(Xt,M,NULL,&G);
533:   MatProductSetType(G,MATPRODUCT_PtAP);
534:   MatProductSetFromOptions(G);
535:   MatProductSymbolic(G);
536:   MatProductNumeric(G);
537:   MatProductClear(G);
538:   MatDestroy(&X);
539:   MatDestroy(&Xt);
540:   MatDestroy(&M);
541:   /* B = B*G */
542:   BVMultInPlace(bv,G,bv->l,bv->k);
543:   MatDestroy(&G);
544:   PetscFree(eig);
545:   DSDestroy(&ds);
546:   PetscLogEventEnd(BV_SetRandom,bv,0,0,0);
547:   PetscObjectStateIncrease((PetscObject)bv);
548:   return(0);
549: }

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

554:    Neighbor-wise Collective on A

556:    Input Parameters:
557: +  V - basis vectors context
558: -  A - the matrix

560:    Output Parameter:
561: .  Y - the result

563:    Notes:
564:    Both V and Y must be distributed in the same manner. Only active columns
565:    (excluding the leading ones) are processed.
566:    In the result Y, columns are overwritten starting from the leading ones.
567:    The number of active columns in V and Y should match, although they need
568:    not be the same columns.

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

573:    Level: beginner

575: .seealso: BVCopy(), BVSetActiveColumns(), BVMatMultColumn(), BVMatMultTranspose(), BVMatMultHermitianTranspose(), BVSetMatMultMethod()
576: @*/
577: PetscErrorCode BVMatMult(BV V,Mat A,BV Y)
578: {

584:   BVCheckSizes(V,1);
585:   BVCheckOp(V,1,matmult);
590:   BVCheckSizes(Y,3);
593:   if (V->n!=Y->n) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension V %D, Y %D",V->n,Y->n);
594:   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);

596:   PetscLogEventBegin(BV_MatMult,V,A,Y,0);
597:   (*V->ops->matmult)(V,A,Y);
598:   PetscLogEventEnd(BV_MatMult,V,A,Y,0);
599:   PetscObjectStateIncrease((PetscObject)Y);
600:   return(0);
601: }

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

607:    Neighbor-wise Collective on A

609:    Input Parameters:
610: +  V - basis vectors context
611: -  A - the matrix

613:    Output Parameter:
614: .  Y - the result

616:    Notes:
617:    Both V and Y must be distributed in the same manner. Only active columns
618:    (excluding the leading ones) are processed.
619:    In the result Y, columns are overwritten starting from the leading ones.
620:    The number of active columns in V and Y should match, although they need
621:    not be the same columns.

623:    Currently implemented via MatCreateTranspose().

625:    Level: beginner

627: .seealso: BVMatMult(), BVMatMultHermitianTranspose()
628: @*/
629: PetscErrorCode BVMatMultTranspose(BV V,Mat A,BV Y)
630: {
632:   Mat            AT;

637:   BVCheckSizes(V,1);
642:   BVCheckSizes(Y,3);
645:   if (V->n!=Y->n) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension V %D, Y %D",V->n,Y->n);
646:   if (V->k-V->l>Y->m-Y->l) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"Y has %D non-leading columns, not enough to store %D columns",Y->m-Y->l,V->k-V->l);

648:   MatCreateTranspose(A,&AT);
649:   BVMatMult(V,AT,Y);
650:   MatDestroy(&AT);
651:   return(0);
652: }

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

658:    Neighbor-wise Collective on A

660:    Input Parameters:
661: +  V - basis vectors context
662: -  A - the matrix

664:    Output Parameter:
665: .  Y - the result

667:    Note:
668:    Both V and Y must be distributed in the same manner. Only active columns
669:    (excluding the leading ones) are processed.
670:    In the result Y, columns are overwritten starting from the leading ones.
671:    The number of active columns in V and Y should match, although they need
672:    not be the same columns.

674:    Currently implemented via MatCreateHermitianTranspose().

676:    Level: beginner

678: .seealso: BVMatMult(), BVMatMultTranspose()
679: @*/
680: PetscErrorCode BVMatMultHermitianTranspose(BV V,Mat A,BV Y)
681: {
683:   Mat            AH;

688:   BVCheckSizes(V,1);
693:   BVCheckSizes(Y,3);
696:   if (V->n!=Y->n) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension V %D, Y %D",V->n,Y->n);
697:   if (V->k-V->l>Y->m-Y->l) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"Y has %D non-leading columns, not enough to store %D columns",Y->m-Y->l,V->k-V->l);

699:   MatCreateHermitianTranspose(A,&AH);
700:   BVMatMult(V,AH,Y);
701:   MatDestroy(&AH);
702:   return(0);
703: }

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

709:    Neighbor-wise Collective on A

711:    Input Parameters:
712: +  V - basis vectors context
713: .  A - the matrix
714: -  j - the column

716:    Level: beginner

718: .seealso: BVMatMult(), BVMatMultTransposeColumn(), BVMatMultHermitianTransposeColumn()
719: @*/
720: PetscErrorCode BVMatMultColumn(BV V,Mat A,PetscInt j)
721: {
723:   Vec            vj,vj1;

728:   BVCheckSizes(V,1);
732:   if (j<0) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
733:   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);

735:   PetscLogEventBegin(BV_MatMultVec,V,A,0,0);
736:   BVGetColumn(V,j,&vj);
737:   BVGetColumn(V,j+1,&vj1);
738:   MatMult(A,vj,vj1);
739:   BVRestoreColumn(V,j,&vj);
740:   BVRestoreColumn(V,j+1,&vj1);
741:   PetscLogEventEnd(BV_MatMultVec,V,A,0,0);
742:   PetscObjectStateIncrease((PetscObject)V);
743:   return(0);
744: }

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

750:    Neighbor-wise Collective on A

752:    Input Parameters:
753: +  V - basis vectors context
754: .  A - the matrix
755: -  j - the column

757:    Level: beginner

759: .seealso: BVMatMult(), BVMatMultColumn(), BVMatMultHermitianTransposeColumn()
760: @*/
761: PetscErrorCode BVMatMultTransposeColumn(BV V,Mat A,PetscInt j)
762: {
764:   Vec            vj,vj1;

769:   BVCheckSizes(V,1);
773:   if (j<0) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
774:   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);

776:   PetscLogEventBegin(BV_MatMultVec,V,A,0,0);
777:   BVGetColumn(V,j,&vj);
778:   BVGetColumn(V,j+1,&vj1);
779:   MatMultTranspose(A,vj,vj1);
780:   BVRestoreColumn(V,j,&vj);
781:   BVRestoreColumn(V,j+1,&vj1);
782:   PetscLogEventEnd(BV_MatMultVec,V,A,0,0);
783:   PetscObjectStateIncrease((PetscObject)V);
784:   return(0);
785: }

787: /*@
788:    BVMatMultHermitianTransposeColumn - Computes the matrix-vector product for a specified
789:    column, storing the result in the next column: v_{j+1}=A*v_j.

791:    Neighbor-wise Collective on A

793:    Input Parameters:
794: +  V - basis vectors context
795: .  A - the matrix
796: -  j - the column

798:    Level: beginner

800: .seealso: BVMatMult(), BVMatMultColumn(), BVMatMultTransposeColumn()
801: @*/
802: PetscErrorCode BVMatMultHermitianTransposeColumn(BV V,Mat A,PetscInt j)
803: {
805:   Vec            vj,vj1;

810:   BVCheckSizes(V,1);
814:   if (j<0) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
815:   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);

817:   PetscLogEventBegin(BV_MatMultVec,V,A,0,0);
818:   BVGetColumn(V,j,&vj);
819:   BVGetColumn(V,j+1,&vj1);
820:   MatMultHermitianTranspose(A,vj,vj1);
821:   BVRestoreColumn(V,j,&vj);
822:   BVRestoreColumn(V,j+1,&vj1);
823:   PetscLogEventEnd(BV_MatMultVec,V,A,0,0);
824:   PetscObjectStateIncrease((PetscObject)V);
825:   return(0);
826: }