{-

A demonstration of LU decomposition and matrix inversion of complex
matrices using exact rational arithmetic.

Copyright 2015 Ken Takusagawa

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or (at
your option) any later version.

This program is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>.
-}
{-# LANGUAGE LambdaCase,ScopedTypeVariables #-}
module Main where {
import Data.Matrix;
import Data.ComplexNum;
import Data.Ratio;
import Data.List;
import System.Random;
import Control.Monad;
import System.Environment(getArgs);
import Data.Maybe;

digitmax :: Integer;
digitmax = 9;

main::IO();
main= getArgs >>= \case {
["example",n] -> example $ read n;
["benchinv",n] -> benchmark_invert $ read n;
["benchlu",n] -> benchmark_lu $ read n;
["inv",n] -> invexample $ read n;
["hilberttest",n] -> hilberttest $ read n;
["hilbertinv",n] -> putStrLn $ showmatrixgeneric showrational $ fromJust $ invert $ hilbert $ read n;
_ -> error "need arguments";
};

type T = Complex Rational;
type M = Matrix T;

{-
m :: M;
m = fromList 3 3 $ map (\x -> x%1:+0) [3,1,4,1,5,9,2,6,5];

q :: M;
q = fromList 3 3 $ map (\x -> 0:+(x%1)) [2,7,1,8,2,8,1,8,3];
-}


showcomplex :: Complex Rational -> String;
showcomplex (x:+y) = if y==0 then showrational x
else if x==0 then showrational y ++ "*i"
else showrational x ++ (if y>0 then "+" else "") ++ showrational y ++"*i";

showrational :: Rational -> String;
showrational 0 = "0";
showrational r = if denominator r == 1 then show $ numerator r
else (show $ numerator r)++"/"++(show $ denominator r);

showmatrixgeneric :: (a -> String) -> Matrix a -> String;
showmatrixgeneric showval x = "[ " ++ (concat $ intersperse " ; " $ map (unwords . map showval) $ toLists x) ++ " ]";

showmatrix :: M -> String;
showmatrix = showmatrixgeneric showcomplex;

outLU :: Maybe (M,M,M,T) -> String;
outLU Nothing = "singular\n";
outLU (Just (u,l,p,d)) = unlines
[ "L= "++showmatrix l
, "U= "++showmatrix u
, "P= "++showmatrix p -- relying on pretty printer to do the right thing for an all real matrix
, "d= "++showcomplex d
];


randomComplexDigit :: RandomGen g => g -> (T, g);
randomComplexDigit g = case randomR (0,digitmax) g of {
(d,g2) -> let { bg = random g2 } in
(if fst bg then d%1:+0 else 0:+d%1, snd bg)
};

randommatrix :: Int -> Int -> IO M;
randommatrix rows cols = replicateM (rows*cols) (getStdRandom randomComplexDigit) >>= return . fromList (fromIntegral rows) (fromIntegral cols);

example :: Int -> IO ();
example size = do {
a <- randommatrix size size;
putStrLn $ "A= "++showmatrix a;
let { inv = fromJust $ invert a
; lu = luDecompWithMag magnitude2 a;
; (u,_,_,d) = fromJust $ lu
; det = d* diagProd u
}
; putStr $ outLU lu;
-- putStrLn $ "Ainv= " ++ showmatrix inv;
putStrLn $ "Ainv= " ++ (showmatrix $ scaleMatrix det inv) ++ " / (" ++ showcomplex det ++ ")";
print $ a*inv == identity (fromIntegral size);
};

benchmark_invert :: Int -> IO();
benchmark_invert size = do {
a <- randommatrix size size;
case invert a of {
Nothing -> putStrLn $ "singular= " ++ showmatrix a;
Just inv -> print $ a*inv == identity (fromIntegral size);
}};

benchmark_lu :: Int -> IO ();
benchmark_lu size = do {
a <- randommatrix size size;
case luDecompWithMag magnitude2 a of {
Nothing -> putStrLn $ "singular= " ++ showmatrix a;
Just (u,l,p,_d) -> print $ p*a == l*u;
}};


invexample_det :: Int -> IO ();
invexample_det size = do {
a <- randommatrix size size;
putStrLn $ "A= "++showmatrix a;
case invert a of {
Nothing -> putStrLn "singular";
Just inv -> do {
let { lu = luDecompWithMag magnitude2 a;
; (u,_,_,d) = fromJust $ lu
; det = d* diagProd u
};
putStrLn $ "Ainv= " ++ (showmatrix $ scaleMatrix det inv) ++ " / (" ++ showcomplex det ++ ")";
print $ a*inv == identity (fromIntegral size);
}}};

invexample :: Int -> IO ();
invexample size = do {
a <- randommatrix size size;
putStrLn $ "A= "++showmatrix a;
case invert a of {
Nothing -> putStrLn "singular";
Just inv -> putStrLn $ "Ainv= " ++ showmatrix inv
}};

hilbert :: Int -> Matrix Rational;
hilbert n = matrix (fromIntegral n) (fromIntegral n) $ \(i,j) -> 1% (fromIntegral $ i+j-1);

hilberttest :: Int -> IO ();
hilberttest n = do {
let { a = hilbert n
; ainv = fromJust $ invert a };
print $ a*ainv == identity (fromIntegral n);
};

} --end;

{-
runtime

$ for i in `seq 50 5 100`; do echo $i
> /usr/bin/time ./matrixlu $i > /dev/null
> done
50
2.51user 0.11system 0:02.63elapsed 99%CPU (0avgtext+0avgdata 406304maxresident)k
3040inputs+0outputs (19major+101706minor)pagefaults 0swaps
55
2.81user 0.13system 0:02.94elapsed 100%CPU (0avgtext+0avgdata 493848maxresident)k
0inputs+0outputs (0major+123602minor)pagefaults 0swaps
60
4.09user 0.20system 0:04.29elapsed 100%CPU (0avgtext+0avgdata 821020maxresident)k
0inputs+0outputs (0major+205405minor)pagefaults 0swaps
65
5.96user 0.27system 0:06.23elapsed 100%CPU (0avgtext+0avgdata 1008932maxresident)k
0inputs+0outputs (0major+252380minor)pagefaults 0swaps
70
8.54user 0.38system 0:08.92elapsed 100%CPU (0avgtext+0avgdata 1603868maxresident)k
0inputs+0outputs (0major+401091minor)pagefaults 0swaps
75
11.19user 0.41system 0:11.59elapsed 100%CPU (0avgtext+0avgdata 1820460maxresident)k
0inputs+0outputs (0major+455263minor)pagefaults 0swaps
80
14.50user 0.53system 0:15.03elapsed 100%CPU (0avgtext+0avgdata 2140968maxresident)k
0inputs+0outputs (0major+535392minor)pagefaults 0swaps
85
19.87user 0.83system 0:20.69elapsed 100%CPU (0avgtext+0avgdata 3476004maxresident)k
0inputs+0outputs (0major+869136minor)pagefaults 0swaps
90
25.24user 0.99system 0:26.22elapsed 100%CPU (0avgtext+0avgdata 4164132maxresident)k
0inputs+0outputs (0major+1041149minor)pagefaults 0swaps
95
31.86user 1.30system 0:33.15elapsed 100%CPU (0avgtext+0avgdata 5196588maxresident)k
0inputs+0outputs (0major+1299297minor)pagefaults 0swaps
100
42.76user 1.20system 0:44.02elapsed 99%CPU (0avgtext+0avgdata 5362476maxresident)k
0inputs+0outputs (0major+1340771minor)pagefaults 0swaps


-}