Cholesky decomposition: Difference between revisions
m
→{{header|Wren}}: Minor tidy
Thundergnat (talk | contribs) m (Automated syntax highlighting fixup (second round - minor fixes)) |
m (→{{header|Wren}}: Minor tidy) |
||
| (7 intermediate revisions by 6 users not shown) | |||
Line 308:
( 9.89949, 1.62455, 1.84971, 1.39262))
</pre>
=={{header|Arturo}}==
<syntaxhighlight lang="arturo">cholesky: function [m][
result: array.of: @[size m, size m] 0.0
loop 0..dec size m\0 'i [
loop 0..i 'j [
s: 0.0
loop 0..j 'k ->
s: s + result\[i]\[k] * result\[j]\[k]
result\[i]\[j]: (i = j)? -> sqrt m\[i]\[i] - s
-> (1.0 // result\[j]\[j]) * (m\[i]\[j] - s)
]
]
return result
]
printMatrix: function [a]->
loop a 'b ->
print to [:string] .format:"8.5f" b
m1: @[
@[25.0, 15.0, neg 5.0]
@[15.0, 18.0, 0.0]
@[neg 5.0, 0.0, 11.0]
]
printMatrix cholesky m1
print ""
m2: [
[18.0, 22.0, 54.0, 42.0]
[22.0, 70.0, 86.0, 62.0]
[54.0, 86.0, 174.0, 134.0]
[42.0, 62.0, 134.0, 106.0]
]
printMatrix cholesky m2</syntaxhighlight>
{{out}}
<pre> 5.00000 0.00000 0.00000
3.00000 3.00000 0.00000
-1.00000 1.00000 3.00000
4.24264 0.00000 0.00000 0.00000
5.18545 6.56591 0.00000 0.00000
12.72792 3.04604 1.64974 0.00000
9.89949 1.62455 1.84971 1.39262</pre>
=={{header|ATS}}==
<syntaxhighlight lang="ats">
%{^
#include <math.h>
#include <float.h>
%}
#include "share/atspre_staload.hats"
macdef NAN = g0f2f ($extval (float, "NAN"))
macdef Zero = g0i2f 0
macdef One = g0i2f 1
(* The sqrt(3) function made part of the ‘g0float’ typekind series.
(The ats2-xprelude package will do this for you, but it is easy
to do if you are not using a lot of math functions. *)
extern fn {tk : tkind} g0float_sqrt : g0float tk -<> g0float tk
overload sqrt with g0float_sqrt
implement g0float_sqrt<fltknd> x = $extfcall (float, "sqrtf", x)
implement g0float_sqrt<dblknd> x = $extfcall (double, "sqrt", x)
implement g0float_sqrt<ldblknd> x = $extfcall (ldouble, "sqrtl", x)
(*------------------------------------------------------------------*)
(* A "very little matrix library" *)
typedef Matrix_Index_Map (m1 : int, n1 : int, m0 : int, n0 : int) =
{i1, j1 : pos | i1 <= m1; j1 <= n1}
(int i1, int j1) -<cloref0>
[i0, j0 : pos | i0 <= m0; j0 <= n0]
@(int i0, int j0)
datatype Real_Matrix (tk : tkind,
m1 : int, n1 : int,
m0 : int, n0 : int) =
| Real_Matrix of (matrixref (g0float tk, m0, n0),
int m1, int n1, int m0, int n0,
Matrix_Index_Map (m1, n1, m0, n0))
typedef Real_Matrix (tk : tkind, m1 : int, n1 : int) =
[m0, n0 : pos] Real_Matrix (tk, m1, n1, m0, n0)
typedef Real_Vector (tk : tkind, m1 : int, n1 : int) =
[m1 == 1 || n1 == 1] Real_Matrix (tk, m1, n1)
typedef Real_Row (tk : tkind, n1 : int) = Real_Vector (tk, 1, n1)
typedef Real_Column (tk : tkind, m1 : int) = Real_Vector (tk, m1, 1)
extern fn {tk : tkind}
Real_Matrix_make_elt :
{m0, n0 : pos}
(int m0, int n0, g0float tk) -< !wrt >
Real_Matrix (tk, m0, n0, m0, n0)
extern fn {tk : tkind}
Real_Matrix_copy :
{m1, n1 : pos}
Real_Matrix (tk, m1, n1) -< !refwrt > Real_Matrix (tk, m1, n1)
extern fn {tk : tkind}
Real_Matrix_copy_to :
{m1, n1 : pos}
(Real_Matrix (tk, m1, n1), (* destination *)
Real_Matrix (tk, m1, n1)) -< !refwrt >
void
extern fn {}
Real_Matrix_dimension :
{tk : tkind}
{m1, n1 : pos}
Real_Matrix (tk, m1, n1) -<> @(int m1, int n1)
extern fn {tk : tkind}
Real_Matrix_get_at :
{m1, n1 : pos}
{i1, j1 : pos | i1 <= m1; j1 <= n1}
(Real_Matrix (tk, m1, n1), int i1, int j1) -< !ref > g0float tk
extern fn {tk : tkind}
Real_Matrix_set_at :
{m1, n1 : pos}
{i1, j1 : pos | i1 <= m1; j1 <= n1}
(Real_Matrix (tk, m1, n1), int i1, int j1, g0float tk) -< !refwrt >
void
extern fn {}
Real_Matrix_reflect_lower_triangle :
(* This operation makes every It is a change in how INDEXING
works. All the storage is still in the lower triangle. *)
{tk : tkind}
{n1 : pos}
{m0, n0 : pos}
Real_Matrix (tk, n1, n1, m0, n0) -<>
Real_Matrix (tk, n1, n1, m0, n0)
extern fn {tk : tkind}
Real_Matrix_fprint :
{m, n : pos}
(FILEref, Real_Matrix (tk, m, n)) -<1> void
overload copy with Real_Matrix_copy
overload copy_to with Real_Matrix_copy_to
overload dimension with Real_Matrix_dimension
overload [] with Real_Matrix_get_at
overload [] with Real_Matrix_set_at
overload reflect_lower_triangle with
Real_Matrix_reflect_lower_triangle
(*------------------------------------------------------------------*)
(* Implementation of the "very little matrix library" *)
implement {tk}
Real_Matrix_make_elt (m0, n0, elt) =
Real_Matrix (matrixref_make_elt<g0float tk> (i2sz m0, i2sz n0, elt),
m0, n0, m0, n0, lam (i1, j1) => @(i1, j1))
implement {}
Real_Matrix_dimension A =
case+ A of Real_Matrix (_, m1, n1, _, _, _) => @(m1, n1)
implement {tk}
Real_Matrix_get_at (A, i1, j1) =
let
val+ Real_Matrix (storage, _, _, _, n0, index_map) = A
val @(i0, j0) = index_map (i1, j1)
in
matrixref_get_at<g0float tk> (storage, pred i0, n0, pred j0)
end
implement {tk}
Real_Matrix_set_at (A, i1, j1, x) =
let
val+ Real_Matrix (storage, _, _, _, n0, index_map) = A
val @(i0, j0) = index_map (i1, j1)
in
matrixref_set_at<g0float tk> (storage, pred i0, n0, pred j0, x)
end
implement {}
Real_Matrix_reflect_lower_triangle {..} {n1} A =
let
typedef t = intBtwe (1, n1)
val+ Real_Matrix (storage, n1, _, m0, n0, index_map) = A
in
Real_Matrix (storage, n1, n1, m0, n0,
lam (i, j) =>
index_map ((if j <= i then i else j) : t,
(if j <= i then j else i) : t))
end
implement {tk}
Real_Matrix_copy A =
let
val @(m1, n1) = dimension A
val C = Real_Matrix_make_elt<tk> (m1, n1, A[1, 1])
val () = copy_to<tk> (C, A)
in
C
end
implement {tk}
Real_Matrix_copy_to (Dst, Src) =
let
val @(m1, n1) = dimension Src
prval [m1 : int] EQINT () = eqint_make_gint m1
prval [n1 : int] EQINT () = eqint_make_gint n1
var i : intGte 1
in
for* {i : pos | i <= m1 + 1} .<(m1 + 1) - i>.
(i : int i) =>
(i := 1; i <> succ m1; i := succ i)
let
var j : intGte 1
in
for* {j : pos | j <= n1 + 1} .<(n1 + 1) - j>.
(j : int j) =>
(j := 1; j <> succ n1; j := succ j)
Dst[i, j] := Src[i, j]
end
end
implement {tk}
Real_Matrix_fprint {m, n} (outf, A) =
let
val @(m, n) = dimension A
var i : intGte 1
in
for* {i : pos | i <= m + 1} .<(m + 1) - i>.
(i : int i) =>
(i := 1; i <> succ m; i := succ i)
let
var j : intGte 1
in
for* {j : pos | j <= n + 1} .<(n + 1) - j>.
(j : int j) =>
(j := 1; j <> succ n; j := succ j)
let
typedef FILEstar = $extype"FILE *"
extern castfn FILEref2star : FILEref -<> FILEstar
val _ = $extfcall (int, "fprintf", FILEref2star outf,
"%16.6g", A[i, j])
in
end;
fprintln! (outf)
end
end
(*------------------------------------------------------------------*)
(* Cholesky-Banachiewicz, in place. See
https://en.wikipedia.org/w/index.php?title=Cholesky_decomposition&oldid=1149960985#The_Cholesky%E2%80%93Banachiewicz_and_Cholesky%E2%80%93Crout_algorithms
I would use Cholesky-Crout if my matrices were stored in column
major order. But it makes little difference. *)
extern fn {tk : tkind}
Real_Matrix_cholesky_decomposition :
(* Only the lower triangle is considered. *)
{n : pos}
Real_Matrix (tk, n, n) -< !refwrt > void
overload cholesky_decomposition with
Real_Matrix_cholesky_decomposition
implement {tk}
Real_Matrix_cholesky_decomposition {n} A =
let
val @(n, _) = dimension A
(* I arrange the nested loops somewhat differently from how it is
done in the Wikipedia article's C snippet. *)
fun
repeat {i, j : pos | j <= i; i <= n + 1} (* <-- allowed values *)
.<(n + 1) - i, i - j>. (* <-- proof of termination *)
(i : int i, j : int j) :<!refwrt> void =
if i = n + 1 then
() (* All done. *)
else
let
fun
_sum {k : pos | k <= j} .<j - k>.
(x : g0float tk, k : int k) :<!refwrt> g0float tk =
if k = j then
x
else
_sum (x + (A[i, k] * A[j, k]), succ k)
val sum = _sum (Zero, 1)
in
if j = i then
begin
A[i, j] := sqrt (A[i, i] - sum);
repeat (succ i, 1)
end
else
begin
A[i, j] := (One / A[j, j]) * (A[i, j] - sum);
repeat (i, succ j)
end
end
in
repeat (1, 1)
end
(*------------------------------------------------------------------*)
fn {tk : tkind} (* We like Fortran, so COLUMN major here. *)
column_major_list_to_square_matrix
{n : pos}
(n : int n,
lst : list (g0float tk, n * n))
: Real_Matrix (tk, n, n) =
let
#define :: list_cons
prval () = mul_gte_gte_gte {n, n} ()
val A = Real_Matrix_make_elt (n, n, NAN)
val lstref : ref (List0 (g0float tk)) = ref lst
var j : intGte 1
in
for* {j : pos | j <= n + 1} .<(n + 1) - j>.
(j : int j) =>
(j := 1; j <> succ n; j := succ j)
let
var i : intGte 1
in
for* {i : pos | i <= n + 1} .<(n + 1) - i>.
(i : int i) =>
(i := 1; i <> succ n; i := succ i)
case- !lstref of
| hd :: tl =>
begin
A[i, j] := hd;
!lstref := tl
end
end;
A
end
implement
main0 () =
let
val _A =
column_major_list_to_square_matrix
(3, $list (25.0, 15.0, ~5.0,
0.0, 18.0, 0.0,
0.0, 0.0, 11.0))
val A = reflect_lower_triangle _A
and B = copy _A
val () =
begin
cholesky_decomposition B;
print! ("\nThe Cholesky decomposition of\n\n");
Real_Matrix_fprint (stdout_ref, A);
print! ("is\n");
Real_Matrix_fprint (stdout_ref, B)
end
val _A =
column_major_list_to_square_matrix
(4, $list (18.0, 22.0, 54.0, 42.0,
0.0, 70.0, 86.0, 62.0,
0.0, 0.0, 174.0, 134.0,
0.0, 0.0, 0.0, 106.0))
val A = reflect_lower_triangle _A
and B = copy _A
val () =
begin
cholesky_decomposition B;
print! ("\nThe Cholesky decomposition of\n\n");
Real_Matrix_fprint (stdout_ref, A);
print! ("is\n");
Real_Matrix_fprint (stdout_ref, B)
end
in
println! ()
end
(*------------------------------------------------------------------*)
</syntaxhighlight>
{{out}}
<pre>$ patscc -std=gnu2x -g -O2 -DATS_MEMALLOC_GCBDW cholesky_decomposition_task.dats -lgc -lm && ./a.out
The Cholesky decomposition of
25 15 -5
15 18 0
-5 0 11
is
5 0 0
3 3 0
-1 1 3
The Cholesky decomposition of
18 22 54 42
22 70 86 62
54 86 174 134
42 62 134 106
is
4.24264 0 0 0
5.18545 6.56591 0 0
12.7279 3.04604 1.64974 0
9.89949 1.62455 1.84971 1.39262
</pre>
=={{header|AutoHotkey}}==
<syntaxhighlight lang="autohotkey">Cholesky_Decomposition(A){
Line 367 ⟶ 780:
,[9.899, 1.625, 1.850, 1.393]]
----</pre>
=={{header|BBC BASIC}}==
{{works with|BBC BASIC for Windows}}
Line 1,209 ⟶ 1,623:
12.72792 3.04604 1.64974 0.00000
9.89949 1.62455 1.84971 1.39262</pre>
=={{header|Frink}}==
Frink's package [https://frinklang.org/fsp/colorize.fsp?f=Matrix.frink Matrix.frink] contains routines for Cholesky-Crout decomposition of a square Hermitian matrix (which can be real or complex.) This code is adapted from that more powerful class to work on raw 2-dimensional arrays. This also demonstrates Frink's layout routines.
<syntaxhighlight lang="frink">Cholesky[array] :=
{
n = length[array]
L = new array[[n,n], 0]
for j = 0 to n-1
{
sum = 0
for k = 0 to j-1
sum = sum + (L@j@k)^2
L@j@j = sqrt[array@j@j - sum]
for i = j+1 to n-1
{
sum = 0
for k = 0 to j-1
sum = sum + L@i@k * L@j@k
L@i@j = (1 / L@j@j * (array@i@j -sum))
}
}
return L
}
A = [[ 25, 15, -5],
[ 15, 18, 0],
[ -5, 0, 11]]
println[formatTable[[[formatMatrix[A], "->", formatMatrix[Cholesky[A]]]]]]
B = [[18, 22, 54, 42],
[22, 70, 86, 62],
[54, 86, 174, 134],
[42, 62, 134, 106]]
println[formatTable[[[formatMatrix[B], "->", formatMatrix[formatFix[Cholesky[B], 1, 5]]]]]]</syntaxhighlight>
{{out}}
<pre>
┌ ┐ ┌ ┐
│25 15 -5│ │ 5 0 0│
│ │ │ │
│15 18 0│ -> │ 3 3 0│
│ │ │ │
│-5 0 11│ │-1 1 3│
└ ┘ └ ┘
┌ ┐ ┌ ┐
│18 22 54 42│ │ 4.24264 0.00000 0.00000 0.00000│
│ │ │ │
│22 70 86 62│ │ 5.18545 6.56591 0.00000 0.00000│
│ │ -> │ │
│54 86 174 134│ │12.72792 3.04604 1.64974 0.00000│
│ │ │ │
│42 62 134 106│ │ 9.89949 1.62455 1.84971 1.39262│
└ ┘ └ ┘
</pre>
=={{header|Go}}==
===Real===
Line 3,228 ⟶ 3,704:
class Matrix
def cholesky_factor
raise ArgumentError, "must provide symmetric matrix" unless symmetric?
Line 3,271 ⟶ 3,737:
[9.899494936611665, 1.6245538642137891, 1.849711005231382, 1.3926212476455924]]
</pre>
=={{header|Rust}}==
Line 3,759 ⟶ 4,226:
End Function
</syntaxhighlight>
=={{header|V (Vlang)}}==
{{trans|go}}
<syntaxhighlight lang="v (vlang)">import math
// Symmetric and Lower use a packed representation that stores only
Line 3,882 ⟶ 4,349:
9.89949 1.62455 1.84971 1.3926
</pre>
=={{header|Wren}}==
{{libheader|Wren-matrix}}
{{libheader|Wren-fmt}}
<syntaxhighlight lang="
import "./fmt" for Fmt
var arrays = [
Line 3,931 ⟶ 4,399:
| 9.89949 1.62455 1.84971 1.39262|
</pre>
=={{header|XPL0}}==
<syntaxhighlight lang "XPL0">real L(4*4);
func real Cholesky(A, N);
real A; int N;
real S;
int I, J, K;
[for I:= 0 to N*N-1 do L(I):= 0.;
for I:= 0 to N-1 do
for J:= 0 to I do
[S:= 0.;
for K:= 0 to J-1 do
S:= S + L(I*N+K) * L(J*N+K);
L(I*N+J):= if I = J then sqrt(A(I*N+I) - S)
else (1.0 / L(J*N+J) * (A(I*N+J) - S));
];
return L;
];
proc ShowMatrix(A, N);
real A; int N;
int I, J;
[for I:= 0 to N-1 do
[for J:= 0 to N-1 do
RlOut(0, A(I*N+J));
CrLf(0);
];
];
int N;
real M1, C1, M2, C2;
[N:= 3;
M1:= [25., 15., -5.,
15., 18., 0.,
-5., 0., 11.];
C1:= Cholesky(M1, N);
ShowMatrix(C1, N);
CrLf(0);
N:= 4;
M2:= [18., 22., 54., 42.,
22., 70., 86., 62.,
54., 86., 174., 134.,
42., 62., 134., 106.];
C2:= Cholesky(M2, N);
ShowMatrix(C2, N);
]</syntaxhighlight>
{{out}}
<pre>
5.00000 0.00000 0.00000
3.00000 3.00000 0.00000
-1.00000 1.00000 3.00000
4.24264 0.00000 0.00000 0.00000
5.18545 6.56591 0.00000 0.00000
12.72792 3.04604 1.64974 0.00000
9.89949 1.62455 1.84971 1.39262
</pre>
=={{header|zkl}}==
Using the GNU Scientific Library:
| |||