Actual source code: bvglobal.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 involving global communication
 12: */

 14: #include <slepc/private/bvimpl.h>      /*I "slepcbv.h" I*/

 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: PETSC_STATIC_INLINE PetscErrorCode BVDot_Private(BV X,BV Y,Mat M)
 21: {
 23:   PetscObjectId  idx,idy;
 24:   PetscInt       i,j,jend,m;
 25:   PetscScalar    *marray;
 26:   PetscBool      symm=PETSC_FALSE;
 27:   Mat            B;
 28:   Vec            z;

 31:   MatGetSize(M,&m,NULL);
 32:   MatDenseGetArray(M,&marray);
 33:   PetscObjectGetId((PetscObject)X,&idx);
 34:   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:     BVGetColumn(X->cached,j,&z);
 42:     (*Y->ops->dotvec)(Y,z,marray+j*m+Y->l);
 43:     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:   MatDenseRestoreArray(M,&marray);
 50:   Y->matrix = B;
 51:   return(0);
 52: }

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

 57:    Collective on BV

 59:    Input Parameters:
 60: +  X, Y - basis vectors
 61: -  M    - Mat object where the result must be placed

 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: {
 89:   PetscBool      match;
 90:   PetscInt       m,n;

 97:   BVCheckSizes(X,1);
 99:   BVCheckSizes(Y,2);
102:   PetscObjectTypeCompare((PetscObject)M,MATSEQDENSE,&match);
103:   if (!match) SETERRQ(PetscObjectComm((PetscObject)X),PETSC_ERR_SUP,"Mat argument must be of type seqdense");

105:   MatGetSize(M,&m,&n);
106:   if (m<Y->k) SETERRQ2(PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_SIZ,"Mat argument has %D rows, should have at least %D",m,Y->k);
107:   if (n<X->k) SETERRQ2(PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_SIZ,"Mat argument has %D columns, should have at least %D",n,X->k);
108:   if (X->n!=Y->n) SETERRQ2(PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension X %D, Y %D",X->n,Y->n);
109:   if (X->matrix!=Y->matrix) SETERRQ(PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_WRONGSTATE,"X and Y must have the same inner product matrix");
110:   if (X->l==X->k || Y->l==Y->k) return(0);

112:   PetscLogEventBegin(BV_Dot,X,Y,0,0);
113:   if (X->matrix) { /* non-standard inner product */
114:     /* compute BX first */
115:     BV_IPMatMultBV(X);
116:     if (X->vmm==BV_MATMULT_VECS) {
117:       /* perform computation column by column */
118:       BVDot_Private(X,Y,M);
119:     } else {
120:       (*X->ops->dot)(X->cached,Y,M);
121:     }
122:   } else {
123:     (*X->ops->dot)(X,Y,M);
124:   }
125:   PetscLogEventEnd(BV_Dot,X,Y,0,0);
126:   return(0);
127: }

129: /*@
130:    BVDotVec - Computes multiple dot products of a vector against all the
131:    column vectors of a BV.

133:    Collective on BV and Vec

135:    Input Parameters:
136: +  X - basis vectors
137: -  y - a vector

139:    Output Parameter:
140: .  m - an array where the result must be placed

142:    Notes:
143:    This is analogue to VecMDot(), but using BV to represent a collection
144:    of vectors. The result is m = X^H*y, so m_i is equal to x_j^H y. Note
145:    that here X is transposed as opposed to BVDot().

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

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

154:    Level: intermediate

156: .seealso: BVDot(), BVDotColumn(), BVSetActiveColumns(), BVSetMatrix()
157: @*/
158: PetscErrorCode BVDotVec(BV X,Vec y,PetscScalar m[])
159: {
161:   PetscInt       n;

167:   BVCheckSizes(X,1);

171:   VecGetLocalSize(y,&n);
172:   if (X->n!=n) SETERRQ2(PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension X %D, y %D",X->n,n);

174:   PetscLogEventBegin(BV_DotVec,X,y,0,0);
175:   (*X->ops->dotvec)(X,y,m);
176:   PetscLogEventEnd(BV_DotVec,X,y,0,0);
177:   return(0);
178: }

180: /*@
181:    BVDotVecBegin - Starts a split phase dot product computation.

183:    Input Parameters:
184: +  X - basis vectors
185: .  y - a vector
186: -  m - an array where the result will go (can be NULL)

188:    Note:
189:    Each call to BVDotVecBegin() should be paired with a call to BVDotVecEnd().

191:    Level: advanced

193: .seealso: BVDotVecEnd(), BVDotVec()
194: @*/
195: PetscErrorCode BVDotVecBegin(BV X,Vec y,PetscScalar *m)
196: {
197:   PetscErrorCode      ierr;
198:   PetscInt            i,n,nv;
199:   PetscSplitReduction *sr;
200:   MPI_Comm            comm;

206:   BVCheckSizes(X,1);

210:   VecGetLocalSize(y,&n);
211:   if (X->n!=n) SETERRQ2(PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension X %D, y %D",X->n,n);

213:   if (X->ops->dotvec_begin) {
214:     (*X->ops->dotvec_begin)(X,y,m);
215:   } else {
216:     nv = X->k-X->l;
217:     PetscObjectGetComm((PetscObject)X,&comm);
218:     PetscSplitReductionGet(comm,&sr);
219:     if (sr->state != STATE_BEGIN) SETERRQ(PetscObjectComm((PetscObject)X),PETSC_ERR_ORDER,"Called before all VecxxxEnd() called");
220:     for (i=0;i<nv;i++) {
221:       if (sr->numopsbegin+i >= sr->maxops) {
222:         PetscSplitReductionExtend(sr);
223:       }
224:       sr->reducetype[sr->numopsbegin+i] = PETSC_SR_REDUCE_SUM;
225:       sr->invecs[sr->numopsbegin+i]     = (void*)X;
226:     }
227:     PetscLogEventBegin(BV_DotVec,X,y,0,0);
228:     (*X->ops->dotvec_local)(X,y,sr->lvalues+sr->numopsbegin);
229:     sr->numopsbegin += nv;
230:     PetscLogEventEnd(BV_DotVec,X,y,0,0);
231:   }
232:   return(0);
233: }

235: /*@
236:    BVDotVecEnd - Ends a split phase dot product computation.

238:    Input Parameters:
239: +  X - basis vectors
240: .  y - a vector
241: -  m - an array where the result will go

243:    Note:
244:    Each call to BVDotVecBegin() should be paired with a call to BVDotVecEnd().

246:    Level: advanced

248: .seealso: BVDotVecBegin(), BVDotVec()
249: @*/
250: PetscErrorCode BVDotVecEnd(BV X,Vec y,PetscScalar *m)
251: {
252:   PetscErrorCode      ierr;
253:   PetscInt            i,nv;
254:   PetscSplitReduction *sr;
255:   MPI_Comm            comm;

260:   BVCheckSizes(X,1);

262:   if (X->ops->dotvec_end) {
263:     (*X->ops->dotvec_end)(X,y,m);
264:   } else {
265:     nv = X->k-X->l;
266:     PetscObjectGetComm((PetscObject)X,&comm);
267:     PetscSplitReductionGet(comm,&sr);
268:     PetscSplitReductionEnd(sr);

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

275:     /* Finished getting all the results so reset to no outstanding requests */
276:     if (sr->numopsend == sr->numopsbegin) {
277:       sr->state       = STATE_BEGIN;
278:       sr->numopsend   = 0;
279:       sr->numopsbegin = 0;
280:     }
281:   }
282:   return(0);
283: }

285: /*@
286:    BVDotColumn - Computes multiple dot products of a column against all the
287:    previous columns of a BV.

289:    Collective on BV

291:    Input Parameters:
292: +  X - basis vectors
293: -  j - the column index

295:    Output Parameter:
296: .  q - an array where the result must be placed

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

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

310:    Level: advanced

312: .seealso: BVDot(), BVDotVec(), BVSetActiveColumns(), BVSetMatrix()
313: @*/
314: PetscErrorCode BVDotColumn(BV X,PetscInt j,PetscScalar *q)
315: {
317:   PetscInt       ksave;
318:   Vec            y;

324:   BVCheckSizes(X,1);

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

329:   PetscLogEventBegin(BV_DotVec,X,0,0,0);
330:   ksave = X->k;
331:   X->k = j;
332:   if (!q && !X->buffer) { BVGetBufferVec(X,&X->buffer); }
333:   BVGetColumn(X,j,&y);
334:   (*X->ops->dotvec)(X,y,q);
335:   BVRestoreColumn(X,j,&y);
336:   X->k = ksave;
337:   PetscLogEventEnd(BV_DotVec,X,0,0,0);
338:   return(0);
339: }

341: /*@
342:    BVDotColumnBegin - Starts a split phase dot product computation.

344:    Input Parameters:
345: +  X - basis vectors
346: -  j - the column index
347: -  m - an array where the result will go (can be NULL)

349:    Note:
350:    Each call to BVDotColumnBegin() should be paired with a call to BVDotColumnEnd().

352:    Level: advanced

354: .seealso: BVDotColumnEnd(), BVDotColumn()
355: @*/
356: PetscErrorCode BVDotColumnBegin(BV X,PetscInt j,PetscScalar *m)
357: {
358:   PetscErrorCode      ierr;
359:   PetscInt            i,nv,ksave;
360:   PetscSplitReduction *sr;
361:   MPI_Comm            comm;
362:   Vec                 y;

368:   BVCheckSizes(X,1);

370:   if (j<0) SETERRQ(PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
371:   if (j>=X->m) SETERRQ2(PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_OUTOFRANGE,"Index j=%D but BV only has %D columns",j,X->m);
372:   ksave = X->k;
373:   X->k = j;
374:   BVGetColumn(X,j,&y);

376:   if (X->ops->dotvec_begin) {
377:     (*X->ops->dotvec_begin)(X,y,m);
378:   } else {
379:     nv = X->k-X->l;
380:     PetscObjectGetComm((PetscObject)X,&comm);
381:     PetscSplitReductionGet(comm,&sr);
382:     if (sr->state != STATE_BEGIN) SETERRQ(PetscObjectComm((PetscObject)X),PETSC_ERR_ORDER,"Called before all VecxxxEnd() called");
383:     for (i=0;i<nv;i++) {
384:       if (sr->numopsbegin+i >= sr->maxops) {
385:         PetscSplitReductionExtend(sr);
386:       }
387:       sr->reducetype[sr->numopsbegin+i] = PETSC_SR_REDUCE_SUM;
388:       sr->invecs[sr->numopsbegin+i]     = (void*)X;
389:     }
390:     PetscLogEventBegin(BV_DotVec,X,0,0,0);
391:     (*X->ops->dotvec_local)(X,y,sr->lvalues+sr->numopsbegin);
392:     sr->numopsbegin += nv;
393:     PetscLogEventEnd(BV_DotVec,X,0,0,0);
394:   }
395:   BVRestoreColumn(X,j,&y);
396:   X->k = ksave;
397:   return(0);
398: }

400: /*@
401:    BVDotColumnEnd - Ends a split phase dot product computation.

403:    Input Parameters:
404: +  X - basis vectors
405: .  j - the column index
406: -  m - an array where the result will go

408:    Notes:
409:    Each call to BVDotColumnBegin() should be paired with a call to BVDotColumnEnd().

411:    Level: advanced

413: .seealso: BVDotColumnBegin(), BVDotColumn()
414: @*/
415: PetscErrorCode BVDotColumnEnd(BV X,PetscInt j,PetscScalar *m)
416: {
417:   PetscErrorCode      ierr;
418:   PetscInt            i,nv,ksave;
419:   PetscSplitReduction *sr;
420:   MPI_Comm            comm;
421:   Vec                 y;

427:   BVCheckSizes(X,1);

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

434:   if (X->ops->dotvec_end) {
435:     BVGetColumn(X,j,&y);
436:     (*X->ops->dotvec_end)(X,y,m);
437:     BVRestoreColumn(X,j,&y);
438:   } else {
439:     nv = X->k-X->l;
440:     PetscObjectGetComm((PetscObject)X,&comm);
441:     PetscSplitReductionGet(comm,&sr);
442:     PetscSplitReductionEnd(sr);

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

449:     /* Finished getting all the results so reset to no outstanding requests */
450:     if (sr->numopsend == sr->numopsbegin) {
451:       sr->state       = STATE_BEGIN;
452:       sr->numopsend   = 0;
453:       sr->numopsbegin = 0;
454:     }
455:   }
456:   X->k = ksave;
457:   return(0);
458: }

460: PETSC_STATIC_INLINE PetscErrorCode BVNorm_Private(BV bv,Vec z,NormType type,PetscReal *val)
461: {
463:   PetscScalar    p;

466:   BV_IPMatMult(bv,z);
467:   VecDot(bv->Bx,z,&p);
468:   BV_SafeSqrt(bv,p,val);
469:   return(0);
470: }

472: PETSC_STATIC_INLINE PetscErrorCode BVNorm_Begin_Private(BV bv,Vec z,NormType type,PetscReal *val)
473: {
475:   PetscScalar    p;

478:   BV_IPMatMult(bv,z);
479:   VecDotBegin(bv->Bx,z,&p);
480:   return(0);
481: }

483: PETSC_STATIC_INLINE PetscErrorCode BVNorm_End_Private(BV bv,Vec z,NormType type,PetscReal *val)
484: {
486:   PetscScalar    p;

489:   VecDotEnd(bv->Bx,z,&p);
490:   BV_SafeSqrt(bv,p,val);
491:   return(0);
492: }

494: /*@
495:    BVNorm - Computes the matrix norm of the BV.

497:    Collective on BV

499:    Input Parameters:
500: +  bv   - basis vectors
501: -  type - the norm type

503:    Output Parameter:
504: .  val  - the norm

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

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

513:    Level: intermediate

515: .seealso: BVNormVec(), BVNormColumn(), BVSetActiveColumns(), BVSetMatrix()
516: @*/
517: PetscErrorCode BVNorm(BV bv,NormType type,PetscReal *val)
518: {

526:   BVCheckSizes(bv,1);

528:   if (type==NORM_2 || type==NORM_1_AND_2) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Requested norm not available");
529:   if (bv->matrix) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Matrix norm not available for non-standard inner product");

531:   PetscLogEventBegin(BV_Norm,bv,0,0,0);
532:   (*bv->ops->norm)(bv,-1,type,val);
533:   PetscLogEventEnd(BV_Norm,bv,0,0,0);
534:   return(0);
535: }

537: /*@
538:    BVNormVec - Computes the norm of a given vector.

540:    Collective on BV

542:    Input Parameters:
543: +  bv   - basis vectors
544: .  v    - the vector
545: -  type - the norm type

547:    Output Parameter:
548: .  val  - the norm

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

556:    Level: developer

558: .seealso: BVNorm(), BVNormColumn(), BVSetMatrix()
559: @*/
560: PetscErrorCode BVNormVec(BV bv,Vec v,NormType type,PetscReal *val)
561: {
563:   PetscInt       n;

571:   BVCheckSizes(bv,1);

575:   if (type==NORM_1_AND_2 && !bv->matrix) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Requested norm not available");

577:   PetscLogEventBegin(BV_NormVec,bv,0,0,0);
578:   if (bv->matrix) { /* non-standard inner product */
579:     VecGetLocalSize(v,&n);
580:     if (bv->n!=n) SETERRQ2(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension bv %D, v %D",bv->n,n);
581:     BVNorm_Private(bv,v,type,val);
582:   } else {
583:     VecNorm(v,type,val);
584:   }
585:   PetscLogEventEnd(BV_NormVec,bv,0,0,0);
586:   return(0);
587: }

589: /*@
590:    BVNormVecBegin - Starts a split phase norm computation.

592:    Input Parameters:
593: +  bv   - basis vectors
594: .  v    - the vector
595: .  type - the norm type
596: -  val  - the norm

598:    Note:
599:    Each call to BVNormVecBegin() should be paired with a call to BVNormVecEnd().

601:    Level: advanced

603: .seealso: BVNormVecEnd(), BVNormVec()
604: @*/
605: PetscErrorCode BVNormVecBegin(BV bv,Vec v,NormType type,PetscReal *val)
606: {
608:   PetscInt       n;

616:   BVCheckSizes(bv,1);

620:   if (type==NORM_1_AND_2 && !bv->matrix) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Requested norm not available");

622:   PetscLogEventBegin(BV_NormVec,bv,0,0,0);
623:   if (bv->matrix) { /* non-standard inner product */
624:     VecGetLocalSize(v,&n);
625:     if (bv->n!=n) SETERRQ2(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension bv %D, v %D",bv->n,n);
626:     BVNorm_Begin_Private(bv,v,type,val);
627:   } else {
628:     VecNormBegin(v,type,val);
629:   }
630:   PetscLogEventEnd(BV_NormVec,bv,0,0,0);
631:   return(0);
632: }

634: /*@
635:    BVNormVecEnd - Ends a split phase norm computation.

637:    Input Parameters:
638: +  bv   - basis vectors
639: .  v    - the vector
640: .  type - the norm type
641: -  val  - the norm

643:    Note:
644:    Each call to BVNormVecBegin() should be paired with a call to BVNormVecEnd().

646:    Level: advanced

648: .seealso: BVNormVecBegin(), BVNormVec()
649: @*/
650: PetscErrorCode BVNormVecEnd(BV bv,Vec v,NormType type,PetscReal *val)
651: {

659:   BVCheckSizes(bv,1);

661:   if (type==NORM_1_AND_2) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Requested norm not available");

663:   if (bv->matrix) { /* non-standard inner product */
664:     BVNorm_End_Private(bv,v,type,val);
665:   } else {
666:     VecNormEnd(v,type,val);
667:   }
668:   return(0);
669: }

671: /*@
672:    BVNormColumn - Computes the vector norm of a selected column.

674:    Collective on BV

676:    Input Parameters:
677: +  bv   - basis vectors
678: .  j    - column number to be used
679: -  type - the norm type

681:    Output Parameter:
682: .  val  - the norm

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

690:    Level: intermediate

692: .seealso: BVNorm(), BVNormVec(), BVSetActiveColumns(), BVSetMatrix()
693: @*/
694: PetscErrorCode BVNormColumn(BV bv,PetscInt j,NormType type,PetscReal *val)
695: {
697:   Vec            z;

705:   BVCheckSizes(bv,1);

707:   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);
708:   if (type==NORM_1_AND_2 && !bv->matrix) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Requested norm not available");

710:   PetscLogEventBegin(BV_NormVec,bv,0,0,0);
711:   if (bv->matrix) { /* non-standard inner product */
712:     BVGetColumn(bv,j,&z);
713:     BVNorm_Private(bv,z,type,val);
714:     BVRestoreColumn(bv,j,&z);
715:   } else {
716:     (*bv->ops->norm)(bv,j,type,val);
717:   }
718:   PetscLogEventEnd(BV_NormVec,bv,0,0,0);
719:   return(0);
720: }

722: /*@
723:    BVNormColumnBegin - Starts a split phase norm computation.

725:    Input Parameters:
726: +  bv   - basis vectors
727: .  j    - column number to be used
728: .  type - the norm type
729: -  val  - the norm

731:    Note:
732:    Each call to BVNormColumnBegin() should be paired with a call to BVNormColumnEnd().

734:    Level: advanced

736: .seealso: BVNormColumnEnd(), BVNormColumn()
737: @*/
738: PetscErrorCode BVNormColumnBegin(BV bv,PetscInt j,NormType type,PetscReal *val)
739: {
740:   PetscErrorCode      ierr;
741:   PetscSplitReduction *sr;
742:   PetscReal           lresult;
743:   MPI_Comm            comm;
744:   Vec                 z;

752:   BVCheckSizes(bv,1);

754:   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);
755:   if (type==NORM_1_AND_2 && !bv->matrix) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Requested norm not available");

757:   PetscLogEventBegin(BV_NormVec,bv,0,0,0);
758:   BVGetColumn(bv,j,&z);
759:   if (bv->matrix) { /* non-standard inner product */
760:     BVNorm_Begin_Private(bv,z,type,val);
761:   } else if (bv->ops->norm_begin) {
762:     (*bv->ops->norm_begin)(bv,j,type,val);
763:   } else {
764:     PetscObjectGetComm((PetscObject)z,&comm);
765:     PetscSplitReductionGet(comm,&sr);
766:     if (sr->state != STATE_BEGIN) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ORDER,"Called before all VecxxxEnd() called");
767:     if (sr->numopsbegin >= sr->maxops) {
768:       PetscSplitReductionExtend(sr);
769:     }
770:     sr->invecs[sr->numopsbegin] = (void*)bv;
771:     (*bv->ops->norm_local)(bv,j,type,&lresult);
772:     if (type == NORM_2) lresult = lresult*lresult;
773:     if (type == NORM_MAX) sr->reducetype[sr->numopsbegin] = PETSC_SR_REDUCE_MAX;
774:     else sr->reducetype[sr->numopsbegin] = PETSC_SR_REDUCE_SUM;
775:     sr->lvalues[sr->numopsbegin++] = lresult;
776:   }
777:   BVRestoreColumn(bv,j,&z);
778:   PetscLogEventEnd(BV_NormVec,bv,0,0,0);
779:   return(0);
780: }

782: /*@
783:    BVNormColumnEnd - Ends a split phase norm computation.

785:    Input Parameters:
786: +  bv   - basis vectors
787: .  j    - column number to be used
788: .  type - the norm type
789: -  val  - the norm

791:    Note:
792:    Each call to BVNormColumnBegin() should be paired with a call to BVNormColumnEnd().

794:    Level: advanced

796: .seealso: BVNormColumnBegin(), BVNormColumn()
797: @*/
798: PetscErrorCode BVNormColumnEnd(BV bv,PetscInt j,NormType type,PetscReal *val)
799: {
800:   PetscErrorCode      ierr;
801:   PetscSplitReduction *sr;
802:   MPI_Comm            comm;
803:   Vec                 z;

811:   BVCheckSizes(bv,1);

813:   if (type==NORM_1_AND_2) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Requested norm not available");

815:   BVGetColumn(bv,j,&z);
816:   if (bv->matrix) { /* non-standard inner product */
817:     BVNorm_End_Private(bv,z,type,val);
818:   } else if (bv->ops->norm_end) {
819:     (*bv->ops->norm_end)(bv,j,type,val);
820:   } else {
821:     PetscObjectGetComm((PetscObject)z,&comm);
822:     PetscSplitReductionGet(comm,&sr);
823:     PetscSplitReductionEnd(sr);

825:     if (sr->numopsend >= sr->numopsbegin) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Called VecxxxEnd() more times then VecxxxBegin()");
826:     if ((void*)bv != sr->invecs[sr->numopsend]) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Called VecxxxEnd() in a different order or with a different vector than VecxxxBegin()");
827:     if (sr->reducetype[sr->numopsend] != PETSC_SR_REDUCE_MAX && type == NORM_MAX) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Called BVNormEnd(,NORM_MAX,) on a reduction started with VecDotBegin() or NORM_1 or NORM_2");
828:     *val = PetscRealPart(sr->gvalues[sr->numopsend++]);
829:     if (type == NORM_2) *val = PetscSqrtReal(*val);
830:     if (sr->numopsend == sr->numopsbegin) {
831:       sr->state       = STATE_BEGIN;
832:       sr->numopsend   = 0;
833:       sr->numopsbegin = 0;
834:     }
835:   }
836:   BVRestoreColumn(bv,j,&z);
837:   return(0);
838: }

840: /*
841:   Compute Y^H*A*X: right part column by column (with MatMult) and bottom
842:   part row by row (with MatMultHermitianTranspose); result placed in marray[*,ldm]
843: */
844: PETSC_STATIC_INLINE PetscErrorCode BVMatProject_Vec(BV X,Mat A,BV Y,PetscScalar *marray,PetscInt ldm,PetscBool symm)
845: {
847:   PetscInt       i,j,lx,ly,kx,ky,ulim;
848:   Vec            z,f;

851:   lx = X->l; kx = X->k;
852:   ly = Y->l; ky = Y->k;
853:   BVCreateVec(X,&f);
854:   for (j=lx;j<kx;j++) {
855:     BVGetColumn(X,j,&z);
856:     MatMult(A,z,f);
857:     BVRestoreColumn(X,j,&z);
858:     ulim = PetscMin(ly+(j-lx)+1,ky);
859:     Y->l = 0; Y->k = ulim;
860:     (*Y->ops->dotvec)(Y,f,marray+j*ldm);
861:     if (symm) {
862:       for (i=0;i<j;i++) marray[j+i*ldm] = PetscConj(marray[i+j*ldm]);
863:     }
864:   }
865:   if (!symm) {
866:     BV_AllocateCoeffs(Y);
867:     for (j=ly;j<ky;j++) {
868:       BVGetColumn(Y,j,&z);
869:       MatMultHermitianTranspose(A,z,f);
870:       BVRestoreColumn(Y,j,&z);
871:       ulim = PetscMin(lx+(j-ly),kx);
872:       X->l = 0; X->k = ulim;
873:       (*X->ops->dotvec)(X,f,Y->h);
874:       for (i=0;i<ulim;i++) marray[j+i*ldm] = PetscConj(Y->h[i]);
875:     }
876:   }
877:   VecDestroy(&f);
878:   X->l = lx; X->k = kx;
879:   Y->l = ly; Y->k = ky;
880:   return(0);
881: }

883: /*
884:   Compute Y^H*A*X= [   --   | Y0'*W1 ]
885:                    [ Y1'*W0 | Y1'*W1 ]
886:   Allocates auxiliary BV to store the result of A*X, then one BVDot
887:   call for top-right part and another one for bottom part;
888:   result placed in marray[*,ldm]
889: */
890: PETSC_STATIC_INLINE PetscErrorCode BVMatProject_MatMult(BV X,Mat A,BV Y,PetscScalar *marray,PetscInt ldm)
891: {
893:   PetscInt       j,lx,ly,kx,ky;
894:   PetscScalar    *harray;
895:   Mat            H;
896:   BV             W;

899:   lx = X->l; kx = X->k;
900:   ly = Y->l; ky = Y->k;
901:   BVDuplicate(X,&W);
902:   X->l = 0; X->k = kx;
903:   W->l = 0; W->k = kx;
904:   BVMatMult(X,A,W);

906:   /* top-right part, Y0'*AX1 */
907:   if (ly>0 && lx<kx) {
908:     MatCreateSeqDense(PETSC_COMM_SELF,ly,kx,NULL,&H);
909:     W->l = lx; W->k = kx;
910:     Y->l = 0;  Y->k = ly;
911:     BVDot(W,Y,H);
912:     MatDenseGetArray(H,&harray);
913:     for (j=lx;j<kx;j++) {
914:       PetscMemcpy(marray+j*ldm,harray+j*ly,ly*sizeof(PetscScalar));
915:     }
916:     MatDenseRestoreArray(H,&harray);
917:     MatDestroy(&H);
918:   }

920:   /* bottom part, Y1'*AX */
921:   if (kx>0 && ly<ky) {
922:     MatCreateSeqDense(PETSC_COMM_SELF,ky,kx,NULL,&H);
923:     W->l = 0;  W->k = kx;
924:     Y->l = ly; Y->k = ky;
925:     BVDot(W,Y,H);
926:     MatDenseGetArray(H,&harray);
927:     for (j=0;j<kx;j++) {
928:       PetscMemcpy(marray+j*ldm+ly,harray+j*ky+ly,(ky-ly)*sizeof(PetscScalar));
929:     }
930:     MatDenseRestoreArray(H,&harray);
931:     MatDestroy(&H);
932:   }
933:   BVDestroy(&W);
934:   X->l = lx; X->k = kx;
935:   Y->l = ly; Y->k = ky;
936:   return(0);
937: }

939: /*
940:   Compute Y^H*A*X= [   --   | Y0'*W1 ]
941:                    [ Y1'*W0 | Y1'*W1 ]
942:   First stage: allocate auxiliary BV to store A*X1, one BVDot for right part;
943:   Second stage: resize BV to accomodate A'*Y1, then call BVDot for transpose of
944:   bottom-left part; result placed in marray[*,ldm]
945: */
946: PETSC_STATIC_INLINE PetscErrorCode BVMatProject_MatMult_2(BV X,Mat A,BV Y,PetscScalar *marray,PetscInt ldm,PetscBool symm)
947: {
949:   PetscInt       i,j,lx,ly,kx,ky;
950:   PetscScalar    *harray;
951:   Mat            H;
952:   BV             W;

955:   lx = X->l; kx = X->k;
956:   ly = Y->l; ky = Y->k;

958:   /* right part, Y'*AX1 */
959:   BVDuplicateResize(X,kx-lx,&W);
960:   if (ky>0 && lx<kx) {
961:     BVMatMult(X,A,W);
962:     MatCreateSeqDense(PETSC_COMM_SELF,ky,kx-lx,NULL,&H);
963:     Y->l = 0; Y->k = ky;
964:     BVDot(W,Y,H);
965:     MatDenseGetArray(H,&harray);
966:     for (j=lx;j<kx;j++) {
967:       PetscMemcpy(marray+j*ldm,harray+(j-lx)*ky,ky*sizeof(PetscScalar));
968:     }
969:     MatDenseRestoreArray(H,&harray);
970:     MatDestroy(&H);
971:   }

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

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

1017:   lx = X->l; kx = X->k;
1018:   ly = Y->l; ky = Y->k;

1020:   /* top-right part, Y0'*X1 */
1021:   if (ly>0 && lx<kx) {
1022:     MatCreateSeqDense(PETSC_COMM_SELF,ly,kx,NULL,&H);
1023:     X->l = lx; X->k = kx;
1024:     Y->l = 0;  Y->k = ly;
1025:     BVDot(X,Y,H);
1026:     MatDenseGetArray(H,&harray);
1027:     for (j=lx;j<kx;j++) {
1028:       PetscMemcpy(marray+j*ldm,harray+j*ly,ly*sizeof(PetscScalar));
1029:     }
1030:     MatDenseRestoreArray(H,&harray);
1031:     MatDestroy(&H);
1032:   }

1034:   /* bottom part, Y1'*X */
1035:   if (kx>0 && ly<ky) {
1036:     MatCreateSeqDense(PETSC_COMM_SELF,ky,kx,NULL,&H);
1037:     X->l = 0;  X->k = kx;
1038:     Y->l = ly; Y->k = ky;
1039:     BVDot(X,Y,H);
1040:     MatDenseGetArray(H,&harray);
1041:     for (j=0;j<kx;j++) {
1042:       PetscMemcpy(marray+j*ldm+ly,harray+j*ky+ly,(ky-ly)*sizeof(PetscScalar));
1043:     }
1044:     MatDenseRestoreArray(H,&harray);
1045:     MatDestroy(&H);
1046:   }
1047:   X->l = lx; X->k = kx;
1048:   Y->l = ly; Y->k = ky;
1049:   return(0);
1050: }

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

1055:    Collective on BV

1057:    Input Parameters:
1058: +  X - basis vectors
1059: .  A - (optional) matrix to be projected
1060: .  Y - left basis vectors, can be equal to X
1061: -  M - Mat object where the result must be placed

1063:    Output Parameter:
1064: .  M - the resulting matrix

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

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

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

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

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

1091:    Level: intermediate

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

1110:   BVCheckSizes(X,1);
1111:   if (A) {
1114:   }
1116:   BVCheckSizes(Y,3);
1119:   PetscObjectTypeCompare((PetscObject)M,MATSEQDENSE,&match);
1120:   if (!match) SETERRQ(PetscObjectComm((PetscObject)X),PETSC_ERR_SUP,"Matrix M must be of type seqdense");

1122:   MatGetSize(M,&m,&n);
1123:   if (m<Y->k) SETERRQ2(PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_SIZ,"Matrix M has %D rows, should have at least %D",m,Y->k);
1124:   if (n<X->k) SETERRQ2(PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_SIZ,"Matrix M has %D columns, should have at least %D",n,X->k);
1125:   if (X->n!=Y->n) SETERRQ2(PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension X %D, Y %D",X->n,Y->n);

1127:   PetscLogEventBegin(BV_MatProject,X,A,Y,0);
1128:   /* temporarily set standard inner product */
1129:   Xmatrix = X->matrix;
1130:   Ymatrix = Y->matrix;
1131:   X->matrix = Y->matrix = NULL;

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

1140:   MatDenseGetArray(M,&marray);

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

1160:   MatDenseRestoreArray(M,&marray);
1161:   PetscLogEventEnd(BV_MatProject,X,A,Y,0);
1162:   /* restore non-standard inner product */
1163:   X->matrix = Xmatrix;
1164:   Y->matrix = Ymatrix;
1165:   return(0);
1166: }