Imagine a typical BLUP program

real x(n,n)
...
x=0
...
x(i,j)=x(i,j)+a
...
call solve(x,n,sol,rhs)

This program is simple for dense matrices but its application is restricted to small problems. Define type of dense matrix in Fortran 90:

type sparsem
integer::n
real (8),pointer::x(:,:) !dimensions to be determined later
end type

and write subroutines operating on that structure so that the same program is:

type (densem):: x
...
call zero(x,n)
...
call addm(a,i,j,x)
...
call solve_iter(x,sol,rs)

Now define structures for sparse matrices and operations on them. Then the same program with hash sparse structure will look virtually identical:

type (sparse_hashm):: x
...
call zero(x,n)
...
call addm(a,i,j,x)
...
call solve_iter(x,sol,rs)

but it will handle millions of equations rather than thousands. Note that the same name applies to many structures, and that details (e.g., n) are hidden whenever possible. Familiar symbols can be used to reduce the learning time. To convert between any two formats, use the "=" operator, as in the example below:

type (sparse_hashm)::x !format good for set up
type (sparse_ija)::y !format good for factorization
...
y=x

Many more operations have been defined.
 

Operation	Description	Comments
call init(x)	Initialize x	Not necessary on systems where pointers are initialized automatically (e.g., Sun, Lahey F90)
call zerom(x,n)	Allocate storage for x as an n*n matrix and zero it	If x was set before, it is reallocated1
call reset(x)	Deallocates storage
call setm(a,i,j,x)	set:  x(i,j)=a	Does not work on SPARSE_IJA
call addm(a,i,j,x)	Add to matrix:  x(i,j)=x(i,j)+a	Does not work on SPARSE_IJA
a=getm(i,j,x)	get one element:  a=x(i,j)
x=y	Conversion between formats	Conversion from sparse to dense formats may require too much storage
call printm(x)	Prints x as square matrix	print(x,'internal') prints sparse matrices in internal format
call solve_iterm(x,rs,sol)	Solves: x sol=rs  iteratively by SOR
call default_iter  (conv,maxround,relax)	Changes default iteration parameters	All parameters are optional; default values are: convergence criterion)=1e-10, max round=1000, relaxation factor=1.0.
x=block(y,i1,i2,j1,j2)	Selects block from y:   x=y(i1:i1,j1:j2)	does not work on dense_symm format
tr=trace(x,y)	Self explanatory	real(8) function; x and y must be in same formats; works on densem and sparse_ija formats only;
q=quadrf(u,x)	q=u'Xu	real(8) function; does not work on dense_symm format

 
In the end, SPARSEM creates easy building blocks for BLUP , REML or MCMC programs. Such programs are being developed separately.