{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}
module Numeric.LAPACK.Vector (
   Vector,
   RealOf,
   ComplexOf,
   Vector.toList,
   Vector.fromList,
   Vector.autoFromList,
   CheckedArray.append, (Vector.+++),
   CheckedArray.take, CheckedArray.drop,
   CheckedArray.takeLeft, CheckedArray.takeRight,
   Vector.swap,
   CheckedArray.singleton,
   Vector.constant,
   Vector.zero,
   Vector.one,
   Vector.unit,
   Vector.dot, Vector.inner, (Vector.-*|),
   Vector.sum,
   Vector.absSum,
   norm1,
   Vector.norm2,
   Vector.norm2Squared,
   normInf,
   Vector.normInf1,
   argAbsMaximum,
   Vector.argAbs1Maximum,
   Vector.product,
   Vector.scale, Vector.scaleReal, (Vector..*|),
   Vector.add, Vector.sub, (Vector.|+|), (Vector.|-|),
   Vector.negate, Vector.raise,
   Vector.mac,
   Vector.mul, Vector.mulConj,
   divide, recip,
   Vector.minimum, Vector.argMinimum,
   Vector.maximum, Vector.argMaximum,
   Vector.limits, Vector.argLimits,
   CheckedArray.foldl,
   CheckedArray.foldl1,
   CheckedArray.foldMap,

   conjugate,
   Vector.fromReal,
   Vector.toComplex,
   Vector.realPart,
   Vector.imaginaryPart,
   Vector.zipComplex,
   Vector.unzipComplex,

   random, RandomDistribution(..),
   ) where

import qualified Numeric.LAPACK.Vector.Private as VectorPriv
import qualified Numeric.LAPACK.Scalar as Scalar
import qualified Numeric.LAPACK.Private as Private
import qualified Numeric.BLAS.Vector as Vector
import Numeric.LAPACK.Linear.Private (withInfo)
import Numeric.LAPACK.Scalar (ComplexOf, RealOf, absolute)
import Numeric.LAPACK.Private (copyConjugate)

import qualified Numeric.LAPACK.FFI.Generic as LapackGen
import qualified Numeric.LAPACK.FFI.Complex as LapackComplex
import qualified Numeric.BLAS.FFI.Real as BlasReal
import qualified Numeric.Netlib.Utility as Call
import qualified Numeric.Netlib.Class as Class

import Foreign.ForeignPtr (withForeignPtr)
import Foreign.Ptr (Ptr)
import Foreign.Storable (Storable, peekElemOff, pokeElemOff)
import Foreign.C.Types (CInt)

import System.IO.Unsafe (unsafePerformIO)

import Control.Monad.Trans.Cont (ContT(ContT), evalContT)
import Control.Monad.IO.Class (liftIO)
import Control.Monad (fmap, return, (=<<))
import Control.Applicative (liftA3, (<$>))

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))

import Data.Function (id, ($), (.))
import Data.Complex (Complex)
import Data.Maybe (Maybe(Nothing,Just), maybe)
import Data.Tuple.HT (mapFst, uncurry3)
import Data.Tuple (snd)
import Data.Word (Word64)
import Data.Bits (shiftR, (.&.))
import Data.Ord (Ord)
import Data.Eq (Eq, (==))
import Data.Bool (Bool(False,True))

import Prelude (Int, fromIntegral, (+), (-), (*), Show, Enum, error, IO)


type Vector = Array


norm1 :: (Shape.C sh, Class.Floating a) => Vector sh a -> RealOf a
norm1 :: forall sh a. (C sh, Floating a) => Vector sh a -> RealOf a
norm1 Vector sh a
arr = IO (RealOf a) -> RealOf a
forall a. IO a -> a
unsafePerformIO (IO (RealOf a) -> RealOf a) -> IO (RealOf a) -> RealOf a
forall a b. (a -> b) -> a -> b
$
   ContT (RealOf a) IO (RealOf a) -> IO (RealOf a)
forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT (ContT (RealOf a) IO (RealOf a) -> IO (RealOf a))
-> ContT (RealOf a) IO (RealOf a) -> IO (RealOf a)
forall a b. (a -> b) -> a -> b
$ IO (RealOf a) -> ContT (RealOf a) IO (RealOf a)
forall a. IO a -> ContT (RealOf a) IO a
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO (RealOf a) -> ContT (RealOf a) IO (RealOf a))
-> ((Ptr CInt, Ptr a, Ptr CInt) -> IO (RealOf a))
-> (Ptr CInt, Ptr a, Ptr CInt)
-> ContT (RealOf a) IO (RealOf a)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a))
-> (Ptr CInt, Ptr a, Ptr CInt) -> IO (RealOf a)
forall a b c d. (a -> b -> c -> d) -> (a, b, c) -> d
uncurry3 Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a)
forall a.
Floating a =>
Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a)
csum1 ((Ptr CInt, Ptr a, Ptr CInt) -> ContT (RealOf a) IO (RealOf a))
-> ContT (RealOf a) IO (Ptr CInt, Ptr a, Ptr CInt)
-> ContT (RealOf a) IO (RealOf a)
forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
=<< Vector sh a -> ContT (RealOf a) IO (Ptr CInt, Ptr a, Ptr CInt)
forall sh a r.
C sh =>
Array sh a -> ContT r IO (Ptr CInt, Ptr a, Ptr CInt)
vectorArgs Vector sh a
arr

csum1 :: Class.Floating a => Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a)
csum1 :: forall a.
Floating a =>
Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a)
csum1 =
   Norm a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a)
forall a. Norm a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a)
getNorm (Norm a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a))
-> Norm a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a)
forall a b. (a -> b) -> a -> b
$
   Norm Float
-> Norm Double
-> Norm (Complex Float)
-> Norm (Complex Double)
-> Norm a
forall a (f :: * -> *).
Floating a =>
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
forall (f :: * -> *).
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
Class.switchFloating
      ((Ptr CInt -> Ptr Float -> Ptr CInt -> IO (RealOf Float))
-> Norm Float
forall a.
(Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a)) -> Norm a
Norm Ptr CInt -> Ptr Float -> Ptr CInt -> IO Float
Ptr CInt -> Ptr Float -> Ptr CInt -> IO (RealOf Float)
forall a. Real a => Ptr CInt -> Ptr a -> Ptr CInt -> IO a
BlasReal.asum)
      ((Ptr CInt -> Ptr Double -> Ptr CInt -> IO (RealOf Double))
-> Norm Double
forall a.
(Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a)) -> Norm a
Norm Ptr CInt -> Ptr Double -> Ptr CInt -> IO Double
Ptr CInt -> Ptr Double -> Ptr CInt -> IO (RealOf Double)
forall a. Real a => Ptr CInt -> Ptr a -> Ptr CInt -> IO a
BlasReal.asum)
      ((Ptr CInt
 -> Ptr (Complex Float) -> Ptr CInt -> IO (RealOf (Complex Float)))
-> Norm (Complex Float)
forall a.
(Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a)) -> Norm a
Norm Ptr CInt -> Ptr (Complex Float) -> Ptr CInt -> IO Float
Ptr CInt
-> Ptr (Complex Float) -> Ptr CInt -> IO (RealOf (Complex Float))
forall a. Real a => Ptr CInt -> Ptr (Complex a) -> Ptr CInt -> IO a
LapackComplex.sum1)
      ((Ptr CInt
 -> Ptr (Complex Double)
 -> Ptr CInt
 -> IO (RealOf (Complex Double)))
-> Norm (Complex Double)
forall a.
(Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a)) -> Norm a
Norm Ptr CInt -> Ptr (Complex Double) -> Ptr CInt -> IO Double
Ptr CInt
-> Ptr (Complex Double) -> Ptr CInt -> IO (RealOf (Complex Double))
forall a. Real a => Ptr CInt -> Ptr (Complex a) -> Ptr CInt -> IO a
LapackComplex.sum1)

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


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 = IO (RealOf a) -> RealOf a
forall a. IO a -> a
unsafePerformIO (IO (RealOf a) -> RealOf a) -> IO (RealOf a) -> RealOf a
forall a b. (a -> b) -> a -> b
$
   ContT (RealOf a) IO (RealOf a) -> IO (RealOf a)
forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT (ContT (RealOf a) IO (RealOf a) -> IO (RealOf a))
-> ContT (RealOf a) IO (RealOf a) -> IO (RealOf a)
forall a b. (a -> b) -> a -> b
$ do
      (Ptr CInt
nPtr, Ptr a
sxPtr, Ptr CInt
incxPtr) <- Vector sh a -> ContT (RealOf a) IO (Ptr CInt, Ptr a, Ptr CInt)
forall sh a r.
C sh =>
Array sh a -> ContT r IO (Ptr CInt, Ptr a, Ptr CInt)
vectorArgs Vector sh a
arr
      IO (RealOf a) -> ContT (RealOf a) IO (RealOf a)
forall a. IO a -> ContT (RealOf a) IO a
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO (RealOf a) -> ContT (RealOf a) IO (RealOf a))
-> IO (RealOf a) -> ContT (RealOf a) IO (RealOf a)
forall a b. (a -> b) -> a -> b
$
         (Maybe (Int, a) -> RealOf a)
-> IO (Maybe (Int, a)) -> IO (RealOf a)
forall a b. (a -> b) -> IO a -> IO b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (a -> RealOf a
forall a. Floating a => a -> RealOf a
absolute (a -> RealOf a)
-> (Maybe (Int, a) -> a) -> Maybe (Int, a) -> RealOf a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. a -> ((Int, a) -> a) -> Maybe (Int, a) -> a
forall b a. b -> (a -> b) -> Maybe a -> b
maybe a
forall a. Floating a => a
Scalar.zero (Int, a) -> a
forall a b. (a, b) -> b
snd) (IO (Maybe (Int, a)) -> IO (RealOf a))
-> IO (Maybe (Int, a)) -> IO (RealOf a)
forall a b. (a -> b) -> a -> b
$
         Ptr a -> CInt -> IO (Maybe (Int, a))
forall a. Storable a => Ptr a -> CInt -> IO (Maybe (Int, a))
peekElemOff1 Ptr a
sxPtr (CInt -> IO (Maybe (Int, a))) -> IO CInt -> IO (Maybe (Int, a))
forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
=<< Ptr CInt -> Ptr a -> Ptr CInt -> IO CInt
forall a. Floating a => Ptr CInt -> Ptr a -> Ptr CInt -> IO CInt
VectorPriv.absMax 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!
-}
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 = IO (Index sh, a) -> (Index sh, a)
forall a. IO a -> a
unsafePerformIO (IO (Index sh, a) -> (Index sh, a))
-> IO (Index sh, a) -> (Index sh, a)
forall a b. (a -> b) -> a -> b
$
   ContT (Index sh, a) IO (Index sh, a) -> IO (Index sh, a)
forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT (ContT (Index sh, a) IO (Index sh, a) -> IO (Index sh, a))
-> ContT (Index sh, a) IO (Index sh, a) -> IO (Index sh, a)
forall a b. (a -> b) -> a -> b
$ do
      (Ptr CInt
nPtr, Ptr a
sxPtr, Ptr CInt
incxPtr) <- Vector sh a -> ContT (Index sh, a) IO (Ptr CInt, Ptr a, Ptr CInt)
forall sh a r.
C sh =>
Array sh a -> ContT r IO (Ptr CInt, Ptr a, Ptr CInt)
vectorArgs Vector sh a
arr
      IO (Index sh, a) -> ContT (Index sh, a) IO (Index sh, a)
forall a. IO a -> ContT (Index sh, a) IO a
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO (Index sh, a) -> ContT (Index sh, a) IO (Index sh, a))
-> IO (Index sh, a) -> ContT (Index sh, a) IO (Index sh, a)
forall a b. (a -> b) -> a -> b
$
         (Maybe (Int, a) -> (Index sh, a))
-> IO (Maybe (Int, a)) -> IO (Index sh, a)
forall a b. (a -> b) -> IO a -> IO b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap
            ((Index sh, a)
-> ((Int, a) -> (Index sh, a)) -> Maybe (Int, a) -> (Index sh, a)
forall b a. b -> (a -> b) -> Maybe a -> b
maybe
               ([Char] -> (Index sh, a)
forall a. HasCallStack => [Char] -> a
error [Char]
"Vector.argAbsMaximum: empty vector")
               ((Int -> Index sh) -> (Int, a) -> (Index sh, a)
forall a c b. (a -> c) -> (a, b) -> (c, b)
mapFst (sh -> Int -> Index sh
forall sh. InvIndexed sh => sh -> Int -> Index sh
Shape.uncheckedIndexFromOffset (sh -> Int -> Index sh) -> sh -> Int -> Index sh
forall a b. (a -> b) -> a -> b
$ Vector sh a -> sh
forall sh a. Array sh a -> sh
Array.shape Vector sh a
arr))) (IO (Maybe (Int, a)) -> IO (Index sh, a))
-> IO (Maybe (Int, a)) -> IO (Index sh, a)
forall a b. (a -> b) -> a -> b
$
         Ptr a -> CInt -> IO (Maybe (Int, a))
forall a. Storable a => Ptr a -> CInt -> IO (Maybe (Int, a))
peekElemOff1 Ptr a
sxPtr (CInt -> IO (Maybe (Int, a))) -> IO CInt -> IO (Maybe (Int, a))
forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
=<< Ptr CInt -> Ptr a -> Ptr CInt -> IO CInt
forall a. Floating a => Ptr CInt -> Ptr a -> Ptr CInt -> IO CInt
VectorPriv.absMax 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) =
   (Ptr CInt -> Ptr a -> Ptr CInt -> (Ptr CInt, Ptr a, Ptr CInt))
-> ContT r IO (Ptr CInt)
-> ContT r IO (Ptr a)
-> ContT r IO (Ptr CInt)
-> ContT r IO (Ptr CInt, Ptr a, Ptr CInt)
forall (f :: * -> *) a b c d.
Applicative f =>
(a -> b -> c -> d) -> f a -> f b -> f c -> f d
liftA3 (,,)
      (Int -> ContT r IO (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint (Int -> ContT r IO (Ptr CInt)) -> Int -> ContT r IO (Ptr CInt)
forall a b. (a -> b) -> a -> b
$ sh -> Int
forall sh. C sh => sh -> Int
Shape.size sh
sh)
      (((Ptr a -> IO r) -> IO r) -> ContT r IO (Ptr a)
forall {k} (r :: k) (m :: k -> *) a.
((a -> m r) -> m r) -> ContT r m a
ContT (((Ptr a -> IO r) -> IO r) -> ContT r IO (Ptr a))
-> ((Ptr a -> IO r) -> IO r) -> ContT r IO (Ptr a)
forall a b. (a -> b) -> a -> b
$ ForeignPtr a -> (Ptr a -> IO r) -> IO r
forall a b. ForeignPtr a -> (Ptr a -> IO b) -> IO b
withForeignPtr ForeignPtr a
x)
      (Int -> ContT r IO (Ptr CInt)
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 = CInt -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral CInt
k1
       ki :: Int
ki = Int
k1iInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1
   in if Int
k1i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0
         then Maybe (Int, a) -> IO (Maybe (Int, a))
forall a. a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return Maybe (Int, a)
forall a. Maybe a
Nothing
         else (Int, a) -> Maybe (Int, a)
forall a. a -> Maybe a
Just ((Int, a) -> Maybe (Int, a))
-> (a -> (Int, a)) -> a -> Maybe (Int, a)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (,) Int
ki (a -> Maybe (Int, a)) -> IO a -> IO (Maybe (Int, a))
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Ptr a -> Int -> IO a
forall a. Storable a => Ptr a -> Int -> IO a
peekElemOff Ptr a
ptr Int
ki


divide ::
   (Shape.C sh, Eq sh, Class.Floating a) =>
   Vector sh a -> Vector sh a -> Vector sh a
divide :: forall sh a.
(C sh, Eq sh, Floating a) =>
Vector sh a -> Vector sh a -> Vector sh a
divide (Array sh
shB ForeignPtr a
b) (Array sh
shA ForeignPtr a
a) =
      sh -> (Int -> Ptr a -> IO ()) -> Array sh a
forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize sh
shB ((Int -> Ptr a -> IO ()) -> Array sh a)
-> (Int -> Ptr a -> IO ()) -> Array sh a
forall a b. (a -> b) -> a -> b
$ \Int
n Ptr a
xPtr -> do
   [Char] -> Bool -> IO ()
Call.assert [Char]
"divide: shapes mismatch" (sh
shA sh -> sh -> Bool
forall a. Eq a => a -> a -> Bool
== sh
shB)
   ContT () IO () -> IO ()
forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT (ContT () IO () -> IO ()) -> ContT () IO () -> IO ()
forall a b. (a -> b) -> a -> b
$ do
      Ptr CInt
nPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
      Ptr CInt
klPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
0
      Ptr CInt
kuPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
0
      Ptr CInt
nrhsPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1
      Ptr a
abPtr <- Int -> ForeignPtr a -> ContT () IO (Ptr a)
forall a r. Storable a => Int -> ForeignPtr a -> ContT r IO (Ptr a)
Private.copyToTemp Int
n ForeignPtr a
a
      Ptr CInt
ldabPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
1
      Ptr CInt
ipivPtr <- Int -> FortranIO () (Ptr CInt)
forall a r. Storable a => Int -> FortranIO r (Ptr a)
Call.allocaArray Int
n
      Ptr a
bPtr <- ((Ptr a -> IO ()) -> IO ()) -> ContT () IO (Ptr a)
forall {k} (r :: k) (m :: k -> *) a.
((a -> m r) -> m r) -> ContT r m a
ContT (((Ptr a -> IO ()) -> IO ()) -> ContT () IO (Ptr a))
-> ((Ptr a -> IO ()) -> IO ()) -> ContT () IO (Ptr a)
forall a b. (a -> b) -> a -> b
$ ForeignPtr a -> (Ptr a -> IO ()) -> IO ()
forall a b. ForeignPtr a -> (Ptr a -> IO b) -> IO b
withForeignPtr ForeignPtr a
b
      Ptr CInt
ldxPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
n
      IO () -> ContT () IO ()
forall a. IO a -> ContT () IO a
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> ContT () IO ()) -> IO () -> ContT () IO ()
forall a b. (a -> b) -> a -> b
$ do
         Int -> Ptr a -> Ptr a -> IO ()
forall a. Floating a => Int -> Ptr a -> Ptr a -> IO ()
Private.copyBlock Int
n Ptr a
bPtr Ptr a
xPtr
         [Char] -> (Ptr CInt -> IO ()) -> IO ()
withInfo [Char]
"gbsv" ((Ptr CInt -> IO ()) -> IO ()) -> (Ptr CInt -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$
            Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> IO ()
forall a.
Floating a =>
Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> IO ()
LapackGen.gbsv Ptr CInt
nPtr Ptr CInt
klPtr Ptr CInt
kuPtr Ptr CInt
nrhsPtr
               Ptr a
abPtr Ptr CInt
ldabPtr Ptr CInt
ipivPtr Ptr a
xPtr Ptr CInt
ldxPtr

recip :: (Shape.C sh, Class.Floating a) => Vector sh a -> Vector sh a
recip :: forall sh a. (C sh, Floating a) => Vector sh a -> Vector sh a
recip Vector sh a
x =
   Array (Unchecked sh) a -> Vector sh a
forall sh a. Array (Unchecked sh) a -> Array sh a
VectorPriv.recheck (Array (Unchecked sh) a -> Vector sh a)
-> Array (Unchecked sh) a -> Vector sh a
forall a b. (a -> b) -> a -> b
$
   Array (Unchecked sh) a
-> Array (Unchecked sh) a -> Array (Unchecked sh) a
forall sh a.
(C sh, Eq sh, Floating a) =>
Vector sh a -> Vector sh a -> Vector sh a
divide
      (Vector sh a -> Array (Unchecked sh) a
forall sh a. Array sh a -> Array (Unchecked sh) a
VectorPriv.uncheck (Vector sh a -> Array (Unchecked sh) a)
-> Vector sh a -> Array (Unchecked sh) a
forall a b. (a -> b) -> a -> b
$ sh -> Vector sh a
forall sh a. (C sh, Floating a) => sh -> Vector sh a
Vector.one (sh -> Vector sh a) -> sh -> Vector sh a
forall a b. (a -> b) -> a -> b
$ Vector sh a -> sh
forall sh a. Array sh a -> sh
Array.shape Vector sh a
x)
      (Vector sh a -> Array (Unchecked sh) a
forall sh a. Array sh a -> Array (Unchecked sh) a
VectorPriv.uncheck Vector sh a
x)


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 =
   Conjugate (Array sh) a -> Array sh a -> Array sh a
forall (f :: * -> *) a. Conjugate f a -> f a -> f a
getConjugate (Conjugate (Array sh) a -> Array sh a -> Array sh a)
-> Conjugate (Array sh) a -> Array sh a -> Array sh a
forall a b. (a -> b) -> a -> b
$
   Conjugate (Array sh) Float
-> Conjugate (Array sh) Double
-> Conjugate (Array sh) (Complex Float)
-> Conjugate (Array sh) (Complex Double)
-> Conjugate (Array sh) a
forall a (f :: * -> *).
Floating a =>
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
forall (f :: * -> *).
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
Class.switchFloating
      ((Array sh Float -> Array sh Float) -> Conjugate (Array sh) Float
forall (f :: * -> *) a. (f a -> f a) -> Conjugate f a
Conjugate Array sh Float -> Array sh Float
forall a. a -> a
id)
      ((Array sh Double -> Array sh Double) -> Conjugate (Array sh) Double
forall (f :: * -> *) a. (f a -> f a) -> Conjugate f a
Conjugate Array sh Double -> Array sh Double
forall a. a -> a
id)
      ((Array sh (Complex Float) -> Array sh (Complex Float))
-> Conjugate (Array sh) (Complex Float)
forall (f :: * -> *) a. (f a -> f a) -> Conjugate f a
Conjugate Array sh (Complex Float) -> Array sh (Complex Float)
forall sh a.
(C sh, Real a) =>
Vector sh (Complex a) -> Vector sh (Complex a)
complexConjugate)
      ((Array sh (Complex Double) -> Array sh (Complex Double))
-> Conjugate (Array sh) (Complex Double)
forall (f :: * -> *) a. (f a -> f a) -> Conjugate f a
Conjugate Array sh (Complex Double) -> Array sh (Complex Double)
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) = sh -> (Int -> Ptr (Complex a) -> IO ()) -> Array sh (Complex a)
forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize sh
sh ((Int -> Ptr (Complex a) -> IO ()) -> Array sh (Complex a))
-> (Int -> Ptr (Complex a) -> IO ()) -> Array sh (Complex a)
forall a b. (a -> b) -> a -> b
$ \Int
n Ptr (Complex a)
syPtr ->
   ContT () IO () -> IO ()
forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT (ContT () IO () -> IO ()) -> ContT () IO () -> IO ()
forall a b. (a -> b) -> a -> b
$ do
      Ptr CInt
nPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
      Ptr (Complex a)
sxPtr <- ((Ptr (Complex a) -> IO ()) -> IO ())
-> ContT () IO (Ptr (Complex a))
forall {k} (r :: k) (m :: k -> *) a.
((a -> m r) -> m r) -> ContT r m a
ContT (((Ptr (Complex a) -> IO ()) -> IO ())
 -> ContT () IO (Ptr (Complex a)))
-> ((Ptr (Complex a) -> IO ()) -> IO ())
-> ContT () IO (Ptr (Complex a))
forall a b. (a -> b) -> a -> b
$ ForeignPtr (Complex a) -> (Ptr (Complex a) -> IO ()) -> IO ()
forall a b. ForeignPtr a -> (Ptr a -> IO b) -> IO b
withForeignPtr ForeignPtr (Complex a)
x
      Ptr CInt
incxPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1
      Ptr CInt
incyPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1
      IO () -> ContT () IO ()
forall a. IO a -> ContT () IO a
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> ContT () IO ()) -> IO () -> ContT () IO ()
forall a b. (a -> b) -> a -> b
$ Ptr CInt
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr CInt
-> IO ()
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


data RandomDistribution =
     UniformBox01
   | UniformBoxPM1
   | Normal
   | UniformDisc
   | UniformCircle
   deriving (RandomDistribution -> RandomDistribution -> Bool
(RandomDistribution -> RandomDistribution -> Bool)
-> (RandomDistribution -> RandomDistribution -> Bool)
-> Eq RandomDistribution
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: RandomDistribution -> RandomDistribution -> Bool
== :: RandomDistribution -> RandomDistribution -> Bool
$c/= :: RandomDistribution -> RandomDistribution -> Bool
/= :: RandomDistribution -> RandomDistribution -> Bool
Eq, Eq RandomDistribution
Eq RandomDistribution
-> (RandomDistribution -> RandomDistribution -> Ordering)
-> (RandomDistribution -> RandomDistribution -> Bool)
-> (RandomDistribution -> RandomDistribution -> Bool)
-> (RandomDistribution -> RandomDistribution -> Bool)
-> (RandomDistribution -> RandomDistribution -> Bool)
-> (RandomDistribution -> RandomDistribution -> RandomDistribution)
-> (RandomDistribution -> RandomDistribution -> RandomDistribution)
-> Ord RandomDistribution
RandomDistribution -> RandomDistribution -> Bool
RandomDistribution -> RandomDistribution -> Ordering
RandomDistribution -> RandomDistribution -> RandomDistribution
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
$ccompare :: RandomDistribution -> RandomDistribution -> Ordering
compare :: RandomDistribution -> RandomDistribution -> Ordering
$c< :: RandomDistribution -> RandomDistribution -> Bool
< :: RandomDistribution -> RandomDistribution -> Bool
$c<= :: RandomDistribution -> RandomDistribution -> Bool
<= :: RandomDistribution -> RandomDistribution -> Bool
$c> :: RandomDistribution -> RandomDistribution -> Bool
> :: RandomDistribution -> RandomDistribution -> Bool
$c>= :: RandomDistribution -> RandomDistribution -> Bool
>= :: RandomDistribution -> RandomDistribution -> Bool
$cmax :: RandomDistribution -> RandomDistribution -> RandomDistribution
max :: RandomDistribution -> RandomDistribution -> RandomDistribution
$cmin :: RandomDistribution -> RandomDistribution -> RandomDistribution
min :: RandomDistribution -> RandomDistribution -> RandomDistribution
Ord, Int -> RandomDistribution -> ShowS
[RandomDistribution] -> ShowS
RandomDistribution -> [Char]
(Int -> RandomDistribution -> ShowS)
-> (RandomDistribution -> [Char])
-> ([RandomDistribution] -> ShowS)
-> Show RandomDistribution
forall a.
(Int -> a -> ShowS) -> (a -> [Char]) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> RandomDistribution -> ShowS
showsPrec :: Int -> RandomDistribution -> ShowS
$cshow :: RandomDistribution -> [Char]
show :: RandomDistribution -> [Char]
$cshowList :: [RandomDistribution] -> ShowS
showList :: [RandomDistribution] -> ShowS
Show, Int -> RandomDistribution
RandomDistribution -> Int
RandomDistribution -> [RandomDistribution]
RandomDistribution -> RandomDistribution
RandomDistribution -> RandomDistribution -> [RandomDistribution]
RandomDistribution
-> RandomDistribution -> RandomDistribution -> [RandomDistribution]
(RandomDistribution -> RandomDistribution)
-> (RandomDistribution -> RandomDistribution)
-> (Int -> RandomDistribution)
-> (RandomDistribution -> Int)
-> (RandomDistribution -> [RandomDistribution])
-> (RandomDistribution
    -> RandomDistribution -> [RandomDistribution])
-> (RandomDistribution
    -> RandomDistribution -> [RandomDistribution])
-> (RandomDistribution
    -> RandomDistribution
    -> RandomDistribution
    -> [RandomDistribution])
-> Enum RandomDistribution
forall a.
(a -> a)
-> (a -> a)
-> (Int -> a)
-> (a -> Int)
-> (a -> [a])
-> (a -> a -> [a])
-> (a -> a -> [a])
-> (a -> a -> a -> [a])
-> Enum a
$csucc :: RandomDistribution -> RandomDistribution
succ :: RandomDistribution -> RandomDistribution
$cpred :: RandomDistribution -> RandomDistribution
pred :: RandomDistribution -> RandomDistribution
$ctoEnum :: Int -> RandomDistribution
toEnum :: Int -> RandomDistribution
$cfromEnum :: RandomDistribution -> Int
fromEnum :: RandomDistribution -> Int
$cenumFrom :: RandomDistribution -> [RandomDistribution]
enumFrom :: RandomDistribution -> [RandomDistribution]
$cenumFromThen :: RandomDistribution -> RandomDistribution -> [RandomDistribution]
enumFromThen :: RandomDistribution -> RandomDistribution -> [RandomDistribution]
$cenumFromTo :: RandomDistribution -> RandomDistribution -> [RandomDistribution]
enumFromTo :: RandomDistribution -> RandomDistribution -> [RandomDistribution]
$cenumFromThenTo :: RandomDistribution
-> RandomDistribution -> RandomDistribution -> [RandomDistribution]
enumFromThenTo :: RandomDistribution
-> RandomDistribution -> RandomDistribution -> [RandomDistribution]
Enum)

{-
@random distribution shape seed@

Only the least significant 47 bits of @seed@ are used.
-}
random ::
   (Shape.C sh, Class.Floating a) =>
   RandomDistribution -> sh -> Word64 -> Vector sh a
random :: forall sh a.
(C sh, Floating a) =>
RandomDistribution -> sh -> Word64 -> Vector sh a
random RandomDistribution
dist sh
sh Word64
seed = sh -> (Int -> Ptr a -> IO ()) -> Array sh a
forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize sh
sh ((Int -> Ptr a -> IO ()) -> Array sh a)
-> (Int -> Ptr a -> IO ()) -> Array sh a
forall a b. (a -> b) -> a -> b
$ \Int
n Ptr a
xPtr ->
   ContT () IO () -> IO ()
forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT (ContT () IO () -> IO ()) -> ContT () IO () -> IO ()
forall a b. (a -> b) -> a -> b
$ do
      Ptr CInt
nPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
      Ptr CInt
distPtr <-
         Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint (Int -> FortranIO () (Ptr CInt)) -> Int -> FortranIO () (Ptr CInt)
forall a b. (a -> b) -> a -> b
$
         case (Ptr a -> Bool -> Bool -> Bool
forall a (f :: * -> *) b. Floating a => f a -> b -> b -> b
Private.caseRealComplexFunc Ptr a
xPtr Bool
False Bool
True, RandomDistribution
dist) of
            (Bool
_, RandomDistribution
UniformBox01) -> Int
1
            (Bool
_, RandomDistribution
UniformBoxPM1) -> Int
2
            (Bool
_, RandomDistribution
Normal) -> Int
3
            (Bool
True, RandomDistribution
UniformDisc) -> Int
4
            (Bool
True, RandomDistribution
UniformCircle) -> Int
5
            (Bool
False, RandomDistribution
UniformDisc) -> Int
2
            (Bool
False, RandomDistribution
UniformCircle) ->
               [Char] -> Int
forall a. HasCallStack => [Char] -> a
error
                  [Char]
"Vector.random: UniformCircle not supported for real numbers"
      Ptr CInt
iseedPtr <- Int -> FortranIO () (Ptr CInt)
forall a r. Storable a => Int -> FortranIO r (Ptr a)
Call.allocaArray Int
4
      IO () -> ContT () IO ()
forall a. IO a -> ContT () IO a
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> ContT () IO ()) -> IO () -> ContT () IO ()
forall a b. (a -> b) -> a -> b
$ do
         Ptr CInt -> Int -> CInt -> IO ()
forall a. Storable a => Ptr a -> Int -> a -> IO ()
pokeElemOff Ptr CInt
iseedPtr Int
0 (CInt -> IO ()) -> CInt -> IO ()
forall a b. (a -> b) -> a -> b
$ Word64 -> CInt
forall a b. (Integral a, Num b) => a -> b
fromIntegral ((Word64
seed Word64 -> Int -> Word64
forall a. Bits a => a -> Int -> a
`shiftR` Int
35) Word64 -> Word64 -> Word64
forall a. Bits a => a -> a -> a
.&. Word64
0xFFF)
         Ptr CInt -> Int -> CInt -> IO ()
forall a. Storable a => Ptr a -> Int -> a -> IO ()
pokeElemOff Ptr CInt
iseedPtr Int
1 (CInt -> IO ()) -> CInt -> IO ()
forall a b. (a -> b) -> a -> b
$ Word64 -> CInt
forall a b. (Integral a, Num b) => a -> b
fromIntegral ((Word64
seed Word64 -> Int -> Word64
forall a. Bits a => a -> Int -> a
`shiftR` Int
23) Word64 -> Word64 -> Word64
forall a. Bits a => a -> a -> a
.&. Word64
0xFFF)
         Ptr CInt -> Int -> CInt -> IO ()
forall a. Storable a => Ptr a -> Int -> a -> IO ()
pokeElemOff Ptr CInt
iseedPtr Int
2 (CInt -> IO ()) -> CInt -> IO ()
forall a b. (a -> b) -> a -> b
$ Word64 -> CInt
forall a b. (Integral a, Num b) => a -> b
fromIntegral ((Word64
seed Word64 -> Int -> Word64
forall a. Bits a => a -> Int -> a
`shiftR` Int
11) Word64 -> Word64 -> Word64
forall a. Bits a => a -> a -> a
.&. Word64
0xFFF)
         Ptr CInt -> Int -> CInt -> IO ()
forall a. Storable a => Ptr a -> Int -> a -> IO ()
pokeElemOff Ptr CInt
iseedPtr Int
3 (CInt -> IO ()) -> CInt -> IO ()
forall a b. (a -> b) -> a -> b
$ Word64 -> CInt
forall a b. (Integral a, Num b) => a -> b
fromIntegral ((Word64
seedWord64 -> Word64 -> Word64
forall a. Bits a => a -> a -> a
.&.Word64
0x7FF)Word64 -> Word64 -> Word64
forall a. Num a => a -> a -> a
*Word64
2Word64 -> Word64 -> Word64
forall a. Num a => a -> a -> a
+Word64
1)
         Ptr CInt -> Ptr CInt -> Ptr CInt -> Ptr a -> IO ()
forall a.
Floating a =>
Ptr CInt -> Ptr CInt -> Ptr CInt -> Ptr a -> IO ()
LapackGen.larnv Ptr CInt
distPtr Ptr CInt
iseedPtr Ptr CInt
nPtr Ptr a
xPtr