9 files changed, 94 insertions, 47 deletions

diff --git a/lib/GSL.hs b/lib/GSL.hs
index d65f8ff..f04cf26 100644
--- a/lib/GSL.hs
+++ b/lib/GSL.hs
@@ -20,6 +20,7 @@ module GSL (
 , module GSL.Minimization
 , module GSL.Special
 , module Complex
+, setErrorHandlerOff
 ) where
 
 import GSL.Integration
@@ -30,3 +31,8 @@ import GSL.Polynomials
 import GSL.Minimization
 import Complex
 import GSL.Special
+
+
+-- | This action removes the GSL default error handler (which aborts the program), so that
+-- GSL errors can be handled by Haskell (using Control.Exception) and ghci doesn't abort.
+foreign import ccall "GSL/gsl-aux.h no_abort_on_error" setErrorHandlerOff :: IO ()
diff --git a/lib/GSL/Matrix.hs b/lib/GSL/Matrix.hs
index ec8ceea..c1bce37 100644
--- a/lib/GSL/Matrix.hs
+++ b/lib/GSL/Matrix.hs
@@ -16,7 +16,7 @@ module GSL.Matrix(
     eigSg, eigHg,
     svdg,
     qr,
-    chol,
+    cholR, -- cholC,
     luSolveR, luSolveC,
     luR, luC,
     fromFile, extractRows
@@ -153,24 +153,29 @@ foreign import ccall "gsl-aux.h QR" c_qr :: TMMM
 
 {- | Cholesky decomposition of a symmetric positive definite real matrix using /gsl_linalg_cholesky_decomp/.
 
-@\> let c = chol $ 'fromLists' [[5,4],[4,5]]
-\
-\> c
-2.236    0.
-1.789 1.342
-\
-\> c \<\> 'trans' c
-5.000 4.000
-4.000 5.000@
+@\> chol $ (2><2) [1,2,
+                   2,9::Double]
+(2><2)
+ [ 1.0,              0.0
+ , 2.0, 2.23606797749979 ]@
 
 -}
-chol :: Matrix Double -> Matrix Double
-chol x = unsafePerformIO $ do
+cholR :: Matrix Double -> Matrix Double
+cholR x = unsafePerformIO $ do
     res <- createMatrix RowMajor r r
-    c_chol // mat cdat x // mat dat res // check "chol" [cdat x]
+    c_cholR // mat cdat x // mat dat res // check "cholR" [cdat x]
     return res
   where r = rows x
-foreign import ccall "gsl-aux.h chol" c_chol :: TMM
+foreign import ccall "gsl-aux.h cholR" c_cholR :: TMM
+
+cholC :: Matrix (Complex Double) -> Matrix (Complex Double)
+cholC x = unsafePerformIO $ do
+    res <- createMatrix RowMajor r r
+    c_cholC // mat cdat x // mat dat res // check "cholC" [cdat x]
+    return res
+  where r = rows x
+foreign import ccall "gsl-aux.h cholC" c_cholC :: TCMCM
+
 
 --------------------------------------------------------
 
diff --git a/lib/GSL/Special.hs b/lib/GSL/Special.hs
index bdee1b2..fa002b9 100644
--- a/lib/GSL/Special.hs
+++ b/lib/GSL/Special.hs
@@ -41,7 +41,6 @@ module GSL.Special (
 , module GSL.Special.Synchrotron
 , module GSL.Special.Trig
 , module GSL.Special.Zeta
-, setErrorHandlerOff
 )
 where
 
@@ -73,8 +72,3 @@ import GSL.Special.Psi
 import GSL.Special.Synchrotron
 import GSL.Special.Trig
 import GSL.Special.Zeta
-
-
--- | This action removes the GSL default error handler which aborts the program, so that
--- GSL errors can be handled by Haskell (using Control.Exception) and ghci doesn't abort.
-foreign import ccall "gsl-aux.h no_abort_on_error" setErrorHandlerOff :: IO ()
diff --git a/lib/GSL/gsl-aux.c b/lib/GSL/gsl-aux.c
index 0e8406c..c602d5e 100644
--- a/lib/GSL/gsl-aux.c
+++ b/lib/GSL/gsl-aux.c
@@ -471,9 +471,9 @@ int QR(KRMAT(x),RMAT(q),RMAT(r)) {
     OK
 }
 
-int chol(KRMAT(x),RMAT(l)) {
+int cholR(KRMAT(x),RMAT(l)) {
     REQUIRES(xr==xc && lr==xr && lr==lc,BAD_SIZE);
-    DEBUGMSG("chol");
+    DEBUGMSG("cholR");
     KDMVIEW(x);
     DMVIEW(l);
     gsl_matrix_memcpy(M(l),M(x));
@@ -488,6 +488,23 @@ int chol(KRMAT(x),RMAT(l)) {
     OK
 }
 
+int cholC(KCMAT(x),CMAT(l)) {
+    REQUIRES(xr==xc && lr==xr && lr==lc,BAD_SIZE);
+    DEBUGMSG("cholC");
+    KCMVIEW(x);
+    CMVIEW(l);
+    gsl_matrix_complex_memcpy(M(l),M(x));
+    int res = 0; // gsl_linalg_complex_cholesky_decomp(M(l));
+    CHECK(res,res);
+    gsl_complex zero = {0.,0.};
+    int r,c;
+    for (r=0; r<xr-1; r++) {
+        for(c=r+1; c<xc; c++) {
+            lp[r*lc+c] = zero;
+        }
+    }
+    OK
+}
 
 int fft(int code, KCVEC(X), CVEC(R)) {
     REQUIRES(Xn == Rn,BAD_SIZE);
diff --git a/lib/GSL/gsl-aux.h b/lib/GSL/gsl-aux.h
index 8f8618a..dd1a247 100644
--- a/lib/GSL/gsl-aux.h
+++ b/lib/GSL/gsl-aux.h
@@ -38,7 +38,9 @@ int eigensystemC(KCMAT(x),RVEC(l),CMAT(v));
 
 int QR(KRMAT(x),RMAT(q),RMAT(r));
 
-int chol(KRMAT(x),RMAT(l));
+int cholR(KRMAT(x),RMAT(l));
+
+int cholC(KCMAT(x),CMAT(l));
 
 int fft(int code, KCVEC(a), CVEC(b));
 
diff --git a/lib/LinearAlgebra.hs b/lib/LinearAlgebra.hs
index a271592..d8f44b4 100644
--- a/lib/LinearAlgebra.hs
+++ b/lib/LinearAlgebra.hs
@@ -8,28 +8,22 @@ Maintainer  :  Alberto Ruiz (aruiz at um dot es)
 Stability   :  provisional
 Portability :  uses ffi
 
-Some linear algebra algorithms, implemented by means of BLAS, LAPACK or GSL.
+Basic matrix computations implemented by BLAS, LAPACK and GSL.
 
 -}
 -----------------------------------------------------------------------------
 module LinearAlgebra (
     module Data.Packed.Vector,
     module Data.Packed.Matrix,
-    module LinearAlgebra.Linear,
     module LinearAlgebra.Instances,
     module LinearAlgebra.Interface,
-    module LAPACK,
-    module GSL.Matrix,
     module LinearAlgebra.Algorithms,
     module Complex
 ) where
 
-import LinearAlgebra.Linear
 import LinearAlgebra.Instances
 import LinearAlgebra.Interface
 import LinearAlgebra.Algorithms
-import LAPACK
-import GSL.Matrix
 import Data.Packed.Matrix
 import Data.Packed.Vector
 import Complex
\ No newline at end of file
diff --git a/lib/LinearAlgebra/Algorithms.hs b/lib/LinearAlgebra/Algorithms.hs
index 7953386..162426f 100644
--- a/lib/LinearAlgebra/Algorithms.hs
+++ b/lib/LinearAlgebra/Algorithms.hs
@@ -9,23 +9,47 @@ Maintainer  :  Alberto Ruiz (aruiz at um dot es)
 Stability   :  provisional
 Portability :  uses ffi
 
+A simple interface to the available matrix computations. We have defined some generic functions 
+on both real and complex matrices in such a way that higher level algorithms and 
+testing properties are valid for both base types.
+
+In any case the specific functions for particular base types can also be explicitly
+imported from the LAPACK and GSL.Matrix modules.
 
 -}
 -----------------------------------------------------------------------------
 
 module LinearAlgebra.Algorithms (
-    GMatrix(..),
+-- * Basic Linear Algebra
+    scale,
+    add,
+    multiply, dot, outer,
+    linearSolve, linearSolveSVD,
+-- * Matrix factorizations
+    svd, lu, eig, eigSH,
+    qr, chol,
+-- * Utilities
     Normed(..), NormType(..),
     det,inv,pinv,full,economy,
     pinvTol,
 --    pinvTolg,
     nullspacePrec,
     nullVector,
-    eps, i
+-- * Conversions
+    toComplex, fromComplex,
+    comp, real, complex,
+    conj, ctrans,
+-- * Misc
+    eps, i,
+    scaleRecip,
+    addConstant,
+    sub,
+    mul,
+    divide,
 ) where
 
 
-import Data.Packed.Internal
+import Data.Packed.Internal hiding (fromComplex, toComplex, comp, conj)
 import Data.Packed.Matrix
 import GSL.Matrix
 import GSL.Vector
@@ -33,7 +57,8 @@ import LAPACK
 import Complex
 import LinearAlgebra.Linear
 
-class (Linear Matrix t) => GMatrix t where
+-- | Base types for which some optimized matrix computations are available
+class (Linear Matrix t) => Optimized t where
     svd         :: Matrix t -> (Matrix t, Vector Double, Matrix t)
     lu          :: Matrix t -> (Matrix t, Matrix t, [Int], t)
     linearSolve :: Matrix t -> Matrix t -> Matrix t
@@ -41,8 +66,9 @@ class (Linear Matrix t) => GMatrix t where
     ctrans :: Matrix t -> Matrix t
     eig         :: Matrix t -> (Vector (Complex Double), Matrix (Complex Double))
     eigSH       :: Matrix t -> (Vector Double, Matrix t)
+    chol        :: Matrix t -> Matrix t
 
-instance GMatrix Double where
+instance Optimized Double where
     svd = svdR
     lu  = luR
     linearSolve = linearSolveR
@@ -50,8 +76,9 @@ instance GMatrix Double where
     ctrans = trans
     eig = eigR
     eigSH = eigS
+    chol = cholR
 
-instance GMatrix (Complex Double) where
+instance Optimized (Complex Double) where
     svd = svdC
     lu  = luC
     linearSolve = linearSolveC
@@ -59,19 +86,20 @@ instance GMatrix (Complex Double) where
     ctrans = conjTrans
     eig = eigC
     eigSH = eigH
+    chol = error "cholC not yet implemented" -- waiting for GSL-1.10
 
 square m = rows m == cols m
 
-det :: GMatrix t => Matrix t -> t
+det :: Optimized t => Matrix t -> t
 det m | square m = s * (product $ toList $ takeDiag $ u)
       | otherwise = error "det of nonsquare matrix"
     where (_,u,_,s) = lu m
 
-inv :: GMatrix t => Matrix t -> Matrix t
+inv :: Optimized t => Matrix t -> Matrix t
 inv m | square m = m `linearSolve` ident (rows m)
       | otherwise = error "inv of nonsquare matrix"
 
-pinv :: GMatrix t => Matrix t -> Matrix t
+pinv :: Optimized t => Matrix t -> Matrix t
 pinv m = linearSolveSVD m (ident (rows m))
 
 
@@ -119,15 +147,15 @@ i = 0:+1
 
 
 -- | matrix product
-mXm :: (Num t, Field t) => Matrix t -> Matrix t -> Matrix t
+mXm :: (Num t, Optimized t) => Matrix t -> Matrix t -> Matrix t
 mXm = multiply
 
 -- | matrix - vector product
-mXv :: (Num t, Field t) => Matrix t -> Vector t -> Vector t
+mXv :: (Num t, Optimized t) => Matrix t -> Vector t -> Vector t
 mXv m v = flatten $ m `mXm` (asColumn v)
 
 -- | vector - matrix product
-vXm :: (Num t, Field t) => Vector t -> Matrix t -> Vector t
+vXm :: (Num t, Optimized t) => Vector t -> Matrix t -> Vector t
 vXm v m = flatten $ (asRow v) `mXm` m
 
 
@@ -171,7 +199,7 @@ pnormCM Infinity m = vectorMax $ liftMatrix (liftVector magnitude) m `mXv` const
 --pnormCM _ _ = error "p norm not yet defined"
 
 -- -- | computes the p-norm of a matrix or vector (with the same definitions as GNU-octave). pnorm 0 denotes \\inf-norm. See also 'norm'.
---pnorm :: (Container t, Field a) => Int -> t a -> Double
+--pnorm :: (Container t, Optimized a) => Int -> t a -> Double
 --pnorm = pnormG
 
 class Normed t where
@@ -194,7 +222,7 @@ instance Normed (Matrix (Complex Double)) where
 -----------------------------------------------------------------------
 
 -- | The nullspace of a matrix from its SVD decomposition.
-nullspacePrec :: GMatrix t
+nullspacePrec :: Optimized t
               => Double          -- ^ relative tolerance in 'eps' units
               -> Matrix t   -- ^ input matrix
               -> [Vector t] -- ^ list of unitary vectors spanning the nullspace
@@ -207,7 +235,7 @@ nullspacePrec t m = ns where
     ns = drop rank $ toRows $ ctrans v
 
 -- | The nullspace of a matrix, assumed to be one-dimensional, with default tolerance (shortcut for @last . nullspacePrec 1@).
-nullVector :: GMatrix t => Matrix t -> Vector t
+nullVector :: Optimized t => Matrix t -> Vector t
 nullVector = last . nullspacePrec 1
 
 ------------------------------------------------------------------------
diff --git a/lib/LinearAlgebra/Interface.hs b/lib/LinearAlgebra/Interface.hs
index 0d2c0a0..3392d54 100644
--- a/lib/LinearAlgebra/Interface.hs
+++ b/lib/LinearAlgebra/Interface.hs
@@ -82,7 +82,7 @@ infixl 3 <->
 
 {- | Horizontal concatenation of matrices and vectors:
 
-@> (ident 3 -&- 3 * ident 3) |&| fromList [1..6.0]
+@> (ident 3 \<-\> 3 * ident 3) \<|\> fromList [1..6.0]
 (6><4)
  [ 1.0, 0.0, 0.0, 1.0
  , 0.0, 1.0, 0.0, 2.0
diff --git a/lib/LinearAlgebra/Linear.hs b/lib/LinearAlgebra/Linear.hs
index 2f1bc6f..85daa4a 100644
--- a/lib/LinearAlgebra/Linear.hs
+++ b/lib/LinearAlgebra/Linear.hs
@@ -12,6 +12,7 @@ Portability :  uses ffi
 
 -}
 -----------------------------------------------------------------------------
+-- #hide
 
 module LinearAlgebra.Linear (
     Linear(..),
@@ -24,7 +25,7 @@ import Data.Packed.Matrix
 import GSL.Vector
 import Complex
 
-
+-- | basic optimized operations
 class (Field e) => Linear c e where
     scale       :: e -> c e -> c e
     scaleRecip  :: e -> c e -> c e