%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. */