Anonymous user
LU decomposition: Difference between revisions
→{{header|Idris}}
Line 1,400:
=={{header|Haskell}}==
''Without elem-at-index modifications; doesn't find maximum but any non-zero element''
<lang Haskell>
import Data.List
import Data.Maybe
import Text.Printf
-- a matrix is represented as a list of columns
mmult :: Num a => [[a]] -> [[a]] -> [[a]]
mmult a b = [ [ sum $ zipWith (*) ak bj | ak <- (transpose a) ] | bj <- b ]
nth mA i j = (mA !! j) !! i
idMatrixPart n m k = [ [if (i==j) then 1 else 0 | i <- [1..n]] | j <- [k..m]]
idMatrix n = idMatrixPart n n 1
permMatrix n ix1 ix2 =
[ [ if ((i==ix1 && j==ix2) || (i==ix2 && j==ix1) || (i==j && j /= ix1 && i /= ix2))
then 1 else 0| i <- [0..n-1]] | j <- [0..n-1]]
permMatrix_inv n ix1 ix2 = permMatrix n ix2 ix1
-- count k from zero
elimColumn :: Int -> [[Rational]] -> Int -> [Rational]
elimMatrix :: Int -> [[Rational]] -> Int -> [[Rational]]
elimMatrix_inv :: Int -> [[Rational]] -> Int -> [[Rational]]
elimColumn n mA k = [(let mAkk = (nth mA k k) in if (i>k) then (-(nth mA i k)/mAkk)
else if (i==k) then 1 else 0) | i <- [0..n-1]]
elimMatrix n mA k = (idMatrixPart n k 1) ++ [elimColumn n mA k] ++ (idMatrixPart n n (k+2))
elimMatrix_inv n mA k = (idMatrixPart n k 1) ++ --mA is elimMatrix there
[let c = (mA!!k) in [if (i==k) then 1 else if (i<k) then 0 else (-(c!!i)) | i <- [0..n-1]]]
++ (idMatrixPart n n (k+2))
swapIndx :: [[Rational]] -> Int -> Int
swapIndx mA k = fromMaybe k (findIndex (>0) (drop k (mA!!k)))
-- LUP; lupStep returns [L:U:P]
paStep_recP :: Int -> [[Rational]] -> [[Rational]] -> [[Rational]] -> Int -> [[[Rational]]]
paStep_recM :: Int -> [[Rational]] -> [[Rational]] -> [[Rational]] -> Int -> [[[Rational]]]
lupStep :: Int -> [[Rational]] -> [[[Rational]]]
paStep_recP n mP mA mL cnt =
let mPt = permMatrix n cnt (swapIndx mA cnt) in
let mPtInv = permMatrix_inv n cnt (swapIndx mA cnt) in
if (cnt >= n) then [(mmult mP mL),mA,mP] else
(paStep_recM n (mmult mPt mP) (mmult mPt mA) (mmult mL mPtInv) cnt)
paStep_recM n mP mA mL cnt =
let mMt = elimMatrix n mA cnt in
let mMtInv = elimMatrix_inv n mMt cnt in
paStep_recP n mP (mmult mMt mA) (mmult mL mMtInv) (cnt + 1)
lupStep n mA = paStep_recP n (idMatrix n) mA (idMatrix n) 0
--IO
matrixFromRationalToString m = concat $ intersperse "\n"
(map (\x -> unwords $ printf "%8.4f" <$> (x::[Double]))
(transpose (matrixFromRational m))) where
matrixFromRational m = map (\x -> map fromRational x) m
solveTask mY = let mLUP = lupStep (length mY) mY in
putStrLn ("A: \n" ++ matrixFromRationalToString mY) >>
putStrLn ("L: \n" ++ matrixFromRationalToString (mLUP!!0)) >>
putStrLn ("U: \n" ++ matrixFromRationalToString (mLUP!!1)) >>
putStrLn ("P: \n" ++ matrixFromRationalToString (mLUP!!2)) >>
putStrLn ("Verify: PA\n" ++ matrixFromRationalToString (mmult (mLUP!!2) mY)) >>
putStrLn ("Verify: LU\n" ++ matrixFromRationalToString (mmult (mLUP!!0) (mLUP!!1)))
mY1 = [[1, 2, 1], [3, 4, 7], [5, 7, 0]] :: [[Rational]]
mY2 = [[11, 1, 3, 2], [9, 5, 17, 5], [24, 2, 18, 7], [2, 6, 1, 1]] :: [[Rational]]
main = putStrLn "Task1: \n" >> solveTask mY1 >>
putStrLn "Task2: \n" >> solveTask mY2
</lang>
{{out}}
<pre>
Task1:
A:
1.0000 3.0000 5.0000
2.0000 4.0000 7.0000
1.0000 7.0000 0.0000
L:
1.0000 0.0000 0.0000
2.0000 1.0000 0.0000
1.0000 -2.0000 1.0000
U:
1.0000 3.0000 5.0000
0.0000 -2.0000 -3.0000
0.0000 0.0000 -11.0000
P:
1.0000 0.0000 0.0000
0.0000 1.0000 0.0000
0.0000 0.0000 1.0000
Verify: PA
1.0000 3.0000 5.0000
2.0000 4.0000 7.0000
1.0000 7.0000 0.0000
Verify: LU
1.0000 3.0000 5.0000
2.0000 4.0000 7.0000
1.0000 7.0000 0.0000
Task2:
A:
11.0000 9.0000 24.0000 2.0000
1.0000 5.0000 2.0000 6.0000
3.0000 17.0000 18.0000 1.0000
2.0000 5.0000 7.0000 1.0000
L:
1.0000 0.5556 0.2317 0.0000
0.0000 1.0000 0.0000 0.0000
0.0000 1.8889 1.0000 0.0000
0.0000 0.0909 0.0000 1.0000
U:
0.0081 0.0000 0.0000 0.5325
11.0000 9.0000 24.0000 2.0000
-17.7778 0.0000 -27.3333 -2.7778
0.0000 4.1818 -0.1818 5.8182
P:
0.0000 0.0000 0.0000 1.0000
1.0000 0.0000 0.0000 0.0000
0.0000 0.0000 1.0000 0.0000
0.0000 1.0000 0.0000 0.0000
Verify: PA
2.0000 5.0000 7.0000 1.0000
11.0000 9.0000 24.0000 2.0000
3.0000 17.0000 18.0000 1.0000
1.0000 5.0000 2.0000 6.0000
Verify: LU
2.0000 5.0000 7.0000 1.0000
11.0000 9.0000 24.0000 2.0000
3.0000 17.0000 18.0000 1.0000
1.0000 5.0000 2.0000 6.0000
</pre>
=={{header|Idris}}==
'''works with Idris 0.10'''
| |||