comfort-blas (empty) → 0.0
raw patch · 29 files changed
+4816/−0 lines, 29 filesdep +QuickCheckdep +basedep +blas-ffisetup-changed
Dependencies added: QuickCheck, base, blas-ffi, comfort-array, comfort-blas, containers, deepseq, doctest-exitcode-stdio, doctest-lib, guarded-allocation, netlib-ffi, transformers, utility-ht
Files
- LICENSE +27/−0
- Makefile +23/−0
- Setup.lhs +3/−0
- comfort-blas.cabal +87/−0
- src/Numeric/BLAS/Matrix/Layout.hs +21/−0
- src/Numeric/BLAS/Matrix/Modifier.hs +49/−0
- src/Numeric/BLAS/Matrix/RowMajor.hs +180/−0
- src/Numeric/BLAS/Private.hs +262/−0
- src/Numeric/BLAS/Scalar.hs +158/−0
- src/Numeric/BLAS/Slice.hs +153/−0
- src/Numeric/BLAS/Vector.hs +801/−0
- src/Numeric/BLAS/Vector/Slice.hs +903/−0
- test/Main.hs +32/−0
- test/Test/ComplexDouble/Numeric/BLAS/Vector.hs +201/−0
- test/Test/ComplexDouble/Numeric/BLAS/Vector/Slice.hs +271/−0
- test/Test/ComplexDouble/Type.hs +5/−0
- test/Test/ComplexFloat/Numeric/BLAS/Vector.hs +201/−0
- test/Test/ComplexFloat/Numeric/BLAS/Vector/Slice.hs +271/−0
- test/Test/ComplexFloat/Type.hs +5/−0
- test/Test/Double/Numeric/BLAS/Vector.hs +201/−0
- test/Test/Double/Numeric/BLAS/Vector/Slice.hs +271/−0
- test/Test/Double/Type.hs +3/−0
- test/Test/Float/Numeric/BLAS/Vector.hs +201/−0
- test/Test/Float/Numeric/BLAS/Vector/Slice.hs +271/−0
- test/Test/Float/Type.hs +3/−0
- test/Test/Generator.hs +28/−0
- test/Test/Numeric/BLAS/Slice.hs +53/−0
- test/Test/Slice.hs +113/−0
- test/Test/Utility.hs +19/−0
+ LICENSE view
@@ -0,0 +1,27 @@+Copyright (c) Henning Thielemann 2023++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions+are met:+1. Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.+2. 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.+3. Neither the name of the author nor the names of his contributors+ may be used to endorse or promote products derived from this software+ without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE REGENTS 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 AUTHORS 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.
+ Makefile view
@@ -0,0 +1,23 @@+TYPES = Float Double ComplexFloat ComplexDouble+TEST_MODULES = test/Test/Numeric/BLAS/Slice.hs \+ $(patsubst %, test/Test/%/Numeric/BLAS/Vector.hs, $(TYPES)) \+ $(patsubst %, test/Test/%/Numeric/BLAS/Vector/Slice.hs, $(TYPES))++run-test: update-test+ runhaskell Setup configure --user --enable-tests+ runhaskell Setup build+ runhaskell Setup haddock+ runhaskell Setup test comfort-blas-test --show-details=streaming++update-test: $(TEST_MODULES)++test/Test/Numeric/BLAS/Slice.hs: src/Numeric/BLAS/Slice.hs+ doctest-extract-0.1 -i src/ -o test/ --module-prefix Test Numeric.BLAS.Slice++test/Test/%/Numeric/BLAS/Vector.hs: src/Numeric/BLAS/Vector.hs+ doctest-extract-0.1 -i src/ -o test/ --module-prefix Test.$* Numeric.BLAS.Vector+ perl -i -p -e "s:NumberModule:$*:" $@++test/Test/%/Numeric/BLAS/Vector/Slice.hs: src/Numeric/BLAS/Vector/Slice.hs+ doctest-extract-0.1 -i src/ -o test/ --module-prefix Test.$* Numeric.BLAS.Vector.Slice+ perl -i -p -e "s:NumberModule:$*:" $@
+ Setup.lhs view
@@ -0,0 +1,3 @@+#! /usr/bin/env runhaskell+> import Distribution.Simple+> main = defaultMain
+ comfort-blas.cabal view
@@ -0,0 +1,87 @@+Cabal-Version: 2.2+Name: comfort-blas+Version: 0.0+License: BSD-3-Clause+License-File: LICENSE+Author: Henning Thielemann <haskell@henning-thielemann.de>+Maintainer: Henning Thielemann <haskell@henning-thielemann.de>+Homepage: https://hub.darcs.net/thielema/comfort-blas/+Category: Math+Synopsis: Numerical Basic Linear Algebra using BLAS+Description:+ This is a high-level interface to BLAS.+ It provides support for working on slices of arrays.+Tested-With: GHC==7.4.2, GHC==7.8.4, GHC==8.6.5, GHC==9.4.5+Build-Type: Simple+Extra-Source-Files:+ Makefile++Source-Repository this+ Tag: 0.0+ Type: darcs+ Location: https://hub.darcs.net/thielema/comfort-blas/++Source-Repository head+ Type: darcs+ Location: https://hub.darcs.net/thielema/comfort-blas/++Library+ Build-Depends:+ blas-ffi >=0.0 && <0.2,+ netlib-ffi >=0.1.1 && <0.2,+ comfort-array >=0.5.2 && <0.6,+ guarded-allocation >=0.0 && <0.1,+ containers >=0.4 && <0.7,+ deepseq >=1.3 && <1.5,+ transformers >=0.4 && <0.7,+ utility-ht >=0.0.13 && <0.1,+ base >=4.5 && <5++ Default-Language: Haskell98+ GHC-Options: -Wall+ Hs-Source-Dirs: src+ Exposed-Modules:+ Numeric.BLAS.Matrix.RowMajor+ Numeric.BLAS.Matrix.Modifier+ Numeric.BLAS.Matrix.Layout+ Numeric.BLAS.Slice+ Numeric.BLAS.Vector.Slice+ Numeric.BLAS.Vector+ Numeric.BLAS.Scalar+ Other-Modules:+ Numeric.BLAS.Private++Test-Suite comfort-blas-test+ Type: exitcode-stdio-1.0+ Build-Depends:+ comfort-blas,+ netlib-ffi,+ comfort-array,+ doctest-exitcode-stdio >=0.0 && <0.1,+ doctest-lib >=0.1 && <0.2,+ QuickCheck >=2.5 && <3,+ containers,+ utility-ht,+ base++ Default-Language: Haskell98+ GHC-Options: -Wall+ Hs-Source-Dirs: test+ Main-Is: Main.hs+ Other-Modules:+ Test.Float.Numeric.BLAS.Vector.Slice+ Test.Double.Numeric.BLAS.Vector.Slice+ Test.ComplexFloat.Numeric.BLAS.Vector.Slice+ Test.ComplexDouble.Numeric.BLAS.Vector.Slice+ Test.Float.Numeric.BLAS.Vector+ Test.Double.Numeric.BLAS.Vector+ Test.ComplexFloat.Numeric.BLAS.Vector+ Test.ComplexDouble.Numeric.BLAS.Vector+ Test.Float.Type+ Test.Double.Type+ Test.ComplexFloat.Type+ Test.ComplexDouble.Type+ Test.Numeric.BLAS.Slice+ Test.Slice+ Test.Generator+ Test.Utility
+ src/Numeric/BLAS/Matrix/Layout.hs view
@@ -0,0 +1,21 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}+module Numeric.BLAS.Matrix.Layout where++import Control.DeepSeq (NFData, rnf)+++data Order = RowMajor | ColumnMajor+ deriving (Eq, Show)++instance NFData Order where+ rnf RowMajor = ()+ rnf ColumnMajor = ()++flipOrder :: Order -> Order+flipOrder RowMajor = ColumnMajor+flipOrder ColumnMajor = RowMajor++transposeFromOrder :: Order -> Char+transposeFromOrder RowMajor = 'T'+transposeFromOrder ColumnMajor = 'N'
+ src/Numeric/BLAS/Matrix/Modifier.hs view
@@ -0,0 +1,49 @@+module Numeric.BLAS.Matrix.Modifier where++import Numeric.BLAS.Matrix.Layout (Order(RowMajor,ColumnMajor), flipOrder)++import Data.Monoid (Monoid, mempty, mappend)+import Data.Semigroup (Semigroup, (<>))++++data Transposition = NonTransposed | Transposed+ deriving (Eq, Show, Enum, Bounded)++instance Semigroup Transposition where+ x<>y = if x==y then NonTransposed else Transposed++instance Monoid Transposition where+ mempty = NonTransposed+ mappend = (<>)++transposeOrder :: Transposition -> Order -> Order+transposeOrder NonTransposed = id+transposeOrder Transposed = flipOrder++++data Conjugation = NonConjugated | Conjugated+ deriving (Eq, Show, Enum, Bounded)++instance Semigroup Conjugation where+ x<>y = if x==y then NonConjugated else Conjugated++instance Monoid Conjugation where+ mempty = NonConjugated+ mappend = (<>)++conjugatedOnRowMajor :: Order -> Conjugation+conjugatedOnRowMajor RowMajor = Conjugated+conjugatedOnRowMajor ColumnMajor = NonConjugated+++data Inversion = NonInverted | Inverted+ deriving (Eq, Show, Enum, Bounded)++instance Semigroup Inversion where+ x<>y = if x==y then NonInverted else Inverted++instance Monoid Inversion where+ mempty = NonInverted+ mappend = (<>)
+ src/Numeric/BLAS/Matrix/RowMajor.hs view
@@ -0,0 +1,180 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}+module Numeric.BLAS.Matrix.RowMajor (+ Matrix,+ Vector,+ takeRow,+ takeColumn,+ fromRows,+ tensorProduct,+ decomplex,+ recomplex,+ scaleRows,+ scaleColumns,+ ) where++import qualified Numeric.BLAS.Private as Private+import Numeric.BLAS.Matrix.Modifier (Conjugation(NonConjugated,Conjugated))+import Numeric.BLAS.Scalar (zero, one)+import Numeric.BLAS.Private (ShapeInt, shapeInt, ComplexShape, pointerSeq)++import qualified Numeric.BLAS.FFI.Generic as Blas+import qualified Numeric.Netlib.Utility as Call+import qualified Numeric.Netlib.Class as Class++import Foreign.Marshal.Array (copyArray, advancePtr)+import Foreign.ForeignPtr (withForeignPtr, castForeignPtr)+import Foreign.Storable (Storable)++import Control.Monad.Trans.Cont (ContT(ContT), evalContT)+import Control.Monad.IO.Class (liftIO)++import qualified Data.Array.Comfort.Storable.Unchecked as Array+import qualified Data.Array.Comfort.Shape as Shape+import Data.Array.Comfort.Storable.Unchecked (Array(Array))++import Data.Complex (Complex)+import Data.Foldable (forM_)+++type Matrix height width = Array (height,width)+type Vector = Array++takeRow ::+ (Shape.Indexed height, Shape.C width, Shape.Index height ~ ix,+ Storable a) =>+ ix -> Matrix height width a -> Vector width a+takeRow ix (Array (height,width) x) =+ Array.unsafeCreateWithSize width $ \n yPtr ->+ withForeignPtr x $ \xPtr ->+ copyArray yPtr (advancePtr xPtr (n * Shape.offset height ix)) n++takeColumn ::+ (Shape.C height, Shape.Indexed width, Shape.Index width ~ ix,+ Class.Floating a) =>+ ix -> Matrix height width a -> Vector height a+takeColumn ix (Array (height,width) x) =+ Array.unsafeCreateWithSize height $ \n yPtr -> evalContT $ do+ let offset = Shape.offset width ix+ nPtr <- Call.cint n+ xPtr <- ContT $ withForeignPtr x+ incxPtr <- Call.cint $ Shape.size width+ incyPtr <- Call.cint 1+ liftIO $ Blas.copy nPtr (advancePtr xPtr offset) incxPtr yPtr incyPtr+++fromRows ::+ (Shape.C width, Eq width, Storable a) =>+ width -> [Vector width a] -> Matrix ShapeInt width a+fromRows width rows =+ Array.unsafeCreate (shapeInt $ length rows, width) $ \dstPtr ->+ let widthSize = Shape.size width+ in forM_ (zip (pointerSeq widthSize dstPtr) rows) $+ \(dstRowPtr, Array.Array rowWidth srcFPtr) ->+ withForeignPtr srcFPtr $ \srcPtr -> do+ Call.assert+ "Matrix.fromRows: non-matching vector size"+ (width == rowWidth)+ copyArray dstRowPtr srcPtr widthSize+++-- ToDo: use lapack:Private.multiplyMatrix+tensorProduct ::+ (Shape.C height, Shape.C width, Class.Floating a) =>+ Either Conjugation Conjugation ->+ Vector height a -> Vector width a -> Matrix height width a+tensorProduct side (Array height x) (Array width y) =+ Array.unsafeCreate (height,width) $ \cPtr -> do+ let m = Shape.size width+ let n = Shape.size height+ let trans conjugated =+ case conjugated of NonConjugated -> 'T'; Conjugated -> 'C'+ let ((transa,transb),(lda,ldb)) =+ case side of+ Left c -> ((trans c, 'N'),(1,1))+ Right c -> (('N', trans c),(m,n))+ evalContT $ do+ transaPtr <- Call.char transa+ transbPtr <- Call.char transb+ mPtr <- Call.cint m+ nPtr <- Call.cint n+ kPtr <- Call.cint 1+ alphaPtr <- Call.number one+ aPtr <- ContT $ withForeignPtr y+ ldaPtr <- Call.leadingDim lda+ bPtr <- ContT $ withForeignPtr x+ ldbPtr <- Call.leadingDim ldb+ betaPtr <- Call.number zero+ ldcPtr <- Call.leadingDim m+ liftIO $+ Blas.gemm+ transaPtr transbPtr mPtr nPtr kPtr alphaPtr+ aPtr ldaPtr bPtr ldbPtr betaPtr cPtr ldcPtr+++decomplex ::+ (Class.Real a) =>+ Matrix height width (Complex a) ->+ Matrix height (width, ComplexShape) a+decomplex (Array (height,width) a) =+ Array (height, (width, Shape.static)) (castForeignPtr a)++recomplex ::+ (Class.Real a) =>+ Matrix height (width, ComplexShape) a ->+ Matrix height width (Complex a)+recomplex (Array (height, (width, Shape.NestedTuple _)) a) =+ Array (height,width) (castForeignPtr a)+++scaleRows ::+ (Shape.C height, Eq height, Shape.C width, Class.Floating a) =>+ Vector height a -> Matrix height width a -> Matrix height width a+scaleRows (Array heightX x) (Array shape@(height,width) a) =+ Array.unsafeCreate shape $ \bPtr -> do+ Call.assert "scaleRows: sizes mismatch" (heightX == height)+ evalContT $ do+ let m = Shape.size height+ let n = Shape.size width+ nPtr <- Call.cint n+ xPtr <- ContT $ withForeignPtr x+ aPtr <- ContT $ withForeignPtr a+ incaPtr <- Call.cint 1+ incbPtr <- Call.cint 1+ liftIO $ sequence_ $ take m $+ zipWith3+ (\xkPtr akPtr bkPtr -> do+ Blas.copy nPtr akPtr incaPtr bkPtr incbPtr+ Blas.scal nPtr xkPtr bkPtr incbPtr)+ (pointerSeq 1 xPtr)+ (pointerSeq n aPtr)+ (pointerSeq n bPtr)++scaleColumns ::+ (Shape.C height, Shape.C width, Eq width, Class.Floating a) =>+ Vector width a -> Matrix height width a -> Matrix height width a+scaleColumns (Array widthX x) (Array shape@(height,width) a) =+ Array.unsafeCreate shape $ \bPtr -> do+ Call.assert "scaleColumns: sizes mismatch" (widthX == width)+ evalContT $ do+ let m = Shape.size height+ let n = Shape.size width+ transPtr <- Call.char 'N'+ nPtr <- Call.cint n+ klPtr <- Call.cint 0+ kuPtr <- Call.cint 0+ alphaPtr <- Call.number one+ xPtr <- ContT $ withForeignPtr x+ ldxPtr <- Call.leadingDim 1+ aPtr <- ContT $ withForeignPtr a+ incaPtr <- Call.cint 1+ betaPtr <- Call.number zero+ incbPtr <- Call.cint 1+ liftIO $ sequence_ $ take m $+ zipWith+ (\akPtr bkPtr ->+ Private.gbmv transPtr+ nPtr nPtr klPtr kuPtr alphaPtr xPtr ldxPtr+ akPtr incaPtr betaPtr bkPtr incbPtr)+ (pointerSeq n aPtr)+ (pointerSeq n bPtr)
+ src/Numeric/BLAS/Private.hs view
@@ -0,0 +1,262 @@+{-# LANGUAGE TypeFamilies #-}+module Numeric.BLAS.Private where++import qualified Numeric.BLAS.FFI.Real as BlasReal+import qualified Numeric.BLAS.FFI.Generic as Blas+import qualified Numeric.Netlib.Utility as Call+import qualified Numeric.Netlib.Class as Class+import Numeric.BLAS.Matrix.Modifier (Conjugation(NonConjugated, Conjugated))+import Numeric.BLAS.Scalar (RealOf, zero, one, minusOne, isZero)++import qualified Foreign.Marshal.Array.Guarded as ForeignArray+import qualified Foreign.Marshal.Utils as Marshal+import qualified Foreign.C.String as CStr+import Foreign.Marshal.Array (advancePtr)+import Foreign.C.Types (CChar, CInt)+import Foreign.Ptr (Ptr, castPtr)+import Foreign.Storable (Storable, peek, pokeElemOff, peekElemOff)++import Control.Monad.Trans.Cont (evalContT)+import Control.Monad.IO.Class (liftIO)+import Control.Monad (when)+import Control.Applicative (liftA2)++import qualified Data.Array.Comfort.Shape as Shape++import qualified Data.Complex as Complex+import Data.Complex (Complex((:+)))++import Prelude hiding (sum)+++type ShapeInt = Shape.ZeroBased Int++shapeInt :: Int -> ShapeInt+shapeInt = Shape.ZeroBased+++realPtr :: Ptr a -> Ptr (RealOf a)+realPtr = castPtr+++pointerSeq :: (Storable a) => Int -> Ptr a -> [Ptr a]+pointerSeq k ptr = iterate (flip advancePtr k) ptr+++fill :: (Class.Floating a) => a -> Int -> Ptr a -> IO ()+fill a n dstPtr = evalContT $ do+ nPtr <- Call.cint n+ srcPtr <- Call.number a+ incxPtr <- Call.cint 0+ incyPtr <- Call.cint 1+ liftIO $ Blas.copy nPtr srcPtr incxPtr dstPtr incyPtr+++copyConjugate ::+ (Class.Floating a) =>+ Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()+copyConjugate nPtr xPtr incxPtr yPtr incyPtr = do+ Blas.copy nPtr xPtr incxPtr yPtr incyPtr+ lacgv nPtr yPtr incyPtr++++newtype Sum a = Sum {runSum :: Int -> Ptr a -> Int -> IO a}++sum :: Class.Floating a => Int -> Ptr a -> Int -> IO a+sum =+ runSum $+ Class.switchFloating+ (Sum sumReal)+ (Sum sumReal)+ (Sum sumComplex)+ (Sum sumComplex)++sumReal :: Class.Real a => Int -> Ptr a -> Int -> IO a+sumReal n xPtr incx =+ evalContT $ do+ nPtr <- Call.cint n+ incxPtr <- Call.cint incx+ yPtr <- Call.real one+ incyPtr <- Call.cint 0+ liftIO $ BlasReal.dot nPtr xPtr incxPtr yPtr incyPtr++sumComplex, sumComplexAlt ::+ Class.Real a => Int -> Ptr (Complex a) -> Int -> IO (Complex a)+sumComplex n xPtr incx =+ evalContT $ do+ nPtr <- Call.cint n+ let sxPtr = realPtr xPtr+ incxPtr <- Call.cint (2*incx)+ yPtr <- Call.real one+ incyPtr <- Call.cint 0+ liftIO $+ liftA2 (Complex.:+)+ (BlasReal.dot nPtr sxPtr incxPtr yPtr incyPtr)+ (BlasReal.dot nPtr (advancePtr sxPtr 1) incxPtr yPtr incyPtr)++sumComplexAlt n aPtr inca =+ evalContT $ do+ transPtr <- Call.char 'N'+ mPtr <- Call.cint 2+ nPtr <- Call.cint n+ onePtr <- Call.number one+ inc0Ptr <- Call.cint 0+ let saPtr = realPtr aPtr+ ldaPtr <- Call.leadingDim (2*inca)+ sxPtr <- Call.allocaArray n+ incxPtr <- Call.cint 1+ betaPtr <- Call.number zero+ yPtr <- Call.alloca+ let syPtr = realPtr yPtr+ incyPtr <- Call.cint 1+ liftIO $ do+ Blas.copy nPtr onePtr inc0Ptr sxPtr incxPtr+ gemv+ transPtr mPtr nPtr onePtr saPtr ldaPtr+ sxPtr incxPtr betaPtr syPtr incyPtr+ peek yPtr+++mul ::+ (Class.Floating a) =>+ Conjugation -> Int -> Ptr a -> Int -> Ptr a -> Int -> Ptr a -> Int -> IO ()+mul conj n aPtr inca xPtr incx yPtr incy =+ mulAdd conj n aPtr inca xPtr incx zero yPtr incy++mulAdd ::+ (Class.Floating a) =>+ Conjugation ->+ Int -> Ptr a -> Int -> Ptr a -> Int -> a -> Ptr a -> Int -> IO ()+mulAdd conj n aPtr inca xPtr incx beta yPtr incy = evalContT $ do+ transPtr <- Call.char $ case conj of NonConjugated -> 'N'; Conjugated -> 'C'+ nPtr <- Call.cint n+ klPtr <- Call.cint 0+ kuPtr <- Call.cint 0+ alphaPtr <- Call.number one+ ldaPtr <- Call.leadingDim inca+ incxPtr <- Call.cint incx+ betaPtr <- Call.number beta+ incyPtr <- Call.cint incy+ liftIO $+ Blas.gbmv transPtr+ nPtr nPtr klPtr kuPtr alphaPtr aPtr ldaPtr+ xPtr incxPtr betaPtr yPtr incyPtr++{- |+Use the foldBalanced trick.+-}+product :: (Class.Floating a) => Int -> Ptr a -> Int -> IO a+product n aPtr inca =+ case compare n 1 of+ LT -> return one+ EQ -> peek aPtr+ GT -> let n2 = div n 2; new = n-n2+ in ForeignArray.alloca (2*new-1) $ \xPtr -> do+ mulPairs n2 aPtr inca xPtr 1+ when (odd n) $ pokeElemOff xPtr n2 =<< peekElemOff aPtr ((n-1)*inca)+ productLoop new xPtr++{- |+If 'mul' would be based on a scalar loop+we would not need to cut the vector into chunks.++The invariance is:+When calling @productLoop n xPtr@,+starting from xPtr there is storage allocated for 2*n-1 elements.+-}+productLoop :: (Class.Floating a) => Int -> Ptr a -> IO a+productLoop n xPtr =+ if n==1+ then peek xPtr+ else do+ let n2 = div n 2+ mulPairs n2 xPtr 1 (advancePtr xPtr n) 1+ productLoop (n-n2) (advancePtr xPtr (2*n2))++mulPairs ::+ (Class.Floating a) =>+ Int -> Ptr a -> Int -> Ptr a -> Int -> IO ()+mulPairs n aPtr inca xPtr incx =+ let inca2 = 2*inca+ in mul NonConjugated n aPtr inca2 (advancePtr aPtr inca) inca2 xPtr incx+++newtype LACGV a = LACGV {getLACGV :: Ptr CInt -> Ptr a -> Ptr CInt -> IO ()}++lacgv :: Class.Floating a => Ptr CInt -> Ptr a -> Ptr CInt -> IO ()+lacgv =+ getLACGV $+ Class.switchFloating+ (LACGV $ const $ const $ const $ return ())+ (LACGV $ const $ const $ const $ return ())+ (LACGV clacgv)+ (LACGV clacgv)++clacgv :: Class.Real a => Ptr CInt -> Ptr (Complex a) -> Ptr CInt -> IO ()+clacgv nPtr xPtr incxPtr =+ Marshal.with minusOne $ \saPtr -> do+ incx <- peek incxPtr+ Marshal.with (2*incx) $ \incyPtr ->+ BlasReal.scal nPtr saPtr (advancePtr (realPtr xPtr) 1) incyPtr+++{-+Work around an inconsistency of BLAS.+In case of a zero-column matrix+BLAS's gemv and gbmv do not initialize the target vector.+In contrast, these work-arounds do.+-}+{-# INLINE gemv #-}+gemv ::+ (Class.Floating a) =>+ Ptr CChar -> Ptr CInt -> Ptr CInt ->+ Ptr a -> Ptr a -> Ptr CInt ->+ Ptr a -> Ptr CInt -> Ptr a -> Ptr a -> Ptr CInt -> IO ()+gemv transPtr mPtr nPtr+ alphaPtr aPtr ldaPtr xPtr incxPtr betaPtr yPtr incyPtr = do+ initializeMV transPtr mPtr nPtr betaPtr yPtr incyPtr+ Blas.gemv transPtr mPtr nPtr+ alphaPtr aPtr ldaPtr xPtr incxPtr betaPtr yPtr incyPtr++{-# INLINE gbmv #-}+gbmv ::+ (Class.Floating a) =>+ Ptr CChar -> Ptr CInt -> Ptr CInt -> Ptr CInt -> Ptr CInt ->+ Ptr a -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt ->+ Ptr a -> Ptr a -> Ptr CInt -> IO ()+gbmv transPtr mPtr nPtr klPtr kuPtr+ alphaPtr aPtr ldaPtr xPtr incxPtr betaPtr yPtr incyPtr = do+ initializeMV transPtr mPtr nPtr betaPtr yPtr incyPtr+ Blas.gbmv transPtr mPtr nPtr klPtr kuPtr+ alphaPtr aPtr ldaPtr xPtr incxPtr betaPtr yPtr incyPtr++initializeMV ::+ Class.Floating a =>+ Ptr CChar -> Ptr CInt -> Ptr CInt -> Ptr a -> Ptr a -> Ptr CInt -> IO ()+initializeMV transPtr mPtr nPtr betaPtr yPtr incyPtr = do+ trans <- peek transPtr+ let (mtPtr,ntPtr) =+ if trans == CStr.castCharToCChar 'N'+ then (mPtr,nPtr) else (nPtr,mPtr)+ n <- peek ntPtr+ beta <- peek betaPtr+ when (n == 0 && isZero beta) $+ Marshal.with 0 $ \incbPtr ->+ Blas.copy mtPtr betaPtr incbPtr yPtr incyPtr+++{-+ToDo:++type ComplexShape =+ Shape.NestedTuple Shape.TupleAccessor (Complex Shape.Element)++This would allow the use of Complex.realPart as accessor,+but it requires GHC>7.6.3 or so, where realPart has no RealFloat constraint.+-}+type ComplexShape = Shape.NestedTuple Shape.TupleIndex (Complex Shape.Element)++ixReal, ixImaginary :: Shape.ElementIndex (Complex Shape.Element)+ixReal :+ ixImaginary =+ Shape.indexTupleFromShape (Shape.static :: ComplexShape)
+ src/Numeric/BLAS/Scalar.hs view
@@ -0,0 +1,158 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE GADTs #-}+module Numeric.BLAS.Scalar (+ RealOf,+ ComplexOf,+ zero,+ one,+ minusOne,+ isZero,++ selectReal,+ selectFloating,+ ComplexSingleton(Real,Complex),+ complexSingleton,+ complexSingletonOf,+ complexSingletonOfFunctor,++ PrecisionSingleton(Float,Double),+ precisionSingleton,+ precisionOf,+ precisionOfFunctor,++ equal,+ fromReal,+ toComplex,+ absolute,+ absoluteSquared,+ norm1,+ realPart,+ conjugate,+ ) where++import qualified Numeric.Netlib.Class as Class++import Data.Functor.Identity (Identity(Identity, runIdentity))++import qualified Data.Complex as Complex+import Data.Complex (Complex((:+)))+++type family RealOf x++type instance RealOf Float = Float+type instance RealOf Double = Double+type instance RealOf (Complex.Complex a) = a+++type ComplexOf x = Complex.Complex (RealOf x)+++data ComplexSingleton a where+ Real :: (Class.Real a, RealOf a ~ a) => ComplexSingleton a+ Complex :: (Class.Real a) => ComplexSingleton (Complex.Complex a)++complexSingleton :: (Class.Floating a) => ComplexSingleton a+complexSingleton = Class.switchFloating Real Real Complex Complex++complexSingletonOf :: (Class.Floating a) => a -> ComplexSingleton a+complexSingletonOf = const complexSingleton++complexSingletonOfFunctor :: (Class.Floating a) => f a -> ComplexSingleton a+complexSingletonOfFunctor = const complexSingleton++withComplexSingleton :: (Class.Floating a) => (ComplexSingleton a -> a) -> a+withComplexSingleton f = f complexSingleton+++data PrecisionSingleton a where+ Float :: PrecisionSingleton Float+ Double :: PrecisionSingleton Double++precisionSingleton :: (Class.Real a) => PrecisionSingleton a+precisionSingleton = Class.switchReal Float Double++precisionOf :: (Class.Real a) => a -> PrecisionSingleton a+precisionOf _ = precisionSingleton++precisionOfFunctor :: (Class.Real a) => f a -> PrecisionSingleton a+precisionOfFunctor _ = precisionSingleton+++-- ToDo: move to netlib-carray:Utility or netlib-ffi:Class+zero, one, minusOne :: Class.Floating a => a+zero = selectFloating 0 0 0 0+one = selectFloating 1 1 1 1+minusOne = selectFloating (-1) (-1) (-1) (-1)++selectReal :: (Class.Real a) => Float -> Double -> a+selectReal rf rd = runIdentity $ Class.switchReal (Identity rf) (Identity rd)++selectFloating ::+ (Class.Floating a) =>+ Float -> Double -> Complex.Complex Float -> Complex.Complex Double -> a+selectFloating rf rd cf cd =+ withComplexSingleton $ \sw ->+ case sw of+ Real -> selectReal rf rd+ Complex -> Class.switchReal cf cd++++equal :: (Class.Floating a) => a -> a -> Bool+equal a b =+ case complexSingletonOf a of+ Real -> a==b+ Complex -> a==b+++isZero :: (Class.Floating a) => a -> Bool+isZero = equal zero+++fromReal :: (Class.Floating a) => RealOf a -> a+fromReal a =+ withComplexSingleton $ \sw ->+ case sw of+ Real -> a+ Complex -> a:+0++toComplex :: (Class.Floating a) => a -> ComplexOf a+toComplex a =+ case complexSingletonOf a of+ Real -> a:+0+ Complex -> a++realPart :: (Class.Floating a) => a -> RealOf a+realPart a =+ case complexSingletonOf a of+ Real -> a+ Complex -> Complex.realPart a++absolute :: (Class.Floating a) => a -> RealOf a+absolute a =+ case complexSingletonOf a of+ Real -> abs a+ Complex -> Complex.magnitude a+++norm1 :: (Class.Floating a) => a -> RealOf a+norm1 a =+ case complexSingletonOf a of+ Real -> abs a+ Complex -> case a of r:+i -> abs r + abs i+++absoluteSquared :: (Class.Floating a) => a -> RealOf a+absoluteSquared a =+ case complexSingletonOf a of+ Real -> a*a+ Complex -> case a of r:+i -> r*r+i*i+++conjugate :: (Class.Floating a) => a -> a+conjugate a =+ case complexSingletonOf a of+ Real -> a+ Complex -> Complex.conjugate a
+ src/Numeric/BLAS/Slice.hs view
@@ -0,0 +1,153 @@+{-# LANGUAGE TypeOperators #-}+module Numeric.BLAS.Slice where++import qualified Data.Array.Comfort.Shape as Shape+import qualified Data.Array.Comfort.Boxed as BoxedArray+import Data.Array.Comfort.Shape ((::+)((::+)))++import qualified Data.Traversable as Trav+import qualified Data.List as List+import Data.Map (Map)+++{- $setup+>>> import qualified Numeric.BLAS.Slice as Slice+>>> import Test.Slice (shapeInt)+>>>+>>> import qualified Data.Array.Comfort.Boxed as Array+>>> import qualified Data.Array.Comfort.Shape as Shape+>>> import qualified Data.Map as Map+>>> import Data.Array.Comfort.Shape ((::+)((::+)))+>>> import Data.Array.Comfort.Boxed ((!))+>>>+>>> import Control.Applicative (liftA3, pure)+>>>+>>> import qualified Test.QuickCheck as QC+>>>+>>> genSlice :: sh -> QC.Gen (Slice.T sh)+>>> genSlice sh =+>>> liftA3 Slice.Cons (QC.choose (0,100)) (QC.choose (1,100)) (pure sh)+>>>+>>> genSlice2 :: shA -> shB -> QC.Gen (Slice.T shA, Slice.T shB)+>>> genSlice2 shA shB = do+>>> s <- QC.choose (0,100)+>>> k <- QC.choose (1,100)+>>> return (Slice.Cons s k shA, Slice.Cons s k shB)+-}+++data T sh = Cons {start, skip :: Int, shape :: sh}+ deriving (Eq, Show)++fromShape :: (Shape.C sh) => sh -> T sh+fromShape = Cons 0 1+++row ::+ (Shape.Indexed sh0, Shape.C sh1) => Shape.Index sh0 -> T (sh0,sh1) -> T sh1+row ix0 (Cons s k (sh0,sh1)) =+ Cons (s + k * Shape.offset sh0 ix0 * Shape.size sh1) k sh1++column ::+ (Shape.C sh0, Shape.Indexed sh1) => Shape.Index sh1 -> T (sh0,sh1) -> T sh0+column ix1 (Cons s k (sh0,sh1)) =+ let (size1, offset1) = Shape.sizeOffset sh1+ in Cons (s + k * offset1 ix1) (k * size1) sh0++{- |+prop> QC.forAll (QC.choose (1,100)) $ \numRows -> QC.forAll (QC.choose (0,100)) $ \numColumns -> QC.forAll (genSlice (shapeInt numRows, shapeInt numColumns)) $ \slice -> QC.forAll (QC.elements $ Shape.indices $ shapeInt numRows) $ \ix -> Slice.row ix slice == Slice.rowArray slice ! ix+-}+rowArray ::+ (Shape.Indexed sh0, Shape.C sh1) =>+ T (sh0,sh1) -> BoxedArray.Array sh0 (T sh1)+rowArray (Cons s k (sh0,sh1)) =+ let step = Shape.size sh1 * k+ in BoxedArray.fromList sh0 $+ List.map (\si -> Cons si k sh1) $+ List.take (Shape.size sh0) $ iterate (step+) s++{- |+prop> QC.forAll (QC.choose (0,100)) $ \numRows -> QC.forAll (QC.choose (1,100)) $ \numColumns -> QC.forAll (genSlice (shapeInt numRows, shapeInt numColumns)) $ \slice -> QC.forAll (QC.elements $ Shape.indices $ shapeInt numColumns) $ \ix -> Slice.column ix slice == Slice.columnArray slice ! ix+-}+columnArray ::+ (Shape.C sh0, Shape.Indexed sh1) =>+ T (sh0,sh1) -> BoxedArray.Array sh1 (T sh0)+columnArray (Cons s k (sh0,sh1)) =+ let step = Shape.size sh1 * k+ in BoxedArray.fromList sh1 $+ List.map (\si -> Cons si step sh0) $+ List.take (Shape.size sh1) $ iterate (k+) s+++topSubmatrix ::+ (Shape.C sh, Shape.C sh0, Shape.C sh1) =>+ T (sh0::+sh1, sh) -> T (sh0,sh)+topSubmatrix (Cons s k (sh0::+_sh1, sh)) =+ Cons s k (sh0,sh)++bottomSubmatrix ::+ (Shape.C sh, Shape.C sh0, Shape.C sh1) =>+ T (sh0::+sh1, sh) -> T (sh1,sh)+bottomSubmatrix (Cons s k (sh0::+sh1, sh)) =+ Cons (s + k * Shape.size sh0) k (sh1,sh)+++cartesianFromSquare :: T (Shape.Square sh) -> T (sh,sh)+cartesianFromSquare (Cons s k (Shape.Square sh)) = Cons s k (sh,sh)++squareRow ::+ (Shape.Indexed sh) => Shape.Index sh -> T (Shape.Square sh) -> T sh+squareRow ix0 = row ix0 . cartesianFromSquare++squareColumn ::+ (Shape.Indexed sh) => Shape.Index sh -> T (Shape.Square sh) -> T sh+squareColumn ix1 = column ix1 . cartesianFromSquare++++plane12 ::+ (Shape.Indexed sh0, Shape.C sh1, Shape.C sh2) =>+ Shape.Index sh0 -> T (sh0,sh1,sh2) -> T (sh1,sh2)+plane12 ix0 (Cons s k (sh0,sh1,sh2)) =+ Cons (s + k * Shape.offset sh0 ix0 * Shape.size (sh1,sh2)) k (sh1,sh2)++plane01 ::+ (Shape.C sh0, Shape.C sh1, Shape.Indexed sh2) =>+ Shape.Index sh2 -> T (sh0,sh1,sh2) -> T (sh0,sh1)+plane01 ix2 (Cons s k (sh0,sh1,sh2)) =+ let (size2, offset2) = Shape.sizeOffset sh2+ in Cons (s + k * offset2 ix2) (k * size2) (sh0,sh1)++column2of3 ::+ (Shape.Indexed sh0, Shape.Indexed sh1, Shape.C sh2) =>+ Shape.Index sh0 -> Shape.Index sh1 -> T (sh0,sh1,sh2) -> T sh2+column2of3 ix0 ix1 = row ix1 . plane12 ix0++column1of3 ::+ (Shape.Indexed sh0, Shape.C sh1, Shape.Indexed sh2) =>+ Shape.Index sh0 -> Shape.Index sh2 -> T (sh0,sh1,sh2) -> T sh1+column1of3 ix0 ix2 = column ix2 . plane12 ix0++column0of3 ::+ (Shape.C sh0, Shape.Indexed sh1, Shape.Indexed sh2) =>+ Shape.Index sh1 -> Shape.Index sh2 -> T (sh0,sh1,sh2) -> T sh0+column0of3 ix1 ix2 = column ix1 . plane01 ix2+++left :: (Shape.C sh0, Shape.C sh1) => T (sh0::+sh1) -> T sh0+left (Cons s k (sh0::+_sh1)) = Cons s k sh0++right :: (Shape.C sh0, Shape.C sh1) => T (sh0::+sh1) -> T sh1+right (Cons s k (sh0::+sh1)) = Cons (s + k * Shape.size sh0) k sh1+++{- |+prop> QC.forAll (fmap shapeInt $ QC.choose (0,100)) $ \shapeA -> QC.forAll (fmap shapeInt $ QC.choose (0,100)) $ \shapeB -> QC.forAll (fmap shapeInt $ QC.choose (0,100)) $ \shapeC -> QC.forAll (genSlice2 (Map.fromList $ ('a', shapeA) : ('b', shapeB) : ('c', shapeC) : []) (shapeA ::+ shapeB ::+ shapeC)) $ \(sliceMap, sliceParted) -> Slice.map sliceMap Map.! 'b' == Slice.left (Slice.right sliceParted)++prop> QC.forAll (QC.choose (0,100)) $ \numRows -> QC.forAll (QC.choose (0,100)) $ \numColumns -> let rowShape = shapeInt numRows; columnShape = shapeInt numColumns; mapShape = Map.fromList $ map (\k -> (k, columnShape)) (Shape.indices rowShape) in QC.forAll (genSlice2 mapShape (rowShape, columnShape)) $ \(sliceMap, sliceMatrix) -> Map.toAscList (Slice.map sliceMap) == Array.toAssociations (Slice.rowArray sliceMatrix)+-}+map :: (Shape.C sh) => T (Map k sh) -> Map k (T sh)+map (Cons s k m) =+ snd $+ Trav.mapAccumL+ (\offset sh -> (offset + Shape.size sh * k, Cons offset k sh)) s m
+ src/Numeric/BLAS/Vector.hs view
@@ -0,0 +1,801 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}+module Numeric.BLAS.Vector (+ Vector,+ RealOf,+ ComplexOf,+ toList,+ fromList,+ autoFromList,+ CheckedArray.append, (+++),+ CheckedArray.take, CheckedArray.drop,+ CheckedArray.takeLeft, CheckedArray.takeRight,+ swap,+ CheckedArray.singleton,+ constant,+ zero,+ one,+ unit,+ dot, inner, (-*|),+ sum,+ absSum,+ norm2,+ norm2Squared,+ normInf,+ normInf1,+ argAbsMaximum,+ argAbs1Maximum,+ product,+ scale, scaleReal, (.*|),+ add, sub, (|+|), (|-|),+ negate, raise,+ mac,+ mul, mulConj,+ minimum, argMinimum,+ maximum, argMaximum,+ limits, argLimits,+ CheckedArray.foldl,+ CheckedArray.foldl1,+ CheckedArray.foldMap,++ conjugate,+ fromReal,+ toComplex,+ realPart,+ imaginaryPart,+ zipComplex,+ unzipComplex,+ ) where++import qualified Numeric.BLAS.Matrix.RowMajor as RowMajor+import qualified Numeric.BLAS.Scalar as Scalar+import qualified Numeric.BLAS.Private as Private+import Numeric.BLAS.Matrix.Modifier (Conjugation(NonConjugated, Conjugated))+import Numeric.BLAS.Scalar (ComplexOf, RealOf, minusOne)+import Numeric.BLAS.Private+ (ComplexShape, ixReal, ixImaginary, fill, copyConjugate, realPtr)++import qualified Numeric.BLAS.FFI.Generic as Blas+import qualified Numeric.BLAS.FFI.Complex as BlasComplex+import qualified Numeric.BLAS.FFI.Real as BlasReal+import qualified Numeric.Netlib.Utility as Call+import qualified Numeric.Netlib.Class as Class++import qualified Foreign.Marshal.Array.Guarded as ForeignArray+import Foreign.Marshal.Array (advancePtr)+import Foreign.ForeignPtr (withForeignPtr, castForeignPtr)+import Foreign.Ptr (Ptr)+import Foreign.Storable (Storable, peek, peekElemOff, pokeElemOff)+import Foreign.C.Types (CInt)++import System.IO.Unsafe (unsafePerformIO)++import Control.Monad.Trans.Cont (ContT(ContT), evalContT)+import Control.Monad.IO.Class (liftIO)+import Control.Monad.ST (runST)+import Control.Monad (fmap, return, (=<<))+import Control.Applicative (liftA3, (<$>))++import qualified Data.Array.Comfort.Storable.Mutable.Unchecked as UMutArray+import qualified Data.Array.Comfort.Storable.Mutable as MutArray+import qualified Data.Array.Comfort.Storable.Unchecked as Array+import qualified Data.Array.Comfort.Storable as CheckedArray+import qualified Data.Array.Comfort.Shape as Shape+import Data.Array.Comfort.Storable.Unchecked (Array(Array), append, (!))+import Data.Array.Comfort.Shape ((::+))++import Data.Function (id, flip, ($), (.))+import Data.Complex (Complex)+import Data.Maybe (Maybe(Nothing,Just), maybe)+import Data.Tuple.HT (mapFst, uncurry3)+import Data.Tuple (fst, snd)+import Data.Ord ((>=))+import Data.Eq (Eq, (==))++import Prelude (Int, fromIntegral, (-), Char, error, IO)+++{- $setup+>>> import Test.NumberModule.Type (Number_)+>>> import Test.Generator (genNumber)+>>> import Test.Utility (approx)+>>> import qualified Numeric.BLAS.Vector as Vector+>>> import qualified Numeric.Netlib.Class as Class+>>> import qualified Data.Array.Comfort.Shape as Shape+>>> import qualified Data.Array.Comfort.Storable as Array+>>> import qualified Data.List as List+>>> import Numeric.BLAS.Vector ((+++), (|+|), (|-|))+>>> import Numeric.BLAS.Scalar (RealOf, absolute, minusOne)+>>> import Data.Array.Comfort.Storable (Array, (!))+>>> import Data.Complex (Complex((:+)))+>>> import Control.Applicative (liftA2)+>>> import Control.Monad (replicateM)+>>>+>>> import qualified Test.QuickCheck as QC+>>> import Test.QuickCheck ((==>))+>>>+>>> type Vector = Vector.Vector (Shape.ZeroBased Int)+>>> type Real_ = RealOf Number_+>>> type Complex_ = Complex Real_+>>>+>>> maxElem :: Integer+>>> maxElem = 10+>>>+>>> maxDim :: Int+>>> maxDim = 100+>>>+>>> genVector :: (Class.Floating a) => Integer -> Int -> QC.Gen (Vector a)+>>> genVector maxE dim =+>>> fmap (Vector.fromList (Shape.ZeroBased dim)) $+>>> replicateM dim $ genNumber maxE+>>>+>>> isNonEmpty :: Shape.C sh => Array sh a -> Bool+>>> isNonEmpty xs = Shape.size (Array.shape xs) > 0+>>>+>>> newtype ArbRealVector = ArbRealVector (Vector Real_)+>>> deriving (Show)+>>>+>>> instance QC.Arbitrary ArbRealVector where+>>> arbitrary = do+>>> dim <- QC.choose (0,maxDim)+>>> fmap ArbRealVector $ genVector maxElem dim+>>> shrink (ArbRealVector xs) =+>>> map (ArbRealVector . Vector.autoFromList) .+>>> QC.shrink . Vector.toList $ xs+>>>+>>> newtype ArbVector = ArbVector (Vector Number_)+>>> deriving (Show)+>>>+>>> instance QC.Arbitrary ArbVector where+>>> arbitrary = do+>>> dim <- QC.choose (0,maxDim)+>>> fmap ArbVector $ genVector maxElem dim+>>> shrink (ArbVector xs) =+>>> map (ArbVector . Vector.autoFromList) .+>>> QC.shrink . Vector.toList $ xs+>>>+>>> data ArbVector2 =+>>> ArbVector2 (Vector Number_) (Vector Number_)+>>> deriving (Show)+>>>+>>> instance QC.Arbitrary ArbVector2 where+>>> arbitrary = do+>>> dim <- QC.choose (0,maxDim)+>>> liftA2 ArbVector2+>>> (genVector maxElem dim) (genVector maxElem dim)+-}+++type Vector = Array+++toList :: (Shape.C sh, Storable a) => Vector sh a -> [a]+toList = Array.toList++fromList :: (Shape.C sh, Storable a) => sh -> [a] -> Vector sh a+fromList = CheckedArray.fromList++autoFromList :: (Storable a) => [a] -> Vector (Shape.ZeroBased Int) a+autoFromList = Array.vectorFromList+++{- |+> constant () = singleton++However, singleton does not need 'Class.Floating' constraint.+-}+constant :: (Shape.C sh, Class.Floating a) => sh -> a -> Vector sh a+constant sh a = Array.unsafeCreateWithSize sh $ fill a++zero :: (Shape.C sh, Class.Floating a) => sh -> Vector sh a+zero = flip constant Scalar.zero++one :: (Shape.C sh, Class.Floating a) => sh -> Vector sh a+one = flip constant Scalar.one++unit ::+ (Shape.Indexed sh, Class.Floating a) =>+ sh -> Shape.Index sh -> Vector sh a+unit sh ix = Array.unsafeCreateWithSize sh $ \n xPtr -> do+ fill Scalar.zero n xPtr+ pokeElemOff xPtr (Shape.offset sh ix) Scalar.one+++{- |+Precedence and associativity (right) of (List.++).+This also matches '(::+)'.+-}+infixr 5 ++++++{- |+prop> \(ArbVector xs) (ArbVector ys) (ArbVector zs) -> Vector.toList ((xs +++ ys) +++ zs) == Vector.toList (xs +++ (ys +++ zs))+-}+(+++) ::+ (Shape.C shx, Shape.C shy, Storable a) =>+ Vector shx a -> Vector shy a -> Vector (shx::+shy) a+(+++) = append+++{- |+prop> QC.forAll (QC.choose (1,100)) $ \dim -> QC.forAll (QC.choose (0, dim-1)) $ \i -> QC.forAll (QC.choose (0, dim-1)) $ \j -> Vector.unit (Shape.ZeroBased dim) i == (Vector.swap i j (Vector.unit (Shape.ZeroBased dim) j) :: Vector Number_)+-}+swap ::+ (Shape.Indexed sh, Storable a) =>+ Shape.Index sh -> Shape.Index sh -> Vector sh a -> Vector sh a+swap i j x =+ runST (do+ y <- MutArray.thaw x+ xi <- MutArray.read y i+ xj <- MutArray.read y j+ MutArray.write y i xj+ MutArray.write y j xi+ UMutArray.unsafeFreeze y)+++infixl 7 -*|, .*|++newtype Dot f a = Dot {runDot :: f a -> f a -> a}++{- |+> dot x y = Matrix.toScalar (singleRow x -*| singleColumn y)+-}+dot, (-*|) ::+ (Shape.C sh, Eq sh, Class.Floating a) =>+ Vector sh a -> Vector sh a -> a+(-*|) = dot+dot =+ runDot $+ Class.switchFloating+ (Dot dotReal)+ (Dot dotReal)+ (Dot $ dotComplex 'T')+ (Dot $ dotComplex 'T')++{- |+prop> \(ArbVector2 xs ys) -> Vector.inner xs ys == Vector.dot (Vector.conjugate xs) ys+-}+inner ::+ (Shape.C sh, Eq sh, Class.Floating a) =>+ Vector sh a -> Vector sh a -> a+inner =+ runDot $+ Class.switchFloating+ (Dot dotReal)+ (Dot dotReal)+ (Dot $ dotComplex 'C')+ (Dot $ dotComplex 'C')++dotReal ::+ (Shape.C sh, Eq sh, Class.Real a) =>+ Vector sh a -> Vector sh a -> a+dotReal arrX@(Array shX _x) (Array shY y) = unsafePerformIO $ do+ Call.assert "dot: shapes mismatch" (shX == shY)+ evalContT $ do+ (nPtr, sxPtr, incxPtr) <- vectorArgs arrX+ syPtr <- ContT $ withForeignPtr y+ incyPtr <- Call.cint 1+ liftIO $ BlasReal.dot nPtr sxPtr incxPtr syPtr incyPtr++{-+We cannot use 'cdot' because Haskell's FFI+does not support Complex numbers as return values.+-}+dotComplex ::+ (Shape.C sh, Eq sh, Class.Real a) =>+ Char -> Vector sh (Complex a) -> Vector sh (Complex a) -> Complex a+dotComplex trans (Array shX x) (Array shY y) = unsafePerformIO $ do+ Call.assert "dot: shapes mismatch" (shX == shY)+ evalContT $ do+ let m = Shape.size shX+ transPtr <- Call.char trans+ mPtr <- Call.cint m+ nPtr <- Call.cint 1+ alphaPtr <- Call.number Scalar.one+ xPtr <- ContT $ withForeignPtr x+ ldxPtr <- Call.leadingDim m+ yPtr <- ContT $ withForeignPtr y+ incyPtr <- Call.cint 1+ betaPtr <- Call.number Scalar.zero+ zPtr <- Call.alloca+ inczPtr <- Call.cint 1+ liftIO $+ Private.gemv+ transPtr mPtr nPtr alphaPtr xPtr ldxPtr+ yPtr incyPtr betaPtr zPtr inczPtr+ liftIO $ peek zPtr++{- |+prop> \(ArbVector xs) -> Vector.sum xs == List.sum (Vector.toList xs)+-}+sum :: (Shape.C sh, Class.Floating a) => Vector sh a -> a+sum (Array sh x) = unsafePerformIO $+ withForeignPtr x $ \xPtr -> Private.sum (Shape.size sh) xPtr 1+++{- |+Sum of the absolute values of real numbers or components of complex numbers.+For real numbers it is equivalent to 'Numeric.LAPACK.Vector.norm1'.+-}+absSum :: (Shape.C sh, Class.Floating a) => Vector sh a -> RealOf a+absSum arr = unsafePerformIO $+ evalContT $ liftIO . uncurry3 asum =<< vectorArgs arr++asum :: Class.Floating a => Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a)+asum =+ getNrm $+ Class.switchFloating+ (Nrm BlasReal.asum) (Nrm BlasReal.asum)+ (Nrm BlasComplex.casum) (Nrm BlasComplex.casum)+++{- |+Euclidean norm of a vector or Frobenius norm of a matrix.+-}+norm2 :: (Shape.C sh, Class.Floating a) => Vector sh a -> RealOf a+norm2 arr = unsafePerformIO $+ evalContT $ liftIO . uncurry3 nrm2 =<< vectorArgs arr++nrm2 :: Class.Floating a => Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a)+nrm2 =+ getNrm $+ Class.switchFloating+ (Nrm BlasReal.nrm2) (Nrm BlasReal.nrm2)+ (Nrm BlasComplex.cnrm2) (Nrm BlasComplex.cnrm2)++newtype Nrm a =+ Nrm {getNrm :: Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a)}+++newtype Norm f a = Norm {getNorm :: f a -> RealOf a}++norm2Squared :: (Shape.C sh, Class.Floating a) => Vector sh a -> RealOf a+norm2Squared =+ getNorm $+ Class.switchFloating+ (Norm norm2SquaredReal)+ (Norm norm2SquaredReal)+ (Norm $ norm2SquaredReal . decomplex)+ (Norm $ norm2SquaredReal . decomplex)++norm2SquaredReal :: (Shape.C sh, Class.Real a) => Vector sh a -> a+norm2SquaredReal arr =+ unsafePerformIO $ evalContT $ do+ (nPtr, sxPtr, incxPtr) <- vectorArgs arr+ liftIO $ BlasReal.dot nPtr sxPtr incxPtr sxPtr incxPtr+++{- |+prop> \(ArbVector xs) -> Vector.normInf xs == List.maximum (0 : List.map absolute (Vector.toList xs))+-}+normInf :: (Shape.C sh, Class.Floating a) => Vector sh a -> RealOf a+normInf arr = unsafePerformIO $+ fmap (Scalar.absolute . maybe Scalar.zero snd) $ absMax arr++{- |+Computes (almost) the infinity norm of the vector.+For complex numbers every element is replaced+by the sum of the absolute component values first.+-}+normInf1 :: (Shape.C sh, Class.Floating a) => Vector sh a -> RealOf a+normInf1 arr = unsafePerformIO $+ evalContT $ do+ (nPtr, sxPtr, incxPtr) <- vectorArgs arr+ liftIO $+ fmap (Scalar.norm1 . maybe Scalar.zero snd) $+ peekElemOff1 sxPtr =<< Blas.iamax nPtr sxPtr incxPtr+++{- |+Returns the index and value of the element with the maximal absolute value.+Caution: It actually returns the value of the element, not its absolute value!++>>> Vector.argAbsMaximum $ Vector.autoFromList [1:+2, 3:+4, 5, 6 :: Complex_]+(3,6.0 :+ 0.0)++prop> \(ArbVector xs) -> isNonEmpty xs ==> let (xi,xm) = Vector.argAbsMaximum xs in xs!xi == xm+prop> \(ArbVector xs) -> isNonEmpty xs ==> let (_xi,xm) = Vector.argAbsMaximum xs in List.all (\x -> absolute x <= absolute xm) $ Vector.toList xs+prop> \(ArbVector xs, ArbVector ys) -> isNonEmpty xs && isNonEmpty ys ==> let (_xi,xm) = Vector.argAbsMaximum xs; (_yi,ym) = Vector.argAbsMaximum ys; (zi,zm) = Vector.argAbsMaximum (xs+++ys) in case zi of Left _ -> xm==zm && absolute xm >= absolute ym; Right _ -> ym==zm && absolute xm < absolute ym+-}+argAbsMaximum ::+ (Shape.InvIndexed sh, Class.Floating a) =>+ Vector sh a -> (Shape.Index sh, a)+argAbsMaximum arr = unsafePerformIO $+ fmap+ (maybe+ (error "Vector.argAbsMaximum: empty vector")+ (mapFst (Shape.uncheckedIndexFromOffset $ Array.shape arr))) $+ absMax arr++absMax ::+ (Shape.C sh, Class.Floating a) =>+ Vector sh a -> IO (Maybe (Int, a))+absMax arr@(Array sh x) =+ case Scalar.complexSingletonOfFunctor arr of+ Scalar.Real -> evalContT $ do+ (nPtr, sxPtr, incxPtr) <- vectorArgs arr+ liftIO $ peekElemOff1 sxPtr =<< Blas.iamax nPtr sxPtr incxPtr+ Scalar.Complex -> evalContT $ do+ let n = Shape.size sh+ sxPtr <- ContT $ withForeignPtr x+ liftIO $ peekElemOff1 sxPtr =<< absMaxComplex n sxPtr++absMaxComplex :: (Class.Real a) => Int -> Ptr (Complex a) -> IO CInt+absMaxComplex n sxPtr =+ ForeignArray.alloca n $ \syPtr -> do+ let xrPtr = realPtr sxPtr+ Private.mul NonConjugated n xrPtr 2 xrPtr 2 syPtr 1+ let xiPtr = advancePtr xrPtr 1+ Private.mulAdd NonConjugated n xiPtr 2 xiPtr 2 Scalar.one syPtr 1+ evalContT $ do+ nPtr <- Call.cint n+ incyPtr <- Call.cint 1+ liftIO $ Blas.iamax nPtr syPtr incyPtr+++{- |+Returns the index and value of the element with the maximal absolute value.+The function does not strictly compare the absolute value of a complex number+but the sum of the absolute complex components.+Caution: It actually returns the value of the element, not its absolute value!++>>> Vector.argAbs1Maximum $ Vector.autoFromList [1:+2, 3:+4, 5, 6 :: Complex_]+(1,3.0 :+ 4.0)++prop> \(ArbRealVector xs) -> isNonEmpty xs ==> Vector.argAbsMaximum xs == Vector.argAbs1Maximum xs+-}+argAbs1Maximum ::+ (Shape.InvIndexed sh, Class.Floating a) =>+ Vector sh a -> (Shape.Index sh, a)+argAbs1Maximum arr = unsafePerformIO $+ evalContT $ do+ (nPtr, sxPtr, incxPtr) <- vectorArgs arr+ liftIO $+ fmap+ (maybe+ (error "Vector.argAbs1Maximum: empty vector")+ (mapFst (Shape.uncheckedIndexFromOffset $ Array.shape arr))) $+ peekElemOff1 sxPtr =<< Blas.iamax nPtr sxPtr incxPtr++vectorArgs ::+ (Shape.C sh) => Array sh a -> ContT r IO (Ptr CInt, Ptr a, Ptr CInt)+vectorArgs (Array sh x) =+ liftA3 (,,)+ (Call.cint $ Shape.size sh)+ (ContT $ withForeignPtr x)+ (Call.cint 1)++peekElemOff1 :: (Storable a) => Ptr a -> CInt -> IO (Maybe (Int, a))+peekElemOff1 ptr k1 =+ let k1i = fromIntegral k1+ ki = k1i-1+ in if k1i == 0+ then return Nothing+ else Just . (,) ki <$> peekElemOff ptr ki+++{- |+prop> QC.forAll (QC.choose (0,10)) $ \dim -> QC.forAll (genVector 3 dim) $ \xs -> approx 1e-2 (Vector.product xs) (List.product (Vector.toList (xs :: Vector Number_)))+-}+product :: (Shape.C sh, Class.Floating a) => Vector sh a -> a+product (Array sh x) = unsafePerformIO $+ withForeignPtr x $ \xPtr -> Private.product (Shape.size sh) xPtr 1+++{- |+For restrictions see 'limits'.++prop> \(ArbRealVector xs) -> isNonEmpty xs ==> Vector.minimum xs == List.minimum (Vector.toList xs)+prop> \(ArbRealVector xs) -> isNonEmpty xs ==> Vector.maximum xs == List.maximum (Vector.toList xs)+prop> \(ArbRealVector xs) -> isNonEmpty xs ==> - Vector.maximum xs == Vector.minimum (Vector.negate xs)+-}+minimum, maximum :: (Shape.C sh, Class.Real a) => Vector sh a -> a+minimum = fst . limits+maximum = snd . limits++{- |+For restrictions see 'limits'.+-}+argMinimum, argMaximum ::+ (Shape.InvIndexed sh, Shape.Index sh ~ ix, Class.Real a) =>+ Vector sh a -> (ix,a)+argMinimum = fst . argLimits+argMaximum = snd . argLimits++{- |+prop> \(ArbRealVector xs) -> isNonEmpty xs ==> Vector.limits xs == Array.limits xs++In contrast to 'Array.limits'+this implementation is based on fast BLAS functions.+It should be faster than @Array.minimum@ and @Array.maximum@+although it is certainly not as fast as possible.+Boths limits share the precision of the limit with the larger absolute value.+This implies for example that you cannot rely on the property+that @raise (- minimum x) x@ has only non-negative elements.+-}+limits :: (Shape.C sh, Class.Real a) => Vector sh a -> (a,a)+limits xs0 =+ let xs = Array.mapShape Shape.Deferred xs0+ x0 = snd $ argAbs1Maximum xs+ x1 = xs ! fst (argAbs1Maximum (raise (-x0) xs))+ in if x0>=0 then (x1,x0) else (x0,x1)++{- |+For restrictions see 'limits'.+-}+argLimits ::+ (Shape.InvIndexed sh, Shape.Index sh ~ ix, Class.Real a) =>+ Vector sh a -> ((ix,a),(ix,a))+argLimits xs =+ let p0@(_i0,x0) = argAbs1Maximum xs+ p1 = (i1,xs!i1); i1 = fst $ argAbs1Maximum $ raise (-x0) xs+ in if x0>=0 then (p1,p0) else (p0,p1)+++{- |+prop> \(ArbVector xs) -> Vector.negate xs == Vector.scale minusOne xs+prop> \(ArbVector xs) -> Vector.scale 2 xs == xs |+| xs+-}+scale, _scale, (.*|) ::+ (Shape.C sh, Class.Floating a) =>+ a -> Vector sh a -> Vector sh a+(.*|) = scale++scale alpha (Array sh x) = Array.unsafeCreateWithSize sh $ \n syPtr -> do+ evalContT $ do+ alphaPtr <- Call.number alpha+ nPtr <- Call.cint n+ sxPtr <- ContT $ withForeignPtr x+ incxPtr <- Call.cint 1+ incyPtr <- Call.cint 1+ liftIO $ Blas.copy nPtr sxPtr incxPtr syPtr incyPtr+ liftIO $ Blas.scal nPtr alphaPtr syPtr incyPtr++_scale a (Array sh b) = Array.unsafeCreateWithSize sh $ \n cPtr -> do+ let m = 1+ let k = 1+ evalContT $ do+ transaPtr <- Call.char 'N'+ transbPtr <- Call.char 'N'+ mPtr <- Call.cint m+ kPtr <- Call.cint k+ nPtr <- Call.cint n+ alphaPtr <- Call.number Scalar.one+ aPtr <- Call.number a+ ldaPtr <- Call.leadingDim m+ bPtr <- ContT $ withForeignPtr b+ ldbPtr <- Call.leadingDim k+ betaPtr <- Call.number Scalar.zero+ ldcPtr <- Call.leadingDim m+ liftIO $+ Blas.gemm+ transaPtr transbPtr mPtr nPtr kPtr alphaPtr+ aPtr ldaPtr bPtr ldbPtr betaPtr cPtr ldcPtr+++scaleReal ::+ (Shape.C sh, Class.Floating a) =>+ RealOf a -> Vector sh a -> Vector sh a+scaleReal =+ getScaleReal $+ Class.switchFloating+ (ScaleReal scale)+ (ScaleReal scale)+ (ScaleReal $ \x -> recomplex . scale x . decomplex)+ (ScaleReal $ \x -> recomplex . scale x . decomplex)++newtype ScaleReal f a = ScaleReal {getScaleReal :: RealOf a -> f a -> f a}+++decomplex ::+ (Class.Real a) =>+ Vector sh (Complex a) -> Vector (sh, ComplexShape) a+decomplex (Array sh a) = Array (sh, Shape.static) (castForeignPtr a)++recomplex ::+ (Class.Real a) =>+ Vector (sh, ComplexShape) a -> Vector sh (Complex a)+recomplex (Array (sh, Shape.NestedTuple _) a) = Array sh (castForeignPtr a)++++infixl 6 |+|, |-|, `add`, `sub`+++{- |+prop> \(ArbVector2 xs ys) -> xs |+| ys == ys |+| xs+prop> \(ArbVector2 xs ys) -> xs == xs |-| ys |+| ys+-}+add, sub, (|+|), (|-|) ::+ (Shape.C sh, Eq sh, Class.Floating a) =>+ Vector sh a -> Vector sh a -> Vector sh a+add = mac Scalar.one+sub x y = mac minusOne y x++(|+|) = add+(|-|) = sub++mac ::+ (Shape.C sh, Eq sh, Class.Floating a) =>+ a -> Vector sh a -> Vector sh a -> Vector sh a+mac alpha (Array shX x) (Array shY y) =+ Array.unsafeCreateWithSize shX $ \n szPtr -> do+ Call.assert "mac: shapes mismatch" (shX == shY)+ evalContT $ do+ nPtr <- Call.cint n+ saPtr <- Call.number alpha+ sxPtr <- ContT $ withForeignPtr x+ incxPtr <- Call.cint 1+ syPtr <- ContT $ withForeignPtr y+ incyPtr <- Call.cint 1+ inczPtr <- Call.cint 1+ liftIO $ Blas.copy nPtr syPtr incyPtr szPtr inczPtr+ liftIO $ Blas.axpy nPtr saPtr sxPtr incxPtr szPtr inczPtr+++{- |+prop> \(ArbVector xs) -> xs == Vector.negate (Vector.negate xs)+-}+negate :: (Shape.C sh, Class.Floating a) => Vector sh a -> Vector sh a+negate =+ getConjugate $+ Class.switchFloating+ (Conjugate $ scaleReal Scalar.minusOne)+ (Conjugate $ scaleReal Scalar.minusOne)+ (Conjugate $ scaleReal Scalar.minusOne)+ (Conjugate $ scaleReal Scalar.minusOne)+++{- |+prop> QC.forAll (genNumber maxElem) $ \d (ArbVector xs) -> xs == Vector.raise (-d) (Vector.raise d xs)+-}+raise :: (Shape.C sh, Class.Floating a) => a -> Vector sh a -> Vector sh a+raise c (Array shX x) =+ Array.unsafeCreateWithSize shX $ \n yPtr -> evalContT $ do+ nPtr <- Call.cint n+ cPtr <- Call.number c+ onePtr <- Call.number Scalar.one+ inccPtr <- Call.cint 0+ xPtr <- ContT $ withForeignPtr x+ inc1Ptr <- Call.cint 1+ liftIO $ do+ Blas.copy nPtr xPtr inc1Ptr yPtr inc1Ptr+ Blas.axpy nPtr onePtr cPtr inccPtr yPtr inc1Ptr+++{- |+prop> \(ArbVector2 xs ys) -> Vector.mul xs ys == Vector.mul ys xs+-}+mul ::+ (Shape.C sh, Eq sh, Class.Floating a) =>+ Vector sh a -> Vector sh a -> Vector sh a+mul = mulConjugation NonConjugated++{- |+prop> \(ArbVector2 xs ys) -> Vector.mulConj xs ys == Vector.mul (Vector.conjugate xs) ys+-}+mulConj ::+ (Shape.C sh, Eq sh, Class.Floating a) =>+ Vector sh a -> Vector sh a -> Vector sh a+mulConj = mulConjugation Conjugated++mulConjugation ::+ (Shape.C sh, Eq sh, Class.Floating a) =>+ Conjugation -> Vector sh a -> Vector sh a -> Vector sh a+mulConjugation conj (Array shA a) (Array shX x) =+ Array.unsafeCreateWithSize shX $ \n yPtr -> do+ Call.assert "mul: shapes mismatch" (shA == shX)+ evalContT $ do+ aPtr <- ContT $ withForeignPtr a+ xPtr <- ContT $ withForeignPtr x+ liftIO $ Private.mul conj n aPtr 1 xPtr 1 yPtr 1+++newtype Conjugate f a = Conjugate {getConjugate :: f a -> f a}++conjugate ::+ (Shape.C sh, Class.Floating a) =>+ Vector sh a -> Vector sh a+conjugate =+ getConjugate $+ Class.switchFloating+ (Conjugate id)+ (Conjugate id)+ (Conjugate complexConjugate)+ (Conjugate complexConjugate)++complexConjugate ::+ (Shape.C sh, Class.Real a) =>+ Vector sh (Complex a) -> Vector sh (Complex a)+complexConjugate (Array sh x) = Array.unsafeCreateWithSize sh $ \n syPtr ->+ evalContT $ do+ nPtr <- Call.cint n+ sxPtr <- ContT $ withForeignPtr x+ incxPtr <- Call.cint 1+ incyPtr <- Call.cint 1+ liftIO $ copyConjugate nPtr sxPtr incxPtr syPtr incyPtr+++fromReal ::+ (Shape.C sh, Class.Floating a) => Vector sh (RealOf a) -> Vector sh a+fromReal =+ getFromReal $+ Class.switchFloating+ (FromReal id)+ (FromReal id)+ (FromReal complexFromReal)+ (FromReal complexFromReal)++newtype FromReal f a = FromReal {getFromReal :: f (RealOf a) -> f a}++toComplex ::+ (Shape.C sh, Class.Floating a) => Vector sh a -> Vector sh (ComplexOf a)+toComplex =+ getToComplex $+ Class.switchFloating+ (ToComplex complexFromReal)+ (ToComplex complexFromReal)+ (ToComplex id)+ (ToComplex id)++newtype ToComplex f a = ToComplex {getToComplex :: f a -> f (ComplexOf a)}++complexFromReal ::+ (Shape.C sh, Class.Real a) => Vector sh a -> Vector sh (Complex a)+complexFromReal (Array sh x) =+ Array.unsafeCreateWithSize sh $ \n yPtr ->+ case realPtr yPtr of+ yrPtr -> evalContT $ do+ nPtr <- Call.cint n+ xPtr <- ContT $ withForeignPtr x+ incxPtr <- Call.cint 1+ incyPtr <- Call.cint 2+ inczPtr <- Call.cint 0+ zPtr <- Call.number Scalar.zero+ liftIO $ do+ Blas.copy nPtr xPtr incxPtr yrPtr incyPtr+ Blas.copy nPtr zPtr inczPtr (advancePtr yrPtr 1) incyPtr+++realPart ::+ (Shape.C sh, Class.Floating a) => Vector sh a -> Vector sh (RealOf a)+realPart =+ getToReal $+ Class.switchFloating+ (ToReal id)+ (ToReal id)+ (ToReal $ RowMajor.takeColumn ixReal . decomplex)+ (ToReal $ RowMajor.takeColumn ixReal . decomplex)+{-+ (ToReal $ RowMajor.takeColumn Complex.realPart . decomplex)+ (ToReal $ RowMajor.takeColumn Complex.realPart . decomplex)+-}++newtype ToReal f a = ToReal {getToReal :: f a -> f (RealOf a)}++imaginaryPart ::+ (Shape.C sh, Class.Real a) => Vector sh (Complex a) -> Vector sh a+imaginaryPart = RowMajor.takeColumn ixImaginary . decomplex+-- imaginaryPart = RowMajor.takeColumn Complex.imagPart . decomplex+++zipComplex ::+ (Shape.C sh, Eq sh, Class.Real a) =>+ Vector sh a -> Vector sh a -> Vector sh (Complex a)+zipComplex (Array shr xr) (Array shi xi) =+ Array.unsafeCreateWithSize shr $ \n yPtr -> evalContT $ do+ liftIO $ Call.assert "zipComplex: shapes mismatch" (shr==shi)+ nPtr <- Call.cint n+ xrPtr <- ContT $ withForeignPtr xr+ xiPtr <- ContT $ withForeignPtr xi+ let yrPtr = realPtr yPtr+ incxPtr <- Call.cint 1+ incyPtr <- Call.cint 2+ liftIO $ do+ Blas.copy nPtr xrPtr incxPtr yrPtr incyPtr+ Blas.copy nPtr xiPtr incxPtr (advancePtr yrPtr 1) incyPtr++unzipComplex ::+ (Shape.C sh, Class.Real a) =>+ Vector sh (Complex a) -> (Vector sh a, Vector sh a)+unzipComplex x = (realPart x, imaginaryPart x)
+ src/Numeric/BLAS/Vector/Slice.hs view
@@ -0,0 +1,903 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}+module Numeric.BLAS.Vector.Slice (+ T,+ shape,+ Vector,+ RealOf,+ ComplexOf,+ slice,+ fromVector,+ toVector,+ extract,+ access,+ dot, inner,+ sum,+ absSum,+ norm2,+ norm2Squared,+ normInf,+ normInf1,+ argAbsMaximum,+ argAbs1Maximum,+ product,+ scale, scaleReal,+ add, sub,+ negate, raise,+ mac,+ mul, mulConj,+ minimum, argMinimum,+ maximum, argMaximum,+ limits, argLimits,++ conjugate,+ fromReal,+ toComplex,+ realFromComplexVector,+ realPart,+ imaginaryPart,+ zipComplex,+ unzipComplex,+ ) where++import qualified Numeric.BLAS.Slice as Slice++import qualified Numeric.BLAS.Scalar as Scalar+import qualified Numeric.BLAS.Private as Private+import Numeric.BLAS.Matrix.Modifier (Conjugation(NonConjugated, Conjugated))+import Numeric.BLAS.Scalar (ComplexOf, RealOf)+import Numeric.BLAS.Private (ComplexShape, copyConjugate, realPtr)++import qualified Numeric.BLAS.FFI.Generic as Blas+import qualified Numeric.BLAS.FFI.Complex as BlasComplex+import qualified Numeric.BLAS.FFI.Real as BlasReal+import qualified Numeric.Netlib.Utility as Call+import qualified Numeric.Netlib.Class as Class++import qualified Foreign.Marshal.Array.Guarded as ForeignArray+import Foreign.Marshal.Array (advancePtr)+import Foreign.ForeignPtr (withForeignPtr, castForeignPtr)+import Foreign.Ptr (Ptr)+import Foreign.Storable (Storable, peek, peekElemOff)+import Foreign.C.Types (CInt)++import System.IO.Unsafe (unsafePerformIO)++import Control.Monad.Trans.Cont (ContT(ContT), evalContT)+import Control.Monad.IO.Class (liftIO)+import Control.Monad (fmap, return, (=<<))+import Control.Applicative (liftA2, pure, (<$>), (<*>))++import qualified Data.Array.Comfort.Storable.Unchecked as Array+import qualified Data.Array.Comfort.Shape as Shape+import Data.Array.Comfort.Storable.Unchecked (Array(Array))++import Text.Show (Show)+import Data.Function (($), (.))+import Data.Complex (Complex)+import Data.Maybe (Maybe(Nothing,Just), maybe)+import Data.Tuple.HT (mapFst, uncurry3)+import Data.Tuple (fst, snd, uncurry)+import Data.Ord ((>=))+import Data.Eq (Eq, (==))++import Prelude (Int, fromIntegral, (-), (+), (*), error, IO)+++{- $setup+>>> :set -XTypeOperators+>>> import qualified Test.Slice as TestSlice+>>> import Test.NumberModule.Type (Number_)+>>> import Test.Generator (genNumber)+>>> import Test.Utility (approx, approxReal)+>>>+>>> import qualified Numeric.BLAS.Matrix.RowMajor as Matrix+>>> import qualified Numeric.BLAS.Vector.Slice as VectorSlice+>>> import qualified Numeric.BLAS.Vector as Vector+>>> import qualified Numeric.BLAS.Slice as Slice+>>> import qualified Numeric.BLAS.Scalar as Scalar+>>> import qualified Numeric.Netlib.Class as Class+>>> import qualified Data.Array.Comfort.Storable as Array+>>> import qualified Data.Array.Comfort.Shape as Shape+>>> import qualified Data.List as List+>>> import Numeric.BLAS.Vector ((+++), (|+|))+>>> import Numeric.BLAS.Scalar (RealOf, absolute, minusOne)+>>> import Data.Array.Comfort.Shape ((::+)((::+)))+>>> import Data.Tuple.HT (mapPair)+>>> import Data.Complex (Complex)+>>> import Control.Monad (replicateM)+>>> import Control.Applicative (liftA2)+>>>+>>> import qualified Test.QuickCheck as QC+>>> import Test.QuickCheck ((==>))+>>>+>>> type Real_ = RealOf Number_+>>> type Complex_ = Complex Real_+>>>+>>> maxElem :: Integer+>>> maxElem = 10+>>>+>>> maxDim :: Int+>>> maxDim = 100+>>>+>>> maxDim1 :: Int+>>> maxDim1 = 10+>>>+>>> type ShapeInt = Shape.ZeroBased Int+>>> type Shape = ShapeInt ::+ (ShapeInt, ShapeInt) ::+ ShapeInt+>>> type Vector = Vector.Vector Shape+>>> type Sliced = VectorSlice.T Shape ShapeInt+>>>+>>> genDim :: QC.Gen Int+>>> genDim = QC.choose (0,maxDim)+>>>+>>> genShapeDim :: Int -> QC.Gen Shape+>>> genShapeDim numRows = do+>>> left <- fmap Shape.ZeroBased $ QC.choose (0,maxDim)+>>> right <- fmap Shape.ZeroBased $ QC.choose (0,maxDim)+>>> columns <- fmap Shape.ZeroBased $ QC.choose (1,maxDim1)+>>> return (left ::+ (Shape.ZeroBased numRows, columns) ::+ right)+>>>+>>> genShape :: QC.Gen Shape+>>> genShape = genShapeDim =<< QC.choose (0,maxDim1)+>>>+>>> genVector ::+>>> (Shape.C sh, Class.Floating a) =>+>>> sh -> QC.Gen a -> QC.Gen (Vector.Vector sh a)+>>> genVector shape genElem =+>>> fmap (Vector.fromList shape) $+>>> replicateM (Shape.size shape) genElem+>>>+>>> forAll_ :: (Show a) => QC.Gen a -> (a -> QC.Property) -> QC.Property+>>> forAll_ = QC.forAll+>>>+>>> real_ :: QC.Gen Real_+>>> real_ = genNumber maxElem+>>>+>>> complex_ :: QC.Gen Complex_+>>> complex_ = genNumber maxElem+>>>+>>> number_ :: QC.Gen Number_+>>> number_ = genNumber maxElem+>>>+>>> isNonEmpty :: Shape.C sh => VectorSlice.T shA sh a -> Bool+>>> isNonEmpty xs = Shape.size (VectorSlice.shape xs) > 0+>>>+>>> takeColumn ::+>>> (Shape.Indexed sh1, Shape.C sh, Shape.C sh0, Shape.C sh2) =>+>>> Shape.Index sh1 ->+>>> Vector.Vector (sh0 ::+ (sh, sh1) ::+ sh2) a ->+>>> VectorSlice.T (sh0 ::+ (sh, sh1) ::+ sh2) sh a+>>> takeColumn c =+>>> VectorSlice.slice+>>> (Slice.column c . Slice.left . Slice.right) . VectorSlice.fromVector+>>>+>>> listFromSlice ::+>>> (Shape.C sh, Class.Floating a) => VectorSlice.T shA sh a -> [a]+>>> listFromSlice = Vector.toList . VectorSlice.toVector+>>>+>>> genSlicedDim ::+>>> (Class.Floating a) =>+>>> Int -> QC.Gen a -> QC.Gen (Int, Vector a)+>>> genSlicedDim numRows genElem = do+>>> shape@(_::+(_rows,columns)::+_) <- genShapeDim numRows+>>> c <- QC.elements (Shape.indices columns)+>>> fmap ((,) c) $ genVector shape genElem+>>>+>>> genSliced ::+>>> (Class.Floating a) =>+>>> QC.Gen a -> QC.Gen (Int, Vector a)+>>> genSliced genElem = flip genSlicedDim genElem =<< genDim+>>>+>>> shrinkSliced ::+>>> (Shape.C sh0, Shape.Indexed sh1, QC.Arbitrary a, Class.Floating a) =>+>>> (Shape.Index sh1,+>>> Vector.Vector (ShapeInt ::+ (sh0, sh1) ::+ ShapeInt) a) ->+>>> [(Shape.Index sh1,+>>> Vector.Vector (ShapeInt ::+ (sh0, sh1) ::+ ShapeInt) a)]+>>> shrinkSliced (c,xs) =+>>> let xs0 = Vector.takeLeft xs in+>>> let xs1 = Vector.takeRight xs in+>>> let xs10 = Vector.takeLeft xs1 in+>>> let xs11 = Vector.takeRight xs1 in+>>> map (\(ysl,ysr) ->+>>> (c,+>>> Vector.autoFromList ysl +++ xs10 +++ Vector.autoFromList ysr)) $+>>> QC.shrink (Vector.toList xs0, Vector.toList xs11)+>>>+>>> forSliced ::+>>> (QC.Testable prop, QC.Arbitrary a, Class.Floating a, Show a) =>+>>> QC.Gen a -> (Sliced a -> prop) -> QC.Property+>>> forSliced genElem prop =+>>> QC.forAllShrink (genSliced genElem) shrinkSliced+>>> (prop . uncurry takeColumn)+>>>+>>> genSliced2 ::+>>> (Class.Floating a) =>+>>> QC.Gen a -> QC.Gen ((Int, Vector a), (Int, Vector a))+>>> genSliced2 genElem = do+>>> dim <- genDim+>>> liftA2 (,) (genSlicedDim dim genElem) (genSlicedDim dim genElem)+>>>+>>> forSliced2 ::+>>> (QC.Testable prop, Class.Floating a, Show a) =>+>>> QC.Gen a -> (Sliced a -> Sliced a -> prop) -> QC.Property+>>> forSliced2 genElem prop =+>>> QC.forAll (genSliced2 genElem)+>>> (uncurry prop . mapPair (uncurry takeColumn, uncurry takeColumn))+-}+++type Vector = Array+++shape :: T sh slice a -> slice+shape (Cons (Slice.Cons _s _k slc) _arr) = slc++mapShape :: (slice0 -> slice1) -> T sh slice0 a -> T sh slice1 a+mapShape f (Cons (Slice.Cons s k slc) arr) =+ Cons (Slice.Cons s k (f slc)) arr+++increment :: T sh slice a -> Int+increment (Cons (Slice.Cons _s k _slc) _arr) = k+++startArg ::+ (Storable a) =>+ T sh slice a -> Call.FortranIO r (Ptr a)+startArg (Cons (Slice.Cons s _k _slice) (Array _sh x)) = do+ sxPtr <- ContT $ withForeignPtr x+ return (advancePtr sxPtr s)++sliceArg ::+ (Storable a) =>+ T sh slice a -> Call.FortranIO r (Ptr a, Ptr CInt)+sliceArg x =+ liftA2 (,) (startArg x) (Call.cint $ increment x)++sizeSliceArg ::+ (Shape.C sh, Storable a) =>+ T shA sh a -> ContT r IO (Ptr CInt, Ptr a, Ptr CInt)+sizeSliceArg x =+ liftA2+ (\nPtr (xPtr,incxPtr) -> (nPtr, xPtr,incxPtr))+ (Call.cint $ Shape.size $ shape x)+ (sliceArg x)++infixl 4 <*|>++(<*|>) ::+ (Storable a) =>+ Call.FortranIO r (Ptr a -> Ptr CInt -> b) ->+ T sh slice a ->+ Call.FortranIO r b+f <*|> x = fmap uncurry f <*> sliceArg x+++data T sh slice a = Cons (Slice.T slice) (Vector sh a)+ deriving (Show)++toVector ::+ (Shape.C slice, Class.Floating a) =>+ T sh slice a -> Vector slice a+toVector x =+ Array.unsafeCreateWithSize (shape x) $ \n syPtr ->+ evalContT $ Call.run $+ pure Blas.copy+ <*> Call.cint n+ <*|> x+ <*> pure syPtr+ <*> Call.cint 1++fromVector :: (Shape.C sh) => Vector sh a -> T sh sh a+fromVector xs = Cons (Slice.fromShape $ Array.shape xs) xs++slice :: (Slice.T shA -> Slice.T shB) -> T sh shA a -> T sh shB a+slice f (Cons slc xs) = Cons (f slc) xs++{- |+prop> QC.forAll genShape $ \shape@(_::+(_rows,columns)::+_) -> QC.forAll (QC.elements (Shape.indices columns)) $ \c -> QC.forAll (genVector shape number_) $ \xs -> VectorSlice.extract (Slice.column c . Slice.left . Slice.right) xs == Matrix.takeColumn c (Vector.takeLeft (Vector.takeRight (xs :: Vector Number_)))++prop> forAll_ (TestSlice.genShapeSelect 4 100) $ \(TestSlice.ShapeSelect sh select) -> QC.forAll (genVector sh number_) $ \xs -> case TestSlice.instantiate sh select of TestSlice.Extraction slice cut -> VectorSlice.extract slice xs == cut xs+-}+extract ::+ (Shape.C slice, Shape.C sh, Class.Floating a) =>+ (Slice.T sh -> Slice.T slice) -> Vector sh a -> Vector slice a+extract slc xs = toVector $ slice slc $ fromVector xs+++access, (!) ::+ (Shape.C shA, Shape.Indexed sh, Storable a) =>+ T shA sh a -> Shape.Index sh -> a+access (Cons (Slice.Cons s k ssh) (Array sh x)) ix =+ Array (Shape.Deferred sh) x+ Array.!+ Shape.DeferredIndex (s + k * Shape.offset ssh ix)+(!) = access++++dot ::+ (Shape.C sh, Eq sh, Class.Floating a) =>+ T shA sh a -> T shB sh a -> a+dot =+ runDot $+ Class.switchFloating+ (Dot dotReal)+ (Dot dotReal)+ (Dot dotComplex)+ (Dot dotComplex)+++{- |+prop> forSliced2 number_ $ \xs ys -> VectorSlice.inner xs ys == Vector.dot (VectorSlice.conjugate xs) (VectorSlice.toVector ys)+prop> forSliced number_ $ \xs -> VectorSlice.inner xs xs == Scalar.fromReal (VectorSlice.norm2Squared xs)+-}+inner ::+ (Shape.C sh, Eq sh, Class.Floating a) =>+ T shA sh a -> T shB sh a -> a+inner =+ runDot $+ Class.switchFloating+ (Dot dotReal)+ (Dot dotReal)+ (Dot $ innerComplex . toVector)+ (Dot $ innerComplex . toVector)+++newtype Dot f g a = Dot {runDot :: f a -> g a -> a}++dotReal ::+ (Shape.C sh, Eq sh, Class.Real a) =>+ T shA sh a -> T shB sh a -> a+dotReal x y = unsafePerformIO $ do+ let shX = shape x+ let shY = shape y+ Call.assert "dot: shapes mismatch" (shX == shY)+ evalContT $ do+ nPtr <- Call.cint $ Shape.size shX+ (sxPtr, incxPtr) <- sliceArg x+ (syPtr, incyPtr) <- sliceArg y+ liftIO $ BlasReal.dot nPtr sxPtr incxPtr syPtr incyPtr++{-+We cannot use 'cdot' because Haskell's FFI+does not support Complex numbers as return values.+-}+dotComplex ::+ (Shape.C sh, Eq sh, Class.Real a) =>+ T shA sh (Complex a) -> T shB sh (Complex a) -> Complex a+dotComplex x y = unsafePerformIO $ do+ Call.assert "dot: shapes mismatch" (shape x == shape y)+ evalContT $ do+ transPtr <- Call.char 'N'+ mPtr <- Call.cint 1+ nPtr <- Call.cint $ Shape.size $ shape x+ alphaPtr <- Call.number Scalar.one+ (xPtr, ldxPtr) <- sliceArg x+ (yPtr, incyPtr) <- sliceArg y+ betaPtr <- Call.number Scalar.zero+ zPtr <- Call.alloca+ inczPtr <- Call.cint 1+ liftIO $+ Private.gemv+ transPtr mPtr nPtr alphaPtr xPtr ldxPtr+ yPtr incyPtr betaPtr zPtr inczPtr+ liftIO $ peek zPtr++{-+In contrast to Vector.inner+we cannot use gemv with 'C' because we need leading dimension > 1.+-}+innerComplex ::+ (Shape.C sh, Eq sh, Class.Real a) =>+ Vector sh (Complex a) -> T shB sh (Complex a) -> Complex a+innerComplex (Array shX x) y = unsafePerformIO $ do+ Call.assert "dot: shapes mismatch" (shX == shape y)+ evalContT $ do+ let m = Shape.size shX+ transPtr <- Call.char 'C'+ mPtr <- Call.cint m+ nPtr <- Call.cint 1+ alphaPtr <- Call.number Scalar.one+ xPtr <- ContT $ withForeignPtr x+ ldxPtr <- Call.leadingDim m+ (yPtr, incyPtr) <- sliceArg y+ betaPtr <- Call.number Scalar.zero+ zPtr <- Call.alloca+ inczPtr <- Call.cint 1+ liftIO $+ Private.gemv+ transPtr mPtr nPtr alphaPtr xPtr ldxPtr+ yPtr incyPtr betaPtr zPtr inczPtr+ liftIO $ peek zPtr+++{- |+prop> forSliced number_ $ \xs -> VectorSlice.sum xs == List.sum (listFromSlice xs)+-}+sum :: (Shape.C sh, Class.Floating a) => T shA sh a -> a+sum x = unsafePerformIO $ evalContT $ do+ xPtr <- startArg x+ liftIO $ Private.sum (Shape.size $ shape x) xPtr (increment x)+++{- |+Sum of the absolute values of real numbers or components of complex numbers.+For real numbers it is equivalent to 'Numeric.LAPACK.Vector.norm1'.+-}+absSum :: (Shape.C sh, Class.Floating a) => T shA sh a -> RealOf a+absSum arr = unsafePerformIO $+ evalContT $ liftIO . uncurry3 asum =<< sizeSliceArg arr++asum :: Class.Floating a => Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a)+asum =+ getNrm $+ Class.switchFloating+ (Nrm BlasReal.asum) (Nrm BlasReal.asum)+ (Nrm BlasComplex.casum) (Nrm BlasComplex.casum)+++{- |+Euclidean norm of a vector or Frobenius norm of a matrix.+-}+norm2 :: (Shape.C sh, Class.Floating a) => T shA sh a -> RealOf a+norm2 arr = unsafePerformIO $+ evalContT $ liftIO . uncurry3 nrm2 =<< sizeSliceArg arr++nrm2 :: Class.Floating a => Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a)+nrm2 =+ getNrm $+ Class.switchFloating+ (Nrm BlasReal.nrm2) (Nrm BlasReal.nrm2)+ (Nrm BlasComplex.cnrm2) (Nrm BlasComplex.cnrm2)++newtype Nrm a = Nrm {getNrm :: Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a)}+++newtype Norm f a = Norm {getNorm :: f a -> RealOf a}++norm2Squared :: (Shape.C sh, Class.Floating a) => T shA sh a -> RealOf a+norm2Squared =+ getNorm $+ Class.switchFloating+ (Norm norm2SquaredReal)+ (Norm norm2SquaredReal)+ (Norm norm2SquaredComplex)+ (Norm norm2SquaredComplex)++norm2SquaredReal :: (Shape.C sh, Class.Real a) => T shA sh a -> a+norm2SquaredReal x =+ unsafePerformIO $ evalContT $ do+ (nPtr, sxPtr, incxPtr) <- sizeSliceArg x+ liftIO $ BlasReal.dot nPtr sxPtr incxPtr sxPtr incxPtr++norm2SquaredComplex :: (Shape.C sh, Class.Real a) => T shA sh (Complex a) -> a+norm2SquaredComplex x =+ unsafePerformIO $ evalContT $ do+ nPtr <- Call.cint $ Shape.size $ shape x+ xPtr <- startArg x+ let xrPtr = realPtr xPtr+ let xiPtr = advancePtr xrPtr 1+ incxPtr <- Call.cint (increment x * 2)+ liftIO $+ liftA2 (+)+ (BlasReal.dot nPtr xrPtr incxPtr xrPtr incxPtr)+ (BlasReal.dot nPtr xiPtr incxPtr xiPtr incxPtr)+++{- |+prop> forSliced number_ $ \xs -> VectorSlice.normInf xs == List.maximum (0 : List.map absolute (listFromSlice xs))+-}+normInf :: (Shape.C sh, Class.Floating a) => T shA sh a -> RealOf a+normInf arr = unsafePerformIO $+ fmap (Scalar.absolute . maybe Scalar.zero snd) $ absMax arr++{- |+Computes (almost) the infinity norm of the vector.+For complex numbers every element is replaced+by the sum of the absolute component values first.++prop> forSliced number_ $ \xs -> VectorSlice.normInf1 xs == List.maximum (0 : List.map Scalar.norm1 (listFromSlice xs))+-}+normInf1 :: (Shape.C sh, Class.Floating a) => T shA sh a -> RealOf a+normInf1 x = unsafePerformIO $+ evalContT $ do+ (nPtr, sxPtr, incxPtr) <- sizeSliceArg x+ liftIO $+ fmap (Scalar.norm1 . maybe Scalar.zero snd) $+ peekElemOff1 sxPtr (increment x) =<< Blas.iamax nPtr sxPtr incxPtr+++{- |+Returns the index and value of the element with the maximal absolute value.+Caution: It actually returns the value of the element, not its absolute value!++prop> forSliced number_ $ \xs -> isNonEmpty xs ==> let (xi,xm) = VectorSlice.argAbsMaximum xs in VectorSlice.access xs xi == xm+prop> forSliced number_ $ \xs -> isNonEmpty xs ==> let (_xi,xm) = VectorSlice.argAbsMaximum xs in List.all (\x -> absolute x <= absolute xm) $ listFromSlice xs+prop> forSliced number_ $ \xs -> forSliced number_ $ \ys -> isNonEmpty xs && isNonEmpty ys ==> let (_xi,xm) = VectorSlice.argAbsMaximum xs; (_yi,ym) = VectorSlice.argAbsMaximum ys; (zi,zm) = Vector.argAbsMaximum (VectorSlice.toVector xs +++ VectorSlice.toVector ys) in case zi of Left _ -> xm==zm && absolute xm >= absolute ym; Right _ -> ym==zm && absolute xm < absolute ym+-}+argAbsMaximum ::+ (Shape.InvIndexed sh, Class.Floating a) =>+ T shA sh a -> (Shape.Index sh, a)+argAbsMaximum arr = unsafePerformIO $+ fmap+ (maybe+ (error "Vector.argAbsMaximum: empty vector")+ (mapFst (Shape.uncheckedIndexFromOffset $ shape arr))) $+ absMax arr++absMax ::+ (Shape.C sh, Class.Floating a) =>+ T shA sh a -> IO (Maybe (Int, a))+absMax x =+ case Scalar.complexSingletonOfFunctor x of+ Scalar.Real -> evalContT $ do+ (nPtr, sxPtr, incxPtr) <- sizeSliceArg x+ liftIO $+ peekElemOff1 sxPtr (increment x) =<< Blas.iamax nPtr sxPtr incxPtr+ Scalar.Complex -> evalContT $ do+ let n = Shape.size $ shape x+ sxPtr <- startArg x+ let incx = increment x+ liftIO $ peekElemOff1 sxPtr incx =<< absMaxComplex n sxPtr incx++absMaxComplex :: (Class.Real a) => Int -> Ptr (Complex a) -> Int -> IO CInt+absMaxComplex n sxPtr incx =+ ForeignArray.alloca n $ \syPtr -> do+ let xrPtr = realPtr sxPtr+ let incx2 = 2*incx+ Private.mul NonConjugated n xrPtr incx2 xrPtr incx2 syPtr 1+ let xiPtr = advancePtr xrPtr 1+ Private.mulAdd NonConjugated n xiPtr incx2 xiPtr incx2 Scalar.one syPtr 1+ evalContT $ do+ nPtr <- Call.cint n+ incyPtr <- Call.cint 1+ liftIO $ Blas.iamax nPtr syPtr incyPtr+++{- |+Returns the index and value of the element with the maximal absolute value.+The function does not strictly compare the absolute value of a complex number+but the sum of the absolute complex components.+Caution: It actually returns the value of the element, not its absolute value!++prop> forSliced real_ $ \xs -> isNonEmpty xs ==> VectorSlice.argAbsMaximum xs == VectorSlice.argAbs1Maximum xs+-}+argAbs1Maximum ::+ (Shape.InvIndexed sh, Class.Floating a) =>+ T shA sh a -> (Shape.Index sh, a)+argAbs1Maximum x = unsafePerformIO $+ evalContT $ do+ (nPtr, sxPtr, incxPtr) <- sizeSliceArg x+ liftIO $+ fmap+ (maybe+ (error "Vector.argAbs1Maximum: empty vector")+ (mapFst (Shape.uncheckedIndexFromOffset $ shape x))) $+ peekElemOff1 sxPtr (increment x) =<< Blas.iamax nPtr sxPtr incxPtr++peekElemOff1 :: (Storable a) => Ptr a -> Int -> CInt -> IO (Maybe (Int, a))+peekElemOff1 ptr inc k1 =+ let k1i = fromIntegral k1+ ki = k1i-1+ in if k1i == 0+ then return Nothing+ else Just . (,) ki <$> peekElemOff ptr (ki*inc)+++{- |+prop> QC.forAll genShape $ \sh@(_::+(_rows,columns)::+_) -> QC.forAll (QC.elements (Shape.indices columns)) $ \c -> QC.forAll (genVector sh $ genNumber 3) $ \xt -> let xs = takeColumn c xt in approx 1e-2 (VectorSlice.product xs) (List.product (listFromSlice (xs :: Sliced Number_)))+-}+product :: (Shape.C sh, Class.Floating a) => T shA sh a -> a+product x = unsafePerformIO $ evalContT $ do+ xPtr <- startArg x+ liftIO $ Private.product (Shape.size $ shape x) xPtr (increment x)+++{- |+For restrictions see 'limits'.++prop> forSliced real_ $ \xs -> isNonEmpty xs ==> VectorSlice.minimum xs == List.minimum (listFromSlice xs)+prop> forSliced real_ $ \xs -> isNonEmpty xs ==> VectorSlice.maximum xs == List.maximum (listFromSlice xs)+-}+minimum, maximum :: (Shape.C shA, Shape.C sh, Class.Real a) => T shA sh a -> a+minimum = fst . limits+maximum = snd . limits++{- |+For restrictions see 'limits'.+-}+argMinimum, argMaximum ::+ (Shape.C shA, Shape.InvIndexed sh, Shape.Index sh ~ ix, Class.Real a) =>+ T shA sh a -> (ix,a)+argMinimum = fst . argLimits+argMaximum = snd . argLimits++{- |+prop> forSliced real_ $ \xs -> isNonEmpty xs ==> VectorSlice.limits xs == Array.limits (VectorSlice.toVector xs)++In contrast to 'Array.limits'+this implementation is based on fast BLAS functions.+It should be faster than @Array.minimum@ and @Array.maximum@+although it is certainly not as fast as possible.+Boths limits share the precision of the limit with the larger absolute value.+This implies for example that you cannot rely on the property+that @raise (- minimum x) x@ has only non-negative elements.+-}+limits :: (Shape.C shA, Shape.C sh, Class.Real a) => T shA sh a -> (a,a)+limits xs0 =+ let xs = mapShape Shape.Deferred xs0+ x0 = snd $ argAbs1Maximum xs+ x1 = xs ! fst (argAbs1Maximum (fromVector (raise (-x0) xs)))+ in if x0>=0 then (x1,x0) else (x0,x1)++{- |+For restrictions see 'limits'.+-}+argLimits ::+ (Shape.C shA, Shape.InvIndexed sh, Shape.Index sh ~ ix, Class.Real a) =>+ T shA sh a -> ((ix,a),(ix,a))+argLimits xs =+ let p0@(_i0,x0) = argAbs1Maximum xs+ p1 = (i1,xs!i1); i1 = fst $ argAbs1Maximum $ fromVector $ raise (-x0) xs+ in if x0>=0 then (p1,p0) else (p0,p1)+++{- |+prop> forSliced number_ $ \xs -> VectorSlice.negate xs == VectorSlice.scale minusOne xs+prop> forSliced number_ $ \xs -> VectorSlice.scale 2 xs == VectorSlice.add xs xs+-}+scale, _scale ::+ (Shape.C sh, Class.Floating a) =>+ a -> T shA sh a -> Vector sh a++scale alpha x = Array.unsafeCreateWithSize (shape x) $ \n syPtr -> do+ evalContT $ do+ alphaPtr <- Call.number alpha+ nPtr <- Call.cint n+ (sxPtr, incxPtr) <- sliceArg x+ incyPtr <- Call.cint 1+ liftIO $ Blas.copy nPtr sxPtr incxPtr syPtr incyPtr+ liftIO $ Blas.scal nPtr alphaPtr syPtr incyPtr++_scale a b = Array.unsafeCreateWithSize (shape b) $ \n cPtr -> do+ let m = 1+ let k = 1+ evalContT $ do+ transaPtr <- Call.char 'N'+ transbPtr <- Call.char 'N'+ mPtr <- Call.cint m+ kPtr <- Call.cint k+ nPtr <- Call.cint n+ alphaPtr <- Call.number Scalar.one+ aPtr <- Call.number a+ ldaPtr <- Call.leadingDim m+ (bPtr, ldbPtr) <- sliceArg b+ betaPtr <- Call.number Scalar.zero+ ldcPtr <- Call.leadingDim m+ liftIO $+ Blas.gemm+ transaPtr transbPtr mPtr nPtr kPtr alphaPtr+ aPtr ldaPtr bPtr ldbPtr betaPtr cPtr ldcPtr+++{- |+Complex implementation requires double number of multiplications+compared to 'Numeric.BLAS.Vector.scaleReal'.+-}+scaleReal ::+ (Shape.C sh, Class.Floating a) =>+ RealOf a -> T shA sh a -> Vector sh a+scaleReal =+ getScaleReal $+ Class.switchFloating+ (ScaleReal scale)+ (ScaleReal scale)+ (ScaleReal $ scale . Scalar.fromReal)+ (ScaleReal $ scale . Scalar.fromReal)++newtype ScaleReal f g a = ScaleReal {getScaleReal :: RealOf a -> f a -> g a}++++infixl 6 `add`, `sub`+++{- |+prop> forSliced2 number_ $ \xs ys -> VectorSlice.add xs ys == VectorSlice.add ys xs+prop> forSliced2 number_ $ \xs ys -> VectorSlice.toVector xs == VectorSlice.sub xs ys |+| VectorSlice.toVector ys+-}+add, sub ::+ (Shape.C sh, Eq sh, Class.Floating a) =>+ T shA sh a -> T shB sh a -> Vector sh a+add = mac Scalar.one+sub x y = mac Scalar.minusOne y x++mac ::+ (Shape.C sh, Eq sh, Class.Floating a) =>+ a -> T shA sh a -> T shB sh a -> Vector sh a+mac alpha x y =+ Array.unsafeCreateWithSize (shape x) $ \n szPtr -> do+ Call.assert "mac: shapes mismatch" (shape x == shape y)+ evalContT $ do+ nPtr <- Call.cint n+ saPtr <- Call.number alpha+ (sxPtr, incxPtr) <- sliceArg x+ (syPtr, incyPtr) <- sliceArg y+ inczPtr <- Call.cint 1+ liftIO $ Blas.copy nPtr syPtr incyPtr szPtr inczPtr+ liftIO $ Blas.axpy nPtr saPtr sxPtr incxPtr szPtr inczPtr+++{- |+prop> forSliced number_ $ \xs -> VectorSlice.toVector xs == Vector.negate (VectorSlice.negate xs)+-}+negate :: (Shape.C sh, Class.Floating a) => T shA sh a -> Vector sh a+negate =+ getConjugate $+ Class.switchFloating+ (Conjugate $ scaleReal Scalar.minusOne)+ (Conjugate $ scaleReal Scalar.minusOne)+ (Conjugate $ scaleReal Scalar.minusOne)+ (Conjugate $ scaleReal Scalar.minusOne)+++{- |+prop> QC.forAll (genNumber maxElem) $ \d -> forSliced number_ $ \xs -> VectorSlice.toVector xs == Vector.raise (-d) (VectorSlice.raise d xs)+-}+raise :: (Shape.C sh, Class.Floating a) => a -> T shA sh a -> Vector sh a+raise c x =+ Array.unsafeCreateWithSize (shape x) $ \n yPtr -> evalContT $ do+ nPtr <- Call.cint n+ cPtr <- Call.number c+ onePtr <- Call.number Scalar.one+ inccPtr <- Call.cint 0+ (xPtr, incxPtr) <- sliceArg x+ inc1Ptr <- Call.cint 1+ liftIO $ do+ Blas.copy nPtr xPtr incxPtr yPtr inc1Ptr+ Blas.axpy nPtr onePtr cPtr inccPtr yPtr inc1Ptr+++{- |+prop> forSliced2 number_ $ \xs ys -> VectorSlice.mul xs ys == VectorSlice.mul ys xs+-}+mul ::+ (Shape.C sh, Eq sh, Class.Floating a) =>+ T shA sh a -> T shB sh a -> Vector sh a+mul = mulConjugation NonConjugated++{- |+prop> forSliced2 number_ $ \xs ys -> VectorSlice.mulConj xs ys == Vector.mul (VectorSlice.conjugate xs) (VectorSlice.toVector ys)+-}+mulConj ::+ (Shape.C sh, Eq sh, Class.Floating a) =>+ T shA sh a -> T shB sh a -> Vector sh a+mulConj = mulConjugation Conjugated++mulConjugation ::+ (Shape.C sh, Eq sh, Class.Floating a) =>+ Conjugation -> T shA sh a -> T shB sh a -> Vector sh a+mulConjugation conj a x =+ Array.unsafeCreateWithSize (shape x) $ \n yPtr -> do+ Call.assert "mul: shapes mismatch" (shape a == shape x)+ evalContT $ do+ aPtr <- startArg a+ xPtr <- startArg x+ liftIO $ Private.mul conj n aPtr (increment a) xPtr (increment x) yPtr 1+++newtype Conjugate f g a = Conjugate {getConjugate :: f a -> g a}++conjugate ::+ (Shape.C sh, Class.Floating a) =>+ T shA sh a -> Vector sh a+conjugate =+ getConjugate $+ Class.switchFloating+ (Conjugate toVector)+ (Conjugate toVector)+ (Conjugate complexConjugate)+ (Conjugate complexConjugate)++complexConjugate ::+ (Shape.C sh, Class.Real a) =>+ T shA sh (Complex a) -> Vector sh (Complex a)+complexConjugate x = Array.unsafeCreateWithSize (shape x) $ \n syPtr ->+ evalContT $ do+ nPtr <- Call.cint n+ (sxPtr, incxPtr) <- sliceArg x+ incyPtr <- Call.cint 1+ liftIO $ copyConjugate nPtr sxPtr incxPtr syPtr incyPtr+++fromReal ::+ (Shape.C sh, Class.Floating a) => T shA sh (RealOf a) -> Vector sh a+fromReal =+ getFromReal $+ Class.switchFloating+ (FromReal toVector)+ (FromReal toVector)+ (FromReal complexFromReal)+ (FromReal complexFromReal)++newtype FromReal f g a = FromReal {getFromReal :: f (RealOf a) -> g a}++toComplex ::+ (Shape.C sh, Class.Floating a) => T shA sh a -> Vector sh (ComplexOf a)+toComplex =+ getToComplex $+ Class.switchFloating+ (ToComplex complexFromReal)+ (ToComplex complexFromReal)+ (ToComplex toVector)+ (ToComplex toVector)++newtype ToComplex f g a = ToComplex {getToComplex :: f a -> g (ComplexOf a)}++complexFromReal ::+ (Shape.C sh, Class.Real a) => T shA sh a -> Vector sh (Complex a)+complexFromReal x =+ Array.unsafeCreateWithSize (shape x) $ \n yPtr ->+ case realPtr yPtr of+ yrPtr -> evalContT $ do+ nPtr <- Call.cint n+ (xPtr, incxPtr) <- sliceArg x+ incyPtr <- Call.cint 2+ inczPtr <- Call.cint 0+ zPtr <- Call.number Scalar.zero+ liftIO $ do+ Blas.copy nPtr xPtr incxPtr yrPtr incyPtr+ Blas.copy nPtr zPtr inczPtr (advancePtr yrPtr 1) incyPtr+++realFromComplexVector ::+ (Shape.C sh) =>+ Vector sh (Complex a) -> T (sh, ComplexShape) (sh, ComplexShape) a+realFromComplexVector (Array sh a) =+ let csh = (sh, Shape.static) in+ Cons (Slice.fromShape csh) (Array csh (castForeignPtr a))+++realPart ::+ (Shape.C sh, Class.Real a) =>+ T shA sh (Complex a) -> T (shA, ComplexShape) sh a+realPart (Cons (Slice.Cons s k slc) (Array sh a)) =+ Cons+ (Slice.Cons (2*s) (2*k) slc)+ (Array (sh, Shape.static) (castForeignPtr a))++imaginaryPart ::+ (Shape.C sh, Class.Real a) =>+ T shA sh (Complex a) -> T (shA, ComplexShape) sh a+imaginaryPart (Cons (Slice.Cons s k slc) (Array sh a)) =+ Cons+ (Slice.Cons (2*s+1) (2*k) slc)+ (Array (sh, Shape.static) (castForeignPtr a))+++zipComplex ::+ (Shape.C sh, Eq sh, Class.Real a) =>+ T shA sh a -> T shB sh a -> Vector sh (Complex a)+zipComplex xr xi =+ Array.unsafeCreateWithSize (shape xr) $ \n yPtr -> evalContT $ do+ liftIO $ Call.assert "zipComplex: shapes mismatch" (shape xr == shape xi)+ nPtr <- Call.cint n+ (xrPtr, incxrPtr) <- sliceArg xr+ (xiPtr, incxiPtr) <- sliceArg xi+ let yrPtr = realPtr yPtr+ incyPtr <- Call.cint 2+ liftIO $ do+ Blas.copy nPtr xrPtr incxrPtr yrPtr incyPtr+ Blas.copy nPtr xiPtr incxiPtr (advancePtr yrPtr 1) incyPtr++{- |+prop> forSliced complex_ $ \xs -> approxReal 1e-2 (VectorSlice.norm2 xs) $ let (xrs,xis) = VectorSlice.unzipComplex xs in sqrt $ VectorSlice.norm2Squared xrs + VectorSlice.norm2Squared xis+-}+unzipComplex ::+ (Shape.C sh, Class.Real a) =>+ T shA sh (Complex a) ->+ (T (shA,ComplexShape) sh a, T (shA,ComplexShape) sh a)+unzipComplex x = (realPart x, imaginaryPart x)
+ test/Main.hs view
@@ -0,0 +1,32 @@+module Main where++import qualified Test.Numeric.BLAS.Slice as Slice++import qualified Test.Float.Numeric.BLAS.Vector as VectorFloat+import qualified Test.Double.Numeric.BLAS.Vector as VectorDouble+import qualified Test.ComplexFloat.Numeric.BLAS.Vector as VectorComplexFloat+import qualified Test.ComplexDouble.Numeric.BLAS.Vector as VectorComplexDouble++import qualified Test.Float.Numeric.BLAS.Vector.Slice as SliceFloat+import qualified Test.Double.Numeric.BLAS.Vector.Slice as SliceDouble+import qualified Test.ComplexFloat.Numeric.BLAS.Vector.Slice as SliceComplexFloat+import qualified Test.ComplexDouble.Numeric.BLAS.Vector.Slice as SliceComplexDouble++import qualified Test.DocTest.Driver as DocTest++main :: IO ()+main = DocTest.run $ do+ DocTest.printLine "\nSlice"+ Slice.test+ DocTest.printLine "\nFloat"+ VectorFloat.test+ SliceFloat.test+ DocTest.printLine "\nDouble"+ VectorDouble.test+ SliceDouble.test+ DocTest.printLine "\nComplex Float"+ VectorComplexFloat.test+ SliceComplexFloat.test+ DocTest.printLine "\nComplex Double"+ VectorComplexDouble.test+ SliceComplexDouble.test
+ test/Test/ComplexDouble/Numeric/BLAS/Vector.hs view
@@ -0,0 +1,201 @@+-- Do not edit! Automatically created with doctest-extract from src/Numeric/BLAS/Vector.hs+{-# LINE 98 "src/Numeric/BLAS/Vector.hs" #-}++module Test.ComplexDouble.Numeric.BLAS.Vector where++import Test.DocTest.Base+import qualified Test.DocTest.Driver as DocTest++{-# LINE 99 "src/Numeric/BLAS/Vector.hs" #-}+import Test.ComplexDouble.Type (Number_)+import Test.Generator (genNumber)+import Test.Utility (approx)+import qualified Numeric.BLAS.Vector as Vector+import qualified Numeric.Netlib.Class as Class+import qualified Data.Array.Comfort.Shape as Shape+import qualified Data.Array.Comfort.Storable as Array+import qualified Data.List as List+import Numeric.BLAS.Vector ((+++), (|+|), (|-|))+import Numeric.BLAS.Scalar (RealOf, absolute, minusOne)+import Data.Array.Comfort.Storable (Array, (!))+import Data.Complex (Complex((:+)))+import Control.Applicative (liftA2)+import Control.Monad (replicateM)++import qualified Test.QuickCheck as QC+import Test.QuickCheck ((==>))++type Vector = Vector.Vector (Shape.ZeroBased Int)+type Real_ = RealOf Number_+type Complex_ = Complex Real_++maxElem :: Integer+maxElem = 10++maxDim :: Int+maxDim = 100++genVector :: (Class.Floating a) => Integer -> Int -> QC.Gen (Vector a)+genVector maxE dim =+ fmap (Vector.fromList (Shape.ZeroBased dim)) $+ replicateM dim $ genNumber maxE++isNonEmpty :: Shape.C sh => Array sh a -> Bool+isNonEmpty xs = Shape.size (Array.shape xs) > 0++newtype ArbRealVector = ArbRealVector (Vector Real_)+ deriving (Show)++instance QC.Arbitrary ArbRealVector where+ arbitrary = do+ dim <- QC.choose (0,maxDim)+ fmap ArbRealVector $ genVector maxElem dim+ shrink (ArbRealVector xs) =+ map (ArbRealVector . Vector.autoFromList) .+ QC.shrink . Vector.toList $ xs++newtype ArbVector = ArbVector (Vector Number_)+ deriving (Show)++instance QC.Arbitrary ArbVector where+ arbitrary = do+ dim <- QC.choose (0,maxDim)+ fmap ArbVector $ genVector maxElem dim+ shrink (ArbVector xs) =+ map (ArbVector . Vector.autoFromList) .+ QC.shrink . Vector.toList $ xs++data ArbVector2 =+ ArbVector2 (Vector Number_) (Vector Number_)+ deriving (Show)++instance QC.Arbitrary ArbVector2 where+ arbitrary = do+ dim <- QC.choose (0,maxDim)+ liftA2 ArbVector2+ (genVector maxElem dim) (genVector maxElem dim)++test :: DocTest.T ()+test = do+ DocTest.printPrefix "Numeric.BLAS.Vector:212: "+{-# LINE 212 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 212 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector xs) (ArbVector ys) (ArbVector zs) -> Vector.toList ((xs +++ ys) +++ zs) == Vector.toList (xs +++ (ys +++ zs)))+ DocTest.printPrefix "Numeric.BLAS.Vector:221: "+{-# LINE 221 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 221 "src/Numeric/BLAS/Vector.hs" #-}+ (QC.forAll (QC.choose (1,100)) $ \dim -> QC.forAll (QC.choose (0, dim-1)) $ \i -> QC.forAll (QC.choose (0, dim-1)) $ \j -> Vector.unit (Shape.ZeroBased dim) i == (Vector.swap i j (Vector.unit (Shape.ZeroBased dim) j) :: Vector Number_))+ DocTest.printPrefix "Numeric.BLAS.Vector:256: "+{-# LINE 256 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 256 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector2 xs ys) -> Vector.inner xs ys == Vector.dot (Vector.conjugate xs) ys)+ DocTest.printPrefix "Numeric.BLAS.Vector:309: "+{-# LINE 309 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 309 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector xs) -> Vector.sum xs == List.sum (Vector.toList xs))+ DocTest.printPrefix "Numeric.BLAS.Vector:369: "+{-# LINE 369 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 369 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector xs) -> Vector.normInf xs == List.maximum (0 : List.map absolute (Vector.toList xs)))+ DocTest.printPrefix "Numeric.BLAS.Vector:396: "+{-# LINE 396 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 396 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector xs) -> isNonEmpty xs ==> let (xi,xm) = Vector.argAbsMaximum xs in xs!xi == xm)+ DocTest.printPrefix "Numeric.BLAS.Vector:397: "+{-# LINE 397 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 397 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector xs) -> isNonEmpty xs ==> let (_xi,xm) = Vector.argAbsMaximum xs in List.all (\x -> absolute x <= absolute xm) $ Vector.toList xs)+ DocTest.printPrefix "Numeric.BLAS.Vector:398: "+{-# LINE 398 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 398 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector xs, ArbVector ys) -> isNonEmpty xs && isNonEmpty ys ==> let (_xi,xm) = Vector.argAbsMaximum xs; (_yi,ym) = Vector.argAbsMaximum ys; (zi,zm) = Vector.argAbsMaximum (xs+++ys) in case zi of Left _ -> xm==zm && absolute xm >= absolute ym; Right _ -> ym==zm && absolute xm < absolute ym)+ DocTest.printPrefix "Numeric.BLAS.Vector:393: "+{-# LINE 393 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.example+{-# LINE 393 "src/Numeric/BLAS/Vector.hs" #-}+ (Vector.argAbsMaximum $ Vector.autoFromList [1:+2, 3:+4, 5, 6 :: Complex_])+ [ExpectedLine [LineChunk "(3,6.0 :+ 0.0)"]]+ DocTest.printPrefix "Numeric.BLAS.Vector:445: "+{-# LINE 445 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 445 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbRealVector xs) -> isNonEmpty xs ==> Vector.argAbsMaximum xs == Vector.argAbs1Maximum xs)+ DocTest.printPrefix "Numeric.BLAS.Vector:442: "+{-# LINE 442 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.example+{-# LINE 442 "src/Numeric/BLAS/Vector.hs" #-}+ (Vector.argAbs1Maximum $ Vector.autoFromList [1:+2, 3:+4, 5, 6 :: Complex_])+ [ExpectedLine [LineChunk "(1,3.0 :+ 4.0)"]]+ DocTest.printPrefix "Numeric.BLAS.Vector:478: "+{-# LINE 478 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 478 "src/Numeric/BLAS/Vector.hs" #-}+ (QC.forAll (QC.choose (0,10)) $ \dim -> QC.forAll (genVector 3 dim) $ \xs -> approx 1e-2 (Vector.product xs) (List.product (Vector.toList (xs :: Vector Number_))))+ DocTest.printPrefix "Numeric.BLAS.Vector:488: "+{-# LINE 488 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 488 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbRealVector xs) -> isNonEmpty xs ==> Vector.minimum xs == List.minimum (Vector.toList xs))+ DocTest.printPrefix "Numeric.BLAS.Vector:489: "+{-# LINE 489 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 489 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbRealVector xs) -> isNonEmpty xs ==> Vector.maximum xs == List.maximum (Vector.toList xs))+ DocTest.printPrefix "Numeric.BLAS.Vector:490: "+{-# LINE 490 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 490 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbRealVector xs) -> isNonEmpty xs ==> - Vector.maximum xs == Vector.minimum (Vector.negate xs))+ DocTest.printPrefix "Numeric.BLAS.Vector:506: "+{-# LINE 506 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 506 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbRealVector xs) -> isNonEmpty xs ==> Vector.limits xs == Array.limits xs)+ DocTest.printPrefix "Numeric.BLAS.Vector:536: "+{-# LINE 536 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 536 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector xs) -> Vector.negate xs == Vector.scale minusOne xs)+ DocTest.printPrefix "Numeric.BLAS.Vector:537: "+{-# LINE 537 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 537 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector xs) -> Vector.scale 2 xs == xs |+| xs)+ DocTest.printPrefix "Numeric.BLAS.Vector:606: "+{-# LINE 606 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 606 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector2 xs ys) -> xs |+| ys == ys |+| xs)+ DocTest.printPrefix "Numeric.BLAS.Vector:607: "+{-# LINE 607 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 607 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector2 xs ys) -> xs == xs |-| ys |+| ys)+ DocTest.printPrefix "Numeric.BLAS.Vector:637: "+{-# LINE 637 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 637 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector xs) -> xs == Vector.negate (Vector.negate xs))+ DocTest.printPrefix "Numeric.BLAS.Vector:650: "+{-# LINE 650 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 650 "src/Numeric/BLAS/Vector.hs" #-}+ (QC.forAll (genNumber maxElem) $ \d (ArbVector xs) -> xs == Vector.raise (-d) (Vector.raise d xs))+ DocTest.printPrefix "Numeric.BLAS.Vector:667: "+{-# LINE 667 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 667 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector2 xs ys) -> Vector.mul xs ys == Vector.mul ys xs)+ DocTest.printPrefix "Numeric.BLAS.Vector:675: "+{-# LINE 675 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 675 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector2 xs ys) -> Vector.mulConj xs ys == Vector.mul (Vector.conjugate xs) ys)
+ test/Test/ComplexDouble/Numeric/BLAS/Vector/Slice.hs view
@@ -0,0 +1,271 @@+-- Do not edit! Automatically created with doctest-extract from src/Numeric/BLAS/Vector/Slice.hs+{-# LINE 87 "src/Numeric/BLAS/Vector/Slice.hs" #-}++{-# OPTIONS_GHC -XTypeOperators #-}+module Test.ComplexDouble.Numeric.BLAS.Vector.Slice where++import qualified Test.DocTest.Driver as DocTest++{-# LINE 89 "src/Numeric/BLAS/Vector/Slice.hs" #-}+import qualified Test.Slice as TestSlice+import Test.ComplexDouble.Type (Number_)+import Test.Generator (genNumber)+import Test.Utility (approx, approxReal)++import qualified Numeric.BLAS.Matrix.RowMajor as Matrix+import qualified Numeric.BLAS.Vector.Slice as VectorSlice+import qualified Numeric.BLAS.Vector as Vector+import qualified Numeric.BLAS.Slice as Slice+import qualified Numeric.BLAS.Scalar as Scalar+import qualified Numeric.Netlib.Class as Class+import qualified Data.Array.Comfort.Storable as Array+import qualified Data.Array.Comfort.Shape as Shape+import qualified Data.List as List+import Numeric.BLAS.Vector ((+++), (|+|))+import Numeric.BLAS.Scalar (RealOf, absolute, minusOne)+import Data.Array.Comfort.Shape ((::+)((::+)))+import Data.Tuple.HT (mapPair)+import Data.Complex (Complex)+import Control.Monad (replicateM)+import Control.Applicative (liftA2)++import qualified Test.QuickCheck as QC+import Test.QuickCheck ((==>))++type Real_ = RealOf Number_+type Complex_ = Complex Real_++maxElem :: Integer+maxElem = 10++maxDim :: Int+maxDim = 100++maxDim1 :: Int+maxDim1 = 10++type ShapeInt = Shape.ZeroBased Int+type Shape = ShapeInt ::+ (ShapeInt, ShapeInt) ::+ ShapeInt+type Vector = Vector.Vector Shape+type Sliced = VectorSlice.T Shape ShapeInt++genDim :: QC.Gen Int+genDim = QC.choose (0,maxDim)++genShapeDim :: Int -> QC.Gen Shape+genShapeDim numRows = do+ left <- fmap Shape.ZeroBased $ QC.choose (0,maxDim)+ right <- fmap Shape.ZeroBased $ QC.choose (0,maxDim)+ columns <- fmap Shape.ZeroBased $ QC.choose (1,maxDim1)+ return (left ::+ (Shape.ZeroBased numRows, columns) ::+ right)++genShape :: QC.Gen Shape+genShape = genShapeDim =<< QC.choose (0,maxDim1)++genVector ::+ (Shape.C sh, Class.Floating a) =>+ sh -> QC.Gen a -> QC.Gen (Vector.Vector sh a)+genVector shape genElem =+ fmap (Vector.fromList shape) $+ replicateM (Shape.size shape) genElem++forAll_ :: (Show a) => QC.Gen a -> (a -> QC.Property) -> QC.Property+forAll_ = QC.forAll++real_ :: QC.Gen Real_+real_ = genNumber maxElem++complex_ :: QC.Gen Complex_+complex_ = genNumber maxElem++number_ :: QC.Gen Number_+number_ = genNumber maxElem++isNonEmpty :: Shape.C sh => VectorSlice.T shA sh a -> Bool+isNonEmpty xs = Shape.size (VectorSlice.shape xs) > 0++takeColumn ::+ (Shape.Indexed sh1, Shape.C sh, Shape.C sh0, Shape.C sh2) =>+ Shape.Index sh1 ->+ Vector.Vector (sh0 ::+ (sh, sh1) ::+ sh2) a ->+ VectorSlice.T (sh0 ::+ (sh, sh1) ::+ sh2) sh a+takeColumn c =+ VectorSlice.slice+ (Slice.column c . Slice.left . Slice.right) . VectorSlice.fromVector++listFromSlice ::+ (Shape.C sh, Class.Floating a) => VectorSlice.T shA sh a -> [a]+listFromSlice = Vector.toList . VectorSlice.toVector++genSlicedDim ::+ (Class.Floating a) =>+ Int -> QC.Gen a -> QC.Gen (Int, Vector a)+genSlicedDim numRows genElem = do+ shape@(_::+(_rows,columns)::+_) <- genShapeDim numRows+ c <- QC.elements (Shape.indices columns)+ fmap ((,) c) $ genVector shape genElem++genSliced ::+ (Class.Floating a) =>+ QC.Gen a -> QC.Gen (Int, Vector a)+genSliced genElem = flip genSlicedDim genElem =<< genDim++shrinkSliced ::+ (Shape.C sh0, Shape.Indexed sh1, QC.Arbitrary a, Class.Floating a) =>+ (Shape.Index sh1,+ Vector.Vector (ShapeInt ::+ (sh0, sh1) ::+ ShapeInt) a) ->+ [(Shape.Index sh1,+ Vector.Vector (ShapeInt ::+ (sh0, sh1) ::+ ShapeInt) a)]+shrinkSliced (c,xs) =+ let xs0 = Vector.takeLeft xs in+ let xs1 = Vector.takeRight xs in+ let xs10 = Vector.takeLeft xs1 in+ let xs11 = Vector.takeRight xs1 in+ map (\(ysl,ysr) ->+ (c,+ Vector.autoFromList ysl +++ xs10 +++ Vector.autoFromList ysr)) $+ QC.shrink (Vector.toList xs0, Vector.toList xs11)++forSliced ::+ (QC.Testable prop, QC.Arbitrary a, Class.Floating a, Show a) =>+ QC.Gen a -> (Sliced a -> prop) -> QC.Property+forSliced genElem prop =+ QC.forAllShrink (genSliced genElem) shrinkSliced+ (prop . uncurry takeColumn)++genSliced2 ::+ (Class.Floating a) =>+ QC.Gen a -> QC.Gen ((Int, Vector a), (Int, Vector a))+genSliced2 genElem = do+ dim <- genDim+ liftA2 (,) (genSlicedDim dim genElem) (genSlicedDim dim genElem)++forSliced2 ::+ (QC.Testable prop, Class.Floating a, Show a) =>+ QC.Gen a -> (Sliced a -> Sliced a -> prop) -> QC.Property+forSliced2 genElem prop =+ QC.forAll (genSliced2 genElem)+ (uncurry prop . mapPair (uncurry takeColumn, uncurry takeColumn))++test :: DocTest.T ()+test = do+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:300: "+{-# LINE 300 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 300 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (QC.forAll genShape $ \shape@(_::+(_rows,columns)::+_) -> QC.forAll (QC.elements (Shape.indices columns)) $ \c -> QC.forAll (genVector shape number_) $ \xs -> VectorSlice.extract (Slice.column c . Slice.left . Slice.right) xs == Matrix.takeColumn c (Vector.takeLeft (Vector.takeRight (xs :: Vector Number_))))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:302: "+{-# LINE 302 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 302 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forAll_ (TestSlice.genShapeSelect 4 100) $ \(TestSlice.ShapeSelect sh select) -> QC.forAll (genVector sh number_) $ \xs -> case TestSlice.instantiate sh select of TestSlice.Extraction slice cut -> VectorSlice.extract slice xs == cut xs)+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:334: "+{-# LINE 334 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 334 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced2 number_ $ \xs ys -> VectorSlice.inner xs ys == Vector.dot (VectorSlice.conjugate xs) (VectorSlice.toVector ys))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:335: "+{-# LINE 335 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 335 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced number_ $ \xs -> VectorSlice.inner xs xs == Scalar.fromReal (VectorSlice.norm2Squared xs))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:418: "+{-# LINE 418 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 418 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced number_ $ \xs -> VectorSlice.sum xs == List.sum (listFromSlice xs))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:491: "+{-# LINE 491 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 491 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced number_ $ \xs -> VectorSlice.normInf xs == List.maximum (0 : List.map absolute (listFromSlice xs)))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:502: "+{-# LINE 502 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 502 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced number_ $ \xs -> VectorSlice.normInf1 xs == List.maximum (0 : List.map Scalar.norm1 (listFromSlice xs)))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:517: "+{-# LINE 517 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 517 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced number_ $ \xs -> isNonEmpty xs ==> let (xi,xm) = VectorSlice.argAbsMaximum xs in VectorSlice.access xs xi == xm)+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:518: "+{-# LINE 518 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 518 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced number_ $ \xs -> isNonEmpty xs ==> let (_xi,xm) = VectorSlice.argAbsMaximum xs in List.all (\x -> absolute x <= absolute xm) $ listFromSlice xs)+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:519: "+{-# LINE 519 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 519 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced number_ $ \xs -> forSliced number_ $ \ys -> isNonEmpty xs && isNonEmpty ys ==> let (_xi,xm) = VectorSlice.argAbsMaximum xs; (_yi,ym) = VectorSlice.argAbsMaximum ys; (zi,zm) = Vector.argAbsMaximum (VectorSlice.toVector xs +++ VectorSlice.toVector ys) in case zi of Left _ -> xm==zm && absolute xm >= absolute ym; Right _ -> ym==zm && absolute xm < absolute ym)+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:566: "+{-# LINE 566 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 566 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced real_ $ \xs -> isNonEmpty xs ==> VectorSlice.argAbsMaximum xs == VectorSlice.argAbs1Maximum xs)+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:591: "+{-# LINE 591 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 591 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (QC.forAll genShape $ \sh@(_::+(_rows,columns)::+_) -> QC.forAll (QC.elements (Shape.indices columns)) $ \c -> QC.forAll (genVector sh $ genNumber 3) $ \xt -> let xs = takeColumn c xt in approx 1e-2 (VectorSlice.product xs) (List.product (listFromSlice (xs :: Sliced Number_))))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:602: "+{-# LINE 602 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 602 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced real_ $ \xs -> isNonEmpty xs ==> VectorSlice.minimum xs == List.minimum (listFromSlice xs))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:603: "+{-# LINE 603 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 603 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced real_ $ \xs -> isNonEmpty xs ==> VectorSlice.maximum xs == List.maximum (listFromSlice xs))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:619: "+{-# LINE 619 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 619 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced real_ $ \xs -> isNonEmpty xs ==> VectorSlice.limits xs == Array.limits (VectorSlice.toVector xs))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:649: "+{-# LINE 649 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 649 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced number_ $ \xs -> VectorSlice.negate xs == VectorSlice.scale minusOne xs)+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:650: "+{-# LINE 650 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 650 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced number_ $ \xs -> VectorSlice.scale 2 xs == VectorSlice.add xs xs)+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:709: "+{-# LINE 709 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 709 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced2 number_ $ \xs ys -> VectorSlice.add xs ys == VectorSlice.add ys xs)+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:710: "+{-# LINE 710 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 710 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced2 number_ $ \xs ys -> VectorSlice.toVector xs == VectorSlice.sub xs ys |+| VectorSlice.toVector ys)+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:735: "+{-# LINE 735 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 735 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced number_ $ \xs -> VectorSlice.toVector xs == Vector.negate (VectorSlice.negate xs))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:748: "+{-# LINE 748 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 748 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (QC.forAll (genNumber maxElem) $ \d -> forSliced number_ $ \xs -> VectorSlice.toVector xs == Vector.raise (-d) (VectorSlice.raise d xs))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:765: "+{-# LINE 765 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 765 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced2 number_ $ \xs ys -> VectorSlice.mul xs ys == VectorSlice.mul ys xs)+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:773: "+{-# LINE 773 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 773 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced2 number_ $ \xs ys -> VectorSlice.mulConj xs ys == Vector.mul (VectorSlice.conjugate xs) (VectorSlice.toVector ys))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:897: "+{-# LINE 897 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 897 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced complex_ $ \xs -> approxReal 1e-2 (VectorSlice.norm2 xs) $ let (xrs,xis) = VectorSlice.unzipComplex xs in sqrt $ VectorSlice.norm2Squared xrs + VectorSlice.norm2Squared xis)
+ test/Test/ComplexDouble/Type.hs view
@@ -0,0 +1,5 @@+module Test.ComplexDouble.Type where++import Data.Complex (Complex)++type Number_ = Complex Double
+ test/Test/ComplexFloat/Numeric/BLAS/Vector.hs view
@@ -0,0 +1,201 @@+-- Do not edit! Automatically created with doctest-extract from src/Numeric/BLAS/Vector.hs+{-# LINE 98 "src/Numeric/BLAS/Vector.hs" #-}++module Test.ComplexFloat.Numeric.BLAS.Vector where++import Test.DocTest.Base+import qualified Test.DocTest.Driver as DocTest++{-# LINE 99 "src/Numeric/BLAS/Vector.hs" #-}+import Test.ComplexFloat.Type (Number_)+import Test.Generator (genNumber)+import Test.Utility (approx)+import qualified Numeric.BLAS.Vector as Vector+import qualified Numeric.Netlib.Class as Class+import qualified Data.Array.Comfort.Shape as Shape+import qualified Data.Array.Comfort.Storable as Array+import qualified Data.List as List+import Numeric.BLAS.Vector ((+++), (|+|), (|-|))+import Numeric.BLAS.Scalar (RealOf, absolute, minusOne)+import Data.Array.Comfort.Storable (Array, (!))+import Data.Complex (Complex((:+)))+import Control.Applicative (liftA2)+import Control.Monad (replicateM)++import qualified Test.QuickCheck as QC+import Test.QuickCheck ((==>))++type Vector = Vector.Vector (Shape.ZeroBased Int)+type Real_ = RealOf Number_+type Complex_ = Complex Real_++maxElem :: Integer+maxElem = 10++maxDim :: Int+maxDim = 100++genVector :: (Class.Floating a) => Integer -> Int -> QC.Gen (Vector a)+genVector maxE dim =+ fmap (Vector.fromList (Shape.ZeroBased dim)) $+ replicateM dim $ genNumber maxE++isNonEmpty :: Shape.C sh => Array sh a -> Bool+isNonEmpty xs = Shape.size (Array.shape xs) > 0++newtype ArbRealVector = ArbRealVector (Vector Real_)+ deriving (Show)++instance QC.Arbitrary ArbRealVector where+ arbitrary = do+ dim <- QC.choose (0,maxDim)+ fmap ArbRealVector $ genVector maxElem dim+ shrink (ArbRealVector xs) =+ map (ArbRealVector . Vector.autoFromList) .+ QC.shrink . Vector.toList $ xs++newtype ArbVector = ArbVector (Vector Number_)+ deriving (Show)++instance QC.Arbitrary ArbVector where+ arbitrary = do+ dim <- QC.choose (0,maxDim)+ fmap ArbVector $ genVector maxElem dim+ shrink (ArbVector xs) =+ map (ArbVector . Vector.autoFromList) .+ QC.shrink . Vector.toList $ xs++data ArbVector2 =+ ArbVector2 (Vector Number_) (Vector Number_)+ deriving (Show)++instance QC.Arbitrary ArbVector2 where+ arbitrary = do+ dim <- QC.choose (0,maxDim)+ liftA2 ArbVector2+ (genVector maxElem dim) (genVector maxElem dim)++test :: DocTest.T ()+test = do+ DocTest.printPrefix "Numeric.BLAS.Vector:212: "+{-# LINE 212 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 212 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector xs) (ArbVector ys) (ArbVector zs) -> Vector.toList ((xs +++ ys) +++ zs) == Vector.toList (xs +++ (ys +++ zs)))+ DocTest.printPrefix "Numeric.BLAS.Vector:221: "+{-# LINE 221 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 221 "src/Numeric/BLAS/Vector.hs" #-}+ (QC.forAll (QC.choose (1,100)) $ \dim -> QC.forAll (QC.choose (0, dim-1)) $ \i -> QC.forAll (QC.choose (0, dim-1)) $ \j -> Vector.unit (Shape.ZeroBased dim) i == (Vector.swap i j (Vector.unit (Shape.ZeroBased dim) j) :: Vector Number_))+ DocTest.printPrefix "Numeric.BLAS.Vector:256: "+{-# LINE 256 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 256 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector2 xs ys) -> Vector.inner xs ys == Vector.dot (Vector.conjugate xs) ys)+ DocTest.printPrefix "Numeric.BLAS.Vector:309: "+{-# LINE 309 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 309 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector xs) -> Vector.sum xs == List.sum (Vector.toList xs))+ DocTest.printPrefix "Numeric.BLAS.Vector:369: "+{-# LINE 369 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 369 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector xs) -> Vector.normInf xs == List.maximum (0 : List.map absolute (Vector.toList xs)))+ DocTest.printPrefix "Numeric.BLAS.Vector:396: "+{-# LINE 396 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 396 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector xs) -> isNonEmpty xs ==> let (xi,xm) = Vector.argAbsMaximum xs in xs!xi == xm)+ DocTest.printPrefix "Numeric.BLAS.Vector:397: "+{-# LINE 397 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 397 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector xs) -> isNonEmpty xs ==> let (_xi,xm) = Vector.argAbsMaximum xs in List.all (\x -> absolute x <= absolute xm) $ Vector.toList xs)+ DocTest.printPrefix "Numeric.BLAS.Vector:398: "+{-# LINE 398 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 398 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector xs, ArbVector ys) -> isNonEmpty xs && isNonEmpty ys ==> let (_xi,xm) = Vector.argAbsMaximum xs; (_yi,ym) = Vector.argAbsMaximum ys; (zi,zm) = Vector.argAbsMaximum (xs+++ys) in case zi of Left _ -> xm==zm && absolute xm >= absolute ym; Right _ -> ym==zm && absolute xm < absolute ym)+ DocTest.printPrefix "Numeric.BLAS.Vector:393: "+{-# LINE 393 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.example+{-# LINE 393 "src/Numeric/BLAS/Vector.hs" #-}+ (Vector.argAbsMaximum $ Vector.autoFromList [1:+2, 3:+4, 5, 6 :: Complex_])+ [ExpectedLine [LineChunk "(3,6.0 :+ 0.0)"]]+ DocTest.printPrefix "Numeric.BLAS.Vector:445: "+{-# LINE 445 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 445 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbRealVector xs) -> isNonEmpty xs ==> Vector.argAbsMaximum xs == Vector.argAbs1Maximum xs)+ DocTest.printPrefix "Numeric.BLAS.Vector:442: "+{-# LINE 442 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.example+{-# LINE 442 "src/Numeric/BLAS/Vector.hs" #-}+ (Vector.argAbs1Maximum $ Vector.autoFromList [1:+2, 3:+4, 5, 6 :: Complex_])+ [ExpectedLine [LineChunk "(1,3.0 :+ 4.0)"]]+ DocTest.printPrefix "Numeric.BLAS.Vector:478: "+{-# LINE 478 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 478 "src/Numeric/BLAS/Vector.hs" #-}+ (QC.forAll (QC.choose (0,10)) $ \dim -> QC.forAll (genVector 3 dim) $ \xs -> approx 1e-2 (Vector.product xs) (List.product (Vector.toList (xs :: Vector Number_))))+ DocTest.printPrefix "Numeric.BLAS.Vector:488: "+{-# LINE 488 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 488 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbRealVector xs) -> isNonEmpty xs ==> Vector.minimum xs == List.minimum (Vector.toList xs))+ DocTest.printPrefix "Numeric.BLAS.Vector:489: "+{-# LINE 489 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 489 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbRealVector xs) -> isNonEmpty xs ==> Vector.maximum xs == List.maximum (Vector.toList xs))+ DocTest.printPrefix "Numeric.BLAS.Vector:490: "+{-# LINE 490 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 490 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbRealVector xs) -> isNonEmpty xs ==> - Vector.maximum xs == Vector.minimum (Vector.negate xs))+ DocTest.printPrefix "Numeric.BLAS.Vector:506: "+{-# LINE 506 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 506 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbRealVector xs) -> isNonEmpty xs ==> Vector.limits xs == Array.limits xs)+ DocTest.printPrefix "Numeric.BLAS.Vector:536: "+{-# LINE 536 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 536 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector xs) -> Vector.negate xs == Vector.scale minusOne xs)+ DocTest.printPrefix "Numeric.BLAS.Vector:537: "+{-# LINE 537 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 537 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector xs) -> Vector.scale 2 xs == xs |+| xs)+ DocTest.printPrefix "Numeric.BLAS.Vector:606: "+{-# LINE 606 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 606 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector2 xs ys) -> xs |+| ys == ys |+| xs)+ DocTest.printPrefix "Numeric.BLAS.Vector:607: "+{-# LINE 607 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 607 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector2 xs ys) -> xs == xs |-| ys |+| ys)+ DocTest.printPrefix "Numeric.BLAS.Vector:637: "+{-# LINE 637 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 637 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector xs) -> xs == Vector.negate (Vector.negate xs))+ DocTest.printPrefix "Numeric.BLAS.Vector:650: "+{-# LINE 650 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 650 "src/Numeric/BLAS/Vector.hs" #-}+ (QC.forAll (genNumber maxElem) $ \d (ArbVector xs) -> xs == Vector.raise (-d) (Vector.raise d xs))+ DocTest.printPrefix "Numeric.BLAS.Vector:667: "+{-# LINE 667 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 667 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector2 xs ys) -> Vector.mul xs ys == Vector.mul ys xs)+ DocTest.printPrefix "Numeric.BLAS.Vector:675: "+{-# LINE 675 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 675 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector2 xs ys) -> Vector.mulConj xs ys == Vector.mul (Vector.conjugate xs) ys)
+ test/Test/ComplexFloat/Numeric/BLAS/Vector/Slice.hs view
@@ -0,0 +1,271 @@+-- Do not edit! Automatically created with doctest-extract from src/Numeric/BLAS/Vector/Slice.hs+{-# LINE 87 "src/Numeric/BLAS/Vector/Slice.hs" #-}++{-# OPTIONS_GHC -XTypeOperators #-}+module Test.ComplexFloat.Numeric.BLAS.Vector.Slice where++import qualified Test.DocTest.Driver as DocTest++{-# LINE 89 "src/Numeric/BLAS/Vector/Slice.hs" #-}+import qualified Test.Slice as TestSlice+import Test.ComplexFloat.Type (Number_)+import Test.Generator (genNumber)+import Test.Utility (approx, approxReal)++import qualified Numeric.BLAS.Matrix.RowMajor as Matrix+import qualified Numeric.BLAS.Vector.Slice as VectorSlice+import qualified Numeric.BLAS.Vector as Vector+import qualified Numeric.BLAS.Slice as Slice+import qualified Numeric.BLAS.Scalar as Scalar+import qualified Numeric.Netlib.Class as Class+import qualified Data.Array.Comfort.Storable as Array+import qualified Data.Array.Comfort.Shape as Shape+import qualified Data.List as List+import Numeric.BLAS.Vector ((+++), (|+|))+import Numeric.BLAS.Scalar (RealOf, absolute, minusOne)+import Data.Array.Comfort.Shape ((::+)((::+)))+import Data.Tuple.HT (mapPair)+import Data.Complex (Complex)+import Control.Monad (replicateM)+import Control.Applicative (liftA2)++import qualified Test.QuickCheck as QC+import Test.QuickCheck ((==>))++type Real_ = RealOf Number_+type Complex_ = Complex Real_++maxElem :: Integer+maxElem = 10++maxDim :: Int+maxDim = 100++maxDim1 :: Int+maxDim1 = 10++type ShapeInt = Shape.ZeroBased Int+type Shape = ShapeInt ::+ (ShapeInt, ShapeInt) ::+ ShapeInt+type Vector = Vector.Vector Shape+type Sliced = VectorSlice.T Shape ShapeInt++genDim :: QC.Gen Int+genDim = QC.choose (0,maxDim)++genShapeDim :: Int -> QC.Gen Shape+genShapeDim numRows = do+ left <- fmap Shape.ZeroBased $ QC.choose (0,maxDim)+ right <- fmap Shape.ZeroBased $ QC.choose (0,maxDim)+ columns <- fmap Shape.ZeroBased $ QC.choose (1,maxDim1)+ return (left ::+ (Shape.ZeroBased numRows, columns) ::+ right)++genShape :: QC.Gen Shape+genShape = genShapeDim =<< QC.choose (0,maxDim1)++genVector ::+ (Shape.C sh, Class.Floating a) =>+ sh -> QC.Gen a -> QC.Gen (Vector.Vector sh a)+genVector shape genElem =+ fmap (Vector.fromList shape) $+ replicateM (Shape.size shape) genElem++forAll_ :: (Show a) => QC.Gen a -> (a -> QC.Property) -> QC.Property+forAll_ = QC.forAll++real_ :: QC.Gen Real_+real_ = genNumber maxElem++complex_ :: QC.Gen Complex_+complex_ = genNumber maxElem++number_ :: QC.Gen Number_+number_ = genNumber maxElem++isNonEmpty :: Shape.C sh => VectorSlice.T shA sh a -> Bool+isNonEmpty xs = Shape.size (VectorSlice.shape xs) > 0++takeColumn ::+ (Shape.Indexed sh1, Shape.C sh, Shape.C sh0, Shape.C sh2) =>+ Shape.Index sh1 ->+ Vector.Vector (sh0 ::+ (sh, sh1) ::+ sh2) a ->+ VectorSlice.T (sh0 ::+ (sh, sh1) ::+ sh2) sh a+takeColumn c =+ VectorSlice.slice+ (Slice.column c . Slice.left . Slice.right) . VectorSlice.fromVector++listFromSlice ::+ (Shape.C sh, Class.Floating a) => VectorSlice.T shA sh a -> [a]+listFromSlice = Vector.toList . VectorSlice.toVector++genSlicedDim ::+ (Class.Floating a) =>+ Int -> QC.Gen a -> QC.Gen (Int, Vector a)+genSlicedDim numRows genElem = do+ shape@(_::+(_rows,columns)::+_) <- genShapeDim numRows+ c <- QC.elements (Shape.indices columns)+ fmap ((,) c) $ genVector shape genElem++genSliced ::+ (Class.Floating a) =>+ QC.Gen a -> QC.Gen (Int, Vector a)+genSliced genElem = flip genSlicedDim genElem =<< genDim++shrinkSliced ::+ (Shape.C sh0, Shape.Indexed sh1, QC.Arbitrary a, Class.Floating a) =>+ (Shape.Index sh1,+ Vector.Vector (ShapeInt ::+ (sh0, sh1) ::+ ShapeInt) a) ->+ [(Shape.Index sh1,+ Vector.Vector (ShapeInt ::+ (sh0, sh1) ::+ ShapeInt) a)]+shrinkSliced (c,xs) =+ let xs0 = Vector.takeLeft xs in+ let xs1 = Vector.takeRight xs in+ let xs10 = Vector.takeLeft xs1 in+ let xs11 = Vector.takeRight xs1 in+ map (\(ysl,ysr) ->+ (c,+ Vector.autoFromList ysl +++ xs10 +++ Vector.autoFromList ysr)) $+ QC.shrink (Vector.toList xs0, Vector.toList xs11)++forSliced ::+ (QC.Testable prop, QC.Arbitrary a, Class.Floating a, Show a) =>+ QC.Gen a -> (Sliced a -> prop) -> QC.Property+forSliced genElem prop =+ QC.forAllShrink (genSliced genElem) shrinkSliced+ (prop . uncurry takeColumn)++genSliced2 ::+ (Class.Floating a) =>+ QC.Gen a -> QC.Gen ((Int, Vector a), (Int, Vector a))+genSliced2 genElem = do+ dim <- genDim+ liftA2 (,) (genSlicedDim dim genElem) (genSlicedDim dim genElem)++forSliced2 ::+ (QC.Testable prop, Class.Floating a, Show a) =>+ QC.Gen a -> (Sliced a -> Sliced a -> prop) -> QC.Property+forSliced2 genElem prop =+ QC.forAll (genSliced2 genElem)+ (uncurry prop . mapPair (uncurry takeColumn, uncurry takeColumn))++test :: DocTest.T ()+test = do+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:300: "+{-# LINE 300 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 300 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (QC.forAll genShape $ \shape@(_::+(_rows,columns)::+_) -> QC.forAll (QC.elements (Shape.indices columns)) $ \c -> QC.forAll (genVector shape number_) $ \xs -> VectorSlice.extract (Slice.column c . Slice.left . Slice.right) xs == Matrix.takeColumn c (Vector.takeLeft (Vector.takeRight (xs :: Vector Number_))))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:302: "+{-# LINE 302 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 302 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forAll_ (TestSlice.genShapeSelect 4 100) $ \(TestSlice.ShapeSelect sh select) -> QC.forAll (genVector sh number_) $ \xs -> case TestSlice.instantiate sh select of TestSlice.Extraction slice cut -> VectorSlice.extract slice xs == cut xs)+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:334: "+{-# LINE 334 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 334 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced2 number_ $ \xs ys -> VectorSlice.inner xs ys == Vector.dot (VectorSlice.conjugate xs) (VectorSlice.toVector ys))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:335: "+{-# LINE 335 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 335 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced number_ $ \xs -> VectorSlice.inner xs xs == Scalar.fromReal (VectorSlice.norm2Squared xs))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:418: "+{-# LINE 418 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 418 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced number_ $ \xs -> VectorSlice.sum xs == List.sum (listFromSlice xs))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:491: "+{-# LINE 491 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 491 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced number_ $ \xs -> VectorSlice.normInf xs == List.maximum (0 : List.map absolute (listFromSlice xs)))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:502: "+{-# LINE 502 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 502 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced number_ $ \xs -> VectorSlice.normInf1 xs == List.maximum (0 : List.map Scalar.norm1 (listFromSlice xs)))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:517: "+{-# LINE 517 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 517 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced number_ $ \xs -> isNonEmpty xs ==> let (xi,xm) = VectorSlice.argAbsMaximum xs in VectorSlice.access xs xi == xm)+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:518: "+{-# LINE 518 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 518 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced number_ $ \xs -> isNonEmpty xs ==> let (_xi,xm) = VectorSlice.argAbsMaximum xs in List.all (\x -> absolute x <= absolute xm) $ listFromSlice xs)+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:519: "+{-# LINE 519 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 519 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced number_ $ \xs -> forSliced number_ $ \ys -> isNonEmpty xs && isNonEmpty ys ==> let (_xi,xm) = VectorSlice.argAbsMaximum xs; (_yi,ym) = VectorSlice.argAbsMaximum ys; (zi,zm) = Vector.argAbsMaximum (VectorSlice.toVector xs +++ VectorSlice.toVector ys) in case zi of Left _ -> xm==zm && absolute xm >= absolute ym; Right _ -> ym==zm && absolute xm < absolute ym)+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:566: "+{-# LINE 566 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 566 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced real_ $ \xs -> isNonEmpty xs ==> VectorSlice.argAbsMaximum xs == VectorSlice.argAbs1Maximum xs)+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:591: "+{-# LINE 591 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 591 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (QC.forAll genShape $ \sh@(_::+(_rows,columns)::+_) -> QC.forAll (QC.elements (Shape.indices columns)) $ \c -> QC.forAll (genVector sh $ genNumber 3) $ \xt -> let xs = takeColumn c xt in approx 1e-2 (VectorSlice.product xs) (List.product (listFromSlice (xs :: Sliced Number_))))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:602: "+{-# LINE 602 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 602 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced real_ $ \xs -> isNonEmpty xs ==> VectorSlice.minimum xs == List.minimum (listFromSlice xs))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:603: "+{-# LINE 603 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 603 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced real_ $ \xs -> isNonEmpty xs ==> VectorSlice.maximum xs == List.maximum (listFromSlice xs))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:619: "+{-# LINE 619 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 619 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced real_ $ \xs -> isNonEmpty xs ==> VectorSlice.limits xs == Array.limits (VectorSlice.toVector xs))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:649: "+{-# LINE 649 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 649 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced number_ $ \xs -> VectorSlice.negate xs == VectorSlice.scale minusOne xs)+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:650: "+{-# LINE 650 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 650 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced number_ $ \xs -> VectorSlice.scale 2 xs == VectorSlice.add xs xs)+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:709: "+{-# LINE 709 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 709 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced2 number_ $ \xs ys -> VectorSlice.add xs ys == VectorSlice.add ys xs)+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:710: "+{-# LINE 710 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 710 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced2 number_ $ \xs ys -> VectorSlice.toVector xs == VectorSlice.sub xs ys |+| VectorSlice.toVector ys)+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:735: "+{-# LINE 735 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 735 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced number_ $ \xs -> VectorSlice.toVector xs == Vector.negate (VectorSlice.negate xs))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:748: "+{-# LINE 748 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 748 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (QC.forAll (genNumber maxElem) $ \d -> forSliced number_ $ \xs -> VectorSlice.toVector xs == Vector.raise (-d) (VectorSlice.raise d xs))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:765: "+{-# LINE 765 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 765 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced2 number_ $ \xs ys -> VectorSlice.mul xs ys == VectorSlice.mul ys xs)+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:773: "+{-# LINE 773 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 773 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced2 number_ $ \xs ys -> VectorSlice.mulConj xs ys == Vector.mul (VectorSlice.conjugate xs) (VectorSlice.toVector ys))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:897: "+{-# LINE 897 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 897 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced complex_ $ \xs -> approxReal 1e-2 (VectorSlice.norm2 xs) $ let (xrs,xis) = VectorSlice.unzipComplex xs in sqrt $ VectorSlice.norm2Squared xrs + VectorSlice.norm2Squared xis)
+ test/Test/ComplexFloat/Type.hs view
@@ -0,0 +1,5 @@+module Test.ComplexFloat.Type where++import Data.Complex (Complex)++type Number_ = Complex Float
+ test/Test/Double/Numeric/BLAS/Vector.hs view
@@ -0,0 +1,201 @@+-- Do not edit! Automatically created with doctest-extract from src/Numeric/BLAS/Vector.hs+{-# LINE 98 "src/Numeric/BLAS/Vector.hs" #-}++module Test.Double.Numeric.BLAS.Vector where++import Test.DocTest.Base+import qualified Test.DocTest.Driver as DocTest++{-# LINE 99 "src/Numeric/BLAS/Vector.hs" #-}+import Test.Double.Type (Number_)+import Test.Generator (genNumber)+import Test.Utility (approx)+import qualified Numeric.BLAS.Vector as Vector+import qualified Numeric.Netlib.Class as Class+import qualified Data.Array.Comfort.Shape as Shape+import qualified Data.Array.Comfort.Storable as Array+import qualified Data.List as List+import Numeric.BLAS.Vector ((+++), (|+|), (|-|))+import Numeric.BLAS.Scalar (RealOf, absolute, minusOne)+import Data.Array.Comfort.Storable (Array, (!))+import Data.Complex (Complex((:+)))+import Control.Applicative (liftA2)+import Control.Monad (replicateM)++import qualified Test.QuickCheck as QC+import Test.QuickCheck ((==>))++type Vector = Vector.Vector (Shape.ZeroBased Int)+type Real_ = RealOf Number_+type Complex_ = Complex Real_++maxElem :: Integer+maxElem = 10++maxDim :: Int+maxDim = 100++genVector :: (Class.Floating a) => Integer -> Int -> QC.Gen (Vector a)+genVector maxE dim =+ fmap (Vector.fromList (Shape.ZeroBased dim)) $+ replicateM dim $ genNumber maxE++isNonEmpty :: Shape.C sh => Array sh a -> Bool+isNonEmpty xs = Shape.size (Array.shape xs) > 0++newtype ArbRealVector = ArbRealVector (Vector Real_)+ deriving (Show)++instance QC.Arbitrary ArbRealVector where+ arbitrary = do+ dim <- QC.choose (0,maxDim)+ fmap ArbRealVector $ genVector maxElem dim+ shrink (ArbRealVector xs) =+ map (ArbRealVector . Vector.autoFromList) .+ QC.shrink . Vector.toList $ xs++newtype ArbVector = ArbVector (Vector Number_)+ deriving (Show)++instance QC.Arbitrary ArbVector where+ arbitrary = do+ dim <- QC.choose (0,maxDim)+ fmap ArbVector $ genVector maxElem dim+ shrink (ArbVector xs) =+ map (ArbVector . Vector.autoFromList) .+ QC.shrink . Vector.toList $ xs++data ArbVector2 =+ ArbVector2 (Vector Number_) (Vector Number_)+ deriving (Show)++instance QC.Arbitrary ArbVector2 where+ arbitrary = do+ dim <- QC.choose (0,maxDim)+ liftA2 ArbVector2+ (genVector maxElem dim) (genVector maxElem dim)++test :: DocTest.T ()+test = do+ DocTest.printPrefix "Numeric.BLAS.Vector:212: "+{-# LINE 212 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 212 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector xs) (ArbVector ys) (ArbVector zs) -> Vector.toList ((xs +++ ys) +++ zs) == Vector.toList (xs +++ (ys +++ zs)))+ DocTest.printPrefix "Numeric.BLAS.Vector:221: "+{-# LINE 221 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 221 "src/Numeric/BLAS/Vector.hs" #-}+ (QC.forAll (QC.choose (1,100)) $ \dim -> QC.forAll (QC.choose (0, dim-1)) $ \i -> QC.forAll (QC.choose (0, dim-1)) $ \j -> Vector.unit (Shape.ZeroBased dim) i == (Vector.swap i j (Vector.unit (Shape.ZeroBased dim) j) :: Vector Number_))+ DocTest.printPrefix "Numeric.BLAS.Vector:256: "+{-# LINE 256 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 256 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector2 xs ys) -> Vector.inner xs ys == Vector.dot (Vector.conjugate xs) ys)+ DocTest.printPrefix "Numeric.BLAS.Vector:309: "+{-# LINE 309 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 309 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector xs) -> Vector.sum xs == List.sum (Vector.toList xs))+ DocTest.printPrefix "Numeric.BLAS.Vector:369: "+{-# LINE 369 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 369 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector xs) -> Vector.normInf xs == List.maximum (0 : List.map absolute (Vector.toList xs)))+ DocTest.printPrefix "Numeric.BLAS.Vector:396: "+{-# LINE 396 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 396 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector xs) -> isNonEmpty xs ==> let (xi,xm) = Vector.argAbsMaximum xs in xs!xi == xm)+ DocTest.printPrefix "Numeric.BLAS.Vector:397: "+{-# LINE 397 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 397 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector xs) -> isNonEmpty xs ==> let (_xi,xm) = Vector.argAbsMaximum xs in List.all (\x -> absolute x <= absolute xm) $ Vector.toList xs)+ DocTest.printPrefix "Numeric.BLAS.Vector:398: "+{-# LINE 398 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 398 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector xs, ArbVector ys) -> isNonEmpty xs && isNonEmpty ys ==> let (_xi,xm) = Vector.argAbsMaximum xs; (_yi,ym) = Vector.argAbsMaximum ys; (zi,zm) = Vector.argAbsMaximum (xs+++ys) in case zi of Left _ -> xm==zm && absolute xm >= absolute ym; Right _ -> ym==zm && absolute xm < absolute ym)+ DocTest.printPrefix "Numeric.BLAS.Vector:393: "+{-# LINE 393 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.example+{-# LINE 393 "src/Numeric/BLAS/Vector.hs" #-}+ (Vector.argAbsMaximum $ Vector.autoFromList [1:+2, 3:+4, 5, 6 :: Complex_])+ [ExpectedLine [LineChunk "(3,6.0 :+ 0.0)"]]+ DocTest.printPrefix "Numeric.BLAS.Vector:445: "+{-# LINE 445 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 445 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbRealVector xs) -> isNonEmpty xs ==> Vector.argAbsMaximum xs == Vector.argAbs1Maximum xs)+ DocTest.printPrefix "Numeric.BLAS.Vector:442: "+{-# LINE 442 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.example+{-# LINE 442 "src/Numeric/BLAS/Vector.hs" #-}+ (Vector.argAbs1Maximum $ Vector.autoFromList [1:+2, 3:+4, 5, 6 :: Complex_])+ [ExpectedLine [LineChunk "(1,3.0 :+ 4.0)"]]+ DocTest.printPrefix "Numeric.BLAS.Vector:478: "+{-# LINE 478 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 478 "src/Numeric/BLAS/Vector.hs" #-}+ (QC.forAll (QC.choose (0,10)) $ \dim -> QC.forAll (genVector 3 dim) $ \xs -> approx 1e-2 (Vector.product xs) (List.product (Vector.toList (xs :: Vector Number_))))+ DocTest.printPrefix "Numeric.BLAS.Vector:488: "+{-# LINE 488 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 488 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbRealVector xs) -> isNonEmpty xs ==> Vector.minimum xs == List.minimum (Vector.toList xs))+ DocTest.printPrefix "Numeric.BLAS.Vector:489: "+{-# LINE 489 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 489 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbRealVector xs) -> isNonEmpty xs ==> Vector.maximum xs == List.maximum (Vector.toList xs))+ DocTest.printPrefix "Numeric.BLAS.Vector:490: "+{-# LINE 490 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 490 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbRealVector xs) -> isNonEmpty xs ==> - Vector.maximum xs == Vector.minimum (Vector.negate xs))+ DocTest.printPrefix "Numeric.BLAS.Vector:506: "+{-# LINE 506 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 506 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbRealVector xs) -> isNonEmpty xs ==> Vector.limits xs == Array.limits xs)+ DocTest.printPrefix "Numeric.BLAS.Vector:536: "+{-# LINE 536 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 536 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector xs) -> Vector.negate xs == Vector.scale minusOne xs)+ DocTest.printPrefix "Numeric.BLAS.Vector:537: "+{-# LINE 537 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 537 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector xs) -> Vector.scale 2 xs == xs |+| xs)+ DocTest.printPrefix "Numeric.BLAS.Vector:606: "+{-# LINE 606 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 606 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector2 xs ys) -> xs |+| ys == ys |+| xs)+ DocTest.printPrefix "Numeric.BLAS.Vector:607: "+{-# LINE 607 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 607 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector2 xs ys) -> xs == xs |-| ys |+| ys)+ DocTest.printPrefix "Numeric.BLAS.Vector:637: "+{-# LINE 637 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 637 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector xs) -> xs == Vector.negate (Vector.negate xs))+ DocTest.printPrefix "Numeric.BLAS.Vector:650: "+{-# LINE 650 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 650 "src/Numeric/BLAS/Vector.hs" #-}+ (QC.forAll (genNumber maxElem) $ \d (ArbVector xs) -> xs == Vector.raise (-d) (Vector.raise d xs))+ DocTest.printPrefix "Numeric.BLAS.Vector:667: "+{-# LINE 667 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 667 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector2 xs ys) -> Vector.mul xs ys == Vector.mul ys xs)+ DocTest.printPrefix "Numeric.BLAS.Vector:675: "+{-# LINE 675 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 675 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector2 xs ys) -> Vector.mulConj xs ys == Vector.mul (Vector.conjugate xs) ys)
+ test/Test/Double/Numeric/BLAS/Vector/Slice.hs view
@@ -0,0 +1,271 @@+-- Do not edit! Automatically created with doctest-extract from src/Numeric/BLAS/Vector/Slice.hs+{-# LINE 87 "src/Numeric/BLAS/Vector/Slice.hs" #-}++{-# OPTIONS_GHC -XTypeOperators #-}+module Test.Double.Numeric.BLAS.Vector.Slice where++import qualified Test.DocTest.Driver as DocTest++{-# LINE 89 "src/Numeric/BLAS/Vector/Slice.hs" #-}+import qualified Test.Slice as TestSlice+import Test.Double.Type (Number_)+import Test.Generator (genNumber)+import Test.Utility (approx, approxReal)++import qualified Numeric.BLAS.Matrix.RowMajor as Matrix+import qualified Numeric.BLAS.Vector.Slice as VectorSlice+import qualified Numeric.BLAS.Vector as Vector+import qualified Numeric.BLAS.Slice as Slice+import qualified Numeric.BLAS.Scalar as Scalar+import qualified Numeric.Netlib.Class as Class+import qualified Data.Array.Comfort.Storable as Array+import qualified Data.Array.Comfort.Shape as Shape+import qualified Data.List as List+import Numeric.BLAS.Vector ((+++), (|+|))+import Numeric.BLAS.Scalar (RealOf, absolute, minusOne)+import Data.Array.Comfort.Shape ((::+)((::+)))+import Data.Tuple.HT (mapPair)+import Data.Complex (Complex)+import Control.Monad (replicateM)+import Control.Applicative (liftA2)++import qualified Test.QuickCheck as QC+import Test.QuickCheck ((==>))++type Real_ = RealOf Number_+type Complex_ = Complex Real_++maxElem :: Integer+maxElem = 10++maxDim :: Int+maxDim = 100++maxDim1 :: Int+maxDim1 = 10++type ShapeInt = Shape.ZeroBased Int+type Shape = ShapeInt ::+ (ShapeInt, ShapeInt) ::+ ShapeInt+type Vector = Vector.Vector Shape+type Sliced = VectorSlice.T Shape ShapeInt++genDim :: QC.Gen Int+genDim = QC.choose (0,maxDim)++genShapeDim :: Int -> QC.Gen Shape+genShapeDim numRows = do+ left <- fmap Shape.ZeroBased $ QC.choose (0,maxDim)+ right <- fmap Shape.ZeroBased $ QC.choose (0,maxDim)+ columns <- fmap Shape.ZeroBased $ QC.choose (1,maxDim1)+ return (left ::+ (Shape.ZeroBased numRows, columns) ::+ right)++genShape :: QC.Gen Shape+genShape = genShapeDim =<< QC.choose (0,maxDim1)++genVector ::+ (Shape.C sh, Class.Floating a) =>+ sh -> QC.Gen a -> QC.Gen (Vector.Vector sh a)+genVector shape genElem =+ fmap (Vector.fromList shape) $+ replicateM (Shape.size shape) genElem++forAll_ :: (Show a) => QC.Gen a -> (a -> QC.Property) -> QC.Property+forAll_ = QC.forAll++real_ :: QC.Gen Real_+real_ = genNumber maxElem++complex_ :: QC.Gen Complex_+complex_ = genNumber maxElem++number_ :: QC.Gen Number_+number_ = genNumber maxElem++isNonEmpty :: Shape.C sh => VectorSlice.T shA sh a -> Bool+isNonEmpty xs = Shape.size (VectorSlice.shape xs) > 0++takeColumn ::+ (Shape.Indexed sh1, Shape.C sh, Shape.C sh0, Shape.C sh2) =>+ Shape.Index sh1 ->+ Vector.Vector (sh0 ::+ (sh, sh1) ::+ sh2) a ->+ VectorSlice.T (sh0 ::+ (sh, sh1) ::+ sh2) sh a+takeColumn c =+ VectorSlice.slice+ (Slice.column c . Slice.left . Slice.right) . VectorSlice.fromVector++listFromSlice ::+ (Shape.C sh, Class.Floating a) => VectorSlice.T shA sh a -> [a]+listFromSlice = Vector.toList . VectorSlice.toVector++genSlicedDim ::+ (Class.Floating a) =>+ Int -> QC.Gen a -> QC.Gen (Int, Vector a)+genSlicedDim numRows genElem = do+ shape@(_::+(_rows,columns)::+_) <- genShapeDim numRows+ c <- QC.elements (Shape.indices columns)+ fmap ((,) c) $ genVector shape genElem++genSliced ::+ (Class.Floating a) =>+ QC.Gen a -> QC.Gen (Int, Vector a)+genSliced genElem = flip genSlicedDim genElem =<< genDim++shrinkSliced ::+ (Shape.C sh0, Shape.Indexed sh1, QC.Arbitrary a, Class.Floating a) =>+ (Shape.Index sh1,+ Vector.Vector (ShapeInt ::+ (sh0, sh1) ::+ ShapeInt) a) ->+ [(Shape.Index sh1,+ Vector.Vector (ShapeInt ::+ (sh0, sh1) ::+ ShapeInt) a)]+shrinkSliced (c,xs) =+ let xs0 = Vector.takeLeft xs in+ let xs1 = Vector.takeRight xs in+ let xs10 = Vector.takeLeft xs1 in+ let xs11 = Vector.takeRight xs1 in+ map (\(ysl,ysr) ->+ (c,+ Vector.autoFromList ysl +++ xs10 +++ Vector.autoFromList ysr)) $+ QC.shrink (Vector.toList xs0, Vector.toList xs11)++forSliced ::+ (QC.Testable prop, QC.Arbitrary a, Class.Floating a, Show a) =>+ QC.Gen a -> (Sliced a -> prop) -> QC.Property+forSliced genElem prop =+ QC.forAllShrink (genSliced genElem) shrinkSliced+ (prop . uncurry takeColumn)++genSliced2 ::+ (Class.Floating a) =>+ QC.Gen a -> QC.Gen ((Int, Vector a), (Int, Vector a))+genSliced2 genElem = do+ dim <- genDim+ liftA2 (,) (genSlicedDim dim genElem) (genSlicedDim dim genElem)++forSliced2 ::+ (QC.Testable prop, Class.Floating a, Show a) =>+ QC.Gen a -> (Sliced a -> Sliced a -> prop) -> QC.Property+forSliced2 genElem prop =+ QC.forAll (genSliced2 genElem)+ (uncurry prop . mapPair (uncurry takeColumn, uncurry takeColumn))++test :: DocTest.T ()+test = do+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:300: "+{-# LINE 300 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 300 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (QC.forAll genShape $ \shape@(_::+(_rows,columns)::+_) -> QC.forAll (QC.elements (Shape.indices columns)) $ \c -> QC.forAll (genVector shape number_) $ \xs -> VectorSlice.extract (Slice.column c . Slice.left . Slice.right) xs == Matrix.takeColumn c (Vector.takeLeft (Vector.takeRight (xs :: Vector Number_))))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:302: "+{-# LINE 302 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 302 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forAll_ (TestSlice.genShapeSelect 4 100) $ \(TestSlice.ShapeSelect sh select) -> QC.forAll (genVector sh number_) $ \xs -> case TestSlice.instantiate sh select of TestSlice.Extraction slice cut -> VectorSlice.extract slice xs == cut xs)+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:334: "+{-# LINE 334 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 334 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced2 number_ $ \xs ys -> VectorSlice.inner xs ys == Vector.dot (VectorSlice.conjugate xs) (VectorSlice.toVector ys))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:335: "+{-# LINE 335 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 335 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced number_ $ \xs -> VectorSlice.inner xs xs == Scalar.fromReal (VectorSlice.norm2Squared xs))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:418: "+{-# LINE 418 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 418 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced number_ $ \xs -> VectorSlice.sum xs == List.sum (listFromSlice xs))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:491: "+{-# LINE 491 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 491 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced number_ $ \xs -> VectorSlice.normInf xs == List.maximum (0 : List.map absolute (listFromSlice xs)))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:502: "+{-# LINE 502 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 502 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced number_ $ \xs -> VectorSlice.normInf1 xs == List.maximum (0 : List.map Scalar.norm1 (listFromSlice xs)))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:517: "+{-# LINE 517 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 517 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced number_ $ \xs -> isNonEmpty xs ==> let (xi,xm) = VectorSlice.argAbsMaximum xs in VectorSlice.access xs xi == xm)+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:518: "+{-# LINE 518 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 518 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced number_ $ \xs -> isNonEmpty xs ==> let (_xi,xm) = VectorSlice.argAbsMaximum xs in List.all (\x -> absolute x <= absolute xm) $ listFromSlice xs)+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:519: "+{-# LINE 519 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 519 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced number_ $ \xs -> forSliced number_ $ \ys -> isNonEmpty xs && isNonEmpty ys ==> let (_xi,xm) = VectorSlice.argAbsMaximum xs; (_yi,ym) = VectorSlice.argAbsMaximum ys; (zi,zm) = Vector.argAbsMaximum (VectorSlice.toVector xs +++ VectorSlice.toVector ys) in case zi of Left _ -> xm==zm && absolute xm >= absolute ym; Right _ -> ym==zm && absolute xm < absolute ym)+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:566: "+{-# LINE 566 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 566 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced real_ $ \xs -> isNonEmpty xs ==> VectorSlice.argAbsMaximum xs == VectorSlice.argAbs1Maximum xs)+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:591: "+{-# LINE 591 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 591 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (QC.forAll genShape $ \sh@(_::+(_rows,columns)::+_) -> QC.forAll (QC.elements (Shape.indices columns)) $ \c -> QC.forAll (genVector sh $ genNumber 3) $ \xt -> let xs = takeColumn c xt in approx 1e-2 (VectorSlice.product xs) (List.product (listFromSlice (xs :: Sliced Number_))))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:602: "+{-# LINE 602 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 602 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced real_ $ \xs -> isNonEmpty xs ==> VectorSlice.minimum xs == List.minimum (listFromSlice xs))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:603: "+{-# LINE 603 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 603 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced real_ $ \xs -> isNonEmpty xs ==> VectorSlice.maximum xs == List.maximum (listFromSlice xs))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:619: "+{-# LINE 619 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 619 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced real_ $ \xs -> isNonEmpty xs ==> VectorSlice.limits xs == Array.limits (VectorSlice.toVector xs))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:649: "+{-# LINE 649 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 649 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced number_ $ \xs -> VectorSlice.negate xs == VectorSlice.scale minusOne xs)+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:650: "+{-# LINE 650 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 650 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced number_ $ \xs -> VectorSlice.scale 2 xs == VectorSlice.add xs xs)+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:709: "+{-# LINE 709 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 709 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced2 number_ $ \xs ys -> VectorSlice.add xs ys == VectorSlice.add ys xs)+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:710: "+{-# LINE 710 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 710 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced2 number_ $ \xs ys -> VectorSlice.toVector xs == VectorSlice.sub xs ys |+| VectorSlice.toVector ys)+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:735: "+{-# LINE 735 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 735 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced number_ $ \xs -> VectorSlice.toVector xs == Vector.negate (VectorSlice.negate xs))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:748: "+{-# LINE 748 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 748 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (QC.forAll (genNumber maxElem) $ \d -> forSliced number_ $ \xs -> VectorSlice.toVector xs == Vector.raise (-d) (VectorSlice.raise d xs))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:765: "+{-# LINE 765 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 765 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced2 number_ $ \xs ys -> VectorSlice.mul xs ys == VectorSlice.mul ys xs)+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:773: "+{-# LINE 773 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 773 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced2 number_ $ \xs ys -> VectorSlice.mulConj xs ys == Vector.mul (VectorSlice.conjugate xs) (VectorSlice.toVector ys))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:897: "+{-# LINE 897 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 897 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced complex_ $ \xs -> approxReal 1e-2 (VectorSlice.norm2 xs) $ let (xrs,xis) = VectorSlice.unzipComplex xs in sqrt $ VectorSlice.norm2Squared xrs + VectorSlice.norm2Squared xis)
+ test/Test/Double/Type.hs view
@@ -0,0 +1,3 @@+module Test.Double.Type where++type Number_ = Double
+ test/Test/Float/Numeric/BLAS/Vector.hs view
@@ -0,0 +1,201 @@+-- Do not edit! Automatically created with doctest-extract from src/Numeric/BLAS/Vector.hs+{-# LINE 98 "src/Numeric/BLAS/Vector.hs" #-}++module Test.Float.Numeric.BLAS.Vector where++import Test.DocTest.Base+import qualified Test.DocTest.Driver as DocTest++{-# LINE 99 "src/Numeric/BLAS/Vector.hs" #-}+import Test.Float.Type (Number_)+import Test.Generator (genNumber)+import Test.Utility (approx)+import qualified Numeric.BLAS.Vector as Vector+import qualified Numeric.Netlib.Class as Class+import qualified Data.Array.Comfort.Shape as Shape+import qualified Data.Array.Comfort.Storable as Array+import qualified Data.List as List+import Numeric.BLAS.Vector ((+++), (|+|), (|-|))+import Numeric.BLAS.Scalar (RealOf, absolute, minusOne)+import Data.Array.Comfort.Storable (Array, (!))+import Data.Complex (Complex((:+)))+import Control.Applicative (liftA2)+import Control.Monad (replicateM)++import qualified Test.QuickCheck as QC+import Test.QuickCheck ((==>))++type Vector = Vector.Vector (Shape.ZeroBased Int)+type Real_ = RealOf Number_+type Complex_ = Complex Real_++maxElem :: Integer+maxElem = 10++maxDim :: Int+maxDim = 100++genVector :: (Class.Floating a) => Integer -> Int -> QC.Gen (Vector a)+genVector maxE dim =+ fmap (Vector.fromList (Shape.ZeroBased dim)) $+ replicateM dim $ genNumber maxE++isNonEmpty :: Shape.C sh => Array sh a -> Bool+isNonEmpty xs = Shape.size (Array.shape xs) > 0++newtype ArbRealVector = ArbRealVector (Vector Real_)+ deriving (Show)++instance QC.Arbitrary ArbRealVector where+ arbitrary = do+ dim <- QC.choose (0,maxDim)+ fmap ArbRealVector $ genVector maxElem dim+ shrink (ArbRealVector xs) =+ map (ArbRealVector . Vector.autoFromList) .+ QC.shrink . Vector.toList $ xs++newtype ArbVector = ArbVector (Vector Number_)+ deriving (Show)++instance QC.Arbitrary ArbVector where+ arbitrary = do+ dim <- QC.choose (0,maxDim)+ fmap ArbVector $ genVector maxElem dim+ shrink (ArbVector xs) =+ map (ArbVector . Vector.autoFromList) .+ QC.shrink . Vector.toList $ xs++data ArbVector2 =+ ArbVector2 (Vector Number_) (Vector Number_)+ deriving (Show)++instance QC.Arbitrary ArbVector2 where+ arbitrary = do+ dim <- QC.choose (0,maxDim)+ liftA2 ArbVector2+ (genVector maxElem dim) (genVector maxElem dim)++test :: DocTest.T ()+test = do+ DocTest.printPrefix "Numeric.BLAS.Vector:212: "+{-# LINE 212 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 212 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector xs) (ArbVector ys) (ArbVector zs) -> Vector.toList ((xs +++ ys) +++ zs) == Vector.toList (xs +++ (ys +++ zs)))+ DocTest.printPrefix "Numeric.BLAS.Vector:221: "+{-# LINE 221 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 221 "src/Numeric/BLAS/Vector.hs" #-}+ (QC.forAll (QC.choose (1,100)) $ \dim -> QC.forAll (QC.choose (0, dim-1)) $ \i -> QC.forAll (QC.choose (0, dim-1)) $ \j -> Vector.unit (Shape.ZeroBased dim) i == (Vector.swap i j (Vector.unit (Shape.ZeroBased dim) j) :: Vector Number_))+ DocTest.printPrefix "Numeric.BLAS.Vector:256: "+{-# LINE 256 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 256 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector2 xs ys) -> Vector.inner xs ys == Vector.dot (Vector.conjugate xs) ys)+ DocTest.printPrefix "Numeric.BLAS.Vector:309: "+{-# LINE 309 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 309 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector xs) -> Vector.sum xs == List.sum (Vector.toList xs))+ DocTest.printPrefix "Numeric.BLAS.Vector:369: "+{-# LINE 369 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 369 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector xs) -> Vector.normInf xs == List.maximum (0 : List.map absolute (Vector.toList xs)))+ DocTest.printPrefix "Numeric.BLAS.Vector:396: "+{-# LINE 396 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 396 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector xs) -> isNonEmpty xs ==> let (xi,xm) = Vector.argAbsMaximum xs in xs!xi == xm)+ DocTest.printPrefix "Numeric.BLAS.Vector:397: "+{-# LINE 397 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 397 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector xs) -> isNonEmpty xs ==> let (_xi,xm) = Vector.argAbsMaximum xs in List.all (\x -> absolute x <= absolute xm) $ Vector.toList xs)+ DocTest.printPrefix "Numeric.BLAS.Vector:398: "+{-# LINE 398 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 398 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector xs, ArbVector ys) -> isNonEmpty xs && isNonEmpty ys ==> let (_xi,xm) = Vector.argAbsMaximum xs; (_yi,ym) = Vector.argAbsMaximum ys; (zi,zm) = Vector.argAbsMaximum (xs+++ys) in case zi of Left _ -> xm==zm && absolute xm >= absolute ym; Right _ -> ym==zm && absolute xm < absolute ym)+ DocTest.printPrefix "Numeric.BLAS.Vector:393: "+{-# LINE 393 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.example+{-# LINE 393 "src/Numeric/BLAS/Vector.hs" #-}+ (Vector.argAbsMaximum $ Vector.autoFromList [1:+2, 3:+4, 5, 6 :: Complex_])+ [ExpectedLine [LineChunk "(3,6.0 :+ 0.0)"]]+ DocTest.printPrefix "Numeric.BLAS.Vector:445: "+{-# LINE 445 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 445 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbRealVector xs) -> isNonEmpty xs ==> Vector.argAbsMaximum xs == Vector.argAbs1Maximum xs)+ DocTest.printPrefix "Numeric.BLAS.Vector:442: "+{-# LINE 442 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.example+{-# LINE 442 "src/Numeric/BLAS/Vector.hs" #-}+ (Vector.argAbs1Maximum $ Vector.autoFromList [1:+2, 3:+4, 5, 6 :: Complex_])+ [ExpectedLine [LineChunk "(1,3.0 :+ 4.0)"]]+ DocTest.printPrefix "Numeric.BLAS.Vector:478: "+{-# LINE 478 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 478 "src/Numeric/BLAS/Vector.hs" #-}+ (QC.forAll (QC.choose (0,10)) $ \dim -> QC.forAll (genVector 3 dim) $ \xs -> approx 1e-2 (Vector.product xs) (List.product (Vector.toList (xs :: Vector Number_))))+ DocTest.printPrefix "Numeric.BLAS.Vector:488: "+{-# LINE 488 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 488 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbRealVector xs) -> isNonEmpty xs ==> Vector.minimum xs == List.minimum (Vector.toList xs))+ DocTest.printPrefix "Numeric.BLAS.Vector:489: "+{-# LINE 489 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 489 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbRealVector xs) -> isNonEmpty xs ==> Vector.maximum xs == List.maximum (Vector.toList xs))+ DocTest.printPrefix "Numeric.BLAS.Vector:490: "+{-# LINE 490 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 490 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbRealVector xs) -> isNonEmpty xs ==> - Vector.maximum xs == Vector.minimum (Vector.negate xs))+ DocTest.printPrefix "Numeric.BLAS.Vector:506: "+{-# LINE 506 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 506 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbRealVector xs) -> isNonEmpty xs ==> Vector.limits xs == Array.limits xs)+ DocTest.printPrefix "Numeric.BLAS.Vector:536: "+{-# LINE 536 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 536 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector xs) -> Vector.negate xs == Vector.scale minusOne xs)+ DocTest.printPrefix "Numeric.BLAS.Vector:537: "+{-# LINE 537 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 537 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector xs) -> Vector.scale 2 xs == xs |+| xs)+ DocTest.printPrefix "Numeric.BLAS.Vector:606: "+{-# LINE 606 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 606 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector2 xs ys) -> xs |+| ys == ys |+| xs)+ DocTest.printPrefix "Numeric.BLAS.Vector:607: "+{-# LINE 607 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 607 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector2 xs ys) -> xs == xs |-| ys |+| ys)+ DocTest.printPrefix "Numeric.BLAS.Vector:637: "+{-# LINE 637 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 637 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector xs) -> xs == Vector.negate (Vector.negate xs))+ DocTest.printPrefix "Numeric.BLAS.Vector:650: "+{-# LINE 650 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 650 "src/Numeric/BLAS/Vector.hs" #-}+ (QC.forAll (genNumber maxElem) $ \d (ArbVector xs) -> xs == Vector.raise (-d) (Vector.raise d xs))+ DocTest.printPrefix "Numeric.BLAS.Vector:667: "+{-# LINE 667 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 667 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector2 xs ys) -> Vector.mul xs ys == Vector.mul ys xs)+ DocTest.printPrefix "Numeric.BLAS.Vector:675: "+{-# LINE 675 "src/Numeric/BLAS/Vector.hs" #-}+ DocTest.property+{-# LINE 675 "src/Numeric/BLAS/Vector.hs" #-}+ (\(ArbVector2 xs ys) -> Vector.mulConj xs ys == Vector.mul (Vector.conjugate xs) ys)
+ test/Test/Float/Numeric/BLAS/Vector/Slice.hs view
@@ -0,0 +1,271 @@+-- Do not edit! Automatically created with doctest-extract from src/Numeric/BLAS/Vector/Slice.hs+{-# LINE 87 "src/Numeric/BLAS/Vector/Slice.hs" #-}++{-# OPTIONS_GHC -XTypeOperators #-}+module Test.Float.Numeric.BLAS.Vector.Slice where++import qualified Test.DocTest.Driver as DocTest++{-# LINE 89 "src/Numeric/BLAS/Vector/Slice.hs" #-}+import qualified Test.Slice as TestSlice+import Test.Float.Type (Number_)+import Test.Generator (genNumber)+import Test.Utility (approx, approxReal)++import qualified Numeric.BLAS.Matrix.RowMajor as Matrix+import qualified Numeric.BLAS.Vector.Slice as VectorSlice+import qualified Numeric.BLAS.Vector as Vector+import qualified Numeric.BLAS.Slice as Slice+import qualified Numeric.BLAS.Scalar as Scalar+import qualified Numeric.Netlib.Class as Class+import qualified Data.Array.Comfort.Storable as Array+import qualified Data.Array.Comfort.Shape as Shape+import qualified Data.List as List+import Numeric.BLAS.Vector ((+++), (|+|))+import Numeric.BLAS.Scalar (RealOf, absolute, minusOne)+import Data.Array.Comfort.Shape ((::+)((::+)))+import Data.Tuple.HT (mapPair)+import Data.Complex (Complex)+import Control.Monad (replicateM)+import Control.Applicative (liftA2)++import qualified Test.QuickCheck as QC+import Test.QuickCheck ((==>))++type Real_ = RealOf Number_+type Complex_ = Complex Real_++maxElem :: Integer+maxElem = 10++maxDim :: Int+maxDim = 100++maxDim1 :: Int+maxDim1 = 10++type ShapeInt = Shape.ZeroBased Int+type Shape = ShapeInt ::+ (ShapeInt, ShapeInt) ::+ ShapeInt+type Vector = Vector.Vector Shape+type Sliced = VectorSlice.T Shape ShapeInt++genDim :: QC.Gen Int+genDim = QC.choose (0,maxDim)++genShapeDim :: Int -> QC.Gen Shape+genShapeDim numRows = do+ left <- fmap Shape.ZeroBased $ QC.choose (0,maxDim)+ right <- fmap Shape.ZeroBased $ QC.choose (0,maxDim)+ columns <- fmap Shape.ZeroBased $ QC.choose (1,maxDim1)+ return (left ::+ (Shape.ZeroBased numRows, columns) ::+ right)++genShape :: QC.Gen Shape+genShape = genShapeDim =<< QC.choose (0,maxDim1)++genVector ::+ (Shape.C sh, Class.Floating a) =>+ sh -> QC.Gen a -> QC.Gen (Vector.Vector sh a)+genVector shape genElem =+ fmap (Vector.fromList shape) $+ replicateM (Shape.size shape) genElem++forAll_ :: (Show a) => QC.Gen a -> (a -> QC.Property) -> QC.Property+forAll_ = QC.forAll++real_ :: QC.Gen Real_+real_ = genNumber maxElem++complex_ :: QC.Gen Complex_+complex_ = genNumber maxElem++number_ :: QC.Gen Number_+number_ = genNumber maxElem++isNonEmpty :: Shape.C sh => VectorSlice.T shA sh a -> Bool+isNonEmpty xs = Shape.size (VectorSlice.shape xs) > 0++takeColumn ::+ (Shape.Indexed sh1, Shape.C sh, Shape.C sh0, Shape.C sh2) =>+ Shape.Index sh1 ->+ Vector.Vector (sh0 ::+ (sh, sh1) ::+ sh2) a ->+ VectorSlice.T (sh0 ::+ (sh, sh1) ::+ sh2) sh a+takeColumn c =+ VectorSlice.slice+ (Slice.column c . Slice.left . Slice.right) . VectorSlice.fromVector++listFromSlice ::+ (Shape.C sh, Class.Floating a) => VectorSlice.T shA sh a -> [a]+listFromSlice = Vector.toList . VectorSlice.toVector++genSlicedDim ::+ (Class.Floating a) =>+ Int -> QC.Gen a -> QC.Gen (Int, Vector a)+genSlicedDim numRows genElem = do+ shape@(_::+(_rows,columns)::+_) <- genShapeDim numRows+ c <- QC.elements (Shape.indices columns)+ fmap ((,) c) $ genVector shape genElem++genSliced ::+ (Class.Floating a) =>+ QC.Gen a -> QC.Gen (Int, Vector a)+genSliced genElem = flip genSlicedDim genElem =<< genDim++shrinkSliced ::+ (Shape.C sh0, Shape.Indexed sh1, QC.Arbitrary a, Class.Floating a) =>+ (Shape.Index sh1,+ Vector.Vector (ShapeInt ::+ (sh0, sh1) ::+ ShapeInt) a) ->+ [(Shape.Index sh1,+ Vector.Vector (ShapeInt ::+ (sh0, sh1) ::+ ShapeInt) a)]+shrinkSliced (c,xs) =+ let xs0 = Vector.takeLeft xs in+ let xs1 = Vector.takeRight xs in+ let xs10 = Vector.takeLeft xs1 in+ let xs11 = Vector.takeRight xs1 in+ map (\(ysl,ysr) ->+ (c,+ Vector.autoFromList ysl +++ xs10 +++ Vector.autoFromList ysr)) $+ QC.shrink (Vector.toList xs0, Vector.toList xs11)++forSliced ::+ (QC.Testable prop, QC.Arbitrary a, Class.Floating a, Show a) =>+ QC.Gen a -> (Sliced a -> prop) -> QC.Property+forSliced genElem prop =+ QC.forAllShrink (genSliced genElem) shrinkSliced+ (prop . uncurry takeColumn)++genSliced2 ::+ (Class.Floating a) =>+ QC.Gen a -> QC.Gen ((Int, Vector a), (Int, Vector a))+genSliced2 genElem = do+ dim <- genDim+ liftA2 (,) (genSlicedDim dim genElem) (genSlicedDim dim genElem)++forSliced2 ::+ (QC.Testable prop, Class.Floating a, Show a) =>+ QC.Gen a -> (Sliced a -> Sliced a -> prop) -> QC.Property+forSliced2 genElem prop =+ QC.forAll (genSliced2 genElem)+ (uncurry prop . mapPair (uncurry takeColumn, uncurry takeColumn))++test :: DocTest.T ()+test = do+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:300: "+{-# LINE 300 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 300 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (QC.forAll genShape $ \shape@(_::+(_rows,columns)::+_) -> QC.forAll (QC.elements (Shape.indices columns)) $ \c -> QC.forAll (genVector shape number_) $ \xs -> VectorSlice.extract (Slice.column c . Slice.left . Slice.right) xs == Matrix.takeColumn c (Vector.takeLeft (Vector.takeRight (xs :: Vector Number_))))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:302: "+{-# LINE 302 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 302 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forAll_ (TestSlice.genShapeSelect 4 100) $ \(TestSlice.ShapeSelect sh select) -> QC.forAll (genVector sh number_) $ \xs -> case TestSlice.instantiate sh select of TestSlice.Extraction slice cut -> VectorSlice.extract slice xs == cut xs)+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:334: "+{-# LINE 334 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 334 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced2 number_ $ \xs ys -> VectorSlice.inner xs ys == Vector.dot (VectorSlice.conjugate xs) (VectorSlice.toVector ys))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:335: "+{-# LINE 335 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 335 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced number_ $ \xs -> VectorSlice.inner xs xs == Scalar.fromReal (VectorSlice.norm2Squared xs))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:418: "+{-# LINE 418 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 418 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced number_ $ \xs -> VectorSlice.sum xs == List.sum (listFromSlice xs))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:491: "+{-# LINE 491 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 491 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced number_ $ \xs -> VectorSlice.normInf xs == List.maximum (0 : List.map absolute (listFromSlice xs)))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:502: "+{-# LINE 502 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 502 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced number_ $ \xs -> VectorSlice.normInf1 xs == List.maximum (0 : List.map Scalar.norm1 (listFromSlice xs)))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:517: "+{-# LINE 517 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 517 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced number_ $ \xs -> isNonEmpty xs ==> let (xi,xm) = VectorSlice.argAbsMaximum xs in VectorSlice.access xs xi == xm)+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:518: "+{-# LINE 518 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 518 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced number_ $ \xs -> isNonEmpty xs ==> let (_xi,xm) = VectorSlice.argAbsMaximum xs in List.all (\x -> absolute x <= absolute xm) $ listFromSlice xs)+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:519: "+{-# LINE 519 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 519 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced number_ $ \xs -> forSliced number_ $ \ys -> isNonEmpty xs && isNonEmpty ys ==> let (_xi,xm) = VectorSlice.argAbsMaximum xs; (_yi,ym) = VectorSlice.argAbsMaximum ys; (zi,zm) = Vector.argAbsMaximum (VectorSlice.toVector xs +++ VectorSlice.toVector ys) in case zi of Left _ -> xm==zm && absolute xm >= absolute ym; Right _ -> ym==zm && absolute xm < absolute ym)+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:566: "+{-# LINE 566 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 566 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced real_ $ \xs -> isNonEmpty xs ==> VectorSlice.argAbsMaximum xs == VectorSlice.argAbs1Maximum xs)+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:591: "+{-# LINE 591 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 591 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (QC.forAll genShape $ \sh@(_::+(_rows,columns)::+_) -> QC.forAll (QC.elements (Shape.indices columns)) $ \c -> QC.forAll (genVector sh $ genNumber 3) $ \xt -> let xs = takeColumn c xt in approx 1e-2 (VectorSlice.product xs) (List.product (listFromSlice (xs :: Sliced Number_))))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:602: "+{-# LINE 602 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 602 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced real_ $ \xs -> isNonEmpty xs ==> VectorSlice.minimum xs == List.minimum (listFromSlice xs))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:603: "+{-# LINE 603 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 603 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced real_ $ \xs -> isNonEmpty xs ==> VectorSlice.maximum xs == List.maximum (listFromSlice xs))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:619: "+{-# LINE 619 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 619 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced real_ $ \xs -> isNonEmpty xs ==> VectorSlice.limits xs == Array.limits (VectorSlice.toVector xs))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:649: "+{-# LINE 649 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 649 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced number_ $ \xs -> VectorSlice.negate xs == VectorSlice.scale minusOne xs)+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:650: "+{-# LINE 650 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 650 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced number_ $ \xs -> VectorSlice.scale 2 xs == VectorSlice.add xs xs)+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:709: "+{-# LINE 709 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 709 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced2 number_ $ \xs ys -> VectorSlice.add xs ys == VectorSlice.add ys xs)+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:710: "+{-# LINE 710 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 710 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced2 number_ $ \xs ys -> VectorSlice.toVector xs == VectorSlice.sub xs ys |+| VectorSlice.toVector ys)+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:735: "+{-# LINE 735 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 735 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced number_ $ \xs -> VectorSlice.toVector xs == Vector.negate (VectorSlice.negate xs))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:748: "+{-# LINE 748 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 748 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (QC.forAll (genNumber maxElem) $ \d -> forSliced number_ $ \xs -> VectorSlice.toVector xs == Vector.raise (-d) (VectorSlice.raise d xs))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:765: "+{-# LINE 765 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 765 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced2 number_ $ \xs ys -> VectorSlice.mul xs ys == VectorSlice.mul ys xs)+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:773: "+{-# LINE 773 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 773 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced2 number_ $ \xs ys -> VectorSlice.mulConj xs ys == Vector.mul (VectorSlice.conjugate xs) (VectorSlice.toVector ys))+ DocTest.printPrefix "Numeric.BLAS.Vector.Slice:897: "+{-# LINE 897 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ DocTest.property+{-# LINE 897 "src/Numeric/BLAS/Vector/Slice.hs" #-}+ (forSliced complex_ $ \xs -> approxReal 1e-2 (VectorSlice.norm2 xs) $ let (xrs,xis) = VectorSlice.unzipComplex xs in sqrt $ VectorSlice.norm2Squared xrs + VectorSlice.norm2Squared xis)
+ test/Test/Float/Type.hs view
@@ -0,0 +1,3 @@+module Test.Float.Type where++type Number_ = Float
+ test/Test/Generator.hs view
@@ -0,0 +1,28 @@+{-# LANGUAGE GADTs #-}+module Test.Generator where++import qualified Numeric.Netlib.Class as Class+import qualified Numeric.BLAS.Scalar as Scalar+import Data.Complex (Complex((:+)))++import Control.Applicative (liftA2)++import qualified Test.QuickCheck as QC+++genReal :: (Class.Real a) => Integer -> QC.Gen a+genReal maxElem = fmap fromIntegral $ QC.choose (-maxElem,maxElem)++genComplex :: (Class.Real a) => Integer -> QC.Gen (Complex a)+genComplex maxElem = liftA2 (:+) (genReal maxElem) (genReal maxElem)++genNumber :: (Class.Floating a) => Integer -> QC.Gen a+genNumber = genNumberAux Scalar.complexSingleton++genNumberAux ::+ (Class.Floating a) =>+ Scalar.ComplexSingleton a -> Integer -> QC.Gen a+genNumberAux singleton =+ case singleton of+ Scalar.Real -> genReal+ Scalar.Complex -> genComplex
+ test/Test/Numeric/BLAS/Slice.hs view
@@ -0,0 +1,53 @@+-- Do not edit! Automatically created with doctest-extract from src/Numeric/BLAS/Slice.hs+{-# LINE 13 "src/Numeric/BLAS/Slice.hs" #-}++module Test.Numeric.BLAS.Slice where++import qualified Test.DocTest.Driver as DocTest++{-# LINE 14 "src/Numeric/BLAS/Slice.hs" #-}+import qualified Numeric.BLAS.Slice as Slice+import Test.Slice (shapeInt)++import qualified Data.Array.Comfort.Boxed as Array+import qualified Data.Array.Comfort.Shape as Shape+import qualified Data.Map as Map+import Data.Array.Comfort.Shape ((::+)((::+)))+import Data.Array.Comfort.Boxed ((!))++import Control.Applicative (liftA3, pure)++import qualified Test.QuickCheck as QC++genSlice :: sh -> QC.Gen (Slice.T sh)+genSlice sh =+ liftA3 Slice.Cons (QC.choose (0,100)) (QC.choose (1,100)) (pure sh)++genSlice2 :: shA -> shB -> QC.Gen (Slice.T shA, Slice.T shB)+genSlice2 shA shB = do+ s <- QC.choose (0,100)+ k <- QC.choose (1,100)+ return (Slice.Cons s k shA, Slice.Cons s k shB)++test :: DocTest.T ()+test = do+ DocTest.printPrefix "Numeric.BLAS.Slice:58: "+{-# LINE 58 "src/Numeric/BLAS/Slice.hs" #-}+ DocTest.property+{-# LINE 58 "src/Numeric/BLAS/Slice.hs" #-}+ (QC.forAll (QC.choose (1,100)) $ \numRows -> QC.forAll (QC.choose (0,100)) $ \numColumns -> QC.forAll (genSlice (shapeInt numRows, shapeInt numColumns)) $ \slice -> QC.forAll (QC.elements $ Shape.indices $ shapeInt numRows) $ \ix -> Slice.row ix slice == Slice.rowArray slice ! ix)+ DocTest.printPrefix "Numeric.BLAS.Slice:70: "+{-# LINE 70 "src/Numeric/BLAS/Slice.hs" #-}+ DocTest.property+{-# LINE 70 "src/Numeric/BLAS/Slice.hs" #-}+ (QC.forAll (QC.choose (0,100)) $ \numRows -> QC.forAll (QC.choose (1,100)) $ \numColumns -> QC.forAll (genSlice (shapeInt numRows, shapeInt numColumns)) $ \slice -> QC.forAll (QC.elements $ Shape.indices $ shapeInt numColumns) $ \ix -> Slice.column ix slice == Slice.columnArray slice ! ix)+ DocTest.printPrefix "Numeric.BLAS.Slice:145: "+{-# LINE 145 "src/Numeric/BLAS/Slice.hs" #-}+ DocTest.property+{-# LINE 145 "src/Numeric/BLAS/Slice.hs" #-}+ (QC.forAll (fmap shapeInt $ QC.choose (0,100)) $ \shapeA -> QC.forAll (fmap shapeInt $ QC.choose (0,100)) $ \shapeB -> QC.forAll (fmap shapeInt $ QC.choose (0,100)) $ \shapeC -> QC.forAll (genSlice2 (Map.fromList $ ('a', shapeA) : ('b', shapeB) : ('c', shapeC) : []) (shapeA ::+ shapeB ::+ shapeC)) $ \(sliceMap, sliceParted) -> Slice.map sliceMap Map.! 'b' == Slice.left (Slice.right sliceParted))+ DocTest.printPrefix "Numeric.BLAS.Slice:147: "+{-# LINE 147 "src/Numeric/BLAS/Slice.hs" #-}+ DocTest.property+{-# LINE 147 "src/Numeric/BLAS/Slice.hs" #-}+ (QC.forAll (QC.choose (0,100)) $ \numRows -> QC.forAll (QC.choose (0,100)) $ \numColumns -> let rowShape = shapeInt numRows; columnShape = shapeInt numColumns; mapShape = Map.fromList $ map (\k -> (k, columnShape)) (Shape.indices rowShape) in QC.forAll (genSlice2 mapShape (rowShape, columnShape)) $ \(sliceMap, sliceMatrix) -> Map.toAscList (Slice.map sliceMap) == Array.toAssociations (Slice.rowArray sliceMatrix))
+ test/Test/Slice.hs view
@@ -0,0 +1,113 @@+{-# LANGUAGE ExistentialQuantification #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE StandaloneDeriving #-}+module Test.Slice where++import qualified Numeric.BLAS.Matrix.RowMajor as Matrix+import qualified Numeric.BLAS.Vector as Vector+import qualified Numeric.BLAS.Slice as Slice+import Numeric.BLAS.Vector (Vector)++import qualified Data.Array.Comfort.Shape as Shape+import Data.Array.Comfort.Shape ((::+)((::+)))++import qualified Numeric.Netlib.Class as Class++import qualified Test.QuickCheck as QC++import Prelude hiding (Either(Left, Right))+++data Extraction sh0 sh1 a =+ Extraction+ (Slice.T sh0 -> Slice.T sh1)+ (Vector sh0 a -> Vector sh1 a)+++type ShapeInt = Shape.ZeroBased Int++shapeInt :: Int -> ShapeInt+shapeInt = Shape.ZeroBased++data Select sh where+ All :: Select ShapeInt+ Left :: (Shape.Indexed shL, Shape.C shR) => Select shL -> Select (shL::+shR)+ Right :: (Shape.Indexed shR, Shape.C shL) => Select shR -> Select (shL::+shR)+ Row ::+ (Shape.Indexed shR, Shape.Index shR ~ ixR, Show ixR, Shape.C shC) =>+ ixR -> Select shC -> Select (shR,shC)+ Column ::+ (Shape.Indexed shC, Shape.Index shC ~ ixC, Show ixC, Shape.C shR) =>+ ixC -> Select shR -> Select (shR,shC)++data ShapeSelect =+ forall sh ix.+ (Shape.Indexed sh, Show sh, Shape.Index sh ~ ix, Show ix) =>+ ShapeSelect sh (Select sh)++deriving instance Show (Select sh)+deriving instance Show ShapeSelect+++instantiate :: (Class.Floating a) => sh -> Select sh -> Extraction sh ShapeInt a+instantiate sh select =+ case select of+ All -> Extraction id id+ Left sel ->+ case instantiate (case sh of shL::+_shR -> shL) sel of+ Extraction fs fv ->+ Extraction (fs . Slice.left) (fv . Vector.takeLeft)+ Right sel ->+ case instantiate (case sh of _shL::+shR -> shR) sel of+ Extraction fs fv ->+ Extraction (fs . Slice.right) (fv . Vector.takeRight)+ Row ix sel ->+ case instantiate (snd sh) sel of+ Extraction fs fv ->+ Extraction (fs . Slice.row ix) (fv . Matrix.takeRow ix)+ Column ix sel ->+ case instantiate (fst sh) sel of+ Extraction fs fv ->+ Extraction (fs . Slice.column ix) (fv . Matrix.takeColumn ix)+++genPrimShape :: Int -> QC.Gen ShapeInt+genPrimShape maxSize =+ fmap Shape.ZeroBased $ QC.choose (1, maxSize)++genShapeSelect :: Int -> Int -> QC.Gen ShapeSelect+genShapeSelect maxDepth maxSize =+ if maxDepth <= 1 || maxSize < 2+ then fmap (flip ShapeSelect All) $ genPrimShape maxSize+ else+ QC.oneof $+ fmap (flip ShapeSelect All) (genPrimShape maxSize) :+ (do+ pivot <- QC.choose (1, maxSize)+ left <- genShapeSelect (maxDepth-1) pivot+ right <- genShapeSelect (maxDepth-1) (maxSize-pivot)+ case (left, right) of+ (ShapeSelect shl sell, ShapeSelect shr selr) ->+ fmap (ShapeSelect (shl::+shr)) $+ QC.elements [Left sell, Right selr]) :+ (do+ let sizeReal :: Double+ sizeReal = fromIntegral maxSize+ pivot <- QC.choose (0, logBase 2 sizeReal)+ let maxRows = 2**pivot+ let maxColumns = sizeReal / maxRows+ rows <- genShapeSelect (maxDepth-1) (floor maxRows)+ columns <- genShapeSelect (maxDepth-1) (floor maxColumns)+ case (rows, columns) of+ (ShapeSelect shr selr, ShapeSelect shc selc) ->+ fmap (ShapeSelect (shr,shc)) $+ QC.oneof $+ (do+ ix <- QC.elements (Shape.indices shr)+ return $ Row ix selc) :+ (do+ ix <- QC.elements (Shape.indices shc)+ return $ Column ix selr) :+ []) :+ []
+ test/Test/Utility.hs view
@@ -0,0 +1,19 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}+module Test.Utility where++import qualified Numeric.Netlib.Class as Class+import Numeric.BLAS.Scalar (RealOf, absolute)++import qualified Test.QuickCheck as QC+++approx ::+ (Class.Floating a, RealOf a ~ ar, Class.Real ar, Show a) =>+ ar -> a -> a -> QC.Property+approx tol x y =+ QC.counterexample (show (x,y)) $ absolute (x-y) <= tol++approxReal :: (Class.Real a, Show a) => a -> a -> a -> QC.Property+approxReal tol x y =+ QC.counterexample (show (x,y)) $ abs (x-y) <= tol