F112 Manipulation of Triangular and Symmetric Matrices

At CERN, matrices are often stored row-wise (TC-convention); furthermore, symmetric matrices are stored packed as the lower left triangular part only, i.e., the Ith diagonal element is found in position $I(I+1)/2$ . The TR-package performs many of the frequently required operations associated with such matrices without resorting to expanding into the unpacked square form. In all the following routines an $M\; xM$ symmetric matrix is taken to be stored in the packed form with $M(M+1)/2$ elements.

Some of these operations produce and require the manipulation of lower triangular matrices which have all elements zero above the leading diagonal. These are also stored in the packed form with all the zeros dropped; therefore, care has to be taken in the interpretation of a packed matrix as to whether it represents a symmetric or lower triangular array. To facilitate this distinction in the Write-up, the following nomenclature has been adopted:

A,B,C

unpacked rectangular matrices (row-wise storage)

Q,R,S,T

packed symmetric matrices

V,W

packed lower triangular matrices

Structure:

SUBROUTINE subprograms
User Entry Names:

TRCHUL,	TRCHLU,	TRSMUL,	TRSMLU,	TRINV,	TRSINV,	TRLA,	TRLTA,
TRAL,	TRALT,	TRSA,	TRAS,	TRSAT,	TRATS,	TRAAT,	TRATA,
TRASAT,	TRATSA,	TRQSQ,	TRPCK,	TRUPCK

Usage:

Choleski Decomposition

CALL TRCHUL(S,W,M)	$S\; =\; W\text{'}W$
CALL TRCHLU(S,V,M)	$S\; =\; VV\text{'}$

S is an $M\; xM$ positive semi-definite symmetric matrix (e.g., error or weight matrix) and the routines calculate the complementary lower triangular Choleski factors. It is allowed to overwrite S by W or V.

Symmetric Multiplication of Lower Triangular Matrices

CALL TRSMUL(W,S,M)	$W\text{'}W$	S
CALL TRSMLU(W,R,M)	$WW\text{'}$	R

W is an $M\; xM$ lower triangular matrix and S, R the two symmetric products of the multiplication of W by its transpose. It is allowed to overwrite W by either S or R.

Lower Triangular Matrix Inversion

CALL TRINV(W,V,M)

$W-1$

V

W is an $M\; xM$ lower triangular matrix which is inverted into V (the inverse of a lower triangular matrix is lower triangular). W may have rows and columns of zeros as produced by the Choleski decomposition of a weight matrix with unmeasured variables. It is allowed to overwrite W by V.

Symmetric Matrix Inversion

CALL TRSINV(S,R,M)

$S-1$

R

S is an $M\; xM$ positive semi-definite symmetric matrix which is inverted into R (also stored packed). It is permissible to overwrite S by R.

Triangular -- Rectangular Multiplication

CALL TRLA (W,A,B,M,N)	$WA$	B
CALL TRLTA(W,A,B,M,N)	$W\text{'}A$	B
CALL TRAL (A,V,B,M,N)	AV	B
CALL TRALT(A,V,B,M,N)	$AV\text{'}$	B

A and B are $M\; xN$ rectangular matrices, W is an $M\; xM$ lower triangular matrix, and V is an $N\; xN$ lower triangular matrix. In each call it is allowed to overwrite A by B.

Symmetric - Rectangular Multiplication

CALL TRSA (S,A,C,M,N)	SA	C
CALL TRAS (A,R,C,M,N)	AR	C
CALL TRSAT(S,B,C,M,N)	$SB\text{'}$	C
CALL TRATS(B,R,C,M,N)	$B\text{'}R$	C

A and C are $M\; xN$ rectangular matrices, B is an $N\; xM$ matrix, S is an $M\; xM$ symmetrix matrix, and R is an $N\; xN$ symmetric matrix. It is not allowed to overwrite A or B by the product matrix C.

Symmetric Multiplication of Rectangular Matrices

CALL TRAAT(A,S,M,N)	$AA\text{'}$	S
CALL TRATA(B,R,M,N)	$B\text{'}B$	R

A is an $M\; xN$ matrix, B is an $N\; xM$ matrix, S is an $M\; xM$

symmetric matrix, and R is an $M\; xM$ symmetric matrix. No overwriting is allowed.

Transformation of Symmetric Matrix

CALL TRASAT(A,S,R,M,N)	$ASA\text{'}$	R
CALL TRATSA(B,S,R,M,N)	$B\text{'}SB$	R
CALL TRQSQ (Q,T,R,M)	QTQ	R

A is an $M\; xN$ matrix, B is an $N\; xM$ matrix, S is an $N\; xN$

symmetric matrix, and R, Q, T are $M\; xM$ symmetric matrices. No overwriting is allowed.

Packing and Unpacking a Symmetric Matrix

CALL TRPCK (A,S,M)	A	S
CALL TRUPCK(S,A,M)	S	A

A is an $M\; xM$ unpacked symmetric matrix (all $M2$ elements) and S is the same matrix stored packed. Overwriting is allowed for both TRPCK and TRUPCK.
$•$

F116