packages feed

arpack (empty) → 0.1.0.0

raw patch · 12 files changed

+830/−0 lines, 12 filesdep +QuickCheckdep +arpackdep +basesetup-changed

Dependencies added: QuickCheck, arpack, base, concurrent-extra, containers, control-monad-loop, data-default, hmatrix, hspec, ieee754, storable-complex, transformers, vector, vector-algorithms

Files

+ LICENSE view
@@ -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.
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ arpack.cabal view
@@ -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
+ src/Arpack/Exceptions.hs view
@@ -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
+ src/Arpack/Foreign.hs view
@@ -0,0 +1,5 @@+module Arpack.Foreign ( Arpack(..) ) where++import Arpack.Foreign.Class+import Arpack.Foreign.Complex ()+import Arpack.Foreign.Real ()
+ src/Arpack/Foreign/Class.hs view
@@ -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)
+ src/Arpack/Foreign/Complex.hs view
@@ -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)
+ src/Arpack/Foreign/Real.hs view
@@ -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)
+ src/Arpack/Lock.hs view
@@ -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
+ src/Arpack/Options.hs view
@@ -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)+                         }
+ src/Numeric/LinearAlgebra/Arnoldi.hs view
@@ -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)
+ tests/tests.hs view
@@ -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)