{-# 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
>>> :set -XGADTs
>>> import qualified Test.Slice as TestSlice
>>> import Test.NumberModule.Numeric.BLAS.Vector
>>>    (maxElem, maxDim, genVector, number_, real_, complex_)
>>> 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.Applicative (liftA2)
>>>
>>> import qualified Test.QuickCheck as QC
>>> import Test.QuickCheck ((==>))
>>>
>>> type Real_ = RealOf Number_
>>> type Complex_ = Complex Real_
>>>
>>> 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)
>>>
>>> forAll_ :: (Show a) => QC.Gen a -> (a -> QC.Property) -> QC.Property
>>> forAll_ = QC.forAll
>>>
>>> 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 :: forall sh slice a. T sh slice a -> slice
shape (Cons (Slice.Cons Int
_s Int
_k slice
slc) Vector sh a
_arr) = slice
slc

mapShape :: (slice0 -> slice1) -> T sh slice0 a -> T sh slice1 a
mapShape :: forall slice0 slice1 sh a.
(slice0 -> slice1) -> T sh slice0 a -> T sh slice1 a
mapShape slice0 -> slice1
f (Cons (Slice.Cons Int
s Int
k slice0
slc) Vector sh a
arr) =
   forall sh slice a. T slice -> Vector sh a -> T sh slice a
Cons (forall sh. Int -> Int -> sh -> T sh
Slice.Cons Int
s Int
k (slice0 -> slice1
f slice0
slc)) Vector sh a
arr


increment :: T sh slice a -> Int
increment :: forall sh slice a. T sh slice a -> Int
increment (Cons (Slice.Cons Int
_s Int
k slice
_slc) Vector sh a
_arr) = Int
k


startArg ::
   (Storable a) =>
   T sh slice a -> Call.FortranIO r (Ptr a)
startArg :: forall a sh slice r.
Storable a =>
T sh slice a -> FortranIO r (Ptr a)
startArg (Cons (Slice.Cons Int
s Int
_k slice
_slice) (Array sh
_sh ForeignPtr a
x)) = do
   Ptr a
sxPtr <- forall {k} (r :: k) (m :: k -> *) a.
((a -> m r) -> m r) -> ContT r m a
ContT forall a b. (a -> b) -> a -> b
$ forall a b. ForeignPtr a -> (Ptr a -> IO b) -> IO b
withForeignPtr ForeignPtr a
x
   forall (m :: * -> *) a. Monad m => a -> m a
return (forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
sxPtr Int
s)

sliceArg ::
   (Storable a) =>
   T sh slice a -> Call.FortranIO r (Ptr a, Ptr CInt)
sliceArg :: forall a sh slice r.
Storable a =>
T sh slice a -> FortranIO r (Ptr a, Ptr CInt)
sliceArg T sh slice a
x =
   forall (f :: * -> *) a b c.
Applicative f =>
(a -> b -> c) -> f a -> f b -> f c
liftA2 (,) (forall a sh slice r.
Storable a =>
T sh slice a -> FortranIO r (Ptr a)
startArg T sh slice a
x) (forall r. Int -> FortranIO r (Ptr CInt)
Call.cint forall a b. (a -> b) -> a -> b
$ forall sh slice a. T sh slice a -> Int
increment T sh slice a
x)

sizeSliceArg ::
   (Shape.C sh, Storable a) =>
   T shA sh a -> ContT r IO (Ptr CInt, Ptr a, Ptr CInt)
sizeSliceArg :: forall sh a shA r.
(C sh, Storable a) =>
T shA sh a -> ContT r IO (Ptr CInt, Ptr a, Ptr CInt)
sizeSliceArg T shA sh a
x =
   forall (f :: * -> *) a b c.
Applicative f =>
(a -> b -> c) -> f a -> f b -> f c
liftA2
      (\Ptr CInt
nPtr (Ptr a
xPtr,Ptr CInt
incxPtr) -> (Ptr CInt
nPtr, Ptr a
xPtr,Ptr CInt
incxPtr))
      (forall r. Int -> FortranIO r (Ptr CInt)
Call.cint forall a b. (a -> b) -> a -> b
$ forall sh. C sh => sh -> Int
Shape.size forall a b. (a -> b) -> a -> b
$ forall sh slice a. T sh slice a -> slice
shape T shA sh a
x)
      (forall a sh slice r.
Storable a =>
T sh slice a -> FortranIO r (Ptr a, Ptr CInt)
sliceArg T shA sh a
x)

infixl 4 <*|>

(<*|>) ::
   (Storable a) =>
   Call.FortranIO r (Ptr a -> Ptr CInt -> b) ->
   T sh slice a ->
   Call.FortranIO r b
FortranIO r (Ptr a -> Ptr CInt -> b)
f <*|> :: forall a r b sh slice.
Storable a =>
FortranIO r (Ptr a -> Ptr CInt -> b)
-> T sh slice a -> FortranIO r b
<*|> T sh slice a
x  =  forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry FortranIO r (Ptr a -> Ptr CInt -> b)
f forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*> forall a sh slice r.
Storable a =>
T sh slice a -> FortranIO r (Ptr a, Ptr CInt)
sliceArg T sh slice a
x


data T sh slice a = Cons (Slice.T slice) (Vector sh a)
   deriving (Int -> T sh slice a -> ShowS
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
forall sh slice a.
(C sh, Storable a, Show slice, Show sh, Show a) =>
Int -> T sh slice a -> ShowS
forall sh slice a.
(C sh, Storable a, Show slice, Show sh, Show a) =>
[T sh slice a] -> ShowS
forall sh slice a.
(C sh, Storable a, Show slice, Show sh, Show a) =>
T sh slice a -> String
showList :: [T sh slice a] -> ShowS
$cshowList :: forall sh slice a.
(C sh, Storable a, Show slice, Show sh, Show a) =>
[T sh slice a] -> ShowS
show :: T sh slice a -> String
$cshow :: forall sh slice a.
(C sh, Storable a, Show slice, Show sh, Show a) =>
T sh slice a -> String
showsPrec :: Int -> T sh slice a -> ShowS
$cshowsPrec :: forall sh slice a.
(C sh, Storable a, Show slice, Show sh, Show a) =>
Int -> T sh slice a -> ShowS
Show)

toVector ::
   (Shape.C slice, Class.Floating a) =>
   T sh slice a -> Vector slice a
toVector :: forall slice a sh.
(C slice, Floating a) =>
T sh slice a -> Vector slice a
toVector T sh slice a
x =
   forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize (forall sh slice a. T sh slice a -> slice
shape T sh slice a
x) forall a b. (a -> b) -> a -> b
$ \Int
n Ptr a
syPtr ->
   forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT forall a b. (a -> b) -> a -> b
$ forall r a. FortranIO r (IO a) -> FortranIO r a
Call.run forall a b. (a -> b) -> a -> b
$
      forall (f :: * -> *) a. Applicative f => a -> f a
pure forall a.
Floating a =>
Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
Blas.copy
         forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*> forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
         forall a r b sh slice.
Storable a =>
FortranIO r (Ptr a -> Ptr CInt -> b)
-> T sh slice a -> FortranIO r b
<*|> T sh slice a
x
         forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*> forall (f :: * -> *) a. Applicative f => a -> f a
pure Ptr a
syPtr
         forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*> forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1

fromVector :: (Shape.C sh) => Vector sh a -> T sh sh a
fromVector :: forall sh a. C sh => Vector sh a -> T sh sh a
fromVector Vector sh a
xs = forall sh slice a. T slice -> Vector sh a -> T sh slice a
Cons (forall sh. C sh => sh -> T sh
Slice.fromShape forall a b. (a -> b) -> a -> b
$ forall sh a. Array sh a -> sh
Array.shape Vector sh a
xs) Vector sh a
xs

slice :: (Slice.T shA -> Slice.T shB) -> T sh shA a -> T sh shB a
slice :: forall shA shB sh a. (T shA -> T shB) -> T sh shA a -> T sh shB a
slice T shA -> T shB
f (Cons T shA
slc Vector sh a
xs) = forall sh slice a. T slice -> Vector sh a -> T sh slice a
Cons (T shA -> T shB
f T shA
slc) Vector sh a
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 :: forall slice sh a.
(C slice, C sh, Floating a) =>
(T sh -> T slice) -> Vector sh a -> Vector slice a
extract T sh -> T slice
slc Vector sh a
xs = forall slice a sh.
(C slice, Floating a) =>
T sh slice a -> Vector slice a
toVector forall a b. (a -> b) -> a -> b
$ forall shA shB sh a. (T shA -> T shB) -> T sh shA a -> T sh shB a
slice T sh -> T slice
slc forall a b. (a -> b) -> a -> b
$ forall sh a. C sh => Vector sh a -> T sh sh a
fromVector Vector sh a
xs


access, (!) ::
   (Shape.C shA, Shape.Indexed sh, Storable a) =>
   T shA sh a -> Shape.Index sh -> a
access :: forall shA sh a.
(C shA, Indexed sh, Storable a) =>
T shA sh a -> Index sh -> a
access (Cons (Slice.Cons Int
s Int
k sh
ssh) (Array shA
sh ForeignPtr a
x)) Index sh
ix =
   forall sh a. sh -> ForeignPtr a -> Array sh a
Array (forall sh. sh -> Deferred sh
Shape.Deferred shA
sh) ForeignPtr a
x
   forall sh a.
(Indexed sh, Storable a) =>
Array sh a -> Index sh -> a
Array.!
   forall sh. Int -> DeferredIndex sh
Shape.DeferredIndex (Int
s forall a. Num a => a -> a -> a
+ Int
k forall a. Num a => a -> a -> a
* forall sh. Indexed sh => sh -> Index sh -> Int
Shape.offset sh
ssh Index sh
ix)
! :: forall shA sh a.
(C shA, Indexed sh, Storable a) =>
T shA sh a -> Index sh -> a
(!) = forall shA sh a.
(C shA, Indexed sh, Storable a) =>
T shA sh a -> Index sh -> a
access



dot ::
   (Shape.C sh, Eq sh, Class.Floating a) =>
   T shA sh a -> T shB sh a -> a
dot :: forall sh a shA shB.
(C sh, Eq sh, Floating a) =>
T shA sh a -> T shB sh a -> a
dot =
   forall (f :: * -> *) (g :: * -> *) a. Dot f g a -> f a -> g a -> a
runDot forall a b. (a -> b) -> a -> b
$
   forall a (f :: * -> *).
Floating a =>
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
Class.switchFloating
      (forall (f :: * -> *) (g :: * -> *) a.
(f a -> g a -> a) -> Dot f g a
Dot forall sh a shA shB.
(C sh, Eq sh, Real a) =>
T shA sh a -> T shB sh a -> a
dotReal)
      (forall (f :: * -> *) (g :: * -> *) a.
(f a -> g a -> a) -> Dot f g a
Dot forall sh a shA shB.
(C sh, Eq sh, Real a) =>
T shA sh a -> T shB sh a -> a
dotReal)
      (forall (f :: * -> *) (g :: * -> *) a.
(f a -> g a -> a) -> Dot f g a
Dot forall sh a shA shB.
(C sh, Eq sh, Real a) =>
T shA sh (Complex a) -> T shB sh (Complex a) -> Complex a
dotComplex)
      (forall (f :: * -> *) (g :: * -> *) a.
(f a -> g a -> a) -> Dot f g a
Dot forall sh a shA shB.
(C sh, Eq sh, Real a) =>
T shA sh (Complex a) -> T shB sh (Complex a) -> Complex a
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 :: forall sh a shA shB.
(C sh, Eq sh, Floating a) =>
T shA sh a -> T shB sh a -> a
inner =
   forall (f :: * -> *) (g :: * -> *) a. Dot f g a -> f a -> g a -> a
runDot forall a b. (a -> b) -> a -> b
$
   forall a (f :: * -> *).
Floating a =>
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
Class.switchFloating
      (forall (f :: * -> *) (g :: * -> *) a.
(f a -> g a -> a) -> Dot f g a
Dot forall sh a shA shB.
(C sh, Eq sh, Real a) =>
T shA sh a -> T shB sh a -> a
dotReal)
      (forall (f :: * -> *) (g :: * -> *) a.
(f a -> g a -> a) -> Dot f g a
Dot forall sh a shA shB.
(C sh, Eq sh, Real a) =>
T shA sh a -> T shB sh a -> a
dotReal)
      (forall (f :: * -> *) (g :: * -> *) a.
(f a -> g a -> a) -> Dot f g a
Dot forall a b. (a -> b) -> a -> b
$ forall sh a shB.
(C sh, Eq sh, Real a) =>
Vector sh (Complex a) -> T shB sh (Complex a) -> Complex a
innerComplex forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall slice a sh.
(C slice, Floating a) =>
T sh slice a -> Vector slice a
toVector)
      (forall (f :: * -> *) (g :: * -> *) a.
(f a -> g a -> a) -> Dot f g a
Dot forall a b. (a -> b) -> a -> b
$ forall sh a shB.
(C sh, Eq sh, Real a) =>
Vector sh (Complex a) -> T shB sh (Complex a) -> Complex a
innerComplex forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall slice a sh.
(C slice, Floating a) =>
T sh slice a -> Vector slice a
toVector)


newtype Dot f g a = Dot {forall (f :: * -> *) (g :: * -> *) a. Dot f g a -> f a -> g a -> a
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 :: forall sh a shA shB.
(C sh, Eq sh, Real a) =>
T shA sh a -> T shB sh a -> a
dotReal T shA sh a
x T shB sh a
y = forall a. IO a -> a
unsafePerformIO forall a b. (a -> b) -> a -> b
$ do
   let shX :: sh
shX = forall sh slice a. T sh slice a -> slice
shape T shA sh a
x
   let shY :: sh
shY = forall sh slice a. T sh slice a -> slice
shape T shB sh a
y
   String -> Bool -> IO ()
Call.assert String
"dot: shapes mismatch" (sh
shX forall a. Eq a => a -> a -> Bool
== sh
shY)
   forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT forall a b. (a -> b) -> a -> b
$ do
      Ptr CInt
nPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint forall a b. (a -> b) -> a -> b
$ forall sh. C sh => sh -> Int
Shape.size sh
shX
      (Ptr a
sxPtr, Ptr CInt
incxPtr) <- forall a sh slice r.
Storable a =>
T sh slice a -> FortranIO r (Ptr a, Ptr CInt)
sliceArg T shA sh a
x
      (Ptr a
syPtr, Ptr CInt
incyPtr) <- forall a sh slice r.
Storable a =>
T sh slice a -> FortranIO r (Ptr a, Ptr CInt)
sliceArg T shB sh a
y
      forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a.
Real a =>
Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO a
BlasReal.dot Ptr CInt
nPtr Ptr a
sxPtr Ptr CInt
incxPtr Ptr a
syPtr Ptr CInt
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 :: forall sh a shA shB.
(C sh, Eq sh, Real a) =>
T shA sh (Complex a) -> T shB sh (Complex a) -> Complex a
dotComplex T shA sh (Complex a)
x T shB sh (Complex a)
y = forall a. IO a -> a
unsafePerformIO forall a b. (a -> b) -> a -> b
$ do
   String -> Bool -> IO ()
Call.assert String
"dot: shapes mismatch" (forall sh slice a. T sh slice a -> slice
shape T shA sh (Complex a)
x forall a. Eq a => a -> a -> Bool
== forall sh slice a. T sh slice a -> slice
shape T shB sh (Complex a)
y)
   forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT forall a b. (a -> b) -> a -> b
$ do
      Ptr CChar
transPtr <- forall r. Char -> FortranIO r (Ptr CChar)
Call.char Char
'N'
      Ptr CInt
mPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1
      Ptr CInt
nPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint forall a b. (a -> b) -> a -> b
$ forall sh. C sh => sh -> Int
Shape.size forall a b. (a -> b) -> a -> b
$ forall sh slice a. T sh slice a -> slice
shape T shA sh (Complex a)
x
      Ptr (Complex a)
alphaPtr <- forall a r. Floating a => a -> FortranIO r (Ptr a)
Call.number forall a. Floating a => a
Scalar.one
      (Ptr (Complex a)
xPtr, Ptr CInt
ldxPtr) <- forall a sh slice r.
Storable a =>
T sh slice a -> FortranIO r (Ptr a, Ptr CInt)
sliceArg T shA sh (Complex a)
x
      (Ptr (Complex a)
yPtr, Ptr CInt
incyPtr) <- forall a sh slice r.
Storable a =>
T sh slice a -> FortranIO r (Ptr a, Ptr CInt)
sliceArg T shB sh (Complex a)
y
      Ptr (Complex a)
betaPtr <- forall a r. Floating a => a -> FortranIO r (Ptr a)
Call.number forall a. Floating a => a
Scalar.zero
      Ptr (Complex a)
zPtr <- forall a r. Storable a => FortranIO r (Ptr a)
Call.alloca
      Ptr CInt
inczPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1
      forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$
         forall a.
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 ()
Private.gemv
            Ptr CChar
transPtr Ptr CInt
mPtr Ptr CInt
nPtr Ptr (Complex a)
alphaPtr Ptr (Complex a)
xPtr Ptr CInt
ldxPtr
            Ptr (Complex a)
yPtr Ptr CInt
incyPtr Ptr (Complex a)
betaPtr Ptr (Complex a)
zPtr Ptr CInt
inczPtr
      forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a. Storable a => Ptr a -> IO a
peek Ptr (Complex a)
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 :: forall sh a shB.
(C sh, Eq sh, Real a) =>
Vector sh (Complex a) -> T shB sh (Complex a) -> Complex a
innerComplex (Array sh
shX ForeignPtr (Complex a)
x) T shB sh (Complex a)
y = forall a. IO a -> a
unsafePerformIO forall a b. (a -> b) -> a -> b
$ do
   String -> Bool -> IO ()
Call.assert String
"dot: shapes mismatch" (sh
shX forall a. Eq a => a -> a -> Bool
== forall sh slice a. T sh slice a -> slice
shape T shB sh (Complex a)
y)
   forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT forall a b. (a -> b) -> a -> b
$ do
      let m :: Int
m = forall sh. C sh => sh -> Int
Shape.size sh
shX
      Ptr CChar
transPtr <- forall r. Char -> FortranIO r (Ptr CChar)
Call.char Char
'C'
      Ptr CInt
mPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
m
      Ptr CInt
nPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1
      Ptr (Complex a)
alphaPtr <- forall a r. Floating a => a -> FortranIO r (Ptr a)
Call.number forall a. Floating a => a
Scalar.one
      Ptr (Complex a)
xPtr <- forall {k} (r :: k) (m :: k -> *) a.
((a -> m r) -> m r) -> ContT r m a
ContT forall a b. (a -> b) -> a -> b
$ forall a b. ForeignPtr a -> (Ptr a -> IO b) -> IO b
withForeignPtr ForeignPtr (Complex a)
x
      Ptr CInt
ldxPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
m
      (Ptr (Complex a)
yPtr, Ptr CInt
incyPtr) <- forall a sh slice r.
Storable a =>
T sh slice a -> FortranIO r (Ptr a, Ptr CInt)
sliceArg T shB sh (Complex a)
y
      Ptr (Complex a)
betaPtr <- forall a r. Floating a => a -> FortranIO r (Ptr a)
Call.number forall a. Floating a => a
Scalar.zero
      Ptr (Complex a)
zPtr <- forall a r. Storable a => FortranIO r (Ptr a)
Call.alloca
      Ptr CInt
inczPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1
      forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$
         forall a.
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 ()
Private.gemv
            Ptr CChar
transPtr Ptr CInt
mPtr Ptr CInt
nPtr Ptr (Complex a)
alphaPtr Ptr (Complex a)
xPtr Ptr CInt
ldxPtr
            Ptr (Complex a)
yPtr Ptr CInt
incyPtr Ptr (Complex a)
betaPtr Ptr (Complex a)
zPtr Ptr CInt
inczPtr
      forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a. Storable a => Ptr a -> IO a
peek Ptr (Complex a)
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 :: forall sh a shA. (C sh, Floating a) => T shA sh a -> a
sum T shA sh a
x = forall a. IO a -> a
unsafePerformIO forall a b. (a -> b) -> a -> b
$ forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT forall a b. (a -> b) -> a -> b
$ do
   Ptr a
xPtr <- forall a sh slice r.
Storable a =>
T sh slice a -> FortranIO r (Ptr a)
startArg T shA sh a
x
   forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a. Floating a => Int -> Ptr a -> Int -> IO a
Private.sum (forall sh. C sh => sh -> Int
Shape.size forall a b. (a -> b) -> a -> b
$ forall sh slice a. T sh slice a -> slice
shape T shA sh a
x) Ptr a
xPtr (forall sh slice a. T sh slice a -> Int
increment T shA sh a
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 :: forall sh a shA. (C sh, Floating a) => T shA sh a -> RealOf a
absSum T shA sh a
arr = forall a. IO a -> a
unsafePerformIO forall a b. (a -> b) -> a -> b
$
   forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT forall a b. (a -> b) -> a -> b
$ forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b c d. (a -> b -> c -> d) -> (a, b, c) -> d
uncurry3 forall a.
Floating a =>
Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a)
asum forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
=<< forall sh a shA r.
(C sh, Storable a) =>
T shA sh a -> ContT r IO (Ptr CInt, Ptr a, Ptr CInt)
sizeSliceArg T shA sh a
arr

asum :: Class.Floating a => Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a)
asum :: forall a.
Floating a =>
Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a)
asum =
   forall a. Nrm a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a)
getNrm forall a b. (a -> b) -> a -> b
$
   forall a (f :: * -> *).
Floating a =>
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
Class.switchFloating
      (forall a. (Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a)) -> Nrm a
Nrm forall a. Real a => Ptr CInt -> Ptr a -> Ptr CInt -> IO a
BlasReal.asum) (forall a. (Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a)) -> Nrm a
Nrm forall a. Real a => Ptr CInt -> Ptr a -> Ptr CInt -> IO a
BlasReal.asum)
      (forall a. (Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a)) -> Nrm a
Nrm forall a. Real a => Ptr CInt -> Ptr (Complex a) -> Ptr CInt -> IO a
BlasComplex.casum) (forall a. (Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a)) -> Nrm a
Nrm forall a. Real a => Ptr CInt -> Ptr (Complex a) -> Ptr CInt -> IO a
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 :: forall sh a shA. (C sh, Floating a) => T shA sh a -> RealOf a
norm2 T shA sh a
arr = forall a. IO a -> a
unsafePerformIO forall a b. (a -> b) -> a -> b
$
   forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT forall a b. (a -> b) -> a -> b
$ forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b c d. (a -> b -> c -> d) -> (a, b, c) -> d
uncurry3 forall a.
Floating a =>
Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a)
nrm2 forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
=<< forall sh a shA r.
(C sh, Storable a) =>
T shA sh a -> ContT r IO (Ptr CInt, Ptr a, Ptr CInt)
sizeSliceArg T shA sh a
arr

nrm2 :: Class.Floating a => Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a)
nrm2 :: forall a.
Floating a =>
Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a)
nrm2 =
   forall a. Nrm a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a)
getNrm forall a b. (a -> b) -> a -> b
$
   forall a (f :: * -> *).
Floating a =>
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
Class.switchFloating
      (forall a. (Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a)) -> Nrm a
Nrm forall a. Real a => Ptr CInt -> Ptr a -> Ptr CInt -> IO a
BlasReal.nrm2) (forall a. (Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a)) -> Nrm a
Nrm forall a. Real a => Ptr CInt -> Ptr a -> Ptr CInt -> IO a
BlasReal.nrm2)
      (forall a. (Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a)) -> Nrm a
Nrm forall a. Real a => Ptr CInt -> Ptr (Complex a) -> Ptr CInt -> IO a
BlasComplex.cnrm2) (forall a. (Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a)) -> Nrm a
Nrm forall a. Real a => Ptr CInt -> Ptr (Complex a) -> Ptr CInt -> IO a
BlasComplex.cnrm2)

newtype Nrm a = Nrm {forall a. Nrm a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a)
getNrm :: Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a)}


newtype Norm f a = Norm {forall (f :: * -> *) a. Norm f a -> f a -> RealOf a
getNorm :: f a -> RealOf a}

norm2Squared :: (Shape.C sh, Class.Floating a) => T shA sh a -> RealOf a
norm2Squared :: forall sh a shA. (C sh, Floating a) => T shA sh a -> RealOf a
norm2Squared =
   forall (f :: * -> *) a. Norm f a -> f a -> RealOf a
getNorm forall a b. (a -> b) -> a -> b
$
   forall a (f :: * -> *).
Floating a =>
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
Class.switchFloating
      (forall (f :: * -> *) a. (f a -> RealOf a) -> Norm f a
Norm forall sh a shA. (C sh, Real a) => T shA sh a -> a
norm2SquaredReal)
      (forall (f :: * -> *) a. (f a -> RealOf a) -> Norm f a
Norm forall sh a shA. (C sh, Real a) => T shA sh a -> a
norm2SquaredReal)
      (forall (f :: * -> *) a. (f a -> RealOf a) -> Norm f a
Norm forall sh a shA. (C sh, Real a) => T shA sh (Complex a) -> a
norm2SquaredComplex)
      (forall (f :: * -> *) a. (f a -> RealOf a) -> Norm f a
Norm forall sh a shA. (C sh, Real a) => T shA sh (Complex a) -> a
norm2SquaredComplex)

norm2SquaredReal :: (Shape.C sh, Class.Real a) => T shA sh a -> a
norm2SquaredReal :: forall sh a shA. (C sh, Real a) => T shA sh a -> a
norm2SquaredReal T shA sh a
x =
   forall a. IO a -> a
unsafePerformIO forall a b. (a -> b) -> a -> b
$ forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT forall a b. (a -> b) -> a -> b
$ do
      (Ptr CInt
nPtr, Ptr a
sxPtr, Ptr CInt
incxPtr) <- forall sh a shA r.
(C sh, Storable a) =>
T shA sh a -> ContT r IO (Ptr CInt, Ptr a, Ptr CInt)
sizeSliceArg T shA sh a
x
      forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a.
Real a =>
Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO a
BlasReal.dot Ptr CInt
nPtr Ptr a
sxPtr Ptr CInt
incxPtr Ptr a
sxPtr Ptr CInt
incxPtr

norm2SquaredComplex :: (Shape.C sh, Class.Real a) => T shA sh (Complex a) -> a
norm2SquaredComplex :: forall sh a shA. (C sh, Real a) => T shA sh (Complex a) -> a
norm2SquaredComplex T shA sh (Complex a)
x =
   forall a. IO a -> a
unsafePerformIO forall a b. (a -> b) -> a -> b
$ forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT forall a b. (a -> b) -> a -> b
$ do
      Ptr CInt
nPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint forall a b. (a -> b) -> a -> b
$ forall sh. C sh => sh -> Int
Shape.size forall a b. (a -> b) -> a -> b
$ forall sh slice a. T sh slice a -> slice
shape T shA sh (Complex a)
x
      Ptr (Complex a)
xPtr <- forall a sh slice r.
Storable a =>
T sh slice a -> FortranIO r (Ptr a)
startArg T shA sh (Complex a)
x
      let xrPtr :: Ptr (RealOf (Complex a))
xrPtr = forall a. Ptr a -> Ptr (RealOf a)
realPtr Ptr (Complex a)
xPtr
      let xiPtr :: Ptr a
xiPtr = forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr (RealOf (Complex a))
xrPtr Int
1
      Ptr CInt
incxPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint (forall sh slice a. T sh slice a -> Int
increment T shA sh (Complex a)
x forall a. Num a => a -> a -> a
* Int
2)
      forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$
         forall (f :: * -> *) a b c.
Applicative f =>
(a -> b -> c) -> f a -> f b -> f c
liftA2 forall a. Num a => a -> a -> a
(+)
            (forall a.
Real a =>
Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO a
BlasReal.dot Ptr CInt
nPtr Ptr (RealOf (Complex a))
xrPtr Ptr CInt
incxPtr Ptr (RealOf (Complex a))
xrPtr Ptr CInt
incxPtr)
            (forall a.
Real a =>
Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO a
BlasReal.dot Ptr CInt
nPtr Ptr a
xiPtr Ptr CInt
incxPtr Ptr a
xiPtr Ptr CInt
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 :: forall sh a shA. (C sh, Floating a) => T shA sh a -> RealOf a
normInf T shA sh a
arr = forall a. IO a -> a
unsafePerformIO forall a b. (a -> b) -> a -> b
$
   forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (forall a. Floating a => a -> RealOf a
Scalar.absolute forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall b a. b -> (a -> b) -> Maybe a -> b
maybe forall a. Floating a => a
Scalar.zero forall a b. (a, b) -> b
snd) forall a b. (a -> b) -> a -> b
$ forall sh a shA.
(C sh, Floating a) =>
T shA sh a -> IO (Maybe (Int, a))
absMax T shA sh a
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 :: forall sh a shA. (C sh, Floating a) => T shA sh a -> RealOf a
normInf1 T shA sh a
x = forall a. IO a -> a
unsafePerformIO forall a b. (a -> b) -> a -> b
$
   forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT forall a b. (a -> b) -> a -> b
$ do
      (Ptr CInt
nPtr, Ptr a
sxPtr, Ptr CInt
incxPtr) <- forall sh a shA r.
(C sh, Storable a) =>
T shA sh a -> ContT r IO (Ptr CInt, Ptr a, Ptr CInt)
sizeSliceArg T shA sh a
x
      forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$
         forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (forall a. Floating a => a -> RealOf a
Scalar.norm1 forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall b a. b -> (a -> b) -> Maybe a -> b
maybe forall a. Floating a => a
Scalar.zero forall a b. (a, b) -> b
snd) forall a b. (a -> b) -> a -> b
$
         forall a. Storable a => Ptr a -> Int -> CInt -> IO (Maybe (Int, a))
peekElemOff1 Ptr a
sxPtr (forall sh slice a. T sh slice a -> Int
increment T shA sh a
x) forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
=<< forall a. Floating a => Ptr CInt -> Ptr a -> Ptr CInt -> IO CInt
Blas.iamax Ptr CInt
nPtr Ptr a
sxPtr Ptr CInt
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 :: forall sh a shA.
(InvIndexed sh, Floating a) =>
T shA sh a -> (Index sh, a)
argAbsMaximum T shA sh a
arr = forall a. IO a -> a
unsafePerformIO forall a b. (a -> b) -> a -> b
$
   forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap
      (forall b a. b -> (a -> b) -> Maybe a -> b
maybe
         (forall a. HasCallStack => String -> a
error String
"Vector.argAbsMaximum: empty vector")
         (forall a c b. (a -> c) -> (a, b) -> (c, b)
mapFst (forall sh. InvIndexed sh => sh -> Int -> Index sh
Shape.uncheckedIndexFromOffset forall a b. (a -> b) -> a -> b
$ forall sh slice a. T sh slice a -> slice
shape T shA sh a
arr))) forall a b. (a -> b) -> a -> b
$
   forall sh a shA.
(C sh, Floating a) =>
T shA sh a -> IO (Maybe (Int, a))
absMax T shA sh a
arr

absMax ::
   (Shape.C sh, Class.Floating a) =>
   T shA sh a -> IO (Maybe (Int, a))
absMax :: forall sh a shA.
(C sh, Floating a) =>
T shA sh a -> IO (Maybe (Int, a))
absMax T shA sh a
x =
   case forall a (f :: * -> *). Floating a => f a -> ComplexSingleton a
Scalar.complexSingletonOfFunctor T shA sh a
x of
      ComplexSingleton a
Scalar.Real -> forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT forall a b. (a -> b) -> a -> b
$ do
         (Ptr CInt
nPtr, Ptr a
sxPtr, Ptr CInt
incxPtr) <- forall sh a shA r.
(C sh, Storable a) =>
T shA sh a -> ContT r IO (Ptr CInt, Ptr a, Ptr CInt)
sizeSliceArg T shA sh a
x
         forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$
            forall a. Storable a => Ptr a -> Int -> CInt -> IO (Maybe (Int, a))
peekElemOff1 Ptr a
sxPtr (forall sh slice a. T sh slice a -> Int
increment T shA sh a
x) forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
=<< forall a. Floating a => Ptr CInt -> Ptr a -> Ptr CInt -> IO CInt
Blas.iamax Ptr CInt
nPtr Ptr a
sxPtr Ptr CInt
incxPtr
      ComplexSingleton a
Scalar.Complex -> forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT forall a b. (a -> b) -> a -> b
$ do
         let n :: Int
n = forall sh. C sh => sh -> Int
Shape.size forall a b. (a -> b) -> a -> b
$ forall sh slice a. T sh slice a -> slice
shape T shA sh a
x
         Ptr a
sxPtr <- forall a sh slice r.
Storable a =>
T sh slice a -> FortranIO r (Ptr a)
startArg T shA sh a
x
         let incx :: Int
incx = forall sh slice a. T sh slice a -> Int
increment T shA sh a
x
         forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a. Storable a => Ptr a -> Int -> CInt -> IO (Maybe (Int, a))
peekElemOff1 Ptr a
sxPtr Int
incx forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
=<< forall a. Real a => Int -> Ptr (Complex a) -> Int -> IO CInt
absMaxComplex Int
n Ptr a
sxPtr Int
incx

absMaxComplex :: (Class.Real a) => Int -> Ptr (Complex a) -> Int -> IO CInt
absMaxComplex :: forall a. Real a => Int -> Ptr (Complex a) -> Int -> IO CInt
absMaxComplex Int
n Ptr (Complex a)
sxPtr Int
incx =
   forall a b. Storable a => Int -> (Ptr a -> IO b) -> IO b
ForeignArray.alloca Int
n forall a b. (a -> b) -> a -> b
$ \Ptr a
syPtr -> do
      let xrPtr :: Ptr (RealOf (Complex a))
xrPtr = forall a. Ptr a -> Ptr (RealOf a)
realPtr Ptr (Complex a)
sxPtr
      let incx2 :: Int
incx2 = Int
2forall a. Num a => a -> a -> a
*Int
incx
      forall a.
Floating a =>
Conjugation
-> Int -> Ptr a -> Int -> Ptr a -> Int -> Ptr a -> Int -> IO ()
Private.mul    Conjugation
NonConjugated Int
n Ptr (RealOf (Complex a))
xrPtr Int
incx2 Ptr (RealOf (Complex a))
xrPtr Int
incx2 Ptr a
syPtr Int
1
      let xiPtr :: Ptr a
xiPtr = forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr (RealOf (Complex a))
xrPtr Int
1
      forall a.
Floating a =>
Conjugation
-> Int
-> Ptr a
-> Int
-> Ptr a
-> Int
-> a
-> Ptr a
-> Int
-> IO ()
Private.mulAdd Conjugation
NonConjugated Int
n Ptr a
xiPtr Int
incx2 Ptr a
xiPtr Int
incx2 forall a. Floating a => a
Scalar.one Ptr a
syPtr Int
1
      forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT forall a b. (a -> b) -> a -> b
$ do
         Ptr CInt
nPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
         Ptr CInt
incyPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1
         forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a. Floating a => Ptr CInt -> Ptr a -> Ptr CInt -> IO CInt
Blas.iamax Ptr CInt
nPtr Ptr a
syPtr Ptr CInt
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 :: forall sh a shA.
(InvIndexed sh, Floating a) =>
T shA sh a -> (Index sh, a)
argAbs1Maximum T shA sh a
x = forall a. IO a -> a
unsafePerformIO forall a b. (a -> b) -> a -> b
$
   forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT forall a b. (a -> b) -> a -> b
$ do
      (Ptr CInt
nPtr, Ptr a
sxPtr, Ptr CInt
incxPtr) <- forall sh a shA r.
(C sh, Storable a) =>
T shA sh a -> ContT r IO (Ptr CInt, Ptr a, Ptr CInt)
sizeSliceArg T shA sh a
x
      forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$
         forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap
            (forall b a. b -> (a -> b) -> Maybe a -> b
maybe
               (forall a. HasCallStack => String -> a
error String
"Vector.argAbs1Maximum: empty vector")
               (forall a c b. (a -> c) -> (a, b) -> (c, b)
mapFst (forall sh. InvIndexed sh => sh -> Int -> Index sh
Shape.uncheckedIndexFromOffset forall a b. (a -> b) -> a -> b
$ forall sh slice a. T sh slice a -> slice
shape T shA sh a
x))) forall a b. (a -> b) -> a -> b
$
         forall a. Storable a => Ptr a -> Int -> CInt -> IO (Maybe (Int, a))
peekElemOff1 Ptr a
sxPtr (forall sh slice a. T sh slice a -> Int
increment T shA sh a
x) forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
=<< forall a. Floating a => Ptr CInt -> Ptr a -> Ptr CInt -> IO CInt
Blas.iamax Ptr CInt
nPtr Ptr a
sxPtr Ptr CInt
incxPtr

peekElemOff1 :: (Storable a) => Ptr a -> Int -> CInt -> IO (Maybe (Int, a))
peekElemOff1 :: forall a. Storable a => Ptr a -> Int -> CInt -> IO (Maybe (Int, a))
peekElemOff1 Ptr a
ptr Int
inc CInt
k1 =
   let k1i :: Int
k1i = forall a b. (Integral a, Num b) => a -> b
fromIntegral CInt
k1
       ki :: Int
ki = Int
k1iforall a. Num a => a -> a -> a
-Int
1
   in if Int
k1i forall a. Eq a => a -> a -> Bool
== Int
0
         then forall (m :: * -> *) a. Monad m => a -> m a
return forall a. Maybe a
Nothing
         else forall a. a -> Maybe a
Just forall b c a. (b -> c) -> (a -> b) -> a -> c
. (,) Int
ki forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> forall a. Storable a => Ptr a -> Int -> IO a
peekElemOff Ptr a
ptr (Int
kiforall a. Num a => a -> a -> a
*Int
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 :: forall sh a shA. (C sh, Floating a) => T shA sh a -> a
product T shA sh a
x = forall a. IO a -> a
unsafePerformIO forall a b. (a -> b) -> a -> b
$ forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT forall a b. (a -> b) -> a -> b
$ do
   Ptr a
xPtr <- forall a sh slice r.
Storable a =>
T sh slice a -> FortranIO r (Ptr a)
startArg T shA sh a
x
   forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a. Floating a => Int -> Ptr a -> Int -> IO a
Private.product (forall sh. C sh => sh -> Int
Shape.size forall a b. (a -> b) -> a -> b
$ forall sh slice a. T sh slice a -> slice
shape T shA sh a
x) Ptr a
xPtr (forall sh slice a. T sh slice a -> Int
increment T shA sh a
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 :: forall shA sh a. (C shA, C sh, Real a) => T shA sh a -> a
minimum = forall a b. (a, b) -> a
fst forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall shA sh a. (C shA, C sh, Real a) => T shA sh a -> (a, a)
limits
maximum :: forall shA sh a. (C shA, C sh, Real a) => T shA sh a -> a
maximum = forall a b. (a, b) -> b
snd forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall shA sh a. (C shA, C sh, Real a) => T shA sh a -> (a, a)
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 :: forall shA sh ix a.
(C shA, InvIndexed sh, Index sh ~ ix, Real a) =>
T shA sh a -> (ix, a)
argMinimum = forall a b. (a, b) -> a
fst forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall shA sh ix a.
(C shA, InvIndexed sh, Index sh ~ ix, Real a) =>
T shA sh a -> ((ix, a), (ix, a))
argLimits
argMaximum :: forall shA sh ix a.
(C shA, InvIndexed sh, Index sh ~ ix, Real a) =>
T shA sh a -> (ix, a)
argMaximum = forall a b. (a, b) -> b
snd forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall shA sh ix a.
(C shA, InvIndexed sh, Index sh ~ ix, Real a) =>
T shA sh a -> ((ix, a), (ix, a))
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 :: forall shA sh a. (C shA, C sh, Real a) => T shA sh a -> (a, a)
limits T shA sh a
xs0 =
   let xs :: T shA (Deferred sh) a
xs = forall slice0 slice1 sh a.
(slice0 -> slice1) -> T sh slice0 a -> T sh slice1 a
mapShape forall sh. sh -> Deferred sh
Shape.Deferred T shA sh a
xs0
       x0 :: a
x0 = forall a b. (a, b) -> b
snd forall a b. (a -> b) -> a -> b
$ forall sh a shA.
(InvIndexed sh, Floating a) =>
T shA sh a -> (Index sh, a)
argAbs1Maximum T shA (Deferred sh) a
xs
       x1 :: a
x1 = T shA (Deferred sh) a
xs forall shA sh a.
(C shA, Indexed sh, Storable a) =>
T shA sh a -> Index sh -> a
! forall a b. (a, b) -> a
fst (forall sh a shA.
(InvIndexed sh, Floating a) =>
T shA sh a -> (Index sh, a)
argAbs1Maximum (forall sh a. C sh => Vector sh a -> T sh sh a
fromVector (forall sh a shA.
(C sh, Floating a) =>
a -> T shA sh a -> Vector sh a
raise (-a
x0) T shA (Deferred sh) a
xs)))
   in if a
x0forall a. Ord a => a -> a -> Bool
>=a
0 then (a
x1,a
x0) else (a
x0,a
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 :: forall shA sh ix a.
(C shA, InvIndexed sh, Index sh ~ ix, Real a) =>
T shA sh a -> ((ix, a), (ix, a))
argLimits T shA sh a
xs =
   let p0 :: (Index sh, a)
p0@(Index sh
_i0,a
x0) = forall sh a shA.
(InvIndexed sh, Floating a) =>
T shA sh a -> (Index sh, a)
argAbs1Maximum T shA sh a
xs
       p1 :: (ix, a)
p1 = (ix
i1,T shA sh a
xsforall shA sh a.
(C shA, Indexed sh, Storable a) =>
T shA sh a -> Index sh -> a
!ix
i1); i1 :: ix
i1 = forall a b. (a, b) -> a
fst forall a b. (a -> b) -> a -> b
$ forall sh a shA.
(InvIndexed sh, Floating a) =>
T shA sh a -> (Index sh, a)
argAbs1Maximum forall a b. (a -> b) -> a -> b
$ forall sh a. C sh => Vector sh a -> T sh sh a
fromVector forall a b. (a -> b) -> a -> b
$ forall sh a shA.
(C sh, Floating a) =>
a -> T shA sh a -> Vector sh a
raise (-a
x0) T shA sh a
xs
   in if a
x0forall a. Ord a => a -> a -> Bool
>=a
0 then ((ix, a)
p1,(Index sh, a)
p0) else ((Index sh, a)
p0,(ix, a)
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 :: forall sh a shA.
(C sh, Floating a) =>
a -> T shA sh a -> Vector sh a
scale a
alpha T shA sh a
x = forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize (forall sh slice a. T sh slice a -> slice
shape T shA sh a
x) forall a b. (a -> b) -> a -> b
$ \Int
n Ptr a
syPtr -> do
   forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT forall a b. (a -> b) -> a -> b
$ do
      Ptr a
alphaPtr <- forall a r. Floating a => a -> FortranIO r (Ptr a)
Call.number a
alpha
      Ptr CInt
nPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
      (Ptr a
sxPtr, Ptr CInt
incxPtr) <- forall a sh slice r.
Storable a =>
T sh slice a -> FortranIO r (Ptr a, Ptr CInt)
sliceArg T shA sh a
x
      Ptr CInt
incyPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1
      forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a.
Floating a =>
Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
Blas.copy Ptr CInt
nPtr Ptr a
sxPtr Ptr CInt
incxPtr Ptr a
syPtr Ptr CInt
incyPtr
      forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a.
Floating a =>
Ptr CInt -> Ptr a -> Ptr a -> Ptr CInt -> IO ()
Blas.scal Ptr CInt
nPtr Ptr a
alphaPtr Ptr a
syPtr Ptr CInt
incyPtr

_scale :: forall sh a shA.
(C sh, Floating a) =>
a -> T shA sh a -> Vector sh a
_scale a
a T shA sh a
b = forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize (forall sh slice a. T sh slice a -> slice
shape T shA sh a
b) forall a b. (a -> b) -> a -> b
$ \Int
n Ptr a
cPtr -> do
   let m :: Int
m = Int
1
   let k :: Int
k = Int
1
   forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT forall a b. (a -> b) -> a -> b
$ do
      Ptr CChar
transaPtr <- forall r. Char -> FortranIO r (Ptr CChar)
Call.char Char
'N'
      Ptr CChar
transbPtr <- forall r. Char -> FortranIO r (Ptr CChar)
Call.char Char
'N'
      Ptr CInt
mPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
m
      Ptr CInt
kPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
k
      Ptr CInt
nPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
      Ptr a
alphaPtr <- forall a r. Floating a => a -> FortranIO r (Ptr a)
Call.number forall a. Floating a => a
Scalar.one
      Ptr a
aPtr <- forall a r. Floating a => a -> FortranIO r (Ptr a)
Call.number a
a
      Ptr CInt
ldaPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
m
      (Ptr a
bPtr, Ptr CInt
ldbPtr) <- forall a sh slice r.
Storable a =>
T sh slice a -> FortranIO r (Ptr a, Ptr CInt)
sliceArg T shA sh a
b
      Ptr a
betaPtr <- forall a r. Floating a => a -> FortranIO r (Ptr a)
Call.number forall a. Floating a => a
Scalar.zero
      Ptr CInt
ldcPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
m
      forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$
         forall a.
Floating a =>
Ptr CChar
-> Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> IO ()
Blas.gemm
            Ptr CChar
transaPtr Ptr CChar
transbPtr Ptr CInt
mPtr Ptr CInt
nPtr Ptr CInt
kPtr Ptr a
alphaPtr
            Ptr a
aPtr Ptr CInt
ldaPtr Ptr a
bPtr Ptr CInt
ldbPtr Ptr a
betaPtr Ptr a
cPtr Ptr CInt
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 :: forall sh a shA.
(C sh, Floating a) =>
RealOf a -> T shA sh a -> Vector sh a
scaleReal =
   forall (f :: * -> *) (g :: * -> *) a.
ScaleReal f g a -> RealOf a -> f a -> g a
getScaleReal forall a b. (a -> b) -> a -> b
$
   forall a (f :: * -> *).
Floating a =>
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
Class.switchFloating
      (forall (f :: * -> *) (g :: * -> *) a.
(RealOf a -> f a -> g a) -> ScaleReal f g a
ScaleReal forall sh a shA.
(C sh, Floating a) =>
a -> T shA sh a -> Vector sh a
scale)
      (forall (f :: * -> *) (g :: * -> *) a.
(RealOf a -> f a -> g a) -> ScaleReal f g a
ScaleReal forall sh a shA.
(C sh, Floating a) =>
a -> T shA sh a -> Vector sh a
scale)
      (forall (f :: * -> *) (g :: * -> *) a.
(RealOf a -> f a -> g a) -> ScaleReal f g a
ScaleReal forall a b. (a -> b) -> a -> b
$ forall sh a shA.
(C sh, Floating a) =>
a -> T shA sh a -> Vector sh a
scale forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Floating a => RealOf a -> a
Scalar.fromReal)
      (forall (f :: * -> *) (g :: * -> *) a.
(RealOf a -> f a -> g a) -> ScaleReal f g a
ScaleReal forall a b. (a -> b) -> a -> b
$ forall sh a shA.
(C sh, Floating a) =>
a -> T shA sh a -> Vector sh a
scale forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Floating a => RealOf a -> a
Scalar.fromReal)

newtype ScaleReal f g a = ScaleReal {forall (f :: * -> *) (g :: * -> *) a.
ScaleReal f g a -> RealOf a -> f a -> g a
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 :: forall sh a shA shB.
(C sh, Eq sh, Floating a) =>
T shA sh a -> T shB sh a -> Vector sh a
add = forall sh a shA shB.
(C sh, Eq sh, Floating a) =>
a -> T shA sh a -> T shB sh a -> Vector sh a
mac forall a. Floating a => a
Scalar.one
sub :: forall sh a shA shB.
(C sh, Eq sh, Floating a) =>
T shA sh a -> T shB sh a -> Vector sh a
sub T shA sh a
x T shB sh a
y = forall sh a shA shB.
(C sh, Eq sh, Floating a) =>
a -> T shA sh a -> T shB sh a -> Vector sh a
mac forall a. Floating a => a
Scalar.minusOne T shB sh a
y T shA sh a
x

mac ::
   (Shape.C sh, Eq sh, Class.Floating a) =>
   a -> T shA sh a -> T shB sh a -> Vector sh a
mac :: forall sh a shA shB.
(C sh, Eq sh, Floating a) =>
a -> T shA sh a -> T shB sh a -> Vector sh a
mac a
alpha T shA sh a
x T shB sh a
y =
   forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize (forall sh slice a. T sh slice a -> slice
shape T shA sh a
x) forall a b. (a -> b) -> a -> b
$ \Int
n Ptr a
szPtr -> do
   String -> Bool -> IO ()
Call.assert String
"mac: shapes mismatch" (forall sh slice a. T sh slice a -> slice
shape T shA sh a
x forall a. Eq a => a -> a -> Bool
== forall sh slice a. T sh slice a -> slice
shape T shB sh a
y)
   forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT forall a b. (a -> b) -> a -> b
$ do
      Ptr CInt
nPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
      Ptr a
saPtr <- forall a r. Floating a => a -> FortranIO r (Ptr a)
Call.number a
alpha
      (Ptr a
sxPtr, Ptr CInt
incxPtr) <- forall a sh slice r.
Storable a =>
T sh slice a -> FortranIO r (Ptr a, Ptr CInt)
sliceArg T shA sh a
x
      (Ptr a
syPtr, Ptr CInt
incyPtr) <- forall a sh slice r.
Storable a =>
T sh slice a -> FortranIO r (Ptr a, Ptr CInt)
sliceArg T shB sh a
y
      Ptr CInt
inczPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1
      forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a.
Floating a =>
Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
Blas.copy Ptr CInt
nPtr Ptr a
syPtr Ptr CInt
incyPtr Ptr a
szPtr Ptr CInt
inczPtr
      forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a.
Floating a =>
Ptr CInt
-> Ptr a -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
Blas.axpy Ptr CInt
nPtr Ptr a
saPtr Ptr a
sxPtr Ptr CInt
incxPtr Ptr a
szPtr Ptr CInt
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 :: forall slice a sh.
(C slice, Floating a) =>
T sh slice a -> Vector slice a
negate =
   forall (f :: * -> *) (g :: * -> *) a. Conjugate f g a -> f a -> g a
getConjugate forall a b. (a -> b) -> a -> b
$
   forall a (f :: * -> *).
Floating a =>
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
Class.switchFloating
      (forall (f :: * -> *) (g :: * -> *) a.
(f a -> g a) -> Conjugate f g a
Conjugate forall a b. (a -> b) -> a -> b
$ forall sh a shA.
(C sh, Floating a) =>
RealOf a -> T shA sh a -> Vector sh a
scaleReal forall a. Floating a => a
Scalar.minusOne)
      (forall (f :: * -> *) (g :: * -> *) a.
(f a -> g a) -> Conjugate f g a
Conjugate forall a b. (a -> b) -> a -> b
$ forall sh a shA.
(C sh, Floating a) =>
RealOf a -> T shA sh a -> Vector sh a
scaleReal forall a. Floating a => a
Scalar.minusOne)
      (forall (f :: * -> *) (g :: * -> *) a.
(f a -> g a) -> Conjugate f g a
Conjugate forall a b. (a -> b) -> a -> b
$ forall sh a shA.
(C sh, Floating a) =>
RealOf a -> T shA sh a -> Vector sh a
scaleReal forall a. Floating a => a
Scalar.minusOne)
      (forall (f :: * -> *) (g :: * -> *) a.
(f a -> g a) -> Conjugate f g a
Conjugate forall a b. (a -> b) -> a -> b
$ forall sh a shA.
(C sh, Floating a) =>
RealOf a -> T shA sh a -> Vector sh a
scaleReal forall a. Floating a => a
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 :: forall sh a shA.
(C sh, Floating a) =>
a -> T shA sh a -> Vector sh a
raise a
c T shA sh a
x =
   forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize (forall sh slice a. T sh slice a -> slice
shape T shA sh a
x) forall a b. (a -> b) -> a -> b
$ \Int
n Ptr a
yPtr -> forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT forall a b. (a -> b) -> a -> b
$ do
      Ptr CInt
nPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
      Ptr a
cPtr <- forall a r. Floating a => a -> FortranIO r (Ptr a)
Call.number a
c
      Ptr a
onePtr <- forall a r. Floating a => a -> FortranIO r (Ptr a)
Call.number forall a. Floating a => a
Scalar.one
      Ptr CInt
inccPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
0
      (Ptr a
xPtr, Ptr CInt
incxPtr) <- forall a sh slice r.
Storable a =>
T sh slice a -> FortranIO r (Ptr a, Ptr CInt)
sliceArg T shA sh a
x
      Ptr CInt
inc1Ptr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1
      forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ do
         forall a.
Floating a =>
Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
Blas.copy Ptr CInt
nPtr Ptr a
xPtr Ptr CInt
incxPtr Ptr a
yPtr Ptr CInt
inc1Ptr
         forall a.
Floating a =>
Ptr CInt
-> Ptr a -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
Blas.axpy Ptr CInt
nPtr Ptr a
onePtr Ptr a
cPtr Ptr CInt
inccPtr Ptr a
yPtr Ptr CInt
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 :: forall sh a shA shB.
(C sh, Eq sh, Floating a) =>
T shA sh a -> T shB sh a -> Vector sh a
mul = forall sh a shA shB.
(C sh, Eq sh, Floating a) =>
Conjugation -> T shA sh a -> T shB sh a -> Vector sh a
mulConjugation Conjugation
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 :: forall sh a shA shB.
(C sh, Eq sh, Floating a) =>
T shA sh a -> T shB sh a -> Vector sh a
mulConj = forall sh a shA shB.
(C sh, Eq sh, Floating a) =>
Conjugation -> T shA sh a -> T shB sh a -> Vector sh a
mulConjugation Conjugation
Conjugated

mulConjugation ::
   (Shape.C sh, Eq sh, Class.Floating a) =>
   Conjugation -> T shA sh a -> T shB sh a -> Vector sh a
mulConjugation :: forall sh a shA shB.
(C sh, Eq sh, Floating a) =>
Conjugation -> T shA sh a -> T shB sh a -> Vector sh a
mulConjugation Conjugation
conj T shA sh a
a T shB sh a
x =
      forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize (forall sh slice a. T sh slice a -> slice
shape T shB sh a
x) forall a b. (a -> b) -> a -> b
$ \Int
n Ptr a
yPtr -> do
   String -> Bool -> IO ()
Call.assert String
"mul: shapes mismatch" (forall sh slice a. T sh slice a -> slice
shape T shA sh a
a forall a. Eq a => a -> a -> Bool
== forall sh slice a. T sh slice a -> slice
shape T shB sh a
x)
   forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT forall a b. (a -> b) -> a -> b
$ do
      Ptr a
aPtr <- forall a sh slice r.
Storable a =>
T sh slice a -> FortranIO r (Ptr a)
startArg T shA sh a
a
      Ptr a
xPtr <- forall a sh slice r.
Storable a =>
T sh slice a -> FortranIO r (Ptr a)
startArg T shB sh a
x
      forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a.
Floating a =>
Conjugation
-> Int -> Ptr a -> Int -> Ptr a -> Int -> Ptr a -> Int -> IO ()
Private.mul Conjugation
conj Int
n Ptr a
aPtr (forall sh slice a. T sh slice a -> Int
increment T shA sh a
a) Ptr a
xPtr (forall sh slice a. T sh slice a -> Int
increment T shB sh a
x) Ptr a
yPtr Int
1


newtype Conjugate f g a = Conjugate {forall (f :: * -> *) (g :: * -> *) a. Conjugate f g a -> f a -> g a
getConjugate :: f a -> g a}

conjugate ::
   (Shape.C sh, Class.Floating a) =>
   T shA sh a -> Vector sh a
conjugate :: forall slice a sh.
(C slice, Floating a) =>
T sh slice a -> Vector slice a
conjugate =
   forall (f :: * -> *) (g :: * -> *) a. Conjugate f g a -> f a -> g a
getConjugate forall a b. (a -> b) -> a -> b
$
   forall a (f :: * -> *).
Floating a =>
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
Class.switchFloating
      (forall (f :: * -> *) (g :: * -> *) a.
(f a -> g a) -> Conjugate f g a
Conjugate forall slice a sh.
(C slice, Floating a) =>
T sh slice a -> Vector slice a
toVector)
      (forall (f :: * -> *) (g :: * -> *) a.
(f a -> g a) -> Conjugate f g a
Conjugate forall slice a sh.
(C slice, Floating a) =>
T sh slice a -> Vector slice a
toVector)
      (forall (f :: * -> *) (g :: * -> *) a.
(f a -> g a) -> Conjugate f g a
Conjugate forall sh a shA.
(C sh, Real a) =>
T shA sh (Complex a) -> Vector sh (Complex a)
complexConjugate)
      (forall (f :: * -> *) (g :: * -> *) a.
(f a -> g a) -> Conjugate f g a
Conjugate forall sh a shA.
(C sh, Real a) =>
T shA sh (Complex a) -> Vector sh (Complex a)
complexConjugate)

complexConjugate ::
   (Shape.C sh, Class.Real a) =>
   T shA sh (Complex a) -> Vector sh (Complex a)
complexConjugate :: forall sh a shA.
(C sh, Real a) =>
T shA sh (Complex a) -> Vector sh (Complex a)
complexConjugate T shA sh (Complex a)
x = forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize (forall sh slice a. T sh slice a -> slice
shape T shA sh (Complex a)
x) forall a b. (a -> b) -> a -> b
$ \Int
n Ptr (Complex a)
syPtr ->
   forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT forall a b. (a -> b) -> a -> b
$ do
      Ptr CInt
nPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
      (Ptr (Complex a)
sxPtr, Ptr CInt
incxPtr) <- forall a sh slice r.
Storable a =>
T sh slice a -> FortranIO r (Ptr a, Ptr CInt)
sliceArg T shA sh (Complex a)
x
      Ptr CInt
incyPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1
      forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a.
Floating a =>
Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
copyConjugate Ptr CInt
nPtr Ptr (Complex a)
sxPtr Ptr CInt
incxPtr Ptr (Complex a)
syPtr Ptr CInt
incyPtr


fromReal ::
   (Shape.C sh, Class.Floating a) => T shA sh (RealOf a) -> Vector sh a
fromReal :: forall sh a shA.
(C sh, Floating a) =>
T shA sh (RealOf a) -> Vector sh a
fromReal =
   forall (f :: * -> *) (g :: * -> *) a.
FromReal f g a -> f (RealOf a) -> g a
getFromReal forall a b. (a -> b) -> a -> b
$
   forall a (f :: * -> *).
Floating a =>
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
Class.switchFloating
      (forall (f :: * -> *) (g :: * -> *) a.
(f (RealOf a) -> g a) -> FromReal f g a
FromReal forall slice a sh.
(C slice, Floating a) =>
T sh slice a -> Vector slice a
toVector)
      (forall (f :: * -> *) (g :: * -> *) a.
(f (RealOf a) -> g a) -> FromReal f g a
FromReal forall slice a sh.
(C slice, Floating a) =>
T sh slice a -> Vector slice a
toVector)
      (forall (f :: * -> *) (g :: * -> *) a.
(f (RealOf a) -> g a) -> FromReal f g a
FromReal forall sh a shA.
(C sh, Real a) =>
T shA sh a -> Vector sh (Complex a)
complexFromReal)
      (forall (f :: * -> *) (g :: * -> *) a.
(f (RealOf a) -> g a) -> FromReal f g a
FromReal forall sh a shA.
(C sh, Real a) =>
T shA sh a -> Vector sh (Complex a)
complexFromReal)

newtype FromReal f g a = FromReal {forall (f :: * -> *) (g :: * -> *) a.
FromReal f g a -> f (RealOf a) -> g a
getFromReal :: f (RealOf a) -> g a}

toComplex ::
   (Shape.C sh, Class.Floating a) => T shA sh a -> Vector sh (ComplexOf a)
toComplex :: forall sh a shA.
(C sh, Floating a) =>
T shA sh a -> Vector sh (ComplexOf a)
toComplex =
   forall (f :: * -> *) (g :: * -> *) a.
ToComplex f g a -> f a -> g (ComplexOf a)
getToComplex forall a b. (a -> b) -> a -> b
$
   forall a (f :: * -> *).
Floating a =>
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
Class.switchFloating
      (forall (f :: * -> *) (g :: * -> *) a.
(f a -> g (ComplexOf a)) -> ToComplex f g a
ToComplex forall sh a shA.
(C sh, Real a) =>
T shA sh a -> Vector sh (Complex a)
complexFromReal)
      (forall (f :: * -> *) (g :: * -> *) a.
(f a -> g (ComplexOf a)) -> ToComplex f g a
ToComplex forall sh a shA.
(C sh, Real a) =>
T shA sh a -> Vector sh (Complex a)
complexFromReal)
      (forall (f :: * -> *) (g :: * -> *) a.
(f a -> g (ComplexOf a)) -> ToComplex f g a
ToComplex forall slice a sh.
(C slice, Floating a) =>
T sh slice a -> Vector slice a
toVector)
      (forall (f :: * -> *) (g :: * -> *) a.
(f a -> g (ComplexOf a)) -> ToComplex f g a
ToComplex forall slice a sh.
(C slice, Floating a) =>
T sh slice a -> Vector slice a
toVector)

newtype ToComplex f g a = ToComplex {forall (f :: * -> *) (g :: * -> *) a.
ToComplex f g a -> f a -> g (ComplexOf a)
getToComplex :: f a -> g (ComplexOf a)}

complexFromReal ::
   (Shape.C sh, Class.Real a) => T shA sh a -> Vector sh (Complex a)
complexFromReal :: forall sh a shA.
(C sh, Real a) =>
T shA sh a -> Vector sh (Complex a)
complexFromReal T shA sh a
x =
   forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize (forall sh slice a. T sh slice a -> slice
shape T shA sh a
x) forall a b. (a -> b) -> a -> b
$ \Int
n Ptr (Complex a)
yPtr ->
   case forall a. Ptr a -> Ptr (RealOf a)
realPtr Ptr (Complex a)
yPtr of
      Ptr (RealOf (Complex a))
yrPtr -> forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT forall a b. (a -> b) -> a -> b
$ do
         Ptr CInt
nPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
         (Ptr a
xPtr, Ptr CInt
incxPtr) <- forall a sh slice r.
Storable a =>
T sh slice a -> FortranIO r (Ptr a, Ptr CInt)
sliceArg T shA sh a
x
         Ptr CInt
incyPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
2
         Ptr CInt
inczPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
0
         Ptr a
zPtr <- forall a r. Floating a => a -> FortranIO r (Ptr a)
Call.number forall a. Floating a => a
Scalar.zero
         forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ do
            forall a.
Floating a =>
Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
Blas.copy Ptr CInt
nPtr Ptr a
xPtr Ptr CInt
incxPtr Ptr (RealOf (Complex a))
yrPtr Ptr CInt
incyPtr
            forall a.
Floating a =>
Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
Blas.copy Ptr CInt
nPtr Ptr a
zPtr Ptr CInt
inczPtr (forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr (RealOf (Complex a))
yrPtr Int
1) Ptr CInt
incyPtr


realFromComplexVector ::
   (Shape.C sh) =>
   Vector sh (Complex a) -> T (sh, ComplexShape) (sh, ComplexShape) a
realFromComplexVector :: forall sh a.
C sh =>
Vector sh (Complex a) -> T (sh, ComplexShape) (sh, ComplexShape) a
realFromComplexVector (Array sh
sh ForeignPtr (Complex a)
a) =
   let csh :: (sh, ComplexShape)
csh = (sh
sh, forall sh. Static sh => sh
Shape.static) in
   forall sh slice a. T slice -> Vector sh a -> T sh slice a
Cons (forall sh. C sh => sh -> T sh
Slice.fromShape (sh, ComplexShape)
csh) (forall sh a. sh -> ForeignPtr a -> Array sh a
Array (sh, ComplexShape)
csh (forall a b. ForeignPtr a -> ForeignPtr b
castForeignPtr ForeignPtr (Complex a)
a))


realPart ::
   (Shape.C sh, Class.Real a) =>
   T shA sh (Complex a) -> T (shA, ComplexShape) sh a
realPart :: forall sh a shA.
(C sh, Real a) =>
T shA sh (Complex a) -> T (shA, ComplexShape) sh a
realPart (Cons (Slice.Cons Int
s Int
k sh
slc) (Array shA
sh ForeignPtr (Complex a)
a)) =
   forall sh slice a. T slice -> Vector sh a -> T sh slice a
Cons
      (forall sh. Int -> Int -> sh -> T sh
Slice.Cons (Int
2forall a. Num a => a -> a -> a
*Int
s) (Int
2forall a. Num a => a -> a -> a
*Int
k) sh
slc)
      (forall sh a. sh -> ForeignPtr a -> Array sh a
Array (shA
sh, forall sh. Static sh => sh
Shape.static) (forall a b. ForeignPtr a -> ForeignPtr b
castForeignPtr ForeignPtr (Complex a)
a))

imaginaryPart ::
   (Shape.C sh, Class.Real a) =>
   T shA sh (Complex a) -> T (shA, ComplexShape) sh a
imaginaryPart :: forall sh a shA.
(C sh, Real a) =>
T shA sh (Complex a) -> T (shA, ComplexShape) sh a
imaginaryPart (Cons (Slice.Cons Int
s Int
k sh
slc) (Array shA
sh ForeignPtr (Complex a)
a)) =
   forall sh slice a. T slice -> Vector sh a -> T sh slice a
Cons
      (forall sh. Int -> Int -> sh -> T sh
Slice.Cons (Int
2forall a. Num a => a -> a -> a
*Int
sforall a. Num a => a -> a -> a
+Int
1) (Int
2forall a. Num a => a -> a -> a
*Int
k) sh
slc)
      (forall sh a. sh -> ForeignPtr a -> Array sh a
Array (shA
sh, forall sh. Static sh => sh
Shape.static) (forall a b. ForeignPtr a -> ForeignPtr b
castForeignPtr ForeignPtr (Complex a)
a))


zipComplex ::
   (Shape.C sh, Eq sh, Class.Real a) =>
   T shA sh a -> T shB sh a -> Vector sh (Complex a)
zipComplex :: forall sh a shA shB.
(C sh, Eq sh, Real a) =>
T shA sh a -> T shB sh a -> Vector sh (Complex a)
zipComplex T shA sh a
xr T shB sh a
xi =
   forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize (forall sh slice a. T sh slice a -> slice
shape T shA sh a
xr) forall a b. (a -> b) -> a -> b
$ \Int
n Ptr (Complex a)
yPtr -> forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT forall a b. (a -> b) -> a -> b
$ do
      forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ String -> Bool -> IO ()
Call.assert String
"zipComplex: shapes mismatch" (forall sh slice a. T sh slice a -> slice
shape T shA sh a
xr forall a. Eq a => a -> a -> Bool
== forall sh slice a. T sh slice a -> slice
shape T shB sh a
xi)
      Ptr CInt
nPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
      (Ptr a
xrPtr, Ptr CInt
incxrPtr) <- forall a sh slice r.
Storable a =>
T sh slice a -> FortranIO r (Ptr a, Ptr CInt)
sliceArg T shA sh a
xr
      (Ptr a
xiPtr, Ptr CInt
incxiPtr) <- forall a sh slice r.
Storable a =>
T sh slice a -> FortranIO r (Ptr a, Ptr CInt)
sliceArg T shB sh a
xi
      let yrPtr :: Ptr (RealOf (Complex a))
yrPtr = forall a. Ptr a -> Ptr (RealOf a)
realPtr Ptr (Complex a)
yPtr
      Ptr CInt
incyPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
2
      forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ do
         forall a.
Floating a =>
Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
Blas.copy Ptr CInt
nPtr Ptr a
xrPtr Ptr CInt
incxrPtr Ptr (RealOf (Complex a))
yrPtr Ptr CInt
incyPtr
         forall a.
Floating a =>
Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
Blas.copy Ptr CInt
nPtr Ptr a
xiPtr Ptr CInt
incxiPtr (forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr (RealOf (Complex a))
yrPtr Int
1) Ptr CInt
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 :: forall sh a shA.
(C sh, Real a) =>
T shA sh (Complex a)
-> (T (shA, ComplexShape) sh a, T (shA, ComplexShape) sh a)
unzipComplex T shA sh (Complex a)
x = (forall sh a shA.
(C sh, Real a) =>
T shA sh (Complex a) -> T (shA, ComplexShape) sh a
realPart T shA sh (Complex a)
x, forall sh a shA.
(C sh, Real a) =>
T shA sh (Complex a) -> T (shA, ComplexShape) sh a
imaginaryPart T shA sh (Complex a)
x)