diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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. diff --git a/Makefile b/Makefile new file mode 100644 --- /dev/null +++ b/Makefile @@ -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:$*:" $@ diff --git a/Setup.lhs b/Setup.lhs new file mode 100644 --- /dev/null +++ b/Setup.lhs @@ -0,0 +1,3 @@ +#! /usr/bin/env runhaskell +> import Distribution.Simple +> main = defaultMain diff --git a/comfort-blas.cabal b/comfort-blas.cabal new file mode 100644 --- /dev/null +++ b/comfort-blas.cabal @@ -0,0 +1,87 @@ +Cabal-Version: 2.2 +Name: comfort-blas +Version: 0.0 +License: BSD-3-Clause +License-File: LICENSE +Author: Henning Thielemann +Maintainer: Henning Thielemann +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 diff --git a/src/Numeric/BLAS/Matrix/Layout.hs b/src/Numeric/BLAS/Matrix/Layout.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/BLAS/Matrix/Layout.hs @@ -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' diff --git a/src/Numeric/BLAS/Matrix/Modifier.hs b/src/Numeric/BLAS/Matrix/Modifier.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/BLAS/Matrix/Modifier.hs @@ -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 = (<>) diff --git a/src/Numeric/BLAS/Matrix/RowMajor.hs b/src/Numeric/BLAS/Matrix/RowMajor.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/BLAS/Matrix/RowMajor.hs @@ -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) diff --git a/src/Numeric/BLAS/Private.hs b/src/Numeric/BLAS/Private.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/BLAS/Private.hs @@ -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) diff --git a/src/Numeric/BLAS/Scalar.hs b/src/Numeric/BLAS/Scalar.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/BLAS/Scalar.hs @@ -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 diff --git a/src/Numeric/BLAS/Slice.hs b/src/Numeric/BLAS/Slice.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/BLAS/Slice.hs @@ -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 diff --git a/src/Numeric/BLAS/Vector.hs b/src/Numeric/BLAS/Vector.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/BLAS/Vector.hs @@ -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) diff --git a/src/Numeric/BLAS/Vector/Slice.hs b/src/Numeric/BLAS/Vector/Slice.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/BLAS/Vector/Slice.hs @@ -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) diff --git a/test/Main.hs b/test/Main.hs new file mode 100644 --- /dev/null +++ b/test/Main.hs @@ -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 diff --git a/test/Test/ComplexDouble/Numeric/BLAS/Vector.hs b/test/Test/ComplexDouble/Numeric/BLAS/Vector.hs new file mode 100644 --- /dev/null +++ b/test/Test/ComplexDouble/Numeric/BLAS/Vector.hs @@ -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) diff --git a/test/Test/ComplexDouble/Numeric/BLAS/Vector/Slice.hs b/test/Test/ComplexDouble/Numeric/BLAS/Vector/Slice.hs new file mode 100644 --- /dev/null +++ b/test/Test/ComplexDouble/Numeric/BLAS/Vector/Slice.hs @@ -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) diff --git a/test/Test/ComplexDouble/Type.hs b/test/Test/ComplexDouble/Type.hs new file mode 100644 --- /dev/null +++ b/test/Test/ComplexDouble/Type.hs @@ -0,0 +1,5 @@ +module Test.ComplexDouble.Type where + +import Data.Complex (Complex) + +type Number_ = Complex Double diff --git a/test/Test/ComplexFloat/Numeric/BLAS/Vector.hs b/test/Test/ComplexFloat/Numeric/BLAS/Vector.hs new file mode 100644 --- /dev/null +++ b/test/Test/ComplexFloat/Numeric/BLAS/Vector.hs @@ -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) diff --git a/test/Test/ComplexFloat/Numeric/BLAS/Vector/Slice.hs b/test/Test/ComplexFloat/Numeric/BLAS/Vector/Slice.hs new file mode 100644 --- /dev/null +++ b/test/Test/ComplexFloat/Numeric/BLAS/Vector/Slice.hs @@ -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) diff --git a/test/Test/ComplexFloat/Type.hs b/test/Test/ComplexFloat/Type.hs new file mode 100644 --- /dev/null +++ b/test/Test/ComplexFloat/Type.hs @@ -0,0 +1,5 @@ +module Test.ComplexFloat.Type where + +import Data.Complex (Complex) + +type Number_ = Complex Float diff --git a/test/Test/Double/Numeric/BLAS/Vector.hs b/test/Test/Double/Numeric/BLAS/Vector.hs new file mode 100644 --- /dev/null +++ b/test/Test/Double/Numeric/BLAS/Vector.hs @@ -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) diff --git a/test/Test/Double/Numeric/BLAS/Vector/Slice.hs b/test/Test/Double/Numeric/BLAS/Vector/Slice.hs new file mode 100644 --- /dev/null +++ b/test/Test/Double/Numeric/BLAS/Vector/Slice.hs @@ -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) diff --git a/test/Test/Double/Type.hs b/test/Test/Double/Type.hs new file mode 100644 --- /dev/null +++ b/test/Test/Double/Type.hs @@ -0,0 +1,3 @@ +module Test.Double.Type where + +type Number_ = Double diff --git a/test/Test/Float/Numeric/BLAS/Vector.hs b/test/Test/Float/Numeric/BLAS/Vector.hs new file mode 100644 --- /dev/null +++ b/test/Test/Float/Numeric/BLAS/Vector.hs @@ -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) diff --git a/test/Test/Float/Numeric/BLAS/Vector/Slice.hs b/test/Test/Float/Numeric/BLAS/Vector/Slice.hs new file mode 100644 --- /dev/null +++ b/test/Test/Float/Numeric/BLAS/Vector/Slice.hs @@ -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) diff --git a/test/Test/Float/Type.hs b/test/Test/Float/Type.hs new file mode 100644 --- /dev/null +++ b/test/Test/Float/Type.hs @@ -0,0 +1,3 @@ +module Test.Float.Type where + +type Number_ = Float diff --git a/test/Test/Generator.hs b/test/Test/Generator.hs new file mode 100644 --- /dev/null +++ b/test/Test/Generator.hs @@ -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 diff --git a/test/Test/Numeric/BLAS/Slice.hs b/test/Test/Numeric/BLAS/Slice.hs new file mode 100644 --- /dev/null +++ b/test/Test/Numeric/BLAS/Slice.hs @@ -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)) diff --git a/test/Test/Slice.hs b/test/Test/Slice.hs new file mode 100644 --- /dev/null +++ b/test/Test/Slice.hs @@ -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) : + []) : + [] diff --git a/test/Test/Utility.hs b/test/Test/Utility.hs new file mode 100644 --- /dev/null +++ b/test/Test/Utility.hs @@ -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