
MODLIB (CCP4: Library)NAMEmodlib  Subroutine Library for handling mathematical operationsDESCRIPTIONThe library modlib contains subroutines for a number of mathematical operations, mainly operations with vectors and matrices. This document is the first stage in an attempt to document and rationalise the mathematical subroutines found in the CCP4 libraries and programs. CONTENTSThe subroutines have been arranged as much as possible into groups with similar functions or origins.
Overview of the MODLIB routinesThe routines are divided into four sections below: vector operations, matrix operations, Harwell routines, and miscellaneous (i.e. those which don't fit in anywhere else). a. Vector Operations
b. MatrixVector Operations
Notes on MatrixVector operations: The subroutines MATVC4 and TRANSFRM perform the same operation, and differ only in that in the first case the result vector is returned separately from the input vector, whilst in the second the result overwrites the input vector. The transformation of the input 3vector by the 4x4 matrix is effected by treating the input vector as a 4vector, with 1 in the last position, and performing a standard matrix times vector operation. 4x4 matrices are commonly used in CCP4 programs to represent symmetry operations, which may have both rotational and translational components (screw operations). For more information on this, see the documentation for SYMLIB, in particular the section on How Symmetry Operations are Stored and Applied. c. Matrix Operations
d. Harwell SubroutinesThe Harwell Subroutine Library was started at Harwell by the precursor of the Numerical Analysis Group at Rutherford Appleton Laboratory in 1963, to provide high quality numerical software to the scientists and engineers working there. Since then it has developed into a portable fully documented Fortran library. See http://www.cse.clrc.ac.uk/Activity/HSL for more details. Some of the Harwell routines are already been distributed as part of CCP4, in modlib.f (ea06c, ea08c, ea09c, fa01as, fa01bs, fa01cs, fa01ds, fm02ad, mc04b) The modlib header has the following comment on the conditions of use of the routines: The routines ea06c, ea08c, ea09c, fa01as, fa01bs, fa01cs, fa01ds, fm02ad, mc04b, (and possibly others) are from the Harwell Subroutine library. The conditions on their external use, reproduced from the Harwell manual are: * due acknowledgement is made of the use of subroutines in any research publications resulting from their use. * the subroutines may be modified for use in research applications by external users. The nature of such modifiactions should be indicated in writing for information to the liaison officer. At no time however, shall the subroutines or modifications thereof become the property of the external user. The liaison officer for the library's external affairs is listed as: Mr. S. Marlow, Building 8.9, Harwell Laboratory, Didcot, Oxon OX11 0RA, UK. Information on some of the Harwell routines can also be found at http://www.netlib.org/harwell/index.html.
e. Miscellaneous
2. Information on BLAS and LAPACK routines(i) IntroductionFrom version 4.2, installation of the CCP4 suite is configured by default to include the BLAS and LAPACK libraries, which are used by some of the programs (e.g. SCALA, REFMAC5, BEAST). The advantage of using these external libraries is their robustness and efficiency of implementation, which should give greater reliability and higher speed. (ii) BLAS = Basic Linear Algebra SubprogramsBLAS contains routines for vectorvector, vectormatrix and matrixmatrix operations. Precompiled optimised versions of BLAS are distributed and should be used in preference to the versions available from NetLib, since they should be considerably faster. The NetLib BLAS FAQ is at http://www.netlib.org/blas/faq.html which includes a list of links to vendor BLAS libraries for various systems, and to the BLAS Quick Reference Guide. For systems where optimised ("vendor") BLAS libraries are unavailable, investigation of the Automatically Tuned Linear Algebra Software project (ATLAS) is recommended for improved BLAS performance see http://mathatlas.sourceforge.net/ for more information. ATLAS can be also be used on systems which have a vendor BLAS library provided. (iii) LAPACK = Linear Algebra PACKageLAPACK is built on the BLAS routines and provides routines for solving systems of simultaneous equations, leastsquares solutions of linear systems of equations, eigenvalue problems and singular value problems (plus others). The NetLib LAPACK FAQ is at http://www.netlib.org/lapack/faq.html There is also a html version of the LAPACK User Guide. (iv) BLAS and LAPACK in CCP4Installation of the CCP4 suite by default includes the BLAS and LAPACK libraries. This section describes how this done, and may be of interest if you have problems. LAPACK configuration can be disabled at build time by running the main CCP4 configure with the disablelapack option, but this is not advised as some programs will fail to build as a result. Many systems carry their own version of the BLAS and/or LAPACK libraries, and configure will try to use these preferentially. Some examples for common systems are given in the table:
If no system LAPACK library is found then configure will set up to build the `reference' version from Fortran source code originally taken from the Netlib archive. Source code for both BLAS and LAPACK v3.0 is included in the $CCP4/lib/lapack directory. Some comments:
Note that building of the Netlib routines can also be forced by specifying the withnetliblapack when running configure. Since the configuration and build has only been tested on a limited number of platforms, CCP4 welcomes any feedback or information about the LAPACK libraries under the Suite. If you have queries or comments please email then to ccp4@ccp4.ac.uk (v) Installation issues for nonIEEE compliant machinesBy default the NetLib version of the LAPACK routine ILAENV (distributed with CCP4) assumes an IEEEcompliant machine, which implies certain behaviour for NaN and "infinity arithmetic". On nonIEEE compliant machines or O/S the LAPACK library will compile successfully but may crash at runtime when ILAENV is called. The CCP4 distribution contains two versions of the ILAENV routine, in $CCP4/lib/lapack/src/. ilaenvdist.f is the standard Netlib version of the routine, whilst ilaenvnonieee.f is the version modified for a nonIEEE compliant machine. The CCP4 configure will copy the appropriate version of the routine to ilaenv.f on installation. 3. List of modlib routines:SUBROUTINE CROSS(A,B,C)Compute vector product A = B x C .. Array Arguments .. REAL A(3),B(3),C(3) REAL FUNCTION DOT(A,B)dot product of two vectors .. Array Arguments .. REAL A(3),B(3) SUBROUTINE EA06C(A,VALUE,VECTOR,M,IA,IV,W)C C** 18/03/70 LAST LIBRARY UPDATE C C ( Calls EA08C(W,W(M1),VALUE,VECTOR,M,IV,W(M+M1)) C and MC04B(A,W,W(M1),M,IA,W(M+M1)) ) C C Given a real MxM symmetric matrix A = {aij} this routine C finds all its eigenvalues (lambda)i i=1,2,.....,m and C eigenvectors xj i=1,2,...,m. i.e. finds the nontrivial C solutions of Ax=(lambda)x C C The matrix is reduced to tridiagonal form by applying C Householder transformations. The eigenvalue problem for C the reduced problem is then solved using the QR algorithm C by calling EA08C. C C Argument list C  C C IA (I) (integer) should be set to the first dimension C of the array A, i.e. if the allocation C for the array A was specified by C DIMENSION A(100,50) C then IA would be set to 100 C C M (I) (integer) should be set to the order m of the matrix C C IV (I) (integer) should be set to the first dimension C of the 2dimensional array VECTOR C C VECTOR(IV,M) (O) (real) 2dimensional array, with first dimension IV, C containing the eigenvectors. The components C of the eigenvector vector(i) corresponding C to the eigenvalue (lambda)i (in VALUE(I)) C are placed in VECTOR(J,I) J=1,2,...,M. C The eigenvectors are normalized so that C xT(i)x(i)=1 i=1,2,...,m. C C VALUE(M) (O) (real) array in which the routine puts C the eigenvalues (lambda)i, i=1,2,...,m. C These are not necessarily in any order. C C W (I) (real(*)) working array used by the routine for C work space. The dimension must be set C to at least 5*M. SUBROUTINE EA08C(A,B,VALUE,VEC,M,IV,W)C (Calls EA09C(A,B,W(M+1),M,W)) C C This uses QR iteration to find the all the eigenvalues and C eigenvectors of the real symmetric tridiagonal matrix C whose diagonal elements are A(i), i=1,M and offdiagonal C elements are B(i),i=2,M. The eigenvalues will have unit C length. The array W is used for workspace and must have C dimension at least 2*M. We treat VEC as if it had C dimensions (IV,M). C C First EA09, which uses the QR algorithm, is used to find C the eigenvalues; using these as shifts the QR algorithm is C again applied but now using the plane rotations to generate C the eigenvectors. Finally the eigenvalues are refined C by taking Rayleigh quotients of the vectors. C C Argument list C  C C A(M) (I) (real) Diagonal elements C C B(M) (I) (real) Offdiagonal elements C C IV (I) (integer) should be set to the first dimen C sion of the 2dimensional array VEC C C M (I) (integer) should be set to the order m of the C matrix C C VALUE(M) (O) (real) Eigenvalues C C VEC (O) (real) Eigenvectors. The dimensions C should be set to (IV,M). C C W(*) (I) (real) Working array.The dimension must be C set to at least 2*M. SUBROUTINE EA09C(A,B,VALUE,M,OFF)C 18/03/70 LAST LIBRARY UPDATE[No info] FUNCTION FA01AS(I)[No info]SUBROUTINE FA01BS(MAX,NRAND)[No info]SUBROUTINE FA01CS(IL,IR)[No info]SUBROUTINE FA01CS(IL,IR)[No info]SUBROUTINE FA01DS(IL,IR)[No info]DOUBLE PRECISION FUNCTION FM02AD(N,A,IA,B,IB)C Compute the inner product of two double precision real C vectors accumulating the result double precision, when the C elements of each vector are stored at some fixed displacement C from neighbouring elements. Given vectors A={a(j)}, C B={b(j)} of length N, evaluates w=a(j)b(j) summed over C j=1..N. Can be used to evaluate inner products involving C rows of multidimensional arrays. C It can be used as an alternative to the assembler version, C but note that it is likely to be significantly slower in execution. C C Argument list C  C C N (I) (integer) The length of the vectors (if N <= 0 FM02AD = 0) C C A (I) (double precision) The first vector C C IA (I) (integer) Subscript displacement between elements of A C C B (I) (double precision) The second vector C C IB (I) (integer) Subscript displacement between elements of B C C FM02AD the result SUBROUTINE ICROSS(A,B,C)C Cross product (integer version) C C A = B x C C C .. Array Arguments .. INTEGER A(3),B(3),C(3) INTEGER FUNCTION IDOT(A,B)C Dot product (integer version) C C IDOT = A . B C C .. Array Arguments .. INTEGER A(3),B(3) SUBROUTINE IMINV3(A,B,D)C Invert a general 3x3 matrix and return determinant in D C (integer version) C C A = (B)1 C C .. Scalar Arguments .. INTEGER D C .. C .. Array Arguments .. INTEGER A(3,3),B(3,3) SUBROUTINE MATMUL(A,B,C)C Multiply two 3x3 matrices C C A = BC C C .. Array Arguments .. REAL A(3,3),B(3,3),C(3,3) SUBROUTINE MATMULNM(N,M,A,B,C)C Multiply NxM MXN matrices C C A = BC C C .. Array Arguments .. REAL A(N,N),B(N,M),C(M,N) SUBROUTINE MATVEC(V,A,B)C Postmultiply a 3x3 matrix by a vector C C V = AB C C .. Array Arguments .. REAL A(3,3),B(3),V(3) SUBROUTINE MATMULGEN(Nb,Mbc,Nc,A,B,C)C Generalised matrix multiplication subroutine C Multiplies a NbxMbc matrix (B) by a MbcXNc C (C) matrix, so that C C A = BC C IMPLICIT NONE C .. C .. Scalar Arguments .. INTEGER Nb,Mbc,Nc C .. C .. Array Arguments .. REAL A(Nb,Nc),B(Nb,Mbc),C(Mbc,Nc) SUBROUTINE MC04B(A,ALPHA,BETA,M,IA,Q)C Transforms a real symmetric matrix A={a(i,j)}, i, j=1..IA C into a tridiagonal matrix having the same eigenvalues as A C using Householder's method. C C .. Scalar Arguments .. INTEGER IA,M C .. C .. Array Arguments .. REAL A(IA,1),ALPHA(1),BETA(1),Q(1) SUBROUTINE MINVN(A,N,D,L,M)C Purpose C ======= C C invert a matrix C C Usage C ====== C C CALL MINVN(A,N,D,L,M) C C Description of parameters C ========================= C C A  input matrix, destroyed in computation and replaced by C resultant inverse. C C N  order of matrix A C C D  resultant determinant C C L  work vector of length n C C M  work vector of length n C C Remarks C ======= C C Matrix a must be a general matrix C C Subroutines and function subprograms required C ============================================= C C NONE C C Method C ====== C C The standard gaussjordan method is used. the determinant C is also calculated. a determinant of zero indicates that C the matrix is singular. C C C Note C ===== C C If a double precision version of this routine is desired, the C c in column 1 should be removed from the double precision C statement which follows. C C double precision a,d,biga,hold C C the c must also be removed from double precision statements C appearing in other routines used in conjunction with this C routine. C C The double precision version of this subroutine must also C contain double precision fortran functions. abs in statement C 10 must be changed to dabs. C ccc REAL*8 D SUBROUTINE MINV3(A,B,D)C Invert a general 3x3 matrix and return determinant in D C C A = (B)1 C C .. Scalar Arguments .. REAL D C .. C .. Array Arguments .. REAL A(3,3),B(3,3) SUBROUTINE RANMAR(RVEC,LEN)C Universal random number generator proposed by Marsaglia and Zaman C in report FSUSCRI8750 C slightly modified by F. James, 1988 to generate a vector C of pseudorandom numbers RVEC of length LEN C and making the COMMON block include everything needed to C specify completely the state of the generator. C Transcribed from CERN report DD/88/22. C Rather inelegant messing about added by D. Love, Jan. 1989 to C make sure initialisation always occurs. C *** James says that this is the preferred generator. C Gives bitidentical results on all machines with at least C 24bit mantissas in the flotaing point representation (i.e. C all common 32bit computers. Fairly fast, satisfies very C stringent tests, has very long period and makes it very C simple to generate independly disjoint sequences. C See also RANECU. C The state of the generator may be saved/restored using the C whole contents of /RASET1/. C Call RANMAR to get a vector, RMARIN to initialise. C C Argument list C  C C VREC (O) (REAL) Random Vector C C LEN (I) (INTEGER) Length of random vector C C C For ENTRY point RMARIN C  C C Initialisation for RANMAR. The input values should C be in the ranges: 0<=ij<=31328, 0<=kl<=30081 C This shows the correspondence between the simplified input seeds C IJ, KL and the original MarsagliaZaman seeds i,j,k,l C To get standard values in MarsagliaZaman paper, C (I=12, J=34, K=56, L=78) put IJ=1802, KL=9373 C C IJ (I) (INTEGER) Seed for random number generator C C KL (I) (INTEGER) Seed for randon number generator SUBROUTINE SCALEV(A,X,B)C Scale vector B with scalar X and put result in A C C .. Scalar Arguments .. REAL X C .. C .. Array Arguments .. REAL A(3),B(3) SUBROUTINE TRANSP(A,B)C Transpose a 3x3 matrix C C A = BT C C .. Array Arguments .. REAL A(3,3),B(3,3) SUBROUTINE UNIT(V)C Vector V reduced to unit vector C C .. Array Arguments .. REAL V(3) SUBROUTINE VDIF(A,B,C)C Subtract two vectors C C A = B  C C C .. Array Arguments .. REAL A(3),B(3),C(3) SUBROUTINE VSET(A,B)C Copy a vector from B to A C C .. Array Arguments .. REAL A(3),B(3) SUBROUTINE VSUM(A,B,C)C Add two vectors C C A = B + C C C .. Array Arguments .. REAL A(3),B(3),C(3) SUBROUTINE ZIPIN(ID,N,BUF)C Fast binary read on unit ID into real array BUF of length N C C .. Scalar Arguments .. INTEGER ID,N C .. C .. Array Arguments .. REAL BUF(N) SUBROUTINE ZIPOUT(ID,N,BUF)C Fast binary write to unit ID of real array BUF length N C C .. Scalar Arguments .. INTEGER ID,N C .. C .. Array Arguments .. REAL BUF(N) SUBROUTINE GMPRD(A,B,R,N,M,L)C Ssp general matrix product C C R(N,L) = A(N,M) * B(M,L) C C .. Scalar Arguments .. INTEGER L,M,N C .. C .. Array Arguments .. REAL A(N*M),B(M*L),R(N*L) SUBROUTINE MATMUL4(A,B,C)C Multiply two 4x4 matrices, A=B*C C C .. Array Arguments .. REAL A(4,4),B(4,4),C(4,4) subroutine matmln(n,a,b,c)C Multiply two nxn matrices C a = b . c C integer n real a(n,n),b(n,n),c(n,n) SUBROUTINE DETMAT(RMAT,DET)C Calculate determinant DET of 3x3 matrix RMAT C C C .. Scalar Arguments .. REAL DET C .. C .. Array Arguments .. REAL RMAT(3,3) SUBROUTINE MATVC4(A,R,B)C Matrix x Vector C A(3) = R(4x4) . B(3) C REAL A(3), R(4,4), B(3) SUBROUTINE ML3MAT(P,A,Q,B,R,C,S,D)C 26Nov1988 J. W. Pflugrath Cold Spring Harbor Laboratory C Edited to conform to Fortran 77. Renamed from Multiply_3_matrices to C ML3MAT C C ============================================================================== C C! to multiply three matrices C SUBROUTINE ML3MAT C ! input: 1st side of 1st matrix 1 (P C ! input: first matrix 2 ,A C ! input: 2nd side of 1st matrix & 1st side of 2nd matrix 3 ,Q C ! input: second matrix 4 ,B C ! input: 2nd side of 2nd matrix & 1st side of 3rd matrix 5 ,R C ! input: third matrix 6 ,C C ! input: 2nd side of 3rd matrix 7 ,S C ! output: product matrix 8 ,D) C CEE Multiplies three real matrices of any dimensions. It is not optimised C for very large matrices. C Multiply_3_matrices C*** this routine is inefficient! C Multiply_3_matrices Created: 15NOV1985 D.J.Thomas, MRC Laboratory of Molecular Biology, C Hills Road, Cambridge, CB2 2QH, England SUBROUTINE INV44(A,AI)C SUBROUTINE TO INVERT 4*4 MATRICES FOR CONVERSION BETWEEN C FRACTIONAL AND ORTHOGONAL AXES C PARAMETERS C C A (I) 4*4 MATRIX TO BE INVERTED C AI (O) INVERSE MATRIX SUBROUTINE MATMLI(A,B,C)C Integer matrix multiply INTEGER A(3,3),B(3,3),C(3,3) SUBROUTINE MATMULTrans(A,B,C)C A=B*C(transpose) for 3x3 matrices C C ..Array arguments REAL A(3,3),B(3,3),C(3,3) SUBROUTINE IMATVEC(V,A,B)C Postmultiply a 3x3 matrix by a vector C C V=AB C C .. Array Arguments .. INTEGER A(3,3),B(3),V(3) SUBROUTINE TRANSFRM(X,MAT)C Transform vector X(3) by quine matrix MAT(4,4) C Return transformed vector in X. C C ..Array arguments.. REAL X(3),MAT(4,4) SUBROUTINE EIGEN_RS_ASC(A, R, N, MV)C SUBROUTINE TO COMPUTE EIGENVALUES & EIGENVECTORS OF A REAL C SYMMETRIC MATRIX, FROM IBM SSP MANUAL (SEE P165). C DESCRIPTION OF PARAMETERS  C A  ORIGINAL MATRIX STORED COLUMNWISE AS UPPER TRIANGLE ONLY, C I.E. "STORAGE MODE" = 1. EIGENVALUES ARE WRITTEN INTO DIAGONAL C ELEMENTS OF A I.E. A(1) A(3) A(6) FOR A 3*3 MATRIX. C R  RESULTANT MATRIX OF EIGENVECTORS STORED COLUMNWISE IN SAME C ORDER AS EIGENVALUES. C N  ORDER OF MATRICES A & R. C MV = 0 TO COMPUTE EIGENVALUES & EIGENVECTORS. C REAL A(*), R(*) INTEGER N, MV 