{-# 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,
   reverse,
   cyclicReverse,
   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, poke, 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 (when, 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, Integral, fromIntegral, (-), Char, error, IO)


{- $setup
>>> import Test.NumberModule.Type (Number_)
>>> import Test.Generator (genNumber)
>>> import Test.Slice (shapeInt)
>>> 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 Data.Monoid ((<>))
>>> import Data.Tuple.HT (mapPair)
>>> 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 ::
>>>    (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
>>>
>>> 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 => Array sh a -> Bool
>>> isNonEmpty xs = Shape.size (Array.shape xs) > 0
>>>
>>> forVector ::
>>>    (QC.Testable prop, QC.Arbitrary a, Class.Floating a, Show a) =>
>>>    QC.Gen a -> (Vector a -> prop) -> QC.Property
>>> forVector genElem =
>>>    QC.forAllShrink
>>>       (flip genVector genElem . shapeInt =<< QC.choose (0,maxDim))
>>>       (map Vector.autoFromList . QC.shrink . Vector.toList)
>>>
>>> forVector2 ::
>>>    (QC.Testable prop, QC.Arbitrary a, Class.Floating a, Show a) =>
>>>    QC.Gen a -> (Vector a -> Vector a -> prop) -> QC.Property
>>> forVector2 genElem prop =
>>>    QC.forAllShrink
>>>       (do shape <- fmap shapeInt $ QC.choose (0,maxDim)
>>>           liftA2 (,) (genVector shape genElem) (genVector shape genElem))
>>>       (map (mapPair (Vector.autoFromList, Vector.autoFromList) . unzip) .
>>>        QC.shrink .
>>>        uncurry zip . mapPair (Vector.toList, Vector.toList))
>>>       (uncurry prop)
>>>
>>> type CyclicVector = Vector.Vector (Shape.Cyclic Int)
>>>
>>> genCyclicVector ::
>>>    (Class.Floating a) =>
>>>    Integer -> Int -> QC.Gen (CyclicVector a)
>>> genCyclicVector maxE dim =
>>>    fmap (Vector.fromList (Shape.Cyclic dim)) $
>>>    replicateM dim $ genNumber maxE
>>>
>>> cyclicVectorFromListGen :: (Class.Floating a) => [a] -> CyclicVector a
>>> cyclicVectorFromListGen xs = Vector.fromList (Shape.Cyclic $ length xs) xs
>>>
>>> cyclicVectorFromList :: [Number_] -> CyclicVector Number_
>>> cyclicVectorFromList = cyclicVectorFromListGen
>>>
>>> forCyclicVector ::
>>>    (QC.Testable prop, QC.Arbitrary a, Class.Floating a, Show a) =>
>>>    QC.Gen a -> (CyclicVector a -> prop) -> QC.Property
>>> forCyclicVector genElem =
>>>    QC.forAllShrink
>>>       (flip genVector genElem . Shape.Cyclic =<< QC.choose (0,maxDim))
>>>       (map cyclicVectorFromListGen . QC.shrink . Vector.toList)
-}


type Vector = Array


toList :: (Shape.C sh, Storable a) => Vector sh a -> [a]
toList :: forall sh a. (C sh, Storable a) => Vector sh a -> [a]
toList = forall sh a. (C sh, Storable a) => Vector sh a -> [a]
Array.toList

fromList :: (Shape.C sh, Storable a) => sh -> [a] -> Vector sh a
fromList :: forall sh a. (C sh, Storable a) => sh -> [a] -> Vector sh a
fromList = forall sh a. (C sh, Storable a) => sh -> [a] -> Vector sh a
CheckedArray.fromList

autoFromList :: (Storable a) => [a] -> Vector (Shape.ZeroBased Int) a
autoFromList :: forall a. Storable a => [a] -> Vector (ZeroBased Int) a
autoFromList = forall a. Storable a => [a] -> Vector (ZeroBased Int) a
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 :: forall sh a. (C sh, Floating a) => sh -> a -> Vector sh a
constant sh
sh a
a = forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize sh
sh forall a b. (a -> b) -> a -> b
$ forall a. Floating a => a -> Int -> Ptr a -> IO ()
fill a
a

zero :: (Shape.C sh, Class.Floating a) => sh -> Vector sh a
zero :: forall sh a. (C sh, Floating a) => sh -> Vector sh a
zero = forall a b c. (a -> b -> c) -> b -> a -> c
flip forall sh a. (C sh, Floating a) => sh -> a -> Vector sh a
constant forall a. Floating a => a
Scalar.zero

one :: (Shape.C sh, Class.Floating a) => sh -> Vector sh a
one :: forall sh a. (C sh, Floating a) => sh -> Vector sh a
one = forall a b c. (a -> b -> c) -> b -> a -> c
flip forall sh a. (C sh, Floating a) => sh -> a -> Vector sh a
constant forall a. Floating a => a
Scalar.one

unit ::
   (Shape.Indexed sh, Class.Floating a) =>
   sh -> Shape.Index sh -> Vector sh a
unit :: forall sh a.
(Indexed sh, Floating a) =>
sh -> Index sh -> Vector sh a
unit sh
sh Index sh
ix = forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize sh
sh forall a b. (a -> b) -> a -> b
$ \Int
n Ptr a
xPtr -> do
   forall a. Floating a => a -> Int -> Ptr a -> IO ()
fill forall a. Floating a => a
Scalar.zero Int
n Ptr a
xPtr
   forall a. Storable a => Ptr a -> Int -> a -> IO ()
pokeElemOff Ptr a
xPtr (forall sh. Indexed sh => sh -> Index sh -> Int
Shape.offset sh
sh Index sh
ix) forall a. Floating a => a
Scalar.one


{- |
Precedence and associativity (right) of (List.++).
This also matches '(::+)'.
-}
infixr 5 +++


{- |
prop> forVector number_ $ \xs -> forVector number_ $ \ys -> forVector number_ $ \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
+++ :: forall shx shy a.
(C shx, C shy, Storable a) =>
Vector shx a -> Vector shy a -> Vector (shx ::+ shy) a
(+++) = forall shx shy a.
(C shx, C shy, Storable a) =>
Vector shx a -> Vector shy a -> Vector (shx ::+ shy) a
append


{- |
prop> Vector.autoFromList [] == (Vector.reverse $ Vector.autoFromList [] :: Vector Number_)
prop> Vector.autoFromList [1] == (Vector.reverse $ Vector.autoFromList [1] :: Vector Number_)
prop> Vector.autoFromList [3,2,1] == (Vector.reverse $ Vector.autoFromList [1,2,3] :: Vector Number_)

prop> forVector number_ $ \xs -> reverse (Vector.toList xs) == Vector.toList (Vector.reverse xs)
prop> forVector number_ $ \xs -> xs == Vector.reverse (Vector.reverse xs)
prop> forVector number_ $ \xs -> forVector number_ $ \ys -> Vector.reverse (xs <> ys) == Vector.reverse ys <> Vector.reverse xs
-}
reverse ::
   (Integral n, Class.Floating a) =>
   Vector (Shape.ZeroBased n) a -> Vector (Shape.ZeroBased n) a
reverse :: forall n a.
(Integral n, Floating a) =>
Vector (ZeroBased n) a -> Vector (ZeroBased n) a
reverse (Array ZeroBased n
sh ForeignPtr a
x) =
   forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize ZeroBased n
sh 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
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 a
x
      Ptr CInt
incxPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1
      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
xPtr Ptr CInt
incxPtr Ptr a
yPtr Ptr CInt
incyPtr


{- |
prop> cyclicVectorFromList [] == Vector.cyclicReverse (cyclicVectorFromList [])
prop> cyclicVectorFromList [1] == Vector.cyclicReverse (cyclicVectorFromList [1])
prop> cyclicVectorFromList [1,3,2] == Vector.cyclicReverse (cyclicVectorFromList [1,2,3])
prop> cyclicVectorFromList [1,6,5,4,3,2] == Vector.cyclicReverse (cyclicVectorFromList [1,2,3,4,5,6])

prop> forCyclicVector number_ $ \xs -> xs == Vector.cyclicReverse (Vector.cyclicReverse xs)
-}
cyclicReverse ::
   (Integral n, Class.Floating a) =>
   Vector (Shape.Cyclic n) a -> Vector (Shape.Cyclic n) a
cyclicReverse :: forall n a.
(Integral n, Floating a) =>
Vector (Cyclic n) a -> Vector (Cyclic n) a
cyclicReverse (Array Cyclic n
sh ForeignPtr a
x) =
   forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize Cyclic n
sh 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
$ forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Int
nforall a. Ord a => a -> a -> Bool
>Int
0) forall a b. (a -> b) -> a -> b
$ do
      Ptr CInt
nPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint (Int
nforall a. Num a => a -> a -> a
-Int
1)
      Ptr 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 a
x
      Ptr CInt
incxPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1
      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
$ do
         forall a. Storable a => Ptr a -> a -> IO ()
poke Ptr a
yPtr forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
=<< forall a. Storable a => Ptr a -> IO a
peek Ptr a
xPtr
         forall a.
Floating a =>
Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
Blas.copy Ptr CInt
nPtr (forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
xPtr Int
1) Ptr CInt
incxPtr (forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
yPtr Int
1) Ptr CInt
incyPtr


{- |
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 :: forall sh a.
(Indexed sh, Storable a) =>
Index sh -> Index sh -> Vector sh a -> Vector sh a
swap Index sh
i Index sh
j Vector sh a
x =
   forall a. (forall s. ST s a) -> a
runST (do
      Array (ST s) sh a
y <- forall (m :: * -> *) sh a.
(PrimMonad m, C sh, Storable a) =>
Array sh a -> m (Array m sh a)
MutArray.thaw Vector sh a
x
      a
xi <- forall (m :: * -> *) sh a.
(PrimMonad m, Indexed sh, Storable a) =>
Array m sh a -> Index sh -> m a
MutArray.read Array (ST s) sh a
y Index sh
i
      a
xj <- forall (m :: * -> *) sh a.
(PrimMonad m, Indexed sh, Storable a) =>
Array m sh a -> Index sh -> m a
MutArray.read Array (ST s) sh a
y Index sh
j
      forall (m :: * -> *) sh a.
(PrimMonad m, Indexed sh, Storable a) =>
Array m sh a -> Index sh -> a -> m ()
MutArray.write Array (ST s) sh a
y Index sh
i a
xj
      forall (m :: * -> *) sh a.
(PrimMonad m, Indexed sh, Storable a) =>
Array m sh a -> Index sh -> a -> m ()
MutArray.write Array (ST s) sh a
y Index sh
j a
xi
      forall (m :: * -> *) sh a.
(PrimMonad m, C sh, Storable a) =>
Array m sh a -> m (Array sh a)
UMutArray.unsafeFreeze Array (ST s) sh a
y)


infixl 7 -*|, .*|

newtype Dot f a = Dot {forall (f :: * -> *) a. Dot f a -> f a -> f a -> a
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
-*| :: forall sh a.
(C sh, Eq sh, Floating a) =>
Vector sh a -> Vector sh a -> a
(-*|) = forall sh a.
(C sh, Eq sh, Floating a) =>
Vector sh a -> Vector sh a -> a
dot
dot :: forall sh a.
(C sh, Eq sh, Floating a) =>
Vector sh a -> Vector sh a -> a
dot =
   forall (f :: * -> *) a. Dot f a -> f a -> f 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 :: * -> *) a. (f a -> f a -> a) -> Dot f a
Dot forall sh a.
(C sh, Eq sh, Real a) =>
Vector sh a -> Vector sh a -> a
dotReal)
      (forall (f :: * -> *) a. (f a -> f a -> a) -> Dot f a
Dot forall sh a.
(C sh, Eq sh, Real a) =>
Vector sh a -> Vector sh a -> a
dotReal)
      (forall (f :: * -> *) a. (f a -> f a -> a) -> Dot f a
Dot forall a b. (a -> b) -> a -> b
$ forall sh a.
(C sh, Eq sh, Real a) =>
Char -> Vector sh (Complex a) -> Vector sh (Complex a) -> Complex a
dotComplex Char
'T')
      (forall (f :: * -> *) a. (f a -> f a -> a) -> Dot f a
Dot forall a b. (a -> b) -> a -> b
$ forall sh a.
(C sh, Eq sh, Real a) =>
Char -> Vector sh (Complex a) -> Vector sh (Complex a) -> Complex a
dotComplex Char
'T')

{- |
prop> forVector2 number_ $ \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 :: forall sh a.
(C sh, Eq sh, Floating a) =>
Vector sh a -> Vector sh a -> a
inner =
   forall (f :: * -> *) a. Dot f a -> f a -> f 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 :: * -> *) a. (f a -> f a -> a) -> Dot f a
Dot forall sh a.
(C sh, Eq sh, Real a) =>
Vector sh a -> Vector sh a -> a
dotReal)
      (forall (f :: * -> *) a. (f a -> f a -> a) -> Dot f a
Dot forall sh a.
(C sh, Eq sh, Real a) =>
Vector sh a -> Vector sh a -> a
dotReal)
      (forall (f :: * -> *) a. (f a -> f a -> a) -> Dot f a
Dot forall a b. (a -> b) -> a -> b
$ forall sh a.
(C sh, Eq sh, Real a) =>
Char -> Vector sh (Complex a) -> Vector sh (Complex a) -> Complex a
dotComplex Char
'C')
      (forall (f :: * -> *) a. (f a -> f a -> a) -> Dot f a
Dot forall a b. (a -> b) -> a -> b
$ forall sh a.
(C sh, Eq sh, Real a) =>
Char -> Vector sh (Complex a) -> Vector sh (Complex a) -> Complex a
dotComplex Char
'C')

dotReal ::
   (Shape.C sh, Eq sh, Class.Real a) =>
   Vector sh a -> Vector sh a -> a
dotReal :: forall sh a.
(C sh, Eq sh, Real a) =>
Vector sh a -> Vector sh a -> a
dotReal arrX :: Vector sh a
arrX@(Array sh
shX ForeignPtr a
_x) (Array sh
shY ForeignPtr 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
== 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, Ptr a
sxPtr, Ptr CInt
incxPtr) <- forall sh a r.
C sh =>
Array sh a -> ContT r IO (Ptr CInt, Ptr a, Ptr CInt)
vectorArgs Vector sh a
arrX
      Ptr a
syPtr <- 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
y
      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.
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) =>
   Char -> Vector sh (Complex a) -> Vector sh (Complex a) -> Complex a
dotComplex :: forall sh a.
(C sh, Eq sh, Real a) =>
Char -> Vector sh (Complex a) -> Vector sh (Complex a) -> Complex a
dotComplex Char
trans (Array sh
shX ForeignPtr (Complex a)
x) (Array sh
shY ForeignPtr (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
== sh
shY)
   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
trans
      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 <- 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)
y
      Ptr CInt
incyPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1
      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> forVector number_ $ \xs -> Vector.sum xs == List.sum (Vector.toList xs)
-}
sum :: (Shape.C sh, Class.Floating a) => Vector sh a -> a
sum :: forall sh a. (C sh, Floating a) => Vector sh a -> a
sum (Array sh
sh ForeignPtr a
x) = forall a. IO a -> a
unsafePerformIO forall a b. (a -> b) -> a -> b
$
   forall a b. ForeignPtr a -> (Ptr a -> IO b) -> IO b
withForeignPtr ForeignPtr a
x forall a b. (a -> b) -> a -> b
$ \Ptr a
xPtr -> forall a. Floating a => Int -> Ptr a -> Int -> IO a
Private.sum (forall sh. C sh => sh -> Int
Shape.size sh
sh) Ptr a
xPtr Int
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 :: forall sh a. (C sh, Floating a) => Vector sh a -> RealOf a
absSum Vector 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 r.
C sh =>
Array sh a -> ContT r IO (Ptr CInt, Ptr a, Ptr CInt)
vectorArgs Vector 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) => Vector sh a -> RealOf a
norm2 :: forall sh a. (C sh, Floating a) => Vector sh a -> RealOf a
norm2 Vector 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 r.
C sh =>
Array sh a -> ContT r IO (Ptr CInt, Ptr a, Ptr CInt)
vectorArgs Vector 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) => Vector sh a -> RealOf a
norm2Squared :: forall sh a. (C sh, Floating a) => Vector 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. (C sh, Real a) => Vector sh a -> a
norm2SquaredReal)
      (forall (f :: * -> *) a. (f a -> RealOf a) -> Norm f a
Norm forall sh a. (C sh, Real a) => Vector sh a -> a
norm2SquaredReal)
      (forall (f :: * -> *) a. (f a -> RealOf a) -> Norm f a
Norm forall a b. (a -> b) -> a -> b
$ forall sh a. (C sh, Real a) => Vector sh a -> a
norm2SquaredReal forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a sh.
Real a =>
Vector sh (Complex a) -> Vector (sh, ComplexShape) a
decomplex)
      (forall (f :: * -> *) a. (f a -> RealOf a) -> Norm f a
Norm forall a b. (a -> b) -> a -> b
$ forall sh a. (C sh, Real a) => Vector sh a -> a
norm2SquaredReal forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a sh.
Real a =>
Vector sh (Complex a) -> Vector (sh, ComplexShape) a
decomplex)

norm2SquaredReal :: (Shape.C sh, Class.Real a) => Vector sh a -> a
norm2SquaredReal :: forall sh a. (C sh, Real a) => Vector sh a -> a
norm2SquaredReal Vector 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
$ do
      (Ptr CInt
nPtr, Ptr a
sxPtr, Ptr CInt
incxPtr) <- forall sh a r.
C sh =>
Array sh a -> ContT r IO (Ptr CInt, Ptr a, Ptr CInt)
vectorArgs Vector sh a
arr
      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


{- |
prop> forVector number_ $ \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 :: forall sh a. (C sh, Floating a) => Vector sh a -> RealOf a
normInf Vector 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.
(C sh, Floating a) =>
Vector sh a -> IO (Maybe (Int, a))
absMax Vector 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.
-}
normInf1 :: (Shape.C sh, Class.Floating a) => Vector sh a -> RealOf a
normInf1 :: forall sh a. (C sh, Floating a) => Vector sh a -> RealOf a
normInf1 Vector 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
$ do
      (Ptr CInt
nPtr, Ptr a
sxPtr, Ptr CInt
incxPtr) <- forall sh a r.
C sh =>
Array sh a -> ContT r IO (Ptr CInt, Ptr a, Ptr CInt)
vectorArgs Vector sh a
arr
      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 -> CInt -> IO (Maybe (Int, a))
peekElemOff1 Ptr a
sxPtr 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!

>>> Vector.argAbsMaximum $ Vector.autoFromList [1:+2, 3:+4, 5, 6 :: Complex_]
(3,6.0 :+ 0.0)

prop> forVector number_ $ \xs -> isNonEmpty xs ==> let (xi,xm) = Vector.argAbsMaximum xs in xs!xi == xm
prop> forVector number_ $ \xs -> isNonEmpty xs ==> let (_xi,xm) = Vector.argAbsMaximum xs in List.all (\x -> absolute x <= absolute xm) $ Vector.toList xs
prop> forVector number_ $ \xs -> forVector number_ $ \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 :: forall sh a.
(InvIndexed sh, Floating a) =>
Vector sh a -> (Index sh, a)
argAbsMaximum Vector 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 a. Array sh a -> sh
Array.shape Vector sh a
arr))) forall a b. (a -> b) -> a -> b
$
   forall sh a.
(C sh, Floating a) =>
Vector sh a -> IO (Maybe (Int, a))
absMax Vector sh a
arr

absMax ::
   (Shape.C sh, Class.Floating a) =>
   Vector sh a -> IO (Maybe (Int, a))
absMax :: forall sh a.
(C sh, Floating a) =>
Vector sh a -> IO (Maybe (Int, a))
absMax arr :: Vector sh a
arr@(Array sh
sh ForeignPtr a
x) =
   case forall a (f :: * -> *). Floating a => f a -> ComplexSingleton a
Scalar.complexSingletonOfFunctor Vector sh a
arr 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 r.
C sh =>
Array sh a -> ContT r IO (Ptr CInt, Ptr a, Ptr CInt)
vectorArgs Vector sh a
arr
         forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a. Storable a => Ptr a -> CInt -> IO (Maybe (Int, a))
peekElemOff1 Ptr a
sxPtr 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 sh
sh
         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. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a. Storable a => Ptr a -> CInt -> IO (Maybe (Int, a))
peekElemOff1 Ptr a
sxPtr forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
=<< forall a. Real a => Int -> Ptr (Complex a) -> IO CInt
absMaxComplex Int
n Ptr a
sxPtr

absMaxComplex :: (Class.Real a) => Int -> Ptr (Complex a) -> IO CInt
absMaxComplex :: forall a. Real a => Int -> Ptr (Complex a) -> IO CInt
absMaxComplex Int
n Ptr (Complex a)
sxPtr =
   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
      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
2 Ptr (RealOf (Complex a))
xrPtr Int
2 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
2 Ptr a
xiPtr Int
2 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!

>>> Vector.argAbs1Maximum $ Vector.autoFromList [1:+2, 3:+4, 5, 6 :: Complex_]
(1,3.0 :+ 4.0)

prop> forVector real_ $ \xs -> isNonEmpty xs ==> Vector.argAbsMaximum xs == Vector.argAbs1Maximum xs
-}
argAbs1Maximum ::
   (Shape.InvIndexed sh, Class.Floating a) =>
   Vector sh a -> (Shape.Index sh, a)
argAbs1Maximum :: forall sh a.
(InvIndexed sh, Floating a) =>
Vector sh a -> (Index sh, a)
argAbs1Maximum Vector 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
$ do
      (Ptr CInt
nPtr, Ptr a
sxPtr, Ptr CInt
incxPtr) <- forall sh a r.
C sh =>
Array sh a -> ContT r IO (Ptr CInt, Ptr a, Ptr CInt)
vectorArgs Vector sh a
arr
      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 a. Array sh a -> sh
Array.shape Vector sh a
arr))) forall a b. (a -> b) -> a -> b
$
         forall a. Storable a => Ptr a -> CInt -> IO (Maybe (Int, a))
peekElemOff1 Ptr a
sxPtr 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

vectorArgs ::
   (Shape.C sh) => Array sh a -> ContT r IO (Ptr CInt, Ptr a, Ptr CInt)
vectorArgs :: forall sh a r.
C sh =>
Array sh a -> ContT r IO (Ptr CInt, Ptr a, Ptr CInt)
vectorArgs (Array sh
sh ForeignPtr a
x) =
   forall (f :: * -> *) a b c d.
Applicative f =>
(a -> b -> c -> d) -> f a -> f b -> f c -> f d
liftA3 (,,)
      (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
sh)
      (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 r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1)

peekElemOff1 :: (Storable a) => Ptr a -> CInt -> IO (Maybe (Int, a))
peekElemOff1 :: forall a. Storable a => Ptr a -> CInt -> IO (Maybe (Int, a))
peekElemOff1 Ptr a
ptr 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
ki


{- |
prop> QC.forAll (QC.choose (0,10)) $ \dim -> QC.forAll (genVector (shapeInt dim) (genNumber 3)) $ \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 :: forall sh a. (C sh, Floating a) => Vector sh a -> a
product (Array sh
sh ForeignPtr a
x) = forall a. IO a -> a
unsafePerformIO forall a b. (a -> b) -> a -> b
$
   forall a b. ForeignPtr a -> (Ptr a -> IO b) -> IO b
withForeignPtr ForeignPtr a
x forall a b. (a -> b) -> a -> b
$ \Ptr a
xPtr -> forall a. Floating a => Int -> Ptr a -> Int -> IO a
Private.product (forall sh. C sh => sh -> Int
Shape.size sh
sh) Ptr a
xPtr Int
1


{- |
For restrictions see 'limits'.

prop> forVector real_ $ \xs -> isNonEmpty xs ==> Vector.minimum xs == List.minimum (Vector.toList xs)
prop> forVector real_ $ \xs -> isNonEmpty xs ==> Vector.maximum xs == List.maximum (Vector.toList xs)
prop> forVector real_ $ \xs -> isNonEmpty xs ==> - Vector.maximum xs == Vector.minimum (Vector.negate xs)
-}
minimum, maximum :: (Shape.C sh, Class.Real a) => Vector sh a -> a
minimum :: forall sh a. (C sh, Real a) => Vector sh a -> a
minimum = forall a b. (a, b) -> a
fst forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall sh a. (C sh, Real a) => Vector sh a -> (a, a)
limits
maximum :: forall sh a. (C sh, Real a) => Vector sh a -> a
maximum = forall a b. (a, b) -> b
snd forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall sh a. (C sh, Real a) => Vector sh a -> (a, a)
limits

{- |
For restrictions see 'limits'.
-}
argMinimum, argMaximum ::
   (Shape.InvIndexed sh, Shape.Index sh ~ ix, Class.Real a) =>
   Vector sh a -> (ix,a)
argMinimum :: forall sh ix a.
(InvIndexed sh, Index sh ~ ix, Real a) =>
Vector sh a -> (ix, a)
argMinimum = forall a b. (a, b) -> a
fst forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall sh ix a.
(InvIndexed sh, Index sh ~ ix, Real a) =>
Vector sh a -> ((ix, a), (ix, a))
argLimits
argMaximum :: forall sh ix a.
(InvIndexed sh, Index sh ~ ix, Real a) =>
Vector sh a -> (ix, a)
argMaximum = forall a b. (a, b) -> b
snd forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall sh ix a.
(InvIndexed sh, Index sh ~ ix, Real a) =>
Vector sh a -> ((ix, a), (ix, a))
argLimits

{- |
prop> forVector real_ $ \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 :: forall sh a. (C sh, Real a) => Vector sh a -> (a, a)
limits Vector sh a
xs0 =
   let xs :: Array (Deferred sh) a
xs = forall sh0 sh1 a. (sh0 -> sh1) -> Array sh0 a -> Array sh1 a
Array.mapShape forall sh. sh -> Deferred sh
Shape.Deferred Vector sh a
xs0
       x0 :: a
x0 = forall a b. (a, b) -> b
snd forall a b. (a -> b) -> a -> b
$ forall sh a.
(InvIndexed sh, Floating a) =>
Vector sh a -> (Index sh, a)
argAbs1Maximum Array (Deferred sh) a
xs
       x1 :: a
x1 = Array (Deferred sh) a
xs forall sh a.
(Indexed sh, Storable a) =>
Array sh a -> Index sh -> a
! forall a b. (a, b) -> a
fst (forall sh a.
(InvIndexed sh, Floating a) =>
Vector sh a -> (Index sh, a)
argAbs1Maximum (forall sh a. (C sh, Floating a) => a -> Vector sh a -> Vector sh a
raise (-a
x0) Array (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.InvIndexed sh, Shape.Index sh ~ ix, Class.Real a) =>
   Vector sh a -> ((ix,a),(ix,a))
argLimits :: forall sh ix a.
(InvIndexed sh, Index sh ~ ix, Real a) =>
Vector sh a -> ((ix, a), (ix, a))
argLimits Vector sh a
xs =
   let p0 :: (Index sh, a)
p0@(Index sh
_i0,a
x0) = forall sh a.
(InvIndexed sh, Floating a) =>
Vector sh a -> (Index sh, a)
argAbs1Maximum Vector sh a
xs
       p1 :: (ix, a)
p1 = (ix
i1,Vector sh a
xsforall sh a.
(Indexed sh, Storable a) =>
Array 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.
(InvIndexed sh, Floating a) =>
Vector sh a -> (Index sh, a)
argAbs1Maximum forall a b. (a -> b) -> a -> b
$ forall sh a. (C sh, Floating a) => a -> Vector sh a -> Vector sh a
raise (-a
x0) Vector 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> forVector number_ $ \xs -> Vector.negate xs == Vector.scale minusOne xs
prop> forVector number_ $ \xs -> Vector.scale 2 xs == xs |+| xs
-}
scale, _scale, (.*|) ::
   (Shape.C sh, Class.Floating a) =>
   a -> Vector sh a -> Vector sh a
.*| :: forall sh a. (C sh, Floating a) => a -> Vector sh a -> Vector sh a
(.*|) = forall sh a. (C sh, Floating a) => a -> Vector sh a -> Vector sh a
scale

scale :: forall sh a. (C sh, Floating a) => a -> Vector sh a -> Vector sh a
scale a
alpha (Array sh
sh ForeignPtr a
x) = forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize sh
sh 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 <- 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
      Ptr CInt
incxPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1
      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. (C sh, Floating a) => a -> Vector sh a -> Vector sh a
_scale a
a (Array sh
sh ForeignPtr a
b) = forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize sh
sh 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 <- 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
b
      Ptr CInt
ldbPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
k
      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


scaleReal ::
   (Shape.C sh, Class.Floating a) =>
   RealOf a -> Vector sh a -> Vector sh a
scaleReal :: forall sh a.
(C sh, Floating a) =>
RealOf a -> Vector sh a -> Vector sh a
scaleReal =
   forall (f :: * -> *) a. ScaleReal f a -> RealOf a -> f a -> f 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 :: * -> *) a. (RealOf a -> f a -> f a) -> ScaleReal f a
ScaleReal forall sh a. (C sh, Floating a) => a -> Vector sh a -> Vector sh a
scale)
      (forall (f :: * -> *) a. (RealOf a -> f a -> f a) -> ScaleReal f a
ScaleReal forall sh a. (C sh, Floating a) => a -> Vector sh a -> Vector sh a
scale)
      (forall (f :: * -> *) a. (RealOf a -> f a -> f a) -> ScaleReal f a
ScaleReal forall a b. (a -> b) -> a -> b
$ \RealOf (Complex Float)
x -> forall a sh.
Real a =>
Vector (sh, ComplexShape) a -> Vector sh (Complex a)
recomplex forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall sh a. (C sh, Floating a) => a -> Vector sh a -> Vector sh a
scale RealOf (Complex Float)
x forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a sh.
Real a =>
Vector sh (Complex a) -> Vector (sh, ComplexShape) a
decomplex)
      (forall (f :: * -> *) a. (RealOf a -> f a -> f a) -> ScaleReal f a
ScaleReal forall a b. (a -> b) -> a -> b
$ \RealOf (Complex Double)
x -> forall a sh.
Real a =>
Vector (sh, ComplexShape) a -> Vector sh (Complex a)
recomplex forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall sh a. (C sh, Floating a) => a -> Vector sh a -> Vector sh a
scale RealOf (Complex Double)
x forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a sh.
Real a =>
Vector sh (Complex a) -> Vector (sh, ComplexShape) a
decomplex)

newtype ScaleReal f a = ScaleReal {forall (f :: * -> *) a. ScaleReal f a -> RealOf a -> f a -> f a
getScaleReal :: RealOf a -> f a -> f a}


decomplex ::
   (Class.Real a) =>
   Vector sh (Complex a) -> Vector (sh, ComplexShape) a
decomplex :: forall a sh.
Real a =>
Vector sh (Complex a) -> Vector (sh, ComplexShape) a
decomplex (Array sh
sh ForeignPtr (Complex a)
a) = forall sh a. sh -> ForeignPtr a -> Array sh a
Array (sh
sh, forall sh. Static sh => sh
Shape.static) (forall a b. ForeignPtr a -> ForeignPtr b
castForeignPtr ForeignPtr (Complex a)
a)

recomplex ::
   (Class.Real a) =>
   Vector (sh, ComplexShape) a -> Vector sh (Complex a)
recomplex :: forall a sh.
Real a =>
Vector (sh, ComplexShape) a -> Vector sh (Complex a)
recomplex (Array (sh
sh, Shape.NestedTuple Complex Element
_) ForeignPtr a
a) = forall sh a. sh -> ForeignPtr a -> Array sh a
Array sh
sh (forall a b. ForeignPtr a -> ForeignPtr b
castForeignPtr ForeignPtr a
a)



infixl 6 |+|, |-|, `add`, `sub`


{- |
prop> forVector2 number_ $ \xs ys -> xs |+| ys == ys |+| xs
prop> forVector2 number_ $ \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 :: forall sh a.
(C sh, Eq sh, Floating a) =>
Vector sh a -> Vector sh a -> Vector sh a
add = forall sh a.
(C sh, Eq sh, Floating a) =>
a -> Vector sh a -> Vector sh a -> Vector sh a
mac forall a. Floating a => a
Scalar.one
sub :: forall sh a.
(C sh, Eq sh, Floating a) =>
Vector sh a -> Vector sh a -> Vector sh a
sub Vector sh a
x Vector sh a
y = forall sh a.
(C sh, Eq sh, Floating a) =>
a -> Vector sh a -> Vector sh a -> Vector sh a
mac forall a. Floating a => a
minusOne Vector sh a
y Vector sh a
x

|+| :: forall sh a.
(C sh, Eq sh, Floating a) =>
Vector sh a -> Vector sh a -> Vector sh a
(|+|) = forall sh a.
(C sh, Eq sh, Floating a) =>
Vector sh a -> Vector sh a -> Vector sh a
add
|-| :: forall sh a.
(C sh, Eq sh, Floating a) =>
Vector sh a -> Vector sh a -> Vector sh a
(|-|) = forall sh a.
(C sh, Eq sh, Floating a) =>
Vector sh a -> Vector sh a -> Vector sh a
sub

mac ::
   (Shape.C sh, Eq sh, Class.Floating a) =>
   a -> Vector sh a -> Vector sh a -> Vector sh a
mac :: forall sh a.
(C sh, Eq sh, Floating a) =>
a -> Vector sh a -> Vector sh a -> Vector sh a
mac a
alpha (Array sh
shX ForeignPtr a
x) (Array sh
shY ForeignPtr a
y) =
   forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize sh
shX forall a b. (a -> b) -> a -> b
$ \Int
n Ptr a
szPtr -> do
   String -> Bool -> IO ()
Call.assert String
"mac: 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 Int
n
      Ptr a
saPtr <- forall a r. Floating a => a -> FortranIO r (Ptr a)
Call.number a
alpha
      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
      Ptr CInt
incxPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1
      Ptr a
syPtr <- 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
y
      Ptr CInt
incyPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1
      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> forVector number_ $ \xs -> xs == Vector.negate (Vector.negate xs)
-}
negate :: (Shape.C sh, Class.Floating a) => Vector sh a -> Vector sh a
negate :: forall sh a. (C sh, Floating a) => Vector sh a -> Vector sh a
negate =
   forall (f :: * -> *) a. Conjugate f a -> f a -> f 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 :: * -> *) a. (f a -> f a) -> Conjugate f a
Conjugate forall a b. (a -> b) -> a -> b
$ forall sh a.
(C sh, Floating a) =>
RealOf a -> Vector sh a -> Vector sh a
scaleReal forall a. Floating a => a
Scalar.minusOne)
      (forall (f :: * -> *) a. (f a -> f a) -> Conjugate f a
Conjugate forall a b. (a -> b) -> a -> b
$ forall sh a.
(C sh, Floating a) =>
RealOf a -> Vector sh a -> Vector sh a
scaleReal forall a. Floating a => a
Scalar.minusOne)
      (forall (f :: * -> *) a. (f a -> f a) -> Conjugate f a
Conjugate forall a b. (a -> b) -> a -> b
$ forall sh a.
(C sh, Floating a) =>
RealOf a -> Vector sh a -> Vector sh a
scaleReal forall a. Floating a => a
Scalar.minusOne)
      (forall (f :: * -> *) a. (f a -> f a) -> Conjugate f a
Conjugate forall a b. (a -> b) -> a -> b
$ forall sh a.
(C sh, Floating a) =>
RealOf a -> Vector sh a -> Vector sh a
scaleReal forall a. Floating a => a
Scalar.minusOne)


{- |
prop> QC.forAll (genNumber maxElem) $ \d -> forVector number_ $ \xs -> xs == Vector.raise (-d) (Vector.raise d xs)
-}
raise :: (Shape.C sh, Class.Floating a) => a -> Vector sh a -> Vector sh a
raise :: forall sh a. (C sh, Floating a) => a -> Vector sh a -> Vector sh a
raise a
c (Array sh
shX ForeignPtr a
x) =
   forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize sh
shX 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 <- 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
      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
inc1Ptr 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> forVector2 number_ $ \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 :: forall sh a.
(C sh, Eq sh, Floating a) =>
Vector sh a -> Vector sh a -> Vector sh a
mul = forall sh a.
(C sh, Eq sh, Floating a) =>
Conjugation -> Vector sh a -> Vector sh a -> Vector sh a
mulConjugation Conjugation
NonConjugated

{- |
prop> forVector2 number_ $ \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 :: forall sh a.
(C sh, Eq sh, Floating a) =>
Vector sh a -> Vector sh a -> Vector sh a
mulConj = forall sh a.
(C sh, Eq sh, Floating a) =>
Conjugation -> Vector sh a -> Vector sh a -> Vector sh a
mulConjugation Conjugation
Conjugated

mulConjugation ::
   (Shape.C sh, Eq sh, Class.Floating a) =>
   Conjugation -> Vector sh a -> Vector sh a -> Vector sh a
mulConjugation :: forall sh a.
(C sh, Eq sh, Floating a) =>
Conjugation -> Vector sh a -> Vector sh a -> Vector sh a
mulConjugation Conjugation
conj (Array sh
shA ForeignPtr a
a) (Array sh
shX ForeignPtr a
x) =
      forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize sh
shX forall a b. (a -> b) -> a -> b
$ \Int
n Ptr a
yPtr -> do
   String -> Bool -> IO ()
Call.assert String
"mul: shapes mismatch" (sh
shA forall a. Eq a => a -> a -> Bool
== sh
shX)
   forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT forall a b. (a -> b) -> a -> b
$ do
      Ptr a
aPtr <- 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
a
      Ptr 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 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 Int
1 Ptr a
xPtr Int
1 Ptr a
yPtr Int
1


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

conjugate ::
   (Shape.C sh, Class.Floating a) =>
   Vector sh a -> Vector sh a
conjugate :: forall sh a. (C sh, Floating a) => Vector sh a -> Vector sh a
conjugate =
   forall (f :: * -> *) a. Conjugate f a -> f a -> f 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 :: * -> *) a. (f a -> f a) -> Conjugate f a
Conjugate forall a. a -> a
id)
      (forall (f :: * -> *) a. (f a -> f a) -> Conjugate f a
Conjugate forall a. a -> a
id)
      (forall (f :: * -> *) a. (f a -> f a) -> Conjugate f a
Conjugate forall sh a.
(C sh, Real a) =>
Vector sh (Complex a) -> Vector sh (Complex a)
complexConjugate)
      (forall (f :: * -> *) a. (f a -> f a) -> Conjugate f a
Conjugate forall sh a.
(C sh, Real a) =>
Vector sh (Complex a) -> Vector sh (Complex a)
complexConjugate)

complexConjugate ::
   (Shape.C sh, Class.Real a) =>
   Vector sh (Complex a) -> Vector sh (Complex a)
complexConjugate :: forall sh a.
(C sh, Real a) =>
Vector sh (Complex a) -> Vector sh (Complex a)
complexConjugate (Array sh
sh ForeignPtr (Complex a)
x) = forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize sh
sh 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 <- 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
incxPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1
      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) => Vector sh (RealOf a) -> Vector sh a
fromReal :: forall sh a.
(C sh, Floating a) =>
Vector sh (RealOf a) -> Vector sh a
fromReal =
   forall (f :: * -> *) a. FromReal f a -> f (RealOf a) -> f 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 :: * -> *) a. (f (RealOf a) -> f a) -> FromReal f a
FromReal forall a. a -> a
id)
      (forall (f :: * -> *) a. (f (RealOf a) -> f a) -> FromReal f a
FromReal forall a. a -> a
id)
      (forall (f :: * -> *) a. (f (RealOf a) -> f a) -> FromReal f a
FromReal forall sh a. (C sh, Real a) => Vector sh a -> Vector sh (Complex a)
complexFromReal)
      (forall (f :: * -> *) a. (f (RealOf a) -> f a) -> FromReal f a
FromReal forall sh a. (C sh, Real a) => Vector sh a -> Vector sh (Complex a)
complexFromReal)

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

toComplex ::
   (Shape.C sh, Class.Floating a) => Vector sh a -> Vector sh (ComplexOf a)
toComplex :: forall sh a.
(C sh, Floating a) =>
Vector sh a -> Vector sh (ComplexOf a)
toComplex =
   forall (f :: * -> *) a. ToComplex f a -> f a -> f (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 :: * -> *) a. (f a -> f (ComplexOf a)) -> ToComplex f a
ToComplex forall sh a. (C sh, Real a) => Vector sh a -> Vector sh (Complex a)
complexFromReal)
      (forall (f :: * -> *) a. (f a -> f (ComplexOf a)) -> ToComplex f a
ToComplex forall sh a. (C sh, Real a) => Vector sh a -> Vector sh (Complex a)
complexFromReal)
      (forall (f :: * -> *) a. (f a -> f (ComplexOf a)) -> ToComplex f a
ToComplex forall a. a -> a
id)
      (forall (f :: * -> *) a. (f a -> f (ComplexOf a)) -> ToComplex f a
ToComplex forall a. a -> a
id)

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

complexFromReal ::
   (Shape.C sh, Class.Real a) => Vector sh a -> Vector sh (Complex a)
complexFromReal :: forall sh a. (C sh, Real a) => Vector sh a -> Vector sh (Complex a)
complexFromReal (Array sh
sh ForeignPtr a
x) =
   forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize sh
sh 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 <- 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
         Ptr CInt
incxPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1
         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


realPart ::
   (Shape.C sh, Class.Floating a) => Vector sh a -> Vector sh (RealOf a)
realPart :: forall sh a.
(C sh, Floating a) =>
Vector sh a -> Vector sh (RealOf a)
realPart =
   forall (f :: * -> *) a. ToReal f a -> f a -> f (RealOf a)
getToReal 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 -> f (RealOf a)) -> ToReal f a
ToReal forall a. a -> a
id)
      (forall (f :: * -> *) a. (f a -> f (RealOf a)) -> ToReal f a
ToReal forall a. a -> a
id)
      (forall (f :: * -> *) a. (f a -> f (RealOf a)) -> ToReal f a
ToReal forall a b. (a -> b) -> a -> b
$ forall height width ix a.
(C height, Indexed width, Index width ~ ix, Floating a) =>
ix -> Matrix height width a -> Vector height a
RowMajor.takeColumn ElementIndex (Complex Element)
ixReal forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a sh.
Real a =>
Vector sh (Complex a) -> Vector (sh, ComplexShape) a
decomplex)
      (forall (f :: * -> *) a. (f a -> f (RealOf a)) -> ToReal f a
ToReal forall a b. (a -> b) -> a -> b
$ forall height width ix a.
(C height, Indexed width, Index width ~ ix, Floating a) =>
ix -> Matrix height width a -> Vector height a
RowMajor.takeColumn ElementIndex (Complex Element)
ixReal forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a sh.
Real a =>
Vector sh (Complex a) -> Vector (sh, ComplexShape) a
decomplex)
{-
      (ToReal $ RowMajor.takeColumn Complex.realPart . decomplex)
      (ToReal $ RowMajor.takeColumn Complex.realPart . decomplex)
-}

newtype ToReal f a = ToReal {forall (f :: * -> *) a. ToReal f a -> f a -> f (RealOf a)
getToReal :: f a -> f (RealOf a)}

imaginaryPart ::
   (Shape.C sh, Class.Real a) => Vector sh (Complex a) -> Vector sh a
imaginaryPart :: forall sh a. (C sh, Real a) => Vector sh (Complex a) -> Vector sh a
imaginaryPart = forall height width ix a.
(C height, Indexed width, Index width ~ ix, Floating a) =>
ix -> Matrix height width a -> Vector height a
RowMajor.takeColumn ElementIndex (Complex Element)
ixImaginary forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a sh.
Real a =>
Vector sh (Complex a) -> Vector (sh, ComplexShape) a
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 :: forall sh a.
(C sh, Eq sh, Real a) =>
Vector sh a -> Vector sh a -> Vector sh (Complex a)
zipComplex (Array sh
shr ForeignPtr a
xr) (Array sh
shi ForeignPtr a
xi) =
   forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize sh
shr 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" (sh
shrforall a. Eq a => a -> a -> Bool
==sh
shi)
      Ptr CInt
nPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
      Ptr a
xrPtr <- 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
xr
      Ptr a
xiPtr <- 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
xi
      let yrPtr :: Ptr (RealOf (Complex a))
yrPtr = forall a. Ptr a -> Ptr (RealOf a)
realPtr Ptr (Complex a)
yPtr
      Ptr CInt
incxPtr <- forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1
      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
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
xiPtr Ptr CInt
incxPtr (forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr (RealOf (Complex a))
yrPtr Int
1) Ptr CInt
incyPtr

unzipComplex ::
   (Shape.C sh, Class.Real a) =>
   Vector sh (Complex a) -> (Vector sh a, Vector sh a)
unzipComplex :: forall sh a.
(C sh, Real a) =>
Vector sh (Complex a) -> (Vector sh a, Vector sh a)
unzipComplex Vector sh (Complex a)
x = (forall sh a.
(C sh, Floating a) =>
Vector sh a -> Vector sh (RealOf a)
realPart Vector sh (Complex a)
x, forall sh a. (C sh, Real a) => Vector sh (Complex a) -> Vector sh a
imaginaryPart Vector sh (Complex a)
x)