{-# LANGUAGE TypeOperators #-}
module Numeric.LAPACK.Matrix.Square.Basic (
   Square,
   LiberalSquare,
   mapSize,
   toFull,
   fromFull,
   fromFullUnchecked,
   liberalFromFull,
   fromScalar,
   toScalar,
   fromList,
   autoFromList,

   transpose,
   adjoint,

   identity,
   identityOrder,
   identityFrom,
   identityFromWidth,
   identityFromHeight,
   diagonal,
   diagonalOrder,
   takeDiagonal,
   trace,

   stack,

   square,
   power,
   congruence,
   congruenceAdjoint,
   ) where


import qualified Numeric.LAPACK.Matrix.Layout.Private as Layout
import qualified Numeric.LAPACK.Matrix.Extent.Private as ExtentPriv
import qualified Numeric.LAPACK.Matrix.Extent as Extent
import qualified Numeric.LAPACK.Matrix.Basic as Basic
import qualified Numeric.LAPACK.Matrix.Private as Matrix
import qualified Numeric.LAPACK.Vector as Vector
import qualified Numeric.LAPACK.Private as Private
import Numeric.LAPACK.Matrix.Layout.Private
         (Order(RowMajor, ColumnMajor), swapOnRowMajor)
import Numeric.LAPACK.Matrix.Private
         (General, argGeneral, LiberalSquare, Square, argSquare,
          Full, mapExtent, ShapeInt, shapeInt)
import Numeric.LAPACK.Vector (Vector)
import Numeric.LAPACK.Scalar (zero, one)
import Numeric.LAPACK.Shape.Private (Unchecked(Unchecked))
import Numeric.LAPACK.Private (pokeCInt)

import qualified Numeric.LAPACK.FFI.Generic as LapackGen
import qualified Numeric.BLAS.FFI.Generic as BlasGen
import qualified Numeric.Netlib.Utility as Call
import qualified Numeric.Netlib.Class as Class

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.Array.Comfort.Shape ((::+))

import Foreign.ForeignPtr (withForeignPtr)
import Foreign.Storable (Storable, peek, poke)

import System.IO.Unsafe (unsafePerformIO)

import Control.Monad.Trans.Cont (ContT(ContT), evalContT)
import Control.Monad.IO.Class (liftIO)

import Data.Function.HT (powerAssociative)


mapSize :: (sh0 -> sh1) -> Square sh0 a -> Square sh1 a
mapSize :: (sh0 -> sh1) -> Square sh0 a -> Square sh1 a
mapSize = (Extent Shape Small Small sh0 sh0
 -> Extent Shape Small Small sh1 sh1)
-> Square sh0 a -> Square sh1 a
forall measA vertA horizA heightA widthA measB vertB horizB heightB
       widthB a.
(Extent measA vertA horizA heightA widthA
 -> Extent measB vertB horizB heightB widthB)
-> Full measA vertA horizA heightA widthA a
-> Full measB vertB horizB heightB widthB a
Basic.mapExtent ((Extent Shape Small Small sh0 sh0
  -> Extent Shape Small Small sh1 sh1)
 -> Square sh0 a -> Square sh1 a)
-> ((sh0 -> sh1)
    -> Extent Shape Small Small sh0 sh0
    -> Extent Shape Small Small sh1 sh1)
-> (sh0 -> sh1)
-> Square sh0 a
-> Square sh1 a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (sh0 -> sh1)
-> Extent Shape Small Small sh0 sh0
-> Extent Shape Small Small sh1 sh1
forall shA shB. (shA -> shB) -> Square shA -> Square shB
ExtentPriv.mapSquareSize

toFull ::
   (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
   Square sh a -> Full meas vert horiz sh sh a
toFull :: Square sh a -> Full meas vert horiz sh sh a
toFull = Map Shape Small Small meas vert horiz sh sh
-> Square sh a -> Full meas vert horiz sh sh a
forall measA vertA horizA measB vertB horizB height width a.
(Measure measA, C vertA, C horizA, Measure measB, C vertB,
 C horizB) =>
Map measA vertA horizA measB vertB horizB height width
-> Full measA vertA horizA height width a
-> Full measB vertB horizB height width a
mapExtent Map Shape Small Small meas vert horiz sh sh
forall meas vert horiz size.
(Measure meas, C vert, C horiz) =>
Square size -> Extent meas vert horiz size size
ExtentPriv.fromSquare

fromFull ::
   (Extent.Measure meas, Extent.C vert, Extent.C horiz, Eq sh) =>
   Full meas vert horiz sh sh a -> Square sh a
fromFull :: Full meas vert horiz sh sh a -> Square sh a
fromFull = Map meas vert horiz Shape Small Small sh sh
-> Full meas vert horiz sh sh a -> Square sh a
forall measA vertA horizA measB vertB horizB height width a.
(Measure measA, C vertA, C horizA, Measure measB, C vertB,
 C horizB) =>
Map measA vertA horizA measB vertB horizB height width
-> Full measA vertA horizA height width a
-> Full measB vertB horizB height width a
mapExtent Map meas vert horiz Shape Small Small sh sh
forall meas vert horiz size.
(Measure meas, C vert, C horiz, Eq size) =>
Extent meas vert horiz size size -> Square size
ExtentPriv.squareFromFull

fromFullUnchecked ::
   (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
   Full meas vert horiz sh sh a -> Square sh a
fromFullUnchecked :: Full meas vert horiz sh sh a -> Square sh a
fromFullUnchecked =
   Full Shape Small Small (Unchecked sh) (Unchecked sh) a
-> Square sh a
forall meas vert horiz height width a.
(Measure meas, C vert, C horiz) =>
Full meas vert horiz (Unchecked height) (Unchecked width) a
-> Full meas vert horiz height width a
Basic.recheck (Full Shape Small Small (Unchecked sh) (Unchecked sh) a
 -> Square sh a)
-> (Full meas vert horiz sh sh a
    -> Full Shape Small Small (Unchecked sh) (Unchecked sh) a)
-> Full meas vert horiz sh sh a
-> Square sh a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Map meas vert horiz Shape Small Small (Unchecked sh) (Unchecked sh)
-> Full meas vert horiz (Unchecked sh) (Unchecked sh) a
-> Full Shape Small Small (Unchecked sh) (Unchecked sh) a
forall measA vertA horizA measB vertB horizB height width a.
(Measure measA, C vertA, C horizA, Measure measB, C vertB,
 C horizB) =>
Map measA vertA horizA measB vertB horizB height width
-> Full measA vertA horizA height width a
-> Full measB vertB horizB height width a
mapExtent Map meas vert horiz Shape Small Small (Unchecked sh) (Unchecked sh)
forall meas vert horiz size.
(Measure meas, C vert, C horiz, Eq size) =>
Extent meas vert horiz size size -> Square size
ExtentPriv.squareFromFull (Full meas vert horiz (Unchecked sh) (Unchecked sh) a
 -> Full Shape Small Small (Unchecked sh) (Unchecked sh) a)
-> (Full meas vert horiz sh sh a
    -> Full meas vert horiz (Unchecked sh) (Unchecked sh) a)
-> Full meas vert horiz sh sh a
-> Full Shape Small Small (Unchecked sh) (Unchecked sh) a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Full meas vert horiz sh sh a
-> Full meas vert horiz (Unchecked sh) (Unchecked sh) a
forall meas vert horiz height width a.
(Measure meas, C vert, C horiz) =>
Full meas vert horiz height width a
-> Full meas vert horiz (Unchecked height) (Unchecked width) a
Basic.uncheck

liberalFromFull ::
   (Extent.Measure meas, Extent.C vert, Extent.C horiz,
    Shape.C height, Shape.C width) =>
   Full meas vert horiz height width a -> LiberalSquare height width a
liberalFromFull :: Full meas vert horiz height width a -> LiberalSquare height width a
liberalFromFull = Map meas vert horiz Size Small Small height width
-> Full meas vert horiz height width a
-> LiberalSquare height width a
forall measA vertA horizA measB vertB horizB height width a.
(Measure measA, C vertA, C horizA, Measure measB, C vertB,
 C horizB) =>
Map measA vertA horizA measB vertB horizB height width
-> Full measA vertA horizA height width a
-> Full measB vertB horizB height width a
mapExtent Map meas vert horiz Size Small Small height width
forall meas vert horiz height width.
(Measure meas, C vert, C horiz, C height, C width) =>
Extent meas vert horiz height width -> LiberalSquare height width
ExtentPriv.liberalSquareFromFull


fromScalar :: (Storable a) => a -> Square () a
fromScalar :: a -> Square () a
fromScalar a
a =
   Square () -> (Ptr a -> IO ()) -> Square () a
forall sh a.
(C sh, Storable a) =>
sh -> (Ptr a -> IO ()) -> Array sh a
Array.unsafeCreate (Order -> () -> Square ()
forall sh. Order -> sh -> Square sh
Layout.square Order
RowMajor ()) ((Ptr a -> IO ()) -> Square () a)
-> (Ptr a -> IO ()) -> Square () a
forall a b. (a -> b) -> a -> b
$ (Ptr a -> a -> IO ()) -> a -> Ptr a -> IO ()
forall a b c. (a -> b -> c) -> b -> a -> c
flip Ptr a -> a -> IO ()
forall a. Storable a => Ptr a -> a -> IO ()
poke a
a

toScalar :: (Storable a) => Square () a -> a
toScalar :: Square () a -> a
toScalar = (Order -> () -> ForeignPtr a -> a) -> Square () a -> a
forall sh a b.
(Order -> sh -> ForeignPtr a -> b) -> Square sh a -> b
argSquare ((Order -> () -> ForeignPtr a -> a) -> Square () a -> a)
-> (Order -> () -> ForeignPtr a -> a) -> Square () a -> a
forall a b. (a -> b) -> a -> b
$ \Order
_ () ForeignPtr a
a ->
   IO a -> a
forall a. IO a -> a
unsafePerformIO (IO a -> a) -> IO a -> a
forall a b. (a -> b) -> a -> b
$ ForeignPtr a -> (Ptr a -> IO a) -> IO a
forall a b. ForeignPtr a -> (Ptr a -> IO b) -> IO b
withForeignPtr ForeignPtr a
a Ptr a -> IO a
forall a. Storable a => Ptr a -> IO a
peek

fromList :: (Shape.C sh, Storable a) => sh -> [a] -> Square sh a
fromList :: sh -> [a] -> Square sh a
fromList sh
sh =
   Square sh -> [a] -> Square sh a
forall sh a. (C sh, Storable a) => sh -> [a] -> Array sh a
CheckedArray.fromList (Order -> sh -> Square sh
forall sh. Order -> sh -> Square sh
Layout.square Order
RowMajor sh
sh)

autoFromList :: (Storable a) => [a] -> Square ShapeInt a
autoFromList :: [a] -> Square ShapeInt a
autoFromList [a]
xs =
   let n :: Int
n = [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
xs
       m :: Int
m = Double -> Int
forall a b. (RealFrac a, Integral b) => a -> b
round (Double -> Int) -> Double -> Int
forall a b. (a -> b) -> a -> b
$ Double -> Double
forall a. Floating a => a -> a
sqrt (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n :: Double)
   in if Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
mInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
m
        then ShapeInt -> [a] -> Square ShapeInt a
forall sh a. (C sh, Storable a) => sh -> [a] -> Square sh a
fromList (Int -> ShapeInt
shapeInt Int
m) [a]
xs
        else [Char] -> Square ShapeInt a
forall a. HasCallStack => [Char] -> a
error [Char]
"Square.autoFromList: no quadratic number of elements"


transpose :: Square sh a -> Square sh a
transpose :: Square sh a -> Square sh a
transpose = (Full Shape Small Small sh sh -> Full Shape Small Small sh sh)
-> Square sh a -> Square sh a
forall sh0 sh1 a. (sh0 -> sh1) -> Array sh0 a -> Array sh1 a
Array.mapShape Full Shape Small Small sh sh -> Full Shape Small Small sh sh
forall meas vert horiz height width.
(Measure meas, C vert, C horiz) =>
Full meas vert horiz height width
-> Full meas horiz vert width height
Layout.transpose

adjoint :: (Shape.C sh, Class.Floating a) => Square sh a -> Square sh a
adjoint :: Square sh a -> Square sh a
adjoint = Square sh a -> Square sh a
forall sh a. Square sh a -> Square sh a
transpose (Square sh a -> Square sh a)
-> (Square sh a -> Square sh a) -> Square sh a -> Square sh a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Square sh a -> Square sh a
forall sh a. (C sh, Floating a) => Vector sh a -> Vector sh a
Vector.conjugate


identity :: (Shape.C sh, Class.Floating a) => sh -> Square sh a
identity :: sh -> Square sh a
identity = Order -> sh -> Square sh a
forall sh a. (C sh, Floating a) => Order -> sh -> Square sh a
identityOrder Order
RowMajor

identityFrom :: (Shape.C sh, Class.Floating a) => Square sh a -> Square sh a
identityFrom :: Square sh a -> Square sh a
identityFrom = (Order -> sh -> ForeignPtr a -> Square sh a)
-> Square sh a -> Square sh a
forall sh a b.
(Order -> sh -> ForeignPtr a -> b) -> Square sh a -> b
argSquare ((Order -> sh -> ForeignPtr a -> Square sh a)
 -> Square sh a -> Square sh a)
-> (Order -> sh -> ForeignPtr a -> Square sh a)
-> Square sh a
-> Square sh a
forall a b. (a -> b) -> a -> b
$ \Order
order sh
sh ForeignPtr a
_ -> Order -> sh -> Square sh a
forall sh a. (C sh, Floating a) => Order -> sh -> Square sh a
identityOrder Order
order sh
sh

identityFromWidth ::
   (Shape.C height, Shape.C width, Class.Floating a) =>
   General height width a -> Square width a
identityFromWidth :: General height width a -> Square width a
identityFromWidth =
   (Order -> height -> width -> ForeignPtr a -> Square width a)
-> General height width a -> Square width a
forall height width a b.
(Order -> height -> width -> ForeignPtr a -> b)
-> General height width a -> b
argGeneral ((Order -> height -> width -> ForeignPtr a -> Square width a)
 -> General height width a -> Square width a)
-> (Order -> height -> width -> ForeignPtr a -> Square width a)
-> General height width a
-> Square width a
forall a b. (a -> b) -> a -> b
$ \Order
order height
_ width
width ForeignPtr a
_ -> Order -> width -> Square width a
forall sh a. (C sh, Floating a) => Order -> sh -> Square sh a
identityOrder Order
order width
width

identityFromHeight ::
   (Shape.C height, Shape.C width, Class.Floating a) =>
   General height width a -> Square height a
identityFromHeight :: General height width a -> Square height a
identityFromHeight =
   (Order -> height -> width -> ForeignPtr a -> Square height a)
-> General height width a -> Square height a
forall height width a b.
(Order -> height -> width -> ForeignPtr a -> b)
-> General height width a -> b
argGeneral ((Order -> height -> width -> ForeignPtr a -> Square height a)
 -> General height width a -> Square height a)
-> (Order -> height -> width -> ForeignPtr a -> Square height a)
-> General height width a
-> Square height a
forall a b. (a -> b) -> a -> b
$ \Order
order height
height width
_ ForeignPtr a
_ -> Order -> height -> Square height a
forall sh a. (C sh, Floating a) => Order -> sh -> Square sh a
identityOrder Order
order height
height

identityOrder, _identityOrder ::
   (Shape.C sh, Class.Floating a) => Order -> sh -> Square sh a
identityOrder :: Order -> sh -> Square sh a
identityOrder Order
order sh
sh =
   Square sh -> (Ptr a -> IO ()) -> Square sh a
forall sh a.
(C sh, Storable a) =>
sh -> (Ptr a -> IO ()) -> Array sh a
Array.unsafeCreate (Order -> sh -> Square sh
forall sh. Order -> sh -> Square sh
Layout.square Order
order sh
sh) ((Ptr a -> IO ()) -> Square sh a)
-> (Ptr a -> IO ()) -> Square sh a
forall a b. (a -> b) -> a -> b
$ \Ptr a
aPtr ->
   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 CChar
uploPtr <- Char -> FortranIO () (Ptr CChar)
forall r. Char -> FortranIO r (Ptr CChar)
Call.char Char
'A'
      Ptr CInt
nPtr <- 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
$ sh -> Int
forall sh. C sh => sh -> Int
Shape.size sh
sh
      Ptr a
alphaPtr <- a -> FortranIO () (Ptr a)
forall a r. Floating a => a -> FortranIO r (Ptr a)
Call.number a
forall a. Floating a => a
zero
      Ptr a
betaPtr <- a -> FortranIO () (Ptr a)
forall a r. Floating a => a -> FortranIO r (Ptr a)
Call.number a
forall a. Floating a => a
one
      IO () -> ContT () IO ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> ContT () IO ()) -> IO () -> ContT () IO ()
forall a b. (a -> b) -> a -> b
$ Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr a
-> Ptr CInt
-> IO ()
forall a.
Floating a =>
Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr a
-> Ptr CInt
-> IO ()
LapackGen.laset Ptr CChar
uploPtr Ptr CInt
nPtr Ptr CInt
nPtr Ptr a
alphaPtr Ptr a
betaPtr Ptr a
aPtr Ptr CInt
nPtr

_identityOrder :: Order -> sh -> Square sh a
_identityOrder Order
order sh
sh =
   Square sh -> (Int -> Ptr a -> IO ()) -> Square sh a
forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize (Order -> sh -> Square sh
forall sh. Order -> sh -> Square sh
Layout.square Order
order sh
sh) ((Int -> Ptr a -> IO ()) -> Square sh a)
-> (Int -> Ptr a -> IO ()) -> Square sh a
forall a b. (a -> b) -> a -> b
$ \Int
blockSize Ptr a
yPtr ->
   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 <- FortranIO () (Ptr CInt)
forall a r. Storable a => FortranIO r (Ptr a)
Call.alloca
      Ptr a
xPtr <- a -> FortranIO () (Ptr a)
forall a r. Floating a => a -> FortranIO r (Ptr a)
Call.number a
forall a. Floating a => a
zero
      Ptr CInt
incxPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
0
      Ptr CInt
incyPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1
      IO () -> ContT () IO ()
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 -> IO ()
pokeCInt Ptr CInt
nPtr Int
blockSize
         Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
forall a.
Floating a =>
Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
BlasGen.copy Ptr CInt
nPtr Ptr a
xPtr Ptr CInt
incxPtr Ptr a
yPtr Ptr CInt
incyPtr
         let n :: CInt
n = Int -> CInt
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> CInt) -> Int -> CInt
forall a b. (a -> b) -> a -> b
$ sh -> Int
forall sh. C sh => sh -> Int
Shape.size sh
sh
         Ptr CInt -> CInt -> IO ()
forall a. Storable a => Ptr a -> a -> IO ()
poke Ptr CInt
nPtr CInt
n
         Ptr a -> a -> IO ()
forall a. Storable a => Ptr a -> a -> IO ()
poke Ptr a
xPtr a
forall a. Floating a => a
one
         Ptr CInt -> CInt -> IO ()
forall a. Storable a => Ptr a -> a -> IO ()
poke Ptr CInt
incyPtr (CInt
nCInt -> CInt -> CInt
forall a. Num a => a -> a -> a
+CInt
1)
         Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
forall a.
Floating a =>
Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
BlasGen.copy Ptr CInt
nPtr Ptr a
xPtr Ptr CInt
incxPtr Ptr a
yPtr Ptr CInt
incyPtr

diagonal :: (Shape.C sh, Class.Floating a) => Vector sh a -> Square sh a
diagonal :: Vector sh a -> Square sh a
diagonal (Array sh
sh ForeignPtr a
x) =
   Square sh -> (Int -> Ptr a -> IO ()) -> Square sh a
forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize (Order -> sh -> Square sh
forall sh. Order -> sh -> Square sh
Layout.square Order
ColumnMajor sh
sh) ((Int -> Ptr a -> IO ()) -> Square sh a)
-> (Int -> Ptr a -> IO ()) -> Square sh a
forall a b. (a -> b) -> a -> b
$
      \Int
blockSize Ptr a
yPtr ->
   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 <- FortranIO () (Ptr CInt)
forall a r. Storable a => FortranIO r (Ptr a)
Call.alloca
      Ptr a
xPtr <- ((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
x
      Ptr a
zPtr <- a -> ContT () IO (Ptr a)
forall a r. Floating a => a -> FortranIO r (Ptr a)
Call.number a
forall a. Floating a => a
zero
      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
      Ptr CInt
inczPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
0
      IO () -> ContT () IO ()
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 -> IO ()
pokeCInt Ptr CInt
nPtr Int
blockSize
         Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
forall a.
Floating a =>
Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
BlasGen.copy Ptr CInt
nPtr Ptr a
zPtr Ptr CInt
inczPtr Ptr a
yPtr Ptr CInt
incyPtr
         let n :: CInt
n = Int -> CInt
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> CInt) -> Int -> CInt
forall a b. (a -> b) -> a -> b
$ sh -> Int
forall sh. C sh => sh -> Int
Shape.size sh
sh
         Ptr CInt -> CInt -> IO ()
forall a. Storable a => Ptr a -> a -> IO ()
poke Ptr CInt
nPtr CInt
n
         Ptr CInt -> CInt -> IO ()
forall a. Storable a => Ptr a -> a -> IO ()
poke Ptr CInt
incyPtr (CInt
nCInt -> CInt -> CInt
forall a. Num a => a -> a -> a
+CInt
1)
         Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
forall a.
Floating a =>
Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
BlasGen.copy Ptr CInt
nPtr Ptr a
xPtr Ptr CInt
incxPtr Ptr a
yPtr Ptr CInt
incyPtr

diagonalOrder ::
   (Shape.C sh, Class.Floating a) => Order -> Vector sh a -> Square sh a
diagonalOrder :: Order -> Vector sh a -> Square sh a
diagonalOrder Order
order =
   (case Order
order of Order
ColumnMajor -> Square sh a -> Square sh a
forall a. a -> a
id; Order
RowMajor -> Square sh a -> Square sh a
forall sh a. Square sh a -> Square sh a
transpose) (Square sh a -> Square sh a)
-> (Vector sh a -> Square sh a) -> Vector sh a -> Square sh a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector sh a -> Square sh a
forall sh a. (C sh, Floating a) => Vector sh a -> Square sh a
diagonal


takeDiagonal :: (Shape.C sh, Class.Floating a) => Square sh a -> Vector sh a
takeDiagonal :: Square sh a -> Vector sh a
takeDiagonal = (Order -> sh -> ForeignPtr a -> Vector sh a)
-> Square sh a -> Vector sh a
forall sh a b.
(Order -> sh -> ForeignPtr a -> b) -> Square sh a -> b
argSquare ((Order -> sh -> ForeignPtr a -> Vector sh a)
 -> Square sh a -> Vector sh a)
-> (Order -> sh -> ForeignPtr a -> Vector sh a)
-> Square sh a
-> Vector sh a
forall a b. (a -> b) -> a -> b
$ \Order
_ sh
sh ForeignPtr a
x ->
   sh -> (Int -> Ptr a -> IO ()) -> Vector 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 ()) -> Vector sh a)
-> (Int -> Ptr a -> IO ()) -> Vector sh a
forall a b. (a -> b) -> a -> b
$ \Int
n Ptr a
yPtr -> 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 a
xPtr <- ((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
x
      Ptr CInt
incxPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)
      Ptr CInt
incyPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1
      IO () -> ContT () IO ()
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 a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
forall a.
Floating a =>
Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
BlasGen.copy Ptr CInt
nPtr Ptr a
xPtr Ptr CInt
incxPtr Ptr a
yPtr Ptr CInt
incyPtr

trace :: (Shape.C sh, Class.Floating a) => Square sh a -> a
trace :: Square sh a -> a
trace = (Order -> sh -> ForeignPtr a -> a) -> Square sh a -> a
forall sh a b.
(Order -> sh -> ForeignPtr a -> b) -> Square sh a -> b
argSquare ((Order -> sh -> ForeignPtr a -> a) -> Square sh a -> a)
-> (Order -> sh -> ForeignPtr a -> a) -> Square sh a -> a
forall a b. (a -> b) -> a -> b
$ \Order
_ sh
sh ForeignPtr a
x -> IO a -> a
forall a. IO a -> a
unsafePerformIO (IO a -> a) -> IO a -> a
forall a b. (a -> b) -> a -> b
$ do
   let n :: Int
n = sh -> Int
forall sh. C sh => sh -> Int
Shape.size sh
sh
   ForeignPtr a -> (Ptr a -> IO a) -> IO a
forall a b. ForeignPtr a -> (Ptr a -> IO b) -> IO b
withForeignPtr ForeignPtr a
x ((Ptr a -> IO a) -> IO a) -> (Ptr a -> IO a) -> IO a
forall a b. (a -> b) -> a -> b
$ \Ptr a
xPtr -> Int -> Ptr a -> Int -> IO a
forall a. Floating a => Int -> Ptr a -> Int -> IO a
Private.sum Int
n Ptr a
xPtr (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)


stack ::
   (Extent.Measure meas, Extent.C vert, Extent.C horiz,
    Shape.C sizeA, Eq sizeA, Shape.C sizeB, Eq sizeB, Class.Floating a) =>
   Square sizeA a -> Full meas vert horiz sizeA sizeB a ->
   Full meas horiz vert sizeB sizeA a -> Square sizeB a ->
   Square (sizeA::+sizeB) a
stack :: Square sizeA a
-> Full meas vert horiz sizeA sizeB a
-> Full meas horiz vert sizeB sizeA a
-> Square sizeB a
-> Square (sizeA ::+ sizeB) a
stack Square sizeA a
a Full meas vert horiz sizeA sizeB a
b Full meas horiz vert sizeB sizeA a
c Square sizeB a
d = Square sizeA a
-> General sizeA sizeB a
-> General sizeB sizeA a
-> Square sizeB a
-> Square (sizeA ::+ sizeB) a
forall meas vert horiz heightA heightB widthA widthB a.
(Measure meas, C vert, C horiz, C heightA, Eq heightA, C heightB,
 Eq heightB, C widthA, Eq widthA, C widthB, Eq widthB,
 Floating a) =>
Full meas vert horiz heightA widthA a
-> General heightA widthB a
-> General heightB widthA a
-> Full meas vert horiz heightB widthB a
-> Full meas vert horiz (heightA ::+ heightB) (widthA ::+ widthB) a
Basic.stack Square sizeA a
a (Full meas vert horiz sizeA sizeB a -> General sizeA sizeB a
forall meas vert horiz height width a.
(Measure meas, C vert, C horiz) =>
Full meas vert horiz height width a -> General height width a
Matrix.fromFull Full meas vert horiz sizeA sizeB a
b) (Full meas horiz vert sizeB sizeA a -> General sizeB sizeA a
forall meas vert horiz height width a.
(Measure meas, C vert, C horiz) =>
Full meas vert horiz height width a -> General height width a
Matrix.fromFull Full meas horiz vert sizeB sizeA a
c) Square sizeB a
d


square :: (Shape.C sh, Class.Floating a) => Square sh a -> Square sh a
square :: Square sh a -> Square sh a
square Square sh a
a = Square sh a -> Square sh a -> Square sh a
forall sh a.
(C sh, Floating a) =>
Square sh a -> Square sh a -> Square sh a
multiplyCommutativeUnchecked Square sh a
a Square sh a
a

power ::
   (Shape.C sh, Class.Floating a) =>
   Integer -> Square sh a -> Square sh a
power :: Integer -> Square sh a -> Square sh a
power Integer
n Square sh a
a =
   (Square sh a -> Square sh a -> Square sh a)
-> Square sh a -> Square sh a -> Integer -> Square sh a
forall a. (a -> a -> a) -> a -> a -> Integer -> a
powerAssociative Square sh a -> Square sh a -> Square sh a
forall sh a.
(C sh, Floating a) =>
Square sh a -> Square sh a -> Square sh a
multiplyCommutativeUnchecked (Square sh a -> Square sh a
forall sh a. (C sh, Floating a) => Square sh a -> Square sh a
identityFrom Square sh a
a) Square sh a
a Integer
n

congruence ::
   (Shape.C height, Eq height, Shape.C width, Class.Floating a) =>
   Square height a -> General height width a -> Square width a
congruence :: Square height a -> General height width a -> Square width a
congruence Square height a
b General height width a
a0 =
   let a :: Full Size Big Big height (Unchecked width) a
a = (width -> Unchecked width)
-> General height width a
-> Full Size Big Big height (Unchecked width) a
forall vert horiz widthA widthB height a.
(C vert, C horiz) =>
(widthA -> widthB)
-> Full Size vert horiz height widthA a
-> Full Size vert horiz height widthB a
Basic.mapWidth width -> Unchecked width
forall sh. sh -> Unchecked sh
Unchecked General height width a
a0
   in (Unchecked width -> width)
-> Square (Unchecked width) a -> Square width a
forall sh0 sh1 a. (sh0 -> sh1) -> Square sh0 a -> Square sh1 a
mapSize (\(Unchecked width
sh) -> width
sh) (Square (Unchecked width) a -> Square width a)
-> Square (Unchecked width) a -> Square width a
forall a b. (a -> b) -> a -> b
$ Full Size Big Big (Unchecked width) (Unchecked width) a
-> Square (Unchecked width) a
forall meas vert horiz sh a.
(Measure meas, C vert, C horiz, Eq sh) =>
Full meas vert horiz sh sh a -> Square sh a
fromFull (Full Size Big Big (Unchecked width) (Unchecked width) a
 -> Square (Unchecked width) a)
-> Full Size Big Big (Unchecked width) (Unchecked width) a
-> Square (Unchecked width) a
forall a b. (a -> b) -> a -> b
$
      Full Size Big Big (Unchecked width) height a
-> Full Size Big Big height (Unchecked width) a
-> Full Size Big Big (Unchecked width) (Unchecked width) a
forall meas vert horiz height fuse width a.
(Measure meas, C vert, C horiz, C height, C fuse, Eq fuse, C width,
 Floating a) =>
Full meas vert horiz height fuse a
-> Full meas vert horiz fuse width a
-> Full meas vert horiz height width a
Basic.multiply (Full Size Big Big height (Unchecked width) a
-> Full Size Big Big (Unchecked width) height a
forall meas vert horiz height width a.
(Measure meas, C vert, C horiz, C height, C width, Floating a) =>
Full meas vert horiz height width a
-> Full meas horiz vert width height a
Basic.adjoint Full Size Big Big height (Unchecked width) a
a) (Full Size Big Big height (Unchecked width) a
 -> Full Size Big Big (Unchecked width) (Unchecked width) a)
-> Full Size Big Big height (Unchecked width) a
-> Full Size Big Big (Unchecked width) (Unchecked width) a
forall a b. (a -> b) -> a -> b
$ Full Size Big Big height height a
-> Full Size Big Big height (Unchecked width) a
-> Full Size Big Big height (Unchecked width) a
forall meas vert horiz height fuse width a.
(Measure meas, C vert, C horiz, C height, C fuse, Eq fuse, C width,
 Floating a) =>
Full meas vert horiz height fuse a
-> Full meas vert horiz fuse width a
-> Full meas vert horiz height width a
Basic.multiply (Square height a -> Full Size Big Big height height a
forall meas vert horiz sh a.
(Measure meas, C vert, C horiz) =>
Square sh a -> Full meas vert horiz sh sh a
toFull Square height a
b) Full Size Big Big height (Unchecked width) a
a

congruenceAdjoint ::
   (Shape.C height, Shape.C width, Eq width, Class.Floating a) =>
   General height width a -> Square width a -> Square height a
congruenceAdjoint :: General height width a -> Square width a -> Square height a
congruenceAdjoint General height width a
a0 Square width a
b =
   let a :: Full Size Big Big (Unchecked height) width a
a = (height -> Unchecked height)
-> General height width a
-> Full Size Big Big (Unchecked height) width a
forall vert horiz heightA heightB width a.
(C vert, C horiz) =>
(heightA -> heightB)
-> Full Size vert horiz heightA width a
-> Full Size vert horiz heightB width a
Basic.mapHeight height -> Unchecked height
forall sh. sh -> Unchecked sh
Unchecked General height width a
a0
   in (Unchecked height -> height)
-> Square (Unchecked height) a -> Square height a
forall sh0 sh1 a. (sh0 -> sh1) -> Square sh0 a -> Square sh1 a
mapSize (\(Unchecked height
sh) -> height
sh) (Square (Unchecked height) a -> Square height a)
-> Square (Unchecked height) a -> Square height a
forall a b. (a -> b) -> a -> b
$ Full Size Big Big (Unchecked height) (Unchecked height) a
-> Square (Unchecked height) a
forall meas vert horiz sh a.
(Measure meas, C vert, C horiz, Eq sh) =>
Full meas vert horiz sh sh a -> Square sh a
fromFull (Full Size Big Big (Unchecked height) (Unchecked height) a
 -> Square (Unchecked height) a)
-> Full Size Big Big (Unchecked height) (Unchecked height) a
-> Square (Unchecked height) a
forall a b. (a -> b) -> a -> b
$
      Full Size Big Big (Unchecked height) width a
-> Full Size Big Big width (Unchecked height) a
-> Full Size Big Big (Unchecked height) (Unchecked height) a
forall meas vert horiz height fuse width a.
(Measure meas, C vert, C horiz, C height, C fuse, Eq fuse, C width,
 Floating a) =>
Full meas vert horiz height fuse a
-> Full meas vert horiz fuse width a
-> Full meas vert horiz height width a
Basic.multiply Full Size Big Big (Unchecked height) width a
a (Full Size Big Big width (Unchecked height) a
 -> Full Size Big Big (Unchecked height) (Unchecked height) a)
-> Full Size Big Big width (Unchecked height) a
-> Full Size Big Big (Unchecked height) (Unchecked height) a
forall a b. (a -> b) -> a -> b
$ Full Size Big Big width width a
-> Full Size Big Big width (Unchecked height) a
-> Full Size Big Big width (Unchecked height) a
forall meas vert horiz height fuse width a.
(Measure meas, C vert, C horiz, C height, C fuse, Eq fuse, C width,
 Floating a) =>
Full meas vert horiz height fuse a
-> Full meas vert horiz fuse width a
-> Full meas vert horiz height width a
Basic.multiply (Square width a -> Full Size Big Big width width a
forall meas vert horiz sh a.
(Measure meas, C vert, C horiz) =>
Square sh a -> Full meas vert horiz sh sh a
toFull Square width a
b) (Full Size Big Big width (Unchecked height) a
 -> Full Size Big Big width (Unchecked height) a)
-> Full Size Big Big width (Unchecked height) a
-> Full Size Big Big width (Unchecked height) a
forall a b. (a -> b) -> a -> b
$ Full Size Big Big (Unchecked height) width a
-> Full Size Big Big width (Unchecked height) a
forall meas vert horiz height width a.
(Measure meas, C vert, C horiz, C height, C width, Floating a) =>
Full meas vert horiz height width a
-> Full meas horiz vert width height a
Basic.adjoint Full Size Big Big (Unchecked height) width a
a

{-
orderA and orderB must be equal but this is not checked.
-}
multiplyCommutativeUnchecked ::
   (Shape.C sh, Class.Floating a) =>
   Square sh a -> Square sh a -> Square sh a
multiplyCommutativeUnchecked :: Square sh a -> Square sh a -> Square sh a
multiplyCommutativeUnchecked
   (Array shape :: Square sh
shape@(Layout.Full Order
order Extent Shape Small Small sh sh
extent) ForeignPtr a
a)
   (Array Square sh
_ ForeignPtr a
b) =
      Square sh -> (Ptr a -> IO ()) -> Square sh a
forall sh a.
(C sh, Storable a) =>
sh -> (Ptr a -> IO ()) -> Array sh a
Array.unsafeCreate Square sh
shape ((Ptr a -> IO ()) -> Square sh a)
-> (Ptr a -> IO ()) -> Square sh a
forall a b. (a -> b) -> a -> b
$ \Ptr a
cPtr ->
   let n :: Int
n = sh -> Int
forall sh. C sh => sh -> Int
Shape.size (sh -> Int) -> sh -> Int
forall a b. (a -> b) -> a -> b
$ Extent Shape Small Small sh sh -> sh
forall meas vert horiz height width.
(Measure meas, C vert, C horiz) =>
Extent meas vert horiz height width -> height
Extent.height Extent Shape Small Small sh sh
extent
       (ForeignPtr a
at,ForeignPtr a
bt) = Order
-> (ForeignPtr a, ForeignPtr a) -> (ForeignPtr a, ForeignPtr a)
forall a. Order -> (a, a) -> (a, a)
swapOnRowMajor Order
order (ForeignPtr a
a,ForeignPtr a
b)
   in  Order
-> Order
-> Int
-> Int
-> Int
-> ForeignPtr a
-> ForeignPtr a
-> Ptr a
-> IO ()
forall a.
Floating a =>
Order
-> Order
-> Int
-> Int
-> Int
-> ForeignPtr a
-> ForeignPtr a
-> Ptr a
-> IO ()
Private.multiplyMatrix Order
ColumnMajor Order
ColumnMajor Int
n Int
n Int
n ForeignPtr a
at ForeignPtr a
bt Ptr a
cPtr