Gauss-Jordan matrix inversion: Difference between revisions
→{{header|ATS}}
Line 519:
Gauss-Jordan 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
%{^
Line 531 ⟶ 540:
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
(*------------------------------------------------------------------*)
Line 978 ⟶ 1,000:
(j := 2; j <> succ n; j := succ j)
C[i, k] :=
end
end
Line 1,133 ⟶ 1,155:
if lead_val <> Zero then
let
val factor = ~lead_val
var k : intGte 1
in
Line 1,139 ⟶ 1,162:
(k : int k) =>
(k := succ j; k <> succ n; k := succ k)
A[i, k] :=
multiply_and_add (A[j, k], factor, A[i, k]);
for* {k : int | 1 <= k; k <= n + 1} .<(n + 1) - k>.
(k : int k) =>
(k := 1; k <> succ n; k := succ k)
B[i, k] :=
multiply_and_add (B[j, k], factor, B[i, k])
end
end
Line 1,236 ⟶ 1,261:
val @(n, _) = dimension A
in
case+ inverse_by_gauss_jordan<tk> (A) of
| None () =>
let
Line 1,375 ⟶ 1,400:
~0.2461829819586655));
(* The following matrix may or may not be detected as singular,
println! ("\nA 5✕5 singular matrix (rank 4):");▼
depending on how rounding is done. *)
(*
print_example
(5, $list (~0.08006407690254358, ~0.1884825001438613,
~0.8954592676839166, 0.2715238629598956,
~0.2872134789517761, ~0.
~0.
~0.
~0.
~0.
~0.
5.0 * ~0.08006407690254358,
5.0 * ~0.1884825001438613,
5.0 * ~0.8954592676839166,
5.0 * 0.2715238629598956,
5.0 * ~0.2872134789517761));
*)
print_example
(5, $list (1.0, 1.0, 2.0, 2.0, 2.0,
5.0, 1.0, 2.0, 2.0, 2.0,
3.0, 2.0, 3.0, 5.0, 7.0,
1.0, 1.0, 2.0, 2.0, 2.0,
4.0, 1.0, 2.0, 2.0, 2.0));
print! separator
Line 1,401 ⟶ 1,437:
{{out}}
(The following settings on my computer ''will'' result in GCC generating fused-multiply-and-add instructions [and thus the ''-lm'' flag actually is not needed]. See comments in the program for a discussion of the matter.)
<pre style="font-size:90%">$ patscc -std=gnu2x -g -
(The examples are printed here after rounding to 4 decimal places.)
Line 1,425 ⟶ 1,461:
[ -0.1885 -0.4588 0.1517 -0.6105 0.5985 ] [ -0.5604 -0.4588 -0.0195 0.1610 0.6702 ] [ 0.0000 1.0000 -0.0000 0.0000 0.0000 ]
[ -0.8955 -0.0195 -0.3178 0.2984 0.0880 ] ✕ [ -0.2402 0.1517 -0.3178 -0.8979 0.1094 ] = [ 0.0000 -0.0000 1.0000 0.0000 0.0000 ]
[ 0.2715 0.1610 -0.8979 -0.1997 0.2327 ] [ -0.3203 -0.6105 0.2984 -0.1997 -0.6291 ] [
[ -0.2872 0.6702 0.1094 -0.6291 -0.2462 ] [ -0.7206 0.5985 0.0880 0.2327 -0.2462 ] [ -0.0000 0.0000 -0.0000 0.0000 1.0000 ]
A 5✕5 singular matrix
[
[
[
[
[
</pre>
| |||