Gaussian elimination: Difference between revisions
add output of the java code, add output of the scala code
Thundergnat (talk | contribs) m (syntax highlighting fixup automation) |
(add output of the java code, add output of the scala code) |
||
| (5 intermediate revisions by 4 users not shown) | |||
Line 537:
<pre>
( -.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>
Line 1,229 ⟶ 2,271:
{{out}}
<pre>1.00, 0.00, 0.00, 4.00</pre>
=={{header|EasyLang}}==
<syntaxhighlight>
proc gauss_elim . a[][] b[] x[] .
n = len a[][]
for i to n
maxr = i
maxv = abs a[i][i]
for j = i + 1 to n
if abs a[j][i] > maxv
maxr = j
maxv = abs a[j][i]
.
.
if maxr <> i
swap a[maxr][] a[i][]
swap b[maxr] b[i]
.
for j = i + 1 to n
f = a[j][i] / a[i][i]
for k = i to n
a[j][k] -= f * a[i][k]
.
b[j] -= f * b[i]
.
.
x[] = [ ]
len x[] n
for i = n downto 1
rhs = b[i]
for j = i + 1 to n
rhs -= a[i][j] * x[j]
.
x[i] = rhs / a[i][i]
.
.
a[][] = [ [ 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 ] ]
b[] = [ -0.01 0.61 0.91 0.99 0.60 0.02 ]
gauss_elim a[][] b[] x[]
print x[]
</syntaxhighlight>
=={{header|F_Sharp|F#}}==
Line 2,841 ⟶ 3,924:
}
}</syntaxhighlight>
{{out}}
<pre>
det: 780.0
0.09750000 0.0000e+00
0.11000000 0.0000e+00
0.12916667 0.0000e+00
0.12333333 1.3878e-17
0.17750000 2.7756e-17
</pre>
=={{header|JavaScript}}==
Line 5,405 ⟶ 6,498:
}
</syntaxhighlight>
=={{header|Scala}}==
{{trans|Java}}
<syntaxhighlight lang="Scala">
object GaussianElimination {
def solve(a: Array[Array[Double]], b: Array[Array[Double]]): Double = {
if (a == null || b == null || a.length == 0 || b.length == 0) {
throw new IllegalArgumentException("Invalid dimensions")
}
val n = b.length
val p = b(0).length
if (a.length != n || a(0).length != n) {
throw new IllegalArgumentException("Invalid dimensions")
}
var det = 1.0
for (i <- 0 until n - 1) {
var k = i
for (j <- i + 1 until n) {
if (Math.abs(a(j)(i)) > Math.abs(a(k)(i))) {
k = j
}
}
if (k != i) {
det = -det
for (j <- i until n) {
val s = a(i)(j)
a(i)(j) = a(k)(j)
a(k)(j) = s
}
for (j <- 0 until p) {
val s = b(i)(j)
b(i)(j) = b(k)(j)
b(k)(j) = s
}
}
for (j <- i + 1 until n) {
val s = a(j)(i) / a(i)(i)
for (k <- i + 1 until n) {
a(j)(k) -= s * a(i)(k)
}
for (k <- 0 until p) {
b(j)(k) -= s * b(i)(k)
}
}
}
for (i <- n - 1 to 0 by -1) {
for (j <- i + 1 until n) {
val s = a(i)(j)
for (k <- 0 until p) {
b(i)(k) -= s * b(j)(k)
}
}
val s = a(i)(i)
det *= s
for (k <- 0 until p) {
b(i)(k) /= s
}
}
det
}
def main(args: Array[String]): Unit = {
val a = Array(
Array(4.0, 1.0, 0.0, 0.0, 0.0),
Array(1.0, 4.0, 1.0, 0.0, 0.0),
Array(0.0, 1.0, 4.0, 1.0, 0.0),
Array(0.0, 0.0, 1.0, 4.0, 1.0),
Array(0.0, 0.0, 0.0, 1.0, 4.0)
)
val b = Array(
Array(1.0 / 2.0),
Array(2.0 / 3.0),
Array(3.0 / 4.0),
Array(4.0 / 5.0),
Array(5.0 / 6.0)
)
val x = Array(39.0 / 400.0, 11.0 / 100.0, 31.0 / 240.0, 37.0 / 300.0, 71.0 / 400.0)
println("det: " + solve(a, b))
for (i <- 0 until 5) {
printf("%12.8f %12.4e\n", b(i)(0), b(i)(0) - x(i))
}
}
}
</syntaxhighlight>
{{out}}
<pre>
det: 780.0
0.09750000 0.0000e+00
0.11000000 0.0000e+00
0.12916667 0.0000e+00
0.12333333 1.3878e-17
0.17750000 2.7756e-17
</pre>
=={{header|Sidef}}==
Line 5,782 ⟶ 6,988:
=={{header|Wren}}==
{{trans|Kotlin}}
{{libheader|Wren-
<syntaxhighlight lang="
var ta = [
| |||