LCOV - code coverage report
Current view: top level - sys/classes/bv/interface - bvops.c (source / functions) Hit Total Coverage
Test: SLEPc Lines: 403 436 92.4 %
Date: 2024-03-29 00:55:19 Functions: 20 21 95.2 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /*
       2             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       3             :    SLEPc - Scalable Library for Eigenvalue Problem Computations
       4             :    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
       5             : 
       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             : */
      13             : 
      14             : #include <slepc/private/bvimpl.h>      /*I "slepcbv.h" I*/
      15             : #include <slepcds.h>
      16             : 
      17             : /*@
      18             :    BVMult - Computes Y = beta*Y + alpha*X*Q.
      19             : 
      20             :    Logically Collective
      21             : 
      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
      28             : 
      29             :    Notes:
      30             :    X and Y must be different objects. The case X=Y can be addressed with
      31             :    BVMultInPlace().
      32             : 
      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.
      39             : 
      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).
      44             : 
      45             :    Level: intermediate
      46             : 
      47             : .seealso: BVMultVec(), BVMultColumn(), BVMultInPlace(), BVSetActiveColumns()
      48             : @*/
      49        6823 : PetscErrorCode BVMult(BV Y,PetscScalar alpha,PetscScalar beta,BV X,Mat Q)
      50             : {
      51        6823 :   PetscInt       m,n;
      52             : 
      53        6823 :   PetscFunctionBegin;
      54        6823 :   PetscValidHeaderSpecific(Y,BV_CLASSID,1);
      55       27292 :   PetscValidLogicalCollectiveScalar(Y,alpha,2);
      56       27292 :   PetscValidLogicalCollectiveScalar(Y,beta,3);
      57        6823 :   PetscValidHeaderSpecific(X,BV_CLASSID,4);
      58        6823 :   if (Q) PetscValidHeaderSpecific(Q,MAT_CLASSID,5);
      59        6823 :   PetscValidType(Y,1);
      60        6823 :   BVCheckSizes(Y,1);
      61        6823 :   BVCheckOp(Y,1,mult);
      62        6823 :   PetscValidType(X,4);
      63        6823 :   BVCheckSizes(X,4);
      64        6823 :   if (Q) PetscValidType(Q,5);
      65        6823 :   PetscCheckSameTypeAndComm(Y,1,X,4);
      66        6823 :   PetscCheck(X!=Y,PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_WRONG,"X and Y arguments must be different");
      67        6823 :   if (Q) {
      68        5680 :     PetscCheckTypeNames(Q,MATSEQDENSE,MATSEQDENSECUDA);
      69        5680 :     PetscCall(MatGetSize(Q,&m,&n));
      70        5680 :     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        5680 :     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        6823 :   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);
      74             : 
      75        6823 :   PetscCall(PetscLogEventBegin(BV_Mult,X,Y,0,0));
      76        6823 :   PetscUseTypeMethod(Y,mult,alpha,beta,X,Q);
      77        6823 :   PetscCall(PetscLogEventEnd(BV_Mult,X,Y,0,0));
      78        6823 :   PetscCall(PetscObjectStateIncrease((PetscObject)Y));
      79        6823 :   PetscFunctionReturn(PETSC_SUCCESS);
      80             : }
      81             : 
      82             : /*@
      83             :    BVMultVec - Computes y = beta*y + alpha*X*q.
      84             : 
      85             :    Logically Collective
      86             : 
      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
      93             : 
      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().
      98             : 
      99             :    If X has leading columns specified, then these columns do not participate
     100             :    in the computation.
     101             : 
     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.
     105             : 
     106             :    Level: intermediate
     107             : 
     108             : .seealso: BVMult(), BVMultColumn(), BVMultInPlace(), BVSetActiveColumns()
     109             : @*/
     110      257039 : PetscErrorCode BVMultVec(BV X,PetscScalar alpha,PetscScalar beta,Vec y,PetscScalar q[])
     111             : {
     112      257039 :   PetscInt       n,N;
     113             : 
     114      257039 :   PetscFunctionBegin;
     115      257039 :   PetscValidHeaderSpecific(X,BV_CLASSID,1);
     116     1028156 :   PetscValidLogicalCollectiveScalar(X,alpha,2);
     117     1028156 :   PetscValidLogicalCollectiveScalar(X,beta,3);
     118      257039 :   PetscValidHeaderSpecific(y,VEC_CLASSID,4);
     119      257039 :   PetscAssertPointer(q,5);
     120      257039 :   PetscValidType(X,1);
     121      257039 :   BVCheckSizes(X,1);
     122      257039 :   BVCheckOp(X,1,multvec);
     123      257039 :   PetscValidType(y,4);
     124      257039 :   PetscCheckSameComm(X,1,y,4);
     125             : 
     126      257039 :   PetscCall(VecGetSize(y,&N));
     127      257039 :   PetscCall(VecGetLocalSize(y,&n));
     128      257039 :   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);
     129             : 
     130      257039 :   PetscCall(PetscLogEventBegin(BV_MultVec,X,y,0,0));
     131      257039 :   PetscUseTypeMethod(X,multvec,alpha,beta,y,q);
     132      257039 :   PetscCall(PetscLogEventEnd(BV_MultVec,X,y,0,0));
     133      257039 :   PetscFunctionReturn(PETSC_SUCCESS);
     134             : }
     135             : 
     136             : /*@
     137             :    BVMultColumn - Computes y = beta*y + alpha*X*q, where y is the j-th column
     138             :    of X.
     139             : 
     140             :    Logically Collective
     141             : 
     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
     148             : 
     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.
     156             : 
     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().
     160             : 
     161             :    Level: advanced
     162             : 
     163             : .seealso: BVMult(), BVMultVec(), BVMultInPlace(), BVSetActiveColumns()
     164             : @*/
     165      354096 : PetscErrorCode BVMultColumn(BV X,PetscScalar alpha,PetscScalar beta,PetscInt j,PetscScalar *q)
     166             : {
     167      354096 :   PetscInt       ksave;
     168      354096 :   Vec            y;
     169             : 
     170      354096 :   PetscFunctionBegin;
     171      354096 :   PetscValidHeaderSpecific(X,BV_CLASSID,1);
     172     1416384 :   PetscValidLogicalCollectiveScalar(X,alpha,2);
     173     1416384 :   PetscValidLogicalCollectiveScalar(X,beta,3);
     174     1416384 :   PetscValidLogicalCollectiveInt(X,j,4);
     175      354096 :   PetscValidType(X,1);
     176      354096 :   BVCheckSizes(X,1);
     177             : 
     178      354096 :   PetscCheck(j>=0,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
     179      354096 :   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);
     180             : 
     181      354096 :   PetscCall(PetscLogEventBegin(BV_MultVec,X,0,0,0));
     182      354096 :   ksave = X->k;
     183      354096 :   X->k = j;
     184      354096 :   if (!q && !X->buffer) PetscCall(BVGetBufferVec(X,&X->buffer));
     185      354096 :   PetscCall(BVGetColumn(X,j,&y));
     186      354096 :   PetscUseTypeMethod(X,multvec,alpha,beta,y,q);
     187      354096 :   PetscCall(BVRestoreColumn(X,j,&y));
     188      354096 :   X->k = ksave;
     189      354096 :   PetscCall(PetscLogEventEnd(BV_MultVec,X,0,0,0));
     190      354096 :   PetscCall(PetscObjectStateIncrease((PetscObject)X));
     191      354096 :   PetscFunctionReturn(PETSC_SUCCESS);
     192             : }
     193             : 
     194             : /*@
     195             :    BVMultInPlace - Update a set of vectors as V(:,s:e-1) = V*Q(:,s:e-1).
     196             : 
     197             :    Logically Collective
     198             : 
     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
     203             : 
     204             :    Input/Output Parameter:
     205             : .  V - basis vectors
     206             : 
     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).
     210             : 
     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.
     215             : 
     216             :    Level: intermediate
     217             : 
     218             : .seealso: BVMult(), BVMultVec(), BVMultInPlaceHermitianTranspose(), BVSetActiveColumns()
     219             : @*/
     220       28794 : PetscErrorCode BVMultInPlace(BV V,Mat Q,PetscInt s,PetscInt e)
     221             : {
     222       28794 :   PetscInt       m,n;
     223             : 
     224       28794 :   PetscFunctionBegin;
     225       28794 :   PetscValidHeaderSpecific(V,BV_CLASSID,1);
     226       28794 :   PetscValidHeaderSpecific(Q,MAT_CLASSID,2);
     227      115176 :   PetscValidLogicalCollectiveInt(V,s,3);
     228      115176 :   PetscValidLogicalCollectiveInt(V,e,4);
     229       28794 :   PetscValidType(V,1);
     230       28794 :   BVCheckSizes(V,1);
     231       28794 :   PetscValidType(Q,2);
     232       28794 :   PetscCheckTypeNames(Q,MATSEQDENSE,MATSEQDENSECUDA);
     233             : 
     234       28794 :   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       28794 :   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       28794 :   PetscCall(MatGetSize(Q,&m,&n));
     237       28794 :   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       28794 :   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             : 
     240       28794 :   PetscCall(PetscLogEventBegin(BV_MultInPlace,V,Q,0,0));
     241       28794 :   PetscUseTypeMethod(V,multinplace,Q,s,e);
     242       28794 :   PetscCall(PetscLogEventEnd(BV_MultInPlace,V,Q,0,0));
     243       28794 :   PetscCall(PetscObjectStateIncrease((PetscObject)V));
     244       28794 :   PetscFunctionReturn(PETSC_SUCCESS);
     245             : }
     246             : 
     247             : /*@
     248             :    BVMultInPlaceHermitianTranspose - Update a set of vectors as V(:,s:e-1) = V*Q'(:,s:e-1).
     249             : 
     250             :    Logically Collective
     251             : 
     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
     256             : 
     257             :    Input/Output Parameter:
     258             : .  V - basis vectors
     259             : 
     260             :    Notes:
     261             :    This is a variant of BVMultInPlace() where the conjugate transpose
     262             :    of Q is used.
     263             : 
     264             :    Level: intermediate
     265             : 
     266             : .seealso: BVMultInPlace()
     267             : @*/
     268           4 : PetscErrorCode BVMultInPlaceHermitianTranspose(BV V,Mat Q,PetscInt s,PetscInt e)
     269             : {
     270           4 :   PetscInt       m,n;
     271             : 
     272           4 :   PetscFunctionBegin;
     273           4 :   PetscValidHeaderSpecific(V,BV_CLASSID,1);
     274           4 :   PetscValidHeaderSpecific(Q,MAT_CLASSID,2);
     275          16 :   PetscValidLogicalCollectiveInt(V,s,3);
     276          16 :   PetscValidLogicalCollectiveInt(V,e,4);
     277           4 :   PetscValidType(V,1);
     278           4 :   BVCheckSizes(V,1);
     279           4 :   PetscValidType(Q,2);
     280           4 :   PetscCheckTypeNames(Q,MATSEQDENSE,MATSEQDENSECUDA);
     281             : 
     282           4 :   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           4 :   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           4 :   PetscCall(MatGetSize(Q,&m,&n));
     285           4 :   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           4 :   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);
     287             : 
     288           4 :   PetscCall(PetscLogEventBegin(BV_MultInPlace,V,Q,0,0));
     289           4 :   PetscUseTypeMethod(V,multinplacetrans,Q,s,e);
     290           4 :   PetscCall(PetscLogEventEnd(BV_MultInPlace,V,Q,0,0));
     291           4 :   PetscCall(PetscObjectStateIncrease((PetscObject)V));
     292           4 :   PetscFunctionReturn(PETSC_SUCCESS);
     293             : }
     294             : 
     295             : /*@
     296             :    BVScale - Multiply the BV entries by a scalar value.
     297             : 
     298             :    Logically Collective
     299             : 
     300             :    Input Parameters:
     301             : +  bv    - basis vectors
     302             : -  alpha - scaling factor
     303             : 
     304             :    Note:
     305             :    All active columns (except the leading ones) are scaled.
     306             : 
     307             :    Level: intermediate
     308             : 
     309             : .seealso: BVScaleColumn(), BVSetActiveColumns()
     310             : @*/
     311        1230 : PetscErrorCode BVScale(BV bv,PetscScalar alpha)
     312             : {
     313        1230 :   PetscFunctionBegin;
     314        1230 :   PetscValidHeaderSpecific(bv,BV_CLASSID,1);
     315        4920 :   PetscValidLogicalCollectiveScalar(bv,alpha,2);
     316        1230 :   PetscValidType(bv,1);
     317        1230 :   BVCheckSizes(bv,1);
     318        1230 :   if (alpha == (PetscScalar)1.0) PetscFunctionReturn(PETSC_SUCCESS);
     319             : 
     320         304 :   PetscCall(PetscLogEventBegin(BV_Scale,bv,0,0,0));
     321         304 :   PetscUseTypeMethod(bv,scale,-1,alpha);
     322         304 :   PetscCall(PetscLogEventEnd(BV_Scale,bv,0,0,0));
     323         304 :   PetscCall(PetscObjectStateIncrease((PetscObject)bv));
     324         304 :   PetscFunctionReturn(PETSC_SUCCESS);
     325             : }
     326             : 
     327             : /*@
     328             :    BVScaleColumn - Scale one column of a BV.
     329             : 
     330             :    Logically Collective
     331             : 
     332             :    Input Parameters:
     333             : +  bv    - basis vectors
     334             : .  j     - column number to be scaled
     335             : -  alpha - scaling factor
     336             : 
     337             :    Level: intermediate
     338             : 
     339             : .seealso: BVScale(), BVSetActiveColumns()
     340             : @*/
     341      132683 : PetscErrorCode BVScaleColumn(BV bv,PetscInt j,PetscScalar alpha)
     342             : {
     343      132683 :   PetscFunctionBegin;
     344      132683 :   PetscValidHeaderSpecific(bv,BV_CLASSID,1);
     345      530732 :   PetscValidLogicalCollectiveInt(bv,j,2);
     346      530732 :   PetscValidLogicalCollectiveScalar(bv,alpha,3);
     347      132683 :   PetscValidType(bv,1);
     348      132683 :   BVCheckSizes(bv,1);
     349             : 
     350      132683 :   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      132683 :   if (alpha == (PetscScalar)1.0) PetscFunctionReturn(PETSC_SUCCESS);
     352             : 
     353      131263 :   PetscCall(PetscLogEventBegin(BV_Scale,bv,0,0,0));
     354      131263 :   PetscUseTypeMethod(bv,scale,j,alpha);
     355      131263 :   PetscCall(PetscLogEventEnd(BV_Scale,bv,0,0,0));
     356      131263 :   PetscCall(PetscObjectStateIncrease((PetscObject)bv));
     357      131263 :   PetscFunctionReturn(PETSC_SUCCESS);
     358             : }
     359             : 
     360        2898 : static inline PetscErrorCode BVSetRandomColumn_Private(BV bv,PetscInt k)
     361             : {
     362        2898 :   PetscInt       i,low,high;
     363        2898 :   PetscScalar    *px,t;
     364        2898 :   Vec            x;
     365             : 
     366        2898 :   PetscFunctionBegin;
     367        2898 :   PetscCall(BVGetColumn(bv,k,&x));
     368        2898 :   if (bv->rrandom) {  /* generate the same vector irrespective of number of processes */
     369          80 :     PetscCall(VecGetOwnershipRange(x,&low,&high));
     370          80 :     PetscCall(VecGetArray(x,&px));
     371        5040 :     for (i=0;i<bv->N;i++) {
     372        4960 :       PetscCall(PetscRandomGetValue(bv->rand,&t));
     373        4960 :       if (i>=low && i<high) px[i-low] = t;
     374             :     }
     375          80 :     PetscCall(VecRestoreArray(x,&px));
     376        2818 :   } else PetscCall(VecSetRandom(x,bv->rand));
     377        2898 :   PetscCall(BVRestoreColumn(bv,k,&x));
     378        2898 :   PetscFunctionReturn(PETSC_SUCCESS);
     379             : }
     380             : 
     381         118 : static inline PetscErrorCode BVSetRandomNormalColumn_Private(BV bv,PetscInt k,Vec w1,Vec w2)
     382             : {
     383         118 :   PetscInt       i,low,high;
     384         118 :   PetscScalar    *px,s,t;
     385         118 :   Vec            x;
     386             : 
     387         118 :   PetscFunctionBegin;
     388         118 :   PetscCall(BVGetColumn(bv,k,&x));
     389         118 :   if (bv->rrandom) {  /* generate the same vector irrespective of number of processes */
     390           0 :     PetscCall(VecGetOwnershipRange(x,&low,&high));
     391           0 :     PetscCall(VecGetArray(x,&px));
     392           0 :     for (i=0;i<bv->N;i++) {
     393           0 :       PetscCall(PetscRandomGetValue(bv->rand,&s));
     394           0 :       PetscCall(PetscRandomGetValue(bv->rand,&t));
     395           0 :       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           0 :         px[i-low] = PetscSqrtReal(-2.0*PetscLogReal(s))*PetscCosReal(2.0*PETSC_PI*t);
     400             : #endif
     401             :       }
     402             :     }
     403           0 :     PetscCall(VecRestoreArray(x,&px));
     404         118 :   } else PetscCall(VecSetRandomNormal(x,bv->rand,w1,w2));
     405         118 :   PetscCall(BVRestoreColumn(bv,k,&x));
     406         118 :   PetscFunctionReturn(PETSC_SUCCESS);
     407             : }
     408             : 
     409         368 : static inline PetscErrorCode BVSetRandomSignColumn_Private(BV bv,PetscInt k)
     410             : {
     411         368 :   PetscInt       i,low,high;
     412         368 :   PetscScalar    *px,t;
     413         368 :   Vec            x;
     414             : 
     415         368 :   PetscFunctionBegin;
     416         368 :   PetscCall(BVGetColumn(bv,k,&x));
     417         368 :   PetscCall(VecGetOwnershipRange(x,&low,&high));
     418         368 :   if (bv->rrandom) {  /* generate the same vector irrespective of number of processes */
     419           0 :     PetscCall(VecGetArray(x,&px));
     420           0 :     for (i=0;i<bv->N;i++) {
     421           0 :       PetscCall(PetscRandomGetValue(bv->rand,&t));
     422           0 :       if (i>=low && i<high) px[i-low] = (PetscRealPart(t)<0.5)? -1.0: 1.0;
     423             :     }
     424           0 :     PetscCall(VecRestoreArray(x,&px));
     425             :   } else {
     426         368 :     PetscCall(VecSetRandom(x,bv->rand));
     427         368 :     PetscCall(VecGetArray(x,&px));
     428      116874 :     for (i=low;i<high;i++) {
     429      172964 :       px[i-low] = (PetscRealPart(px[i-low])<0.5)? -1.0: 1.0;
     430             :     }
     431         368 :     PetscCall(VecRestoreArray(x,&px));
     432             :   }
     433         368 :   PetscCall(BVRestoreColumn(bv,k,&x));
     434         368 :   PetscFunctionReturn(PETSC_SUCCESS);
     435             : }
     436             : 
     437             : /*@
     438             :    BVSetRandom - Set the columns of a BV to random numbers.
     439             : 
     440             :    Logically Collective
     441             : 
     442             :    Input Parameters:
     443             : .  bv - basis vectors
     444             : 
     445             :    Note:
     446             :    All active columns (except the leading ones) are modified.
     447             : 
     448             :    Level: advanced
     449             : 
     450             : .seealso: BVSetRandomContext(), BVSetRandomColumn(), BVSetRandomNormal(), BVSetRandomSign(), BVSetRandomCond(), BVSetActiveColumns()
     451             : @*/
     452          16 : PetscErrorCode BVSetRandom(BV bv)
     453             : {
     454          16 :   PetscInt       k;
     455             : 
     456          16 :   PetscFunctionBegin;
     457          16 :   PetscValidHeaderSpecific(bv,BV_CLASSID,1);
     458          16 :   PetscValidType(bv,1);
     459          16 :   BVCheckSizes(bv,1);
     460             : 
     461          16 :   PetscCall(BVGetRandomContext(bv,&bv->rand));
     462          16 :   PetscCall(PetscLogEventBegin(BV_SetRandom,bv,0,0,0));
     463          96 :   for (k=bv->l;k<bv->k;k++) PetscCall(BVSetRandomColumn_Private(bv,k));
     464          16 :   PetscCall(PetscLogEventEnd(BV_SetRandom,bv,0,0,0));
     465          16 :   PetscCall(PetscObjectStateIncrease((PetscObject)bv));
     466          16 :   PetscFunctionReturn(PETSC_SUCCESS);
     467             : }
     468             : 
     469             : /*@
     470             :    BVSetRandomColumn - Set one column of a BV to random numbers.
     471             : 
     472             :    Logically Collective
     473             : 
     474             :    Input Parameters:
     475             : +  bv - basis vectors
     476             : -  j  - column number to be set
     477             : 
     478             :    Level: advanced
     479             : 
     480             : .seealso: BVSetRandomContext(), BVSetRandom(), BVSetRandomNormal(), BVSetRandomCond()
     481             : @*/
     482        2786 : PetscErrorCode BVSetRandomColumn(BV bv,PetscInt j)
     483             : {
     484        2786 :   PetscFunctionBegin;
     485        2786 :   PetscValidHeaderSpecific(bv,BV_CLASSID,1);
     486       11144 :   PetscValidLogicalCollectiveInt(bv,j,2);
     487        2786 :   PetscValidType(bv,1);
     488        2786 :   BVCheckSizes(bv,1);
     489        2786 :   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);
     490             : 
     491        2786 :   PetscCall(BVGetRandomContext(bv,&bv->rand));
     492        2786 :   PetscCall(PetscLogEventBegin(BV_SetRandom,bv,0,0,0));
     493        2786 :   PetscCall(BVSetRandomColumn_Private(bv,j));
     494        2786 :   PetscCall(PetscLogEventEnd(BV_SetRandom,bv,0,0,0));
     495        2786 :   PetscCall(PetscObjectStateIncrease((PetscObject)bv));
     496        2786 :   PetscFunctionReturn(PETSC_SUCCESS);
     497             : }
     498             : 
     499             : /*@
     500             :    BVSetRandomNormal - Set the columns of a BV to random numbers with a normal
     501             :    distribution.
     502             : 
     503             :    Logically Collective
     504             : 
     505             :    Input Parameter:
     506             : .  bv - basis vectors
     507             : 
     508             :    Notes:
     509             :    All active columns (except the leading ones) are modified.
     510             : 
     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).
     514             : 
     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.
     518             : 
     519             :    Level: advanced
     520             : 
     521             : .seealso: BVSetRandomContext(), BVSetRandom(), BVSetRandomColumn(), BVSetRandomCond(), BVSetActiveColumns()
     522             : @*/
     523          11 : PetscErrorCode BVSetRandomNormal(BV bv)
     524             : {
     525          11 :   PetscInt       k;
     526          11 :   Vec            w1=NULL,w2=NULL;
     527             : 
     528          11 :   PetscFunctionBegin;
     529          11 :   PetscValidHeaderSpecific(bv,BV_CLASSID,1);
     530          11 :   PetscValidType(bv,1);
     531          11 :   BVCheckSizes(bv,1);
     532             : 
     533          11 :   PetscCall(BVGetRandomContext(bv,&bv->rand));
     534          11 :   if (!bv->rrandom) {
     535          11 :     PetscCall(BVCreateVec(bv,&w1));
     536          11 :     PetscCall(BVCreateVec(bv,&w2));
     537             :   }
     538          11 :   PetscCall(PetscLogEventBegin(BV_SetRandom,bv,0,0,0));
     539         129 :   for (k=bv->l;k<bv->k;k++) PetscCall(BVSetRandomNormalColumn_Private(bv,k,w1,w2));
     540          11 :   PetscCall(PetscLogEventEnd(BV_SetRandom,bv,0,0,0));
     541          11 :   if (!bv->rrandom) {
     542          11 :     PetscCall(VecDestroy(&w1));
     543          11 :     PetscCall(VecDestroy(&w2));
     544             :   }
     545          11 :   PetscCall(PetscObjectStateIncrease((PetscObject)bv));
     546          11 :   PetscFunctionReturn(PETSC_SUCCESS);
     547             : }
     548             : 
     549             : /*@
     550             :    BVSetRandomSign - Set the entries of a BV to values 1 or -1 with equal
     551             :    probability.
     552             : 
     553             :    Logically Collective
     554             : 
     555             :    Input Parameter:
     556             : .  bv - basis vectors
     557             : 
     558             :    Notes:
     559             :    All active columns (except the leading ones) are modified.
     560             : 
     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].
     564             : 
     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.
     568             : 
     569             :    Level: advanced
     570             : 
     571             : .seealso: BVSetRandomContext(), BVSetRandom(), BVSetRandomColumn(), BVSetActiveColumns()
     572             : @*/
     573          26 : PetscErrorCode BVSetRandomSign(BV bv)
     574             : {
     575          26 :   PetscScalar    low,high;
     576          26 :   PetscInt       k;
     577             : 
     578          26 :   PetscFunctionBegin;
     579          26 :   PetscValidHeaderSpecific(bv,BV_CLASSID,1);
     580          26 :   PetscValidType(bv,1);
     581          26 :   BVCheckSizes(bv,1);
     582             : 
     583          26 :   PetscCall(BVGetRandomContext(bv,&bv->rand));
     584          26 :   PetscCall(PetscRandomGetInterval(bv->rand,&low,&high));
     585          26 :   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          26 :   PetscCall(PetscLogEventBegin(BV_SetRandom,bv,0,0,0));
     587         394 :   for (k=bv->l;k<bv->k;k++) PetscCall(BVSetRandomSignColumn_Private(bv,k));
     588          26 :   PetscCall(PetscLogEventEnd(BV_SetRandom,bv,0,0,0));
     589          26 :   PetscCall(PetscObjectStateIncrease((PetscObject)bv));
     590          26 :   PetscFunctionReturn(PETSC_SUCCESS);
     591             : }
     592             : 
     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.
     596             : 
     597             :    Logically Collective
     598             : 
     599             :    Input Parameters:
     600             : +  bv    - basis vectors
     601             : -  condn - condition number
     602             : 
     603             :    Note:
     604             :    All active columns (except the leading ones) are modified.
     605             : 
     606             :    Level: advanced
     607             : 
     608             : .seealso: BVSetRandomContext(), BVSetRandom(), BVSetRandomColumn(), BVSetRandomNormal(), BVSetActiveColumns()
     609             : @*/
     610           4 : PetscErrorCode BVSetRandomCond(BV bv,PetscReal condn)
     611             : {
     612           4 :   PetscInt       k,i;
     613           4 :   PetscScalar    *eig,*d;
     614           4 :   DS             ds;
     615           4 :   Mat            A,X,Xt,M,G;
     616             : 
     617           4 :   PetscFunctionBegin;
     618           4 :   PetscValidHeaderSpecific(bv,BV_CLASSID,1);
     619           4 :   PetscValidType(bv,1);
     620           4 :   BVCheckSizes(bv,1);
     621             : 
     622           4 :   PetscCall(BVGetRandomContext(bv,&bv->rand));
     623           4 :   PetscCall(PetscLogEventBegin(BV_SetRandom,bv,0,0,0));
     624             :   /* B = rand(n,k) */
     625          36 :   for (k=bv->l;k<bv->k;k++) PetscCall(BVSetRandomColumn_Private(bv,k));
     626           4 :   PetscCall(DSCreate(PetscObjectComm((PetscObject)bv),&ds));
     627           4 :   PetscCall(DSSetType(ds,DSHEP));
     628           4 :   PetscCall(DSAllocate(ds,bv->m));
     629           4 :   PetscCall(DSSetDimensions(ds,bv->k,bv->l,bv->k));
     630             :   /* [V,S] = eig(B'*B) */
     631           4 :   PetscCall(DSGetMat(ds,DS_MAT_A,&A));
     632           4 :   PetscCall(BVDot(bv,bv,A));
     633           4 :   PetscCall(DSRestoreMat(ds,DS_MAT_A,&A));
     634           4 :   PetscCall(PetscMalloc1(bv->m,&eig));
     635           4 :   PetscCall(DSSolve(ds,eig,NULL));
     636           4 :   PetscCall(DSSynchronize(ds,eig,NULL));
     637           4 :   PetscCall(DSVectors(ds,DS_MAT_X,NULL,NULL));
     638             :   /* M = diag(linspace(1/condn,1,n)./sqrt(diag(S)))' */
     639           4 :   PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,bv->k,bv->k,NULL,&M));
     640           4 :   PetscCall(MatZeroEntries(M));
     641           4 :   PetscCall(MatDenseGetArray(M,&d));
     642          36 :   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           4 :   PetscCall(MatDenseRestoreArray(M,&d));
     644             :   /* G = X*M*X' */
     645           4 :   PetscCall(DSGetMat(ds,DS_MAT_X,&X));
     646           4 :   PetscCall(MatTranspose(X,MAT_INITIAL_MATRIX,&Xt));
     647           4 :   PetscCall(MatProductCreate(Xt,M,NULL,&G));
     648           4 :   PetscCall(MatProductSetType(G,MATPRODUCT_PtAP));
     649           4 :   PetscCall(MatProductSetFromOptions(G));
     650           4 :   PetscCall(MatProductSymbolic(G));
     651           4 :   PetscCall(MatProductNumeric(G));
     652           4 :   PetscCall(MatProductClear(G));
     653           4 :   PetscCall(DSRestoreMat(ds,DS_MAT_X,&X));
     654           4 :   PetscCall(MatDestroy(&Xt));
     655           4 :   PetscCall(MatDestroy(&M));
     656             :   /* B = B*G */
     657           4 :   PetscCall(BVMultInPlace(bv,G,bv->l,bv->k));
     658           4 :   PetscCall(MatDestroy(&G));
     659           4 :   PetscCall(PetscFree(eig));
     660           4 :   PetscCall(DSDestroy(&ds));
     661           4 :   PetscCall(PetscLogEventEnd(BV_SetRandom,bv,0,0,0));
     662           4 :   PetscCall(PetscObjectStateIncrease((PetscObject)bv));
     663           4 :   PetscFunctionReturn(PETSC_SUCCESS);
     664             : }
     665             : 
     666             : /*@
     667             :    BVMatMult - Computes the matrix-vector product for each column, Y=A*V.
     668             : 
     669             :    Neighbor-wise Collective
     670             : 
     671             :    Input Parameters:
     672             : +  V - basis vectors context
     673             : -  A - the matrix
     674             : 
     675             :    Output Parameter:
     676             : .  Y - the result
     677             : 
     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.
     684             : 
     685             :    It is possible to choose whether the computation is done column by column
     686             :    or as a Mat-Mat product, see BVSetMatMultMethod().
     687             : 
     688             :    Level: beginner
     689             : 
     690             : .seealso: BVCopy(), BVSetActiveColumns(), BVMatMultColumn(), BVMatMultTranspose(), BVMatMultHermitianTranspose(), BVSetMatMultMethod()
     691             : @*/
     692       11588 : PetscErrorCode BVMatMult(BV V,Mat A,BV Y)
     693             : {
     694       11588 :   PetscInt       M,N,m,n;
     695             : 
     696       11588 :   PetscFunctionBegin;
     697       11588 :   PetscValidHeaderSpecific(V,BV_CLASSID,1);
     698       11588 :   PetscValidType(V,1);
     699       11588 :   BVCheckSizes(V,1);
     700       11588 :   BVCheckOp(V,1,matmult);
     701       11588 :   PetscValidHeaderSpecific(A,MAT_CLASSID,2);
     702       11588 :   PetscValidType(A,2);
     703       11588 :   PetscValidHeaderSpecific(Y,BV_CLASSID,3);
     704       11588 :   PetscValidType(Y,3);
     705       11588 :   BVCheckSizes(Y,3);
     706       11588 :   PetscCheckSameComm(V,1,A,2);
     707       11588 :   PetscCheckSameTypeAndComm(V,1,Y,3);
     708             : 
     709       11588 :   PetscCall(MatGetSize(A,&M,&N));
     710       11588 :   PetscCall(MatGetLocalSize(A,&m,&n));
     711       11588 :   PetscCheck(M==Y->N,PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Mismatching row dimension A %" PetscInt_FMT ", Y %" PetscInt_FMT,M,Y->N);
     712       11588 :   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       11588 :   PetscCheck(N==V->N,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching column dimension A %" PetscInt_FMT ", V %" PetscInt_FMT,N,V->N);
     714       11588 :   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       11588 :   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);
     716             : 
     717       11588 :   PetscCall(PetscLogEventBegin(BV_MatMult,V,A,Y,0));
     718       11588 :   PetscUseTypeMethod(V,matmult,A,Y);
     719       11588 :   PetscCall(PetscLogEventEnd(BV_MatMult,V,A,Y,0));
     720       11588 :   PetscCall(PetscObjectStateIncrease((PetscObject)Y));
     721       11588 :   PetscFunctionReturn(PETSC_SUCCESS);
     722             : }
     723             : 
     724             : /*@
     725             :    BVMatMultTranspose - Computes the matrix-vector product with the transpose
     726             :    of a matrix for each column, Y=A^T*V.
     727             : 
     728             :    Neighbor-wise Collective
     729             : 
     730             :    Input Parameters:
     731             : +  V - basis vectors context
     732             : -  A - the matrix
     733             : 
     734             :    Output Parameter:
     735             : .  Y - the result
     736             : 
     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.
     743             : 
     744             :    Currently implemented via MatCreateTranspose().
     745             : 
     746             :    Level: beginner
     747             : 
     748             : .seealso: BVMatMult(), BVMatMultHermitianTranspose()
     749             : @*/
     750          16 : PetscErrorCode BVMatMultTranspose(BV V,Mat A,BV Y)
     751             : {
     752          16 :   PetscInt       M,N,m,n;
     753          16 :   Mat            AT;
     754             : 
     755          16 :   PetscFunctionBegin;
     756          16 :   PetscValidHeaderSpecific(V,BV_CLASSID,1);
     757          16 :   PetscValidType(V,1);
     758          16 :   BVCheckSizes(V,1);
     759          16 :   PetscValidHeaderSpecific(A,MAT_CLASSID,2);
     760          16 :   PetscValidType(A,2);
     761          16 :   PetscValidHeaderSpecific(Y,BV_CLASSID,3);
     762          16 :   PetscValidType(Y,3);
     763          16 :   BVCheckSizes(Y,3);
     764          16 :   PetscCheckSameComm(V,1,A,2);
     765          16 :   PetscCheckSameTypeAndComm(V,1,Y,3);
     766             : 
     767          16 :   PetscCall(MatGetSize(A,&M,&N));
     768          16 :   PetscCall(MatGetLocalSize(A,&m,&n));
     769          16 :   PetscCheck(M==V->N,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching row dimension A %" PetscInt_FMT ", V %" PetscInt_FMT,M,V->N);
     770          16 :   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          16 :   PetscCheck(N==Y->N,PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Mismatching column dimension A %" PetscInt_FMT ", Y %" PetscInt_FMT,N,Y->N);
     772          16 :   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          16 :   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);
     774             : 
     775          16 :   PetscCall(MatCreateTranspose(A,&AT));
     776          16 :   PetscCall(BVMatMult(V,AT,Y));
     777          16 :   PetscCall(MatDestroy(&AT));
     778          16 :   PetscFunctionReturn(PETSC_SUCCESS);
     779             : }
     780             : 
     781             : /*@
     782             :    BVMatMultHermitianTranspose - Computes the matrix-vector product with the
     783             :    conjugate transpose of a matrix for each column, Y=A^H*V.
     784             : 
     785             :    Neighbor-wise Collective
     786             : 
     787             :    Input Parameters:
     788             : +  V - basis vectors context
     789             : -  A - the matrix
     790             : 
     791             :    Output Parameter:
     792             : .  Y - the result
     793             : 
     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.
     800             : 
     801             :    Currently implemented via MatCreateHermitianTranspose().
     802             : 
     803             :    Level: beginner
     804             : 
     805             : .seealso: BVMatMult(), BVMatMultTranspose()
     806             : @*/
     807          67 : PetscErrorCode BVMatMultHermitianTranspose(BV V,Mat A,BV Y)
     808             : {
     809          67 :   PetscInt       j,M,N,m,n;
     810          67 :   Vec            v,y;
     811             : 
     812          67 :   PetscFunctionBegin;
     813          67 :   PetscValidHeaderSpecific(V,BV_CLASSID,1);
     814          67 :   PetscValidType(V,1);
     815          67 :   BVCheckSizes(V,1);
     816          67 :   PetscValidHeaderSpecific(A,MAT_CLASSID,2);
     817          67 :   PetscValidType(A,2);
     818          67 :   PetscValidHeaderSpecific(Y,BV_CLASSID,3);
     819          67 :   PetscValidType(Y,3);
     820          67 :   BVCheckSizes(Y,3);
     821          67 :   PetscCheckSameComm(V,1,A,2);
     822          67 :   PetscCheckSameTypeAndComm(V,1,Y,3);
     823             : 
     824          67 :   PetscCall(MatGetSize(A,&M,&N));
     825          67 :   PetscCall(MatGetLocalSize(A,&m,&n));
     826          67 :   PetscCheck(M==V->N,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching row dimension A %" PetscInt_FMT ", V %" PetscInt_FMT,M,V->N);
     827          67 :   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          67 :   PetscCheck(N==Y->N,PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Mismatching column dimension A %" PetscInt_FMT ", Y %" PetscInt_FMT,N,Y->N);
     829          67 :   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          67 :   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);
     831             : 
     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         174 :   for (j=0;j<V->k-V->l;j++) {
     838         107 :     PetscCall(BVGetColumn(V,V->l+j,&v));
     839         107 :     PetscCall(BVGetColumn(Y,Y->l+j,&y));
     840         107 :     PetscCall(MatMultHermitianTranspose(A,v,y));
     841         107 :     PetscCall(BVRestoreColumn(V,V->l+j,&v));
     842         107 :     PetscCall(BVRestoreColumn(Y,Y->l+j,&y));
     843             :   }
     844          67 :   PetscFunctionReturn(PETSC_SUCCESS);
     845             : }
     846             : 
     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.
     850             : 
     851             :    Neighbor-wise Collective
     852             : 
     853             :    Input Parameters:
     854             : +  V - basis vectors context
     855             : .  A - the matrix
     856             : -  j - the column
     857             : 
     858             :    Level: beginner
     859             : 
     860             : .seealso: BVMatMult(), BVMatMultTransposeColumn(), BVMatMultHermitianTransposeColumn()
     861             : @*/
     862       68314 : PetscErrorCode BVMatMultColumn(BV V,Mat A,PetscInt j)
     863             : {
     864       68314 :   Vec            vj,vj1;
     865             : 
     866       68314 :   PetscFunctionBegin;
     867       68314 :   PetscValidHeaderSpecific(V,BV_CLASSID,1);
     868       68314 :   PetscValidType(V,1);
     869       68314 :   BVCheckSizes(V,1);
     870       68314 :   PetscValidHeaderSpecific(A,MAT_CLASSID,2);
     871       68314 :   PetscCheckSameComm(V,1,A,2);
     872      273256 :   PetscValidLogicalCollectiveInt(V,j,3);
     873       68314 :   PetscCheck(j>=0,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
     874       68314 :   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);
     875             : 
     876       68314 :   PetscCall(PetscLogEventBegin(BV_MatMultVec,V,A,0,0));
     877       68314 :   PetscCall(BVGetColumn(V,j,&vj));
     878       68314 :   PetscCall(BVGetColumn(V,j+1,&vj1));
     879       68314 :   PetscCall(MatMult(A,vj,vj1));
     880       68314 :   PetscCall(BVRestoreColumn(V,j,&vj));
     881       68314 :   PetscCall(BVRestoreColumn(V,j+1,&vj1));
     882       68314 :   PetscCall(PetscLogEventEnd(BV_MatMultVec,V,A,0,0));
     883       68314 :   PetscCall(PetscObjectStateIncrease((PetscObject)V));
     884       68314 :   PetscFunctionReturn(PETSC_SUCCESS);
     885             : }
     886             : 
     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.
     890             : 
     891             :    Neighbor-wise Collective
     892             : 
     893             :    Input Parameters:
     894             : +  V - basis vectors context
     895             : .  A - the matrix
     896             : -  j - the column
     897             : 
     898             :    Level: beginner
     899             : 
     900             : .seealso: BVMatMult(), BVMatMultColumn(), BVMatMultHermitianTransposeColumn()
     901             : @*/
     902          16 : PetscErrorCode BVMatMultTransposeColumn(BV V,Mat A,PetscInt j)
     903             : {
     904          16 :   Vec            vj,vj1;
     905             : 
     906          16 :   PetscFunctionBegin;
     907          16 :   PetscValidHeaderSpecific(V,BV_CLASSID,1);
     908          16 :   PetscValidType(V,1);
     909          16 :   BVCheckSizes(V,1);
     910          16 :   PetscValidHeaderSpecific(A,MAT_CLASSID,2);
     911          16 :   PetscCheckSameComm(V,1,A,2);
     912          64 :   PetscValidLogicalCollectiveInt(V,j,3);
     913          16 :   PetscCheck(j>=0,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
     914          16 :   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);
     915             : 
     916          16 :   PetscCall(PetscLogEventBegin(BV_MatMultVec,V,A,0,0));
     917          16 :   PetscCall(BVGetColumn(V,j,&vj));
     918          16 :   PetscCall(BVGetColumn(V,j+1,&vj1));
     919          16 :   PetscCall(MatMultTranspose(A,vj,vj1));
     920          16 :   PetscCall(BVRestoreColumn(V,j,&vj));
     921          16 :   PetscCall(BVRestoreColumn(V,j+1,&vj1));
     922          16 :   PetscCall(PetscLogEventEnd(BV_MatMultVec,V,A,0,0));
     923          16 :   PetscCall(PetscObjectStateIncrease((PetscObject)V));
     924          16 :   PetscFunctionReturn(PETSC_SUCCESS);
     925             : }
     926             : 
     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.
     930             : 
     931             :    Neighbor-wise Collective
     932             : 
     933             :    Input Parameters:
     934             : +  V - basis vectors context
     935             : .  A - the matrix
     936             : -  j - the column
     937             : 
     938             :    Level: beginner
     939             : 
     940             : .seealso: BVMatMult(), BVMatMultColumn(), BVMatMultTransposeColumn()
     941             : @*/
     942           0 : PetscErrorCode BVMatMultHermitianTransposeColumn(BV V,Mat A,PetscInt j)
     943             : {
     944           0 :   Vec            vj,vj1;
     945             : 
     946           0 :   PetscFunctionBegin;
     947           0 :   PetscValidHeaderSpecific(V,BV_CLASSID,1);
     948           0 :   PetscValidType(V,1);
     949           0 :   BVCheckSizes(V,1);
     950           0 :   PetscValidHeaderSpecific(A,MAT_CLASSID,2);
     951           0 :   PetscCheckSameComm(V,1,A,2);
     952           0 :   PetscValidLogicalCollectiveInt(V,j,3);
     953           0 :   PetscCheck(j>=0,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
     954           0 :   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);
     955             : 
     956           0 :   PetscCall(PetscLogEventBegin(BV_MatMultVec,V,A,0,0));
     957           0 :   PetscCall(BVGetColumn(V,j,&vj));
     958           0 :   PetscCall(BVGetColumn(V,j+1,&vj1));
     959           0 :   PetscCall(MatMultHermitianTranspose(A,vj,vj1));
     960           0 :   PetscCall(BVRestoreColumn(V,j,&vj));
     961           0 :   PetscCall(BVRestoreColumn(V,j+1,&vj1));
     962           0 :   PetscCall(PetscLogEventEnd(BV_MatMultVec,V,A,0,0));
     963           0 :   PetscCall(PetscObjectStateIncrease((PetscObject)V));
     964           0 :   PetscFunctionReturn(PETSC_SUCCESS);
     965             : }

Generated by: LCOV version 1.14