Actual source code: bvops.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2019, 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>      /*I "slepcbv.h" I*/
 15: #include <slepcds.h>

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

 20:    Logically Collective on BV

 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);
 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) {
 70:     PetscObjectTypeCompare((PetscObject)Q,MATSEQDENSE,&match);
 71:     if (!match) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_SUP,"Mat argument must be of type seqdense");
 72:     MatGetSize(Q,&m,&n);
 73:     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);
 74:     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);
 75:   }
 76:   if (X->n!=Y->n) SETERRQ2(PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension X %D, Y %D",X->n,Y->n);

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

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

 88:    Logically Collective on BV and Vec

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

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

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

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

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

111:    Level: intermediate

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

127:   BVCheckSizes(X,1);

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 BV

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 BV

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(), BVMultInPlaceTranspose(), BVSetActiveColumns()
224: @*/
225: PetscErrorCode BVMultInPlace(BV V,Mat Q,PetscInt s,PetscInt e)
226: {
228:   PetscBool      match;
229:   PetscInt       m,n;

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

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

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

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

259:    Logically Collective on BV

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

266:    Input/Output Parameter:
267: .  V - basis vectors

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

273:    Level: intermediate

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

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

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

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

308: /*@
309:    BVScale - Multiply the BV entries by a scalar value.

311:    Logically Collective on BV

313:    Input Parameters:
314: +  bv    - basis vectors
315: -  alpha - scaling factor

317:    Note:
318:    All active columns (except the leading ones) are scaled.

320:    Level: intermediate

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

332:   BVCheckSizes(bv,1);
333:   if (alpha == (PetscScalar)1.0) return(0);

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

344: /*@
345:    BVScaleColumn - Scale one column of a BV.

347:    Logically Collective on BV

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

354:    Level: intermediate

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

367:   BVCheckSizes(bv,1);

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

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

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

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

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

408:    Logically Collective on BV

410:    Input Parameters:
411: .  bv - basis vectors

413:    Note:
414:    All active columns (except the leading ones) are modified.

416:    Level: advanced

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

428:   BVCheckSizes(bv,1);

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

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

443:    Logically Collective on BV

445:    Input Parameters:
446: +  bv - basis vectors
447: -  j  - column number to be set

449:    Level: advanced

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

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

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

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

476:    Logically Collective on BV

478:    Input Parameters:
479: +  bv    - basis vectors
480: -  condn - condition number

482:    Note:
483:    All active columns (except the leading ones) are modified.

485:    Level: advanced

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

500:   BVCheckSizes(bv,1);

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

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

549:    Neighbor-wise Collective on Mat and BV

551:    Input Parameters:
552: +  V - basis vectors context
553: -  A - the matrix

555:    Output Parameter:
556: .  Y - the result

558:    Note:
559:    Both V and Y must be distributed in the same manner. Only active columns
560:    (excluding the leading ones) are processed.
561:    In the result Y, columns are overwritten starting from the leading ones.
562:    The number of active columns in V and Y should match, although they need
563:    not be the same columns.

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

568:    Level: beginner

570: .seealso: BVCopy(), BVSetActiveColumns(), BVMatMultColumn(), BVMatMultHermitianTranspose(), BVSetMatMultMethod()
571: @*/
572: PetscErrorCode BVMatMult(BV V,Mat A,BV Y)
573: {

579:   BVCheckSizes(V,1);
584:   BVCheckSizes(Y,3);
587:   if (V->n!=Y->n) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension V %D, Y %D",V->n,Y->n);
588:   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);

590:   PetscLogEventBegin(BV_MatMult,V,A,Y,0);
591:   (*V->ops->matmult)(V,A,Y);
592:   PetscLogEventEnd(BV_MatMult,V,A,Y,0);
593:   PetscObjectStateIncrease((PetscObject)Y);
594:   return(0);
595: }

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

601:    Neighbor-wise Collective on Mat and BV

603:    Input Parameters:
604: +  V - basis vectors context
605: -  A - the matrix

607:    Output Parameter:
608: .  Y - the result

610:    Note:
611:    Both V and Y must be distributed in the same manner. Only active columns
612:    (excluding the leading ones) are processed.
613:    In the result Y, columns are overwritten starting from the leading ones.

615:    As opposed to BVMatMult(), this operation is always done column by column,
616:    with a sequence of calls to MatMultHermitianTranspose().

618:    Level: beginner

620: .seealso: BVCopy(), BVSetActiveColumns(), BVMatMult(), BVMatMultColumn()
621: @*/
622: PetscErrorCode BVMatMultHermitianTranspose(BV V,Mat A,BV Y)
623: {
625:   PetscInt       j;
626:   Vec            z,f;

631:   BVCheckSizes(V,1);
636:   BVCheckSizes(Y,3);
639:   if (V->n!=Y->n) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension V %D, Y %D",V->n,Y->n);
640:   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);

642:   PetscLogEventBegin(BV_MatMult,V,A,Y,0);
643:   for (j=0;j<V->k-V->l;j++) {
644:     BVGetColumn(V,V->l+j,&z);
645:     BVGetColumn(Y,Y->l+j,&f);
646:     MatMultHermitianTranspose(A,z,f);
647:     BVRestoreColumn(V,V->l+j,&z);
648:     BVRestoreColumn(Y,Y->l+j,&f);
649:   }
650:   PetscLogEventEnd(BV_MatMult,V,A,Y,0);
651:   PetscObjectStateIncrease((PetscObject)Y);
652:   return(0);
653: }

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

659:    Neighbor-wise Collective on Mat and BV

661:    Input Parameters:
662: +  V - basis vectors context
663: .  A - the matrix
664: -  j - the column

666:    Output Parameter:
667: .  Y - the result

669:    Level: beginner

671: .seealso: BVMatMult()
672: @*/
673: PetscErrorCode BVMatMultColumn(BV V,Mat A,PetscInt j)
674: {
676:   Vec            vj,vj1;

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

687:   PetscLogEventBegin(BV_MatMultVec,V,A,0,0);
688:   BVGetColumn(V,j,&vj);
689:   BVGetColumn(V,j+1,&vj1);
690:   MatMult(A,vj,vj1);
691:   BVRestoreColumn(V,j,&vj);
692:   BVRestoreColumn(V,j+1,&vj1);
693:   PetscLogEventEnd(BV_MatMultVec,V,A,0,0);
694:   PetscObjectStateIncrease((PetscObject)V);
695:   return(0);
696: }