Table of Contents
Karin's FORTRAN gems
This page gathers some small programs or routines, with limited documentation.
If you like it, cite it!
B-splines
Simple program to generate B-splines functions for given range or given values of covariable, with choice of equi-distant splines or user-defined knots.
- Source code Bsplines.f90
- Executables bsplines.zip
Generalised inverse of symmetric matrix
Fortran 90 rewrite of subroutine DKMWHF (derived from Henderson's DJNVHF) to obtain a generalised inverse of a real, symmetric matrix.
Simple front-end program to DKMXHF which reads in the matrix from a file & writes out the inverse to another file (DKMXHF not included).
Complete program for downloading (ginverse.f90 with dkmxhf.f90) ginverse.zip
Hash tables
Hash tables are a standard computational tool to facilitate fast `look-up' in a table or structured list.
A common task in analyzing quantitative genetic data is to identify the unique levels of an effect (with many different levels), e.g. animal IDs, in large data sets. Hash tables provide an effective way to carry out this task, and recode levels in running order (1 to N) when no numerical ordering of the original codes is required.
Here is a simple FORTRAN 95 module (derived from original code by Bruce Tier) which provides all the necessary routines to use one or more hash lists for such task:
And here is a small example program illustrating how to use it:
Download both: hash.zip
Approximate inverse of a sparse matrix
An approximate inverse of a sparse (symmetric and positive definite) matrix – or selected elements thereof - can be obtained using Gibbs sampling, as described by Harville (1999).
Here is a FORTRAN 95 module which contains the necessary subroutines to store the non-zero elements of the sparse matrix in linked list format, sort it and obtain the approximate inverse. The version given determines the diagonal elements of the inverse only, but is readily extended to return diagonal blocks or a complete inverse. It uses routines from the ranlibf
package (Texas Medical Center) to sample random normal deviates (not included; other RNG software can readily be substituted):
And here is a small test program illustrating its use:
Download both: harville.zip
@ARTICLE{harv99, author = {Harville, D. A.}, title = {Use of the Gibbs sampler to invert large, possibly sparse, positive definite matrices}, journal = {Lin. Alg. Appl.}, volume = {289}, year = {1999}, pages = {203--224}}
Updating the inverse of a dense matrix
"Super-nodal" sparse matrix factorisation and inversion
For symmetric, positive definite matrices stored in compressed format (see George and Liu, 1981) & assuming that matrices have been re-ordered to reduce fill-in and the symbolic factorisation has been carried out.
Computations require BLAS and LAPACK libraries. Fortran subroutines:
- Routine “super_nodes” works out the node boundaries and number of nodes, starting from the last row
- Routine “super_gsfct” performs the Cholesky factorisation; it replaces George & Liu's routine GSFCT, but uses a right- rather than left-looking algorithm
- Routine “super_sparsinv” takes the Cholesky factor as input and returns the subset of elements of the sparse inverse corresponding to non-zero elements in the Cholesky factor (Only those! This means that the product of the `inverse' obtained & the original matrix usually is not quite an identity matrix). It uses a block-wise adaptation of Takahashi's method.