LU decomposition: Difference between revisions
Content added Content deleted
Thundergnat (talk | contribs) m (syntax highlighting fixup automation) |
(Added ATS.) |
||
Line 463: | Line 463: | ||
0.00 1.00 0.00 0.00 |
0.00 1.00 0.00 0.00 |
||
0.00 0.00 0.00 1.00</pre> |
0.00 0.00 0.00 1.00</pre> |
||
=={{header|ATS}}== |
|||
<syntaxhighlight lang="ATS"> |
|||
(* There is a "little matrix library" included below. Not all of it is |
|||
used, though unused parts may prove useful for playing with the |
|||
code. |
|||
One might, by the way, find interesting how I get the P matrix from |
|||
a permutation vector. *) |
|||
%{^ |
|||
#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 |
|||
(* You can substitute an "fma" function for this definition: *) |
|||
macdef multiply_and_add (x, y, z) = (,(x) * ,(y)) + ,(z) |
|||
exception Exc_degenerate_problem of string |
|||
(*------------------------------------------------------------------*) |
|||
(* A "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 {tk : tkind} |
|||
Real_Matrix_fill_with_elt : |
|||
{m1, n1 : pos} |
|||
(Real_Matrix (tk, m1, n1), g0float tk) -< !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_apply_index_map : |
|||
{tk : tkind} |
|||
{m1, n1 : pos} |
|||
{m0, n0 : pos} |
|||
(Real_Matrix (tk, m0, n0), int m1, int n1, |
|||
Matrix_Index_Map (m1, n1, m0, n0)) -<> |
|||
Real_Matrix (tk, m1, n1) |
|||
extern fn {} |
|||
Real_Matrix_transpose : |
|||
(* This is transposed INDEXING. It does NOT copy the data. *) |
|||
{tk : tkind} |
|||
{m1, n1 : pos} |
|||
{m0, n0 : pos} |
|||
Real_Matrix (tk, m1, n1, m0, n0) -<> |
|||
Real_Matrix (tk, n1, m1, m0, n0) |
|||
extern fn {} |
|||
Real_Matrix_block : |
|||
(* This is block (submatrix) INDEXING. It does NOT copy the data. *) |
|||
{tk : tkind} |
|||
{p0, p1 : pos | p0 <= p1} |
|||
{q0, q1 : pos | q0 <= q1} |
|||
{m1, n1 : pos | p1 <= m1; q1 <= n1} |
|||
{m0, n0 : pos} |
|||
(Real_Matrix (tk, m1, n1, m0, n0), |
|||
int p0, int p1, int q0, int q1) -<> |
|||
Real_Matrix (tk, p1 - p0 + 1, q1 - q0 + 1, m0, n0) |
|||
extern fn {tk : tkind} |
|||
Real_Matrix_unit_matrix : |
|||
{m : pos} |
|||
int m -< !refwrt > Real_Matrix (tk, m, m) |
|||
extern fn {tk : tkind} |
|||
Real_Matrix_unit_matrix_to : |
|||
{m : pos} |
|||
Real_Matrix (tk, m, m) -< !refwrt > void |
|||
extern fn {tk : tkind} |
|||
Real_Matrix_matrix_sum : |
|||
{m, n : pos} |
|||
(Real_Matrix (tk, m, n), Real_Matrix (tk, m, n)) -< !refwrt > |
|||
Real_Matrix (tk, m, n) |
|||
extern fn {tk : tkind} |
|||
Real_Matrix_matrix_sum_to : |
|||
{m, n : pos} |
|||
(Real_Matrix (tk, m, n), (* destination*) |
|||
Real_Matrix (tk, m, n), |
|||
Real_Matrix (tk, m, n)) -< !refwrt > |
|||
void |
|||
extern fn {tk : tkind} |
|||
Real_Matrix_matrix_difference : |
|||
{m, n : pos} |
|||
(Real_Matrix (tk, m, n), Real_Matrix (tk, m, n)) -< !refwrt > |
|||
Real_Matrix (tk, m, n) |
|||
extern fn {tk : tkind} |
|||
Real_Matrix_matrix_difference_to : |
|||
{m, n : pos} |
|||
(Real_Matrix (tk, m, n), (* destination*) |
|||
Real_Matrix (tk, m, n), |
|||
Real_Matrix (tk, m, n)) -< !refwrt > |
|||
void |
|||
extern fn {tk : tkind} |
|||
Real_Matrix_matrix_product : |
|||
{m, n, p : pos} |
|||
(Real_Matrix (tk, m, n), Real_Matrix (tk, n, p)) -< !refwrt > |
|||
Real_Matrix (tk, m, p) |
|||
extern fn {tk : tkind} |
|||
Real_Matrix_matrix_product_to : |
|||
(* For the matrix product, the destination should not be the same as |
|||
either of the other matrices. *) |
|||
{m, n, p : pos} |
|||
(Real_Matrix (tk, m, p), (* destination*) |
|||
Real_Matrix (tk, m, n), |
|||
Real_Matrix (tk, n, p)) -< !refwrt > |
|||
void |
|||
extern fn {tk : tkind} |
|||
Real_Matrix_scalar_product : |
|||
{m, n : pos} |
|||
(Real_Matrix (tk, m, n), g0float tk) -< !refwrt > |
|||
Real_Matrix (tk, m, n) |
|||
extern fn {tk : tkind} |
|||
Real_Matrix_scalar_product_2 : |
|||
{m, n : pos} |
|||
(g0float tk, Real_Matrix (tk, m, n)) -< !refwrt > |
|||
Real_Matrix (tk, m, n) |
|||
extern fn {tk : tkind} |
|||
Real_Matrix_scalar_product : |
|||
{m, n : pos} |
|||
(Real_Matrix (tk, m, n), g0float tk) -< !refwrt > |
|||
Real_Matrix (tk, m, n) |
|||
extern fn {tk : tkind} |
|||
Real_Matrix_scalar_product_2 : |
|||
{m, n : pos} |
|||
(g0float tk, Real_Matrix (tk, m, n)) -< !refwrt > |
|||
Real_Matrix (tk, m, n) |
|||
extern fn {tk : tkind} |
|||
Real_Matrix_scalar_product_to : |
|||
{m, n : pos} |
|||
(Real_Matrix (tk, m, n), (* destination*) |
|||
Real_Matrix (tk, m, n), |
|||
g0float tk) -< !refwrt > |
|||
void |
|||
extern fn {tk : tkind} (* Useful for debugging. *) |
|||
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 fill_with_elt with Real_Matrix_fill_with_elt |
|||
overload dimension with Real_Matrix_dimension |
|||
overload [] with Real_Matrix_get_at |
|||
overload [] with Real_Matrix_set_at |
|||
overload apply_index_map with Real_Matrix_apply_index_map |
|||
overload transpose with Real_Matrix_transpose |
|||
overload block with Real_Matrix_block |
|||
overload unit_matrix with Real_Matrix_unit_matrix |
|||
overload unit_matrix_to with Real_Matrix_unit_matrix_to |
|||
overload matrix_sum with Real_Matrix_matrix_sum |
|||
overload matrix_sum_to with Real_Matrix_matrix_sum_to |
|||
overload matrix_difference with Real_Matrix_matrix_difference |
|||
overload matrix_difference_to with Real_Matrix_matrix_difference_to |
|||
overload matrix_product with Real_Matrix_matrix_product |
|||
overload matrix_product_to with Real_Matrix_matrix_product_to |
|||
overload scalar_product with Real_Matrix_scalar_product |
|||
overload scalar_product with Real_Matrix_scalar_product_2 |
|||
overload scalar_product_to with Real_Matrix_scalar_product_to |
|||
overload + with matrix_sum |
|||
overload - with matrix_difference |
|||
overload * with matrix_product |
|||
overload * with scalar_product |
|||
(*------------------------------------------------------------------*) |
|||
(* Implementation of the "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_apply_index_map (A, m1, n1, index_map) = |
|||
(* This is not the most efficient way to acquire new indexing, but |
|||
it will work. It requires three closures, instead of the two |
|||
needed by our implementations of "transpose" and "block". *) |
|||
let |
|||
val+ Real_Matrix (storage, m1a, n1a, m0, n0, index_map_1a) = A |
|||
in |
|||
Real_Matrix (storage, m1, n1, m0, n0, |
|||
lam (i1, j1) => |
|||
index_map_1a (i1a, j1a) where |
|||
{ val @(i1a, j1a) = index_map (i1, j1) }) |
|||
end |
|||
implement {} |
|||
Real_Matrix_transpose A = |
|||
let |
|||
val+ Real_Matrix (storage, m1, n1, m0, n0, index_map) = A |
|||
in |
|||
Real_Matrix (storage, n1, m1, m0, n0, |
|||
lam (i1, j1) => index_map (j1, i1)) |
|||
end |
|||
implement {} |
|||
Real_Matrix_block (A, p0, p1, q0, q1) = |
|||
let |
|||
val+ Real_Matrix (storage, m1, n1, m0, n0, index_map) = A |
|||
in |
|||
Real_Matrix (storage, succ (p1 - p0), succ (q1 - q0), m0, n0, |
|||
lam (i1, j1) => |
|||
index_map (p0 + pred i1, q0 + pred j1)) |
|||
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_fill_with_elt (A, elt) = |
|||
let |
|||
val @(m1, n1) = dimension A |
|||
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) |
|||
A[i, j] := elt |
|||
end |
|||
end |
|||
implement {tk} |
|||
Real_Matrix_unit_matrix {m} m = |
|||
let |
|||
val A = Real_Matrix_make_elt<tk> (m, m, Zero) |
|||
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) |
|||
A[i, i] := One; |
|||
A |
|||
end |
|||
implement {tk} |
|||
Real_Matrix_unit_matrix_to A = |
|||
let |
|||
val @(m, _) = dimension A |
|||
prval [m : int] EQINT () = eqint_make_gint m |
|||
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 <= m + 1} .<(m + 1) - j>. |
|||
(j : int j) => |
|||
(j := 1; j <> succ m; j := succ j) |
|||
A[i, j] := (if i = j then One else Zero) |
|||
end |
|||
end |
|||
implement {tk} |
|||
Real_Matrix_matrix_sum (A, B) = |
|||
let |
|||
val @(m, n) = dimension A |
|||
val C = Real_Matrix_make_elt<tk> (m, n, NAN) |
|||
val () = matrix_sum_to<tk> (C, A, B) |
|||
in |
|||
C |
|||
end |
|||
implement {tk} |
|||
Real_Matrix_matrix_sum_to (C, A, B) = |
|||
let |
|||
val @(m, n) = dimension A |
|||
prval [m : int] EQINT () = eqint_make_gint m |
|||
prval [n : int] EQINT () = eqint_make_gint n |
|||
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) |
|||
C[i, j] := A[i, j] + B[i, j] |
|||
end |
|||
end |
|||
implement {tk} |
|||
Real_Matrix_matrix_difference (A, B) = |
|||
let |
|||
val @(m, n) = dimension A |
|||
val C = Real_Matrix_make_elt<tk> (m, n, NAN) |
|||
val () = matrix_difference_to<tk> (C, A, B) |
|||
in |
|||
C |
|||
end |
|||
implement {tk} |
|||
Real_Matrix_matrix_difference_to (C, A, B) = |
|||
let |
|||
val @(m, n) = dimension A |
|||
prval [m : int] EQINT () = eqint_make_gint m |
|||
prval [n : int] EQINT () = eqint_make_gint n |
|||
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) |
|||
C[i, j] := A[i, j] - B[i, j] |
|||
end |
|||
end |
|||
implement {tk} |
|||
Real_Matrix_matrix_product (A, B) = |
|||
let |
|||
val @(m, n) = dimension A and @(_, p) = dimension B |
|||
val C = Real_Matrix_make_elt<tk> (m, p, NAN) |
|||
val () = matrix_product_to<tk> (C, A, B) |
|||
in |
|||
C |
|||
end |
|||
implement {tk} |
|||
Real_Matrix_matrix_product_to (C, A, B) = |
|||
let |
|||
val @(m, n) = dimension A and @(_, p) = dimension B |
|||
prval [m : int] EQINT () = eqint_make_gint m |
|||
prval [n : int] EQINT () = eqint_make_gint n |
|||
prval [p : int] EQINT () = eqint_make_gint p |
|||
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 k : intGte 1 |
|||
in |
|||
for* {k : pos | k <= p + 1} .<(p + 1) - k>. |
|||
(k : int k) => |
|||
(k := 1; k <> succ p; k := succ k) |
|||
let |
|||
var j : intGte 1 |
|||
in |
|||
C[i, k] := A[i, 1] * B[1, k]; |
|||
for* {j : pos | j <= n + 1} .<(n + 1) - j>. |
|||
(j : int j) => |
|||
(j := 2; j <> succ n; j := succ j) |
|||
C[i, k] := |
|||
multiply_and_add (A[i, j], B[j, k], C[i, k]) |
|||
end |
|||
end |
|||
end |
|||
implement {tk} |
|||
Real_Matrix_scalar_product (A, r) = |
|||
let |
|||
val @(m, n) = dimension A |
|||
val C = Real_Matrix_make_elt<tk> (m, n, NAN) |
|||
val () = scalar_product_to<tk> (C, A, r) |
|||
in |
|||
C |
|||
end |
|||
implement {tk} |
|||
Real_Matrix_scalar_product_2 (r, A) = |
|||
Real_Matrix_scalar_product<tk> (A, r) |
|||
implement {tk} |
|||
Real_Matrix_scalar_product_to (C, A, r) = |
|||
let |
|||
val @(m, n) = dimension A |
|||
prval [m : int] EQINT () = eqint_make_gint m |
|||
prval [n : int] EQINT () = eqint_make_gint n |
|||
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) |
|||
C[i, j] := A[i, j] * r |
|||
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, |
|||
"%12.6lf", A[i, j]) |
|||
in |
|||
end; |
|||
fprintln! (outf) |
|||
end |
|||
end |
|||
(*------------------------------------------------------------------*) |
|||
(* LUP decomposition. Based on |
|||
https://en.wikipedia.org/w/index.php?title=LU_decomposition&oldid=1146366204#C_code_example |
|||
*) |
|||
extern fn {tk : tkind} |
|||
Real_Matrix_LUP_decomposition : |
|||
{n : pos} |
|||
(Real_Matrix (tk, n, n), |
|||
g0float tk (* tolerance *) ) -< !exnrefwrt > |
|||
@(Real_Matrix (tk, n, n), |
|||
Real_Matrix (tk, n, n), |
|||
Real_Matrix (tk, n, n)) |
|||
overload LUP_decomposition with Real_Matrix_LUP_decomposition |
|||
implement {tk} |
|||
Real_Matrix_LUP_decomposition {n} (A, tol) = |
|||
let |
|||
val @(n, _) = dimension A |
|||
typedef one_to_n = intBtwe (1, n) |
|||
(* The initial permutation is [1,2,3,...,n]. *) |
|||
implement |
|||
array_tabulate$fopr<one_to_n> i = |
|||
let |
|||
val i = g1ofg0 (sz2i (succ i)) |
|||
val () = assertloc ((1 <= i) * (i <= n)) |
|||
in |
|||
i |
|||
end |
|||
val permutation = |
|||
$effmask_all arrayref_tabulate<one_to_n> (i2sz n) |
|||
fn |
|||
index_map : Matrix_Index_Map (n, n, n, n) = |
|||
lam (i1, j1) => $effmask_ref |
|||
(@(i0, j1) where { val i0 = permutation[i1 - 1] }) |
|||
val A = apply_index_map (copy<tk> A, n, n, index_map) |
|||
fun |
|||
select_pivot {i, k : pos | i <= k; k <= n + 1} |
|||
.<(n + 1) - k>. |
|||
(i : int i, |
|||
k : int k, |
|||
max_abs : g0float tk, |
|||
k_max_abs : intBtwe (i, n)) |
|||
:<!ref> @(g0float tk, intBtwe (i, n)) = |
|||
if k = succ n then |
|||
@(max_abs, k_max_abs) |
|||
else |
|||
let |
|||
val absval = A[k, i] |
|||
in |
|||
if absval > max_abs then |
|||
select_pivot (i, succ k, absval, k) |
|||
else |
|||
select_pivot (i, succ k, max_abs, k_max_abs) |
|||
end |
|||
fn {} |
|||
exchange_rows (i1 : one_to_n, |
|||
i2 : one_to_n) :<!refwrt> void = |
|||
if i1 <> i2 then |
|||
let |
|||
val k1 = permutation[pred i1] |
|||
and k2 = permutation[pred i2] |
|||
in |
|||
permutation[pred i1] := k2; |
|||
permutation[pred i2] := k1 |
|||
end |
|||
val () = |
|||
let |
|||
var i : Int |
|||
in |
|||
for* {i : pos | i <= n + 1} .<(n + 1) - i>. |
|||
(i : int i) => |
|||
(i := 1; i <> succ n; i := succ i) |
|||
let |
|||
val @(maxabs, i_pivot) = select_pivot (i, i, Zero, i) |
|||
prval [i_pivot : int] EQINT () = eqint_make_gint i_pivot |
|||
var j : Int |
|||
in |
|||
if maxabs < tol then |
|||
$raise Exc_degenerate_problem |
|||
("Real_Matrix_LUP_decomposition"); |
|||
exchange_rows (i_pivot, i); |
|||
for* {j : int | i + 1 <= j; j <= n + 1} |
|||
.<(n + 1) - j>. |
|||
(j : int j) => |
|||
(j := succ i; j <> succ n; j := succ j) |
|||
let |
|||
var k : Int |
|||
in |
|||
A[j, i] := A[j, i] / A[i, i]; |
|||
for* {k : int | i + 1 <= k; k <= n + 1} |
|||
.<(n + 1) - k>. |
|||
(k : int k) => |
|||
(k := succ i; k <> succ n; k := succ k) |
|||
A[j, k] := |
|||
multiply_and_add |
|||
(~A[j, i], A[i, k], A[j, k]) |
|||
end |
|||
end |
|||
end |
|||
val U = A |
|||
val L = Real_Matrix_unit_matrix<tk> n |
|||
val () = |
|||
let |
|||
var i : Int |
|||
in |
|||
for* {i : int | 2 <= i; i <= n + 1} .<(n + 1) - i>. |
|||
(i : int i) => |
|||
(i := 2; i <> succ n; i := succ i) |
|||
let |
|||
var j : Int |
|||
in |
|||
for* {j : pos | j <= i} .<i - j>. |
|||
(j : int j) => |
|||
(j := 1; j <> i; j := succ j) |
|||
begin |
|||
L[i, j] := U[i, j]; |
|||
U[i, j] := Zero |
|||
end |
|||
end |
|||
end |
|||
val P = apply_index_map (Real_Matrix_unit_matrix<tk> n, |
|||
n, n, index_map) |
|||
in |
|||
@(L, U, P) |
|||
end |
|||
(*------------------------------------------------------------------*) |
|||
implement |
|||
main0 () = |
|||
(* I use tolerances of zero, secure in the knowledge that IEEE |
|||
floating point will not crash the program just because a matrix |
|||
was singular. :) *) |
|||
let |
|||
val A = Real_Matrix_make_elt<dblknd> (3, 3, NAN) |
|||
val () = |
|||
(A[1, 1] := 1.0; A[1, 2] := 3.0; A[1, 3] := 5.0; |
|||
A[2, 1] := 2.0; A[2, 2] := 4.0; A[2, 3] := 7.0; |
|||
A[3, 1] := 1.0; A[3, 2] := 1.0; A[3, 3] := 0.0) |
|||
val @(L, U, P) = LUP_decomposition (A, 0.0) |
|||
val () = println! "A" |
|||
val () = Real_Matrix_fprint (stdout_ref, A) |
|||
val () = println! "L" |
|||
val () = Real_Matrix_fprint (stdout_ref, L) |
|||
val () = println! "U" |
|||
val () = Real_Matrix_fprint (stdout_ref, U) |
|||
val () = println! "P" |
|||
val () = Real_Matrix_fprint (stdout_ref, P) |
|||
val () = println! "PA - LU" |
|||
val () = Real_Matrix_fprint (stdout_ref, P * A - L * U) |
|||
val () = println! "\n------------------------------------------\n" |
|||
val A = Real_Matrix_make_elt<dblknd> (4, 4, NAN) |
|||
val () = |
|||
(A[1, 1] := 11.0; A[1, 2] := 9.0; A[1, 3] := 24.0; A[1, 4] := 2.0; |
|||
A[2, 1] := 1.0; A[2, 2] := 5.0; A[2, 3] := 2.0; A[2, 4] := 6.0; |
|||
A[3, 1] := 3.0; A[3, 2] := 17.0; A[3, 3] := 18.0; A[3, 4] := 1.0; |
|||
A[4, 1] := 2.0; A[4, 2] := 5.0; A[4, 3] := 7.0; A[4, 4] := 1.0) |
|||
val @(L, U, P) = LUP_decomposition (A, 0.0) |
|||
val () = println! "A" |
|||
val () = Real_Matrix_fprint (stdout_ref, A) |
|||
val () = println! "L" |
|||
val () = Real_Matrix_fprint (stdout_ref, L) |
|||
val () = println! "U" |
|||
val () = Real_Matrix_fprint (stdout_ref, U) |
|||
val () = println! "P" |
|||
val () = Real_Matrix_fprint (stdout_ref, P) |
|||
val () = println! "PA - LU" |
|||
val () = Real_Matrix_fprint (stdout_ref, P * A - L * U) |
|||
val () = println! "\n------------------------------------------\n" |
|||
val () = println! ("I have added an example having an asymmetric", |
|||
" permutation\nmatrix, ", |
|||
"because I needed one for testing:") |
|||
val () = println! () |
|||
val A = Real_Matrix_make_elt<dblknd> (4, 4, NAN) |
|||
val () = |
|||
(A[1, 1] := 11.0; A[1, 2] := 9.0; A[1, 3] := 24.0; A[1, 4] := 2.0; |
|||
A[2, 1] := 1.0; A[2, 2] := 5.0; A[2, 3] := 2.0; A[2, 4] := 6.0; |
|||
A[3, 1] := 3.0; A[3, 2] := 175.0; A[3, 3] := 18.0; A[3, 4] := 1.0; |
|||
A[4, 1] := 2.0; A[4, 2] := 5.0; A[4, 3] := 7.0; A[4, 4] := 1.0) |
|||
val @(L, U, P) = LUP_decomposition (A, 0.0) |
|||
val () = println! "A" |
|||
val () = Real_Matrix_fprint (stdout_ref, A) |
|||
val () = println! "L" |
|||
val () = Real_Matrix_fprint (stdout_ref, L) |
|||
val () = println! "U" |
|||
val () = Real_Matrix_fprint (stdout_ref, U) |
|||
val () = println! "P" |
|||
val () = Real_Matrix_fprint (stdout_ref, P) |
|||
val () = println! "PA - LU" |
|||
val () = Real_Matrix_fprint (stdout_ref, P * A - L * U) |
|||
in |
|||
end |
|||
(*------------------------------------------------------------------*) |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre style="font-size:90%">$ patscc -std=gnu2x -g -O2 -march=native -DATS_MEMALLOC_GCBDW lu_decomposition_task.dats -lgc && ./a.out |
|||
A |
|||
1.000000 3.000000 5.000000 |
|||
2.000000 4.000000 7.000000 |
|||
1.000000 1.000000 0.000000 |
|||
L |
|||
1.000000 0.000000 0.000000 |
|||
0.500000 1.000000 0.000000 |
|||
0.500000 -1.000000 1.000000 |
|||
U |
|||
2.000000 4.000000 7.000000 |
|||
0.000000 1.000000 1.500000 |
|||
0.000000 0.000000 -2.000000 |
|||
P |
|||
0.000000 1.000000 0.000000 |
|||
1.000000 0.000000 0.000000 |
|||
0.000000 0.000000 1.000000 |
|||
PA - LU |
|||
0.000000 0.000000 0.000000 |
|||
0.000000 0.000000 0.000000 |
|||
0.000000 0.000000 0.000000 |
|||
------------------------------------------ |
|||
A |
|||
11.000000 9.000000 24.000000 2.000000 |
|||
1.000000 5.000000 2.000000 6.000000 |
|||
3.000000 17.000000 18.000000 1.000000 |
|||
2.000000 5.000000 7.000000 1.000000 |
|||
L |
|||
1.000000 0.000000 0.000000 0.000000 |
|||
0.272727 1.000000 0.000000 0.000000 |
|||
0.090909 0.287500 1.000000 0.000000 |
|||
0.181818 0.231250 0.003597 1.000000 |
|||
U |
|||
11.000000 9.000000 24.000000 2.000000 |
|||
0.000000 14.545455 11.454545 0.454545 |
|||
0.000000 0.000000 -3.475000 5.687500 |
|||
0.000000 0.000000 0.000000 0.510791 |
|||
P |
|||
1.000000 0.000000 0.000000 0.000000 |
|||
0.000000 0.000000 1.000000 0.000000 |
|||
0.000000 1.000000 0.000000 0.000000 |
|||
0.000000 0.000000 0.000000 1.000000 |
|||
PA - LU |
|||
0.000000 0.000000 0.000000 0.000000 |
|||
0.000000 0.000000 0.000000 0.000000 |
|||
0.000000 0.000000 0.000000 0.000000 |
|||
0.000000 0.000000 0.000000 0.000000 |
|||
------------------------------------------ |
|||
I have added an example having an asymmetric permutation |
|||
matrix, because I needed one for testing: |
|||
A |
|||
11.000000 9.000000 24.000000 2.000000 |
|||
1.000000 5.000000 2.000000 6.000000 |
|||
3.000000 175.000000 18.000000 1.000000 |
|||
2.000000 5.000000 7.000000 1.000000 |
|||
L |
|||
1.000000 0.000000 0.000000 0.000000 |
|||
0.272727 1.000000 0.000000 0.000000 |
|||
0.181818 0.019494 1.000000 0.000000 |
|||
0.090909 0.024236 -0.190393 1.000000 |
|||
U |
|||
11.000000 9.000000 24.000000 2.000000 |
|||
0.000000 172.545455 11.454545 0.454545 |
|||
0.000000 0.000000 2.413066 0.627503 |
|||
0.000000 0.000000 0.000000 5.926638 |
|||
P |
|||
1.000000 0.000000 0.000000 0.000000 |
|||
0.000000 0.000000 1.000000 0.000000 |
|||
0.000000 0.000000 0.000000 1.000000 |
|||
0.000000 1.000000 0.000000 0.000000 |
|||
PA - LU |
|||
0.000000 0.000000 0.000000 0.000000 |
|||
0.000000 0.000000 0.000000 0.000000 |
|||
0.000000 0.000000 0.000000 0.000000 |
|||
0.000000 0.000000 0.000000 0.000000 |
|||
</pre> |
|||
=={{header|AutoHotkey}}== |
=={{header|AutoHotkey}}== |