procedure gresolpiv (A,B,X,N,IMPOSSIBLE);
value N;
array A,B,X;
integer N;
label IMPOSSIBLE;
Solution of AX = B by Gaussian elimination with row pivoting. N is the order of A. Exit to IMPOSSIBLE if A is singular. A and B are preserved. Declare A[1 :N, 1:N], B,X[1:N].
procedure tridiagreson (N,A1,B1,C1,F,IMPOS);
value N;
integer N;
real array A1,B1,C1,F;
label IMPOS;
Solves AX = B where A is tridiagonal, with diagonal elements in A1[1:N] and sub- and super-diagonal elements respectively in B1,C1[2:N]. Solution in F[1:N]. Exit to IMPOS if any minor is less than 10-6 in modulus. A1,B1, and C1 are destroyed.
procedure deter (A,N);
value N;
integer N;
real array A;
Determinant of A[1:N,1:N] by Gaussian elimination with row pivoting. A is preserved.
procedure trace (A,T,N,P);
value N,P;
integer N,P;
real array A,T;
Traces of powers of A [1:N,1:N] up to AP are produced in T[1:P]. A is preserved.
procedure pivotmax (A,B,N,IMPOSSIBLE);
value N;
array A,B;
integer N;
label IMPOSSIBLE;
The inverse of A[1:N,1:N], computed by Gaussian elimination with row interchanges, is stored in B[1:N,1:N]. Exit to IMPOSSIBLE if A is singular. A is preserved.
Cholesky decomposition of a positive definite symmetric matrix.
procedure detsym 1(A,N);
value N;
array A;
integer N;
Determinant of the +ve def. symmetric matrix M whose upper triangle is input in successive columns to A[1:N×(N+1)/2]. The Cholesky matrix of M is overwritten on A.
procedure solsym 1(A,N,B);
value N;
array A,B;
integer N;
Called after detsym 1 and solves Mx = B. The solution x is overwritten on B[1:N]. The Cholesky matrix of M is preserved in A so the procedure can be called repeatedly for various B.
procedure detinvsym 1(A,N);
value N;
array A;
integer N;
Determinant of M. Also the upper triangle of M-1 is overwritten on A (if det(M) ≥ 0). The procedure calls detsym 1. Note that in detsym 1 the decomposition is discontinued if M is not positive definite. In that case detsym 1:= -k, where k is the stage number at which a fault occurs.
procedure detbnd(A,N,LW,RW,AUX,M,P);
value N,LW,RW;
integer N,LW,RW;
integer array P;
array A,M,AUX;
Determinant of a band matrix order N with LW non-zero left codiagonal elements and RW right codiagonal elements. The elements of the matrix are input into A [1 : (N-1 )×(LW + RW) + N] successively by rows. A relative tolerance must be specified as AUX[0]. Method is standard LU decomposition. U is delivered in A and L in M[1:(N-2)×LW+1]. Pivotal indices are declared in P [1:N].
procedure solbnd(A,N,LW,RW,M,P,B);
value N,LW,RW;
integer N,LW,RW;
integer array P;
array A,B,M;
Solves the system Rx = B where R is the original band matrix described above. The solution vector x is overwritten on B. The procedure should be called after detbnd (if the determinant is non-zero). Since solbnd leaves A and M unaltered, it can be called successively for different right-hand sides B.
procedure lsqdec(A,N,M,AUX,AID,CI);
value N,M;
integer N,M;
array A,AUX,AID;
integer array CI;
Rank of the N × M matrix A [1:N,1:M]. AUX [0] is set in AUX [0:2].
procedure lsqsol(A,N,M,AID,CI,B);
value N,M;
integer N,M;
array A,AID,B;
integer array CI;
Least squares solution x of Ax = B is written on the first M elements of B. The Householder-triangularised form of A is left in A, except for the diagonal elements which are output in AID [1:M]. The pivotal indices are left in CI[1:M]. The procedure should only be called if Rank = M. A, AID and CI are left intact so the procedure can be repeatedly called for different RHS. B.
procedure lsqdglinv(A,M,AID,CI,DIAG);
value M;
integer M;
array A,AID,DIAG;
integer array CI;
Can be called after lsqdec if Rank = M, and calculates the diagonal of the inverse of A'A in DIAG [1:M].
procedure lsqdecsol(A,N,M,AUX,DIAG,B);
value N,M;
integer N,M;
array A,AUX,DIAG,B;
Combines the 3 preceding procedures to compute the least squares solution of Ax = B. The elements of A are altered.
procedure pseudoinv(B,M,N);
array B;
integer M,N;
Pseudoinverse (transposed) of a rectangular or non-singular square matrix B, based on the relation BI = (B'B)IB'. B is M × N where M ≤ N. Note that if rank(B) = N then BI = (B'B)-1B'. Declare B[1:M,1:N]. B is destroyed.
procedure pseudoinverse(A,M,N,IMPOSSIBLE);
array A;
integer M,N;
label IMPOSSIBLE;
Pseudoinverse (transposed) of a rectangular or (possibly singular) M × N matrix A. The procedure borders A to obtain a non-singular matrix hence we most declare A [1:M+N, 1:M+N]. An exit to IMPOSSIBLE is only possible if there b a machine failure, or the M+N - order matrix happens to be nearly singular. A is destroyed.