How to resolve the algorithm Cholesky decomposition step by step in the PL/I programming language
How to resolve the algorithm Cholesky decomposition step by step in the PL/I programming language
Table of Contents
Problem Statement
Every symmetric, positive definite matrix A can be decomposed into a product of a unique lower triangular matrix L and its transpose:
L
{\displaystyle L}
is called the Cholesky factor of
A
{\displaystyle A}
, and can be interpreted as a generalized square root of
A
{\displaystyle A}
, as described in Cholesky decomposition. In a 3x3 example, we have to solve the following system of equations: We can see that for the diagonal elements (
l
k k
{\displaystyle l_{kk}}
) of
L
{\displaystyle L}
there is a calculation pattern: or in general: For the elements below the diagonal (
l
i k
{\displaystyle l_{ik}}
, where
i
k
{\displaystyle i>k}
) there is also a calculation pattern: which can also be expressed in a general formula: Task description The task is to implement a routine which will return a lower Cholesky factor
L
{\displaystyle L}
for every given symmetric, positive definite nxn matrix
A
{\displaystyle A}
. You should then test it on the following two examples and include your output. Example 1: Example 2:
Let's start with the solution:
Step by Step solution about How to resolve the algorithm Cholesky decomposition step by step in the PL/I programming language
Source code in the pl/i programming language
(subscriptrange):
decompose: procedure options (main); /* 31 October 2013 */
declare a(*,*) float controlled;
allocate a(3,3) initial (25, 15, -5,
15, 18, 0,
-5, 0, 11);
put skip list ('Original matrix:');
put edit (a) (skip, 3 f(4));
call cholesky(a);
put skip list ('Decomposed matrix');
put edit (a) (skip, 3 f(4));
free a;
allocate a(4,4) initial (18, 22, 54, 42,
22, 70, 86, 62,
54, 86, 174, 134,
42, 62, 134, 106);
put skip list ('Original matrix:');
put edit (a) (skip, (hbound(a,1)) f(12) );
call cholesky(a);
put skip list ('Decomposed matrix');
put edit (a) (skip, (hbound(a,1)) f(12,5) );
cholesky: procedure(a);
declare a(*,*) float;
declare L(hbound(a,1), hbound(a,2)) float;
declare s float;
declare (i, j, k) fixed binary;
L = 0;
do i = lbound(a,1) to hbound(a,1);
do j = lbound(a,2) to i;
s = 0;
do k = lbound(a,2) to j-1;
s = s + L(i,k) * L(j,k);
end;
if i = j then
L(i,j) = sqrt(a(i,i) - s);
else
L(i,j) = (a(i,j) - s) / L(j,j);
end;
end;
a = L;
end cholesky;
end decompose;
You may also check:How to resolve the algorithm Roots of a quadratic function step by step in the Nim programming language
You may also check:How to resolve the algorithm Apply a callback to an array step by step in the FunL programming language
You may also check:How to resolve the algorithm Binary digits step by step in the FALSE programming language
You may also check:How to resolve the algorithm World Cup group stage step by step in the Raku programming language
You may also check:How to resolve the algorithm Convert seconds to compound duration step by step in the Erlang programming language