diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,30 @@ +Copyright (c) 2013, Thomas Tuegel + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. + + * Neither the name of Thomas Tuegel nor the names of other + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/arpack.cabal b/arpack.cabal new file mode 100644 --- /dev/null +++ b/arpack.cabal @@ -0,0 +1,61 @@ +name: arpack +version: 0.1.0.0 +synopsis: Solve large scale eigenvalue problems +description: + Bindings to ARPACK, a Fortran package implementing the implicitly restarted + Arnoldi (Lanczos) method for solving large scale eigenvalue problems. +license: BSD3 +license-file: LICENSE +author: Thomas Tuegel +maintainer: ttuegel@gmail.com +copyright: (c) 2016 Thomas Tuegel +category: Science +build-type: Simple +cabal-version: >=1.10 + +source-repository head + Type: git + Location: https://github.com/ttuegel/arpack.git + +library + exposed-modules: + Numeric.LinearAlgebra.Arnoldi + other-modules: + Arpack.Exceptions, + Arpack.Foreign, + Arpack.Foreign.Class, + Arpack.Foreign.Complex, + Arpack.Foreign.Real, + Arpack.Lock, + Arpack.Options + build-depends: + base >= 4.6 && < 5, + concurrent-extra >= 0.7, + containers >= 0.5, + control-monad-loop >= 0.1, + data-default >= 0.5, + hmatrix >= 0.16, + ieee754 >= 0.7, + storable-complex >= 0.2, + transformers >= 0.3, + vector >= 0.10, + vector-algorithms >= 0.6 + hs-source-dirs: src + default-language: Haskell2010 + ghc-options: -Wall + pkgconfig-depends: arpack + +test-suite tests + type: exitcode-stdio-1.0 + main-is: tests.hs + hs-source-dirs: tests + build-depends: + arpack + , base + , hmatrix + , hspec >= 2.2 + , QuickCheck >= 2.8 + , vector + , vector-algorithms >= 0.7 + default-language: Haskell2010 + ghc-options: -Wall diff --git a/src/Arpack/Exceptions.hs b/src/Arpack/Exceptions.hs new file mode 100644 --- /dev/null +++ b/src/Arpack/Exceptions.hs @@ -0,0 +1,31 @@ +{-# LANGUAGE DeriveDataTypeable #-} + +module Arpack.Exceptions where + +import Control.Exception +import Data.Int (Int32) +import Data.Typeable + +data MaxIterations = MaxIterations deriving (Show, Typeable) + +instance Exception MaxIterations + +data NoShifts = NoShifts deriving (Show, Typeable) + +instance Exception NoShifts + +data Reallocate = Reallocate deriving (Show, Typeable) + +instance Exception Reallocate + +data XYAUPD = XYAUPD Int32 deriving (Show, Typeable) + +instance Exception XYAUPD + +data XYEUPD = XYEUPD Int32 deriving (Show, Typeable) + +instance Exception XYEUPD + +data Unimplemented = Unimplemented Int32 deriving (Show, Typeable) + +instance Exception Unimplemented diff --git a/src/Arpack/Foreign.hs b/src/Arpack/Foreign.hs new file mode 100644 --- /dev/null +++ b/src/Arpack/Foreign.hs @@ -0,0 +1,5 @@ +module Arpack.Foreign ( Arpack(..) ) where + +import Arpack.Foreign.Class +import Arpack.Foreign.Complex () +import Arpack.Foreign.Real () diff --git a/src/Arpack/Foreign/Class.hs b/src/Arpack/Foreign/Class.hs new file mode 100644 --- /dev/null +++ b/src/Arpack/Foreign/Class.hs @@ -0,0 +1,12 @@ +{-# LANGUAGE RecordWildCards #-} + +module Arpack.Foreign.Class where + +import Data.Vector.Storable.Mutable (IOVector) +import Numeric.LinearAlgebra (Matrix, Vector) + +import Arpack.Options + +class Arpack t where + arpack :: Options t -> Int -> (IOVector t -> IOVector t -> IO ()) + -> IO (Vector t, Matrix t) diff --git a/src/Arpack/Foreign/Complex.hs b/src/Arpack/Foreign/Complex.hs new file mode 100644 --- /dev/null +++ b/src/Arpack/Foreign/Complex.hs @@ -0,0 +1,276 @@ +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE ForeignFunctionInterface #-} +{-# LANGUAGE LambdaCase #-} +{-# LANGUAGE RecordWildCards #-} +{-# LANGUAGE TypeFamilies #-} + +{-# OPTIONS_GHC -fno-warn-orphans #-} + +module Arpack.Foreign.Complex () where + +import qualified Control.Concurrent.Lock as Lock +import Control.Exception (bracket, throwIO) +import Data.Complex +import Data.Maybe (fromMaybe) +import qualified Data.Vector.Storable as VS +import Data.Vector.Storable.Mutable (IOVector) +import qualified Data.Vector.Storable.Mutable as VSM +import Foreign +import Foreign.C.String +import Foreign.C.Types (CChar) +import Foreign.Storable.Complex () +import qualified Numeric.LinearAlgebra.Devel as Dense + +import Arpack.Exceptions +import Arpack.Foreign.Class +import qualified Arpack.Lock as Arpack +import Arpack.Options + +-- +-- * Foreign functions +-- + +foreign import ccall unsafe "znaupd_" + znaupd_ :: Ptr Int32 -- ido + -> Ptr CChar -- bmat + -> Ptr Int32 -- n + -> Ptr CChar -- which + -> Ptr Int32 -- nev + -> Ptr Double -- tol + -> Ptr (Complex Double) -- resid + -> Ptr Int32 -- ncv + -> Ptr (Complex Double) -- v + -> Ptr Int32 -- ldv + -> Ptr Int32 -- iparam + -> Ptr Int32 -- ipntr + -> Ptr (Complex Double) -- workd + -> Ptr (Complex Double) -- workl + -> Ptr Int32 -- lworkl + -> Ptr Double -- rwork + -> Ptr Int32 -- info + -> IO () + +foreign import ccall unsafe "zneupd_" + zneupd_ :: Ptr Int32 -- rvec + -> Ptr CChar -- howmny + -> Ptr Int32 -- select + -> Ptr (Complex Double) -- d + -> Ptr (Complex Double) -- z + -> Ptr Int32 -- ldz + -> Ptr (Complex Double) -- sigma + -> Ptr (Complex Double) -- workev + -> Ptr CChar -- bmat + -> Ptr Int32 -- n + -> Ptr CChar -- which + -> Ptr Int32 -- nev + -> Ptr Double -- tol + -> Ptr (Complex Double) -- resid + -> Ptr Int32 -- ncv + -> Ptr (Complex Double) -- v + -> Ptr Int32 -- ldv + -> Ptr Int32 -- iparam + -> Ptr Int32 -- ipntr + -> Ptr (Complex Double) -- workd + -> Ptr (Complex Double) -- workl + -> Ptr Int32 -- lworkl + -> Ptr Double -- rwork + -> Ptr Int32 -- info + -> IO () + +-- +-- * Wrappers +-- + + +data AUPD + = AUPD + { ido :: {-# UNPACK #-} !(Ptr Int32) + , bmat :: {-# UNPACK #-} !CString + , n :: {-# UNPACK #-} !(Ptr Int32) + , which :: {-# UNPACK #-} !CString + , nev :: {-# UNPACK #-} !(Ptr Int32) + , tol :: {-# UNPACK #-} !(Ptr Double) + , resid :: {-# UNPACK #-} !(IOVector (Complex Double)) + , ncv :: {-# UNPACK #-} !(Ptr Int32) + , v :: {-# UNPACK #-} !(IOVector (Complex Double)) + , ldv :: {-# UNPACK #-} !(Ptr Int32) + , iparam :: {-# UNPACK #-} !(IOVector Int32) + , ipntr :: {-# UNPACK #-} !(IOVector Int32) + , workd :: {-# UNPACK #-} !(IOVector (Complex Double)) + , workl :: {-# UNPACK #-} !(IOVector (Complex Double)) + , lworkl :: {-# UNPACK #-} !(Ptr Int32) + , rwork :: {-# UNPACK #-} !(IOVector Double) + , info :: {-# UNPACK #-} !(Ptr Int32) + } + +withAUPD :: Options (Complex Double) -> Int -> (AUPD -> IO a) -> IO a +withAUPD options dim = bracket initAUPD freeAUPD where + + initAUPD = do + let + _nev = number options + + -- Largest number of basis vectors to use. + -- Work per iteration is O(dim * ncv ^ 2). + + _ncv = min dim (4 * _nev) + + _lworkl = 3 * _ncv * _ncv + 5 * _ncv + + ido <- new 0 + bmat <- newCString "I" + n <- new (fromIntegral dim) + which <- newCString "SR" + nev <- new (fromIntegral _nev) + tol <- new 0 + resid <- VSM.new dim + ncv <- new (fromIntegral _ncv) + v <- VSM.new (dim * _ncv) + ldv <- new (fromIntegral dim) + + iparam <- VSM.new 11 + -- shift strategy + VSM.write iparam (1 - 1) 1 + -- maximum number of iterations + VSM.write iparam (3 - 1) + (fromIntegral (fromMaybe (3 * dim) (maxIterations options))) + -- block size + VSM.write iparam (4 - 1) 1 + -- eigenproblem type + VSM.write iparam (7 - 1) 1 + + ipntr <- VSM.new 14 + workd <- VSM.new (3 * dim) + workl <- VSM.new _lworkl + lworkl <- new (fromIntegral _lworkl) + rwork <- VSM.new _ncv + info <- new 0 + + pure AUPD {..} + + freeAUPD AUPD {..} = do + free ido + free bmat + free n + free which + free nev + free tol + free ncv + free ldv + free lworkl + free info + +aupd :: AUPD -> IO () +aupd (AUPD {..}) = + VSM.unsafeWith resid $ \_resid -> + VSM.unsafeWith v $ \_v -> + VSM.unsafeWith iparam $ \_iparam -> + VSM.unsafeWith ipntr $ \_ipntr -> + VSM.unsafeWith workd $ \_workd -> + VSM.unsafeWith workl $ \_workl -> + VSM.unsafeWith rwork $ \_rwork -> + + znaupd_ ido bmat n which nev tol _resid ncv _v ldv _iparam _ipntr + _workd _workl lworkl _rwork info + +data EUPD + = EUPD + { rvec :: {-# UNPACK #-} !(Ptr Int32) + , howmny :: {-# UNPACK #-} !CString + , select :: {-# UNPACK #-} !(IOVector Int32) + , d :: {-# UNPACK #-} !(IOVector (Complex Double)) + , z :: {-# UNPACK #-} !(IOVector (Complex Double)) + , ldz :: {-# UNPACK #-} !(Ptr Int32) + , sigma :: {-# UNPACK #-} !(Ptr (Complex Double)) + , workev :: {-# UNPACK #-} !(IOVector (Complex Double)) + } + +withEUPD :: AUPD -> (EUPD -> IO a) -> IO a +withEUPD (AUPD {..}) = bracket initEUPD freeEUPD where + + initEUPD = do + _nev <- fromIntegral <$> peek nev + _ncv <- fromIntegral <$> peek ncv + dim <- fromIntegral <$> peek n + + rvec <- new 1 + howmny <- newCString "A" + select <- VSM.new _ncv + d <- VSM.new (_nev + 1) + z <- VSM.new (dim * _nev) + ldz <- new (fromIntegral dim) + sigma <- malloc + workev <- VSM.new (2 * _ncv) + pure EUPD {..} + + freeEUPD (EUPD {..}) = do + free rvec + free howmny + free ldz + free sigma + +eupd :: EUPD -> AUPD -> IO () +eupd (EUPD {..}) (AUPD {..}) = + VSM.unsafeWith select $ \_select -> + VSM.unsafeWith d $ \_d -> + VSM.unsafeWith z $ \_z -> + VSM.unsafeWith workev $ \_workev -> + VSM.unsafeWith resid $ \_resid -> + VSM.unsafeWith v $ \_v -> + VSM.unsafeWith iparam $ \_iparam -> + VSM.unsafeWith ipntr $ \_ipntr -> + VSM.unsafeWith workd $ \_workd -> + VSM.unsafeWith workl $ \_workl -> + VSM.unsafeWith rwork $ \_rwork -> + + zneupd_ rvec howmny _select _d _z ldz sigma _workev + bmat n which nev tol _resid ncv _v ldv _iparam _ipntr + _workd _workl lworkl _rwork info + +instance Arpack (Complex Double) where + + arpack !opts !dim !multiply + -- These variables are all banged because we need to be strict + -- in them _before_ we enter the locked segment of code! If we + -- wait until we're inside the lock, and evaluating one of these + -- variables invokes 'arpack' again, the program will deadlock! + = withAUPD opts dim $ \stateA@(AUPD {..}) -> do + + let + loop = do + aupd stateA + peek ido >>= \case + 99 -> do + peek info >>= \case + 0 -> pure () + 1 -> throwIO MaxIterations + 3 -> throwIO NoShifts + i -> throwIO (XYAUPD i) + + i | abs i == 1 -> do + xi <- fromIntegral <$> VSM.read ipntr 0 + let + x = VSM.slice (xi - 1) dim workd + yi <- fromIntegral <$> VSM.read ipntr 1 + let + y = VSM.slice (yi - 1) dim workd + multiply y x + loop + + | otherwise -> throwIO (Unimplemented i) + + extract = withEUPD stateA $ \stateE@(EUPD {..}) -> do + eupd stateE stateA + + peek info >>= \case + 0 -> pure () + 1 -> throwIO Reallocate + i -> throwIO (XYEUPD i) + + evals <- VS.unsafeFreeze (VSM.slice 0 (number opts) d) + evecs <- VS.unsafeFreeze (VSM.slice 0 (number opts * dim) z) + let matrixFromVector = Dense.matrixFromVector Dense.ColumnMajor + pure (evals, matrixFromVector dim (number opts) evecs) + + Lock.with Arpack.lock (loop >> extract) diff --git a/src/Arpack/Foreign/Real.hs b/src/Arpack/Foreign/Real.hs new file mode 100644 --- /dev/null +++ b/src/Arpack/Foreign/Real.hs @@ -0,0 +1,279 @@ +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE ForeignFunctionInterface #-} +{-# LANGUAGE LambdaCase #-} +{-# LANGUAGE RecordWildCards #-} +{-# LANGUAGE TypeFamilies #-} + +{-# OPTIONS_GHC -fno-warn-orphans #-} + +module Arpack.Foreign.Real () where + +import qualified Control.Concurrent.Lock as Lock +import Control.Exception (bracket, throwIO) +import Data.Maybe (fromMaybe) +import qualified Data.Vector.Storable as VS +import Data.Vector.Storable.Mutable (IOVector) +import qualified Data.Vector.Storable.Mutable as VSM +import Foreign +import Foreign.C.String +import Foreign.C.Types (CChar) +import Foreign.Storable.Complex () +import qualified Numeric.LinearAlgebra.Devel as Dense + +import Arpack.Exceptions +import Arpack.Foreign.Class +import qualified Arpack.Lock as Arpack +import Arpack.Options + +-- +-- * Foreign functions +-- + +foreign import ccall unsafe "dnaupd_" + dnaupd_ :: Ptr Int32 -- ido + -> Ptr CChar -- bmat + -> Ptr Int32 -- n + -> Ptr CChar -- which + -> Ptr Int32 -- nev + -> Ptr Double -- tol + -> Ptr Double -- resid + -> Ptr Int32 -- ncv + -> Ptr Double -- v + -> Ptr Int32 -- ldv + -> Ptr Int32 -- iparam + -> Ptr Int32 -- ipntr + -> Ptr Double -- workd + -> Ptr Double -- workl + -> Ptr Int32 -- lworkl + -> Ptr Int32 -- info + -> IO () + +foreign import ccall unsafe "dneupd_" + dneupd_ :: Ptr Int32 -- rvec + -> Ptr CChar -- howmny + -> Ptr Int32 -- select + -> Ptr Double -- dr + -> Ptr Double -- di + -> Ptr Double -- z + -> Ptr Int32 -- ldz + -> Ptr Double -- sigmar + -> Ptr Double -- sigmai + -> Ptr Double -- workev + -> CString -- bmat + -> Ptr Int32 -- n + -> Ptr CChar -- which + -> Ptr Int32 -- nev + -> Ptr Double -- tol + -> Ptr Double -- resid + -> Ptr Int32 -- ncv + -> Ptr Double -- v + -> Ptr Int32 -- ldv + -> Ptr Int32 -- iparam + -> Ptr Int32 -- ipntr + -> Ptr Double -- workd + -> Ptr Double -- workl + -> Ptr Int32 -- lworkl + -> Ptr Int32 -- info + -> IO () + +-- +-- * Types +-- + +data AUPD + = AUPD + { ido :: {-# UNPACK #-} !(Ptr Int32) + , bmat :: {-# UNPACK #-} !(Ptr CChar) + , n :: {-# UNPACK #-} !(Ptr Int32) + , which :: {-# UNPACK #-} !(Ptr CChar) + , nev :: {-# UNPACK #-} !(Ptr Int32) + , tol :: {-# UNPACK #-} !(Ptr Double) + , resid :: {-# UNPACK #-} !(IOVector Double) + , ncv :: {-# UNPACK #-} !(Ptr Int32) + , v :: {-# UNPACK #-} !(IOVector Double) + , ldv :: {-# UNPACK #-} !(Ptr Int32) + , iparam :: {-# UNPACK #-} !(IOVector Int32) + , ipntr :: {-# UNPACK #-} !(IOVector Int32) + , workd :: {-# UNPACK #-} !(IOVector Double) + , workl :: {-# UNPACK #-} !(IOVector Double) + , lworkl :: {-# UNPACK #-} !(Ptr Int32) + , info :: {-# UNPACK #-} !(Ptr Int32) + } + +withAUPD :: Options Double -> Int -> (AUPD -> IO a) -> IO a +withAUPD options dim = bracket initAUPD freeAUPD where + + initAUPD = do + let + _nev = number options + + -- Largest number of basis vectors to use. + -- Work per iteration is O(dim * ncv ^ 2). + + _ncv = min dim (4 * _nev) + + _lworkl = 3 * _ncv * _ncv + 6 * _ncv + + ido <- new 0 + bmat <- newCString "I" + n <- new (fromIntegral dim) + which <- newCString "SR" + nev <- new (fromIntegral _nev) + tol <- new 0 + resid <- VSM.new dim + ncv <- new (fromIntegral _ncv) + v <- VSM.new (dim * _ncv) + ldv <- new (fromIntegral dim) + + iparam <- VSM.new 11 + -- shift strategy + VSM.write iparam (1 - 1) 1 + -- maximum number of iterations + VSM.write iparam (3 - 1) + (fromIntegral (fromMaybe (3 * dim) (maxIterations options))) + -- block size + VSM.write iparam (4 - 1) 1 + -- eigenproblem type + VSM.write iparam (7 - 1) 1 + + ipntr <- VSM.new 14 + workd <- VSM.new (3 * dim) + workl <- VSM.new _lworkl + lworkl <- new (fromIntegral _lworkl) + info <- new 0 + + pure AUPD {..} + + freeAUPD AUPD {..} = do + free ido + free bmat + free n + free which + free nev + free tol + free ncv + free ldv + free lworkl + free info + +data EUPD + = EUPD + { rvec :: {-# UNPACK #-} !(Ptr Int32) + , howmny :: {-# UNPACK #-} !(Ptr CChar) + , select :: {-# UNPACK #-} !(IOVector Int32) + , dr :: {-# UNPACK #-} !(IOVector Double) + , di :: {-# UNPACK #-} !(IOVector Double) + , z :: {-# UNPACK #-} !(IOVector Double) + , ldz :: {-# UNPACK #-} !(Ptr Int32) + , sigmar :: {-# UNPACK #-} !(Ptr Double) + , sigmai :: {-# UNPACK #-} !(Ptr Double) + , workev :: {-# UNPACK #-} !(IOVector Double) + } + +withEUPD :: AUPD -> (EUPD -> IO a) -> IO a +withEUPD (AUPD {..}) = bracket initEUPD freeEUPD where + + initEUPD = do + _nev <- fromIntegral <$> peek nev + _ncv <- fromIntegral <$> peek ncv + dim <- fromIntegral <$> peek n + + rvec <- new 1 + howmny <- newCString "A" + select <- VSM.new _ncv + dr <- VSM.new (_nev + 1) + di <- VSM.new (_nev + 1) + z <- VSM.new (dim * _nev) + ldz <- new (fromIntegral dim) + sigmar <- malloc + sigmai <- malloc + workev <- VSM.new (2 * _ncv) + pure EUPD {..} + + freeEUPD (EUPD {..}) = do + free rvec + free howmny + free ldz + free sigmar + free sigmai + +-- +-- * Wrappers +-- + +aupd :: AUPD -> IO () +aupd (AUPD {..}) = + VSM.unsafeWith resid $ \_resid -> + VSM.unsafeWith v $ \_v -> + VSM.unsafeWith iparam $ \_iparam -> + VSM.unsafeWith ipntr $ \_ipntr -> + VSM.unsafeWith workd $ \_workd -> + VSM.unsafeWith workl $ \_workl -> + + dnaupd_ ido bmat n which nev tol _resid ncv _v ldv + _iparam _ipntr _workd _workl lworkl info + +eupd :: EUPD -> AUPD -> IO () +eupd (EUPD {..}) (AUPD {..}) = + VSM.unsafeWith select $ \_select -> + VSM.unsafeWith dr $ \_dr -> + VSM.unsafeWith di $ \_di -> + VSM.unsafeWith z $ \_z -> + VSM.unsafeWith workev $ \_workev -> + VSM.unsafeWith resid $ \_resid -> + VSM.unsafeWith v $ \_v -> + VSM.unsafeWith iparam $ \_iparam -> + VSM.unsafeWith ipntr $ \_ipntr -> + VSM.unsafeWith workd $ \_workd -> + VSM.unsafeWith workl $ \_workl -> + + dneupd_ rvec howmny _select _dr _di _z ldz sigmar sigmai _workev + bmat n which nev tol _resid ncv _v ldv _iparam _ipntr + _workd _workl lworkl info + +instance Arpack Double where + + arpack !opts !dim !multiply + -- These variables are all banged because we need to be strict + -- in them _before_ we enter the locked segment of code! If we + -- wait until we're inside the lock, and evaluating one of these + -- variables invokes 'arpack' again, the program will deadlock! + = withAUPD opts dim $ \stateA@(AUPD {..}) -> do + + let + loop = do + aupd stateA + peek ido >>= \case + 99 -> do + peek info >>= \case + 0 -> pure () + 1 -> throwIO MaxIterations + 3 -> throwIO NoShifts + i -> throwIO (XYAUPD i) + + i | abs i == 1 -> do + xi <- fromIntegral <$> VSM.read ipntr 0 + let + x = VSM.slice (xi - 1) dim workd + yi <- fromIntegral <$> VSM.read ipntr 1 + let + y = VSM.slice (yi - 1) dim workd + multiply y x + loop + + | otherwise -> throwIO (Unimplemented i) + + extract = withEUPD stateA $ \stateE@(EUPD {..}) -> do + eupd stateE stateA + + peek info >>= \case + 0 -> pure () + 1 -> throwIO Reallocate + i -> throwIO (XYEUPD i) + + evals <- VS.unsafeFreeze (VSM.slice 0 (number opts) dr) + evecs <- VS.unsafeFreeze (VSM.slice 0 (number opts * dim) z) + let matrixFromVector = Dense.matrixFromVector Dense.ColumnMajor + pure (evals, matrixFromVector dim (number opts) evecs) + + Lock.with Arpack.lock (loop >> extract) diff --git a/src/Arpack/Lock.hs b/src/Arpack/Lock.hs new file mode 100644 --- /dev/null +++ b/src/Arpack/Lock.hs @@ -0,0 +1,11 @@ +{-# OPTIONS_GHC -fno-cse -fno-full-laziness #-} -- `unsafePerformIO` + +module Arpack.Lock ( lock ) where + +import Control.Concurrent.Lock (Lock) +import qualified Control.Concurrent.Lock as Lock +import System.IO.Unsafe (unsafePerformIO) + +lock :: Lock +{-# NOINLINE lock #-} +lock = unsafePerformIO $ Lock.new diff --git a/src/Arpack/Options.hs b/src/Arpack/Options.hs new file mode 100644 --- /dev/null +++ b/src/Arpack/Options.hs @@ -0,0 +1,27 @@ +{-# LANGUAGE GADTs #-} +{-# LANGUAGE KindSignatures #-} + +module Arpack.Options where + +import Data.Complex + +data Which :: * -> * where + LM :: Which a + SM :: Which a + LR :: Which a + SR :: Which a + LI :: Which (Complex a) + SI :: Which (Complex a) + +showWhich :: Which a -> String +showWhich LM = "LM" +showWhich SM = "SM" +showWhich LR = "LR" +showWhich SR = "SR" +showWhich LI = "LI" +showWhich SI = "SI" + +data Options a = Options { which :: !(Which a) + , number :: {-# UNPACK #-} !Int + , maxIterations :: !(Maybe Int) + } diff --git a/src/Numeric/LinearAlgebra/Arnoldi.hs b/src/Numeric/LinearAlgebra/Arnoldi.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/LinearAlgebra/Arnoldi.hs @@ -0,0 +1,30 @@ +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE LambdaCase #-} +{-# LANGUAGE RecordWildCards #-} + +module Numeric.LinearAlgebra.Arnoldi + ( Arpack, eig + -- * Options + , Options(..), Which(..) + -- * Exceptions + , MaxIterations(..) + , NoShifts(..) + , Reallocate(..) + , XYAUPD(..) + , XYEUPD(..) + , Unimplemented(..) + ) where + +import Data.Vector.Storable (Vector) +import Data.Vector.Storable.Mutable (IOVector) +import Numeric.LinearAlgebra (Matrix) +import System.IO.Unsafe (unsafePerformIO) + +import Arpack.Exceptions +import Arpack.Foreign +import Arpack.Options + +eig :: Arpack t => Options t -> Int -> (IOVector t -> IOVector t -> IO ()) + -> (Vector t, Matrix t) +eig !options !dim !multiply = unsafePerformIO (arpack options dim multiply) diff --git a/tests/tests.hs b/tests/tests.hs new file mode 100644 --- /dev/null +++ b/tests/tests.hs @@ -0,0 +1,66 @@ +module Main where + +import Data.Ord (comparing) +import qualified Data.Vector.Algorithms.Intro as Intro +import qualified Data.Vector.Storable as VS +import Data.Vector.Storable.Mutable (IOVector) +import Numeric.LinearAlgebra hiding (eig, vector) +import Numeric.LinearAlgebra.Arnoldi +import Test.Hspec +import Test.QuickCheck + +main :: IO () +main = hspec $ do + describe "Numeric.LinearAlgebra.Arnoldi.eig" $ do + + it "returns the diagonal of an arbitrary diagonal matrix" $ property $ do + dim <- arbitrary `suchThat` (> 3) + diagonal <- VS.fromList <$> vector dim + let + _matrix = diag diagonal :: Matrix Double + nev = dim - 3 + options = Options { which = SR + , number = nev + , maxIterations = Nothing + } + sort = VS.modify Intro.sort + actual = (sort . fst) (eig options dim (multiply _matrix)) + expected = (VS.take nev . sort) diagonal + relative = VS.maximum (VS.zipWith (%) expected actual) + pure (counterexample (show (expected, actual)) (relative < 1E-4)) + + it "returns the diagonal of an arbitrary diagonal matrix" $ property $ do + dim <- arbitrary `suchThat` (> 3) + diagonal <- VS.fromList <$> vector dim + let + _matrix = diag diagonal :: Matrix (Complex Double) + nev = dim - 3 + options = Options { which = SR + , number = nev + , maxIterations = Nothing + } + sort = VS.modify (Intro.sortBy (comparing realPart)) + actual = (sort . fst) (eig options dim (multiply _matrix)) + expected = (VS.take nev . sort) diagonal + relative = VS.maximum (VS.zipWith (%%) expected actual) + pure (counterexample (show (expected, actual)) (relative < 1E-4)) + +multiply :: Numeric a => Matrix a -> IOVector a -> IOVector a -> IO () +multiply _matrix dst src = do + x <- VS.freeze src + let y = _matrix #> x + VS.copy dst y + +(%) :: Double -> Double -> Double +(%) a b = + let a_plus_b = a + b + in if a_plus_b /= 0 + then abs ((a - b) / a_plus_b) + else a - b + +(%%) :: Complex Double -> Complex Double -> Double +(%%) a b = + let a_plus_b = a + b + in if a_plus_b /= 0 + then magnitude ((a - b) / a_plus_b) + else magnitude (a - b)