# SpMV solvers I

Back to reservoir simulation which I’ve enjoyed previously, lots of physics and numerics. The heart of the problem is the linear solver; after discretization of partial differential equation & linearization, one has a rather neat looking matrix equation A.x = b, where A are the coefficients, b the RHS vector and x the unknowns. Rather trivial for a small non-singular matrix, for a practical system composed of >1e7 non-zero elements, and matrices of height ~ 1e6 elements, one is faced with using a direct or iterative solution, finding inv(A) is just impossible. Consideration must also be given to storage format (eg., compressed sparse row) and memory usage; direct methods tend to blow up (like n^3) due to fill rates when factoring, and for other reasons. Iterative methods scale better but are more susceptible to numerical problems, and are a strong function of the matrix structure. As I trawl through the literature and iterate through possible solvers, a few collections standout eg., LIS, ILUPACK, and those in the intel MKL. Octave has useful formats, functions and solvers, including a nice polymorphic (adaptive) direct solver. What follows is a little bit ‘o’ script for converting output from the simulation, converting from CSR to octave’s format. Too bad the direct method doesn’t scale well, this worked nicely…

```
# output
fid =fopen('matrix_data.txt','r');
# unpack CSR details; output
[cA,b] =fscanf(fid,"%f,",Inf);
[rA,b] =fscanf(fid,"%f,",Inf);
[dA,b] =fscanf(fid,"%e,",Inf);
[rhs,b] =fscanf(fid,"%e,",Inf);
fclose(fid);
# matrix reconstruction A
rInd =zeros(size(cA));
rInd(rA(1:end-1),:) = 1;
#skip zero rows
rInd(rA(1:end-1),:) +=((rA(2:end)-rA(1:end-1))==0);
rInd = cumsum(rInd);
#store in octave sparse
sp = sparse(rInd, cA, dA);
#rhs vector
b=rhs;
#direct method
tic; x = b' / sp; timeA = toc;
#lu decomp for precond
tic; [l,u,p,q]=lu(sp); timeB=toc;
spy(u);
```