Actual source code: bvops.c

slepc-3.19.2 2023-09-05
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:     PetscCheckTypeNames(Q,MATSEQDENSE,MATSEQDENSECUDA);
 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;
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:   PetscCheckTypeNames(Q,MATSEQDENSE,MATSEQDENSECUDA);

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);
239:   if (s>=e) PetscFunctionReturn(PETSC_SUCCESS);

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

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

251:    Logically Collective

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

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

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

265:    Level: intermediate

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

273:   PetscFunctionBegin;
279:   BVCheckSizes(V,1);
281:   PetscCheckTypeNames(Q,MATSEQDENSE,MATSEQDENSECUDA);

283:   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);
284:   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);
285:   PetscCall(MatGetSize(Q,&m,&n));
286:   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);
287:   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:   if (s>=e || !V->n) PetscFunctionReturn(PETSC_SUCCESS);

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

297: /*@
298:    BVScale - Multiply the BV entries by a scalar value.

300:    Logically Collective

302:    Input Parameters:
303: +  bv    - basis vectors
304: -  alpha - scaling factor

306:    Note:
307:    All active columns (except the leading ones) are scaled.

309:    Level: intermediate

311: .seealso: BVScaleColumn(), BVSetActiveColumns()
312: @*/
313: PetscErrorCode BVScale(BV bv,PetscScalar alpha)
314: {
315:   PetscFunctionBegin;
319:   BVCheckSizes(bv,1);
320:   if (alpha == (PetscScalar)1.0) PetscFunctionReturn(PETSC_SUCCESS);

322:   PetscCall(PetscLogEventBegin(BV_Scale,bv,0,0,0));
323:   if (bv->n) PetscUseTypeMethod(bv,scale,-1,alpha);
324:   PetscCall(PetscLogEventEnd(BV_Scale,bv,0,0,0));
325:   PetscCall(PetscObjectStateIncrease((PetscObject)bv));
326:   PetscFunctionReturn(PETSC_SUCCESS);
327: }

329: /*@
330:    BVScaleColumn - Scale one column of a BV.

332:    Logically Collective

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

339:    Level: intermediate

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

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

355:   PetscCall(PetscLogEventBegin(BV_Scale,bv,0,0,0));
356:   if (bv->n) PetscUseTypeMethod(bv,scale,j,alpha);
357:   PetscCall(PetscLogEventEnd(BV_Scale,bv,0,0,0));
358:   PetscCall(PetscObjectStateIncrease((PetscObject)bv));
359:   PetscFunctionReturn(PETSC_SUCCESS);
360: }

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

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

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

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

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

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

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

442:    Logically Collective

444:    Input Parameters:
445: .  bv - basis vectors

447:    Note:
448:    All active columns (except the leading ones) are modified.

450:    Level: advanced

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

458:   PetscFunctionBegin;
461:   BVCheckSizes(bv,1);

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

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

474:    Logically Collective

476:    Input Parameters:
477: +  bv - basis vectors
478: -  j  - column number to be set

480:    Level: advanced

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

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

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

505:    Logically Collective

507:    Input Parameter:
508: .  bv - basis vectors

510:    Notes:
511:    All active columns (except the leading ones) are modified.

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

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

521:    Level: advanced

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

530:   PetscFunctionBegin;
533:   BVCheckSizes(bv,1);

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

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

555:    Logically Collective

557:    Input Parameter:
558: .  bv - basis vectors

560:    Notes:
561:    All active columns (except the leading ones) are modified.

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

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

571:    Level: advanced

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

580:   PetscFunctionBegin;
583:   BVCheckSizes(bv,1);

585:   PetscCall(BVGetRandomContext(bv,&bv->rand));
586:   PetscCall(PetscRandomGetInterval(bv->rand,&low,&high));
587:   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]");
588:   PetscCall(PetscLogEventBegin(BV_SetRandom,bv,0,0,0));
589:   for (k=bv->l;k<bv->k;k++) PetscCall(BVSetRandomSignColumn_Private(bv,k));
590:   PetscCall(PetscLogEventEnd(BV_SetRandom,bv,0,0,0));
591:   PetscCall(PetscObjectStateIncrease((PetscObject)bv));
592:   PetscFunctionReturn(PETSC_SUCCESS);
593: }

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

599:    Logically Collective

601:    Input Parameters:
602: +  bv    - basis vectors
603: -  condn - condition number

605:    Note:
606:    All active columns (except the leading ones) are modified.

608:    Level: advanced

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

619:   PetscFunctionBegin;
622:   BVCheckSizes(bv,1);

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

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

671:    Neighbor-wise Collective

673:    Input Parameters:
674: +  V - basis vectors context
675: -  A - the matrix

677:    Output Parameter:
678: .  Y - the result

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

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

690:    Level: beginner

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

698:   PetscFunctionBegin;
701:   BVCheckSizes(V,1);
702:   BVCheckOp(V,1,matmult);
707:   BVCheckSizes(Y,3);
708:   PetscCheckSameComm(V,1,A,2);
709:   PetscCheckSameTypeAndComm(V,1,Y,3);

711:   PetscCall(MatGetSize(A,&M,&N));
712:   PetscCall(MatGetLocalSize(A,&m,&n));
713:   PetscCheck(M==Y->N,PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Mismatching row dimension A %" PetscInt_FMT ", Y %" PetscInt_FMT,M,Y->N);
714:   PetscCheck(m==Y->n,PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Mismatching local row dimension A %" PetscInt_FMT ", Y %" PetscInt_FMT,m,Y->n);
715:   PetscCheck(N==V->N,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching column dimension A %" PetscInt_FMT ", V %" PetscInt_FMT,N,V->N);
716:   PetscCheck(n==V->n,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching local column dimension A %" PetscInt_FMT ", V %" PetscInt_FMT,n,V->n);
717:   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);

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

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

730:    Neighbor-wise Collective

732:    Input Parameters:
733: +  V - basis vectors context
734: -  A - the matrix

736:    Output Parameter:
737: .  Y - the result

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

746:    Currently implemented via MatCreateTranspose().

748:    Level: beginner

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

757:   PetscFunctionBegin;
760:   BVCheckSizes(V,1);
765:   BVCheckSizes(Y,3);
766:   PetscCheckSameComm(V,1,A,2);
767:   PetscCheckSameTypeAndComm(V,1,Y,3);

769:   PetscCall(MatGetSize(A,&M,&N));
770:   PetscCall(MatGetLocalSize(A,&m,&n));
771:   PetscCheck(M==V->N,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching row dimension A %" PetscInt_FMT ", V %" PetscInt_FMT,M,V->N);
772:   PetscCheck(m==V->n,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching local row dimension A %" PetscInt_FMT ", V %" PetscInt_FMT,m,V->n);
773:   PetscCheck(N==Y->N,PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Mismatching column dimension A %" PetscInt_FMT ", Y %" PetscInt_FMT,N,Y->N);
774:   PetscCheck(n==Y->n,PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Mismatching local column dimension A %" PetscInt_FMT ", Y %" PetscInt_FMT,n,Y->n);
775:   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);

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

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

787:    Neighbor-wise Collective

789:    Input Parameters:
790: +  V - basis vectors context
791: -  A - the matrix

793:    Output Parameter:
794: .  Y - the result

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

803:    Currently implemented via MatCreateHermitianTranspose().

805:    Level: beginner

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

814:   PetscFunctionBegin;
817:   BVCheckSizes(V,1);
822:   BVCheckSizes(Y,3);
823:   PetscCheckSameComm(V,1,A,2);
824:   PetscCheckSameTypeAndComm(V,1,Y,3);

826:   PetscCall(MatGetSize(A,&M,&N));
827:   PetscCall(MatGetLocalSize(A,&m,&n));
828:   PetscCheck(M==V->N,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching row dimension A %" PetscInt_FMT ", V %" PetscInt_FMT,M,V->N);
829:   PetscCheck(m==V->n,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching local row dimension A %" PetscInt_FMT ", V %" PetscInt_FMT,m,V->n);
830:   PetscCheck(N==Y->N,PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Mismatching column dimension A %" PetscInt_FMT ", Y %" PetscInt_FMT,N,Y->N);
831:   PetscCheck(n==Y->n,PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Mismatching local column dimension A %" PetscInt_FMT ", Y %" PetscInt_FMT,n,Y->n);
832:   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);

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

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

853:    Neighbor-wise Collective

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

860:    Level: beginner

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

868:   PetscFunctionBegin;
871:   BVCheckSizes(V,1);
873:   PetscCheckSameComm(V,1,A,2);
875:   PetscCheck(j>=0,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
876:   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);

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

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

893:    Neighbor-wise Collective

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

900:    Level: beginner

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

908:   PetscFunctionBegin;
911:   BVCheckSizes(V,1);
913:   PetscCheckSameComm(V,1,A,2);
915:   PetscCheck(j>=0,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
916:   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);

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

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

933:    Neighbor-wise Collective

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

940:    Level: beginner

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

948:   PetscFunctionBegin;
951:   BVCheckSizes(V,1);
953:   PetscCheckSameComm(V,1,A,2);
955:   PetscCheck(j>=0,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
956:   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);

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