initial support for sparse matrix

8 files changed, 278 insertions, 23 deletions

diff --git a/packages/base/hmatrix.cabal b/packages/base/hmatrix.cabal
index 1948897..06c1a3c 100644
--- a/packages/base/hmatrix.cabal
+++ b/packages/base/hmatrix.cabal
@@ -48,6 +48,7 @@ library
                         Numeric.LinearAlgebra.Data
 
                         Numeric.LinearAlgebra.Compat
+                        Data.Packed.Internal.Sparse
                         
     other-modules:      Data.Packed.Internal,
                         Data.Packed.Internal.Common,
diff --git a/packages/base/src/C/lapack-aux.h b/packages/base/src/C/lapack-aux.h
index a3f1899..c95a2a3 100644
--- a/packages/base/src/C/lapack-aux.h
+++ b/packages/base/src/C/lapack-aux.h
@@ -36,6 +36,7 @@ typedef short ftnlen;
 
 /********************************************************/
 
+#define IVEC(A) int A##n, int*A##p
 #define FVEC(A) int A##n, float*A##p
 #define DVEC(A) int A##n, double*A##p
 #define QVEC(A) int A##n, complex*A##p
@@ -47,6 +48,7 @@ typedef short ftnlen;
 #define CMAT(A) int A##r, int A##c, doublecomplex* A##p
 #define PMAT(A) int A##r, int A##c, void* A##p, int A##s
 
+#define KIVEC(A) int A##n, const int*A##p
 #define KFVEC(A) int A##n, const float*A##p
 #define KDVEC(A) int A##n, const double*A##p
 #define KQVEC(A) int A##n, const complex*A##p
diff --git a/packages/base/src/C/vector-aux.c b/packages/base/src/C/vector-aux.c
index 5b9c171..53b56aa 100644
--- a/packages/base/src/C/vector-aux.c
+++ b/packages/base/src/C/vector-aux.c
@@ -744,3 +744,29 @@ int random_vector(int seed, int code, DVEC(r)) {
     }
 }
 
+////////////////////////////////////////////////////////////////////////////////
+
+int smXv(KDVEC(vals),KIVEC(cols),KIVEC(rows),KDVEC(x),DVEC(r)) {
+    int r, c;
+    for (r = 0; r < rowsn - 1; r++) {
+        rp[r] = 0;
+        for (c = rowsp[r]; c < rowsp[r+1]; c++) {
+            rp[r] += valsp[c-1] * xp[colsp[c-1]-1];
+        }
+    }
+    OK
+}
+
+int smTXv(KDVEC(vals),KIVEC(cols),KIVEC(rows),KDVEC(x),DVEC(r)) {
+    int r,c;
+    for (c = 0; c < rn; c++) {
+        rp[c] = 0;
+    }
+    for (r = 0; r < rowsn - 1; r++) {
+        for (c = rowsp[r]; c < rowsp[r+1]; c++) {
+            rp[colsp[c-1]-1] += valsp[c-1] * xp[r];
+        }
+    }
+    OK
+}
+
diff --git a/packages/base/src/Data/Packed/Internal/Numeric.hs b/packages/base/src/Data/Packed/Internal/Numeric.hs
index 81a8083..3528e96 100644
--- a/packages/base/src/Data/Packed/Internal/Numeric.hs
+++ b/packages/base/src/Data/Packed/Internal/Numeric.hs
@@ -20,6 +20,7 @@ module Data.Packed.Internal.Numeric (
     ident, diag, ctrans,
     -- * Generic operations
     Container(..),
+    Transposable(..), Linear(..), Testable(..),
     -- * Matrix product and related functions
     Product(..), udot,
     mXm,mXv,vXm,
@@ -605,3 +606,25 @@ condV f a b l e t = f a' b' l' e' t'
   where
     [a', b', l', e', t'] = conformVs [a,b,l,e,t]
 
+--------------------------------------------------------------------------------
+
+class Transposable t
+  where
+    tr :: t -> t
+
+
+class Linear t v
+  where
+    scalarL :: t -> v
+    addL    :: v -> v -> v
+    scaleL  :: t -> v -> v
+
+
+class Testable t
+  where
+    checkT   :: t -> (Bool, IO())
+    ioCheckT :: t -> IO (Bool, IO())
+    ioCheckT = return . checkT
+
+--------------------------------------------------------------------------------
+
diff --git a/packages/base/src/Data/Packed/Internal/Sparse.hs b/packages/base/src/Data/Packed/Internal/Sparse.hs
new file mode 100644
index 0000000..544c913
--- /dev/null
+++ b/packages/base/src/Data/Packed/Internal/Sparse.hs
@@ -0,0 +1,191 @@
+{-# LANGUAGE RecordWildCards #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE FlexibleInstances #-}
+
+
+
+module Data.Packed.Internal.Sparse(
+    SMatrix(..),
+    mkCSR, mkDiag,
+    AssocMatrix,
+    toDense,
+    smXv
+)where
+
+import Numeric.Container
+import qualified Data.Vector.Storable as V
+import Data.Function(on)
+import Control.Arrow((***))
+import Control.Monad(when)
+import Data.List(groupBy, sort)
+import Foreign.C.Types(CInt(..))
+import Numeric.LinearAlgebra.Devel
+import System.IO.Unsafe(unsafePerformIO)
+import Foreign(Ptr)
+import Text.Printf(printf)
+
+infixl 0 ~!~
+c ~!~ msg = when c (error msg)
+
+type AssocMatrix = [((Int,Int),Double)]
+
+data SMatrix
+    = CSR
+        { csrVals :: Vector Double
+        , csrCols :: Vector CInt
+        , csrRows :: Vector CInt
+        , nRows   :: Int
+        , nCols   :: Int
+        }
+    | CSC
+        { cscVals :: Vector Double
+        , cscRows :: Vector CInt
+        , cscCols :: Vector CInt
+        , nRows   :: Int
+        , nCols   :: Int
+        }
+    | Diag
+        { diagVals :: Vector Double
+        , nRows   :: Int
+        , nCols   :: Int
+        }
+--    | Banded
+    deriving Show
+
+mkCSR :: AssocMatrix -> SMatrix
+mkCSR sm' = CSR{..}
+  where
+    sm = sort sm'
+    rws = map ((fromList *** fromList)
+              . unzip
+              . map ((succ.fi.snd) *** id)
+              )
+        . groupBy ((==) `on` (fst.fst))
+        $ sm
+    rszs = map (fi . dim . fst) rws
+    csrRows = fromList (scanl (+) 1 rszs)
+    csrVals = vjoin (map snd rws)
+    csrCols = vjoin (map fst rws)
+    nRows = dim csrRows - 1
+    nCols = fromIntegral (V.maximum csrCols)
+
+
+mkDiagR r c v
+    | dim v <= min r c = Diag{..}
+    | otherwise = error $ printf "mkDiagR: incorrect sizes (%d,%d) [%d]" r c (dim v)
+  where
+    nRows = r
+    nCols = c
+    diagVals = v
+
+mkDiag v = mkDiagR (dim v) (dim v) v
+
+
+type IV t = CInt -> Ptr CInt   -> t
+type  V t = CInt -> Ptr Double -> t
+type SMxV = V (IV (IV (V (V (IO CInt)))))
+
+smXv :: SMatrix -> Vector Double -> Vector Double
+smXv CSR{..} v = unsafePerformIO $ do
+    dim v /= nCols ~!~ printf "smXv (CSR): incorrect sizes: (%d,%d) x %d" nRows nCols (dim v)
+    r <- createVector nRows
+    app5 c_smXv vec csrVals vec csrCols vec csrRows vec v vec r "CSRXv"
+    return r
+
+smXv CSC{..} v = unsafePerformIO $ do
+    dim v /= nCols ~!~ printf "smXv (CSC): incorrect sizes: (%d,%d) x %d" nRows nCols (dim v)
+    r <- createVector nRows
+    app5 c_smTXv vec cscVals vec cscRows vec cscCols vec v vec r "CSCXv"
+    return r
+
+smXv Diag{..} v
+    | dim v == nCols
+        = vjoin [ subVector 0 (dim diagVals) v `mul` diagVals
+                , konst 0 (nRows - dim diagVals) ]
+    | otherwise = error $ printf "smXv (Diag): incorrect sizes: (%d,%d) [%d] x %d"
+                                 nRows nCols (dim diagVals) (dim v)
+
+
+instance Contraction SMatrix (Vector Double) (Vector Double)
+  where
+    contraction = smXv
+
+--------------------------------------------------------------------------------
+
+foreign import ccall unsafe "smXv"
+  c_smXv :: SMxV
+
+foreign import ccall unsafe "smTXv"
+  c_smTXv :: SMxV
+
+--------------------------------------------------------------------------------
+
+toDense :: AssocMatrix -> Matrix Double
+toDense asm = assoc (r+1,c+1) 0 asm
+  where
+    (r,c) = (maximum *** maximum) . unzip . map fst $ asm
+
+
+
+instance Transposable (SMatrix)
+  where
+    tr (CSR vs cs rs n m) = CSC vs cs rs m n
+    tr (CSC vs rs cs n m) = CSR vs rs cs m n
+    tr (Diag v n m) = Diag v m n
+
+instance Transposable (Matrix Double)
+  where
+    tr = trans
+
+
+
+--------------------------------------------------------------------------------
+
+instance Testable SMatrix
+  where
+    checkT _ = (ok,info)
+      where
+        sma = convo2 20 3
+        x1 = vect [1..20]
+        x2 = vect [1..40]
+        sm = mkCSR sma
+
+        s1 = sm ◇ x1
+        d1 = toDense sma  ◇  x1
+
+        s2 = tr sm  ◇  x2
+        d2 = tr (toDense sma) ◇  x2
+
+        sdia = mkDiagR 40 20 (vect [1..10])
+        s3 =    sdia  ◇  x1
+        s4 = tr sdia  ◇  x2
+        ddia = diagRect 0 (vect [1..10])  40 20
+        d3 = ddia  ◇  x1
+        d4 = tr ddia  ◇  x2
+
+        info = do
+            print sm
+            disp (toDense sma)
+            print s1; print d1
+            print s2; print d2
+            print s3; print d3
+            print s4; print d4
+
+        ok = s1==d1
+          && s2==d2
+          && s3==d3
+          && s4==d4
+
+        disp = putStr . dispf 2
+
+        vect = fromList :: [Double] -> Vector Double
+
+        convomat :: Int -> Int -> AssocMatrix
+        convomat n k = [ ((i,j `mod` n),1) | i<-[0..n-1], j <- [i..i+k-1]]
+
+        convo2 :: Int -> Int -> AssocMatrix
+        convo2 n k = m1 ++ m2
+          where
+            m1 = convomat n k
+            m2 = map (((+n) *** id) *** id) m1
+
diff --git a/packages/base/src/Numeric/Container.hs b/packages/base/src/Numeric/Container.hs
index 264a619..0633640 100644
--- a/packages/base/src/Numeric/Container.hs
+++ b/packages/base/src/Numeric/Container.hs
@@ -32,11 +32,11 @@ module Numeric.Container (
     diag, ident,
     ctrans,
     -- * Generic operations
-    Container(..),
+    Container(..), Transposable(..), Linear(..),
     -- * Matrix product
     Product(..), udot, dot, (◇),
     Mul(..),
-    Contraction(..),
+    Contraction(..),(<.>),
     optimiseMult,
     mXm,mXv,vXm,LSDiv(..),
     outer, kronecker,
@@ -55,7 +55,9 @@ module Numeric.Container (
     IndexOf,
     module Data.Complex,
     -- * IO
-    module Data.Packed.IO
+    module Data.Packed.IO,
+    -- * Misc
+    Testable(..)
 ) where
 
 import Data.Packed hiding (stepD, stepF, condD, condF, conjugateC, conjugateQ)
@@ -87,10 +89,9 @@ linspace n (a,b) = addConstant a $ scale s $ fromList $ map fromIntegral [0 .. n
 
 --------------------------------------------------------
 
-class Contraction a b c | a b -> c
-  where
-    infixl 7 <.>
-    {- | Matrix product, matrix - vector product, and dot product
+{- | Matrix product, matrix - vector product, and dot product (equivalent to 'contraction')
+
+(This operator can also be written using the unicode symbol ◇ (25c7).)
 
 Examples:
 
@@ -129,22 +130,28 @@ For complex vectors the first argument is conjugated:
 
 >>> fromList [1,i,1-i] <.> complex a
 fromList [10.0 :+ 4.0,12.0 :+ 4.0,14.0 :+ 4.0,16.0 :+ 4.0]
-
 -}
-    (<.>) :: a -> b -> c
+infixl 7 <.>
+(<.>) :: Contraction a b c => a -> b -> c
+(<.>) = contraction
 
 
+class Contraction a b c | a b -> c
+  where
+    -- | Matrix product, matrix - vector product, and dot product
+    contraction :: a -> b -> c
+
 instance (Product t, Container Vector t) => Contraction (Vector t) (Vector t) t where
-    u <.> v = conj u `udot` v
+    u `contraction` v = conj u `udot` v
 
 instance Product t => Contraction (Matrix t) (Vector t) (Vector t) where
-    (<.>) = mXv
+    contraction = mXv
 
 instance (Container Vector t, Product t) => Contraction (Vector t) (Matrix t) (Vector t) where
-    (<.>) v m = (conj v) `vXm` m
+    contraction v m = (conj v) `vXm` m
 
 instance Product t => Contraction (Matrix t) (Matrix t) (Matrix t) where
-    (<.>) = mXm
+    contraction = mXm
 
 
 --------------------------------------------------------------------------------
@@ -229,10 +236,10 @@ instance Container Matrix e => Build (Int,Int) (e -> e -> e) Matrix e
 
 --------------------------------------------------------------------------------
 
--- | alternative unicode symbol (25c7) for the contraction operator '(\<.\>)'
+-- | alternative unicode symbol (25c7) for 'contraction'
 (◇) :: Contraction a b c => a -> b -> c
 infixl 7 ◇
-(◇) = (<.>)
+(◇) = contraction
 
 -- | dot product: @cdot u v = 'udot' ('conj' u) v@
 dot :: (Container Vector t, Product t) => Vector t -> Vector t -> t
diff --git a/packages/base/src/Numeric/LinearAlgebra.hs b/packages/base/src/Numeric/LinearAlgebra.hs
index 96bf29f..549ebd0 100644
--- a/packages/base/src/Numeric/LinearAlgebra.hs
+++ b/packages/base/src/Numeric/LinearAlgebra.hs
@@ -37,11 +37,11 @@ module Numeric.LinearAlgebra (
 
     -- * Matrix product
     (<.>),
-    
-    -- | This operator can also be written using the unicode symbol ◇ (25c7).
+
+    -- | The overloaded multiplication operator may need type annotations to remove
+    -- ambiguity. In those cases we can also use the specific functions 'mXm', 'mXv', and 'dot'.
     --
-    
-    -- | The matrix x matrix product is also implemented in the "Data.Monoid" instance, where
+    -- The matrix x matrix product is also implemented in the "Data.Monoid" instance, where
     -- single-element matrices (created from numeric literals or using 'scalar')
     -- are used for scaling.
     --
@@ -52,6 +52,7 @@ module Numeric.LinearAlgebra (
     -- , 4.0, 10.0, 0.0 ]
     --
     -- 'mconcat' uses 'optimiseMult' to get the optimal association order.
+
      
     -- * Other products
     outer, kronecker, cross,
@@ -125,7 +126,7 @@ module Numeric.LinearAlgebra (
     RandDist(..), randomVector, rand, randn, gaussianSample, uniformSample,
     
     -- * Misc
-    meanCov, peps, relativeError, haussholder, optimiseMult, udot, Seed, (◇)
+    meanCov, peps, relativeError, haussholder, optimiseMult, dot, udot, mXm, mXv, smXv, (<>), (◇), Seed, checkT
 ) where
 
 import Numeric.LinearAlgebra.Data
@@ -136,6 +137,5 @@ import Numeric.Container
 import Numeric.LinearAlgebra.Algorithms
 import Numeric.LinearAlgebra.Util
 import Numeric.LinearAlgebra.Random
-
-
+import Data.Packed.Internal.Sparse(smXv)
 
diff --git a/packages/base/src/Numeric/LinearAlgebra/Data.hs b/packages/base/src/Numeric/LinearAlgebra/Data.hs
index 7e8af03..49bc1c0 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Data.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Data.hs
@@ -48,6 +48,10 @@ module Numeric.LinearAlgebra.Data(
     latexFormat,
     dispf, disps, dispcf, format,
 
+    -- * Sparse
+    SMatrix, AssocMatrix, mkCSR, toDense,
+    mkDiag,
+
 -- * Conversion
     Convert(..),
     
@@ -56,7 +60,7 @@ module Numeric.LinearAlgebra.Data(
     rows, cols,
     separable,
 
-    module Data.Complex
+    module Data.Complex,
 
 ) where
 
@@ -65,4 +69,5 @@ import Data.Packed.Matrix
 import Numeric.Container
 import Numeric.LinearAlgebra.Util
 import Data.Complex
+import Data.Packed.Internal.Sparse