Gaussian elimination: Difference between revisions
Content deleted Content added
→version 3: added REXX version 3, changed indentation of the Gauss_elim subroutine. |
→{{header|Go}}: add version with scaling |
||
Line 680: | Line 680: | ||
=={{header|Go}}== |
=={{header|Go}}== |
||
===Partial pivoting, no scaling=== |
|||
Gaussian elimination with partial pivoting by pseudocode on WP page "Gaussian elimination." Scaling not implemented. |
|||
Gaussian elimination with partial pivoting by [https://en.wikipedia.org/wiki/Gaussian_elimination#Pseudocode pseudocode] on WP page [https://en.wikipedia.org/wiki/Gaussian_elimination Gaussian elimination]." |
|||
<lang go>package main |
<lang go>package main |
||
Line 781: | Line 782: | ||
<pre> |
<pre> |
||
[-0.01 1.6027903945020987 -1.613203059905494 1.245494121371364 -0.49098971958462834 0.06576069617522803] |
[-0.01 1.6027903945020987 -1.613203059905494 1.245494121371364 -0.49098971958462834 0.06576069617522803] |
||
</pre> |
|||
===Scaled partial pivoting=== |
|||
Changes from above version noted with comments. For the example data scaling does help a bit. |
|||
<lang go>package main |
|||
import ( |
|||
"errors" |
|||
"fmt" |
|||
"log" |
|||
"math" |
|||
) |
|||
type testCase struct { |
|||
a [][]float64 |
|||
b []float64 |
|||
x []float64 |
|||
} |
|||
var tc = testCase{ |
|||
a: [][]float64{ |
|||
{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: []float64{-0.01, 0.61, 0.91, 0.99, 0.60, 0.02}, |
|||
x: []float64{-0.01, 1.602790394502114, -1.6132030599055613, |
|||
1.2454941213714368, -0.4909897195846576, 0.065760696175232}, |
|||
} |
|||
// result from above test case turns out to be correct to this tolerance. |
|||
const ε = 1e-14 |
|||
func main() { |
|||
x, err := GaussPartial(tc.a, tc.b) |
|||
if err != nil { |
|||
log.Fatal(err) |
|||
} |
|||
fmt.Println(x) |
|||
for i, xi := range x { |
|||
if math.Abs(tc.x[i]-xi) > ε { |
|||
log.Println("out of tolerance") |
|||
log.Fatal("expected", tc.x) |
|||
} |
|||
} |
|||
} |
|||
func GaussPartial(a0 [][]float64, b0 []float64) ([]float64, error) { |
|||
m := len(b0) |
|||
a := make([][]float64, m) |
|||
for i, ai := range a0 { |
|||
row := make([]float64, m+1) |
|||
copy(row, ai) |
|||
row[m] = b0[i] |
|||
a[i] = row |
|||
} |
|||
for k := range a { |
|||
iMax := 0 |
|||
max := -1. |
|||
for i := k; i < m; i++ { |
|||
row := a[i] |
|||
// compute scale factor s = max abs in row |
|||
s := -1. |
|||
for j := k; j < m; j++ { |
|||
x := math.Abs(row[j]) |
|||
if x > s { |
|||
s = x |
|||
} |
|||
} |
|||
// scale the abs used to pick the pivot. |
|||
if abs := math.Abs(row[k]) / s; abs > max { |
|||
iMax = i |
|||
max = abs |
|||
} |
|||
} |
|||
if a[iMax][k] == 0 { |
|||
return nil, errors.New("singular") |
|||
} |
|||
a[k], a[iMax] = a[iMax], a[k] |
|||
for i := k + 1; i < m; i++ { |
|||
for j := k + 1; j <= m; j++ { |
|||
a[i][j] -= a[k][j] * (a[i][k] / a[k][k]) |
|||
} |
|||
a[i][k] = 0 |
|||
} |
|||
} |
|||
x := make([]float64, m) |
|||
for i := m - 1; i >= 0; i-- { |
|||
x[i] = a[i][m] |
|||
for j := i + 1; j < m; j++ { |
|||
x[i] -= a[i][j] * x[j] |
|||
} |
|||
x[i] /= a[i][i] |
|||
} |
|||
return x, nil |
|||
}</lang> |
|||
{{out}} |
|||
<pre> |
|||
[-0.01 1.6027903945021131 -1.6132030599055596 1.245494121371436 -0.49098971958465754 0.065760696175232] |
|||
</pre> |
</pre> |
||