introduce sparse matrix multiplication for big models

[stemangiola]

e.g.

data {
  int<lower=1> N;
  int<lower=1> K;
  vector[N] y;

  int<lower=0> nnz;
  vector[nnz] X_w;          // non-zero values
  array[nnz] int X_v;       // column indices
  array[N + 1] int X_u;     // row-start indices
}

parameters {
  real alpha;
  vector[K] beta;
  real<lower=0> sigma;
}

transformed parameters {
  vector[N] mu;

  mu = alpha
       + csr_matrix_times_vector(
           N,
           K,
           X_w,
           X_v,
           X_u,
           beta
         );
}

model {
  alpha ~ normal(0, 2);
  beta ~ normal(0, 1);
  sigma ~ normal(0, 1);

  y ~ normal(mu, sigma);
}