title: Solving large system of coupled differential equations
body: <p>I have a system of coupled ordinary differential equations<br>
  <code>dx/dt = (A + C_d(t) * B) * x</code>,<br>
  where A and B are constant matrices and C_d is a diagonal coefficient matrix which smoothly varies depending on the current value of the integration variable. </p>
  
  <p>The square matrices A and B are built up from smaller 60*60 upper triangular or zero matrices. The dimension of the full system is around 2500*2500. A and B are sparse with ~10% non-zero elements. The diagonal elements are negative or zero. The main (physical) constraint is that elements of x(t) are not allowed to become negative during integration.</p>
  
  <p>Currently, I employ a ‘naïve’ step solver<br>
  <code>x_(i+1) = A * x_i * dt_i  + B * (C_d(t_i) * x_i) * dt_i + x_i</code><br>
  or in the CPU/GPU versions</p>
  
  <pre><code>def solve_CPU(nsteps, dt, c_d, x):
      for step in xrange(nsteps): 
          x += (A.dot(x) + B.dot(x * c_d[step])) * dt[step]
  
  def solve_GPU(m, n, nsteps, dt, c_d, cu_curr_x, cu_delta_x, cu_A, cu_B):
      for step in xrange(nsteps): 
          cubl.gemv(trans='T', m=m, n=n, alpha=1.0, A=cu_A,
                    x=cu_curr_x, beta=0.0, y=cu_delta_x)
          cubl.gemv(trans='T', m=m, n=n, alpha=c_d[step], A=cu_B,
                    x=cu_curr_x, beta=1.0, y=cu_delta_x)
          cubl.axpy(alpha=dt[step], x=cu_delta_x, y=cu_curr_x)
  </code></pre>
  
  <p>and make use of a feature, that the step sizes <code>dt_i</code>this can be computed <em>a priory</em> in a way that the elements of x are always >=0 during integration. Depending on the amount of approximations and the settings the number of integration steps varies between 25k and 10M. </p>
  
  <p>I  have tried several methods to optimize performance on general purpose hardware:</p>
  
  <ul>
  <li>(unknown) When using ODEPACK’s VODE solver, I do not know how to express the x>=0 constraint</li>
  <li>(slowest) Dense BLAS 2 dot-product using Intel MKL</li>
  <li>(medium) Dense BLAS using single precision cuBLAS on NVIDIA GPU</li>
  <li>(fastest) SCIPY sparse module using CSR/CSC formats</li>
  </ul>
  
  <p>The code is written in Python and has access to the above listed libraries via Anaconda, Numba, Accelerate, Numpy etc. SCIPY's sparse BLAS routines are not properly linked to MKL in Anaconda and Python wrappers around cuSPARSE are to my knowledge not available, yet. I would know how to squeeze out a little bit more performance by directly interfacing to cuSPARSE/C-MKL sparse dot product, but that’s it. This exercise has to be solved dozens of times, again and again if models change, so performance is always an issue. I’m not an expert in this matter, so I don’t know much about preconditioners, factorization theorems etc. what brings me to my question:</p>
  
  <p>Is there a more elegant or better way to solve such a linear-algebra task?</p>