Gauss-Jordan matrix inversion: Difference between revisions
Put Generic at the right place
(→{{header|Go}}: Trans PowerShell) |
(Put Generic at the right place) |
||
Line 480:
</pre>
=={{header|Generic}}==
<lang cpp>
Line 1,416 ⟶ 1,154:
}
</lang>
=={{header|Go}}==
{{trans|Kotlin}}
<lang go>package main
import "fmt"
type vector = []float64
type matrix []vector
func (m matrix) inverse() matrix {
le := len(m)
for _, v := range m {
if len(v) != le {
panic("Not a square matrix")
}
}
aug := make(matrix, le)
for i := 0; i < le; i++ {
aug[i] = make(vector, 2*le)
copy(aug[i], m[i])
// augment by identity matrix to right
aug[i][i+le] = 1
}
aug.toReducedRowEchelonForm()
inv := make(matrix, le)
// remove identity matrix to left
for i := 0; i < le; i++ {
inv[i] = make(vector, le)
copy(inv[i], aug[i][le:])
}
return inv
}
// note: this mutates the matrix in place
func (m matrix) toReducedRowEchelonForm() {
lead := 0
rowCount, colCount := len(m), len(m[0])
for r := 0; r < rowCount; r++ {
if colCount <= lead {
return
}
i := r
for m[i][lead] == 0 {
i++
if rowCount == i {
i = r
lead++
if colCount == lead {
return
}
}
}
m[i], m[r] = m[r], m[i]
if div := m[r][lead]; div != 0 {
for j := 0; j < colCount; j++ {
m[r][j] /= div
}
}
for k := 0; k < rowCount; k++ {
if k != r {
mult := m[k][lead]
for j := 0; j < colCount; j++ {
m[k][j] -= m[r][j] * mult
}
}
}
lead++
}
}
func (m matrix) print(title string) {
fmt.Println(title)
for _, v := range m {
fmt.Printf("% f\n", v)
}
fmt.Println()
}
func main() {
a := matrix{{1, 2, 3}, {4, 1, 6}, {7, 8, 9}}
a.inverse().print("Inverse of A is:\n")
b := matrix{{2, -1, 0}, {-1, 2, -1}, {0, -1, 2}}
b.inverse().print("Inverse of B is:\n")
}</lang>
{{out}}
<pre>
Inverse of A is:
[-0.812500 0.125000 0.187500]
[ 0.125000 -0.250000 0.125000]
[ 0.520833 0.125000 -0.145833]
Inverse of B is:
[ 0.750000 0.500000 0.250000]
[ 0.500000 1.000000 0.500000]
[ 0.250000 0.500000 0.750000]
</pre>
===Alternative===
{{trans|PoweShell}}
<lang go>package main
import (
"fmt"
"math"
)
type vector = []float64
type matrix []vector
func (m matrix) print(title string) {
fmt.Println(title)
for _, v := range m {
fmt.Printf("% f\n", v)
}
}
func I(n int) matrix {
b := make(matrix, n)
for i := 0; i < n; i++ {
a := make(vector, n)
a[i] = 1
b[i] = a
}
return b
}
func (m matrix) inverse() matrix {
n := len(m)
for _, v := range m {
if (n != len(v)){
panic("Not a square matrix")
}
}
b := I(n)
a := make(matrix, n)
for i, v := range m {
row := make([]float64, n)
copy(row, v)
a[i] = row
}
for k := range a {
iMax := 0
max := -1.
for i := k; i < n; i++ {
row := a[i]
// compute scale factor s = max abs in row
s := -1.
for j := k; j < n; 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 {
panic("Irregular matrix")
}
if k != iMax {
a[k], a[iMax] = a[iMax], a[k]
b[k], b[iMax] = b[iMax], b[k]
}
akk := a[k][k]
for j := 0; j < n; j++ {
a[k][j] /= akk
b[k][j] /= akk
}
for i := 0; i < n; i++ {
if (i != k) {
aik := a[i][k]
for j := 0; j < n; j++ {
a[i][j] -= a[k][j] * aik
b[i][j] -= b[k][j] * aik
}
}
}
}
return b
}
func main() {
a := matrix{{1, 2, 3}, {4, 1, 6}, {7, 8, 9}}
a.inverse().print("Inverse of A is:\n")
b := matrix{{2, -1, 0}, {-1, 2, -1}, {0, -1, 2}}
b.inverse().print("Inverse of B is:\n")
}</lang>
{{out}}
<pre>Same output than above</pre>
=={{header|Haskell}}==
<lang Haskell>isMatrix xs = null xs || all ((== (length.head $ xs)).length) xs
isSquareMatrix xs = null xs || all ((== (length xs)).length) xs
mult:: Num a => [[a]] -> [[a]] -> [[a]]
mult uss vss = map ((\xs -> if null xs then [] else foldl1 (zipWith (+)) xs). zipWith (\vs u -> map (u*) vs) vss) uss
matI::(Num a) => Int -> [[a]]
matI n = [ [fromIntegral.fromEnum $ i == j | j <- [1..n]] | i <- [1..n]]
inversion xs = gauss xs (matI $ length xs)
gauss::[[Double]] -> [[Double]] -> [[Double]]
gauss xs bs = map (map fromRational) $ solveGauss (toR xs) (toR bs)
where toR = map $ map toRational
solveGauss:: (Fractional a, Ord a) => [[a]] -> [[a]] -> [[a]]
solveGauss xs bs | null xs || null bs || length xs /= length bs || (not $ isSquareMatrix xs) || (not $ isMatrix bs) = []
| otherwise = uncurry solveTriangle $ triangle xs bs
solveTriangle::(Fractional a,Eq a) => [[a]] -> [[a]] -> [[a]]
solveTriangle us _ | not.null.dropWhile ((/= 0).head) $ us = []
solveTriangle ([c]:as) (b:bs) = go as bs [map (/c) b]
where
val us vs ws = let u = head us in map (/u) $ zipWith (-) vs (head $ mult [tail us] ws)
go [] _ zs = zs
go _ [] zs = zs
go (x:xs) (y:ys) zs = go xs ys $ (val x y zs):zs
triangle::(Num a, Ord a) => [[a]] -> [[a]] -> ([[a]],[[a]])
triangle xs bs = triang ([],[]) (xs,bs)
where
triang ts (_,[]) = ts
triang ts ([],_) = ts
triang (os,ps) zs = triang (us:os,cs:ps).unzip $ [(fun tus vs, fun cs es) | (v:vs,es) <- zip uss css,let fun = zipWith (\x y -> v*x - u*y)]
where ((us@(u:tus)):uss,cs:css) = bubble zs
bubble::(Num a, Ord a) => ([[a]],[[a]]) -> ([[a]],[[a]])
bubble (xs,bs) = (go xs, go bs)
where
idmax = snd.maximum.flip zip [0..].map (abs.head) $ xs
go ys = let (us,vs) = splitAt idmax ys in vs ++ us
main = do
let a = [[1, 2, 3], [4, 1, 6], [7, 8, 9]]
let b = [[2, -1, 0], [-1, 2, -1], [0, -1, 2]]
putStrLn "inversion a ="
mapM_ print $ inversion a
putStrLn "\ninversion b ="
mapM_ print $ inversion b</lang>
{{out}}
<pre>
inversion a =
[-0.8125,0.125,0.1875]
[0.125,-0.25,0.125]
[0.5208333333333334,0.125,-0.14583333333333334]
inversion b =
[0.75,0.5,0.25]
[0.5,1.0,0.5]
[0.25,0.5,0.75]
</pre>
=={{header|J}}==
| |||