Actual source code: bvops.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, except those involving global communication
 12: */

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

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

 20:    Logically Collective

 22:    Input Parameters:
 23: +  Y     - first basis vectors context (modified on output)
 24: .  alpha - first scalar
 25: .  beta  - second scalar
 26: .  X     - second basis vectors context
 27: -  Q     - (optional) sequential dense matrix

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

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

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

 45:    Level: intermediate

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

 53:   PetscFunctionBegin;
 60:   BVCheckSizes(Y,1);
 61:   BVCheckOp(Y,1,mult);
 63:   BVCheckSizes(X,4);
 65:   PetscCheckSameTypeAndComm(Y,1,X,4);
 66:   PetscCheck(X!=Y,PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_WRONG,"X and Y arguments must be different");
 67:   if (Q) {
 68:     SlepcMatCheckSeq(Q);
 69:     PetscCall(MatGetSize(Q,&m,&n));
 70:     PetscCheck(m>=X->k,PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_SIZ,"Mat argument has %" PetscInt_FMT " rows, should have at least %" PetscInt_FMT,m,X->k);
 71:     PetscCheck(n>=Y->k,PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_SIZ,"Mat argument has %" PetscInt_FMT " columns, should have at least %" PetscInt_FMT,n,Y->k);
 72:   }
 73:   PetscCheck(X->n==Y->n,PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension X %" PetscInt_FMT ", Y %" PetscInt_FMT,X->n,Y->n);

 75:   PetscCall(PetscLogEventBegin(BV_Mult,X,Y,0,0));
 76:   PetscUseTypeMethod(Y,mult,alpha,beta,X,Q);
 77:   PetscCall(PetscLogEventEnd(BV_Mult,X,Y,0,0));
 78:   PetscCall(PetscObjectStateIncrease((PetscObject)Y));
 79:   PetscFunctionReturn(PETSC_SUCCESS);
 80: }

 82: /*@
 83:    BVMultVec - Computes y = beta*y + alpha*X*q.

 85:    Logically Collective

 87:    Input Parameters:
 88: +  X     - a basis vectors object
 89: .  alpha - first scalar
 90: .  beta  - second scalar
 91: .  y     - a vector (modified on output)
 92: -  q     - an array of scalars

 94:    Notes:
 95:    This operation is the analogue of BVMult() but with a BV and a Vec,
 96:    instead of two BV. Note that arguments are listed in different order
 97:    with respect to BVMult().

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

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

106:    Level: intermediate

108: .seealso: BVMult(), BVMultColumn(), BVMultInPlace(), BVSetActiveColumns()
109: @*/
110: PetscErrorCode BVMultVec(BV X,PetscScalar alpha,PetscScalar beta,Vec y,PetscScalar q[])
111: {
112:   PetscInt       n,N;

114:   PetscFunctionBegin;
119:   PetscAssertPointer(q,5);
121:   BVCheckSizes(X,1);
122:   BVCheckOp(X,1,multvec);
124:   PetscCheckSameComm(X,1,y,4);

126:   PetscCall(VecGetSize(y,&N));
127:   PetscCall(VecGetLocalSize(y,&n));
128:   PetscCheck(N==X->N && n==X->n,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_INCOMP,"Vec sizes (global %" PetscInt_FMT ", local %" PetscInt_FMT ") do not match BV sizes (global %" PetscInt_FMT ", local %" PetscInt_FMT ")",N,n,X->N,X->n);

130:   PetscCall(PetscLogEventBegin(BV_MultVec,X,y,0,0));
131:   PetscUseTypeMethod(X,multvec,alpha,beta,y,q);
132:   PetscCall(PetscLogEventEnd(BV_MultVec,X,y,0,0));
133:   PetscFunctionReturn(PETSC_SUCCESS);
134: }

136: /*@
137:    BVMultColumn - Computes y = beta*y + alpha*X*q, where y is the j-th column
138:    of X.

140:    Logically Collective

142:    Input Parameters:
143: +  X     - a basis vectors object
144: .  alpha - first scalar
145: .  beta  - second scalar
146: .  j     - the column index
147: -  q     - an array of scalars

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

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

161:    Level: advanced

163: .seealso: BVMult(), BVMultVec(), BVMultInPlace(), BVSetActiveColumns()
164: @*/
165: PetscErrorCode BVMultColumn(BV X,PetscScalar alpha,PetscScalar beta,PetscInt j,PetscScalar *q)
166: {
167:   PetscInt       ksave;
168:   Vec            y;

170:   PetscFunctionBegin;
176:   BVCheckSizes(X,1);

178:   PetscCheck(j>=0,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
179:   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);

181:   PetscCall(PetscLogEventBegin(BV_MultVec,X,0,0,0));
182:   ksave = X->k;
183:   X->k = j;
184:   if (!q && !X->buffer) PetscCall(BVGetBufferVec(X,&X->buffer));
185:   PetscCall(BVGetColumn(X,j,&y));
186:   PetscUseTypeMethod(X,multvec,alpha,beta,y,q);
187:   PetscCall(BVRestoreColumn(X,j,&y));
188:   X->k = ksave;
189:   PetscCall(PetscLogEventEnd(BV_MultVec,X,0,0,0));
190:   PetscCall(PetscObjectStateIncrease((PetscObject)X));
191:   PetscFunctionReturn(PETSC_SUCCESS);
192: }

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

197:    Logically Collective

199:    Input Parameters:
200: +  Q - a sequential dense matrix
201: .  s - first column of V to be overwritten
202: -  e - first column of V not to be overwritten

204:    Input/Output Parameter:
205: .  V - basis vectors

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

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

216:    Level: intermediate

218: .seealso: BVMult(), BVMultVec(), BVMultInPlaceHermitianTranspose(), BVSetActiveColumns()
219: @*/
220: PetscErrorCode BVMultInPlace(BV V,Mat Q,PetscInt s,PetscInt e)
221: {
222:   PetscInt       m,n;

224:   PetscFunctionBegin;
230:   BVCheckSizes(V,1);
232:   SlepcMatCheckSeq(Q);

234:   PetscCheck(s>=V->l && s<=V->m,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Argument s has wrong value %" PetscInt_FMT ", should be between %" PetscInt_FMT " and %" PetscInt_FMT,s,V->l,V->m);
235:   PetscCheck(e>=V->l && e<=V->m,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Argument e has wrong value %" PetscInt_FMT ", should be between %" PetscInt_FMT " and %" PetscInt_FMT,e,V->l,V->m);
236:   PetscCall(MatGetSize(Q,&m,&n));
237:   PetscCheck(m>=V->k,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"Mat argument has %" PetscInt_FMT " rows, should have at least %" PetscInt_FMT,m,V->k);
238:   PetscCheck(e<=n,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"Mat argument only has %" PetscInt_FMT " columns, the requested value of e is larger: %" PetscInt_FMT,n,e);

240:   PetscCall(PetscLogEventBegin(BV_MultInPlace,V,Q,0,0));
241:   PetscUseTypeMethod(V,multinplace,Q,s,e);
242:   PetscCall(PetscLogEventEnd(BV_MultInPlace,V,Q,0,0));
243:   PetscCall(PetscObjectStateIncrease((PetscObject)V));
244:   PetscFunctionReturn(PETSC_SUCCESS);
245: }

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

250:    Logically Collective

252:    Input Parameters:
253: +  Q - a sequential dense matrix
254: .  s - first column of V to be overwritten
255: -  e - first column of V not to be overwritten

257:    Input/Output Parameter:
258: .  V - basis vectors

260:    Notes:
261:    This is a variant of BVMultInPlace() where the conjugate transpose
262:    of Q is used.

264:    Level: intermediate

266: .seealso: BVMultInPlace()
267: @*/
268: PetscErrorCode BVMultInPlaceHermitianTranspose(BV V,Mat Q,PetscInt s,PetscInt e)
269: {
270:   PetscInt       m,n;

272:   PetscFunctionBegin;
278:   BVCheckSizes(V,1);
280:   SlepcMatCheckSeq(Q);

282:   PetscCheck(s>=V->l && s<=V->m,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Argument s has wrong value %" PetscInt_FMT ", should be between %" PetscInt_FMT " and %" PetscInt_FMT,s,V->l,V->m);
283:   PetscCheck(e>=V->l && e<=V->m,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Argument e has wrong value %" PetscInt_FMT ", should be between %" PetscInt_FMT " and %" PetscInt_FMT,e,V->l,V->m);
284:   PetscCall(MatGetSize(Q,&m,&n));
285:   PetscCheck(n>=V->k,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"Mat argument has %" PetscInt_FMT " columns, should have at least %" PetscInt_FMT,n,V->k);
286:   PetscCheck(e<=m,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"Mat argument only has %" PetscInt_FMT " rows, the requested value of e is larger: %" PetscInt_FMT,m,e);

288:   PetscCall(PetscLogEventBegin(BV_MultInPlace,V,Q,0,0));
289:   PetscUseTypeMethod(V,multinplacetrans,Q,s,e);
290:   PetscCall(PetscLogEventEnd(BV_MultInPlace,V,Q,0,0));
291:   PetscCall(PetscObjectStateIncrease((PetscObject)V));
292:   PetscFunctionReturn(PETSC_SUCCESS);
293: }

295: /*@
296:    BVScale - Multiply the BV entries by a scalar value.

298:    Logically Collective

300:    Input Parameters:
301: +  bv    - basis vectors
302: -  alpha - scaling factor

304:    Note:
305:    All active columns (except the leading ones) are scaled.

307:    Level: intermediate

309: .seealso: BVScaleColumn(), BVSetActiveColumns()
310: @*/
311: PetscErrorCode BVScale(BV bv,PetscScalar alpha)
312: {
313:   PetscFunctionBegin;
317:   BVCheckSizes(bv,1);
318:   if (alpha == (PetscScalar)1.0) PetscFunctionReturn(PETSC_SUCCESS);

320:   PetscCall(PetscLogEventBegin(BV_Scale,bv,0,0,0));
321:   PetscUseTypeMethod(bv,scale,-1,alpha);
322:   PetscCall(PetscLogEventEnd(BV_Scale,bv,0,0,0));
323:   PetscCall(PetscObjectStateIncrease((PetscObject)bv));
324:   PetscFunctionReturn(PETSC_SUCCESS);
325: }

327: /*@
328:    BVScaleColumn - Scale one column of a BV.

330:    Logically Collective

332:    Input Parameters:
333: +  bv    - basis vectors
334: .  j     - column number to be scaled
335: -  alpha - scaling factor

337:    Level: intermediate

339: .seealso: BVScale(), BVSetActiveColumns()
340: @*/
341: PetscErrorCode BVScaleColumn(BV bv,PetscInt j,PetscScalar alpha)
342: {
343:   PetscFunctionBegin;
348:   BVCheckSizes(bv,1);

350:   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);
351:   if (alpha == (PetscScalar)1.0) PetscFunctionReturn(PETSC_SUCCESS);

353:   PetscCall(PetscLogEventBegin(BV_Scale,bv,0,0,0));
354:   PetscUseTypeMethod(bv,scale,j,alpha);
355:   PetscCall(PetscLogEventEnd(BV_Scale,bv,0,0,0));
356:   PetscCall(PetscObjectStateIncrease((PetscObject)bv));
357:   PetscFunctionReturn(PETSC_SUCCESS);
358: }

360: static inline PetscErrorCode BVSetRandomColumn_Private(BV bv,PetscInt k)
361: {
362:   PetscInt       i,low,high;
363:   PetscScalar    *px,t;
364:   Vec            x;

366:   PetscFunctionBegin;
367:   PetscCall(BVGetColumn(bv,k,&x));
368:   if (bv->rrandom) {  /* generate the same vector irrespective of number of processes */
369:     PetscCall(VecGetOwnershipRange(x,&low,&high));
370:     PetscCall(VecGetArray(x,&px));
371:     for (i=0;i<bv->N;i++) {
372:       PetscCall(PetscRandomGetValue(bv->rand,&t));
373:       if (i>=low && i<high) px[i-low] = t;
374:     }
375:     PetscCall(VecRestoreArray(x,&px));
376:   } else PetscCall(VecSetRandom(x,bv->rand));
377:   PetscCall(BVRestoreColumn(bv,k,&x));
378:   PetscFunctionReturn(PETSC_SUCCESS);
379: }

381: static inline PetscErrorCode BVSetRandomNormalColumn_Private(BV bv,PetscInt k,Vec w1,Vec w2)
382: {
383:   PetscInt       i,low,high;
384:   PetscScalar    *px,s,t;
385:   Vec            x;

387:   PetscFunctionBegin;
388:   PetscCall(BVGetColumn(bv,k,&x));
389:   if (bv->rrandom) {  /* generate the same vector irrespective of number of processes */
390:     PetscCall(VecGetOwnershipRange(x,&low,&high));
391:     PetscCall(VecGetArray(x,&px));
392:     for (i=0;i<bv->N;i++) {
393:       PetscCall(PetscRandomGetValue(bv->rand,&s));
394:       PetscCall(PetscRandomGetValue(bv->rand,&t));
395:       if (i>=low && i<high) {
396: #if defined(PETSC_USE_COMPLEX)
397:         px[i-low] = PetscCMPLX(PetscSqrtReal(-2.0*PetscLogReal(PetscRealPart(s)))*PetscCosReal(2.0*PETSC_PI*PetscRealPart(t)),PetscSqrtReal(-2.0*PetscLogReal(PetscImaginaryPart(s)))*PetscCosReal(2.0*PETSC_PI*PetscImaginaryPart(t)));
398: #else
399:         px[i-low] = PetscSqrtReal(-2.0*PetscLogReal(s))*PetscCosReal(2.0*PETSC_PI*t);
400: #endif
401:       }
402:     }
403:     PetscCall(VecRestoreArray(x,&px));
404:   } else PetscCall(VecSetRandomNormal(x,bv->rand,w1,w2));
405:   PetscCall(BVRestoreColumn(bv,k,&x));
406:   PetscFunctionReturn(PETSC_SUCCESS);
407: }

409: static inline PetscErrorCode BVSetRandomSignColumn_Private(BV bv,PetscInt k)
410: {
411:   PetscInt       i,low,high;
412:   PetscScalar    *px,t;
413:   Vec            x;

415:   PetscFunctionBegin;
416:   PetscCall(BVGetColumn(bv,k,&x));
417:   PetscCall(VecGetOwnershipRange(x,&low,&high));
418:   if (bv->rrandom) {  /* generate the same vector irrespective of number of processes */
419:     PetscCall(VecGetArray(x,&px));
420:     for (i=0;i<bv->N;i++) {
421:       PetscCall(PetscRandomGetValue(bv->rand,&t));
422:       if (i>=low && i<high) px[i-low] = (PetscRealPart(t)<0.5)? -1.0: 1.0;
423:     }
424:     PetscCall(VecRestoreArray(x,&px));
425:   } else {
426:     PetscCall(VecSetRandom(x,bv->rand));
427:     PetscCall(VecGetArray(x,&px));
428:     for (i=low;i<high;i++) {
429:       px[i-low] = (PetscRealPart(px[i-low])<0.5)? -1.0: 1.0;
430:     }
431:     PetscCall(VecRestoreArray(x,&px));
432:   }
433:   PetscCall(BVRestoreColumn(bv,k,&x));
434:   PetscFunctionReturn(PETSC_SUCCESS);
435: }

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

440:    Logically Collective

442:    Input Parameters:
443: .  bv - basis vectors

445:    Note:
446:    All active columns (except the leading ones) are modified.

448:    Level: advanced

450: .seealso: BVSetRandomContext(), BVSetRandomColumn(), BVSetRandomNormal(), BVSetRandomSign(), BVSetRandomCond(), BVSetActiveColumns()
451: @*/
452: PetscErrorCode BVSetRandom(BV bv)
453: {
454:   PetscInt       k;

456:   PetscFunctionBegin;
459:   BVCheckSizes(bv,1);

461:   PetscCall(BVGetRandomContext(bv,&bv->rand));
462:   PetscCall(PetscLogEventBegin(BV_SetRandom,bv,0,0,0));
463:   for (k=bv->l;k<bv->k;k++) PetscCall(BVSetRandomColumn_Private(bv,k));
464:   PetscCall(PetscLogEventEnd(BV_SetRandom,bv,0,0,0));
465:   PetscCall(PetscObjectStateIncrease((PetscObject)bv));
466:   PetscFunctionReturn(PETSC_SUCCESS);
467: }

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

472:    Logically Collective

474:    Input Parameters:
475: +  bv - basis vectors
476: -  j  - column number to be set

478:    Level: advanced

480: .seealso: BVSetRandomContext(), BVSetRandom(), BVSetRandomNormal(), BVSetRandomCond()
481: @*/
482: PetscErrorCode BVSetRandomColumn(BV bv,PetscInt j)
483: {
484:   PetscFunctionBegin;
488:   BVCheckSizes(bv,1);
489:   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);

491:   PetscCall(BVGetRandomContext(bv,&bv->rand));
492:   PetscCall(PetscLogEventBegin(BV_SetRandom,bv,0,0,0));
493:   PetscCall(BVSetRandomColumn_Private(bv,j));
494:   PetscCall(PetscLogEventEnd(BV_SetRandom,bv,0,0,0));
495:   PetscCall(PetscObjectStateIncrease((PetscObject)bv));
496:   PetscFunctionReturn(PETSC_SUCCESS);
497: }

499: /*@
500:    BVSetRandomNormal - Set the columns of a BV to random numbers with a normal
501:    distribution.

503:    Logically Collective

505:    Input Parameter:
506: .  bv - basis vectors

508:    Notes:
509:    All active columns (except the leading ones) are modified.

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

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

519:    Level: advanced

521: .seealso: BVSetRandomContext(), BVSetRandom(), BVSetRandomColumn(), BVSetRandomCond(), BVSetActiveColumns()
522: @*/
523: PetscErrorCode BVSetRandomNormal(BV bv)
524: {
525:   PetscInt       k;
526:   Vec            w1=NULL,w2=NULL;

528:   PetscFunctionBegin;
531:   BVCheckSizes(bv,1);

533:   PetscCall(BVGetRandomContext(bv,&bv->rand));
534:   if (!bv->rrandom) {
535:     PetscCall(BVCreateVec(bv,&w1));
536:     PetscCall(BVCreateVec(bv,&w2));
537:   }
538:   PetscCall(PetscLogEventBegin(BV_SetRandom,bv,0,0,0));
539:   for (k=bv->l;k<bv->k;k++) PetscCall(BVSetRandomNormalColumn_Private(bv,k,w1,w2));
540:   PetscCall(PetscLogEventEnd(BV_SetRandom,bv,0,0,0));
541:   if (!bv->rrandom) {
542:     PetscCall(VecDestroy(&w1));
543:     PetscCall(VecDestroy(&w2));
544:   }
545:   PetscCall(PetscObjectStateIncrease((PetscObject)bv));
546:   PetscFunctionReturn(PETSC_SUCCESS);
547: }

549: /*@
550:    BVSetRandomSign - Set the entries of a BV to values 1 or -1 with equal
551:    probability.

553:    Logically Collective

555:    Input Parameter:
556: .  bv - basis vectors

558:    Notes:
559:    All active columns (except the leading ones) are modified.

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

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

569:    Level: advanced

571: .seealso: BVSetRandomContext(), BVSetRandom(), BVSetRandomColumn(), BVSetActiveColumns()
572: @*/
573: PetscErrorCode BVSetRandomSign(BV bv)
574: {
575:   PetscScalar    low,high;
576:   PetscInt       k;

578:   PetscFunctionBegin;
581:   BVCheckSizes(bv,1);

583:   PetscCall(BVGetRandomContext(bv,&bv->rand));
584:   PetscCall(PetscRandomGetInterval(bv->rand,&low,&high));
585:   PetscCheck(PetscRealPart(low)==0.0 && PetscRealPart(high)==1.0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"The PetscRandom object in the BV must have interval [0,1]");
586:   PetscCall(PetscLogEventBegin(BV_SetRandom,bv,0,0,0));
587:   for (k=bv->l;k<bv->k;k++) PetscCall(BVSetRandomSignColumn_Private(bv,k));
588:   PetscCall(PetscLogEventEnd(BV_SetRandom,bv,0,0,0));
589:   PetscCall(PetscObjectStateIncrease((PetscObject)bv));
590:   PetscFunctionReturn(PETSC_SUCCESS);
591: }

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

597:    Logically Collective

599:    Input Parameters:
600: +  bv    - basis vectors
601: -  condn - condition number

603:    Note:
604:    All active columns (except the leading ones) are modified.

606:    Level: advanced

608: .seealso: BVSetRandomContext(), BVSetRandom(), BVSetRandomColumn(), BVSetRandomNormal(), BVSetActiveColumns()
609: @*/
610: PetscErrorCode BVSetRandomCond(BV bv,PetscReal condn)
611: {
612:   PetscInt       k,i;
613:   PetscScalar    *eig,*d;
614:   DS             ds;
615:   Mat            A,X,Xt,M,G;

617:   PetscFunctionBegin;
620:   BVCheckSizes(bv,1);

622:   PetscCall(BVGetRandomContext(bv,&bv->rand));
623:   PetscCall(PetscLogEventBegin(BV_SetRandom,bv,0,0,0));
624:   /* B = rand(n,k) */
625:   for (k=bv->l;k<bv->k;k++) PetscCall(BVSetRandomColumn_Private(bv,k));
626:   PetscCall(DSCreate(PetscObjectComm((PetscObject)bv),&ds));
627:   PetscCall(DSSetType(ds,DSHEP));
628:   PetscCall(DSAllocate(ds,bv->m));
629:   PetscCall(DSSetDimensions(ds,bv->k,bv->l,bv->k));
630:   /* [V,S] = eig(B'*B) */
631:   PetscCall(DSGetMat(ds,DS_MAT_A,&A));
632:   PetscCall(BVDot(bv,bv,A));
633:   PetscCall(DSRestoreMat(ds,DS_MAT_A,&A));
634:   PetscCall(PetscMalloc1(bv->m,&eig));
635:   PetscCall(DSSolve(ds,eig,NULL));
636:   PetscCall(DSSynchronize(ds,eig,NULL));
637:   PetscCall(DSVectors(ds,DS_MAT_X,NULL,NULL));
638:   /* M = diag(linspace(1/condn,1,n)./sqrt(diag(S)))' */
639:   PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,bv->k,bv->k,NULL,&M));
640:   PetscCall(MatZeroEntries(M));
641:   PetscCall(MatDenseGetArray(M,&d));
642:   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]);
643:   PetscCall(MatDenseRestoreArray(M,&d));
644:   /* G = X*M*X' */
645:   PetscCall(DSGetMat(ds,DS_MAT_X,&X));
646:   PetscCall(MatTranspose(X,MAT_INITIAL_MATRIX,&Xt));
647:   PetscCall(MatProductCreate(Xt,M,NULL,&G));
648:   PetscCall(MatProductSetType(G,MATPRODUCT_PtAP));
649:   PetscCall(MatProductSetFromOptions(G));
650:   PetscCall(MatProductSymbolic(G));
651:   PetscCall(MatProductNumeric(G));
652:   PetscCall(MatProductClear(G));
653:   PetscCall(DSRestoreMat(ds,DS_MAT_X,&X));
654:   PetscCall(MatDestroy(&Xt));
655:   PetscCall(MatDestroy(&M));
656:   /* B = B*G */
657:   PetscCall(BVMultInPlace(bv,G,bv->l,bv->k));
658:   PetscCall(MatDestroy(&G));
659:   PetscCall(PetscFree(eig));
660:   PetscCall(DSDestroy(&ds));
661:   PetscCall(PetscLogEventEnd(BV_SetRandom,bv,0,0,0));
662:   PetscCall(PetscObjectStateIncrease((PetscObject)bv));
663:   PetscFunctionReturn(PETSC_SUCCESS);
664: }

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

669:    Neighbor-wise Collective

671:    Input Parameters:
672: +  V - basis vectors context
673: -  A - the matrix

675:    Output Parameter:
676: .  Y - the result

678:    Notes:
679:    Both V and Y must be distributed in the same manner. Only active columns
680:    (excluding the leading ones) are processed.
681:    In the result Y, columns are overwritten starting from the leading ones.
682:    The number of active columns in V and Y should match, although they need
683:    not be the same columns.

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

688:    Level: beginner

690: .seealso: BVCopy(), BVSetActiveColumns(), BVMatMultColumn(), BVMatMultTranspose(), BVMatMultHermitianTranspose(), BVSetMatMultMethod()
691: @*/
692: PetscErrorCode BVMatMult(BV V,Mat A,BV Y)
693: {
694:   PetscInt       M,N,m,n;

696:   PetscFunctionBegin;
699:   BVCheckSizes(V,1);
700:   BVCheckOp(V,1,matmult);
705:   BVCheckSizes(Y,3);
706:   PetscCheckSameComm(V,1,A,2);
707:   PetscCheckSameTypeAndComm(V,1,Y,3);

709:   PetscCall(MatGetSize(A,&M,&N));
710:   PetscCall(MatGetLocalSize(A,&m,&n));
711:   PetscCheck(M==Y->N,PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Mismatching row dimension A %" PetscInt_FMT ", Y %" PetscInt_FMT,M,Y->N);
712:   PetscCheck(m==Y->n,PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Mismatching local row dimension A %" PetscInt_FMT ", Y %" PetscInt_FMT,m,Y->n);
713:   PetscCheck(N==V->N,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching column dimension A %" PetscInt_FMT ", V %" PetscInt_FMT,N,V->N);
714:   PetscCheck(n==V->n,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching local column dimension A %" PetscInt_FMT ", V %" PetscInt_FMT,n,V->n);
715:   PetscCheck(V->k-V->l==Y->k-Y->l,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"Y has %" PetscInt_FMT " active columns, should match %" PetscInt_FMT " active columns in V",Y->k-Y->l,V->k-V->l);

717:   PetscCall(PetscLogEventBegin(BV_MatMult,V,A,Y,0));
718:   PetscUseTypeMethod(V,matmult,A,Y);
719:   PetscCall(PetscLogEventEnd(BV_MatMult,V,A,Y,0));
720:   PetscCall(PetscObjectStateIncrease((PetscObject)Y));
721:   PetscFunctionReturn(PETSC_SUCCESS);
722: }

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

728:    Neighbor-wise Collective

730:    Input Parameters:
731: +  V - basis vectors context
732: -  A - the matrix

734:    Output Parameter:
735: .  Y - the result

737:    Notes:
738:    Both V and Y must be distributed in the same manner. Only active columns
739:    (excluding the leading ones) are processed.
740:    In the result Y, columns are overwritten starting from the leading ones.
741:    The number of active columns in V and Y should match, although they need
742:    not be the same columns.

744:    Currently implemented via MatCreateTranspose().

746:    Level: beginner

748: .seealso: BVMatMult(), BVMatMultHermitianTranspose()
749: @*/
750: PetscErrorCode BVMatMultTranspose(BV V,Mat A,BV Y)
751: {
752:   PetscInt       M,N,m,n;
753:   Mat            AT;

755:   PetscFunctionBegin;
758:   BVCheckSizes(V,1);
763:   BVCheckSizes(Y,3);
764:   PetscCheckSameComm(V,1,A,2);
765:   PetscCheckSameTypeAndComm(V,1,Y,3);

767:   PetscCall(MatGetSize(A,&M,&N));
768:   PetscCall(MatGetLocalSize(A,&m,&n));
769:   PetscCheck(M==V->N,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching row dimension A %" PetscInt_FMT ", V %" PetscInt_FMT,M,V->N);
770:   PetscCheck(m==V->n,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching local row dimension A %" PetscInt_FMT ", V %" PetscInt_FMT,m,V->n);
771:   PetscCheck(N==Y->N,PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Mismatching column dimension A %" PetscInt_FMT ", Y %" PetscInt_FMT,N,Y->N);
772:   PetscCheck(n==Y->n,PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Mismatching local column dimension A %" PetscInt_FMT ", Y %" PetscInt_FMT,n,Y->n);
773:   PetscCheck(V->k-V->l==Y->k-Y->l,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"Y has %" PetscInt_FMT " active columns, should match %" PetscInt_FMT " active columns in V",Y->k-Y->l,V->k-V->l);

775:   PetscCall(MatCreateTranspose(A,&AT));
776:   PetscCall(BVMatMult(V,AT,Y));
777:   PetscCall(MatDestroy(&AT));
778:   PetscFunctionReturn(PETSC_SUCCESS);
779: }

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

785:    Neighbor-wise Collective

787:    Input Parameters:
788: +  V - basis vectors context
789: -  A - the matrix

791:    Output Parameter:
792: .  Y - the result

794:    Note:
795:    Both V and Y must be distributed in the same manner. Only active columns
796:    (excluding the leading ones) are processed.
797:    In the result Y, columns are overwritten starting from the leading ones.
798:    The number of active columns in V and Y should match, although they need
799:    not be the same columns.

801:    Currently implemented via MatCreateHermitianTranspose().

803:    Level: beginner

805: .seealso: BVMatMult(), BVMatMultTranspose()
806: @*/
807: PetscErrorCode BVMatMultHermitianTranspose(BV V,Mat A,BV Y)
808: {
809:   PetscInt       j,M,N,m,n;
810:   Vec            v,y;

812:   PetscFunctionBegin;
815:   BVCheckSizes(V,1);
820:   BVCheckSizes(Y,3);
821:   PetscCheckSameComm(V,1,A,2);
822:   PetscCheckSameTypeAndComm(V,1,Y,3);

824:   PetscCall(MatGetSize(A,&M,&N));
825:   PetscCall(MatGetLocalSize(A,&m,&n));
826:   PetscCheck(M==V->N,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching row dimension A %" PetscInt_FMT ", V %" PetscInt_FMT,M,V->N);
827:   PetscCheck(m==V->n,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching local row dimension A %" PetscInt_FMT ", V %" PetscInt_FMT,m,V->n);
828:   PetscCheck(N==Y->N,PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Mismatching column dimension A %" PetscInt_FMT ", Y %" PetscInt_FMT,N,Y->N);
829:   PetscCheck(n==Y->n,PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Mismatching local column dimension A %" PetscInt_FMT ", Y %" PetscInt_FMT,n,Y->n);
830:   PetscCheck(V->k-V->l==Y->k-Y->l,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"Y has %" PetscInt_FMT " active columns, should match %" PetscInt_FMT " active columns in V",Y->k-Y->l,V->k-V->l);

832:   /* TODO: recover this code if PETSc ever gets MATPRODUCT_AhB done
833:   PetscCall(MatCreateHermitianTranspose(A,&AH));
834:   PetscCall(BVMatMult(V,AH,Y));
835:   PetscCall(MatDestroy(&AH));
836:   */
837:   for (j=0;j<V->k-V->l;j++) {
838:     PetscCall(BVGetColumn(V,V->l+j,&v));
839:     PetscCall(BVGetColumn(Y,Y->l+j,&y));
840:     PetscCall(MatMultHermitianTranspose(A,v,y));
841:     PetscCall(BVRestoreColumn(V,V->l+j,&v));
842:     PetscCall(BVRestoreColumn(Y,Y->l+j,&y));
843:   }
844:   PetscFunctionReturn(PETSC_SUCCESS);
845: }

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

851:    Neighbor-wise Collective

853:    Input Parameters:
854: +  V - basis vectors context
855: .  A - the matrix
856: -  j - the column

858:    Level: beginner

860: .seealso: BVMatMult(), BVMatMultTransposeColumn(), BVMatMultHermitianTransposeColumn()
861: @*/
862: PetscErrorCode BVMatMultColumn(BV V,Mat A,PetscInt j)
863: {
864:   Vec            vj,vj1;

866:   PetscFunctionBegin;
869:   BVCheckSizes(V,1);
871:   PetscCheckSameComm(V,1,A,2);
873:   PetscCheck(j>=0,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
874:   PetscCheck(j+1<V->m,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Result should go in index j+1=%" PetscInt_FMT " but BV only has %" PetscInt_FMT " columns",j+1,V->m);

876:   PetscCall(PetscLogEventBegin(BV_MatMultVec,V,A,0,0));
877:   PetscCall(BVGetColumn(V,j,&vj));
878:   PetscCall(BVGetColumn(V,j+1,&vj1));
879:   PetscCall(MatMult(A,vj,vj1));
880:   PetscCall(BVRestoreColumn(V,j,&vj));
881:   PetscCall(BVRestoreColumn(V,j+1,&vj1));
882:   PetscCall(PetscLogEventEnd(BV_MatMultVec,V,A,0,0));
883:   PetscCall(PetscObjectStateIncrease((PetscObject)V));
884:   PetscFunctionReturn(PETSC_SUCCESS);
885: }

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

891:    Neighbor-wise Collective

893:    Input Parameters:
894: +  V - basis vectors context
895: .  A - the matrix
896: -  j - the column

898:    Level: beginner

900: .seealso: BVMatMult(), BVMatMultColumn(), BVMatMultHermitianTransposeColumn()
901: @*/
902: PetscErrorCode BVMatMultTransposeColumn(BV V,Mat A,PetscInt j)
903: {
904:   Vec            vj,vj1;

906:   PetscFunctionBegin;
909:   BVCheckSizes(V,1);
911:   PetscCheckSameComm(V,1,A,2);
913:   PetscCheck(j>=0,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
914:   PetscCheck(j+1<V->m,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Result should go in index j+1=%" PetscInt_FMT " but BV only has %" PetscInt_FMT " columns",j+1,V->m);

916:   PetscCall(PetscLogEventBegin(BV_MatMultVec,V,A,0,0));
917:   PetscCall(BVGetColumn(V,j,&vj));
918:   PetscCall(BVGetColumn(V,j+1,&vj1));
919:   PetscCall(MatMultTranspose(A,vj,vj1));
920:   PetscCall(BVRestoreColumn(V,j,&vj));
921:   PetscCall(BVRestoreColumn(V,j+1,&vj1));
922:   PetscCall(PetscLogEventEnd(BV_MatMultVec,V,A,0,0));
923:   PetscCall(PetscObjectStateIncrease((PetscObject)V));
924:   PetscFunctionReturn(PETSC_SUCCESS);
925: }

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

931:    Neighbor-wise Collective

933:    Input Parameters:
934: +  V - basis vectors context
935: .  A - the matrix
936: -  j - the column

938:    Level: beginner

940: .seealso: BVMatMult(), BVMatMultColumn(), BVMatMultTransposeColumn()
941: @*/
942: PetscErrorCode BVMatMultHermitianTransposeColumn(BV V,Mat A,PetscInt j)
943: {
944:   Vec            vj,vj1;

946:   PetscFunctionBegin;
949:   BVCheckSizes(V,1);
951:   PetscCheckSameComm(V,1,A,2);
953:   PetscCheck(j>=0,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
954:   PetscCheck(j+1<V->m,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Result should go in index j+1=%" PetscInt_FMT " but BV only has %" PetscInt_FMT " columns",j+1,V->m);

956:   PetscCall(PetscLogEventBegin(BV_MatMultVec,V,A,0,0));
957:   PetscCall(BVGetColumn(V,j,&vj));
958:   PetscCall(BVGetColumn(V,j+1,&vj1));
959:   PetscCall(MatMultHermitianTranspose(A,vj,vj1));
960:   PetscCall(BVRestoreColumn(V,j,&vj));
961:   PetscCall(BVRestoreColumn(V,j+1,&vj1));
962:   PetscCall(PetscLogEventEnd(BV_MatMultVec,V,A,0,0));
963:   PetscCall(PetscObjectStateIncrease((PetscObject)V));
964:   PetscFunctionReturn(PETSC_SUCCESS);
965: }