9/4/97 - 10/1/97
Traditionally, programming in animal breeding is done in 2 stages: in a matrix language and in a regular programming language. Programs in a matrix language such as IML SAS, Matlab, Mathematica or APL are reasonably simple and useful for creating examples but inefficient for large problems. Programs in a regular programming language such as Fortran or C/C++ are much more efficient but could take much longer to write and require substantial training.
Matrix languages are easy to deal with matrices partly because usually only one format is usually supported: dense rectangular. Operations on such matrices are easy to specify and program, but large matrices require large memory and long running time. Also, memory and computations are equal whether matrices are sparse (contain very few nonzero elements) or not. In animal breeding, many matrices are sparse. If that sparsity is taken into account, the memory requirements and computations can decrease dramatically. Unfortunately, there is more than one format for storing sparse matrices, and some computations are fast with one format and but not with another one. Also, the storage formats and operations are considerably more complicated than dense rectangular matrices. A library to handle multiple matrix formats and multiple operations would contain many subroutines, each with a long list of arguments. Such a library would involve considerable learning, and many details associated with the library would create many opportunities for making a mistake.
One matrix package, Matlab, has some forms of sparse-matrix storage and operations included.
Modern programming languages with "object-oriented" features, such as C++ or Fortran 90, have abilities to create classes/modules, where many implementation details on specific data structures can be hidden. A technique called overloading allows single function/subroutine to work with different formats of its arguments. Therefore, the number of details to remember can be drastically reduced. Subsequently, programming can be done much easier and quicker.
SPARSEM is a module for Fortran 90 that enables programming common sparse matrix operations almost as easily as with dense matrices. It supports two dense matrix formats, useful for testing, and two sparse matrix formats. Changing a program from dense- to sparse-matrix format using DENSEM can be as simple as changing one declaration line. SPARSEM incorporates an interface to FSPAK, which enables efficient sparse matrix factorization, solving, sparse inversion and calculation of determinant on matrices much larger than possible with dense matrix structures.
Four matrix formats are available.
DENSEM - dense square matrix.
DENSE_SYMM -dense symmetric upper-stored.
It has approximately only half memory requirements of the dense square matrix.
SPARSE_HASHM - sparse triple accessed by hash algorithm.
This is a very efficient format for set-up and for iterative-solving of sparse matrices.
SPARSE_IJA - Sparse IJA.
This is a memory-efficient format for sparse matrices used by sparsse matrix packages. Format IJA cannot easily be set up directly but can be derived by conversion from the hash format.
For more information on all these formats see Duff et al, George and Liu, or my class notes.
Currently, all the formats use double precision except SPARSE_HASM, which was chosen to use single precision for efficiency. real(4) except in SPARSE_IJA, where they are real(8).
A popular format that is not included here is linked list. That format is reasonably efficient for creating and computing with sparse matrices if the number of nonzero elements per row is not too high and the matrix is not too large. However, the combination of hash plus ija is generally more efficient.
The following subroutines/functions are supported. All real scalars and vectors are single precision unless indicated otherwise.
|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 addm(a,i,j,x)||Add to matrix:
|Does not work on SPARSE_IJA|
|call setm(a,i,j,x)||sets element of matrix:
|Does not work on SPARSE_IJA|
|y=getm(i,j,x)||find element of matrix:
|real(4) function; returns lower-diagonal elements of upper-stored matrix|
|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
|Changes default iteration parameters||All parameters are optional; default values are: conv(ergence criterion)=1e-10, max round(s)=1000, relax(ation factor)=1.0.|
|x=block(y,i1,i2,j1,j2)||Selects block from y:
|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|
All operations assume that the densem type is general while all the other types are upper-stored.
Operations tr, quadf work with both upper- and full-stored matrices but the block operation works literally, i.e., selecting a lower block of an upper-stored matrix would return an empty matrix and selecting an upper block would return only an upper-stored matrix. This could be a source of incompatibility between densem and other formats that use the block operation without taking its limitations into consideration. This problem is easy to notice in examples by printing matrices of interest.
FSPAK is a sparse matrix package written in F77 that performs operations on sparse matrices in format SPARSE_IJA. Operations include solving a system of linear equations by factorization, calculating a (log)determinant or finding a sparse inverse of a matrix. A sparse inverse is such a matrix that contains inverse values only for those elements that were nonzero in the original matrix. For sparse matrices, FSPAK is very efficient computationally. FSPAK90 is a F90 interface written to simplify the use of FSPAK.
A complete call to FSPAK90 is:
operation= "factorize" - calculate sparse factorization
"invert" - calculate sparse inverse
"solve" - solve a system of equation
"reset" - reset the storage
"det" - calculate determinant
"stat" - print statistics
ija = matrix in SPARSE_IJA form
rs = real (4) vector of right hand side
sol = real (4) vector of solutions
det = real(8) determinant or log-determinant
msglev= message level from 0 (minimum) to 3 (maximum); default=0
maxmem=maximum memory available in the system; default=infinite
rank=rank of matrix
rs8,sol8 = real(8) equivalents of rs and sol
All the arguments of fspak90 except "operation" and "ija" are optional except when they are needed in the specific "operation". Thus, rs and sol are needed for solving and det for "det" or "ldet".
All preceding steps are done automatically.
To solve using double precision right hand side and solutions:
To sparse invert:
To obtain the determinant d:
To obtain the log determinant ld:
To obtain rank r with any operation:
To force new factorization, when the input matrix has changed:
To deallocate the internal memory:
To limit memory to a maximum od maxmem, e.g., 20,000k, with any
Note that only relevant arguments for each step need to be included in calling FSPAK90. Reordering is performed the first time when FSPAK90 is called. Subsequent factorization except after the option "reset" will reuse the ordering. Subsequent solves will reuse the factorization.
Hints on using SPARSEM
Initially all the matrices can be implemented in DENSEM format. After the program works well with an example, convert all data structures for potentially large matrices to sparse formats and verify that same results are obtained.
Matrix types and functions subroutine are defined in module sparsem. Subroutine fspak90 is in module sparseop. In general, an arbitrary program xx.f90 can be compiled as on the Sun workstation as
f90 -Maa xx.f90 aa/sparsem.a
where aa is the directory containing the modules and the library. The conventions on compiling are slightly different in different systems.
SPARSEM has been successfully compiled under Unix on HP, RS/6000 (with
some modifications to makefiles), on SGI (with slight modification to code
to circumvent a compiler bug), under DOS using Lahey LF90 and under Linux
using Absoft f90.
Dense matrix solution program
integer :: i,j
real :: rs(n),sol(n)
! set up a sample matrix
call addm(10.0*i/j,i,j,x); call addm(10.0*i/j,j,i,x)
print*,'matrix' ; call printm(y)
call solve_iterm(y,rs,sol) !solve iteratively
Triangular dense matrix iterative-solution program
(The rest of the program remains identical)
Sparse hash matrix iterative-solution program
Sparse IJA matrix iterative-solution program
Matrix in ija form cannot be set up directly but can be converted from
call reset(x) ! Optional statement to release storage
print*,'matrix' ; call printm(y)
Sparse IJA matrix finite-solution and inversion program with FSPAK90
use sparseop !fspak90 is in module sparseop
George, A. and Liu, J.W.H. (1981) Computer solution of large sparse positive definite systems. Prentice-Hall, Englewood Cliffs, N.J.
Definitions of structures (types)
type densem !traditional dense square matrix
integer :: n
end type densem
type dense_symm !upper stored symmetric dense matrix
end type dense_symm
integer:: n,& ! for compatibility mainly
nel,& ! number of elements
filled,& ! number of filled elements
status ! 1 if ready to hash, 2 if in sorted
real , pointer :: x(:,:)
end type sparse_hashm
integer :: n,& ! number of equations
nel ! number of nonzeroes
integer, pointer::ia(:),ja(:) !will be ia(n+1), ja(m)
real (8), pointer::a(:) !will be a(m)
Structures can be accessed within the application program using the "%" symbol. This is useful, e.g., when using Fortran 77 programs. The example below shows how to use a determinant program written in F77.
type (densem):: z
! initialize z
call addm(i**j/10., i,j,z)
!calculate determinant for a 2x2 matrix
real :: x(n,n),det
The SPARSEM library includes the following files:
sparse.f90 - type definitions + main subroutines,
sparse2.f - supporting subroutines (in f77),
fspak.f90 - f90 interface to fspak
fspak.f - main fspak subroutine (in f77),
fspaksub.f - supporting fspak subroutines (in f77),
sparssub.f - low-level subroutines from the book of George and Liu (in f77),
second.f - timing subroutine specific to Sun (in f77).
Subroutines second() specific to other computers can be found in the FSPAK manual.