Gaussian elimination: Difference between revisions
Content added Content deleted
m (→{{header|Wren}}: Wren-trait -> Wren-iterate) |
(Added ATS.) |
||
Line 537: | Line 537: | ||
<pre> |
<pre> |
||
( -.010000000000002, 1.602790394502130, -1.613203059905640, 1.245494121371510, -.490989719584686, .065760696175236) |
( -.010000000000002, 1.602790394502130, -1.613203059905640, 1.245494121371510, -.490989719584686, .065760696175236) |
||
</pre> |
|||
=={{header|ATS}}== |
|||
This program was written by modifying [[Gauss-Jordan_matrix_inversion#ATS]]. There is a commented out portion of the code, whose removal makes this "Gaussian" elimination (with back substitution) rather than "Gauss-Jordan" elimination (without the need for back substitution). |
|||
<syntaxhighlight lang="ats"> |
|||
(* There is a "little matrix library" in the code below. Not all of it |
|||
will be used, but it travels from task to task. Furthermore, the |
|||
"unused" parts are useful during the debugging phase. Also, reading |
|||
them may make it easier to understand other parts of the code. (For |
|||
instance, seeing how "block" and "transpose" work may help one |
|||
understand how "apply_index_map" works.) *) |
|||
(* Set to 1 for debugging: to fill in ones and zeros and other values |
|||
that are not actually used but are part of the theory of the |
|||
Gaussian elimination algorithm. *) |
|||
#define DO_THINGS_THAT_DO_NOT_NEED_TO_BE_DONE 0 |
|||
(* Setting this to 1 may cause rounding to change, and the change in |
|||
rounding is not unlikely to cause detection of singularity of a |
|||
matrix to change. (To invert a matrix that might be nearly |
|||
singular, the SVD seems a popular method.) The |
|||
-fexpensive-optimizations option to GCC also may cause the same |
|||
rounding changes (due to fused-multiply-and-add instructions being |
|||
generated). *) |
|||
#define USE_MULTIPLY_AND_ADD 0 |
|||
%{^ |
|||
#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 |
|||
macdef Two = g0i2f 2 |
|||
#if USE_MULTIPLY_AND_ADD #then |
|||
(* "fma" from the C math library, although your system may have it as |
|||
a built-in. *) |
|||
extern fn {tk : tkind} g0float_fma : (g0float tk, g0float tk, g0float tk) -<> g0float tk |
|||
implement g0float_fma<fltknd> (x, y, z) = $extfcall (float, "fmaf", x, y, z) |
|||
implement g0float_fma<dblknd> (x, y, z) = $extfcall (double, "fma", x, y, z) |
|||
implement g0float_fma<ldblknd> (x, y, z) = $extfcall (ldouble, "fmal", x, y, z) |
|||
overload fma with g0float_fma |
|||
macdef multiply_and_add = fma |
|||
#else |
|||
macdef multiply_and_add (x, y, z) = (,(x) * ,(y)) + ,(z) |
|||
#endif |
|||
(*------------------------------------------------------------------*) |
|||
(* 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} |
|||
Real_Matrix_scalar_multiply_and_add_to : |
|||
{m, n : pos} |
|||
(Real_Matrix (tk, m, n), (* destination*) |
|||
Real_Matrix (tk, m, n), |
|||
g0float tk, |
|||
Real_Matrix (tk, m, n)) -< !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 scalar_multiply_and_add_to with |
|||
Real_Matrix_scalar_multiply_and_add_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_scalar_multiply_and_add_to (C, A, r, 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] := multiply_and_add (A[i, j], r, B[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 |
|||
(*------------------------------------------------------------------*) |
|||
(* Gaussian elimination *) |
|||
extern fn {tk : tkind} |
|||
Real_Matrix_gaussian_elimination : |
|||
(* Solve k systems of linear equations in n unknowns. (A special |
|||
case is if the second argument is a unit matrix. In that case, |
|||
the solution columns constitute the matrix inverse of the first |
|||
argument. See also the Gauss-Jordan matrix inversion task.) *) |
|||
{n, k : pos} |
|||
(Real_Matrix (tk, n, n), Real_Matrix (tk, n, k)) -< !refwrt > |
|||
Option (Real_Matrix (tk, n, k)) |
|||
#if DO_THINGS_THAT_DO_NOT_NEED_TO_BE_DONE #then |
|||
macdef do_needless_things = true |
|||
#else |
|||
macdef do_needless_things = false |
|||
#endif |
|||
fn {tk : tkind} |
|||
needlessly_set_to_value |
|||
{n : pos} |
|||
{i, j : pos | i <= n; j <= n} |
|||
(A : Real_Matrix (tk, n, n), |
|||
i : int i, |
|||
j : int j, |
|||
x : g0float tk) :<!refwrt> void = |
|||
if do_needless_things then A[i, j] := x |
|||
implement {tk} |
|||
Real_Matrix_gaussian_elimination {n, k} (A, B) = |
|||
let |
|||
val @(n, k) = dimension B |
|||
typedef one_to_n = intBtwe (1, n) |
|||
(* Partial pivoting. *) |
|||
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 rows_permutation = |
|||
$effmask_all arrayref_tabulate<one_to_n> (i2sz n) |
|||
fn |
|||
index_map_A : Matrix_Index_Map (n, n, n, n) = |
|||
lam (i1, j1) => $effmask_ref |
|||
(@(i0, j1) where { val i0 = rows_permutation[i1 - 1] }) |
|||
fn |
|||
index_map_B : Matrix_Index_Map (n, k, n, k) = |
|||
lam (i1, j1) => $effmask_ref |
|||
(@(i0, j1) where { val i0 = rows_permutation[i1 - 1] }) |
|||
val A = apply_index_map (copy<tk> A, n, n, index_map_A) |
|||
and B = apply_index_map (copy<tk> B, n, k, index_map_B) |
|||
fn {} |
|||
exchange_rows (i1 : one_to_n, |
|||
i2 : one_to_n) :<!refwrt> void = |
|||
if i1 <> i2 then |
|||
let |
|||
val k1 = rows_permutation[pred i1] |
|||
and k2 = rows_permutation[pred i2] |
|||
in |
|||
rows_permutation[pred i1] := k2; |
|||
rows_permutation[pred i2] := k1 |
|||
end |
|||
fn {} |
|||
normalize_pivot_row (j : one_to_n) :<!refwrt> void = |
|||
let |
|||
prval [j : int] EQINT () = eqint_make_gint j |
|||
val pivot_val = A[j, j] |
|||
var p : intGte 1 |
|||
in |
|||
needlessly_set_to_value (A, j, j, One); |
|||
for* {p : int | j + 1 <= p; p <= n + 1} .<(n + 1) - p>. |
|||
(p : int p) => |
|||
(p := succ j; p <> succ n; p := succ p) |
|||
A[j, p] := A[j, p] / pivot_val; |
|||
for* {p : int | 1 <= p; p <= k + 1} .<(k + 1) - p>. |
|||
(p : int p) => |
|||
(p := 1; p <> succ k; p := succ p) |
|||
B[j, p] := B[j, p] / pivot_val; |
|||
end |
|||
fn |
|||
subtract_normalized_pivot_row |
|||
(i : one_to_n, j : one_to_n) :<!refwrt> void = |
|||
let |
|||
prval [j : int] EQINT () = eqint_make_gint j |
|||
val lead_val = A[i, j] |
|||
in |
|||
if lead_val <> Zero then |
|||
let |
|||
val factor = ~lead_val |
|||
var p : intGte 1 |
|||
in |
|||
needlessly_set_to_value (A, i, j, Zero); |
|||
for* {p : int | j + 1 <= p; p <= n + 1} .<(n + 1) - p>. |
|||
(p : int p) => |
|||
(p := succ j; p <> succ n; p := succ p) |
|||
A[i, p] := |
|||
multiply_and_add (A[j, p], factor, A[i, p]); |
|||
for* {p : int | 1 <= p; p <= k + 1} .<(k + 1) - p>. |
|||
(p : int p) => |
|||
(p := 1; p <> succ k; p := succ p) |
|||
B[i, p] := |
|||
multiply_and_add (B[j, p], factor, B[i, p]) |
|||
end |
|||
end |
|||
fun |
|||
main_loop {j : pos | j <= n + 1} .<(n + 1) - j>. |
|||
(j : int j, |
|||
success : &bool? >> bool) :<!refwrt> void = |
|||
if j = succ n then |
|||
success := true |
|||
else |
|||
let |
|||
fun |
|||
select_pivot {i : int | j <= i; i <= n + 1} |
|||
.<(n + 1) - i>. |
|||
(i : int i, |
|||
max_abs : g0float tk, |
|||
i_max_abs : intBtwe (j - 1, n)) |
|||
:<!ref> intBtwe (j - 1, n) = |
|||
if i = succ n then |
|||
i_max_abs |
|||
else |
|||
let |
|||
val abs_aij = abs A[i, j] |
|||
in |
|||
if abs_aij > max_abs then |
|||
select_pivot (succ i, abs_aij, i) |
|||
else |
|||
select_pivot (succ i, max_abs, i_max_abs) |
|||
end |
|||
val i_pivot = select_pivot (j, Zero, pred j) |
|||
prval [i_pivot : int] EQINT () = eqint_make_gint i_pivot |
|||
in |
|||
if i_pivot = pred j then |
|||
success := false |
|||
else |
|||
let |
|||
var i : intGte 1 |
|||
in |
|||
exchange_rows (i_pivot, j); |
|||
normalize_pivot_row (j); |
|||
(* For Gauss-Jordan elimination, we would do this, |
|||
instead of switching to back substitution: |
|||
for* {i : int | 1 <= i; i <= j} |
|||
.<j - i>. |
|||
(i : int i) => |
|||
(i := 1; i <> j; i := succ i) |
|||
subtract_normalized_pivot_row (i, j); *) |
|||
for* {i : int | j + 1 <= i; i <= n + 1} |
|||
.<(n + 1) - i>. |
|||
(i : int i) => |
|||
(i := succ j; i <> succ n; i := succ i) |
|||
subtract_normalized_pivot_row (i, j); |
|||
main_loop (succ j, success) |
|||
end |
|||
end |
|||
var success : bool |
|||
val () = main_loop (1, success) |
|||
in |
|||
if ~success then |
|||
None () |
|||
else |
|||
(* Back substitution. (Doing this with block operations on rows, |
|||
as is done below, is not the most efficient way, but helps |
|||
convey what is going on.) *) |
|||
let |
|||
val bottom_row = block (B, n, n, 1, k) |
|||
(* The rows array will treat the rows of B as |
|||
submatrices. (The zeroth entry will not be used.) *) |
|||
val rows = arrayref_make_elt (i2sz (succ n), bottom_row) |
|||
var i : intGte 0 |
|||
in |
|||
(* Fill in the rows array (ignoring its zeroth entry). *) |
|||
for* {i : nat | i <= n - 1} .<i>. |
|||
(i : int i) => |
|||
(i := pred n; i <> 0; i := pred i) |
|||
rows[i] := block (B, i, i, 1, k); |
|||
(* Now do back substitution, one solution row at a time. *) |
|||
for* {i : nat | i <= n - 1} .<i>. |
|||
(i : int i) => |
|||
(i := pred n; i <> 0; i := pred i) |
|||
let |
|||
var j : intGte 0 |
|||
in |
|||
for* {j : int | i <= j; j <= n} .<j>. |
|||
(j : int j) => |
|||
(j := n; j <> i; j := pred j) |
|||
scalar_multiply_and_add_to |
|||
(rows[i], rows[j], ~A[i, j], rows[i]) |
|||
end; |
|||
(* The returned matrix will "contain" the rows_permutation |
|||
array and some extra closures. If you want a "clean" |
|||
matrix, you can use Real_Matrix_copy to get one. *) |
|||
Some B |
|||
end |
|||
end |
|||
overload gaussian_elimination with Real_Matrix_gaussian_elimination |
|||
(*------------------------------------------------------------------*) |
|||
fn {tk : tkind} |
|||
fprint_matrices_and_solutions |
|||
{n, k : pos} |
|||
(outf : FILEref, |
|||
A : Real_Matrix (tk, n, n), |
|||
B : Real_Matrix (tk, n, k)) : void = |
|||
let |
|||
typedef FILEstar = $extype"FILE *" |
|||
extern castfn FILEref2star : FILEref -<> FILEstar |
|||
macdef fmt = "%9.4f" |
|||
macdef left_bracket = "[" |
|||
macdef right_bracket = " ]" |
|||
macdef product = " ✕ " |
|||
macdef equals = " = " |
|||
macdef spacing = " " |
|||
macdef msg_for_singular = " appears to be singular" |
|||
macdef print_num (x) = |
|||
ignoret ($extfcall (int, "fprintf", FILEref2star outf, |
|||
fmt, ,(x))) |
|||
val @(n, k) = dimension B |
|||
in |
|||
case+ gaussian_elimination<tk> (A, B) of |
|||
| None () => |
|||
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) |
|||
let |
|||
var j : intGte 1 |
|||
in |
|||
fprint! (outf, left_bracket); |
|||
for* {j : pos | j <= n + 1} .<(n + 1) - j>. |
|||
(j : int j) => |
|||
(j := 1; j <> succ n; j := succ j) |
|||
print_num A[i, j]; |
|||
fprint! (outf, right_bracket); |
|||
if pred i = half n then |
|||
fprint! (outf, msg_for_singular); |
|||
fprintln! (outf) |
|||
end |
|||
end |
|||
| Some X => |
|||
let |
|||
val AX = A * X |
|||
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) |
|||
let |
|||
var j : intGte 1 |
|||
in |
|||
fprint! (outf, left_bracket); |
|||
for* {j : pos | j <= n + 1} .<(n + 1) - j>. |
|||
(j : int j) => |
|||
(j := 1; j <> succ n; j := succ j) |
|||
print_num A[i, j]; |
|||
fprint! (outf, right_bracket); |
|||
if pred i = half n then |
|||
fprint! (outf, product) |
|||
else |
|||
fprint! (outf, spacing); |
|||
fprint! (outf, left_bracket); |
|||
for* {j : pos | j <= k + 1} .<(k + 1) - j>. |
|||
(j : int j) => |
|||
(j := 1; j <> succ k; j := succ j) |
|||
print_num X[i, j]; |
|||
fprint! (outf, right_bracket); |
|||
if pred i = half n then |
|||
fprint! (outf, equals) |
|||
else |
|||
fprint! (outf, spacing); |
|||
fprint! (outf, left_bracket); |
|||
for* {j : pos | j <= k + 1} .<(k + 1) - j>. |
|||
(j : int j) => |
|||
(j := 1; j <> succ k; j := succ j) |
|||
print_num AX[i, j]; |
|||
fprint! (outf, right_bracket); |
|||
fprintln! (outf) |
|||
end |
|||
end |
|||
end |
|||
fn {tk : tkind} |
|||
column_major_list_to_matrix |
|||
{m, n : pos} |
|||
(m : int m, |
|||
n : int n, |
|||
lst : list (g0float tk, m * n)) |
|||
: Real_Matrix (tk, m, n) = |
|||
let |
|||
#define :: list_cons |
|||
prval () = mul_gte_gte_gte {m, n} () |
|||
val A = Real_Matrix_make_elt (m, 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 <= m + 1} .<(m + 1) - i>. |
|||
(i : int i) => |
|||
(i := 1; i <> succ m; i := succ i) |
|||
case- !lstref of |
|||
| hd :: tl => |
|||
begin |
|||
A[i, j] := hd; |
|||
!lstref := tl |
|||
end |
|||
end; |
|||
A |
|||
end |
|||
fn {tk : tkind} |
|||
row_major_list_to_matrix |
|||
{m, n : pos} |
|||
(m : int m, |
|||
n : int n, |
|||
lst : list (g0float tk, m * n)) |
|||
: Real_Matrix (tk, m, n) = |
|||
transpose (column_major_list_to_matrix<tk> (n, m, lst)) |
|||
macdef separator = "\n" |
|||
fn |
|||
print_example |
|||
{n, k : pos} |
|||
(n : int n, |
|||
k : int k, |
|||
lst_A : list (double, n * n), |
|||
lst_B : list (double, n * k)) : void = |
|||
let |
|||
val A = row_major_list_to_matrix (n, n, lst_A) |
|||
and B = row_major_list_to_matrix (n, k, lst_B) |
|||
in |
|||
print! separator; |
|||
fprint_matrices_and_solutions<dblknd> (stdout_ref, A, B) |
|||
end |
|||
implement |
|||
main0 () = |
|||
begin |
|||
println! ("\n(The examples are printed here after ", |
|||
"rounding to 4 decimal places.)"); |
|||
println! ("\nSolving Ax=b where ", |
|||
"transpose(b)=[-0.01 0.61 0.91 0.99 0.60 0.02]:"); |
|||
print_example |
|||
(6, 1, $list (1.00, 0.00, 0.00, 0.00, 0.00, 0.00, |
|||
1.00, 0.63, 0.39, 0.25, 0.16, 0.10, |
|||
1.00, 1.26, 1.58, 1.98, 2.49, 3.13, |
|||
1.00, 1.88, 3.55, 6.70, 12.62, 23.80, |
|||
1.00, 2.51, 6.32, 15.88, 39.90, 100.28, |
|||
1.00, 3.14, 9.87, 31.01, 97.41, 306.02), |
|||
$list (~0.01, 0.61, 0.91, 0.99, 0.60, 0.02)); |
|||
println! ("\nAn orthonormal matrix of rank 4,", |
|||
" solving for its inverse:"); |
|||
print_example |
|||
(4, 4, |
|||
$list (~0.1543033499620918, ~0.1307837816808, |
|||
0.8649377296669811, 0.4593134032861146, |
|||
~0.6172133998483675, ~0.2062359634197235, |
|||
0.2410548538544755, ~0.7200047943403952, |
|||
~0.7715167498104593, 0.1911455270719389, |
|||
~0.3658314290169772, 0.4841411548150931, |
|||
0.0, ~0.950697489910433, |
|||
~0.2448318745080428, 0.190346095055507), |
|||
$list (1.0, 0.0, 0.0, 0.0, |
|||
0.0, 1.0, 0.0, 0.0, |
|||
0.0, 0.0, 1.0, 0.0, |
|||
0.0, 0.0, 0.0, 1.0)); |
|||
print! separator |
|||
end |
|||
(*------------------------------------------------------------------*) |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
On my computer, the following compiler options result in "fused multiply-and-add" instructions being generated. See the program text for a brief discussion of how that optimization may affect results, when the '''A''' matrix that is singular or nearly so. |
|||
<pre style="font-size:90%">$ patscc -std=gnu2x -g -O3 -march=native -DATS_MEMALLOC_GCBDW gaussian_elimination_task.dats -lgc && ./a.out |
|||
(The examples are printed here after rounding to 4 decimal places.) |
|||
Solving Ax=b where transpose(b)=[-0.01 0.61 0.91 0.99 0.60 0.02]: |
|||
[ 1.0000 0.0000 0.0000 0.0000 0.0000 0.0000 ] [ -0.0100 ] [ -0.0100 ] |
|||
[ 1.0000 0.6300 0.3900 0.2500 0.1600 0.1000 ] [ 1.6028 ] [ 0.6100 ] |
|||
[ 1.0000 1.2600 1.5800 1.9800 2.4900 3.1300 ] [ -1.6132 ] [ 0.9100 ] |
|||
[ 1.0000 1.8800 3.5500 6.7000 12.6200 23.8000 ] ✕ [ 1.2455 ] = [ 0.9900 ] |
|||
[ 1.0000 2.5100 6.3200 15.8800 39.9000 100.2800 ] [ -0.4910 ] [ 0.6000 ] |
|||
[ 1.0000 3.1400 9.8700 31.0100 97.4100 306.0200 ] [ 0.0658 ] [ 0.0200 ] |
|||
An orthonormal matrix of rank 4, solving for its inverse: |
|||
[ -0.1543 -0.1308 0.8649 0.4593 ] [ -0.1543 -0.6172 -0.7715 -0.0000 ] [ 1.0000 0.0000 -0.0000 0.0000 ] |
|||
[ -0.6172 -0.2062 0.2411 -0.7200 ] [ -0.1308 -0.2062 0.1911 -0.9507 ] [ -0.0000 1.0000 -0.0000 0.0000 ] |
|||
[ -0.7715 0.1911 -0.3658 0.4841 ] ✕ [ 0.8649 0.2411 -0.3658 -0.2448 ] = [ -0.0000 -0.0000 1.0000 0.0000 ] |
|||
[ 0.0000 -0.9507 -0.2448 0.1903 ] [ 0.4593 -0.7200 0.4841 0.1903 ] [ -0.0000 0.0000 0.0000 1.0000 ] |
|||
</pre> |
</pre> |
||