packages feed

hmatrix 0.4.0.0 → 0.5.0.1

raw patch · 28 files changed

+565/−1048 lines, 28 files

Files

+ INSTALL view
@@ -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.
README view
@@ -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.
examples/benchmarks.hs view
@@ -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
+ examples/tests.hs view
@@ -0,0 +1,3 @@+import Numeric.LinearAlgebra.Tests++main = runTests 20
hmatrix.cabal view
@@ -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
lib/Data/Packed/Convert.hs view
@@ -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
lib/Data/Packed/Internal/Common.hs view
@@ -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
lib/Data/Packed/Internal/Matrix.hs view
@@ -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)
lib/Data/Packed/Internal/auxi.c view
@@ -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");+}+//---------------------------------------
lib/Data/Packed/Internal/auxi.h view
@@ -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();
lib/Data/Packed/ST.hs view
@@ -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)
lib/Graphics/Plot.hs view
@@ -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
− lib/Numeric/GSL/Matrix.hs 
@@ -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
lib/Numeric/GSL/Special/Internal.hsc view
@@ -1,4 +1,4 @@-{-# OPTIONS -ffi #-}+ {-# LANGUAGE ForeignFunctionInterface #-} ----------------------------------------------------------------------------- {- | Module      :  Numeric.GSL.Special.Internal
lib/Numeric/GSL/Vector.hs view
@@ -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
lib/Numeric/GSL/gsl-aux.c view
@@ -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);
lib/Numeric/GSL/gsl-aux.h view
@@ -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,
lib/Numeric/LinearAlgebra.hs view
@@ -16,7 +16,6 @@     module Data.Packed,     module Numeric.LinearAlgebra.Linear,     module Numeric.LinearAlgebra.Algorithms,-    module Numeric.LinearAlgebra.Instances,     module Numeric.LinearAlgebra.Interface ) where
lib/Numeric/LinearAlgebra/Algorithms.hs view
@@ -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
lib/Numeric/LinearAlgebra/Instances.hs view
@@ -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+
lib/Numeric/LinearAlgebra/Interface.hs view
@@ -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
lib/Numeric/LinearAlgebra/LAPACK.hs view
@@ -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
lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c view
@@ -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 }
lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h view
@@ -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));
lib/Numeric/LinearAlgebra/Linear.hs view
@@ -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
lib/Numeric/LinearAlgebra/Tests.hs view
@@ -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"
lib/Numeric/LinearAlgebra/Tests/Instances.hs view
@@ -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)
lib/Numeric/LinearAlgebra/Tests/Properties.hs view
@@ -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)