1 files changed, 102 insertions, 0 deletions

diff --git a/lib/Numeric/LinearAlgebra/Linear.hs b/lib/Numeric/LinearAlgebra/Linear.hs
new file mode 100644
index 0000000..a3a0183
--- /dev/null
+++ b/lib/Numeric/LinearAlgebra/Linear.hs
@@ -0,0 +1,102 @@
+{-# OPTIONS_GHC -fglasgow-exts #-}
+-----------------------------------------------------------------------------
+{- |
+Module      :  Numeric.LinearAlgebra.Linear
+Copyright   :  (c) Alberto Ruiz 2006-7
+License     :  GPL-style
+
+Maintainer  :  Alberto Ruiz (aruiz at um dot es)
+Stability   :  provisional
+Portability :  uses ffi
+
+Basic optimized operations on vectors and matrices.
+
+-}
+-----------------------------------------------------------------------------
+
+module Numeric.LinearAlgebra.Linear (
+    Linear(..),
+    multiply, dot, outer
+) where
+
+
+import Data.Packed.Internal
+import Data.Packed
+import Numeric.GSL.Vector
+import Complex
+
+-- | basic optimized operations
+class (Container c e) => Linear c e where
+    scale       :: e -> c e -> c e
+    addConstant :: e -> c e -> c e
+    add         :: c e -> c e -> c e
+    sub         :: c e -> c e -> c e
+    -- | element by element multiplication
+    mul         :: c e -> c e -> c e
+    -- | element by element division
+    divide      :: c e -> c e -> c e
+    -- | scale the element by element reciprocal of the object: @scaleRecip 2 (fromList [5,i]) == 2 |> [0.4 :+ 0.0,0.0 :+ (-2.0)]@
+    scaleRecip  :: e -> c e -> c e
+    equal       :: c e -> c e -> Bool
+--  numequal    :: Double -> c e -> c e -> Bool
+
+instance Linear Vector Double where
+    scale = vectorMapValR Scale
+    scaleRecip = vectorMapValR Recip
+    addConstant = vectorMapValR AddConstant
+    add = vectorZipR Add
+    sub = vectorZipR Sub
+    mul = vectorZipR Mul
+    divide = vectorZipR Div
+    equal u v = dim u == dim v && vectorMax (vectorMapR Abs (sub u v)) == 0.0
+
+instance Linear Vector (Complex Double) where
+    scale = vectorMapValC Scale
+    scaleRecip = vectorMapValC Recip
+    addConstant = vectorMapValC AddConstant
+    add = vectorZipC Add
+    sub = vectorZipC Sub
+    mul = vectorZipC Mul
+    divide = vectorZipC Div
+    equal u v = dim u == dim v && vectorMax (liftVector magnitude (sub u v)) == 0.0
+
+instance Linear Matrix Double where
+    scale x = liftMatrix (scale x)
+    scaleRecip x = liftMatrix (scaleRecip x)
+    addConstant x = liftMatrix (addConstant x)
+    add = liftMatrix2 add
+    sub = liftMatrix2 sub
+    mul = liftMatrix2 mul
+    divide = liftMatrix2 divide
+    equal a b = cols a == cols b && cdat a `equal` cdat b
+
+
+instance Linear Matrix (Complex Double) where
+    scale x = liftMatrix (scale x)
+    scaleRecip x = liftMatrix (scaleRecip x)
+    addConstant x = liftMatrix (addConstant x)
+    add = liftMatrix2 add
+    sub = liftMatrix2 sub
+    mul = liftMatrix2 mul
+    divide = liftMatrix2 divide
+    equal a b = cols a == cols b && cdat a `equal` cdat b
+
+--------------------------------------------------
+
+-- | euclidean inner product
+dot :: (Field t) => Vector t -> Vector t -> t
+dot u v = dat (multiply r c) `at` 0
+    where r = asRow u
+          c = asColumn v
+
+
+{- | Outer product of two vectors.
+
+@\> 'fromList' [1,2,3] \`outer\` 'fromList' [5,2,3]
+(3><3)
+ [  5.0, 2.0, 3.0
+ , 10.0, 4.0, 6.0
+ , 15.0, 6.0, 9.0 ]@
+-}
+outer :: (Field t) => Vector t -> Vector t -> Matrix t
+outer u v = asColumn u `multiply` asRow v