dispf, disps

3 files changed, 54 insertions, 1 deletions

diff --git a/lib/Data/Packed/Matrix.hs b/lib/Data/Packed/Matrix.hs
index 6a6942d..7e197d1 100644
--- a/lib/Data/Packed/Matrix.hs
+++ b/lib/Data/Packed/Matrix.hs
@@ -29,7 +29,7 @@ module Data.Packed.Matrix (
     extractRows,
     ident, diag, diagRect, takeDiag,
     liftMatrix, liftMatrix2,
-    format,
+    format, dispf, disps, vecdisp,
     loadMatrix, saveMatrix, fromFile, fileDimensions,
     readMatrix, fromArray2D
 ) where
@@ -40,6 +40,7 @@ import Data.Packed.Vector
 import Data.List(transpose,intersperse)
 import Data.Array
 import System.Process(readProcess)
+import Text.Printf(printf)
 
 -- | creates a matrix from a vertical list of matrices
 joinVert :: Element t => [Matrix t] -> Matrix t
@@ -245,6 +246,43 @@ dispC :: Int -> Matrix (Complex Double) -> IO ()
 dispC d m = disp m (shfc d)
 -}
 
+-------------------------------------------------------------------
+-- display utilities
+
+-- | Print a matrix with \"autoscaling\" and a given number of decimal places.
+disps :: Int -> Matrix Double -> String
+disps d x = sdims x ++ "  " ++ formatScaled d x
+
+-- | Print a matrix with a given number of decimal places.
+dispf :: Int -> Matrix Double -> String
+dispf d x = sdims x ++ "\n" ++ formatFixed (if isInt x then 0 else d) x
+
+sdims x = show (rows x) ++ "x" ++ show (cols x)
+
+formatFixed d x = format "  " (printf ("%."++show d++"f")) $ x
+
+isInt = all lookslikeInt . toList . flatten where
+    lookslikeInt x = show (round x :: Int) ++".0" == shx || "-0.0" == shx
+        where shx = show x
+
+formatScaled dec t = "E"++show o++"\n" ++ ss
+    where ss = format " " (printf fmt. g) t
+          g x | o >= 0    = x/10^(o::Int)
+              | otherwise = x*10^(-o)
+          o = floor $ maximum $ map (logBase 10 . abs) $ toList $ flatten t
+          fmt = '%':show (dec+3) ++ '.':show dec ++"f"
+
+-- | Print a vector using a function for printing matrices.
+vecdisp :: (Element t) => (Matrix t -> String) -> Vector t -> String
+vecdisp f v
+    = ((show (dim v) ++ " |> ") ++) . (++"\n")
+    . unwords . lines .  tail . dropWhile (not . (`elem` " \n"))
+    . f . trans . reshape 1
+    $ v
+
+
+--------------------------------------------------------------------
+
 -- | reads a matrix from a string containing a table of numbers.
 readMatrix :: String -> Matrix Double
 readMatrix = fromLists . map (map read). map words . filter (not.null) . lines
diff --git a/lib/Data/Packed/Vector.hs b/lib/Data/Packed/Vector.hs
index 457da71..ed0e49c 100644
--- a/lib/Data/Packed/Vector.hs
+++ b/lib/Data/Packed/Vector.hs
@@ -32,6 +32,10 @@ import Numeric.GSL.Vector
 
 @\> linspace 5 (-3,7)
 5 |> [-3.0,-0.5,2.0,4.5,7.0]@
+
+Logarithmic spacing can be defined as follows:
+
+@logspace n (a,b) = 10 ** linspace n (a,b)@
 -}
 linspace :: Int -> (Double, Double) -> Vector Double
 linspace n (a,b) = add a $ scale s  $ fromList [0 .. fromIntegral n-1]
diff --git a/lib/Numeric/LinearAlgebra/Tests.hs b/lib/Numeric/LinearAlgebra/Tests.hs
index 9557ac3..063ae86 100644
--- a/lib/Numeric/LinearAlgebra/Tests.hs
+++ b/lib/Numeric/LinearAlgebra/Tests.hs
@@ -321,6 +321,8 @@ manymult n = foldl1' (<>) (map rot2 angles) where
         where c = cos a
               s = sin a
 
+multb n = foldl1' (<>) (replicate (10^6) (ident n :: Matrix Double))
+
 --------------------------------
 
 multBench = do
@@ -328,6 +330,15 @@ multBench = do
     let b = ident 2000 :: Matrix Double
     a `seq` b `seq` putStrLn ""
     time "product of 1M different 3x3 matrices" (manymult (10^6))
+    putStrLn ""
+    time "product of 1M constant  1x1 matrices" (multb 1)
+    time "product of 1M constant  3x3 matrices" (multb 3)
+    --time "product of 1M constant  5x5 matrices" (multb 5)
+    time "product of 1M const.  10x10 matrices" (multb 10)
+    --time "product of 1M const.  15x15 matrices" (multb 15)
+    time "product of 1M const.  20x20 matrices" (multb 20)
+    --time "product of 1M const.  25x25 matrices" (multb 25)
+    putStrLn ""
     time "product (1000 x 1000)<>(1000 x 1000)" (a<>a)
     time "product (2000 x 2000)<>(2000 x 2000)" (b<>b)