The 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.
The subroutines have been arranged as much as possible into groups with similar functions or origins.
The routines are divided into six sections below: vector operations, matrix-vector operations, matrix operations, Harwell routines, Bessel function routines, and miscellaneous (i.e. those which don't fit in anywhere else).
| Type | Name | Function | Call using | Arguments |
|---|---|---|---|---|
| s/r | cross | compute vector product A = B x C | CROSS(A,B,C) | REAL: A(3), B(3), C(3) |
| fun | dot | dot product of two vectors A.B | DOT(A,B) | REAL: A(3), B(3) |
| s/r | icross | cross product (integer version) A = B x C | ICROSS(A,B,C) | INTEGER: A(3), B(3), C(3) |
| fun | idot | dot product (integer version) A.B | IDOT(A,B) | INTEGER: A(3), B(3) |
| s/r | unit | Vector V reduced to unit vector | UNIT(V) | REAL V(3) |
| s/r | scalev | Scale vector B with scalar X and put result in A | SCALEV(A,X,B) | REAL: A(3), B(3), X |
| s/r | vsum | Add two vectors A = B + C | VSUM(A,B,C) | REAL: A(3),B(3),C(3) |
| s/r | vdif | Subtract two vectors A = B - C | VDIF(A,B,C) | REAL: A(3),B(3),C(3) |
| s/r | vset | Copy a vector from B to A | VSET(A,B) | REAL A(3),B(3) |
| Type | Name | Function | Call using | Arguments |
|---|---|---|---|---|
| s/r | matvec | Post-multiply a 3x3 matrix by a vector V = AB | MATVEC(V,A,B) | REAL: A(3,3), B(3), V(3) |
| s/r | imatvec | Post-multiply a 3x3 matrix by a vector V=AB (integer version) | IMATVEC(V,A,B) | INTEGER A(3,3),B(3),V(3) |
| s/r | matvc4 | Matrix x Vector, A(3) = R(4x4) . B(3) (see notes) | MATVC4(A,R,B) | REAL A(3), R(4,4), B(3) |
| s/r | transfrm | Transform vector X(3) by quine matrix MAT(4,4), return transformed vector in X (see notes) | TRANSFRM(X,MAT) | REAL X(3),MAT(4,4) |
Notes on Matrix-Vector 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 3-vector by the 4x4 matrix is effected by treating the input vector as a 4-vector, 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.
| Type | Name | Function | Call using | Arguments |
|---|---|---|---|---|
| s/r | minv3 | Invert a general 3x3 matrix and return determinant in
D A=(B)-1 | MINV3(A,B,D) | REAL: A(3,3), B(3,3), D |
| s/r | iminv3 | Invert a general 3x3 matrix and return determinant in D (integer version) A = (B)-1 | IMINV3(A,B,D) | INTEGER: A(3,3),B(3,3), D |
| s/r | minvn | Invert a matrix A of order N (overwritten on output) and return determinant D | MINVN(A,N,D,L,M) | INTEGER: N REAL: A(N,N), D, L(N), M(N) L,M are work-vectors |
| s/r | matmul | Multiply two 3x3 matrices A = BC | MATMUL(A,B,C) | REAL: A(3,3),B(3,3),C(3,3) |
| s/r | matmulnm | Multiply NxM MXN matrices A = BC | MATMULNM(N,M,A,B,C) | INTEGER: N, M REAL: A(N,N), B(N,M), C(M,N) |
| s/r | matmulgen | Generalised matrix multiplication subroutine Multiplies a NbxMbc matrix (B) by a MbcXNc (C) matrix, so that A = BC | MATMULGEN(Nb,Mbc,Nc,A,B,C) | INTEGER: Nb, Mbc, Nc REAL: A(Nb,Nc),B(Nb,Mbc),C(Mbc,Nc) |
| s/r | transp | Transpose a 3x3 matrix A = BT | TRANSP(A,B) | REAL: A(3,3),B(3,3) |
| s/r | gmprd | general matrix product R(N,L) = A(N,M) * B(M,L) | GMPRD(A,B,R,N,M,L) | REAL A(N*M),B(M*L),R(N*L) INTEGER L,M,N |
| s/r | matmul4 | Multiply two 4x4 matrices, A=B*C | MATMUL4(A,B,C) | REAL A(4,4),B(4,4),C(4,4) |
| s/r | matmln | Multiply two nxn matrices, a = b . c | MATMLN(N,A,B,C) | INTEGER N REAL A(N,N),B(N,N),C(N,N) |
| s/r | detmat | Calculate determinant DET of 3x3 matrix RMAT | DETMAT(RMAT,DET) | REAL DET,RMAT(3,3) |
| s/r | ml3mat | Multiply together three matrices of any size D = A.B.C | ML3MAT(P,A,Q,B,R,C,S,D) | REAL A(P,Q),B(Q,R),C(R,S),D(P,S) INTEGER P,Q,R,S |
| s/r | inv44 | 4x4 matrix inversion, returns IA as inverse of A | INV44(A,AI) | REAL A(4,4),AI(4,4) |
| s/r | matmli | Multiply two 3x3 matrices A = BC (integer version) | MATMLI(A,B,C) | INTEGER A(3,3),B(3,3),C(3,3) |
| s/r | matmultrans | A=B*CT for 3x3 matrices | MATMULTrans(A,B,C) | REAL A(3,3),B(3,3),C(3,3) |
| s/r | eigen_rs_asc | eigenvalues and eigenvectors of real symmetric matrix | EIGEN_RS_ASC(A, R, N, MV) | REAL A(*), R(*) INTEGER N, MV |
The 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 modifications 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.
| Type | Name | Function |
|---|---|---|
| s/r | ea06c | Given a real MxM symmetric matrix A = {aij} this routine finds all its eigenvalues (lambda)i i=1,2,.....,m and eigenvectors xj i=1,2,...,m. i.e. finds the non-trivial solutions of Ax=(lambda)x |
| s/r | ea08c | This uses QR iteration to find the all the eigenvalues and eigenvectors of a real symmetric tri-diagonal matrix |
| s/r | ea09c | Purpose is unspecified (no comments in header) |
| fun | fa01as | Purpose is unspecified (no comments in header) |
| s/r | fa01bs | Purpose is unspecified (no comments in header) |
| s/r | fa01cs | Purpose is unspecified (no comments in header) |
| s/r | fa01ds | Purpose is unspecified (no comments in header) |
| fun | fm02ad | Compute the inner product of two double precision real vectors accumulating the result double precision, when the elements of each vector are stored at some fixed displacement from neighbouring elements. |
| s/r | mc04b | Transforms a real symmetric matrix A={a(i,j)}, i, j=1..IA into a tri-diagonal matrix having the same eigenvalues as A using Householder's method. |
Subroutine ZJVX is from Ian Tickle. Subroutines JDVX, JVX, GAMMA and functions MSTA1, MSTA2 are from MJYV at http://iris-lee3.ece.uiuc.edu/~jjin/routines/routines.html.
That site states:
All the programs and subroutines contained in this archive are copyrighted. However, we give permission to the user who downloads these routines to incorporate any of these routines into his or her programs provided that the copyright is acknowledged.
| Type | Name | Function |
|---|---|---|
| s/r | ZJVX | Compute zero of Bessel function of 1st kind Jv(x) |
| s/r | JDVX | Compute Bessel functions of 1st kind Jv(x) and their derivatives |
| s/r | JVX | Compute Bessel functions of 1st kind Jv(x) |
| s/r | GAMMA | Compute gamma function GA(x) |
| fun | MSTA1 | Determine the starting point for backward recurrence such that the magnitude of Jn(x) at that point is about 10^(-MP) |
| fun | MSTA2 | Determine the starting point for backward recurrence such that all Jn(x) has MP significant digits |
| Type | Name | Function | Call using | Arguments |
|---|---|---|---|---|
| s/r | ranmar | Universal random number generator; returns a random vector RVEC(LEN) | RANMAR(RVEC,LEN) | INTEGER: LEN REAL: RVEC(LEN) |
| ent | rmarin | initialisation for ranmar | RMARIN(IJ,KL) | INTEGER: IJ, KL (seeds for ranmar) |
| s/r | zipin | Fast binary read on unit ID into real array BUF of length N | ZIPIN(ID,N,BUF) | INTEGER: ID, N REAL: BUF(N) |
| s/r | zipout | Fast binary write to unit ID of real array BUF length N | ZIPOUT(ID,N,BUF) | INTEGER: ID, N REAL: BUF(N) |
From 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.
BLAS contains routines for vector-vector, vector-matrix and matrix-matrix operations. Pre-compiled 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 (warning the list may not be up-to-date), 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://math-atlas.sourceforge.net/ for more information. ATLAS can be also be used on systems which have a vendor BLAS library provided.
LAPACK is built on the BLAS routines and provides routines for solving systems of simultaneous equations, least-squares 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.
Installation 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 --disable-lapack 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:
| System | LAPACK and BLAS | BLAS only | See also: |
|---|---|---|---|
| OSF1 (DEC/Compaq) | libcxml,libdxml | [none known] | http://www5.compaq.com/math/documentation/index.html |
| IRIX (SGI) | libscs | libblas | http://www.sgi.com/software/scsl.html |
| HPUX | HP mlib | libblas | http://h21007.www2.hp.com/dspp/tech/tech_TechSoftwareDetailPage_IDX/1,1703,1204,00.html |
| SUNOS | Sun Performance Library | [none known] | http://docs.sun.com/db/doc/806-7995 |
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:
The Netlib LAPACK routines will normally only be built if configure cannot find a system LAPACK library.
As of CCP4 v5.0, all the versions of the routines are built (i.e. single (real), double precision, complex and complex*16). versions of the routines.
If you are installing on a non-IEEE compliant machine then there is a small but very important change which must be made to prevent run-time crashes; see below. Configure should normally make the change automatically, but if you encounter problems then please let us know.
Note that building of the Netlib routines can also be forced by specifying the --with-netlib-lapack 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 e-mail then to ccp4@ccp4.ac.uk
By default the NetLib version of the LAPACK routine ILAENV (distributed with CCP4) assumes an IEEE-compliant machine, which implies certain behaviour for NaN and "infinity arithmetic". On non-IEEE compliant machines or O/S the LAPACK library will compile successfully but may crash at run-time when ILAENV is called.
The CCP4 distribution contains two versions of the ILAENV routine, in $CCP4/lib/lapack/src/. ilaenv-dist.f is the standard Netlib version of the routine, whilst ilaenv-non-ieee.f is the version modified for a non-IEEE compliant machine.
The CCP4 configure will copy the appropriate version of the routine to ilaenv.f on installation.
Compute vector product A = B x C
.. Array Arguments ..
REAL A(3),B(3),C(3)
dot product of two vectors
.. Array Arguments ..
REAL A(3),B(3)
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 non-trivial
C solutions of Ax=(lambda)x
C
C The matrix is reduced to tri-diagonal 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 2-dimensional array VECTOR
C
C VECTOR(IV,M) (O) (real) 2-dimensional 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.
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 tri-diagonal matrix C whose diagonal elements are A(i), i=1,M and off-diagonal 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) Off-diagonal elements C C IV (I) (integer) should be set to the first dimen- C sion of the 2-dimensional 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.
C 18/03/70 LAST LIBRARY UPDATE[No info]
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 multi-dimensional 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
C Cross product (integer version)
C
C A = B x C
C
C .. Array Arguments ..
INTEGER A(3),B(3),C(3)
C Dot product (integer version)
C
C IDOT = A . B
C
C .. Array Arguments ..
INTEGER A(3),B(3)
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)
C Multiply two 3x3 matrices
C
C A = BC
C
C .. Array Arguments ..
REAL A(3,3),B(3,3),C(3,3)
C Multiply NxM MXN matrices
C
C A = BC
C
C .. Array Arguments ..
REAL A(N,N),B(N,M),C(M,N)
C Post-multiply a 3x3 matrix by a vector
C
C V = AB
C
C .. Array Arguments ..
REAL A(3,3),B(3),V(3)
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)
C Transforms a real symmetric matrix A={a(i,j)}, i, j=1..IA
C into a tri-diagonal 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)
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 gauss-jordan 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
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)
C Universal random number generator proposed by Marsaglia and Zaman C in report FSU-SCRI-87-50 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 bit-identical results on all machines with at least C 24-bit mantissas in the floating point representation (i.e. C all common 32-bit computers. Fairly fast, satisfies very C stringent tests, has very long period and makes it very C simple to generate independently 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 Marsaglia-Zaman seeds i,j,k,l C To get standard values in Marsaglia-Zaman 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 random number generator
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)
C---- Transpose a 3x3 matrix
C
C A = BT
C
C .. Array Arguments ..
REAL A(3,3),B(3,3)
C Vector V reduced to unit vector
C
C .. Array Arguments ..
REAL V(3)
C Subtract two vectors
C
C A = B - C
C
C .. Array Arguments ..
REAL A(3),B(3),C(3)
C Copy a vector from B to A
C
C .. Array Arguments ..
REAL A(3),B(3)
C Add two vectors
C
C A = B + C
C
C .. Array Arguments ..
REAL A(3),B(3),C(3)
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)
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)
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)
C Multiply two 4x4 matrices, A=B*C
C
C .. Array Arguments ..
REAL A(4,4),B(4,4),C(4,4)
C Multiply two nxn matrices
C a = b . c
C
integer n
real a(n,n),b(n,n),c(n,n)
C---- Calculate determinant DET of 3x3 matrix RMAT
C
C
C .. Scalar Arguments ..
REAL DET
C ..
C .. Array Arguments ..
REAL RMAT(3,3)
C Matrix x Vector
C A(3) = R(4x4) . B(3)
C
REAL A(3), R(4,4), B(3)
C 26-Nov-1988 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: 15-NOV-1985 D.J.Thomas, MRC Laboratory of Molecular Biology,
C Hills Road, Cambridge, CB2 2QH, England
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
C Integer matrix multiply
INTEGER A(3,3),B(3,3),C(3,3)
C A=B*C(transpose) for 3x3 matrices
C
C ..Array arguments
REAL A(3,3),B(3,3),C(3,3)
C---- Post-multiply a 3x3 matrix by a vector
C
C V=AB
C
C .. Array Arguments ..
INTEGER A(3,3),B(3),V(3)
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)
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
C Purpose: Compute zero of Bessel function Jv(x) C Input : V --- Order of Jv(x) C IZ --- Index of previous zero of Jv(x) C X --- Value of previous zero of Jv(x) C Output: BJ(N) --- Jv(x) for N = 0...INT(V) C DJ(N) --- J'v(x) for N = 0...INT(V) C X --- Value of next zero of Jv(x)
C Purpose: Compute Bessel functions Jv(x) and their derivatives C Input : x --- Argument of Jv(x) C v --- Order of Jv(x) C ( v = n+v0, 0 <= v0 < 1, n = 0,1,2,... ) C Output: BJ(n) --- Jn+v0(x) C DJ(n) --- Jn+v0'(x) C VM --- Highest order computed C Routines called: C (1) GAMMA for computing gamma function C (2) MSTA1 and MSTA2 for computing the starting C point for backward recurrence
C Purpose: Compute Bessel functions Jv(x) C Input : x --- Argument of Jv(x) C v --- Order of Jv(x) C ( v = n+v0, 0 <= v0 < 1, n = 0,1,2,... ) C Output: BJ(n) --- Jn+v0(x) C VM --- Highest order computed C Routines called: C (1) GAMMA for computing gamma function C (2) MSTA1 and MSTA2 for computing the starting C point for backward recurrence
C Purpose: Compute gamma function GA(x) C Input : x --- Argument of GA(x) C ( x is not equal to 0,-1,-2,... ) C Output: GA --- GA(x)
C Purpose: Determine the starting point for backward C recurrence such that the magnitude of C Jn(x) at that point is about 10^(-MP) C Input : x --- Argument of Jn(x) C MP --- Value of magnitude C Output: MSTA1 --- Starting point
C Purpose: Determine the starting point for backward C recurrence such that all Jn(x) has MP C significant digits C Input : x --- Argument of Jn(x) C n --- Order of Jn(x) C MP --- Significant digit C Output: MSTA2 --- Starting point