diff --git a/INSTALL b/INSTALL
new file mode 100644
--- /dev/null
+++ b/INSTALL
@@ -0,0 +1,58 @@
+-----------------------------------------
+ A simple scientific library for Haskell
+-----------------------------------------
+
+INSTALLATION
+
+Recommended method:
+    $ sudo apt-get install libgsl0-dev refblas3-dev lapack3-dev atlas3-[your arch]-dev
+    $ cabal install hmatrix
+
+Detailed installation instructions:
+    http://www.hmatrix.googlepages.com/installation
+
+INSTALLATION ON WINDOWS ----------------------------------------
+
+1) Download the developer files gsl-1.8-lib.zip from
+   http://gnuwin32.sourceforge.net/packages/gsl.htm
+   and copy the gsl headers folder (under include) to:
+       C:\ghc\ghc.6.x.1\include
+   These headers are also available from:
+       http://perception.inf.um.es/~aruiz/darcs/hmatrix/gsl.zip
+
+2) Copy libgsl.dll, libcblas.dll (from the binaries package gsl-1.8.bin.zip)
+   and liblapack.dll (borrowed from the R system) to the ghc folder, e.g.:
+       C:\ghc\ghc-6.x.x.
+   Rename libcblas.dll to libblas.dll.
+   They are needed to compile programs.
+   These three dlls are available from:
+       http://perception.inf.um.es/~aruiz/darcs/hmatrix/dll1.zip
+
+2.5) Remove the following functions from the export list of
+     lib/Numeric/GSL/Special/Ellint.hs:
+     ellint_Pcomp_e, ellint_Pcomp, ellint_Dcomp_e, ellint_Dcomp
+
+3) Install the package as usual:
+       runhaskell Setup.lhs configure
+       runhaskell Setup.lhs build
+       runhaskell Setup.lhs install
+
+3.5) If configure cannot find ld please see:
+       http://article.gmane.org/gmane.comp.lang.haskell.cafe/32025
+
+4) Copy the dlls available from:
+       http://perception.inf.um.es/~aruiz/darcs/hmatrix/dll2.zip
+   to the working directory or C:\windows\system
+   They are required to run the programs and ghci.
+
+5) run the tests
+
+Unfortunately the lapack dll supplied by the R system does not include
+zgels_, zgelss_, and zgees_, so the functions depending on them
+(linearSolveLS, linearSolveSVD, and schur for complex data)
+will produce a "non supported in this OS" runtime error.
+
+If you find an alternative free and complete lapack.dll which works well
+for this system please let me know.
+
+The examples using graphics do not yet work in windows.
diff --git a/README b/README
--- a/README
+++ b/README
@@ -2,65 +2,16 @@
  A simple scientific library for Haskell
 -----------------------------------------
 
-REQUIREMENTS ----------------------------
-
-1) GNU Scientific Library (http://www.gnu.org/software/gsl).
-   In Ubuntu we need the package "libgsl0-dev".
-
-2) BLAS/LAPACK (http://www.netlib.org/lapack).
-   An optimized implementation is recommended. I have tested:
-
-   - Intel's MKL (http://www.intel.com/cd/software/products).
-     There is a free noncommercial download of MKL for Linux.
-
-   - ATLAS (http://math-atlas.sourceforge.net).
-     In Ubuntu the required packages are "refblas3-dev", "lapack3-dev",
-     and "atlas3-base-dev" (or a version tuned for your machine).
-     Please note that ATLAS currently requires compilation -fviaC in 32bit
-     machines. Otherwise many functions fail, producing strange NaN's.
-     Even with -fvia-C we may get wrong behavior in some cases.
-
-For ghc-6.8.x you may also need:
-
-- libgmp3-dev.
-
-The following packages are used for simple graphics:
-
-- gnuplot
-- imagemagick
-
-GNU-Octave can be used to check if the results
-obtained by this library are correct.
-
-INSTALLATION --------------------------------------
-
-Automatic (using cabal-install and HackageDB):
+INSTALLATION
 
+Recommended method (ok in Ubuntu/Debian systems):
+    $ sudo apt-get install libgsl0-dev refblas3-dev lapack3-dev atlas3-[your_arch]-dev
     $ cabal install hmatrix
 
-Manual:
-
-    Install storable-complex from HackageDB and then
-
-    $ runhaskell Setup.lhs configure --prefix=$HOME --user
-    $ runhaskell Setup.lhs build
-    $ runhaskell Setup.lhs haddock
-    $ runhaskell Setup.lhs install
-
-Using Intel's MKL:
-
-    - add/modify environment variables (e.g. in your .bashrc):
-          export LD_LIBRARY_PATH=/path/to/mkl/lib/arch
-          export    LIBRARY_PATH=/path/to/mkl/lib/arch
-      where arch = "32" or "em64t"
-
-    - add the "-fmkl" flag in the cabal configuration command:
-           $ runhaskell Setup.lhs configure --prefix=$HOME --user -fmkl
-           $ runhaskell Setup.lhs build
-           $ runhaskell Setup.lhs install
-
+Detailed installation instructions:
+    http://www.hmatrix.googlepages.com/installation
 
-See below for installation on Windows.
+For installation in Windows see the companion INSTALL file.
 
 TESTS ---------------------------------------------
 
@@ -70,7 +21,7 @@
 
 Additional tests with big matrices (taking a few minutes):
 
-$ runhaskell examples/experiments bigtests
+$ runhaskell examples/experiments/bigtests
 
 EXAMPLES ------------------------------------------------------
 
@@ -98,12 +49,6 @@
 
 KNOWN PROBLEMS / BUGS -------------------------------
 
-- Compilation with -O -fasm on 32-bit machines produces strange
-  NaN's results on certain blas/lapack calls. In these machines
-  the library is automatically compiled -fvia-C, which apparently
-  solves the problem.
-  On 64-bit, or using MKL, the default and faster -fasm seems to work well.
-
 - On 64-bit machines the example "minimize.hs", when run from ghci,
   produces a segmentation fault. It happens in the call to
   gsl_multimin_fdfminimizer_alloc, inside the C wrapper.
@@ -112,90 +57,12 @@
   program seems to work perfectly well.
 
 - On Ubuntu 6.06 LTS (Dapper) atlas3-sse2-dev (3.6.0-20)
-  produces segmentation faults when working with big matrices 
-  on compiled programs. To expose the problem:
-
-  $ cd examples
-  $ ghc --make -O -fvia-C tests.hs
-  $ ./tests --big
-
-  If this crashes, just uninstall atlas3-sse2 and use atlas3-base-dev instead.
-  Fortunately, atlas3-sse2-dev seems to work well on Ubuntu 7.10 Gutsy.
-  A similar problem was reported at:
-      http://article.gmane.org/gmane.linux.debian.devel.bugs.general/323065
+  produced segmentation faults when working with big matrices
+  on compiled programs.
 
 - On distributions with old GSL versions you should comment out a couple of functions
   in the export lists of Ellint.hs and Debye.hs
 
-CHANGES ---------------------------------------------------------
-
-This is a new version of the library previously known as GSLHaskell.
-It has been renamed to "hmatrix" because only a small part of GSL is actually
-available, and most linear algebra is based on LAPACK.
-
-The code has been extensively refactored. There is a new internal representation
-which admits both C and Fortran matrices and avoids many transposes.
-
-There are only minor API changes:
-
-- The matrix product operator (<>) is now overloaded only for matrix-matrix,
-  matrix-vector and vector-matrix, with the same base type. Dot product and scaling
-  of vectors or matrices is now denoted by `dot` or (<.>) and `scale` or (.*).
-  Conversions from real to complex objects must now be explicit.
-
-- Most linear algebra functions admit both real and complex objects. Utilities such as
-  ident or constant are now polymorphic.
-
-- Runtime errors produced by GSL or LAPACK can be handled using Control.Exeception.catch.
-
-Old GSLHaskell code will work with small modifications.
-
-INSTALLATION ON WINDOWS ----------------------------------------
-
-1) Download the developer files gsl-1.8-lib.zip from
-   http://gnuwin32.sourceforge.net/packages/gsl.htm
-   and copy the gsl headers folder (under include) to:
-       C:\ghc\ghc.6.x.1\include
-   These headers are also available from:
-       http://perception.inf.um.es/~aruiz/darcs/hmatrix/gsl.zip
-
-2) Copy libgsl.dll, libcblas.dll (from the binaries package gsl-1.8.bin.zip)
-   and liblapack.dll (borrowed from the R system) to the ghc folder, e.g.:
-       C:\ghc\ghc-6.x.x.
-   Rename libcblas.dll to libblas.dll.
-   They are needed to compile programs.
-   These three dlls are available from:
-       http://perception.inf.um.es/~aruiz/darcs/hmatrix/dll1.zip
-
-2.5) Remove the following functions from the export list of
-     lib/Numeric/GSL/Special/Ellint.hs:
-     ellint_Pcomp_e, ellint_Pcomp, ellint_Dcomp_e, ellint_Dcomp
-
-3) Install the package as usual:
-       runhaskell Setup.lhs configure
-       runhaskell Setup.lhs build
-       runhaskell Setup.lhs install
-
-3.5) If configure cannot find ld please see:
-       http://article.gmane.org/gmane.comp.lang.haskell.cafe/32025
-
-4) Copy the dlls available from:
-       http://perception.inf.um.es/~aruiz/darcs/hmatrix/dll2.zip
-   to the working directory or C:\windows\system
-   They are required to run the programs and ghci.
-
-5) run the tests
-
-Unfortunately the lapack dll supplied by the R system does not include
-zgels_, zgelss_, and zgees_, so the functions depending on them
-(linearSolveLS, linearSolveSVD, and schur for complex data)
-will produce a "non supported in this OS" runtime error.
-
-If you find an alternative free and complete lapack.dll which works well
-for this system please let me know.
-
-The examples using graphics do not yet work in windows.
-
 ACKNOWLEDGEMENTS -----------------------------------------------------
 
 I thank Don Stewart, Henning Thielemann, Bulat Ziganshin and all the people
@@ -228,3 +95,11 @@
 - Don Stewart fixed the implementation of the internal data structures
   to achieve excellent, C-like performance in Haskell functions which
   explicitly work with the elements of vectors and matrices.
+
+- Dylan Alex Simon improved the numeric instances to allow optimized
+  implementations of signum and abs on Vectors.
+
+- Pedro E. López de Teruel discovered the need of asm("finit") to
+  avoid the wrong NaNs produced by foreign functions.
+
+- Reiner Pope added support for luSolve, based on (d|z)getrs.
diff --git a/examples/benchmarks.hs b/examples/benchmarks.hs
--- a/examples/benchmarks.hs
+++ b/examples/benchmarks.hs
@@ -1,9 +1,8 @@
-{-# OPTIONS -fbang-patterns #-}
+{-# LANGUAGE BangPatterns #-}
 
--- compile as:
--- ghc --make -O2 -optc-O2 -fvia-C benchmarks.hs
--- ghc --make -O benchmarks.hs
+-- $ ghc --make -O2 benchmarks.hs
 
+
 import Numeric.LinearAlgebra
 import System.Time
 import System.CPUTime
@@ -19,18 +18,19 @@
 
 --------------------------------------------------------------------------------
 
-main = sequence_ [bench1,bench2,bench3]
+main = sequence_ [bench1,bench2,bench3,bench4,bench5 1000000 3]
 
 w :: Vector Double
-w = constant 1 30000000
+w = constant 1 5000000
+w2 = 1 * w
 
 bench1 = do
-    putStrLn "Sum of a vector with 30M doubles:"
-    print$ vectorMax w -- evaluate it
-    time $ printf "     BLAS: %.2f: " $ sumVB w
-    time $ printf "  Haskell: %.2f: " $ sumVH w
-    time $ printf "     BLAS: %.2f: " $ sumVB w
-    time $ printf "  Haskell: %.2f: " $ sumVH w
+    putStrLn "Sum of a vector with 5M doubles:"
+    print$ vectorMax (w+w2) -- evaluate it
+    time $ printf "         BLAS: %.2f: " $ sumVB w
+    time $ printf "BLAS only dot: %.2f: " $ w <.> w2
+    time $ printf "      Haskell: %.2f: " $ sumVH w
+    time $ printf "       innerH: %.2f: " $ innerH w w2
 
 sumVB v = constant 1 (dim v) <.> v
 
@@ -41,14 +41,21 @@
        go 0 s = s + (v @> 0)
        go !j !s = go (j - 1) (s + (v @> j))
 
+innerH u v = go (d - 1) 0
+     where
+       d = dim u
+       go :: Int -> Double -> Double
+       go 0 s = s + (u @> 0) * (v @> 0)
+       go !j !s = go (j - 1) (s + (u @> j) * (v @> j))
+
 --------------------------------------------------------------------------------
 
 bench2 = do
     putStrLn "-------------------------------------------------------"
     putStrLn "Multiplication of 1M different 3x3 matrices:"
---     putStrLn "from [[]]"
---     time $ print $ fun (10^6) rot'
---     putStrLn "from []"
+--    putStrLn "from [[]]"
+--    time $ print $ manymult (10^6) rot'
+--    putStrLn "from (3><3) []"
     time $ print $ manymult (10^6) rot
     print $ cos (10^6/2)
 
@@ -81,18 +88,18 @@
     putStrLn "foldVector"
     let v = flatten $ ident 500 :: Vector Double
     print $ vectorMax v  -- evaluate it
-    let getPos k s = if k `mod` 500 < 200 && v@>k > 0 then k:s else s
-    putStrLn "indices extraction, dim=0.25M:"
-    time $ print $ (`divMod` 500) $ maximum $ foldLoop getPos [] (dim v)
-    putStrLn "sum, dim=30M:"
-    --time $ print $ foldLoop (\k s -> w@>k + s) 0.0 (dim w)
-    time $ print $ foldVector (\k v s -> v k + s) 0.0 w
+
+    putStrLn "sum, dim=5M:"
+    -- time $ print $ foldLoop (\k s -> w@>k + s) 0.0 (dim w)
+    time $ print $ sumVector w
+
     putStrLn "sum, dim=0.25M:"
-    --time $ print $ foldVector (\k v s -> v k + s) 0.0 v
-    time $ print $ foldLoop (\k s -> v@>k + s) 0.0 (dim v)
+    --time $ print $ foldLoop (\k s -> v@>k + s) 0.0 (dim v)
+    time $ print $ sumVector v
 
--- foldVector is slower if it is used in two places. (!?)
--- this does not happen with foldLoop
+    let getPos k s = if k `mod` 500 < 200 && v@>k > 0 then k:s else s
+    putStrLn "foldLoop for element selection, dim=0.25M:"
+    time $ print $ (`divMod` 500) $ maximum $ foldLoop getPos [] (dim v)
 
 foldLoop f s d = go (d - 1) s
      where
@@ -101,3 +108,38 @@
 
 foldVector f s v = foldLoop g s (dim v)
     where g !k !s = f k (v@>) s
+          {-# INLINE g #-} -- Thanks Ryan Ingram (http://permalink.gmane.org/gmane.comp.lang.haskell.cafe/46479)
+
+sumVector = foldVector (\k v s -> v k + s) 0.0
+
+-- foldVector is slower if used in two places unless we use the above INLINE
+-- this does not happen with foldLoop
+--------------------------------------------------------------------------------
+
+bench4 = do
+    putStrLn "-------------------------------------------------------"
+    putStrLn "1000x1000 inverse"
+    let a = ident 1000 :: Matrix Double
+    let b = 2*a
+    print $ vectorMax $ flatten (a+b) -- evaluate it
+    time $ print $ vectorMax $ flatten $ linearSolve a b
+
+--------------------------------------------------------------------------------
+
+op1 a b = a <> trans b
+
+op2 a b = a + trans b
+
+timep = time . print . vectorMax . flatten
+
+bench5 n d = do
+    putStrLn "-------------------------------------------------------"
+    putStrLn "transpose in multiply"
+    let ms = replicate n ((ident d :: Matrix Double))
+    let mz = replicate n (diag (constant (0::Double) d))
+    timep $ foldl1' (<>) ms
+    timep $ foldl1' op1  ms
+    putStrLn "-------------------------------------------------------"
+    putStrLn "transpose in add"
+    timep $ foldl1' (+)  ms
+    timep $ foldl1' op2  ms
diff --git a/examples/tests.hs b/examples/tests.hs
new file mode 100644
--- /dev/null
+++ b/examples/tests.hs
@@ -0,0 +1,3 @@
+import Numeric.LinearAlgebra.Tests
+
+main = runTests 20
diff --git a/hmatrix.cabal b/hmatrix.cabal
--- a/hmatrix.cabal
+++ b/hmatrix.cabal
@@ -1,5 +1,5 @@
 Name:               hmatrix
-Version:            0.4.0.0
+Version:            0.5.0.1
 License:            GPL
 License-file:       LICENSE
 Author:             Alberto Ruiz
@@ -7,13 +7,11 @@
 Stability:          provisional
 Homepage:           http://www.hmatrix.googlepages.com
 Synopsis:           Linear algebra and numerical computations
-Description:        A purely functional interface to basic linear algebra computations
-                    and other numerical routines, internally implemented using
+Description:        This library provides a purely functional interface to basic linear algebra
+                    and other numerical computations, internally implemented using
                     GSL, BLAS and LAPACK.
-                    .
-                    More information: <http://www.hmatrix.googlepages.com>
-Category:           Numerical, Math
-tested-with:        GHC ==6.8.3
+Category:           Math
+tested-with:        GHC ==6.10.0
 
 cabal-version:      >=1.2
 build-type:         Simple
@@ -25,24 +23,17 @@
     description:    Link with Intel's MKL optimized libraries.
     default:        False
 
-flag gsl
-    description:    Link with GSL unoptimized blas.
-    default:        False
-
 flag unsafe
     description: Compile the library with bound checking disabled.
     default: False
 
+
 library
     if flag(splitBase)
       build-depends:    base >= 3, array, QuickCheck, HUnit, storable-complex
     else
       build-depends:    base < 3, QuickCheck, HUnit, storable-complex
 
-    if !flag(mkl)
-      if !arch(x86_64)
-        ghc-options:    -fvia-C
-
     Build-Depends:      haskell98
     Extensions:         ForeignFunctionInterface,
                         CPP
@@ -100,26 +91,45 @@
                         Data.Packed.Internal.Vector,
                         Data.Packed.Internal.Matrix,
                         Numeric.GSL.Special.Internal,
-                        Numeric.GSL.Matrix,
                         Numeric.LinearAlgebra.Tests.Instances,
                         Numeric.LinearAlgebra.Tests.Properties
     C-sources:          lib/Data/Packed/Internal/auxi.c,
                         lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c,
                         lib/Numeric/GSL/gsl-aux.c
+
+    ghc-prof-options:   -auto-all
+
+    ghc-options:        -Wall -fno-warn-missing-signatures -fno-warn-orphans -fno-warn-unused-binds
+
+    if flag(unsafe)
+        cpp-options: -DUNSAFE
+
+
     if flag(mkl)
       if arch(x86_64)
         extra-libraries:   gsl mkl_lapack mkl_intel_lp64 mkl_sequential mkl_core
       else
         extra-libraries:   gsl mkl_lapack mkl_intel mkl_sequential mkl_core
     else
-      if flag(gsl)
-        extra-libraries:  gsl gslcblas lapack
-      else
-        extra-libraries:  gsl blas lapack
 
-    cc-options:         -O4
-    ghc-prof-options:   -auto-all
+        extra-libraries: gsl lapack
 
-    if flag(unsafe)
-        cpp-options: -DUNSAFE
+        -- Include additional libraries if they are
+        -- required by your system to link -lgsl -llapack
 
+        -- (In ubuntu/debian cblas is included in blas,
+        --  which is automatically linked by lapack, so
+        -- nothing more is required.)
+
+        -- Examples:
+
+            -------- if blas/cblas are not automatically linked by lapack:
+            -- blas cblas
+
+            -------- Nonoptimized cblas included in gsl:
+            -- gslcblas
+
+            -------- Arch Linux with atlas-lapack:
+            -- f77blas cblas atlas gcc_s
+            -------- Arch Linux with normal blas and lapack:
+            -- blas gslcblas gfortran
diff --git a/lib/Data/Packed/Convert.hs b/lib/Data/Packed/Convert.hs
--- a/lib/Data/Packed/Convert.hs
+++ b/lib/Data/Packed/Convert.hs
@@ -29,7 +29,6 @@
 import Foreign
 import Control.Monad.ST
 import Data.Array.ST
-import Data.Array.IArray
 import Data.Array.Unboxed
 
 -- | Creates a StorableArray indexed from 0 to dim -1.
@@ -88,7 +87,7 @@
 matrixFromArray :: UArray (Int, Int) Double -> Matrix Double
 matrixFromArray m = reshape c . fromList . elems $ m
     where ((r1,c1),(r2,c2)) = bounds m
-          r = r2-r1+1
+          _r = r2-r1+1
           c = c2-c1+1
 
 arrayFromMatrix :: Matrix Double -> UArray (Int, Int) Double
diff --git a/lib/Data/Packed/Internal/Common.hs b/lib/Data/Packed/Internal/Common.hs
--- a/lib/Data/Packed/Internal/Common.hs
+++ b/lib/Data/Packed/Internal/Common.hs
@@ -22,7 +22,7 @@
 import Debug.Trace
 import Foreign.C.String(peekCString)
 import Foreign.C.Types
-import Foreign.Storable.Complex
+import Foreign.Storable.Complex()
 
 
 -- | @debug x = trace (show x) x@
@@ -80,9 +80,14 @@
 errorCode 2007 = "not yet supported in this OS"
 errorCode n    = "code "++show n
 
+
+-- | clear the fpu
+foreign import ccall "auxi.h asm_finit" finit :: IO ()
+
 -- | check the error code
 check :: String -> IO CInt -> IO ()
 check msg f = do
+    finit
     err <- f
     when (err/=0) $ if err > 1024
                       then (error (msg++": "++errorCode err)) -- our errors
diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs
--- a/lib/Data/Packed/Internal/Matrix.hs
+++ b/lib/Data/Packed/Internal/Matrix.hs
@@ -22,10 +22,8 @@
 
 import Foreign hiding (xor)
 import Complex
-import Control.Monad(when)
 import Foreign.C.String
 import Foreign.C.Types
-import Data.List(transpose)
 
 -----------------------------------------------------------------
 
@@ -212,7 +210,6 @@
 class (Storable a, Floating a) => Element a where
     constantD :: a -> Int -> Vector a
     transdata :: Int -> Vector a -> Int -> Vector a
-    multiplyD :: Matrix a -> Matrix a -> Matrix a
     subMatrixD :: (Int,Int) -- ^ (r0,c0) starting position 
                -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix
                -> Matrix a -> Matrix a
@@ -221,14 +218,12 @@
 instance Element Double where
     constantD  = constantR
     transdata = transdataR
-    multiplyD  = multiplyR
     subMatrixD = subMatrixR
     diagD      = diagR
 
 instance Element (Complex Double) where
     constantD  = constantC
     transdata  = transdataC
-    multiplyD  = multiplyC
     subMatrixD = subMatrixC
     diagD      = diagC
 
@@ -266,33 +261,6 @@
 foreign import ccall "auxi.h transR" ctransR :: TMM
 foreign import ccall "auxi.h transC" ctransC :: TCMCM
 
-------------------------------------------------------------------
-
-gmatC MF { rows = r, cols = c } p f = f 1 (fi c) (fi r) p
-gmatC MC { rows = r, cols = c } p f = f 0 (fi r) (fi c) p
-
-dtt MC { cdat = d } = d
-dtt MF { fdat = d } = d
-
-multiplyAux fun a b = unsafePerformIO $ do
-    when (cols a /= rows b) $ error $ "inconsistent dimensions in contraction "++
-                                      show (rows a,cols a) ++ " x " ++ show (rows b, cols b)
-    r <- createMatrix RowMajor (rows a) (cols b)
-    withForeignPtr (fptr (dtt a)) $ \pa -> withForeignPtr (fptr (dtt b)) $ \pb ->
-        withMatrix r $ \r' ->
-            fun // gmatC a pa // gmatC b pb // r' // check "multiplyAux"
-    return r
-
-multiplyR = multiplyAux cmultiplyR
-foreign import ccall "auxi.h multiplyR" cmultiplyR :: TauxMul Double
-
-multiplyC = multiplyAux cmultiplyC
-foreign import ccall "auxi.h multiplyC" cmultiplyC :: TauxMul (Complex Double)
-
--- | matrix product
-multiply :: (Element a) => Matrix a -> Matrix a -> Matrix a
-multiply = multiplyD
-
 ----------------------------------------------------------------------
 
 -- | extraction of a submatrix from a real matrix
@@ -370,7 +338,12 @@
 
 -- | obtains the complex conjugate of a complex vector
 conj :: Vector (Complex Double) -> Vector (Complex Double)
-conj v = asComplex $ flatten $ reshape 2 (asReal v) `multiply` diag (fromList [1,-1])
+conj v = unsafePerformIO $ do
+    r <- createVector (dim v)
+    app2 cconjugate vec v vec r "cconjugate"
+    return r
+foreign import ccall "auxi.h conjugate" cconjugate :: TCVCV
+
 
 -- | creates a complex vector from vectors with real and imaginary parts
 toComplex :: (Vector Double, Vector Double) ->  Vector (Complex Double)
diff --git a/lib/Data/Packed/Internal/auxi.c b/lib/Data/Packed/Internal/auxi.c
--- a/lib/Data/Packed/Internal/auxi.c
+++ b/lib/Data/Packed/Internal/auxi.c
@@ -4,14 +4,9 @@
 #include <gsl/gsl_matrix.h>
 #include <gsl/gsl_math.h>
 #include <gsl/gsl_errno.h>
-#include <gsl/gsl_fft_complex.h>
-#include <gsl/gsl_eigen.h>
-#include <gsl/gsl_integration.h>
-#include <gsl/gsl_deriv.h>
-#include <gsl/gsl_poly.h>
-#include <gsl/gsl_multimin.h>
 #include <gsl/gsl_complex.h>
 #include <gsl/gsl_complex_math.h>
+#include <gsl/gsl_cblas.h>
 #include <string.h>
 #include <stdio.h>
 
@@ -118,78 +113,6 @@
 }
 
 
-int multiplyR(int ta, KRMAT(a), int tb, KRMAT(b),RMAT(r)) {
-    //printf("%d %d %d %d %d %d\n",ar,ac,br,bc,rr,rc);
-    //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);
-    DEBUGMSG("multiplyR (gsl_blas_dgemm)");
-    KDMVIEW(a);
-    KDMVIEW(b);
-    DMVIEW(r);
-    int k;
-    for(k=0;k<rr*rc;k++) rp[k]=0;
-    int debug = 0;
-    if(debug) {
-        printf("---------------------------\n");
-        printf("%p: ",ap); for(k=0;k<ar*ac;k++) printf("%f ",ap[k]); printf("\n");
-        printf("%p: ",bp); for(k=0;k<br*bc;k++) printf("%f ",bp[k]); printf("\n");
-        printf("%p: ",rp); for(k=0;k<rr*rc;k++) printf("%f ",rp[k]); printf("\n");
-    }
-    int res = gsl_blas_dgemm(
-         ta?CblasTrans:CblasNoTrans,
-         tb?CblasTrans:CblasNoTrans,
-         1.0, M(a), M(b),
-         0.0, M(r));
-    if(debug) {
-        printf("--------------\n");
-        printf("%p: ",ap); for(k=0;k<ar*ac;k++) printf("%f ",ap[k]); printf("\n");
-        printf("%p: ",bp); for(k=0;k<br*bc;k++) printf("%f ",bp[k]); printf("\n");
-        printf("%p: ",rp); for(k=0;k<rr*rc;k++) printf("%f ",rp[k]); printf("\n");
-    }
-    CHECK(res,res);
-    OK
-}
-
-int multiplyC(int ta, KCMAT(a), int tb, KCMAT(b),CMAT(r)) {
-    //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);
-    DEBUGMSG("multiplyC (gsl_blas_zgemm)");
-    KCMVIEW(a);
-    KCMVIEW(b);
-    CMVIEW(r);
-    int k;
-    gsl_complex alpha, beta;
-    GSL_SET_COMPLEX(&alpha,1.,0.);
-    GSL_SET_COMPLEX(&beta,0.,0.);
-    //double *TEMP = (double*)malloc(rr*rc*2*sizeof(double));
-    //gsl_matrix_complex_view T = gsl_matrix_complex_view_array(TEMP,rr,rc);
-    for(k=0;k<rr*rc;k++) rp[k]=beta;
-    //for(k=0;k<2*rr*rc;k++) TEMP[k]=0;
-    int debug = 0;
-    if(debug) {
-        printf("---------------------------\n");
-        printf("%p: ",ap); for(k=0;k<2*ar*ac;k++) printf("%f ",((double*)ap)[k]); printf("\n");
-        printf("%p: ",bp); for(k=0;k<2*br*bc;k++) printf("%f ",((double*)bp)[k]); printf("\n");
-        printf("%p: ",rp); for(k=0;k<2*rr*rc;k++) printf("%f ",((double*)rp)[k]); printf("\n");
-        //printf("%p: ",T); for(k=0;k<2*rr*rc;k++) printf("%f ",TEMP[k]); printf("\n");
-    }
-    int res = gsl_blas_zgemm(
-         ta?CblasTrans:CblasNoTrans,
-         tb?CblasTrans:CblasNoTrans,
-         alpha, M(a), M(b),
-         beta, M(r)); 
-         //&T.matrix);
-    //memcpy(rp,TEMP,2*rr*rc*sizeof(double));
-    if(debug) {
-        printf("--------------\n");
-        printf("%p: ",ap); for(k=0;k<2*ar*ac;k++) printf("%f ",((double*)ap)[k]); printf("\n");
-        printf("%p: ",bp); for(k=0;k<2*br*bc;k++) printf("%f ",((double*)bp)[k]); printf("\n");
-        printf("%p: ",rp); for(k=0;k<2*rr*rc;k++) printf("%f ",((double*)rp)[k]); printf("\n");
-        //printf("%p: ",T); for(k=0;k<2*rr*rc;k++) printf("%f ",TEMP[k]); printf("\n");
-    }
-    CHECK(res,res);
-    OK
-}
-
-
 int diagR(KRVEC(d),RMAT(r)) {
     REQUIRES(dn==rr && rr==rc,BAD_SIZE);
     DEBUGMSG("diagR");
@@ -215,3 +138,20 @@
     }
     OK
 }
+
+int conjugate(KCVEC(x),CVEC(t)) {
+    REQUIRES(xn==tn,BAD_SIZE);
+    DEBUGMSG("conjugate");
+    int k;
+    for (k=0; k<xn; k++) {
+        tp[k].dat[0] =   xp[k].dat[0];
+        tp[k].dat[1] = - xp[k].dat[1];
+    }
+    OK
+}
+
+//---------------------------------------
+void asm_finit() {
+    asm("finit");
+}
+//---------------------------------------
diff --git a/lib/Data/Packed/Internal/auxi.h b/lib/Data/Packed/Internal/auxi.h
--- a/lib/Data/Packed/Internal/auxi.h
+++ b/lib/Data/Packed/Internal/auxi.h
@@ -10,16 +10,12 @@
 #define KCVEC(A) int A##n, const gsl_complex*A##p
 #define KCMAT(A) int A##r, int A##c, const gsl_complex* A##p
 
-
 int transR(KRMAT(x),RMAT(t));
 int transC(KCMAT(x),CMAT(t));
 
 int constantR(double *val     , RVEC(r));
 int constantC(gsl_complex *val, CVEC(r));
 
-int multiplyR(int ta, KRMAT(a), int tb,  KRMAT(b),RMAT(r));
-int multiplyC(int ta, KCMAT(a), int tb, KCMAT(b),CMAT(r));
-
 int submatrixR(int r1, int r2, int c1, int c2, KRMAT(x),RMAT(r));
 
 int diagR(KRVEC(d),RMAT(r));
@@ -28,3 +24,7 @@
 const char * gsl_strerror (const int gsl_errno);
 
 int matrix_fscanf(char*filename, RMAT(a));
+
+int conjugate(KCVEC(x),CVEC(t));
+
+void asm_finit();
diff --git a/lib/Data/Packed/ST.hs b/lib/Data/Packed/ST.hs
--- a/lib/Data/Packed/ST.hs
+++ b/lib/Data/Packed/ST.hs
@@ -30,9 +30,7 @@
 ) where
 
 import Data.Packed.Internal
-import Data.Array.Storable
 import Control.Monad.ST
-import Data.Array.ST
 import Foreign
 
 {-# INLINE ioReadV #-}
@@ -97,13 +95,13 @@
 
 {-# INLINE ioReadM #-}
 ioReadM :: Storable t => Matrix t -> Int -> Int -> IO t
-ioReadM (MC nr nc cv) r c = ioReadV cv (r*nc+c)
-ioReadM (MF nr nc fv) r c = ioReadV fv (c*nr+r)
+ioReadM (MC _ nc cv) r c = ioReadV cv (r*nc+c)
+ioReadM (MF nr _ fv) r c = ioReadV fv (c*nr+r)
 
 {-# INLINE ioWriteM #-}
 ioWriteM :: Storable t => Matrix t -> Int -> Int -> t -> IO ()
-ioWriteM (MC nr nc cv) r c val = ioWriteV cv (r*nc+c) val
-ioWriteM (MF nr nc fv) r c val = ioWriteV fv (c*nr+r) val
+ioWriteM (MC _ nc cv) r c val = ioWriteV cv (r*nc+c) val
+ioWriteM (MF nr _ fv) r c val = ioWriteV fv (c*nr+r) val
 
 newtype STMatrix s t = STMatrix (Matrix t)
 
diff --git a/lib/Graphics/Plot.hs b/lib/Graphics/Plot.hs
--- a/lib/Graphics/Plot.hs
+++ b/lib/Graphics/Plot.hs
@@ -9,7 +9,9 @@
 -- Portability :  uses gnuplot and ImageMagick
 --
 -- Very basic (and provisional) drawing tools using gnuplot and imageMagick.
---
+-- 
+-- This module is deprecated. It will be replaced by improved drawing tools based
+-- on the Gnuplot package by Henning Thielemann.
 -----------------------------------------------------------------------------
 
 module Graphics.Plot(
@@ -28,7 +30,7 @@
 
 import Data.Packed.Vector
 import Data.Packed.Matrix
-import Numeric.LinearAlgebra.Linear(outer)
+import Numeric.LinearAlgebra(outer)
 import Numeric.GSL.Vector(FunCodeS(Max,Min),toScalarR)
 import Data.List(intersperse)
 import System
diff --git a/lib/Numeric/GSL/Matrix.hs b/lib/Numeric/GSL/Matrix.hs
deleted file mode 100644
--- a/lib/Numeric/GSL/Matrix.hs
+++ /dev/null
@@ -1,311 +0,0 @@
------------------------------------------------------------------------------
--- |
--- Module      :  Numeric.GSL.Matrix
--- Copyright   :  (c) Alberto Ruiz 2007
--- License     :  GPL-style
---
--- Maintainer  :  Alberto Ruiz <aruiz@um.es>
--- Stability   :  provisional
--- Portability :  portable (uses FFI)
---
--- A few linear algebra computations based on GSL.
---
------------------------------------------------------------------------------
--- #hide
-
-module Numeric.GSL.Matrix(
-    eigSg, eigHg,
-    svdg,
-    qr, qrPacked, unpackQR,
-    cholR, cholC,
-    luSolveR, luSolveC,
-    luR, luC
-) where
-
-import Data.Packed.Internal
-import Data.Packed.Matrix(ident)
-import Numeric.GSL.Vector
-import Foreign
-import Complex
-
-{- | eigendecomposition of a real symmetric matrix using /gsl_eigen_symmv/.
-
-> > let (l,v) = eigS $ 'fromLists' [[1,2],[2,1]]
-> > l
-> 3.000 -1.000
->
-> > v
-> 0.707 -0.707
-> 0.707  0.707
->
-> > v <> diag l <> trans v
-> 1.000 2.000
-> 2.000 1.000
-
--}
-eigSg :: Matrix Double -> (Vector Double, Matrix Double)
-eigSg = eigSg' . cmat
-
-eigSg' m
-    | r == 1 = (fromList [cdat m `at` 0], singleton 1)
-    | otherwise = unsafePerformIO $ do
-        l <- createVector r
-        v <- createMatrix RowMajor r r
-        app3 c_eigS mat m vec l mat v "eigSg"
-        return (l,v)
-  where r = rows m
-foreign import ccall "gsl-aux.h eigensystemR" c_eigS :: TMVM
-
-------------------------------------------------------------------
-
-
-
-{- | eigendecomposition of a complex hermitian matrix using /gsl_eigen_hermv/
-
-> > let (l,v) = eigH $ 'fromLists' [[1,2+i],[2-i,3]]
->
-> > l
-> 4.449 -0.449
->
-> > v
->         -0.544          0.839
-> (-0.751,0.375) (-0.487,0.243)
->
-> > v <> diag l <> (conjTrans) v
->          1.000 (2.000,1.000)
-> (2.000,-1.000)         3.000
-
--}
-eigHg :: Matrix (Complex Double)-> (Vector Double, Matrix (Complex Double))
-eigHg = eigHg' . cmat
-
-eigHg' m
-    | r == 1 = (fromList [realPart $ cdat m `at` 0], singleton 1)
-    | otherwise = unsafePerformIO $ do
-        l <- createVector r
-        v <- createMatrix RowMajor r r
-        app3 c_eigH mat m vec l mat v "eigHg"
-        return (l,v)
-  where r = rows m
-foreign import ccall "gsl-aux.h eigensystemC" c_eigH :: TCMVCM
-
-
-{- | Singular value decomposition of a real matrix, using /gsl_linalg_SV_decomp_mod/:
-
-
--}
-svdg :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
-svdg x = if rows x >= cols x
-    then svd' (cmat x)
-    else (v, s, u) where (u,s,v) = svd' (cmat (trans x))
-
-svd' x = unsafePerformIO $ do
-    u <- createMatrix RowMajor r c
-    s <- createVector c
-    v <- createMatrix RowMajor c c
-    app4 c_svd mat x mat u vec s mat v "svdg"
-    return (u,s,v)
-  where r = rows x
-        c = cols x
-foreign import ccall "gsl-aux.h svd" c_svd :: TMMVM
-
-{- | QR decomposition of a real matrix using /gsl_linalg_QR_decomp/ and /gsl_linalg_QR_unpack/.
-
--}
-qr :: Matrix Double -> (Matrix Double, Matrix Double)
-qr = qr' . cmat
-
-qr' x = unsafePerformIO $ do
-    q <- createMatrix RowMajor r r
-    rot <- createMatrix RowMajor r c
-    app3 c_qr mat x mat q mat rot  "qr"
-    return (q,rot)
-  where r = rows x
-        c = cols x
-foreign import ccall "gsl-aux.h QR" c_qr :: TMMM
-
-qrPacked :: Matrix Double -> (Matrix Double, Vector Double)
-qrPacked = qrPacked' . cmat
-
-qrPacked' x = unsafePerformIO $ do
-    qrp <- createMatrix RowMajor r c
-    tau <- createVector (min r c)
-    app3 c_qrPacked mat x mat qrp vec tau "qrUnpacked"
-    return (qrp,tau)
-  where r = rows x
-        c = cols x
-foreign import ccall "gsl-aux.h QRpacked" c_qrPacked :: TMMV
-
-unpackQR :: (Matrix Double, Vector Double) -> (Matrix Double, Matrix Double)
-unpackQR (qrp,tau) = unpackQR' (cmat qrp, tau)
-
-unpackQR' (qrp,tau) = unsafePerformIO $ do
-    q <- createMatrix RowMajor r r
-    res <- createMatrix RowMajor r c
-    app4 c_qrUnpack mat qrp vec tau mat q mat res "qrUnpack"
-    return (q,res)
-  where r = rows qrp
-        c = cols qrp
-foreign import ccall "gsl-aux.h QRunpack" c_qrUnpack :: TMVMM
-
-{- | Cholesky decomposition of a symmetric positive definite real matrix using /gsl_linalg_cholesky_decomp/.
-
-@\> chol $ (2><2) [1,2,
-                   2,9::Double]
-(2><2)
- [ 1.0,              0.0
- , 2.0, 2.23606797749979 ]@
-
--}
-cholR :: Matrix Double -> Matrix Double
-cholR = cholR' . cmat
-
-cholR' x = unsafePerformIO $ do
-    r <- createMatrix RowMajor n n
-    app2 c_cholR mat x mat r "cholR"
-    return r
-  where n = rows x
-foreign import ccall "gsl-aux.h cholR" c_cholR :: TMM
-
-cholC :: Matrix (Complex Double) -> Matrix (Complex Double)
-cholC = cholC' . cmat
-
-cholC' x = unsafePerformIO $ do
-    r <- createMatrix RowMajor n n
-    app2 c_cholC mat x mat r "cholC"
-    return r
-  where n = rows x
-foreign import ccall "gsl-aux.h cholC" c_cholC :: TCMCM
-
-
---------------------------------------------------------
-
-{- -| efficient multiplication by the inverse of a matrix (for real matrices)
--}
-luSolveR :: Matrix Double -> Matrix Double -> Matrix Double
-luSolveR a b = luSolveR' (cmat a) (cmat b)
-
-luSolveR' a b
-    | n1==n2 && n1==r = unsafePerformIO $ do
-        s <- createMatrix RowMajor r c
-        app3 c_luSolveR mat a mat b mat s "luSolveR"
-        return s
-    | otherwise = error "luSolveR of nonsquare matrix"
-  where n1 = rows a
-        n2 = cols a
-        r  = rows b
-        c  = cols b
-foreign import ccall "gsl-aux.h luSolveR" c_luSolveR ::  TMMM
-
-{- -| efficient multiplication by the inverse of a matrix (for complex matrices). 
--}
-luSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
-luSolveC a b = luSolveC' (cmat a) (cmat b)
-
-luSolveC' a b
-    | n1==n2 && n1==r = unsafePerformIO $ do
-        s <- createMatrix RowMajor r c
-        app3 c_luSolveC mat a mat b mat s "luSolveC"
-        return s
-    | otherwise = error "luSolveC of nonsquare matrix"
-  where n1 = rows a
-        n2 = cols a
-        r  = rows b
-        c  = cols b
-foreign import ccall "gsl-aux.h luSolveC" c_luSolveC ::  TCMCMCM
-
-{- | lu decomposition of real matrix (packed as a vector including l, u, the permutation and sign)
--}
-luRaux  :: Matrix Double -> Vector Double
-luRaux = luRaux' . cmat
-
-luRaux' x = unsafePerformIO $ do
-    res <- createVector (r*r+r+1)
-    app2 c_luRaux mat x vec res "luRaux"
-    return res
-  where r = rows x
-foreign import ccall "gsl-aux.h luRaux" c_luRaux :: TMV
-
-{- | lu decomposition of complex matrix (packed as a vector including l, u, the permutation and sign)
--}
-luCaux  :: Matrix (Complex Double) -> Vector (Complex Double)
-luCaux = luCaux' . cmat
-
-luCaux' x = unsafePerformIO $ do
-    res <- createVector (r*r+r+1)
-    app2 c_luCaux mat x vec res "luCaux"
-    return res
-  where r = rows x
-foreign import ccall "gsl-aux.h luCaux" c_luCaux :: TCMCV
-
-{- | The LU decomposition of a square matrix. Is based on /gsl_linalg_LU_decomp/ and  /gsl_linalg_complex_LU_decomp/ as described in <http://www.gnu.org/software/Numeric.GSL/manual/Numeric.GSL-ref_13.html#SEC223>.
-
-@\> let m = 'fromLists' [[1,2,-3],[2+3*i,-7,0],[1,-i,2*i]]
-\> let (l,u,p,s) = luR m@
-
-L is the lower triangular:
-
-@\> l
-          1.            0.  0.
-0.154-0.231i            1.  0.
-0.154-0.231i  0.624-0.522i  1.@
-
-U is the upper triangular:
-
-@\> u
-2.+3.i           -7.            0.
-    0.  3.077-1.615i           -3.
-    0.            0.  1.873+0.433i@
-
-p is a permutation:
-
-@\> p
-[1,0,2]@
-
-L \* U obtains a permuted version of the original matrix:
-
-@\> extractRows p m
-  2.+3.i   -7.   0.
-      1.    2.  -3.
-      1.  -1.i  2.i
-\ -- CPP
-\> l \<\> u
- 2.+3.i   -7.   0.
-     1.    2.  -3.
-     1.  -1.i  2.i@
-
-s is the sign of the permutation, required to obtain sign of the determinant:
-
-@\> s * product ('toList' $ 'takeDiag' u)
-(-18.0) :+ (-16.000000000000004)
-\> 'LinearAlgebra.Algorithms.det' m
-(-18.0) :+ (-16.000000000000004)@
-
- -}
-luR :: Matrix Double -> (Matrix Double, Matrix Double, [Int], Double)
-luR m = (l,u,p, fromIntegral s') where
-    r = rows m
-    v = luRaux m
-    lu = reshape r $ subVector 0 (r*r) v
-    s':p = map round . toList . subVector (r*r) (r+1) $ v
-    u = triang r r 0 1`mul` lu
-    l = (triang r r 0 0 `mul` lu) `add` ident r
-    add = liftMatrix2 $ vectorZipR Add
-    mul = liftMatrix2 $ vectorZipR Mul
-
--- | Complex version of 'luR'.
-luC :: Matrix (Complex Double) -> (Matrix (Complex Double), Matrix (Complex Double), [Int], Complex Double)
-luC m = (l,u,p, fromIntegral s') where
-    r = rows m
-    v = luCaux m
-    lu = reshape r $ subVector 0 (r*r) v
-    s':p = map (round.realPart) . toList . subVector (r*r) (r+1) $ v
-    u = triang r r 0 1 `mul` lu
-    l = (triang r r 0 0 `mul` lu) `add` liftMatrix comp (ident r)
-    add = liftMatrix2 $ vectorZipC Add
-    mul = liftMatrix2 $ vectorZipC Mul
-
-{- auxiliary function to get triangular matrices
--}
-triang r c h v = reshape c $ fromList [el i j | i<-[0..r-1], j<-[0..c-1]]
-    where el i j = if j-i>=h then v else 1 - v
diff --git a/lib/Numeric/GSL/Special/Internal.hsc b/lib/Numeric/GSL/Special/Internal.hsc
--- a/lib/Numeric/GSL/Special/Internal.hsc
+++ b/lib/Numeric/GSL/Special/Internal.hsc
@@ -1,4 +1,4 @@
-{-# OPTIONS -ffi #-}
+ {-# LANGUAGE ForeignFunctionInterface #-}
 -----------------------------------------------------------------------------
 {- |
 Module      :  Numeric.GSL.Special.Internal
diff --git a/lib/Numeric/GSL/Vector.hs b/lib/Numeric/GSL/Vector.hs
--- a/lib/Numeric/GSL/Vector.hs
+++ b/lib/Numeric/GSL/Vector.hs
@@ -27,7 +27,7 @@
 import Foreign
 import Foreign.C.Types(CInt)
 
-fromei x = fromIntegral (fromEnum x)
+fromei x = fromIntegral (fromEnum x) :: CInt
 
 data FunCodeV = Sin
               | Cos
diff --git a/lib/Numeric/GSL/gsl-aux.c b/lib/Numeric/GSL/gsl-aux.c
--- a/lib/Numeric/GSL/gsl-aux.c
+++ b/lib/Numeric/GSL/gsl-aux.c
@@ -1,11 +1,8 @@
 #include "gsl-aux.h"
 #include <gsl/gsl_blas.h>
-#include <gsl/gsl_linalg.h>
-#include <gsl/gsl_matrix.h>
 #include <gsl/gsl_math.h>
 #include <gsl/gsl_errno.h>
 #include <gsl/gsl_fft_complex.h>
-#include <gsl/gsl_eigen.h>
 #include <gsl/gsl_integration.h>
 #include <gsl/gsl_deriv.h>
 #include <gsl/gsl_poly.h>
@@ -98,7 +95,27 @@
     }
 }
 
+inline gsl_complex complex_abs(gsl_complex z) {
+    gsl_complex r;
+    r.dat[0] = gsl_complex_abs(z);
+    r.dat[1] = 0;
+    return r;
+}
 
+inline gsl_complex complex_signum(gsl_complex z) {
+    gsl_complex r;
+    double mag;
+    if (z.dat[0] == 0 && z.dat[1] == 0) {
+        r.dat[0] = 0;
+        r.dat[1] = 0;
+    } else {
+        mag = gsl_complex_abs(z);
+        r.dat[0] = z.dat[0]/mag;
+        r.dat[1] = z.dat[1]/mag;
+    }
+    return r;
+}
+
 #define OP(C,F) case C: { for(k=0;k<xn;k++) rp[k] = F(xp[k]); OK }
 #define OPV(C,E) case C: { for(k=0;k<xn;k++) rp[k] = E; OK }
 int mapR(int code, KRVEC(x), RVEC(r)) {
@@ -127,6 +144,7 @@
     }
 }
 
+
 int mapCAux(int code, KGCVEC(x), GCVEC(r)) {
     int k;
     REQUIRES(xn == rn,BAD_SIZE);
@@ -135,7 +153,7 @@
         OP(0,gsl_complex_sin)
         OP(1,gsl_complex_cos)
         OP(2,gsl_complex_tan)
-
+        OP(3,complex_abs)
         OP(4,gsl_complex_arcsin)
         OP(5,gsl_complex_arccos)
         OP(6,gsl_complex_arctan)
@@ -147,7 +165,7 @@
         OP(12,gsl_complex_arctanh)
         OP(13,gsl_complex_exp)
         OP(14,gsl_complex_log)
-
+        OP(15,complex_signum)
         OP(16,gsl_complex_sqrt)
 
         // gsl_complex_arg
@@ -161,47 +179,6 @@
 }
 
 
-/*
-int scaleR(double* alpha, KRVEC(x), RVEC(r)) {
-    REQUIRES(xn == rn,BAD_SIZE);
-    DEBUGMSG("scaleR");
-    KDVVIEW(x);
-    DVVIEW(r);
-    CHECK( gsl_vector_memcpy(V(r),V(x)) , MEM);
-    int res = gsl_vector_scale(V(r),*alpha);
-    CHECK(res,res);
-    OK
-}
-
-int scaleC(gsl_complex *alpha, KCVEC(x), CVEC(r)) {
-    REQUIRES(xn == rn,BAD_SIZE);
-    DEBUGMSG("scaleC");
-    //KCVVIEW(x);
-    CVVIEW(r);
-    gsl_vector_const_view vrx = gsl_vector_const_view_array((double*)xp,xn*2);
-    gsl_vector_view vrr = gsl_vector_view_array((double*)rp,rn*2);
-    CHECK(gsl_vector_memcpy(V(vrr),V(vrx)) , MEM);
-    gsl_blas_zscal(*alpha,V(r)); // void !
-    int res = 0; 
-    CHECK(res,res);
-    OK
-}
-
-int addConstantR(double offs, KRVEC(x), RVEC(r)) { 
-    REQUIRES(xn == rn,BAD_SIZE); 
-    DEBUGMSG("addConstantR");
-    KDVVIEW(x);
-    DVVIEW(r);
-    CHECK(gsl_vector_memcpy(V(r),V(x)), MEM);
-    int res = gsl_vector_add_constant(V(r),offs);
-    CHECK(res,res);
-    OK
-}
-
-*/
-
-
-
 int mapValR(int code, double* pval, KRVEC(x), RVEC(r)) {
     int k;
     double val = *pval;
@@ -290,248 +267,6 @@
 }
 
 
-
-
-int luSolveR(KRMAT(a),KRMAT(b),RMAT(r)) {
-    REQUIRES(ar==ac && ac==br && ar==rr && bc==rc,BAD_SIZE);
-    DEBUGMSG("luSolveR");
-    KDMVIEW(a);
-    KDMVIEW(b);
-    DMVIEW(r);
-    int res;
-    gsl_matrix *LU = gsl_matrix_alloc(ar,ar);
-    CHECK(!LU,MEM);
-    int s;
-    gsl_permutation * p = gsl_permutation_alloc (ar);
-    CHECK(!p,MEM);
-    CHECK(gsl_matrix_memcpy(LU,M(a)),MEM);
-    res = gsl_linalg_LU_decomp(LU, p, &s);
-    CHECK(res,res);
-    int c;
-
-    for (c=0; c<bc; c++) {
-        gsl_vector_const_view colb = gsl_matrix_const_column (M(b), c);
-        gsl_vector_view colr = gsl_matrix_column (M(r), c);
-        res = gsl_linalg_LU_solve (LU, p, V(colb), V(colr));
-        CHECK(res,res);
-    }
-    gsl_permutation_free(p);
-    gsl_matrix_free(LU);
-    OK
-}
-
-
-int luSolveC(KCMAT(a),KCMAT(b),CMAT(r)) {
-    REQUIRES(ar==ac && ac==br && ar==rr && bc==rc,BAD_SIZE);
-    DEBUGMSG("luSolveC");
-    KCMVIEW(a);
-    KCMVIEW(b);
-    CMVIEW(r);
-    gsl_matrix_complex *LU = gsl_matrix_complex_alloc(ar,ar);
-    CHECK(!LU,MEM);
-    int s;
-    gsl_permutation * p = gsl_permutation_alloc (ar);
-    CHECK(!p,MEM);
-    CHECK(gsl_matrix_complex_memcpy(LU,M(a)),MEM);
-    int res;
-    res = gsl_linalg_complex_LU_decomp(LU, p, &s);
-    CHECK(res,res);
-    int c;
-    for (c=0; c<bc; c++) {
-        gsl_vector_complex_const_view colb = gsl_matrix_complex_const_column (M(b), c);
-        gsl_vector_complex_view colr = gsl_matrix_complex_column (M(r), c);
-        res = gsl_linalg_complex_LU_solve (LU, p, V(colb), V(colr));
-        CHECK(res,res);
-    }
-    gsl_permutation_free(p);
-    gsl_matrix_complex_free(LU);
-    OK
-}
-
-
-int luRaux(KRMAT(a),RVEC(b)) {
-    REQUIRES(ar==ac && bn==ar*ar+ar+1,BAD_SIZE);
-    DEBUGMSG("luRaux");
-    KDMVIEW(a);
-    //DVVIEW(b);
-    gsl_matrix_view LU = gsl_matrix_view_array(bp,ar,ac);
-    int s;
-    gsl_permutation * p = gsl_permutation_alloc (ar);
-    CHECK(!p,MEM);
-    CHECK(gsl_matrix_memcpy(M(LU),M(a)),MEM);
-    gsl_linalg_LU_decomp(M(LU), p, &s);
-    bp[ar*ar] = s;
-    int k;
-    for (k=0; k<ar; k++) {
-        bp[ar*ar+k+1] = gsl_permutation_get(p,k);
-    }
-    gsl_permutation_free(p);
-    OK
-}
-
-int luCaux(KCMAT(a),CVEC(b)) {
-    REQUIRES(ar==ac && bn==ar*ar+ar+1,BAD_SIZE);
-    DEBUGMSG("luCaux");
-    KCMVIEW(a);
-    //DVVIEW(b);
-    gsl_matrix_complex_view LU = gsl_matrix_complex_view_array((double*)bp,ar,ac);
-    int s;
-    gsl_permutation * p = gsl_permutation_alloc (ar);
-    CHECK(!p,MEM);
-    CHECK(gsl_matrix_complex_memcpy(M(LU),M(a)),MEM);
-    int res;
-    res = gsl_linalg_complex_LU_decomp(M(LU), p, &s);
-    CHECK(res,res);
-    ((double*)bp)[2*ar*ar] = s;
-    ((double*)bp)[2*ar*ar+1] = 0;
-    int k;
-    for (k=0; k<ar; k++) {
-        ((double*)bp)[2*ar*ar+2*k+2] = gsl_permutation_get(p,k);
-        ((double*)bp)[2*ar*ar+2*k+2+1] = 0;
-    }
-    gsl_permutation_free(p);
-    OK
-}
-
-int svd(KRMAT(a),RMAT(u), RVEC(s),RMAT(v)) {
-    REQUIRES(ar==ur && ac==uc && ac==sn && ac==vr && ac==vc,BAD_SIZE);
-    DEBUGMSG("svd");
-    KDMVIEW(a);
-    DMVIEW(u);
-    DVVIEW(s);
-    DMVIEW(v);
-    gsl_vector  *workv = gsl_vector_alloc(ac);
-    CHECK(!workv,MEM);
-    gsl_matrix  *workm = gsl_matrix_alloc(ac,ac);
-    CHECK(!workm,MEM);
-    CHECK(gsl_matrix_memcpy(M(u),M(a)),MEM);
-    // int res = gsl_linalg_SV_decomp_jacobi (&U.matrix, &V.matrix, &S.vector);
-    // doesn't work 
-    //int res = gsl_linalg_SV_decomp (&U.matrix, &V.matrix, &S.vector, workv);
-    int res = gsl_linalg_SV_decomp_mod (M(u), workm, M(v), V(s), workv);
-    CHECK(res,res);
-    //gsl_matrix_transpose(M(v));
-    gsl_vector_free(workv);
-    gsl_matrix_free(workm);
-    OK
-}
-
-
-// for real symmetric matrices
-int eigensystemR(KRMAT(x),RVEC(l),RMAT(v)) {
-    REQUIRES(xr==xc && xr==ln && xr==vr && vr==vc,BAD_SIZE);
-    DEBUGMSG("eigensystemR (gsl_eigen_symmv)");
-    KDMVIEW(x);
-    DVVIEW(l);
-    DMVIEW(v);
-    gsl_matrix *XC = gsl_matrix_alloc(xr,xr);
-    gsl_matrix_memcpy(XC,M(x)); // needed because the argument is destroyed
-                                // many thanks to Nico Mahlo for the bug report 
-    gsl_eigen_symmv_workspace * w = gsl_eigen_symmv_alloc (xc);
-    int res = gsl_eigen_symmv (XC, V(l), M(v), w);
-    CHECK(res,res);
-    gsl_eigen_symmv_free (w);
-    gsl_matrix_free(XC);
-    gsl_eigen_symmv_sort (V(l), M(v), GSL_EIGEN_SORT_ABS_DESC);
-    OK
-}
-
-// for hermitian matrices
-int eigensystemC(KCMAT(x),RVEC(l),CMAT(v)) {
-    REQUIRES(xr==xc && xr==ln && xr==vr && vr==vc,BAD_SIZE);
-    DEBUGMSG("eigensystemC");
-    KCMVIEW(x);
-    DVVIEW(l);
-    CMVIEW(v);
-    gsl_matrix_complex *XC = gsl_matrix_complex_alloc(xr,xr);
-    gsl_matrix_complex_memcpy(XC,M(x)); // again needed because the argument is destroyed
-    gsl_eigen_hermv_workspace * w = gsl_eigen_hermv_alloc (xc);
-    int res = gsl_eigen_hermv (XC, V(l), M(v), w);
-    CHECK(res,res);
-    gsl_eigen_hermv_free (w);
-    gsl_matrix_complex_free(XC);
-    gsl_eigen_hermv_sort (V(l), M(v), GSL_EIGEN_SORT_ABS_DESC);
-    OK
-}
-
-int QR(KRMAT(x),RMAT(q),RMAT(r)) {
-    REQUIRES(xr==rr && xc==rc && qr==qc && xr==qr,BAD_SIZE);
-    DEBUGMSG("QR");
-    KDMVIEW(x);
-    DMVIEW(q);
-    DMVIEW(r);
-    gsl_matrix * a = gsl_matrix_alloc(xr,xc);
-    gsl_vector * tau = gsl_vector_alloc(MIN(xr,xc));
-    gsl_matrix_memcpy(a,M(x));
-    int res = gsl_linalg_QR_decomp(a,tau);
-    CHECK(res,res);
-    gsl_linalg_QR_unpack(a,tau,M(q),M(r));
-    gsl_vector_free(tau);
-    gsl_matrix_free(a);
-    OK
-}
-
-int QRpacked(KRMAT(x),RMAT(qr),RVEC(tau)) {
-    //REQUIRES(xr==rr && xc==rc && qr==qc && xr==qr,BAD_SIZE);
-    DEBUGMSG("QRpacked");
-    KDMVIEW(x);
-    DMVIEW(qr);
-    DVVIEW(tau);
-    //gsl_matrix * a = gsl_matrix_alloc(xr,xc);
-    //gsl_vector * tau = gsl_vector_alloc(MIN(xr,xc));
-    gsl_matrix_memcpy(M(qr),M(x));
-    int res = gsl_linalg_QR_decomp(M(qr),V(tau));
-    CHECK(res,res);
-    OK
-}
-
-
-int QRunpack(KRMAT(xqr),KRVEC(tau),RMAT(q),RMAT(r)) {
-    //REQUIRES(xr==rr && xc==rc && qr==qc && xr==qr,BAD_SIZE);
-    DEBUGMSG("QRunpack");
-    KDMVIEW(xqr);
-    KDVVIEW(tau);
-    DMVIEW(q);
-    DMVIEW(r);
-    gsl_linalg_QR_unpack(M(xqr),V(tau),M(q),M(r));
-    OK
-}
-
-
-int cholR(KRMAT(x),RMAT(l)) {
-    REQUIRES(xr==xc && lr==xr && lr==lc,BAD_SIZE);
-    DEBUGMSG("cholR");
-    KDMVIEW(x);
-    DMVIEW(l);
-    gsl_matrix_memcpy(M(l),M(x));
-    int res = gsl_linalg_cholesky_decomp(M(l));
-    CHECK(res,res);
-    int r,c;
-    for (r=0; r<xr-1; r++) {
-        for(c=r+1; c<xc; c++) {
-            lp[r*lc+c] = 0.;
-        }
-    }
-    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/Numeric/GSL/gsl-aux.h b/lib/Numeric/GSL/gsl-aux.h
--- a/lib/Numeric/GSL/gsl-aux.h
+++ b/lib/Numeric/GSL/gsl-aux.h
@@ -26,25 +26,6 @@
 int zipC(int code, KCVEC(a), KCVEC(b), CVEC(r));
 
 
-int luSolveR(KRMAT(a),KRMAT(b),RMAT(r)); 
-int luSolveC(KCMAT(a),KCMAT(b),CMAT(r));
-int luRaux(KRMAT(a),RVEC(b));
-int luCaux(KCMAT(a),CVEC(b));
-
-int svd(KRMAT(x),RMAT(u), RVEC(s),RMAT(v));
-
-int eigensystemR(KRMAT(x),RVEC(l),RMAT(v));
-int eigensystemC(KCMAT(x),RVEC(l),CMAT(v));
-
-int QR(KRMAT(x),RMAT(q),RMAT(r));
-
-int QRpacked(KRMAT(x),RMAT(qr),RVEC(tau));
-int QRunpack(KRMAT(qr),KRVEC(tau),RMAT(q),RMAT(r));
-
-int cholR(KRMAT(x),RMAT(l));
-
-int cholC(KCMAT(x),CMAT(l));
-
 int fft(int code, KCVEC(a), CVEC(b));
 
 int integrate_qng(double f(double, void*), double a, double b, double prec,
diff --git a/lib/Numeric/LinearAlgebra.hs b/lib/Numeric/LinearAlgebra.hs
--- a/lib/Numeric/LinearAlgebra.hs
+++ b/lib/Numeric/LinearAlgebra.hs
@@ -16,7 +16,6 @@
     module Data.Packed,
     module Numeric.LinearAlgebra.Linear,
     module Numeric.LinearAlgebra.Algorithms,
-    module Numeric.LinearAlgebra.Instances,
     module Numeric.LinearAlgebra.Interface
 ) where
 
diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs
--- a/lib/Numeric/LinearAlgebra/Algorithms.hs
+++ b/lib/Numeric/LinearAlgebra/Algorithms.hs
@@ -1,4 +1,5 @@
-{-# OPTIONS_GHC -fglasgow-exts #-}
+{-# OPTIONS_GHC -XFlexibleContexts -XFlexibleInstances #-}
+{-# LANGUAGE CPP #-}
 -----------------------------------------------------------------------------
 {- |
 Module      :  Numeric.LinearAlgebra.Algorithms
@@ -20,6 +21,7 @@
 
 module Numeric.LinearAlgebra.Algorithms (
 -- * Linear Systems
+    multiply, dot,
     linearSolve,
     inv, pinv,
     pinvTol, det, rank, rcond,
@@ -38,7 +40,7 @@
 -- ** Schur
     schur,
 -- ** LU
-    lu,
+    lu, luPacked, luSolve,
 -- * Matrix functions
     expm,
     sqrtm,
@@ -51,6 +53,7 @@
 -- * Misc
     ctrans,
     eps, i,
+    outer, kronecker,
 -- * Util
     haussholder,
     unpackQR, unpackHess,
@@ -60,7 +63,6 @@
 
 import Data.Packed.Internal hiding (fromComplex, toComplex, comp, conj, (//))
 import Data.Packed
-import qualified Numeric.GSL.Matrix as GSL
 import Numeric.GSL.Vector
 import Numeric.LinearAlgebra.LAPACK as LAPACK
 import Complex
@@ -68,12 +70,19 @@
 import Data.List(foldl1')
 import Data.Array
 
+
+
 -- | Auxiliary typeclass used to define generic computations for both real and complex matrices.
 class (Normed (Matrix t), Linear Vector t, Linear Matrix t) => Field t where
     -- | Singular value decomposition using lapack's dgesvd or zgesvd.
     svd         :: Matrix t -> (Matrix t, Vector Double, Matrix t)
+    -- | Obtains the LU decomposition of a matrix in a compact data structure suitable for 'luSolve'.
     luPacked    :: Matrix t -> (Matrix t, [Int])
-    -- | Solution of a general linear system (for several right-hand sides) using lapacks' dgesv and zgesv.
+    -- | Solution of a linear system (for several right hand sides) from the precomputed LU factorization
+    --   obtained by 'luPacked'.
+    luSolve     :: (Matrix t, [Int]) -> Matrix t -> Matrix t
+    -- | Solution of a general linear system (for several right-hand sides) using lapacks' dgesv or zgesv.
+    -- It is similar to 'luSolve' . 'luPacked', but @linearSolve@ raises an error if called on a singular system.
     --  See also other versions of linearSolve in "Numeric.LinearAlgebra.LAPACK".
     linearSolve :: Matrix t -> Matrix t -> Matrix t
     linearSolveSVD :: Matrix t -> Matrix t -> Matrix t
@@ -105,24 +114,29 @@
     schur       :: Matrix t -> (Matrix t, Matrix t)
     -- | Conjugate transpose.
     ctrans :: Matrix t -> Matrix t
+    -- | Matrix product.
+    multiply :: Matrix t -> Matrix t -> Matrix t
 
 
 instance Field Double where
     svd = svdR
     luPacked = luR
-    linearSolve = linearSolveR
+    luSolve (l_u,perm) = lusR l_u perm
+    linearSolve = linearSolveR                 -- (luSolve . luPacked) ??
     linearSolveSVD = linearSolveSVDR Nothing
     ctrans = trans
     eig = eigR
     eigSH' = eigS
     cholSH = cholS
-    qr = GSL.unpackQR . qrR
+    qr = unpackQR . qrR
     hess = unpackHess hessR
     schur = schurR
+    multiply = multiplyR
 
 instance Field (Complex Double) where
     svd = svdC
     luPacked = luC
+    luSolve (l_u,perm) = lusC l_u perm
     linearSolve = linearSolveC
     linearSolveSVD = linearSolveSVDC Nothing
     ctrans = conj . trans
@@ -132,7 +146,9 @@
     qr = unpackQR . qrC
     hess = unpackHess hessC
     schur = schurC
+    multiply = multiplyC
 
+
 -- | Eigenvalues and Eigenvectors of a complex hermitian or real symmetric matrix using lapack's dsyev or zheev.
 --
 -- If @(s,v) = eigSH m@ then @m == v \<> diag s \<> ctrans v@
@@ -151,12 +167,12 @@
 
 -- | determinant of a square matrix, computed from the LU decomposition.
 det :: Field t => Matrix t -> t
-det m | square m = s * (product $ toList $ takeDiag $ lu)
+det m | square m = s * (product $ toList $ takeDiag $ lup)
       | otherwise = error "det of nonsquare matrix"
-    where (lu,perm) = luPacked m
+    where (lup,perm) = luPacked m
           s = signlp (rows m) perm
 
--- | LU factorization of a general matrix using lapack's dgetrf or zgetrf.
+-- | Explicit LU factorization of a general matrix using lapack's dgetrf or zgetrf.
 --
 -- If @(l,u,p,s) = lu m@ then @m == p \<> l \<> u@, where l is lower triangular,
 -- u is upper triangular, p is a permutation matrix and s is the signature of the permutation.
@@ -484,20 +500,69 @@
           s = toColumns (ident r)
           (res,sign) = foldl swap (listArray (0,r-1) s, 1) (zip v vals)
 
-triang r c h v = reshape c $ fromList [el i j | i<-[0..r-1], j<-[0..c-1]]
-    where el i j = if j-i>=h then v else 1 - v
+triang r c h v = (r><c) [el s t | s<-[0..r-1], t<-[0..c-1]]
+    where el p q = if q-p>=h then v else 1 - v
 
-luFact (lu,perm) | r <= c    = (l ,u ,p, s)
-                 | otherwise = (l',u',p, s)
+luFact (l_u,perm) | r <= c    = (l ,u ,p, s)
+                  | otherwise = (l',u',p, s)
   where
-    r = rows lu
-    c = cols lu
+    r = rows l_u
+    c = cols l_u
     tu = triang r c 0 1
     tl = triang r c 0 0
-    l = takeColumns r (lu |*| tl) |+| diagRect (constant 1 r) r r
-    u = lu |*| tu
+    l = takeColumns r (l_u |*| tl) |+| diagRect (constant 1 r) r r
+    u = l_u |*| tu
     (p,s) = fixPerm r perm
-    l' = (lu |*| tl) |+| diagRect (constant 1 c) r c
-    u' = takeRows c (lu |*| tu)
+    l' = (l_u |*| tl) |+| diagRect (constant 1 c) r c
+    u' = takeRows c (l_u |*| tu)
     (|+|) = add
     (|*|) = mul
+
+--------------------------------------------------
+
+-- | Euclidean inner product.
+dot :: (Field t) => Vector t -> Vector t -> t
+dot u v = multiply r c  @@> (0,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
+
+{- | Kronecker product of two matrices.
+
+@m1=(2><3)
+ [ 1.0,  2.0, 0.0
+ , 0.0, -1.0, 3.0 ]
+m2=(4><3)
+ [  1.0,  2.0,  3.0
+ ,  4.0,  5.0,  6.0
+ ,  7.0,  8.0,  9.0
+ , 10.0, 11.0, 12.0 ]@
+
+@\> kronecker m1 m2
+(8><9)
+ [  1.0,  2.0,  3.0,   2.0,   4.0,   6.0,  0.0,  0.0,  0.0
+ ,  4.0,  5.0,  6.0,   8.0,  10.0,  12.0,  0.0,  0.0,  0.0
+ ,  7.0,  8.0,  9.0,  14.0,  16.0,  18.0,  0.0,  0.0,  0.0
+ , 10.0, 11.0, 12.0,  20.0,  22.0,  24.0,  0.0,  0.0,  0.0
+ ,  0.0,  0.0,  0.0,  -1.0,  -2.0,  -3.0,  3.0,  6.0,  9.0
+ ,  0.0,  0.0,  0.0,  -4.0,  -5.0,  -6.0, 12.0, 15.0, 18.0
+ ,  0.0,  0.0,  0.0,  -7.0,  -8.0,  -9.0, 21.0, 24.0, 27.0
+ ,  0.0,  0.0,  0.0, -10.0, -11.0, -12.0, 30.0, 33.0, 36.0 ]@
+-}
+kronecker :: (Field t) => Matrix t -> Matrix t -> Matrix t
+kronecker a b = fromBlocks
+              . partit (cols a)
+              . map (reshape (cols b))
+              . toRows
+              $ flatten a `outer` flatten b
diff --git a/lib/Numeric/LinearAlgebra/Instances.hs b/lib/Numeric/LinearAlgebra/Instances.hs
--- a/lib/Numeric/LinearAlgebra/Instances.hs
+++ b/lib/Numeric/LinearAlgebra/Instances.hs
@@ -1,4 +1,4 @@
-{-# OPTIONS_GHC -fglasgow-exts -fallow-undecidable-instances #-}
+{-# LANGUAGE UndecidableInstances, FlexibleInstances #-}
 -----------------------------------------------------------------------------
 {- |
 Module      :  Numeric.LinearAlgebra.Instances
@@ -22,7 +22,6 @@
 import Numeric.LinearAlgebra.Linear
 import Numeric.GSL.Vector
 import Data.Packed.Matrix
-import Data.Packed.Vector
 import Complex
 import Data.List(transpose,intersperse)
 import Foreign(Storable)
@@ -49,11 +48,11 @@
 ------------------------------------------------------------------
 
 instance (Element a, Read a) => Read (Matrix a) where
-    readsPrec _ s = [((rows><cols) . read $ listnums, rest)]
+    readsPrec _ s = [((rs><cs) . read $ listnums, rest)]
         where (thing,rest) = breakAt ']' s
               (dims,listnums) = breakAt ')' thing
-              cols = read . init . fst. breakAt ')' . snd . breakAt '<' $ dims
-              rows = read . snd . breakAt '(' .init . fst . breakAt '>' $ dims
+              cs = read . init . fst. breakAt ')' . snd . breakAt '<' $ dims
+              rs = read . snd . breakAt '(' .init . fst . breakAt '>' $ dims
 
 instance (Element a, Read a) => Read (Vector a) where
     readsPrec _ s = [((d |>) . read $ listnums, rest)]
@@ -83,18 +82,26 @@
 instance (Eq a, Element a) => Eq (Vector a) where
     a == b = dim a == dim b && toList a == toList b
 
-instance (Linear Vector a) => Num (Vector a) where
+instance Num (Vector Double) where
     (+) = adaptScalar addConstant add (flip addConstant)
     negate = scale (-1)
     (*) = adaptScalar scale mul (flip scale)
-    signum = liftVector signum
-    abs = liftVector abs
+    signum = vectorMapR Sign
+    abs = vectorMapR Abs
     fromInteger = fromList . return . fromInteger
 
+instance Num (Vector (Complex Double)) where
+    (+) = adaptScalar addConstant add (flip addConstant)
+    negate = scale (-1)
+    (*) = adaptScalar scale mul (flip scale)
+    signum = vectorMapC Sign
+    abs = vectorMapC Abs
+    fromInteger = fromList . return . fromInteger
+
 instance (Eq a, Element a) => Eq (Matrix a) where
     a == b = cols a == cols b && flatten a == flatten b
 
-instance (Linear Vector a) => Num (Matrix a) where
+instance (Linear Matrix a, Num (Vector a)) => Num (Matrix a) where
     (+) = liftMatrix2' (+)
     (-) = liftMatrix2' (-)
     negate = liftMatrix negate
@@ -105,7 +112,7 @@
 
 ---------------------------------------------------
 
-instance (Linear Vector a) => Fractional (Vector a) where
+instance (Linear Vector a, Num (Vector a)) => Fractional (Vector a) where
     fromRational n = fromList [fromRational n]
     (/) = adaptScalar f divide g where
         r `f` v = scaleRecip r v
@@ -113,7 +120,7 @@
 
 -------------------------------------------------------
 
-instance (Linear Vector a, Fractional (Vector a)) => Fractional (Matrix a) where
+instance (Linear Vector a, Fractional (Vector a), Num (Matrix a)) => Fractional (Matrix a) where
     fromRational n = (1><1) [fromRational n]
     (/) = liftMatrix2' (/)
 
@@ -161,7 +168,7 @@
 
 -----------------------------------------------------------
 
-instance (Linear Vector a, Floating (Vector a)) => Floating (Matrix a) where
+instance (Linear Vector a, Floating (Vector a), Fractional (Matrix a)) => Floating (Matrix a) where
     sin   = liftMatrix sin
     cos   = liftMatrix cos
     tan   = liftMatrix tan
@@ -188,3 +195,4 @@
     mconcat = j . filter ((>0).dim)
         where j [] = mempty
               j l  = join l
+
diff --git a/lib/Numeric/LinearAlgebra/Interface.hs b/lib/Numeric/LinearAlgebra/Interface.hs
--- a/lib/Numeric/LinearAlgebra/Interface.hs
+++ b/lib/Numeric/LinearAlgebra/Interface.hs
@@ -29,7 +29,7 @@
 class Mul a b c | a b -> c where
  infixl 7 <>
  -- | matrix product
- (<>) :: Element t => a t -> b t -> c t
+ (<>) :: Field t => a t -> b t -> c t
 
 instance Mul Matrix Matrix Matrix where
     (<>) = multiply
@@ -43,7 +43,7 @@
 ---------------------------------------------------
 
 -- | @u \<.\> v = dot u v@
-(<.>) :: (Element t) => Vector t -> Vector t -> t
+(<.>) :: (Field t) => Vector t -> Vector t -> t
 infixl 7 <.>
 (<.>) = dot
 
diff --git a/lib/Numeric/LinearAlgebra/LAPACK.hs b/lib/Numeric/LinearAlgebra/LAPACK.hs
--- a/lib/Numeric/LinearAlgebra/LAPACK.hs
+++ b/lib/Numeric/LinearAlgebra/LAPACK.hs
@@ -14,12 +14,13 @@
 -----------------------------------------------------------------------------
 
 module Numeric.LinearAlgebra.LAPACK (
+    multiplyR, multiplyC,
     svdR, svdRdd, svdC,
     eigC, eigR, eigS, eigH, eigS', eigH',
     linearSolveR, linearSolveC,
     linearSolveLSR, linearSolveLSC,
     linearSolveSVDR, linearSolveSVDC,
-    luR, luC,
+    luR, luC, lusR, lusC,
     cholS, cholH,
     qrR, qrC,
     hessR, hessC,
@@ -34,7 +35,35 @@
 import Numeric.GSL.Vector(vectorMapValR, FunCodeSV(Scale))
 import Complex
 import Foreign
+import Foreign.C.Types (CInt)
+import Control.Monad(when)
 
+-----------------------------------------------------------------------------------
+
+foreign import ccall "LAPACK/lapack-aux.h multiplyR" dgemmc :: CInt -> CInt -> TMMM
+foreign import ccall "LAPACK/lapack-aux.h multiplyC" zgemmc :: CInt -> CInt -> TCMCMCM
+
+isT MF{} = 0
+isT MC{} = 1
+
+tt x@MF{} = x
+tt x@MC{} = trans x
+
+multiplyAux f st a b = unsafePerformIO $ do
+    when (cols a /= rows b) $ error $ "inconsistent dimensions in matrix product "++
+                                       show (rows a,cols a) ++ " x " ++ show (rows b, cols b)
+    s <- createMatrix ColumnMajor (rows a) (cols b)
+    app3 (f (isT a) (isT b)) mat (tt a) mat (tt b) mat s st
+    return s
+
+-- | Matrix product based on BLAS's /dgemm/.
+multiplyR :: Matrix Double -> Matrix Double -> Matrix Double
+multiplyR a b = multiplyAux dgemmc "dgemmc" a b
+
+-- | Matrix product based on BLAS's /zgemm/.
+multiplyC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
+multiplyC a b = multiplyAux zgemmc "zgemmc" a b
+
 -----------------------------------------------------------------------------
 foreign import ccall "LAPACK/lapack-aux.h svd_l_R" dgesvd :: TMMVM
 foreign import ccall "LAPACK/lapack-aux.h svd_l_C" zgesvd :: TCMCMVCM
@@ -337,3 +366,30 @@
     return (lu, map (pred.round) (toList piv))
   where n = rows a
         m = cols a
+
+-----------------------------------------------------------------------------------
+type TW a = CInt -> PD -> a
+type TQ a = CInt -> CInt -> PC -> a
+
+foreign import ccall "LAPACK/lapack-aux.h luS_l_R" dgetrs :: TMVMM
+foreign import ccall "LAPACK/lapack-aux.h luS_l_C" zgetrs :: TQ (TW (TQ (TQ (IO CInt))))
+
+-- | Wrapper for LAPACK's /dgetrs/, which solves a general real linear system (for several right-hand sides) from a precomputed LU decomposition.
+lusR :: Matrix Double -> [Int] -> Matrix Double -> Matrix Double
+lusR a piv b = lusAux dgetrs "lusR" (fmat a) piv (fmat b)
+
+-- | Wrapper for LAPACK's /zgetrs/, which solves a general real linear system (for several right-hand sides) from a precomputed LU decomposition.
+lusC :: Matrix (Complex Double) -> [Int] -> Matrix (Complex Double) -> Matrix (Complex Double)
+lusC a piv b = lusAux zgetrs "lusC" (fmat a) piv (fmat b)
+
+lusAux f st a piv b
+    | n1==n2 && n2==n =unsafePerformIO $ do
+         x <- createMatrix ColumnMajor n m
+         app4 f mat a vec piv' mat b mat x st
+         return x
+    | otherwise = error $ st ++ " on LU factorization of nonsquare matrix"
+  where n1 = rows a
+        n2 = cols a
+        n = rows b
+        m = cols b
+        piv' = fromList (map (fromIntegral.succ) piv) :: Vector Double
diff --git a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
--- a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
+++ b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
@@ -20,6 +20,9 @@
 #define OK return 0;
 #endif
 
+#define TRACEMAT(M) {int q; printf(" %d x %d: ",M##r,M##c); \
+                     for(q=0;q<M##r*M##c;q++) printf("%.1f ",M##p[q]); printf("\n");}
+
 #define CHECK(RES,CODE) MACRO(if(RES) return CODE;)
 
 #define BAD_SIZE 2000
@@ -812,5 +815,89 @@
         ipivp[k] = auxipiv[k];
     }
     free(auxipiv);
+    OK
+}
+
+
+//////////////////// LU substitution /////////////////////////
+
+int luS_l_R(KDMAT(a), KDVEC(ipiv), KDMAT(b), DMAT(x)) {
+  integer m = ar;
+  integer n = ac;
+  integer mrhs = br;
+  integer nrhs = bc;
+
+  REQUIRES(m==n && m==mrhs && m==ipivn,BAD_SIZE);
+  integer* auxipiv = (integer*)malloc(n*sizeof(integer));
+  int k;
+  for (k=0; k<n; k++) {
+    auxipiv[k] = (integer)ipivp[k];
+  }
+  integer res;
+  memcpy(xp,bp,mrhs*nrhs*sizeof(double));
+  dgetrs_ ("N",&n,&nrhs,(/*no const (!?)*/ double*)ap,&m,auxipiv,xp,&mrhs,&res);
+  CHECK(res,res);
+  free(auxipiv);
+  OK
+}
+
+int luS_l_C(KCMAT(a), KDVEC(ipiv), KCMAT(b), CMAT(x)) {
+    integer m = ar;
+    integer n = ac;
+    integer mrhs = br;
+    integer nrhs = bc;
+
+    REQUIRES(m==n && m==mrhs && m==ipivn,BAD_SIZE);
+    integer* auxipiv = (integer*)malloc(n*sizeof(integer));
+    int k;
+    for (k=0; k<n; k++) {
+        auxipiv[k] = (integer)ipivp[k];
+    }
+    integer res;
+    memcpy(xp,bp,mrhs*nrhs*sizeof(doublecomplex));
+    zgetrs_ ("N",&n,&nrhs,(doublecomplex*)ap,&m,auxipiv,(doublecomplex*)xp,&mrhs,&res);
+    CHECK(res,res);
+    free(auxipiv);
+    OK
+}
+
+//////////////////// Matrix Product /////////////////////////
+
+void dgemm_(char *, char *, integer *, integer *, integer *,
+           double *, const double *, integer *, const double *,
+           integer *, double *, double *, integer *);
+
+int multiplyR(int ta, int tb, KDMAT(a),KDMAT(b),DMAT(r)) {
+    //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);
+    integer m = ta?ac:ar;
+    integer n = tb?br:bc;
+    integer k = ta?ar:ac;
+    integer lda = ar;
+    integer ldb = br;
+    integer ldc = rr;
+    double alpha = 1;
+    double beta = 0;
+    dgemm_(ta?"T":"N",tb?"T":"N",&m,&n,&k,&alpha,ap,&lda,bp,&ldb,&beta,rp,&ldc);
+    OK
+}
+
+void zgemm_(char *, char *, integer *, integer *, integer *,
+           doublecomplex *, const doublecomplex *, integer *, const doublecomplex *,
+           integer *, doublecomplex *, doublecomplex *, integer *);
+
+int multiplyC(int ta, int tb, KCMAT(a),KCMAT(b),CMAT(r)) {
+    //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);
+    integer m = ta?ac:ar;
+    integer n = tb?br:bc;
+    integer k = ta?ar:ac;
+    integer lda = ar;
+    integer ldb = br;
+    integer ldc = rr;
+    doublecomplex alpha = {1,0};
+    doublecomplex beta = {0,0};
+    zgemm_(ta?"T":"N",tb?"T":"N",&m,&n,&k,&alpha,
+           (doublecomplex*)ap,&lda,
+           (doublecomplex*)bp,&ldb,&beta,
+           (doublecomplex*)rp,&ldc);
     OK
 }
diff --git a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h
--- a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h
+++ b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h
@@ -45,12 +45,16 @@
 #define DMAT(A) int A##r, int A##c, double* A##p
 #define CMAT(A) int A##r, int A##c, double* A##p
 
-// const pointer versions for the parameters 
 #define KDVEC(A) int A##n, const double*A##p
 #define KCVEC(A) int A##n, const double*A##p
 #define KDMAT(A) int A##r, int A##c, const double* A##p
-#define KCMAT(A) int A##r, int A##c, const double* A##p 
+#define KCMAT(A) int A##r, int A##c, const double* A##p
 
+/********************************************************/
+
+int multiplyR(int ta, int tb, KDMAT(a),KDMAT(b),DMAT(r));
+int multiplyC(int ta, int tb, KCMAT(a),KCMAT(b),CMAT(r));
+
 int svd_l_R(KDMAT(x),DMAT(u),DVEC(s),DMAT(v));
 int svd_l_Rdd(KDMAT(x),DMAT(u),DVEC(s),DMAT(v));
 int svd_l_C(KCMAT(a),CMAT(u),DVEC(s),CMAT(v));
@@ -84,3 +88,6 @@
 
 int lu_l_R(KDMAT(a), DVEC(ipiv), DMAT(r));
 int lu_l_C(KCMAT(a), DVEC(ipiv), CMAT(r));
+
+int luS_l_R(KDMAT(a), KDVEC(ipiv), KDMAT(b), DMAT(x));
+int luS_l_C(KCMAT(a), KDVEC(ipiv), KCMAT(b), CMAT(x));
diff --git a/lib/Numeric/LinearAlgebra/Linear.hs b/lib/Numeric/LinearAlgebra/Linear.hs
--- a/lib/Numeric/LinearAlgebra/Linear.hs
+++ b/lib/Numeric/LinearAlgebra/Linear.hs
@@ -1,4 +1,4 @@
-{-# OPTIONS_GHC -fglasgow-exts -fallow-undecidable-instances #-}
+{-# LANGUAGE UndecidableInstances, MultiParamTypeClasses, FlexibleInstances #-}
 -----------------------------------------------------------------------------
 {- |
 Module      :  Numeric.LinearAlgebra.Linear
@@ -15,12 +15,9 @@
 -----------------------------------------------------------------------------
 
 module Numeric.LinearAlgebra.Linear (
-    Linear(..),
-    multiply, dot, outer, kronecker
+    Linear(..)
 ) where
 
-
-import Data.Packed.Internal(multiply,partit)
 import Data.Packed
 import Numeric.GSL.Vector
 import Complex
@@ -69,52 +66,3 @@
     mul = liftMatrix2 mul
     divide = liftMatrix2 divide
     equal a b = cols a == cols b && flatten a `equal` flatten b
-
---------------------------------------------------
-
--- | euclidean inner product
-dot :: (Element t) => Vector t -> Vector t -> t
-dot u v = multiply r c  @@> (0,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 :: (Element t) => Vector t -> Vector t -> Matrix t
-outer u v = asColumn u `multiply` asRow v
-
-{- | Kronecker product of two matrices.
-
-@m1=(2><3)
- [ 1.0,  2.0, 0.0
- , 0.0, -1.0, 3.0 ]
-m2=(4><3)
- [  1.0,  2.0,  3.0
- ,  4.0,  5.0,  6.0
- ,  7.0,  8.0,  9.0
- , 10.0, 11.0, 12.0 ]@
-
-@\> kronecker m1 m2
-(8><9)
- [  1.0,  2.0,  3.0,   2.0,   4.0,   6.0,  0.0,  0.0,  0.0
- ,  4.0,  5.0,  6.0,   8.0,  10.0,  12.0,  0.0,  0.0,  0.0
- ,  7.0,  8.0,  9.0,  14.0,  16.0,  18.0,  0.0,  0.0,  0.0
- , 10.0, 11.0, 12.0,  20.0,  22.0,  24.0,  0.0,  0.0,  0.0
- ,  0.0,  0.0,  0.0,  -1.0,  -2.0,  -3.0,  3.0,  6.0,  9.0
- ,  0.0,  0.0,  0.0,  -4.0,  -5.0,  -6.0, 12.0, 15.0, 18.0
- ,  0.0,  0.0,  0.0,  -7.0,  -8.0,  -9.0, 21.0, 24.0, 27.0
- ,  0.0,  0.0,  0.0, -10.0, -11.0, -12.0, 30.0, 33.0, 36.0 ]@
--}
-kronecker :: (Element t) => Matrix t -> Matrix t -> Matrix t
-kronecker a b = fromBlocks
-              . partit (cols a)
-              . map (reshape (cols b))
-              . toRows
-              $ flatten a `outer` flatten b
diff --git a/lib/Numeric/LinearAlgebra/Tests.hs b/lib/Numeric/LinearAlgebra/Tests.hs
--- a/lib/Numeric/LinearAlgebra/Tests.hs
+++ b/lib/Numeric/LinearAlgebra/Tests.hs
@@ -22,12 +22,16 @@
 import Numeric.LinearAlgebra
 import Numeric.LinearAlgebra.Tests.Instances
 import Numeric.LinearAlgebra.Tests.Properties
-import Test.QuickCheck
+import Test.QuickCheck hiding (test)
 import Test.HUnit hiding ((~:),test)
 import System.Info
 import Data.List(foldl1')
 import Numeric.GSL hiding (sin,cos,exp,choose)
+import Prelude hiding ((^))
+import qualified Prelude
 
+a ^ b = a Prelude.^ (b :: Int)
+
 qCheck n = check defaultConfig {configSize = const n}
 
 utest str b = TestCase $ assertBool str b
@@ -73,7 +77,7 @@
           expected = -0.17759677131433830434739701
 
 exponentialTest = utest "exp_e10_e" ( abs (v*10^e - expected) < 4E-2 )
-    where (v,e,err) = exp_e10_e 30.0
+    where (v,e,_err) = exp_e10_e 30.0
           expected = exp 30.0
 
 ---------------------------------------------------------------------
@@ -115,7 +119,6 @@
     where fun n = foldl1' (<>) (map rot angles)
               where angles = toList $ linspace n (0,1)
 
-
 -- | All tests must pass with a maximum dimension of about 20
 --  (some tests may fail with bigger sizes due to precision loss).
 runTests :: Int  -- ^ maximum dimension
@@ -123,13 +126,21 @@
 runTests n = do
     setErrorHandlerOff
     let test p = qCheck n p
+    putStrLn "------ mult"
+    test (multProp1  . rConsist)
+    test (multProp1  . cConsist)
+    test (multProp2  . rConsist)
+    test (multProp2  . cConsist)
     putStrLn "------ lu"
     test (luProp    . rM)
     test (luProp    . cM)
-    putStrLn "------ inv"
+    putStrLn "------ inv (linearSolve)"
     test (invProp   . rSqWC)
     test (invProp   . cSqWC)
-    putStrLn "------ pinv"
+    putStrLn "------ luSolve"
+    test (linearSolveProp (luSolve.luPacked) . rSqWC)
+    test (linearSolveProp (luSolve.luPacked) . cSqWC)
+    putStrLn "------ pinv (linearSolveSVD)"
     test (pinvProp  . rM)
     if os == "mingw32"
         then putStrLn "complex pinvTest skipped in this OS"
diff --git a/lib/Numeric/LinearAlgebra/Tests/Instances.hs b/lib/Numeric/LinearAlgebra/Tests/Instances.hs
--- a/lib/Numeric/LinearAlgebra/Tests/Instances.hs
+++ b/lib/Numeric/LinearAlgebra/Tests/Instances.hs
@@ -1,4 +1,4 @@
-{-# OPTIONS #-}
+{-# LANGUAGE FlexibleContexts, UndecidableInstances #-}
 -----------------------------------------------------------------------------
 {- |
 Module      :  Numeric.LinearAlgebra.Tests.Instances
@@ -20,6 +20,7 @@
     WC(..),     rWC,cWC,
     SqWC(..),   rSqWC, cSqWC,
     PosDef(..), rPosDef, cPosDef,
+    Consistent(..), rConsist, cConsist,
     RM,CM, rM,cM
 ) where
 
@@ -29,9 +30,9 @@
 
 instance (Arbitrary a, RealFloat a) => Arbitrary (Complex a) where
     arbitrary = do
-        r <- arbitrary
-        i <- arbitrary
-        return (r:+i)
+        re <- arbitrary
+        im <- arbitrary
+        return (re :+ im)
     coarbitrary = undefined
 
 chooseDim = sized $ \m -> choose (1,max 1 m)
@@ -70,7 +71,7 @@
 
 -- a complex hermitian or real symmetric matrix
 newtype (Her a) = Her (Matrix a) deriving Show
-instance (Field a, Arbitrary a) => Arbitrary (Her a) where
+instance (Field a, Arbitrary a, Num (Vector a)) => Arbitrary (Her a) where
     arbitrary = do
         Sq m <- arbitrary
         let m' = m/2
@@ -105,7 +106,7 @@
 
 -- a positive definite square matrix (the eigenvalues are between 0 and 100)
 newtype (PosDef a) = PosDef (Matrix a) deriving Show
-instance (Field a, Arbitrary a) => Arbitrary (PosDef a) where
+instance (Field a, Arbitrary a, Num (Vector a)) => Arbitrary (PosDef a) where
     arbitrary = do
         Her m <- arbitrary
         let (_,v) = eigSH m
@@ -116,6 +117,19 @@
         return $ PosDef (0.5 .* p + 0.5 .* ctrans p)
     coarbitrary = undefined
 
+-- a pair of matrices that can be multiplied
+newtype (Consistent a) = Consistent (Matrix a, Matrix a) deriving Show
+instance (Field a, Arbitrary a) => Arbitrary (Consistent a) where
+    arbitrary = do
+        n <- chooseDim
+        k <- chooseDim
+        m <- chooseDim
+        la <- vector (n*k)
+        lb <- vector (k*m)
+        return $ Consistent ((n><k) la, (k><m) lb)
+    coarbitrary = undefined
+
+
 type RM = Matrix Double
 type CM = Matrix (Complex Double)
 
@@ -140,3 +154,5 @@
 rPosDef (PosDef m) = m :: RM
 cPosDef (PosDef m) = m :: CM
 
+rConsist (Consistent (a,b)) = (a,b::RM)
+cConsist (Consistent (a,b)) = (a,b::CM)
diff --git a/lib/Numeric/LinearAlgebra/Tests/Properties.hs b/lib/Numeric/LinearAlgebra/Tests/Properties.hs
--- a/lib/Numeric/LinearAlgebra/Tests/Properties.hs
+++ b/lib/Numeric/LinearAlgebra/Tests/Properties.hs
@@ -34,15 +34,16 @@
     hessProp,
     schurProp1, schurProp2,
     cholProp,
-    expmDiagProp
+    expmDiagProp,
+    multProp1, multProp2,
+    linearSolveProp
 ) where
 
 import Numeric.LinearAlgebra
-import Numeric.LinearAlgebra.Tests.Instances(Sq(..),Her(..),Rot(..))
 import Test.QuickCheck
-import Debug.Trace
+-- import Debug.Trace
 
-debug x = trace (show x) x
+-- debug x = trace (show x) x
 
 -- relative error
 dist :: (Normed t, Num t) => t -> t -> Double
@@ -76,7 +77,7 @@
 wellCond m = rcond m > 1/100
 
 positiveDefinite m = minimum (toList e) > 0
-    where (e,v) = eigSH m
+    where (e,_v) = eigSH m
 
 upperTriang m = rows m == 1 || down == z
     where down = fromList $ concat $ zipWith drop [1..] (toLists (ctrans m))
@@ -106,8 +107,8 @@
 
 detProp m = s d1 |~| s d2
     where d1 = det m
-          d2 = det' m * det q
-          det' m = product $ toList $ takeDiag r
+          d2 = det' * det q
+          det' = product $ toList $ takeDiag r
           (q,r) = qr m
           s x = fromList [x]
 
@@ -146,8 +147,17 @@
 
 cholProp m = m |~| ctrans c <> c && upperTriang c
     where c = chol m
-          pos = positiveDefinite m
+          -- pos = positiveDefinite m
 
 expmDiagProp m = expm (logm m) :~ 7 ~: complex m
-    where logm m = matFunc log m
+    where logm = matFunc log
 
+-- reference multiply
+mulH a b = fromLists [[ doth ai bj | bj <- toColumns b] | ai <- toRows a ]
+    where doth u v = sum $ zipWith (*) (toList u) (toList v)
+
+multProp1 (a,b) = a <> b |~| mulH a b
+
+multProp2 (a,b) = ctrans (a <> b) |~| ctrans b <> ctrans a
+
+linearSolveProp f m = f m m |~| ident (rows m)