*Peter Briggs
CCP4 Daresbury Laboratory*

As of version 4.1, CCP4 now includes the option of using the LAPACK linear algebra package. This article is intended to give some background information about LAPACK and the functionality it provides, as well some details of how CCP4 can be configured at installation time to use it. Links to relevant documentation are given, and a simple example using LAPACK in the CCP4 program MLPHARE is provided as a practical demonstration.

LAPACK stands for **L**inear **A**lgebra **Pack**age
[1], and is a collection of standardised subroutines for
various mathematical operations, for example solving systems of linear
equations or eigenvector problems. Some systems or system configurations
include libraries which incorporate the LAPACK routines; for those systems
which don't, the nessary source code can be obtain from the NetLib archive
[2].

LAPACK is based on the BLAS (**B**asic **L**inear **Al**gebra
**S**ubprograms) [3], another set of standardised
subroutines and functions which deal with vector and matrix operations.
Again, many systems provide a vendor version of the BLAS routines; if not then
a ``reference'' BLAS can be obtained from NetLib. Since the vendor BLAS have
been optimised for a specific system these are expected to be more
efficient than the reference version. On systems without a vendor BLAS (or
for those people not wishing for whatever reason to use vendor BLAS) freely
available projects such as ATLAS [4] can be used to
optimise BLAS performance.

LAPACK is a freely-available software package which can be incorporated into other software packages (including commercial ones) provided that proper credit is given to the authors. To reference LAPACK in a scientific publication you should cite the LAPACK Users' Guide [1]. For more details it is recommended that you consult the LAPACK FAQ [5].

As of version 4.1, CCP4 now includes the option of linking in the
LAPACK libraries. This has to be done at install time, by running the
main configure with an extra option **--with-lapack**.

This latest release includes the source code for both LAPACK 3.0 and the reference BLAS, both obtained from NetLib - but configure should first search for ``native'' LAPACK libraries on your system, and will use these preferentially if found. If it cannot find a LAPACK library then it will search for native/vendor BLAS libraries and use these in preference to the reference BLAS, since they should be more efficient. Only if neither are found will the reference BLAS be built.

This behaviour can be over-ridden using the **--with-lapack=FORCE**
option with configure. This forces building of the reference BLAS and
LAPACK routines. Using the reference BLAS can be expected to result in
poorer performance (in terms of speed) than when using vendor BLAS.

More information can be found in the CCP4 installation document [6] and in the MODLIB documentation [7] (which includes links to vendor BLAS and LAPACK libraries for some systems).

It should be noted that LAPACK provides routines for real, double
precision, complex and double complex variables, and while the source code
for all four precisions is supplied with CCP4 4.1, only real, double and
complex are actually built under the **--with-lapack** option. For
non-IEEE compliant machines there is one further issue which is addressed
below.

Once CCP4 has been configured to use the LAPACK libraries, the
`XLAPACK_LIB` variable in the Makefiles should show how to reference
them for linking purposes - for example:

XLAPACK_LIB = -L/ccpdisk/xtal/ccp4-4.1/lib/lapack -llapack -L/usr/lib -lblas

None of the existing CCP4 libraries or programs currently use LAPACK - its inclusion at this point is to allow developers to take advantage of the routines in their programs in future. It is expected that future releases of the suite will use LAPACK routines, for example to replace some of the existing MODLIB routines, and at this point LAPACK will become a standard part of the CCP4 installation.

LAPACK is capable of solving systems of linear equations, linear least squares problems, eigenvalue problems and singular value problems. It is also able to handle many associated computations, such as matrix factorizations or estimating condition numbers.

The routines themselves are divided into three sets:

**Driver routines**for solving standard types of problem**Computational routines**for distinct computational tasks**Auxiliary routines**to perform subtasks or low-level computations.

Typically driver routines will call a sequence of computational routines, and for standard problems it is likely that there will be a suitable driver routine which can be used directly without needing to deal with the computational routines underneath.

There are various ways of quickly identifying a suitable driver routine for a particular problem. I would recommend the LAPACK User Guide [1] as the best starting point:

- The user guide contents can be used to find lists of suitable routines for each problem quickly and easily; see http://www.netlib.org/lapack/lug/node1.html
- The user guide also contains a combined index of the driver and computational routines at http://www.netlib.org/lapack/lug/node142.html, which can be used as a quick reference.

However the User Guide appears to stop short of giving a full specification
of the subroutine arguments, so at this stage it is necessary to refer to the
comments in the source code for the specific routine, for example by looking in
the directory *$CCP4/lib/lapack/src/*, or by directly accessing the source
code on the web (which can be done using the links from the LAPACK FAQ at
http://www.netlib.org/lapack/faq.html#1.15.)

It may be that you have a particular problem which is not addressed by an existing driver routine. It should be noted however that taken together the computational routines can perform a much wider range of tasks than represented by the driver routines alone - so in these cases it should be possible to create your own "custom" driver. The computational routines are documented with the same level of detail as the drivers, and this documentation can be accessed in the same ways as outlined above.

In principle many programs in the CCP4 suite contain subroutines which could be replaced by calls to LAPACK routines; for example MLPHARE (heavy atom refinement and phase calculation program) uses the following:

- EIGN1M: find the eigenvalues and vectors of a real symmetric matrix
- MATSOL: uses a modified Cholesky method to solve the matrix system
__A__**x**=**b** - MATIN1: matrix inversion with accompanying solution of linear equations

In practice replacing these routines means creating wrappers to prepare input for the appropriate LAPACK routine(s) and then performing any necessary conversion of the LAPACK output into the form expected by the calling subprogram. In many cases this is likely to result in a disproportionate amount of additional code.

For demonstration purposes however a version of MLPHARE was prepared in which the subroutine EIGN1M was replaced by a new routine using the LAPACK routine DSYEVR to perform the same task. DSYEVR computes selected eigenvalues and, optionally, eigenvectors of a real symmetric matrix using "Reasonably Robust Representations" (RRR).

The full code of the altered program is not reproduced here, but a patch with the replacement routine EIGN1M_W_LAPACK can be accessed at http://www.ccp4.ac.uk/newsletter39/04_mlphare_patch.f. It should be possible to compile and run a version of MLPHARE with this patch under CCP4 4.1, provided that it has been configured to use LAPACK as described previously (though see the comments below regarding installation on non-IEEE compliant machines).

Eigenvalues and eigenvectors of the B-factor matrix are only calculated in MLPHARE when anisotropic B factors are supplied. Comparison of the output from the two versions in all the test cases gave the same sets of eigenvalues. For the simple case where the anisotropic B factor coefficients are estimated from the isotropic value <B> using:

B_{1} = B_{4} = B_{6} = <B>

B_{2} = B_{3} = B_{5} = 0.0.

the eigenvalues are three-fold degenerate and equal to <B> - so any set of three linearly independent vectors forms a valid set of eigenvectors. In this case the sets of eigenvectors from the two versions differed simply as a consequence of using different algorithms.

Once the anisotropic B estimate is perturbed slightly (for example by
setting B_{3} = 1.0) to break the degeneracy, then the two
versions reassuringly give identical sets of eigenvalues and vectors.

There is at least one installation issue which isn't yet automatically resolved by configure: the routine ILAENV, which tests for (amongst other things) IEEE-754 compliance for Nan and infinity arithmetic. If LAPACK is being installed on a non-IEEE compliant machine then these tests will fail at run-time and cause a program crash.

Once it has been established that you are installing on a non-IEEE
compliant machine then the following change needs to be made to the
ILAENV source code, *$CCP4/lib/lapack/src/ilaenv.f*:

diff ilaenv-dist.f ilaenv-non-ieee.f 527,528c527,528 < C ILAENV = 0 < ILAENV = 1 --- > ILAENV = 0 > C ILAENV = 1 538,539c538,539 < C ILAENV = 0 < ILAENV = 1 --- > ILAENV = 0 > C ILAENV = 1

This is necessary at least for on a DEC Alpha running OSF1 v4.0 - for instance, on this platform the call to DSYEVR fails in the MLPHARE example above unless the change is made.

Note that this is not strictly a bug: more background information about this issue can be found in section 6.1.4 of the LAPACK "Quick Installation Guide for Unix Systems" [8]; in future the CCP4 install should take care of this problem automatically.

As of CCP4 4.1 it possible to make use of the LAPACK linear algebra
library as part of the suite. To enable LAPACK within CCP4, run the main
configure script with the **--with-lapack** option.

LAPACK offers a robust and comprehensive set of routines which can be used to solve a number of standard numerical problems. In future it is envisaged that LAPACK will be a standard part of the CCP4 installation, and it is hoped that developers will take advantage of the routines to develop new code more quickly.

Finally, the LAPACK installation under CCP4 is still relatively untested - so I would be grateful to hear from anyone who can give me feedback (good or bad) about installing and using the routines. Please e-mail me at p.j.briggs@ccp4.dl.ac.uk if you have any comments.

[1] Anderson, E. and Bai, Z. and
Bischof, C. and Blackford, S. and Demmel, J. and Dongarra, J. and Du Croz,
J. and Greenbaum, A. and Hammarling, S. and McKenney, A. and Sorensen, D.,
*LAPACK Users' Guide, Third Edition, 1999, pub. Society for
Industrial and Applied Mathematics, Philadelphia, PA, ISBN 0-89871-447-8
(paperback)
A HTML version of the LAPACK 3.0 Users Guide can be found at http://www.netlib.org/lapack/lug/lapack_lug.html
*

*
[2] Netlib web address: http://www.netlib.org
*

*
[3] BLAS: http://www.netlib.org/blas/
*

*
[4] Automatically Tuned Linear Algebra Software (ATLAS): http://www.netlib.org/atlas/
*

*
[5] LAPACK FAQ: http://www.netlib.org/lapack/faq.html
*

*
[6] CCP4 installation document: http://www.ccp4.ac.uk/dist/INSTALL.html), section E "Configure Options"
The local version is in $CCP4/INSTALL.html.
*

*
[7] CCP4 MODLIB documentation: http://www.ccp4.ac.uk/dist/html/modlib.html), section 2 "Information on BLAS and LAPACK routines".
The local version is in $CCP4/html/modlib.html.
*

*
[8] LAPACK Working Note 81: Quick Installation Guide for LAPACK on Unix Systems
can be downloaded as a postscript file from
http://www.netlib.org/lapack/lawns/lawn81.ps, but it should be noted that the routines distributed with CCP4 4.1 form only a subset of the full LAPACK distribution available from NetLib.
The full set of working notes are available via http://www.netlib.org/lapack/lawns/index.html
*

*
*