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 reallocated^{1} 
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)=1e10, 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 