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