MODULE SPARSEM

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 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 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.