MATRIX revised by R N Caffin from the original by H W Holdaway CSIRO Division of Textile Physics 338 Blaxland Road Ryde NSW 2112 Australia Page 2 Index| Introduction 3 Installation and Test Programs 4 Routines 5 References 6 Page 3 Introduction This matrix package comes in two versions: MATRIX, for REAL*4 variables, and CMTRIX, for COMPLEX*8 variables. It performs the usual simple operations such as addition, subtraction and multiplication (with and without in-line transposes), and two more complex operations: the Moore-Penrose Generalised Inverse (or psuedo inverse) and the Singular Value Decomposition. For a theoretical treatment of the subject (this sort of stuff is heavy) the reader should consult books on matrices, linear equations and numerical analysis. A couple are listed at the end. The ordinary inverse of a square matrix is well known, provided that the matrix is not singular. If the matrix is singular the determinant is zero and there is no solution, and if the matrix is not square there is normally no determinant. However, if the problem is treated as involving a set of linear equations with noise, a least squares solution to the intersection of the equations can at least be found. That is, a solution to A.x=c where A and c are given, (even though there are more rows of data in A than variables in x,) can be found. This is done by finding the Generalised Inverse of A, denoted hereafter by A'. The singular value decomposition of A gives the arrays U, Q and V where A=U'.Q.V, where U' is the transpose of U, U.U'=V.V'=I and Q is the diagonal ordered array of singular values of A'.A . The number of non-zero elements in Q gives the rank of A, allowing for round-off errors and so on. In a suitable case the singular values are the same as the eigenvalues of the equation Ax=0. It should be noted that the package uses several (internal) routines for the calculation of dot products. In the REAL*4 case the calculations in these dot product routines are done in ordinary double precision. In the COMPLEX*8 case the calculations are done again in double precision by breaking out the real and imaginary parts in this version, as COMPLEX*16 is not supported in the all versions of Fortran. At present the complex SVD routine does an immediate return with an error message, as I am not too sure of the theoretical meaning of some of the steps in the reductions in the complex plane. In particular, a>b has no clear mean- ing and is not supported in the complex case. Advice from mathematicians would be appreciated. It is interesting that this field is relatively modern: some of my text books dating back to the mid-60's do not cover the GINV and SVD in any real detail at all. Others from that period reflect only too painfully techniques relevant to a manual (ie mechanical, like Facit) calculator! It should be noted that while the original author (HWH) well understood what he was doing, the present author is NOT a professional mathematician. The routines appear to work reliably, but no guarantees are offered. Page 4 Installation Installation requires no more than the compiling of the MATRIX or CMTRIX file. For greater efficiency it is very useful to convert the resulting object file into a library: this saves both link time and memory space. Command files MATRIX.COM and CMTRIX.COM are on the distribution for these operations. Test Programs Also included on the distribution are two test programs MATEST and CMTEST and their associated COM files for compilation and execution. MATEST is for the real case and CMTEST is for the complex case. These programs require the appropiately named DAT files to be on DK: for execution. Running the test programs requires only that you enter the size of the arrays to be used: they are m rows by n columns. You can pick square arrays by making m=n, or non-square by making m>n. The case of m A, A(m,n)=0, m>=n MXUNT(A,m,n) I -> A, A(m,n)=I, m=n MXMOV(A,m,n,R) A -> R, R(m,n)=A(m,n) MXADD(A,B,m,n,R) A+B -> R, A(m,n), B(m,n), R(m,n) MXSUB(A,B,m,n,R) A-B -> R, A(m,n), B(m,n), R(m,n) MXSCM(A,m,n,s) s*A -> A, A(m,n), scalar multiply MXTRIM(A,m,n,e) A -> A, zero all elements of A(m,n) < e MXWRT(A,m,n) TYPE out A(m,n) MXPRT(A,m,n) PRINT out A(m,n) MXTRAP(A,m,n) A' -> A, transpose of A(m,n) MXNORM(A,m,n,s) |A| -> s, A(m,n), Frobenius norm s MXMLT(A,B,m,l,n,R) A.B -> R, A(m,l), B(l,n), R(m,n) MXMAT(A,B,m,l,n,R) A'.B-> R, A(l,m), B(l,n), R(M,n) MXMBT(A,B,m,l,n,R) A.B'-> R, A(m,l), B(n,l), R(m,n) GINV(A,m,n,U,AT,NF,C,X,ix,CQ) Solves A.x=C to give A' A(m,n), U(n,m), AT(n-1), NF(n-1), CQ(n), C(m), X(n) for ix=0 returns A', the generalised inverse for ix=2 returns x, the eigenvector solution for ix=3 returns A' and x. SVD(A,m,n,U,V,Q) A -> U.Q.V', where U.U'=V.V'=I and Q(1,n) contains the singular values of A'.A st Q(i)>Q(i+1)(real only) MLT(A,B,m,l,n,R) Internal routine for A.B DOT(A,m,n,j,k) Internal routine DOT1(A,m,n,j,k) Internal routine DOT2(A,m,n,j,k) Internal routine DOT3(A,m,n,j,k) Internal routine Notes: 1. Variables l,m,n etc are Integer*2 (not *4). 2. Variables e & s are either real*4 or complex*8. Use of a real variable where a complex variable is required can/will have the predicable strange consequences. 3. MXTRIM "trims" elements of A to zero when they are less than some user-specified error limit e. In the complex case the real and imaginary components are dealt with separately. Page 6 References Stewart G W, Introduction to Matrix Computations, AP, 1973