Merge pull request #10 from voduchuy/improve_sensitivity

Extend the sensitivity FSP to allow for the case where the time-indep…

src/Matrix/FspMatrixBase.cpp

@@ -38,6 +38,8 @@ int FspMatrixBase::Action(PetscReal t, Vec x, Vec y) {
 
   ierr = VecSet(y, 0.0);
   CHKERRQ(ierr);
+  if (has_values_ == PETSC_FALSE) return 0;
+
   if (!tv_reactions_.empty()) {
     ierr = t_fun_(t, num_reactions_, time_coefficients_.memptr(), t_fun_args_);
     PACMENSLCHKERRQ(ierr);
@@ -244,6 +246,7 @@ PacmenslErrorCode FspMatrixBase::GenerateValues(const StateSetBase &fsp,
     ierr = MatAssemblyEnd(ti_mat_, MAT_FINAL_ASSEMBLY);
     CHKERRQ(ierr);
   }
+  has_values_ = PETSC_TRUE;
   return 0;
 }
 
@@ -266,7 +269,7 @@ int FspMatrixBase::Destroy() {
     ierr = VecDestroy(work_.mem());
     CHKERRQ(ierr);
   }
-
+  has_values_ = PETSC_FALSE;
   return 0;
 }
 

src/Matrix/FspMatrixBase.h

@@ -58,7 +58,7 @@ class FspMatrixBase {
    * @param comm Communicator context for the processes that own the matrix.
    */
   explicit FspMatrixBase(MPI_Comm comm);
-  
+
   /**
    * @brief Populate the fields of the FSP matrix object.
    * This method is __collective__, meaning that it must be called by all processes that own this object.
@@ -69,13 +69,13 @@ class FspMatrixBase {
   virtual PacmenslErrorCode
   GenerateValues(const StateSetBase &fsp,
                  const Model &model);
-  
+
   /**
    * @brief Populate the fields of the FSP matrix object.
    * This method is __collective__, meaning that it must be called by all processes that own this object.
    * @param fsp (in) The FSP-truncated state space from which the transition rate matrix is built.
    * @param SM (in) The stoichiometry matrix of the reaction network.
-   * @param time_varying (in) List of reactions whose propensities are time-varying.
+   * @param time_varying (in) List of reactions whose propensities are time-varying. If empty, all reaction propensities are time-independent.
    * @param new_prop_t (in) Function pointer to evaluate the vector of time-varying coefficients $c_r(t,\theta)$ at time $t$.
    * @param new_prop_x (in) Function pointer to evaluate the time-independent component of the propensities.
    * @param enable_reactions (in) Vector of reactions that are active. If left empty, we assume all reactions are active. Propensity functions of inactive reactions are not added to the data structure.
@@ -92,11 +92,15 @@ class FspMatrixBase {
                  const std::vector<int> &enable_reactions,
                  void *prop_t_args,
                  void *prop_x_args);
-  
+
   PacmenslErrorCode SetTimeFun(TcoefFun new_t_fun, void *new_t_fun_args);
-
+
+  /**
+   * @brief Clear internal state of the object.
+   * @return Error code: 0 (success) or -1 (failure).
+   */
   virtual int Destroy();
-  
+
   /**
   * @brief Compute y = A*x.
   * This method is __collective__, meaning that it must be called by all processes that own this object.
@@ -106,15 +110,15 @@ class FspMatrixBase {
   * @return Error code: 0 (success) or -1 (failure)
   */
   virtual PacmenslErrorCode Action(PetscReal t, Vec x, Vec y);
-  
+
   /**
    * @brief Generate a PETSc Matrix object with the same sparsity structure as the FSP matrix. This serves as the Jacobian matrix required by ODE solvers.
    * This method is __collective__, meaning that it must be called by all processes that own this object.
    * @param A (out) pointer to PETSc Mat object.
    * @return Error code. 0 (success) or -1 (failure)
    */
   virtual PacmenslErrorCode CreateRHSJacobian(Mat *A);
-  
+
   /**
    * @brief Populate the values of a PETSc Mat object with the entries of the FSP matrix at any given time.
    * This method is __collective__, meaning that it must be called by all processes that own this object.
@@ -124,22 +128,22 @@ class FspMatrixBase {
    * @return Error code: 0 (success) or -1 (failure)
    */
   virtual PacmenslErrorCode ComputeRHSJacobian(PetscReal t, Mat A);
-  
+
   /**
    * @brief Get an estimate for the number of floating-point operations needed by the calling process per matrix-vector multiplication (i.e., per call to the \ref Action() method).
    * This method is __local__, meaning that it can be called independently by the local process.
    * @param nflops (out) number of flops.
    * @return Error code: 0 (success) or -1 (failure)
    */
   virtual PacmenslErrorCode GetLocalMVFlops(PetscInt *nflops);
-  
+
   /**
    * @brief Get the number of rows owned by the calling process.
    * This method is __local__, meaning that it can be called independently by the local process.
    * @return Number of rows owned by the calling process.
    */
   int GetNumLocalRows() const { return num_rows_local_; };
-  
+
   /**
    * @brief Object destructor.
    */
@@ -148,26 +152,26 @@ class FspMatrixBase {
   MPI_Comm comm_ = MPI_COMM_NULL; ///< Communicator context
   int rank_; ///< Rank of local process
   int comm_size_; ///< Total number of processes that own this object
-  
+
   Int num_reactions_ = 0; ///< Number of reactions
   Int num_rows_global_ = 0; ///< Global number of rows
   Int num_rows_local_ = 0; ///< Number of matrix rows owned by this process
-  
+
   std::vector<int> enable_reactions_ = std::vector<int>();
   std::vector<int> tv_reactions_ = std::vector<int>();
   std::vector<int> ti_reactions_ = std::vector<int>();
-  
+
   // Local data of the matrix
   std::vector<Petsc<Mat>> tv_mats_; ///< List of PETSc Mat objects for storing the $A_r$ factors for reactions $r$ with time-varying propensities
   Petsc<Mat> ti_mat_; ///< PETsc Mat object to store the time-independent term of the matrix sum decomposition
-  
+
   // Data for computing the matrix action
   Petsc<Vec> work_; ///< Work vector for computing operator times vector
-  
+
   TcoefFun t_fun_ = nullptr;  ///< Pointer to external function that evaluates time-varying factors of reaction propensities.
   void *t_fun_args_ = nullptr; ///< Pointer for extra data needed in \ref t_fun_ evaluation.
   arma::Row<Real> time_coefficients_; ///< Vector to store the time-varying factors of reaction propensities
-  
+
   /**
    * @brief Populate the data members with information about the inter-process layout of the object.
    * This method is __collective__, meaning that it must be called by all processes that own this object.
@@ -176,7 +180,7 @@ class FspMatrixBase {
    * Modified fields: \ref num_rows_local_, \ref num_rows_global_, \ref work_.
    */
   virtual int DetermineLayout_(const StateSetBase &fsp);
-  
+
   // arrays for counting nonzero entries on the diagonal and off-diagonal blocks
   arma::Mat<Int> dblock_nz_, oblock_nz_;
   arma::Col<Int> ti_dblock_nz_, ti_oblock_nz_;
@@ -185,6 +189,8 @@ class FspMatrixBase {
   // array of matrix values
   arma::Mat<PetscReal> offdiag_vals_;
   arma::Mat<PetscReal> diag_vals_;
+
+  PetscBool has_values_ = PETSC_FALSE; ///< whether the matrix has values.
 };
 
 }