/************************************************************************/
/* */
/* Algorithm 298. SQUARE ROOT OF A POSITIVE DEFINITE MATRIX. */
/* */
/* H. Spath, Institut fur Neutronenphysik und Reaktortechnik */
/* Kernforschungszentrum Karlsruhe, Germany */
/* */
/************************************************************************/
/* Translated to PL/I by R. A. Vowels, 10 February 2006. */
/* Let A be a symmetric positive-definite matrix of the order N.
Further let Lambda(min) be the smallest end Lambda(max) be the
greatest eigenvalue of A.
It is known [I. Babusea, M. Prager, E. Vitasek, Numerical Processes
in Differential Equations, John Wiley & Sons, London, 1966,
p. 31 ff.] that for all theta with 0 < theta < 1 the sequence
B(i+1) = B(k) + c(A-B(k**2), B(0) = 2cA (1)
with c = theta / (2sqrt(Lambda(max)))
converges to sqrt(A). The rate of convergence of the above
sequence is given by the rate of convergence to zero of the
sequence
x(k) = (1 - theta.sqrt(Lambda(min)/Lambda(max))) ** k (2)
As || A || = alpha Lambda(max) with alpha >= 1, we set
c(1) = theta(1) / (2sqrt( ||A|| ) .
(In the program we choose || A || = max((i) || SIGMA(k) |a(i,k) | ).
Then the sequence (1) with c = c(1) converges for all theta(1) with
0 < theta(1) > sqrt(alpha) = sqrt( ||A|| / Lambda(max))
and therefore in any case for theta(1) with 0 < theta(1) < 1. Because
of (2) it is favorable to choose theta close to 1 and theta(1) close
to sqrt(alpha), respectively. If nothing at all is known about alpha,
the optimum is to choose theta(1) close to 1. The computing time is
proportional to f(theta(1))N**2, where f(theta(1)) decreases as
theta(1) increases.
Meaning of symbols in the formal parameter list:
A = A(1:N, 1:N) must be symnmnetric and positive-definite.
where N is the order of A and B.
A is not destroyed after leaving Matrix_Square_Root.
B = B(1:N, 1:N) contains sqrt(A) when Matrix_Square_Root is left.
theta = theta(1) is an input parameter as described above;
eps is an accuracy parameter. The iteration is stopped when
k+1 k
max | b = b | < eps
i,j i,j i,j
*/
/* The procedure was executed with the 3 x 3 test data given with the */
/* algorithm. Single precision arithmetic was used, and the results */
/* agreed to 5 significant figures. */
Matrix_Square_Root:
procedure (A, B, theta, eps);
declare (A(*,*), B(*,*), theta, eps) float (18);
declare (i, j, k, nits, N) fixed binary, (delta, s, c, bb(hbound(a,1))) float (18);
N = hbound(a,1);
if N ^= hbound(b,2) | N ^= hbound(b,1) | N ^= hbound(b,2) then
do;
put skip list ('The bounds of A and/or B are not identical');
signal error;
end;
/* Determination of c. */
c = 0;
do i = 1 to N;
s = 0;
do j = 1 to N;
s = s + abs(A(i,j));
end;
if c < s then c = s;
end;
c = 0.5 * theta / sqrt(c);
/* now B(0) is set */
do i = 1 to N;
do j = i to N;
B(i,j) = 2 * c * A(i,j); B(j,i) = B(i,j);
end;
end;
/* start of iteration. */
nits = 0;
REPEAT: delta = 0;
nits = nits + 1;
if nits > 100 then return;
put skip list ('Iteration=', nits);
do i = 1 to N;
do j = i to N;
s = 0;
do k = 1 to N;
s = s - B(i,k) * B(k,j);
end;
bb(j) = B(i,j) + c * (A(i,j) + s);
end;
do j = i to N;
s = abs(B(i,j) - bb(j));
if s > delta then delta = s;
B(i,j) = bb(j);
end;
end;
do i = 1 to N-1;
do j = i+1 to N;
B(j,i) = B(i,j);
end;
end;
do i = 1 to N;
put skip;
do j = 1 to N;
put list (B(i,j));
end;
end;
put skip list ('delta=', delta); put skip list (' eps=', eps);
if delta > eps then go to REPEAT;
end Matrix_Square_Root;