1 #ifndef STAN_MATH_PRIM_FUN_CSR_MATRIX_TIMES_VECTOR_HPP
2 #define STAN_MATH_PRIM_FUN_CSR_MATRIX_TIMES_VECTOR_HPP
3 
4 #include <stan/math/prim/err.hpp>
5 #include <stan/math/prim/fun/csr_u_to_z.hpp>
6 #include <stan/math/prim/fun/Eigen.hpp>
7 #include <stan/math/prim/fun/dot_product.hpp>
8 #include <vector>
9 
10 namespace stan {
11 namespace math {
12 
13 /**
14  * @defgroup csr_format Compressed Sparse Row matrix format.
15  *  A compressed Sparse Row (CSR) sparse matrix is defined by four
16  *  component vectors labeled w, v, and u.  They are defined as:
17  *    - w: the non-zero values in the sparse matrix.
18  *    - v: column index for each value in w,  as a result this
19  *      is the same length as w.
20  *    - u: index of where each row starts in w,  length
21  *      is equal to the number of rows plus one.  Last entry is
22  *      one-past-the-end in w, following the Eigen spec.
23  *  Indexing is either zero-based or one-based depending on the value of
24  *  stan::error_index::value.  Following the definition of the format in
25  *  Eigen, we allow for unused garbage values in w/v which are never read.
26  *  All indexing _internal_ to a given function is zero-based.
27  *
28  *  With only m/n/w/v/u in hand, it is possible to check all
29  *  dimensions are sane _except_ the column dimension since it is
30  *  implicit.  The error-checking strategy is to check all dimensions
31  *  except the column dimension before any work is done inside a
32  *  function.  The column index is checked as it is constructed and
33  *  used for each entry.  If the column index is not needed it is
34  *  not checked.  As a result indexing mistakes might produce non-sensical
35  *  operations but out-of-bounds indexing will be caught.
36  *
37  *  Except for possible garbage values in w/v/u, memory usage is
38  *  calculated from the number of non-zero entries (NNZE) and the number of
39  *  rows (NR):  2*NNZE + 2*NR + 1.
40  */
41 
42 /**
43  * \addtogroup csr_format
44  * Return the multiplication of the sparse matrix (specified by
45  * by values and indexing) by the specified dense vector.
46  *
47  * The sparse matrix X of dimension m by n is represented by the
48  * vector w (of values), the integer array v (containing one-based
49  * column index of each value), the integer array u (containing
50  * one-based indexes of where each row starts in w).
51  *
52  * @tparam T1 type of the sparse matrix
53  * @tparam T2 type of the dense vector
54  * @param m Number of rows in matrix.
55  * @param n Number of columns in matrix.
56  * @param w Vector of non-zero values in matrix.
57  * @param v Column index of each non-zero value, same
58  *          length as w.
59  * @param u Index of where each row starts in w, length equal to
60  *          the number of rows plus one.
61  * @param b Eigen vector which the matrix is multiplied by.
62  * @return Dense vector for the product.
63  * @throw std::domain_error if m and n are not positive or are nan.
64  * @throw std::domain_error if the implied sparse matrix and b are
65  *                          not multiplicable.
66  * @throw std::invalid_argument if m/n/w/v/u are not internally
67  *   consistent, as defined by the indexing scheme.  Extractors are
68  *   defined in Stan which guarantee a consistent set of m/n/w/v/u
69  *   for a given sparse matrix.
70  * @throw std::out_of_range if any of the indexes are out of range.
71  */
72 template <typename T1, typename T2,
73           require_all_not_rev_matrix_t<T1, T2>* = nullptr>
74 inline Eigen::Matrix<return_type_t<T1, T2>, Eigen::Dynamic, 1>
csr_matrix_times_vector(int m,int n,const T1 & w,const std::vector<int> & v,const std::vector<int> & u,const T2 & b)75 csr_matrix_times_vector(int m, int n, const T1& w, const std::vector<int>& v,
76                         const std::vector<int>& u, const T2& b) {
77   check_positive("csr_matrix_times_vector", "m", m);
78   check_positive("csr_matrix_times_vector", "n", n);
79   check_size_match("csr_matrix_times_vector", "n", n, "b", b.size());
80   check_size_match("csr_matrix_times_vector", "w", w.size(), "v", v.size());
81   check_size_match("csr_matrix_times_vector", "m", m, "u", u.size() - 1);
82   check_size_match("csr_matrix_times_vector", "u/z",
83                    u[m - 1] + csr_u_to_z(u, m - 1) - 1, "v", v.size());
84   for (int i : v) {
85     check_range("csr_matrix_times_vector", "v[]", n, i);
86   }
87   std::vector<int> v_zero_based(v.size());
88   std::transform(v.begin(), v.end(), v_zero_based.begin(),
89                  [](auto&& x) { return x - 1; });
90   std::vector<int> u_zero_based(u.size());
91   std::transform(u.begin(), u.end(), u_zero_based.begin(),
92                  [](auto&& x) { return x - 1; });
93   // u_zero_based[u.size()] = w.size();
94   auto&& w_ref = to_ref(w);
95   Eigen::Map<const Eigen::SparseMatrix<scalar_type_t<T1>, Eigen::RowMajor>>
96       w_mat(m, n, w_ref.size(), u_zero_based.data(), v_zero_based.data(),
97             w_ref.data());
98   return w_mat * b;
99 }
100 
101 }  // namespace math
102 }  // namespace stan
103 
104 #endif
105