Future Plans for MODLIB

Peter Briggs
January 2000

Contents of this document:

Warning: this is an evolving document - the work is in progress! Things may change - apologies for any inconsistencies.

• Background
• Plan for rationalisation
• Mathematical routines in the CCP4 Suite
  • Comments and suggested plan of action
• Using BLAS and LAPACK
  • Overview
  • Vendor BLAS for various systems
  • Implementation within CCP4
• Appendix A: Mathematical subroutines in CCP4 Supported Programs
• Appendix B: Mathematical subroutines in CCP4 Libraries

Background

MODLIB is a CCP4 subroutine library containing a motley collection of mathematical-type routines. The operations catered for are rather arbitrary. Furthermore, not all the mathematical routines present in the suite are collected in MODLIB - some are in other libraries (for example, in SYMLIB), and there are any number of other routines in various programs.

Plan for rationalisation

There are four basic stages in the rationalisation programme:

1. Document the contents of MODLIB

   This has now been done; the resulting document is available at http://www.dl.ac.uk/CCP/CCP4/ccp4/html/modlib.html, and will be also be included as part of future CCP4 releases (v4.0 onwards).

2. Compile a list of mathematical subroutines in the suite

   This will mean sorting through the programs and other libraries and listing the ``mathematical'' subroutines.

   Progress with this stage is below

3. Rationalise MODLIB and the suite

   This stage follows on directly from (2), and would consist in principle of two sub-stages:
   1. Replacing subroutines in programs with calls to equivalent routines in MODLIB,
   2. Transferring useful routines from programs or other libraries into MODLIB, and removing where possible duplicated functionality.

   It is also possible at this stage that omissions might be identified in the available routines in the library, or that a single generalised routine could replace a set of specific versions.

   Comments and a plan of action is below

4. Replace routines in MODLIB with optimised versions from BLAS etc

   BLAS and LAPACK are freely-available libraries of mathematical subroutines, and versions optimised for specific systems can be obtained. It is possible that other external libraries might also be used.

   BLAS and LAPACK are now available within CCP4 (v4.1 onwards). Some comments and ideas about BLAS and LAPACK are given below.

Mathematical Routines in the CCP4 Suite

This is work in progress. Appendix A lists the routines identified in various programs as possible candidates for rationalisation, either by relocation to the library, or by replacement with calls to library routines.

Programs examined so far:

areaimol, act, almn, absurd, baverage, bones2pdb, bplot, cad, contact, coordconv, crossec, detwin, distang, dm, dmmulti, dyndom, ecalc, f2mtz, fft, fhscal, findncs, freerflag, gensym, geomcalc, getax, hgen, hklplot, icoefl, lsqkab, makedict, mama2ccp4, mapdump, mapmask, maprot, mapsig, maptona4, matthews_coef, mlphare, molrep, mtz2various, mtzMADmod, mtzdump, mtzmnf, mtztona4, mtzutils, ncsmask, npo, oasis, omit, overlapmap, pdbset, peakmax, phistats, polarrfn, polypose, protin, rantan, rebatch, refmac, reindex, restrain, revise, rfcorr, rotmat, rsearch, rsps, rstats, rwcontents, scaleit, scalepack2mtz, sc, scala, sfall, sigmaa, sortmtz, sortwater, stnet, surface, tffc, tracer, truncate, undupl, unique, vecref, vectors, volume, watertidy, watncs, watpeak, wilson

(Programs in bold have routines which should be examined for rationalisation.)

Details of routines in the libraries can be found in Appendix B.

Comments and suggested plan of action

This exercise has revealed a substantial amount of duplication of functionality, and in many cases duplication of actual code, in routines found in different programs.

It is also apparent that there is a deficiency in the range of functions offered by MODLIB, in particular routines capable of dealing with:

• 4-matrices and vectors,
• eigenvalue and eigenvector problems,
• sorting routines, and
• generating rotation matrices from various conventions, and extracting rotation angles from existing matrices.

My suggested next set of steps is as follows:

1. Identify duplicated code and where possible bring these routines into the library.
   Possible candidates include GMPRD, DCOSFD, MATMUL4.

2. Check for duplication of functionality between routines in different programs, and between those in the programs and in the library.
   Replace routines in programs with calls to equivalent routines in MODLIB (if routines with the required functionality don't exist in the library, transfer them from the programs first).

3. Check for duplication of functionality between routines within the library, (for example GMPRD probably duplicates an existing library routine) and rationalise.
   Essentially this should be a cleaning-up exercise after the previous steps.

4. Identify routines which are not duplicated, but which seem to have useful functionality which is not covered by the library, and bring these into MODLIB too. A program like rsps has its own "library" which might be a useful model for identifying gaps in MODLIB at this stage.

It will not be appropriate to do this for every single routine or program; in some cases it will left up to the program authors/maintainers to decide whether to make the switch. After all, this job will be hard enough as it is!

Using BLAS and LAPACK

NOTE: since the time of writing CCP4 4.1 has been released, which has the option of linking in BLAS and LAPACK at install-time, and includes the source code for BLAS and LAPACK taken from NetLib.
Some of the comments below are therefore now out-of-date; furthermore information about LAPACK within CCP4 has become rather dispersed:
BLAS and LAPACK in CCP4 (from the MODLIB documentation).
Configure options from the current CCP4 INSTALL document (4.1 at the time of writing) - see the --with-lapack option.
LAPACK in CCP4 4.1 from CCP4 newsletter 39 (may not be available yet).

Overview

This would be the last step in rationalisation, and consists of replacing CCP4 code with calls to external libraries to perform various maths functions. The motivation for doing this is that we can access routines which are robust and efficient.

BLAS (Basic Linear Algebra Subprograms) contains routines for "basic" vector and matrix operations, for example dot product. Architecture-specific optimised versions of BLAS are available for different systems (see list). Generic (non-optimised) versions of BLAS are also available from NetLib, though these will be slower (and in some cases may give significantly worse performance).

LAPACK (Linear Algebra PACKage) contains routines for: solving systems of linear equations, linear least squares problems, eigenvalue problems and singular value problems. LAPACK can also handle many associated computations such as matrix factorizations or estimating condition numbers.

The complete LAPACK can be freely obtained from NetLib. The routines are written so that as much as possible computations are performed by calls to the BLAS. Having optimised BLAS routines is therefore an advantage speed-wise; the generic version of BLAS allows LAPACK to be run on machines which do not have any other implementation of the BLAS. Some vendors provide maths libraries which include BLAS and LAPACK.

For more information about the BLAS, visit http://www.netlib.org/blas/faq.html; this site also offers links to vendor-supplied BLAS (though some of these seem out of date now).

For information about LAPACK see http://www.netlib.org/lapack/faq.html.

Vendor BLAS for various systems

The BLAS libraries are provided in a number of different ways for different platforms, and it may be difficult to rationalise.

System	Name of library	Link	Comments
Digital/Compq	CXML (Compaq eXtended Math Library), libcxml	http://www.compaq.com/math/	Freely downloadable, includes BLAS 1-3, LAPACK + others.
SGI	SCSL	http://www.sgi.com/software/scsl.html	Freely available (only for IRIX 6.4 onwards however), includes BLAS and LAPACK.
HP	MLIB (Mathemetical Software Library)	http://www.hp.com/rsn/mlib/mlibhome.html	Free evaluation copy available, otherwise must be paid for(?), includes BLAS and LAPACK.
NT	Intel Math Kernal Library	http://developer.intel.com/vtune/perflibst/mkl/	Licensed binaries are freely redistributable, includes BLAS and a "substantial subset" of LAPACK.
Sun	?	Try http://www.sun.com/workshop/fortran/	Possibly part of Forte?
Linux	Intel BLAS	http://www.cs.utk.edu/~ghenry/distrib/	BLAS only.

Implementation within CCP4

It will necessary for those wishing to use the BLAS etc routines to obtain the libraries themselves; they would not be supplied with the suite. This shouldn't be a problem since most vendors seem to supply a version of BLAS, either as part of the standard compilers or else freely downloadable, whilst Netlib also provides an (un-optimised) version. Netlib also provides the LAPACK routines.

A clone of MODLIB would be written which contains the same subroutines and functions as the CCP4 version, but replaces the code within those subroutines with calls to the BLAS etc libraries instead. A switch within configure (eg. --with-blas) would be necessary to tell the compilation which version of MODLIB it should use.

This stage is still some way off.

Appendix A: Mathematical subroutines in CCP4 Supported Programs

I have decided to organise the subroutines into groups based on similar/related functions, and within these groups to classify the routines with codes which identify those with (near-)identical functionality.

The groups are:

• Matrix operations
• Vector operations
• Rotations etc
• Sorting routines
• Eigen values/vectors
• Testing/setting matrices
• General array & matrix manipulations
• Solving sets of linear equations
• Miscellanous routines

Matrix Operations

The codes in this group are:

     MM  = matrix-matrix multiplication
     MV  = matrix-vector multiplication
     IM  = matrix inversion
     DET = calculate matrix determinant

Vector Operations

The codes in this group are:

     VNORM  = normalise vector (i.e. reduce to unit vector)
     VDOT   = dot product
     VCROSS = cross product
     VDIFF  = difference between two vectors
     VSUM   = sum of two vectors
     VTRIP  = vector triple product
     VPERM  = permute vector elements
     VUNIT  = set to unit vector (same as VNORM??)

     ANG2MAT = rotation matrix from angles
     MAT2ANG = angles from rotation matrix
     GETDIR  = direction cosines from angles etc
     ANGDIR  = angles from direction cosines etc

Sorting routines

The codes in this group are:

     SORT = sorting routine (details vary)

Eigenvalues/eigenvector routines

The codes in this group are:

    EIGEN = obtain eigen values and/or vectors

Testing/setting matrices

The codes in this group are:

    SETIDENT  = set matrix to identity
    TESTIDENT = test if matrix is identity
    TESTDET   =	test determinant of matrix
    TESTROT   = test if matrix is a rotation
    TESTORTH  = test if rotation is orthonormal

General array & matrix manipulations

The codes in this group are:

    COPY  = copy one array into another
    ZERO  = set an array to zero
    STATS = get statistics (e.g. mean) from an array of data

Solving sets of linear equations

The codes in this group are:

    SOLVE = solve sets of linear equations

Miscellaneous Routines

These are the ones left over, for the time being. Some of the codes are:

    HASH = hashing routines
    INTERP = interpolation of functions
    MTRANS = transform matrices
    ANGLE  = get angles from sets of vectors/points/etc

Other programs

These are (generally) the programs which are made up from several files. Often they use internal routines with similar functions to those listed above but they have standard headers peculiar to that particular program, which makes them unsuitable at present for inclusion in the library. They are listed here for the sake of completeness.

dm, dmmulti, mapmask, maprot, ncsmask

These programs contain a number of generalised routines for manipulating vectors/matrices (symuvw, symmatmul etc) and for generating rotation matrices (eulertomatrix, polartomatrix).

molrep

Contains matrix operations, rotation-matrix generation, sorting, finding eigenvalues and eigenvectors (amongst others).

refmac

refmac has its own set of routines for matrix operations etc in the linalgebra.f file. This also contains routines for linear equation solution by Gauss-Jordan method (``GAUSSJ'') and solving the normal equation (a*dv=v) by conjugate gradient method (``CGSOLVE'').

restrain

restrain contains a number of routines for the solution of linear equations (some derived from LINPACK?).

Subroutines	Description
congra	This subroutine solves the normal equations by the Hestenes & Stiefel Conjugate Gradient method.
detmat	THIS SUBROUTINE CHECKS DETERMINANT OF 3*3 ROTATION MATRIX
orthog	THIS SUBROUTINE CHECKS IF 3*3 ROTATION MATRIX IS ORTHOGONAL
renorm	THIS SUBROUTINE RENORMALISES 3*3 ORTHOGONAL ROTATION MATRIX
eigen	SUBROUTINE TO COMPUTE EIGENVALUES & EIGENVECTORS OF A REAL SYMMETRIC MATRIX, FROM IBM SSP MANUAL [nb copy of eigen found in other progs?]
matpol	Given rotation matrix R, returns direction D of rotation axis, and polar angles A (deg.).
invmat	Invert 3x3 matrix.
condir	Given a starting point XP that is a vector of length NP, Fletcher-Reeves-Polak-Ribiere conjugate direction minimisation is performed on a function defined by DELVEC(x), using its gradient calculated by SCGRAD(x,g)
linmin	Minimisation of FX=FUNC(X) by line search along X
EIGSOL	Solution of linear equations by eigenvalue filtering.
solve	SOLVE NORMAL EQUATIONS BY CHOLESKY FACTORISATION.
solved	(LINPACK) Single precision positive-definite packed real symmetric matrix. Factorization step.
solver	(LINPACK (SPPSL)) Single precision positive-definite packed real symmetric matrix. Solve linear equations step.
solvei	SOLVE NORMAL EQUATIONS BY CHOLESKY FACTORISATION - INVERSION STEP
deigen	SUBROUTINE TO COMPUTE EIGENVALUES & EIGENVECTORS OF A REAL SYMMETRIC MATRIX, FROM IBM SSP MANUAL (SEE P165).
saxpy	(LINPACK) Single precision constant*vector + vector.
MATRIX	A SUBROUTINE TO TRANSFORM A TENSOR (TEN) TO THE TENSOR (CAR) WRT TO ANOTHER SET OF AXES BY MULTIPLYING BY THE MATRIX OF EIGENVECTORS (VEC) AND ITS TRANSPOSE: CAR = VEC~.TEN.VEC .
MATMLT	A SUBROUTINE FOR MULTIPLYING TWO 3*3 MATRICES TOGETHER.
CALCRS	THIS SUBROUTINE CALCULATE EIGENVALUES (EVAL) AND VECTORS (EVEC) OF MATRIX R OF ORDER N
EIGEN	SUBROUTINE TO COMPUTE EIGENVALUES & EIGENVECTORS OF A REAL SYMMETRIC MATRIX, FROM IBM SSP MANUAL

rsps

Matrix and vector operations can be found in veclib.f. This appears to be quite a comprehensive set of routines:

*       s/r matcmp     compare two n x m real matrices
*       s/r rmatcmp    = matcmp
*       s/r imatcmp    compare two n x m integer matrices
*
*       s/r matcop     copy n x m real matrix
*       s/r rmatcop    = matcop
*       s/r imatcop    copy n x m integer matrix
*       s/r drmatcop   copy n x m double precision matrix to real
*
*       s/r matclr     clear n x m real matrix
*       s/r rmatclr    = matclr
*       s/r imatclr    clear n x m integer matrix
*       s/r dmatclr    clear n x m double precision matrix
*       s/r cmatclr    clear n x m character matrix
*
*       s/r rmatset    set all elements in n x m real matrix
*       s/r imatset    set all elements in n x m integer matrix
*       s/r dmatset    set all elements in n x m double precision matrix
*       s/r cmatset    set all elements in n x m character matrix
*
*       s/r vecdif     subtract two real vectors
*       s/r vecadd     add two real vectors
*       s/r vecabs     abs of real vector elements
*       s/r rvecdif    = vecdif
*       s/r rvecadd    = vecadd
*       s/r rvecdif    = vecdif
*       s/r rvecmod    modulus of real vector
*       s/r ivecdif    subtract two integer vectors
*       s/r ivecadd    add two integer vectors
*       s/r ivecabs    abs of integer vector elements
*
*       s/r ivecf      convert integer vector to real
*       s/r fveci      convert real vector to integer
*
*       s/r mtxtps     transpose n x m real matrix
*       s/r mtxinv     invert n x n matrix
*       s/r mtxmlt     multiply n x m matrix by m x n matrix
*       s/r vmxmlt     mulitply n x 1 vector by m x n matrix
*
*       s/r transpose  transpose n x n real matrix
*       s/r matinv     invert 3 x 3 real matrix
*       s/r vmtply     multiply 3 x 1 vector by 3 x 3 matrix
*       s/r mmtply     multiply 3 x 3 matrix by 3 x 3 matrix
*       s/r cosabg     cosine of three angles
*       s/r ortho      get (de)orthogoanlization matrix from cell
*       s/r pnt2pln    get plane equation from three points
*       fn  plndst     distance between point and plane
*       s/r polar2mtx  convert polar angles to matrix
*       s/r nfold2mtx  convert direction cosines and rot angle to matrix
*       s/r msk        .and. two logical vectors, count and point
*       s/r mskcmp     compare two logical vectors
*
*       s/r isort      sort real array, co-sort integer vector
*       s/r fsort      sort real array, co-sort real vector
*       s/r qsortfi    quicksort real array, co-sort integer vector
*       s/r qsortff    quicksort real array, co-sort real vector
*       s/r qsortffi   quicksort real array, co-sort real and integer vector
*       s/r qsortpf    quicksort pointers
*       s/r nsorta     Sort N columns on values in key column, acending
*
*       fun det        determinant of 3 x 3 matrix
*       fun vlen       length (A) of fractional vector in ]-0.5,0.5[
*       fun vlen2      length (A) of fractional vector
*
*       fun imaxvc     maximum element in vector (integer)
*       fun iminvc     minimum element in vector (integer)
*       fun fmaxvc     maximum element in vector (real)
*       fun fminvc     minimum element in vector (real)
*       fun random     pseudo-random number generator (real)

Appendix B: Mathematical subroutines in CCP4 Libraries