1010 lines
30 KiB
Text
1010 lines
30 KiB
Text
(* 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
|
|
|
|
(*------------------------------------------------------------------*)
|