Gauss-Jordan matrix inversion: Difference between revisions
Content added Content deleted
Line 519: | Line 519: | ||
Gauss-Jordan algorithm. *) |
Gauss-Jordan algorithm. *) |
||
#define DO_THINGS_THAT_DO_NOT_NEED_TO_BE_DONE 0 |
#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: | Line 540: | ||
macdef One = g0i2f 1 |
macdef One = g0i2f 1 |
||
macdef Two = g0i2f 2 |
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: | Line 1,000: | ||
(j := 2; j <> succ n; j := succ j) |
(j := 2; j <> succ n; j := succ j) |
||
C[i, k] := |
C[i, k] := |
||
multiply_and_add (A[i, j], B[j, k], C[i, k]) |
|||
end |
end |
||
end |
end |
||
Line 1,133: | Line 1,155: | ||
if lead_val <> Zero then |
if lead_val <> Zero then |
||
let |
let |
||
val factor = ~lead_val |
|||
var k : intGte 1 |
var k : intGte 1 |
||
in |
in |
||
Line 1,139: | Line 1,162: | ||
(k : int k) => |
(k : int k) => |
||
(k := succ j; k <> succ n; k := succ k) |
(k := succ j; k <> succ n; k := succ k) |
||
A[i, 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>. |
for* {k : int | 1 <= k; k <= n + 1} .<(n + 1) - k>. |
||
(k : int k) => |
(k : int k) => |
||
(k := 1; k <> succ n; k := succ k) |
(k := 1; k <> succ n; k := succ k) |
||
B[i, k] := |
B[i, k] := |
||
multiply_and_add (B[j, k], factor, B[i, k]) |
|||
end |
end |
||
end |
end |
||
Line 1,236: | Line 1,261: | ||
val @(n, _) = dimension A |
val @(n, _) = dimension A |
||
in |
in |
||
case+ inverse_by_gauss_jordan A of |
case+ inverse_by_gauss_jordan<tk> (A) of |
||
| None () => |
| None () => |
||
let |
let |
||
Line 1,375: | Line 1,400: | ||
~0.2461829819586655)); |
~0.2461829819586655)); |
||
(* The following matrix may or may not be detected as singular, |
|||
⚫ | |||
depending on how rounding is done. *) |
|||
(* |
|||
print_example |
print_example |
||
(5, $list (~0.08006407690254358, ~0.1884825001438613, |
(5, $list (~0.08006407690254358, ~0.1884825001438613, |
||
~0.8954592676839166, 0.2715238629598956, |
~0.8954592676839166, 0.2715238629598956, |
||
~0.2872134789517761, ~0. |
~0.2872134789517761, ~0.5604485383178048, |
||
0. |
~0.4587939881550582, ~0.01946650581921516, |
||
0.1609525690123523, 0.6701647842208114, |
|||
~0. |
~0.2401922307076307, 0.1517054269450593, |
||
0. |
~0.3178281430868086, ~0.8979459113320678, |
||
0.1094146586482966, ~0.3202563076101743, |
|||
0. |
~0.6104994151001173, 0.2983616372675922, |
||
0. |
~0.1997416922923375, ~0.6291342872277008, |
||
5.0 * ~0.08006407690254358, |
5.0 * ~0.08006407690254358, |
||
5.0 * ~0.1884825001438613, |
5.0 * ~0.1884825001438613, |
||
5.0 * ~0.8954592676839166, |
5.0 * ~0.8954592676839166, |
||
5.0 * 0.2715238629598956, |
5.0 * 0.2715238629598956, |
||
5.0 * ~0.2872134789517761)); |
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 |
print! separator |
||
Line 1,401: | Line 1,437: | ||
{{out}} |
{{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.) |
|||
(When compiling specifically for Opteron 6344 with GCC 12 (and nearby versions), I had trouble unless I added <code>-fno-expensive-optimizations</code>. It is some problem involving fused multiply-add instructions. The same problem occurs if you use the built-in <code>fma</code> function of GCC. Other projects have had similar difficulties.) |
|||
<pre style="font-size:90%">$ patscc -std=gnu2x -g - |
<pre style="font-size:90%">$ patscc -std=gnu2x -g -O3 -march=native -DATS_MEMALLOC_GCBDW gauss_jordan_task.dats -lgc -lm && ./a.out |
||
(The examples are printed here after rounding to 4 decimal places.) |
(The examples are printed here after rounding to 4 decimal places.) |
||
Line 1,425: | Line 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.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.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.2715 0.1610 -0.8979 -0.1997 0.2327 ] [ -0.3203 -0.6105 0.2984 -0.1997 -0.6291 ] [ -0.0000 0.0000 -0.0000 1.0000 -0.0000 ] |
||
[ -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 ] |
[ -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 |
A 5✕5 singular matrix: |
||
[ |
[ 1.0000 5.0000 3.0000 1.0000 4.0000 ] |
||
[ |
[ 1.0000 1.0000 2.0000 1.0000 1.0000 ] |
||
[ |
[ 2.0000 2.0000 3.0000 2.0000 2.0000 ] appears to be singular |
||
[ |
[ 2.0000 2.0000 5.0000 2.0000 2.0000 ] |
||
[ |
[ 2.0000 2.0000 7.0000 2.0000 2.0000 ] |
||
</pre> |
</pre> |