%PROCESS MARGINS(1,100);
/* WR-21 */
/* SYMMETRIC DECOMPOSITION OF A POSITIVE DEFINITE MATRIX - LDLt. */
/* From J. H. Wilkinson & C. Reinsch, "Handbook for Automatic Computation",
Vol. II, "Linear Algebra", Springer-Verlag, New York, 1971, pp. 21-24. */
/* Input:
n = order of the matrix a.
a = elements of the positive definite matrix A, in upper triangular form.
Only the diagonal and above-diagonal elements are required.
Output:
a = elements of the lower triangle L of the LDLt decomposition of A
(excepting the diagonal elements). These are stored below the diagonal.
p = a vector consisting of the reciprocals of the diagonal elements of D.
d1 and d2 = two numbers which together give the determinant of A. */
/* The upper triangle of a positive definite symmetric matrix, A, is
stored in the upper triangle of an n x n array a[i,j], i = 1(1)n, j = 1(1)n.
The decomposition A = LDU, where U is the tranpose of L, which is
a unit lower triangular matrix, and D is a diagonal matrix, is
performed and L is stored in the remainder of the array a, omitting the
unit diagonal, with the reciprocals of the elements of D stored in the
array p[i], i = 1(1)n. A is retained so that the solution obtained can
be subsequently improved. The determinant, d1 x 2^d2, of A is also
computed. The procedure will fail if A, modified by the rounding
errors, is singular. */
/* P. 22 */
/* Input:
n = order of the matrix a.
r = number of right-hand sides for which Ax = b is to be solved.
a = elements of L the lower triangle of the LDLt decomposition of a
positive definite matrix A as produced by procedure symdet.
p = a vector consisting of the reciprocals of the diagonal elements of L
as produced by procedure symdet.
b = the n x r matrix consisting of the r right-hand sides.
Output:
x = the n x r matrix consisting of the r solution vectors. */
/* Solves Ax = b, where A is a positive definite symmetric matrix and
b is an n x r matrix of r right-hand sides. The procedure symbol must
be preceded by symdet in which L is produced in a[i,j] and D is
produced in p[i], from A. Ax = b is solved in two steps, Ly = b and
DUx = y. The matrix b is retained in order to facilitate the refinement
of x, but x is overwritten on y. However, x and b can be identified in
the call of the procedure. */
/* P. 22. */
/* Input:
n = order of the matrix a.
a = elements of the positive definite matrix A, in upper triangular form.
Only the diagonal and above-diagonal elements are required.
Output:
a = elements of inv(A). */
/* The upper triangle of a positive definite symmetric matrix, A, is
stored in the upper triangle of an (n+1) x n array a[i,j], i = 1(1)n + 1,
j = 1(1)n. The decomposition A = LDU, where U is the tranpose of L,
which is a unit lower triangular matrix, and D is a diagonal matrix,
is performed and L and D together are stored in the remainder of the
array a, omitting the unit diagonal. The reciprocals of the elements
of D stored instead of the elements themselves. L is then replaced
by its inverse and then both D and the inverse of L are then replaced by
the lower triangle of the inverse of A. A is retained so that the inverse
can be subsequently be improved. The procedure will fail if A, modified
by the rounding errors, is singular. */