Actual source code: bvglobal.c

slepc-main 2024-11-15
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 operations involving global communication
 12: */

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

 16: /*
 17:   BVDot for the particular case of non-standard inner product with
 18:   matrix B, which is assumed to be symmetric (or complex Hermitian)
 19: */
 20: static inline PetscErrorCode BVDot_Private(BV X,BV Y,Mat M)
 21: {
 22:   PetscObjectId  idx,idy;
 23:   PetscInt       i,j,jend,m;
 24:   PetscScalar    *marray;
 25:   PetscBool      symm=PETSC_FALSE;
 26:   Mat            B;
 27:   Vec            z;

 29:   PetscFunctionBegin;
 30:   BVCheckOp(Y,1,dotvec);
 31:   PetscCall(MatGetSize(M,&m,NULL));
 32:   PetscCall(MatDenseGetArray(M,&marray));
 33:   PetscCall(PetscObjectGetId((PetscObject)X,&idx));
 34:   PetscCall(PetscObjectGetId((PetscObject)Y,&idy));
 35:   B = Y->matrix;
 36:   Y->matrix = NULL;
 37:   if (idx==idy) symm=PETSC_TRUE;  /* M=X'BX is symmetric */
 38:   jend = X->k;
 39:   for (j=X->l;j<jend;j++) {
 40:     if (symm) Y->k = j+1;
 41:     PetscCall(BVGetColumn(X->cached,j,&z));
 42:     PetscUseTypeMethod(Y,dotvec,z,marray+j*m+Y->l);
 43:     PetscCall(BVRestoreColumn(X->cached,j,&z));
 44:     if (symm) {
 45:       for (i=X->l;i<j;i++)
 46:         marray[j+i*m] = PetscConj(marray[i+j*m]);
 47:     }
 48:   }
 49:   PetscCall(MatDenseRestoreArray(M,&marray));
 50:   Y->matrix = B;
 51:   PetscFunctionReturn(PETSC_SUCCESS);
 52: }

 54: /*@
 55:    BVDot - Computes the 'block-dot' product of two basis vectors objects.

 57:    Collective

 59:    Input Parameters:
 60: +  X - first basis vectors
 61: -  Y - second basis vectors

 63:    Output Parameter:
 64: .  M - the resulting matrix

 66:    Notes:
 67:    This is the generalization of VecDot() for a collection of vectors, M = Y^H*X.
 68:    The result is a matrix M whose entry m_ij is equal to y_i^H x_j (where y_i^H
 69:    denotes the conjugate transpose of y_i).

 71:    If a non-standard inner product has been specified with BVSetMatrix(),
 72:    then the result is M = Y^H*B*X. In this case, both X and Y must have
 73:    the same associated matrix.

 75:    On entry, M must be a sequential dense Mat with dimensions m,n at least, where
 76:    m is the number of active columns of Y and n is the number of active columns of X.
 77:    Only rows (resp. columns) of M starting from ly (resp. lx) are computed,
 78:    where ly (resp. lx) is the number of leading columns of Y (resp. X).

 80:    X and Y need not be different objects.

 82:    Level: intermediate

 84: .seealso: BVDotVec(), BVDotColumn(), BVSetActiveColumns(), BVSetMatrix()
 85: @*/
 86: PetscErrorCode BVDot(BV X,BV Y,Mat M)
 87: {
 88:   PetscInt       m,n;

 90:   PetscFunctionBegin;
 95:   BVCheckSizes(X,1);
 97:   BVCheckSizes(Y,2);
 99:   PetscCheckSameTypeAndComm(X,1,Y,2);
100:   SlepcMatCheckSeq(M);

102:   PetscCall(MatGetSize(M,&m,&n));
103:   PetscCheck(m>=Y->k,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_SIZ,"Mat argument has %" PetscInt_FMT " rows, should have at least %" PetscInt_FMT,m,Y->k);
104:   PetscCheck(n>=X->k,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_SIZ,"Mat argument has %" PetscInt_FMT " columns, should have at least %" PetscInt_FMT,n,X->k);
105:   PetscCheck(X->n==Y->n,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension X %" PetscInt_FMT ", Y %" PetscInt_FMT,X->n,Y->n);
106:   PetscCheck(X->matrix==Y->matrix,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_WRONGSTATE,"X and Y must have the same inner product matrix");
107:   if (X->l==X->k || Y->l==Y->k) PetscFunctionReturn(PETSC_SUCCESS);

109:   PetscCall(PetscLogEventBegin(BV_Dot,X,Y,0,0));
110:   if (X->matrix) { /* non-standard inner product */
111:     /* compute BX first */
112:     PetscCall(BV_IPMatMultBV(X));
113:     if (X->vmm==BV_MATMULT_VECS) {
114:       /* perform computation column by column */
115:       PetscCall(BVDot_Private(X,Y,M));
116:     } else PetscUseTypeMethod(X->cached,dot,Y,M);
117:   } else PetscUseTypeMethod(X,dot,Y,M);
118:   PetscCall(PetscLogEventEnd(BV_Dot,X,Y,0,0));
119:   PetscFunctionReturn(PETSC_SUCCESS);
120: }

122: /*@
123:    BVDotVec - Computes multiple dot products of a vector against all the
124:    column vectors of a BV.

126:    Collective

128:    Input Parameters:
129: +  X - basis vectors
130: -  y - a vector

132:    Output Parameter:
133: .  m - an array where the result must be placed

135:    Notes:
136:    This is analogue to VecMDot(), but using BV to represent a collection
137:    of vectors. The result is m = X^H*y, so m_i is equal to x_i^H y. Note
138:    that here X is transposed as opposed to BVDot().

140:    If a non-standard inner product has been specified with BVSetMatrix(),
141:    then the result is m = X^H*B*y.

143:    The length of array m must be equal to the number of active columns of X
144:    minus the number of leading columns, i.e. the first entry of m is the
145:    product of the first non-leading column with y.

147:    Level: intermediate

149: .seealso: BVDot(), BVDotColumn(), BVSetActiveColumns(), BVSetMatrix()
150: @*/
151: PetscErrorCode BVDotVec(BV X,Vec y,PetscScalar m[])
152: {
153:   PetscInt       n;

155:   PetscFunctionBegin;
159:   BVCheckSizes(X,1);
160:   BVCheckOp(X,1,dotvec);
162:   PetscCheckSameTypeAndComm(X,1,y,2);

164:   PetscCall(VecGetLocalSize(y,&n));
165:   PetscCheck(X->n==n,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension X %" PetscInt_FMT ", y %" PetscInt_FMT,X->n,n);

167:   PetscCall(PetscLogEventBegin(BV_DotVec,X,y,0,0));
168:   PetscUseTypeMethod(X,dotvec,y,m);
169:   PetscCall(PetscLogEventEnd(BV_DotVec,X,y,0,0));
170:   PetscFunctionReturn(PETSC_SUCCESS);
171: }

173: /*@
174:    BVDotVecBegin - Starts a split phase dot product computation.

176:    Input Parameters:
177: +  X - basis vectors
178: .  y - a vector
179: -  m - an array where the result will go (can be NULL)

181:    Note:
182:    Each call to BVDotVecBegin() should be paired with a call to BVDotVecEnd().

184:    Level: advanced

186: .seealso: BVDotVecEnd(), BVDotVec()
187: @*/
188: PetscErrorCode BVDotVecBegin(BV X,Vec y,PetscScalar *m)
189: {
190:   PetscInt            i,n,nv;
191:   PetscSplitReduction *sr;
192:   MPI_Comm            comm;

194:   PetscFunctionBegin;
198:   BVCheckSizes(X,1);
200:   PetscCheckSameTypeAndComm(X,1,y,2);

202:   PetscCall(VecGetLocalSize(y,&n));
203:   PetscCheck(X->n==n,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension X %" PetscInt_FMT ", y %" PetscInt_FMT,X->n,n);

205:   if (X->ops->dotvec_begin) PetscUseTypeMethod(X,dotvec_begin,y,m);
206:   else {
207:     BVCheckOp(X,1,dotvec_local);
208:     nv = X->k-X->l;
209:     PetscCall(PetscObjectGetComm((PetscObject)X,&comm));
210:     PetscCall(PetscSplitReductionGet(comm,&sr));
211:     PetscCheck(sr->state==STATE_BEGIN,PetscObjectComm((PetscObject)X),PETSC_ERR_ORDER,"Called before all VecxxxEnd() called");
212:     for (i=0;i<nv;i++) {
213:       if (sr->numopsbegin+i >= sr->maxops) PetscCall(PetscSplitReductionExtend(sr));
214:       sr->reducetype[sr->numopsbegin+i] = PETSC_SR_REDUCE_SUM;
215:       sr->invecs[sr->numopsbegin+i]     = (void*)X;
216:     }
217:     PetscCall(PetscLogEventBegin(BV_DotVec,X,y,0,0));
218:     PetscUseTypeMethod(X,dotvec_local,y,sr->lvalues+sr->numopsbegin);
219:     sr->numopsbegin += nv;
220:     PetscCall(PetscLogEventEnd(BV_DotVec,X,y,0,0));
221:   }
222:   PetscFunctionReturn(PETSC_SUCCESS);
223: }

225: /*@
226:    BVDotVecEnd - Ends a split phase dot product computation.

228:    Input Parameters:
229: +  X - basis vectors
230: .  y - a vector
231: -  m - an array where the result will go

233:    Note:
234:    Each call to BVDotVecBegin() should be paired with a call to BVDotVecEnd().

236:    Level: advanced

238: .seealso: BVDotVecBegin(), BVDotVec()
239: @*/
240: PetscErrorCode BVDotVecEnd(BV X,Vec y,PetscScalar *m)
241: {
242:   PetscInt            i,nv;
243:   PetscSplitReduction *sr;
244:   MPI_Comm            comm;

246:   PetscFunctionBegin;
249:   BVCheckSizes(X,1);

251:   if (X->ops->dotvec_end) PetscUseTypeMethod(X,dotvec_end,y,m);
252:   else {
253:     nv = X->k-X->l;
254:     PetscCall(PetscObjectGetComm((PetscObject)X,&comm));
255:     PetscCall(PetscSplitReductionGet(comm,&sr));
256:     PetscCall(PetscSplitReductionEnd(sr));

258:     PetscCheck(sr->numopsend<sr->numopsbegin,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_WRONGSTATE,"Called VecxxxEnd() more times than VecxxxBegin()");
259:     PetscCheck((void*)X==sr->invecs[sr->numopsend],PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_WRONGSTATE,"Called BVxxxEnd() in a different order or with a different BV than BVxxxBegin()");
260:     PetscCheck(sr->reducetype[sr->numopsend]==PETSC_SR_REDUCE_SUM,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_WRONGSTATE,"Wrong type of reduction");
261:     for (i=0;i<nv;i++) m[i] = sr->gvalues[sr->numopsend++];

263:     /* Finished getting all the results so reset to no outstanding requests */
264:     if (sr->numopsend == sr->numopsbegin) {
265:       sr->state       = STATE_BEGIN;
266:       sr->numopsend   = 0;
267:       sr->numopsbegin = 0;
268:     }
269:   }
270:   PetscFunctionReturn(PETSC_SUCCESS);
271: }

273: /*@
274:    BVDotColumn - Computes multiple dot products of a column against all the
275:    previous columns of a BV.

277:    Collective

279:    Input Parameters:
280: +  X - basis vectors
281: -  j - the column index

283:    Output Parameter:
284: .  q - an array where the result must be placed

286:    Notes:
287:    This operation is equivalent to BVDotVec() but it uses column j of X
288:    rather than taking a Vec as an argument. The number of active columns of
289:    X is set to j before the computation, and restored afterwards.
290:    If X has leading columns specified, then these columns do not participate
291:    in the computation. Therefore, the length of array q must be equal to j
292:    minus the number of leading columns.

294:    Developer Notes:
295:    If q is NULL, then the result is written in position nc+l of the internal
296:    buffer vector, see BVGetBufferVec().

298:    Level: advanced

300: .seealso: BVDot(), BVDotVec(), BVSetActiveColumns(), BVSetMatrix()
301: @*/
302: PetscErrorCode BVDotColumn(BV X,PetscInt j,PetscScalar *q)
303: {
304:   PetscInt       ksave;
305:   Vec            y;

307:   PetscFunctionBegin;
311:   BVCheckSizes(X,1);
312:   BVCheckOp(X,1,dotvec);

314:   PetscCheck(j>=0,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
315:   PetscCheck(j<X->m,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_OUTOFRANGE,"Index j=%" PetscInt_FMT " but BV only has %" PetscInt_FMT " columns",j,X->m);

317:   PetscCall(PetscLogEventBegin(BV_DotVec,X,0,0,0));
318:   ksave = X->k;
319:   X->k = j;
320:   if (!q && !X->buffer) PetscCall(BVGetBufferVec(X,&X->buffer));
321:   PetscCall(BVGetColumn(X,j,&y));
322:   PetscUseTypeMethod(X,dotvec,y,q);
323:   PetscCall(BVRestoreColumn(X,j,&y));
324:   X->k = ksave;
325:   PetscCall(PetscLogEventEnd(BV_DotVec,X,0,0,0));
326:   PetscFunctionReturn(PETSC_SUCCESS);
327: }

329: /*@
330:    BVDotColumnBegin - Starts a split phase dot product computation.

332:    Input Parameters:
333: +  X - basis vectors
334: .  j - the column index
335: -  m - an array where the result will go (can be NULL)

337:    Note:
338:    Each call to BVDotColumnBegin() should be paired with a call to BVDotColumnEnd().

340:    Level: advanced

342: .seealso: BVDotColumnEnd(), BVDotColumn()
343: @*/
344: PetscErrorCode BVDotColumnBegin(BV X,PetscInt j,PetscScalar *m)
345: {
346:   PetscInt            i,nv,ksave;
347:   PetscSplitReduction *sr;
348:   MPI_Comm            comm;
349:   Vec                 y;

351:   PetscFunctionBegin;
355:   BVCheckSizes(X,1);

357:   PetscCheck(j>=0,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
358:   PetscCheck(j<X->m,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_OUTOFRANGE,"Index j=%" PetscInt_FMT " but BV only has %" PetscInt_FMT " columns",j,X->m);
359:   ksave = X->k;
360:   X->k = j;
361:   PetscCall(BVGetColumn(X,j,&y));

363:   if (X->ops->dotvec_begin) PetscUseTypeMethod(X,dotvec_begin,y,m);
364:   else {
365:     BVCheckOp(X,1,dotvec_local);
366:     nv = X->k-X->l;
367:     PetscCall(PetscObjectGetComm((PetscObject)X,&comm));
368:     PetscCall(PetscSplitReductionGet(comm,&sr));
369:     PetscCheck(sr->state==STATE_BEGIN,PetscObjectComm((PetscObject)X),PETSC_ERR_ORDER,"Called before all VecxxxEnd() called");
370:     for (i=0;i<nv;i++) {
371:       if (sr->numopsbegin+i >= sr->maxops) PetscCall(PetscSplitReductionExtend(sr));
372:       sr->reducetype[sr->numopsbegin+i] = PETSC_SR_REDUCE_SUM;
373:       sr->invecs[sr->numopsbegin+i]     = (void*)X;
374:     }
375:     PetscCall(PetscLogEventBegin(BV_DotVec,X,0,0,0));
376:     PetscUseTypeMethod(X,dotvec_local,y,sr->lvalues+sr->numopsbegin);
377:     sr->numopsbegin += nv;
378:     PetscCall(PetscLogEventEnd(BV_DotVec,X,0,0,0));
379:   }
380:   PetscCall(BVRestoreColumn(X,j,&y));
381:   X->k = ksave;
382:   PetscFunctionReturn(PETSC_SUCCESS);
383: }

385: /*@
386:    BVDotColumnEnd - Ends a split phase dot product computation.

388:    Input Parameters:
389: +  X - basis vectors
390: .  j - the column index
391: -  m - an array where the result will go

393:    Notes:
394:    Each call to BVDotColumnBegin() should be paired with a call to BVDotColumnEnd().

396:    Level: advanced

398: .seealso: BVDotColumnBegin(), BVDotColumn()
399: @*/
400: PetscErrorCode BVDotColumnEnd(BV X,PetscInt j,PetscScalar *m)
401: {
402:   PetscInt            i,nv,ksave;
403:   PetscSplitReduction *sr;
404:   MPI_Comm            comm;
405:   Vec                 y;

407:   PetscFunctionBegin;
411:   BVCheckSizes(X,1);

413:   PetscCheck(j>=0,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
414:   PetscCheck(j<X->m,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_OUTOFRANGE,"Index j=%" PetscInt_FMT " but BV only has %" PetscInt_FMT " columns",j,X->m);
415:   ksave = X->k;
416:   X->k = j;

418:   if (X->ops->dotvec_end) {
419:     PetscCall(BVGetColumn(X,j,&y));
420:     PetscUseTypeMethod(X,dotvec_end,y,m);
421:     PetscCall(BVRestoreColumn(X,j,&y));
422:   } else {
423:     nv = X->k-X->l;
424:     PetscCall(PetscObjectGetComm((PetscObject)X,&comm));
425:     PetscCall(PetscSplitReductionGet(comm,&sr));
426:     PetscCall(PetscSplitReductionEnd(sr));

428:     PetscCheck(sr->numopsend<sr->numopsbegin,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_WRONGSTATE,"Called VecxxxEnd() more times than VecxxxBegin()");
429:     PetscCheck((void*)X==sr->invecs[sr->numopsend],PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_WRONGSTATE,"Called BVxxxEnd() in a different order or with a different BV than BVxxxBegin()");
430:     PetscCheck(sr->reducetype[sr->numopsend]==PETSC_SR_REDUCE_SUM,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_WRONGSTATE,"Wrong type of reduction");
431:     for (i=0;i<nv;i++) m[i] = sr->gvalues[sr->numopsend++];

433:     /* Finished getting all the results so reset to no outstanding requests */
434:     if (sr->numopsend == sr->numopsbegin) {
435:       sr->state       = STATE_BEGIN;
436:       sr->numopsend   = 0;
437:       sr->numopsbegin = 0;
438:     }
439:   }
440:   X->k = ksave;
441:   PetscFunctionReturn(PETSC_SUCCESS);
442: }

444: static inline PetscErrorCode BVNorm_Private(BV bv,Vec z,NormType type,PetscReal *val)
445: {
446:   PetscScalar    p;

448:   PetscFunctionBegin;
449:   PetscCall(BV_IPMatMult(bv,z));
450:   PetscCall(VecDot(bv->Bx,z,&p));
451:   PetscCall(BV_SafeSqrt(bv,p,val));
452:   PetscFunctionReturn(PETSC_SUCCESS);
453: }

455: static inline PetscErrorCode BVNorm_Begin_Private(BV bv,Vec z,NormType type,PetscReal *val)
456: {
457:   PetscScalar    p;

459:   PetscFunctionBegin;
460:   PetscCall(BV_IPMatMult(bv,z));
461:   PetscCall(VecDotBegin(bv->Bx,z,&p));
462:   PetscFunctionReturn(PETSC_SUCCESS);
463: }

465: static inline PetscErrorCode BVNorm_End_Private(BV bv,Vec z,NormType type,PetscReal *val)
466: {
467:   PetscScalar    p;

469:   PetscFunctionBegin;
470:   PetscCall(VecDotEnd(bv->Bx,z,&p));
471:   PetscCall(BV_SafeSqrt(bv,p,val));
472:   PetscFunctionReturn(PETSC_SUCCESS);
473: }

475: /*@
476:    BVNorm - Computes the matrix norm of the BV.

478:    Collective

480:    Input Parameters:
481: +  bv   - basis vectors
482: -  type - the norm type

484:    Output Parameter:
485: .  val  - the norm

487:    Notes:
488:    All active columns (except the leading ones) are considered as a matrix.
489:    The allowed norms are NORM_1, NORM_FROBENIUS, and NORM_INFINITY.

491:    This operation fails if a non-standard inner product has been
492:    specified with BVSetMatrix().

494:    Level: intermediate

496: .seealso: BVNormVec(), BVNormColumn(), BVNormalize(), BVSetActiveColumns(), BVSetMatrix()
497: @*/
498: PetscErrorCode BVNorm(BV bv,NormType type,PetscReal *val)
499: {
500:   PetscFunctionBegin;
503:   PetscAssertPointer(val,3);
505:   BVCheckSizes(bv,1);

507:   PetscCheck(type!=NORM_2 && type!=NORM_1_AND_2,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Requested norm not available");
508:   PetscCheck(!bv->matrix,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Matrix norm not available for non-standard inner product");

510:   PetscCall(PetscLogEventBegin(BV_Norm,bv,0,0,0));
511:   PetscUseTypeMethod(bv,norm,-1,type,val);
512:   PetscCall(PetscLogEventEnd(BV_Norm,bv,0,0,0));
513:   PetscFunctionReturn(PETSC_SUCCESS);
514: }

516: /*@
517:    BVNormVec - Computes the norm of a given vector.

519:    Collective

521:    Input Parameters:
522: +  bv   - basis vectors
523: .  v    - the vector
524: -  type - the norm type

526:    Output Parameter:
527: .  val  - the norm

529:    Notes:
530:    This is the analogue of BVNormColumn() but for a vector that is not in the BV.
531:    If a non-standard inner product has been specified with BVSetMatrix(),
532:    then the returned value is sqrt(v'*B*v), where B is the inner product
533:    matrix (argument 'type' is ignored). Otherwise, VecNorm() is called.

535:    Level: developer

537: .seealso: BVNorm(), BVNormColumn(), BVSetMatrix()
538: @*/
539: PetscErrorCode BVNormVec(BV bv,Vec v,NormType type,PetscReal *val)
540: {
541:   PetscInt       n;

543:   PetscFunctionBegin;
547:   PetscAssertPointer(val,4);
549:   BVCheckSizes(bv,1);
551:   PetscCheckSameComm(bv,1,v,2);

553:   PetscCheck(type!=NORM_1_AND_2 || bv->matrix,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Requested norm not available");

555:   PetscCall(PetscLogEventBegin(BV_NormVec,bv,0,0,0));
556:   if (bv->matrix) { /* non-standard inner product */
557:     PetscCall(VecGetLocalSize(v,&n));
558:     PetscCheck(bv->n==n,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension bv %" PetscInt_FMT ", v %" PetscInt_FMT,bv->n,n);
559:     PetscCall(BVNorm_Private(bv,v,type,val));
560:   } else PetscCall(VecNorm(v,type,val));
561:   PetscCall(PetscLogEventEnd(BV_NormVec,bv,0,0,0));
562:   PetscFunctionReturn(PETSC_SUCCESS);
563: }

565: /*@
566:    BVNormVecBegin - Starts a split phase norm computation.

568:    Input Parameters:
569: +  bv   - basis vectors
570: .  v    - the vector
571: .  type - the norm type
572: -  val  - the norm

574:    Note:
575:    Each call to BVNormVecBegin() should be paired with a call to BVNormVecEnd().

577:    Level: advanced

579: .seealso: BVNormVecEnd(), BVNormVec()
580: @*/
581: PetscErrorCode BVNormVecBegin(BV bv,Vec v,NormType type,PetscReal *val)
582: {
583:   PetscInt       n;

585:   PetscFunctionBegin;
589:   PetscAssertPointer(val,4);
591:   BVCheckSizes(bv,1);
593:   PetscCheckSameTypeAndComm(bv,1,v,2);

595:   PetscCheck(type!=NORM_1_AND_2 || bv->matrix,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Requested norm not available");

597:   PetscCall(PetscLogEventBegin(BV_NormVec,bv,0,0,0));
598:   if (bv->matrix) { /* non-standard inner product */
599:     PetscCall(VecGetLocalSize(v,&n));
600:     PetscCheck(bv->n==n,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension bv %" PetscInt_FMT ", v %" PetscInt_FMT,bv->n,n);
601:     PetscCall(BVNorm_Begin_Private(bv,v,type,val));
602:   } else PetscCall(VecNormBegin(v,type,val));
603:   PetscCall(PetscLogEventEnd(BV_NormVec,bv,0,0,0));
604:   PetscFunctionReturn(PETSC_SUCCESS);
605: }

607: /*@
608:    BVNormVecEnd - Ends a split phase norm computation.

610:    Input Parameters:
611: +  bv   - basis vectors
612: .  v    - the vector
613: .  type - the norm type
614: -  val  - the norm

616:    Note:
617:    Each call to BVNormVecBegin() should be paired with a call to BVNormVecEnd().

619:    Level: advanced

621: .seealso: BVNormVecBegin(), BVNormVec()
622: @*/
623: PetscErrorCode BVNormVecEnd(BV bv,Vec v,NormType type,PetscReal *val)
624: {
625:   PetscFunctionBegin;
628:   PetscAssertPointer(val,4);
630:   BVCheckSizes(bv,1);

632:   PetscCheck(type!=NORM_1_AND_2,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Requested norm not available");

634:   if (bv->matrix) PetscCall(BVNorm_End_Private(bv,v,type,val));  /* non-standard inner product */
635:   else PetscCall(VecNormEnd(v,type,val));
636:   PetscFunctionReturn(PETSC_SUCCESS);
637: }

639: /*@
640:    BVNormColumn - Computes the vector norm of a selected column.

642:    Collective

644:    Input Parameters:
645: +  bv   - basis vectors
646: .  j    - column number to be used
647: -  type - the norm type

649:    Output Parameter:
650: .  val  - the norm

652:    Notes:
653:    The norm of V[j] is computed (NORM_1, NORM_2, or NORM_INFINITY).
654:    If a non-standard inner product has been specified with BVSetMatrix(),
655:    then the returned value is sqrt(V[j]'*B*V[j]),
656:    where B is the inner product matrix (argument 'type' is ignored).

658:    Level: intermediate

660: .seealso: BVNorm(), BVNormVec(), BVNormalize(), BVSetActiveColumns(), BVSetMatrix()
661: @*/
662: PetscErrorCode BVNormColumn(BV bv,PetscInt j,NormType type,PetscReal *val)
663: {
664:   Vec            z;

666:   PetscFunctionBegin;
670:   PetscAssertPointer(val,4);
672:   BVCheckSizes(bv,1);

674:   PetscCheck(j>=0 && j<bv->m,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Argument j has wrong value %" PetscInt_FMT ", the number of columns is %" PetscInt_FMT,j,bv->m);
675:   PetscCheck(type!=NORM_1_AND_2 || bv->matrix,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Requested norm not available");

677:   PetscCall(PetscLogEventBegin(BV_NormVec,bv,0,0,0));
678:   if (bv->matrix) { /* non-standard inner product */
679:     PetscCall(BVGetColumn(bv,j,&z));
680:     PetscCall(BVNorm_Private(bv,z,type,val));
681:     PetscCall(BVRestoreColumn(bv,j,&z));
682:   } else PetscUseTypeMethod(bv,norm,j,type,val);
683:   PetscCall(PetscLogEventEnd(BV_NormVec,bv,0,0,0));
684:   PetscFunctionReturn(PETSC_SUCCESS);
685: }

687: /*@
688:    BVNormColumnBegin - Starts a split phase norm computation.

690:    Input Parameters:
691: +  bv   - basis vectors
692: .  j    - column number to be used
693: .  type - the norm type
694: -  val  - the norm

696:    Note:
697:    Each call to BVNormColumnBegin() should be paired with a call to BVNormColumnEnd().

699:    Level: advanced

701: .seealso: BVNormColumnEnd(), BVNormColumn()
702: @*/
703: PetscErrorCode BVNormColumnBegin(BV bv,PetscInt j,NormType type,PetscReal *val)
704: {
705:   PetscSplitReduction *sr;
706:   PetscReal           lresult;
707:   MPI_Comm            comm;
708:   Vec                 z;

710:   PetscFunctionBegin;
714:   PetscAssertPointer(val,4);
716:   BVCheckSizes(bv,1);

718:   PetscCheck(j>=0 && j<bv->m,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Argument j has wrong value %" PetscInt_FMT ", the number of columns is %" PetscInt_FMT,j,bv->m);
719:   PetscCheck(type!=NORM_1_AND_2 || bv->matrix,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Requested norm not available");

721:   PetscCall(PetscLogEventBegin(BV_NormVec,bv,0,0,0));
722:   PetscCall(BVGetColumn(bv,j,&z));
723:   if (bv->matrix) PetscCall(BVNorm_Begin_Private(bv,z,type,val)); /* non-standard inner product */
724:   else if (bv->ops->norm_begin) PetscUseTypeMethod(bv,norm_begin,j,type,val);
725:   else {
726:     BVCheckOp(bv,1,norm_local);
727:     PetscCall(PetscObjectGetComm((PetscObject)z,&comm));
728:     PetscCall(PetscSplitReductionGet(comm,&sr));
729:     PetscCheck(sr->state==STATE_BEGIN,PetscObjectComm((PetscObject)bv),PETSC_ERR_ORDER,"Called before all VecxxxEnd() called");
730:     if (sr->numopsbegin >= sr->maxops) PetscCall(PetscSplitReductionExtend(sr));
731:     sr->invecs[sr->numopsbegin] = (void*)bv;
732:     PetscUseTypeMethod(bv,norm_local,j,type,&lresult);
733:     if (type == NORM_2) lresult = lresult*lresult;
734:     if (type == NORM_MAX) sr->reducetype[sr->numopsbegin] = PETSC_SR_REDUCE_MAX;
735:     else sr->reducetype[sr->numopsbegin] = PETSC_SR_REDUCE_SUM;
736:     sr->lvalues[sr->numopsbegin++] = lresult;
737:   }
738:   PetscCall(BVRestoreColumn(bv,j,&z));
739:   PetscCall(PetscLogEventEnd(BV_NormVec,bv,0,0,0));
740:   PetscFunctionReturn(PETSC_SUCCESS);
741: }

743: /*@
744:    BVNormColumnEnd - Ends a split phase norm computation.

746:    Input Parameters:
747: +  bv   - basis vectors
748: .  j    - column number to be used
749: .  type - the norm type
750: -  val  - the norm

752:    Note:
753:    Each call to BVNormColumnBegin() should be paired with a call to BVNormColumnEnd().

755:    Level: advanced

757: .seealso: BVNormColumnBegin(), BVNormColumn()
758: @*/
759: PetscErrorCode BVNormColumnEnd(BV bv,PetscInt j,NormType type,PetscReal *val)
760: {
761:   PetscSplitReduction *sr;
762:   MPI_Comm            comm;
763:   Vec                 z;

765:   PetscFunctionBegin;
769:   PetscAssertPointer(val,4);
771:   BVCheckSizes(bv,1);

773:   PetscCheck(type!=NORM_1_AND_2,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Requested norm not available");

775:   PetscCall(BVGetColumn(bv,j,&z));
776:   if (bv->matrix) PetscCall(BVNorm_End_Private(bv,z,type,val)); /* non-standard inner product */
777:   else if (bv->ops->norm_end) PetscUseTypeMethod(bv,norm_end,j,type,val);
778:   else {
779:     PetscCall(PetscObjectGetComm((PetscObject)z,&comm));
780:     PetscCall(PetscSplitReductionGet(comm,&sr));
781:     PetscCall(PetscSplitReductionEnd(sr));

783:     PetscCheck(sr->numopsend<sr->numopsbegin,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Called VecxxxEnd() more times then VecxxxBegin()");
784:     PetscCheck((void*)bv==sr->invecs[sr->numopsend],PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Called VecxxxEnd() in a different order or with a different vector than VecxxxBegin()");
785:     PetscCheck(sr->reducetype[sr->numopsend]==PETSC_SR_REDUCE_MAX || type!=NORM_MAX,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Called BVNormEnd(,NORM_MAX,) on a reduction started with VecDotBegin() or NORM_1 or NORM_2");
786:     *val = PetscRealPart(sr->gvalues[sr->numopsend++]);
787:     if (type == NORM_2) *val = PetscSqrtReal(*val);
788:     if (sr->numopsend == sr->numopsbegin) {
789:       sr->state       = STATE_BEGIN;
790:       sr->numopsend   = 0;
791:       sr->numopsbegin = 0;
792:     }
793:   }
794:   PetscCall(BVRestoreColumn(bv,j,&z));
795:   PetscFunctionReturn(PETSC_SUCCESS);
796: }

798: /*@
799:    BVNormalize - Normalize all columns (starting from the leading ones).

801:    Collective

803:    Input Parameters:
804: +  bv   - basis vectors
805: -  eigi - (optional) imaginary parts of eigenvalues

807:    Notes:
808:    On output, all columns will have unit norm. The normalization is done with
809:    respect to the 2-norm, or to the B-norm if a non-standard inner product has
810:    been specified with BVSetMatrix(), see BVNormColumn().

812:    If the optional argument eigi is passed (taken into account only in real
813:    scalars) it is interpreted as the imaginary parts of the eigenvalues and
814:    the BV is supposed to contain the corresponding eigenvectors. Suppose the
815:    first three values are eigi = { 0, alpha, -alpha }, then the first column
816:    is normalized as usual, but the second and third ones are normalized assuming
817:    that they contain the real and imaginary parts of a complex conjugate pair of
818:    eigenvectors.

820:    If eigi is passed, the inner-product matrix is ignored.

822:    If there are leading columns, they are not modified (are assumed to be already
823:    normalized).

825:    Level: intermediate

827: .seealso: BVNormColumn()
828: @*/
829: PetscErrorCode BVNormalize(BV bv,PetscScalar *eigi)
830: {
831:   PetscReal      norm;
832:   PetscInt       i;

834:   PetscFunctionBegin;
837:   BVCheckSizes(bv,1);

839:   PetscCall(PetscLogEventBegin(BV_Normalize,bv,0,0,0));
840:   if (bv->matrix && !eigi) {
841:     for (i=bv->l;i<bv->k;i++) {
842:       PetscCall(BVNormColumn(bv,i,NORM_2,&norm));
843:       PetscCall(BVScaleColumn(bv,i,1.0/norm));
844:     }
845:   } else PetscTryTypeMethod(bv,normalize,eigi);
846:   PetscCall(PetscLogEventEnd(BV_Normalize,bv,0,0,0));
847:   PetscCall(PetscObjectStateIncrease((PetscObject)bv));
848:   PetscFunctionReturn(PETSC_SUCCESS);
849: }

851: /*
852:   Compute Y^H*A*X: right part column by column (with MatMult) and bottom
853:   part row by row (with MatMultHermitianTranspose); result placed in marray[*,ldm]
854: */
855: static inline PetscErrorCode BVMatProject_Vec(BV X,Mat A,BV Y,PetscScalar *marray,PetscInt ldm,PetscBool symm)
856: {
857:   PetscInt       i,j,lx,ly,kx,ky,ulim;
858:   Vec            z,f;

860:   PetscFunctionBegin;
861:   lx = X->l; kx = X->k;
862:   ly = Y->l; ky = Y->k;
863:   PetscCall(BVCreateVec(X,&f));
864:   BVCheckOp(Y,3,dotvec);
865:   for (j=lx;j<kx;j++) {
866:     PetscCall(BVGetColumn(X,j,&z));
867:     PetscCall(MatMult(A,z,f));
868:     PetscCall(BVRestoreColumn(X,j,&z));
869:     ulim = PetscMin(ly+(j-lx)+1,ky);
870:     Y->l = 0; Y->k = ulim;
871:     PetscUseTypeMethod(Y,dotvec,f,marray+j*ldm);
872:     if (symm) {
873:       for (i=0;i<j;i++) marray[j+i*ldm] = PetscConj(marray[i+j*ldm]);
874:     }
875:   }
876:   if (!symm) {
877:     BVCheckOp(X,1,dotvec);
878:     PetscCall(BV_AllocateCoeffs(Y));
879:     for (j=ly;j<ky;j++) {
880:       PetscCall(BVGetColumn(Y,j,&z));
881:       PetscCall(MatMultHermitianTranspose(A,z,f));
882:       PetscCall(BVRestoreColumn(Y,j,&z));
883:       ulim = PetscMin(lx+(j-ly),kx);
884:       X->l = 0; X->k = ulim;
885:       PetscUseTypeMethod(X,dotvec,f,Y->h);
886:       for (i=0;i<ulim;i++) marray[j+i*ldm] = PetscConj(Y->h[i]);
887:     }
888:   }
889:   PetscCall(VecDestroy(&f));
890:   X->l = lx; X->k = kx;
891:   Y->l = ly; Y->k = ky;
892:   PetscFunctionReturn(PETSC_SUCCESS);
893: }

895: /*
896:   Compute Y^H*A*X= [   --   | Y0'*W1 ]
897:                    [ Y1'*W0 | Y1'*W1 ]
898:   Allocates auxiliary BV to store the result of A*X, then one BVDot
899:   call for top-right part and another one for bottom part;
900:   result placed in marray[*,ldm]
901: */
902: static inline PetscErrorCode BVMatProject_MatMult(BV X,Mat A,BV Y,PetscScalar *marray,PetscInt ldm)
903: {
904:   PetscInt          j,lx,ly,kx,ky;
905:   const PetscScalar *harray;
906:   Mat               H;
907:   BV                W;

909:   PetscFunctionBegin;
910:   lx = X->l; kx = X->k;
911:   ly = Y->l; ky = Y->k;
912:   PetscCall(BVDuplicate(X,&W));
913:   X->l = 0; X->k = kx;
914:   W->l = 0; W->k = kx;
915:   PetscCall(BVMatMult(X,A,W));

917:   /* top-right part, Y0'*AX1 */
918:   if (ly>0 && lx<kx) {
919:     PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,ly,kx,NULL,&H));
920:     W->l = lx; W->k = kx;
921:     Y->l = 0;  Y->k = ly;
922:     PetscCall(BVDot(W,Y,H));
923:     PetscCall(MatDenseGetArrayRead(H,&harray));
924:     for (j=lx;j<kx;j++) PetscCall(PetscArraycpy(marray+j*ldm,harray+j*ly,ly));
925:     PetscCall(MatDenseRestoreArrayRead(H,&harray));
926:     PetscCall(MatDestroy(&H));
927:   }

929:   /* bottom part, Y1'*AX */
930:   if (kx>0 && ly<ky) {
931:     PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,ky,kx,NULL,&H));
932:     W->l = 0;  W->k = kx;
933:     Y->l = ly; Y->k = ky;
934:     PetscCall(BVDot(W,Y,H));
935:     PetscCall(MatDenseGetArrayRead(H,&harray));
936:     for (j=0;j<kx;j++) PetscCall(PetscArraycpy(marray+j*ldm+ly,harray+j*ky+ly,ky-ly));
937:     PetscCall(MatDenseRestoreArrayRead(H,&harray));
938:     PetscCall(MatDestroy(&H));
939:   }
940:   PetscCall(BVDestroy(&W));
941:   X->l = lx; X->k = kx;
942:   Y->l = ly; Y->k = ky;
943:   PetscFunctionReturn(PETSC_SUCCESS);
944: }

946: /*
947:   Compute Y^H*A*X= [   --   | Y0'*W1 ]
948:                    [ Y1'*W0 | Y1'*W1 ]
949:   First stage: allocate auxiliary BV to store A*X1, one BVDot for right part;
950:   Second stage: resize BV to accommodate A'*Y1, then call BVDot for transpose of
951:   bottom-left part; result placed in marray[*,ldm]
952: */
953: static inline PetscErrorCode BVMatProject_MatMult_2(BV X,Mat A,BV Y,PetscScalar *marray,PetscInt ldm,PetscBool symm)
954: {
955:   PetscInt          i,j,lx,ly,kx,ky;
956:   const PetscScalar *harray;
957:   Mat               H;
958:   BV                W;

960:   PetscFunctionBegin;
961:   lx = X->l; kx = X->k;
962:   ly = Y->l; ky = Y->k;

964:   /* right part, Y'*AX1 */
965:   PetscCall(BVDuplicateResize(X,kx-lx,&W));
966:   if (ky>0 && lx<kx) {
967:     PetscCall(BVMatMult(X,A,W));
968:     PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,ky,kx-lx,NULL,&H));
969:     Y->l = 0; Y->k = ky;
970:     PetscCall(BVDot(W,Y,H));
971:     PetscCall(MatDenseGetArrayRead(H,&harray));
972:     for (j=lx;j<kx;j++) PetscCall(PetscArraycpy(marray+j*ldm,harray+(j-lx)*ky,ky));
973:     PetscCall(MatDenseRestoreArrayRead(H,&harray));
974:     PetscCall(MatDestroy(&H));
975:   }

977:   /* bottom-left part, Y1'*AX0 */
978:   if (lx>0 && ly<ky) {
979:     if (symm) {
980:       /* do not compute, just copy symmetric elements */
981:       for (i=ly;i<ky;i++) {
982:         for (j=0;j<lx;j++) marray[i+j*ldm] = PetscConj(marray[j+i*ldm]);
983:       }
984:     } else {
985:       PetscCall(BVResize(W,ky-ly,PETSC_FALSE));
986:       Y->l = ly; Y->k = ky;
987:       PetscCall(BVMatMultHermitianTranspose(Y,A,W));
988:       PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,lx,ky-ly,NULL,&H));
989:       X->l = 0; X->k = lx;
990:       PetscCall(BVDot(W,X,H));
991:       PetscCall(MatDenseGetArrayRead(H,&harray));
992:       for (i=0;i<ky-ly;i++) {
993:         for (j=0;j<lx;j++) {
994:           marray[i+j*ldm+ly] = PetscConj(harray[j+i*(ky-ly)]);
995:         }
996:       }
997:       PetscCall(MatDenseRestoreArrayRead(H,&harray));
998:       PetscCall(MatDestroy(&H));
999:     }
1000:   }
1001:   PetscCall(BVDestroy(&W));
1002:   X->l = lx; X->k = kx;
1003:   Y->l = ly; Y->k = ky;
1004:   PetscFunctionReturn(PETSC_SUCCESS);
1005: }

1007: /*
1008:   Compute Y^H*X = [   --   | Y0'*X1 ]     (X contains A*X):
1009:                   [ Y1'*X0 | Y1'*X1 ]
1010:   one BVDot call for top-right part and another one for bottom part;
1011:   result placed in marray[*,ldm]
1012: */
1013: static inline PetscErrorCode BVMatProject_Dot(BV X,BV Y,PetscScalar *marray,PetscInt ldm)
1014: {
1015:   PetscInt          j,lx,ly,kx,ky;
1016:   const PetscScalar *harray;
1017:   Mat               H;

1019:   PetscFunctionBegin;
1020:   lx = X->l; kx = X->k;
1021:   ly = Y->l; ky = Y->k;

1023:   /* top-right part, Y0'*X1 */
1024:   if (ly>0 && lx<kx) {
1025:     PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,ly,kx,NULL,&H));
1026:     X->l = lx; X->k = kx;
1027:     Y->l = 0;  Y->k = ly;
1028:     PetscCall(BVDot(X,Y,H));
1029:     PetscCall(MatDenseGetArrayRead(H,&harray));
1030:     for (j=lx;j<kx;j++) PetscCall(PetscArraycpy(marray+j*ldm,harray+j*ly,ly));
1031:     PetscCall(MatDenseRestoreArrayRead(H,&harray));
1032:     PetscCall(MatDestroy(&H));
1033:   }

1035:   /* bottom part, Y1'*X */
1036:   if (kx>0 && ly<ky) {
1037:     PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,ky,kx,NULL,&H));
1038:     X->l = 0;  X->k = kx;
1039:     Y->l = ly; Y->k = ky;
1040:     PetscCall(BVDot(X,Y,H));
1041:     PetscCall(MatDenseGetArrayRead(H,&harray));
1042:     for (j=0;j<kx;j++) PetscCall(PetscArraycpy(marray+j*ldm+ly,harray+j*ky+ly,ky-ly));
1043:     PetscCall(MatDenseRestoreArrayRead(H,&harray));
1044:     PetscCall(MatDestroy(&H));
1045:   }
1046:   X->l = lx; X->k = kx;
1047:   Y->l = ly; Y->k = ky;
1048:   PetscFunctionReturn(PETSC_SUCCESS);
1049: }

1051: /*@
1052:    BVMatProject - Computes the projection of a matrix onto a subspace.

1054:    Collective

1056:    Input Parameters:
1057: +  X - basis vectors
1058: .  A - (optional) matrix to be projected
1059: -  Y - left basis vectors, can be equal to X

1061:    Output Parameter:
1062: .  M - the resulting matrix

1064:    Notes:
1065:    If A=NULL, then it is assumed that X already contains A*X.

1067:    This operation is similar to BVDot(), with important differences.
1068:    The goal is to compute the matrix resulting from the orthogonal projection
1069:    of A onto the subspace spanned by the columns of X, M = X^H*A*X, or the
1070:    oblique projection onto X along Y, M = Y^H*A*X.

1072:    A difference with respect to BVDot() is that the standard inner product
1073:    is always used, regardless of a non-standard inner product being specified
1074:    with BVSetMatrix().

1076:    On entry, M must be a sequential dense Mat with dimensions ky,kx at least,
1077:    where ky (resp. kx) is the number of active columns of Y (resp. X).
1078:    Another difference with respect to BVDot() is that all entries of M are
1079:    computed except the leading ly,lx part, where ly (resp. lx) is the
1080:    number of leading columns of Y (resp. X). Hence, the leading columns of
1081:    X and Y participate in the computation, as opposed to BVDot().
1082:    The leading part of M is assumed to be already available from previous
1083:    computations.

1085:    In the orthogonal projection case, Y=X, some computation can be saved if
1086:    A is real symmetric (or complex Hermitian). In order to exploit this
1087:    property, the symmetry flag of A must be set with MatSetOption().

1089:    Level: intermediate

1091: .seealso: BVDot(), BVSetActiveColumns(), BVSetMatrix()
1092: @*/
1093: PetscErrorCode BVMatProject(BV X,Mat A,BV Y,Mat M)
1094: {
1095:   PetscBool      set,flg,symm=PETSC_FALSE;
1096:   PetscInt       m,n,ldm;
1097:   PetscScalar    *marray;
1098:   Mat            Xmatrix,Ymatrix;
1099:   PetscObjectId  idx,idy;

1101:   PetscFunctionBegin;
1107:   BVCheckSizes(X,1);
1108:   if (A) {
1110:     PetscCheckSameComm(X,1,A,2);
1111:   }
1113:   BVCheckSizes(Y,3);
1115:   PetscCheckSameTypeAndComm(X,1,Y,3);
1116:   SlepcMatCheckSeq(M);

1118:   PetscCall(MatGetSize(M,&m,&n));
1119:   PetscCall(MatDenseGetLDA(M,&ldm));
1120:   PetscCheck(m>=Y->k,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_SIZ,"Matrix M has %" PetscInt_FMT " rows, should have at least %" PetscInt_FMT,m,Y->k);
1121:   PetscCheck(n>=X->k,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_SIZ,"Matrix M has %" PetscInt_FMT " columns, should have at least %" PetscInt_FMT,n,X->k);
1122:   PetscCheck(X->n==Y->n,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension X %" PetscInt_FMT ", Y %" PetscInt_FMT,X->n,Y->n);

1124:   PetscCall(PetscLogEventBegin(BV_MatProject,X,A,Y,0));
1125:   /* temporarily set standard inner product */
1126:   Xmatrix = X->matrix;
1127:   Ymatrix = Y->matrix;
1128:   X->matrix = Y->matrix = NULL;

1130:   PetscCall(PetscObjectGetId((PetscObject)X,&idx));
1131:   PetscCall(PetscObjectGetId((PetscObject)Y,&idy));
1132:   if (A && idx==idy) { /* check symmetry of M=X'AX */
1133:     PetscCall(MatIsHermitianKnown(A,&set,&flg));
1134:     symm = set? flg: PETSC_FALSE;
1135:   }

1137:   PetscCall(MatDenseGetArray(M,&marray));

1139:   if (A) {
1140:     if (X->vmm==BV_MATMULT_VECS) {
1141:       /* perform computation column by column */
1142:       PetscCall(BVMatProject_Vec(X,A,Y,marray,ldm,symm));
1143:     } else {
1144:       /* use BVMatMult, then BVDot */
1145:       PetscCall(MatHasOperation(A,MATOP_MULT_TRANSPOSE,&flg));
1146:       if (symm || (flg && X->l>=X->k/2 && Y->l>=Y->k/2)) PetscCall(BVMatProject_MatMult_2(X,A,Y,marray,ldm,symm));
1147:       else PetscCall(BVMatProject_MatMult(X,A,Y,marray,ldm));
1148:     }
1149:   } else {
1150:     /* use BVDot on subblocks */
1151:     PetscCall(BVMatProject_Dot(X,Y,marray,ldm));
1152:   }

1154:   PetscCall(MatDenseRestoreArray(M,&marray));
1155:   PetscCall(PetscLogEventEnd(BV_MatProject,X,A,Y,0));
1156:   /* restore non-standard inner product */
1157:   X->matrix = Xmatrix;
1158:   Y->matrix = Ymatrix;
1159:   PetscFunctionReturn(PETSC_SUCCESS);
1160: }