====== 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 **[[BasicSplines|Bsplines.f90]]** 
    * Executables **{{:fortran:bsplines.zip|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.

**[[GenInverse|dkmxhf.f90]]**

Simple front-end program to DKMXHF which reads in the matrix from a file & writes out the inverse to another file
(DKMXHF not included).

**[[FrontGenInv|ginverse.f90]]**

Complete program for downloading (ginverse.f90 with dkmxhf.f90) {{:fortran:ginverse.zip|}}

===== Hash tables =====

[[wp>Hash_table| 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:

 **[[fortran:hashtabs|hash.f90]]**

And here is a small example program illustrating how to use it:

[[fortran:hashfront|hash_example.f90]]

Download both: {{:fortran:hash.zip|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): 

**[[fortran:harville|harville.f90]]**

And here is a small test program illustrating its use:

[[fortran:harvilletest|harville_test.f90]]

Download both: {{:fortran:harville.zip|harville.zip}}

<code bibtex>
@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}}
</code>

===== Updating the inverse of a dense matrix =====

[[http://didgeridoo.une.edu.au/km/tawny.php|Separate page for "TAWNY"]]


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



{{:fortran:super.tar.gz|}}