69 lines
2.5 KiB
D
69 lines
2.5 KiB
D
import std.stdio, std.math, std.algorithm, std.range, std.numeric,
|
|
std.typecons;
|
|
|
|
Tuple!(double[],"x", string,"err")
|
|
gaussPartial(in double[][] a0, in double[] b0) pure /*nothrow*/
|
|
in {
|
|
assert(a0.length == a0[0].length);
|
|
assert(a0.length == b0.length);
|
|
assert(a0.all!(row => row.length == a0[0].length));
|
|
} body {
|
|
enum eps = 1e-6;
|
|
immutable m = b0.length;
|
|
|
|
// Make augmented matrix.
|
|
//auto a = a0.zip(b0).map!(c => c[0] ~ c[1]).array; // Not mutable.
|
|
auto a = a0.zip(b0).map!(c => [] ~ c[0] ~ c[1]).array;
|
|
|
|
// Wikipedia algorithm from Gaussian elimination page,
|
|
// produces row-eschelon form.
|
|
foreach (immutable k; 0 .. a.length) {
|
|
// Find pivot for column k and swap.
|
|
a[k .. m].minPos!((x, y) => x[k] > y[k]).front.swap(a[k]);
|
|
if (a[k][k].abs < eps)
|
|
return typeof(return)(null, "singular");
|
|
|
|
// Do for all rows below pivot.
|
|
foreach (immutable i; k + 1 .. m) {
|
|
// Do for all remaining elements in current row.
|
|
a[i][k+1 .. m+1] -= a[k][k+1 .. m+1] * (a[i][k] / a[k][k]);
|
|
|
|
a[i][k] = 0; // Fill lower triangular matrix with zeros.
|
|
}
|
|
}
|
|
|
|
// End of WP algorithm. Now back substitute to get result.
|
|
auto x = new double[m];
|
|
foreach_reverse (immutable i; 0 .. m)
|
|
x[i] = (a[i][m] - a[i][i+1 .. m].dotProduct(x[i+1 .. m])) / a[i][i];
|
|
|
|
return typeof(return)(x, null);
|
|
}
|
|
|
|
void main() {
|
|
// The test case result is correct to this tolerance.
|
|
enum eps = 1e-13;
|
|
|
|
// Common RC example. Result computed with rational arithmetic
|
|
// then converted to double, and so should be about as close to
|
|
// correct as double represention allows.
|
|
immutable 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]];
|
|
immutable b = [-0.01, 0.61, 0.91, 0.99, 0.60, 0.02];
|
|
|
|
immutable r = gaussPartial(a, b);
|
|
if (!r.err.empty)
|
|
return writeln("Error: ", r.err);
|
|
r.x.writeln;
|
|
|
|
immutable result = [-0.01, 1.602790394502114,
|
|
-1.6132030599055613, 1.2454941213714368,
|
|
-0.4909897195846576, 0.065760696175232];
|
|
foreach (immutable i, immutable xi; result)
|
|
if (abs(r.x[i] - xi) > eps)
|
|
return writeln("Out of tolerance: ", r.x[i], " ", xi);
|
|
}
|