{-# LANGUAGE TypeFamilies #-}
module Numeric.LAPACK.Orthogonal.Basic (
   leastSquares,
   minimumNorm,
   leastSquaresMinimumNormRCond,
   pseudoInverseRCond,

   leastSquaresConstraint,
   gaussMarkovLinearModel,

   determinantAbsolute,
   complement,
   extractComplement,
   ) where

import qualified Numeric.LAPACK.Orthogonal.Plain as HH

import qualified Numeric.LAPACK.Matrix.Square.Basic as Square
import qualified Numeric.LAPACK.Matrix.Layout.Private as Layout
import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
import qualified Numeric.LAPACK.Matrix.Basic as Basic
import qualified Numeric.LAPACK.Vector as Vector
import Numeric.LAPACK.Matrix.Layout.Private (Order(RowMajor,ColumnMajor))
import Numeric.LAPACK.Matrix.Private
         (Full, General, Tall, Wide, ShapeInt, shapeInt)
import Numeric.LAPACK.Vector (Vector)
import Numeric.LAPACK.Scalar (RealOf, zero, absolute)
import Numeric.LAPACK.Private
         (lacgv, peekCInt,
          copySubMatrix, copyToTemp,
          copyToColumnMajorTemp, copyToSubColumnMajor,
          withAutoWorkspaceInfo, rankMsg, errorCodeMsg, createHigherArray)

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

import qualified Data.Array.Comfort.Storable.Unchecked.Monadic as ArrayIO
import qualified Data.Array.Comfort.Storable.Unchecked as Array
import qualified Data.Array.Comfort.Shape as Shape
import Data.Array.Comfort.Storable.Unchecked (Array(Array))

import System.IO.Unsafe (unsafePerformIO)

import Foreign.Marshal.Array (pokeArray)
import Foreign.C.Types (CInt, CChar)
import Foreign.ForeignPtr (ForeignPtr, withForeignPtr)
import Foreign.Ptr (Ptr)

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

import Data.Complex (Complex)
import Data.Tuple.HT (mapSnd)


leastSquares ::
   (Extent.Measure meas, Extent.C vert, Extent.C horiz,
    Shape.C height, Eq height, Shape.C width, Shape.C nrhs, Class.Floating a) =>
   Full meas horiz Extent.Small height width a ->
   Full meas vert horiz height nrhs a ->
   Full meas vert horiz width nrhs a
leastSquares :: forall meas vert horiz height width nrhs a.
(Measure meas, C vert, C horiz, C height, Eq height, C width,
 C nrhs, Floating a) =>
Full meas horiz Small height width a
-> Full meas vert horiz height nrhs a
-> Full meas vert horiz width nrhs a
leastSquares
   (Array shapeA :: Full meas horiz Small height width
shapeA@(Layout.Full Order
orderA Extent meas horiz Small height width
extentA) ForeignPtr a
a)
   (Array        (Layout.Full Order
orderB Extent meas vert horiz height nrhs
extentB) ForeignPtr a
b) =

 case Extent meas vert horiz width height
-> Extent meas vert horiz height nrhs
-> Maybe (Extent meas vert horiz width nrhs)
forall meas vert horiz fuse height width.
(Measure meas, C vert, C horiz, Eq fuse) =>
Extent meas vert horiz height fuse
-> Extent meas vert horiz fuse width
-> Maybe (Extent meas vert horiz height width)
Extent.fuse (Extent meas Small horiz width height
-> Extent meas vert horiz width height
forall meas vert horiz height width.
(Measure meas, C vert, C horiz) =>
Extent meas Small horiz height width
-> Extent meas vert horiz height width
Extent.weakenWide (Extent meas Small horiz width height
 -> Extent meas vert horiz width height)
-> Extent meas Small horiz width height
-> Extent meas vert horiz width height
forall a b. (a -> b) -> a -> b
$ Extent meas horiz Small height width
-> Extent meas Small horiz width height
forall meas vert horiz height width.
(Measure meas, C vert, C horiz) =>
Extent meas vert horiz height width
-> Extent meas horiz vert width height
Extent.transpose Extent meas horiz Small height width
extentA) Extent meas vert horiz height nrhs
extentB of
  Maybe (Extent meas vert horiz width nrhs)
Nothing -> [Char] -> Full meas vert horiz width nrhs a
forall a. HasCallStack => [Char] -> a
error [Char]
"leastSquares: height shapes mismatch"
  Just Extent meas vert horiz width nrhs
extent ->
      Full meas vert horiz width nrhs
-> (Ptr a -> IO ()) -> Full meas vert horiz width nrhs a
forall sh a.
(C sh, Storable a) =>
sh -> (Ptr a -> IO ()) -> Array sh a
Array.unsafeCreate (Order
-> Extent meas vert horiz width nrhs
-> Full meas vert horiz width nrhs
forall meas vert horiz height width.
Order
-> Extent meas vert horiz height width
-> Full meas vert horiz height width
Layout.Full Order
ColumnMajor Extent meas vert horiz width nrhs
extent) ((Ptr a -> IO ()) -> Full meas vert horiz width nrhs a)
-> (Ptr a -> IO ()) -> Full meas vert horiz width nrhs a
forall a b. (a -> b) -> a -> b
$ \Ptr a
xPtr -> do

   let widthA :: width
widthA = Extent meas horiz Small height width -> width
forall meas vert horiz height width.
(Measure meas, C vert, C horiz) =>
Extent meas vert horiz height width -> width
Extent.width Extent meas horiz Small height width
extentA
   let (height
height,nrhs
widthB) = Extent meas vert horiz height nrhs -> (height, nrhs)
forall meas vert horiz height width.
(Measure meas, C vert, C horiz) =>
Extent meas vert horiz height width -> (height, width)
Extent.dimensions Extent meas vert horiz height nrhs
extentB
   let (Int
m,Int
n) = Full meas horiz Small height width -> (Int, Int)
forall meas vert horiz height width.
(Measure meas, C vert, C horiz, C height, C width) =>
Full meas vert horiz height width -> (Int, Int)
Layout.dimensions Full meas horiz Small height width
shapeA
   let lda :: Int
lda = Int
m
   let nrhs :: Int
nrhs = nrhs -> Int
forall sh. C sh => sh -> Int
Shape.size nrhs
widthB
   let ldb :: Int
ldb = height -> Int
forall sh. C sh => sh -> Int
Shape.size height
height
   let ldx :: Int
ldx = width -> Int
forall sh. C sh => sh -> Int
Shape.size width
widthA
   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
mPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
m
      Ptr CInt
nPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
      Ptr CInt
nrhsPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
nrhs
      (Ptr CChar
transPtr,Ptr a
aPtr) <- Order -> Int -> ForeignPtr a -> ContT () IO (Ptr CChar, Ptr a)
forall a r.
Floating a =>
Order -> Int -> ForeignPtr a -> ContT r IO (Ptr CChar, Ptr a)
adjointA Order
orderA (Int
mInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
n) ForeignPtr a
a
      Ptr CInt
ldaPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
lda
      Ptr CInt
ldbPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
ldb
      Ptr a
bPtr <- Order -> Int -> Int -> ForeignPtr a -> ContT () IO (Ptr a)
forall a r.
Floating a =>
Order -> Int -> Int -> ForeignPtr a -> ContT r IO (Ptr a)
copyToColumnMajorTemp Order
orderB Int
ldb Int
nrhs ForeignPtr a
b
      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
$ [Char]
-> [Char] -> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
forall a.
Floating a =>
[Char]
-> [Char] -> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
withAutoWorkspaceInfo [Char]
rankMsg [Char]
"gels" ((Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ())
-> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$
         Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> IO ()
forall a.
Floating a =>
Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> IO ()
LapackGen.gels Ptr CChar
transPtr
            Ptr CInt
mPtr Ptr CInt
nPtr Ptr CInt
nrhsPtr Ptr a
aPtr Ptr CInt
ldaPtr Ptr a
bPtr Ptr CInt
ldbPtr
      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
$ Int -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
forall a.
Floating a =>
Int -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
copySubMatrix Int
ldx Int
nrhs Int
ldb Ptr a
bPtr Int
ldx Ptr a
xPtr

minimumNorm ::
   (Extent.Measure meas, Extent.C vert, Extent.C horiz,
    Shape.C height, Eq height, Shape.C width, Shape.C nrhs, Class.Floating a) =>
   Full meas Extent.Small vert height width a ->
   Full meas vert horiz height nrhs a ->
   Full meas vert horiz width nrhs a
minimumNorm :: forall meas vert horiz height width nrhs a.
(Measure meas, C vert, C horiz, C height, Eq height, C width,
 C nrhs, Floating a) =>
Full meas Small vert height width a
-> Full meas vert horiz height nrhs a
-> Full meas vert horiz width nrhs a
minimumNorm
   (Array shapeA :: Full meas Small vert height width
shapeA@(Layout.Full Order
orderA Extent meas Small vert height width
extentA) ForeignPtr a
a)
   (Array        (Layout.Full Order
orderB Extent meas vert horiz height nrhs
extentB) ForeignPtr a
b) =

 case Extent meas vert horiz width height
-> Extent meas vert horiz height nrhs
-> Maybe (Extent meas vert horiz width nrhs)
forall meas vert horiz fuse height width.
(Measure meas, C vert, C horiz, Eq fuse) =>
Extent meas vert horiz height fuse
-> Extent meas vert horiz fuse width
-> Maybe (Extent meas vert horiz height width)
Extent.fuse (Extent meas vert Small width height
-> Extent meas vert horiz width height
forall meas vert horiz height width.
(Measure meas, C vert, C horiz) =>
Extent meas vert Small height width
-> Extent meas vert horiz height width
Extent.weakenTall (Extent meas vert Small width height
 -> Extent meas vert horiz width height)
-> Extent meas vert Small width height
-> Extent meas vert horiz width height
forall a b. (a -> b) -> a -> b
$ Extent meas Small vert height width
-> Extent meas vert Small width height
forall meas vert horiz height width.
(Measure meas, C vert, C horiz) =>
Extent meas vert horiz height width
-> Extent meas horiz vert width height
Extent.transpose Extent meas Small vert height width
extentA) Extent meas vert horiz height nrhs
extentB of
  Maybe (Extent meas vert horiz width nrhs)
Nothing -> [Char] -> Full meas vert horiz width nrhs a
forall a. HasCallStack => [Char] -> a
error [Char]
"minimumNorm: height shapes mismatch"
  Just Extent meas vert horiz width nrhs
extent ->
      Full meas vert horiz width nrhs
-> (Ptr a -> IO ()) -> Full meas vert horiz width nrhs a
forall sh a.
(C sh, Storable a) =>
sh -> (Ptr a -> IO ()) -> Array sh a
Array.unsafeCreate (Order
-> Extent meas vert horiz width nrhs
-> Full meas vert horiz width nrhs
forall meas vert horiz height width.
Order
-> Extent meas vert horiz height width
-> Full meas vert horiz height width
Layout.Full Order
ColumnMajor Extent meas vert horiz width nrhs
extent) ((Ptr a -> IO ()) -> Full meas vert horiz width nrhs a)
-> (Ptr a -> IO ()) -> Full meas vert horiz width nrhs a
forall a b. (a -> b) -> a -> b
$ \Ptr a
xPtr -> do

   let widthA :: width
widthA = Extent meas Small vert height width -> width
forall meas vert horiz height width.
(Measure meas, C vert, C horiz) =>
Extent meas vert horiz height width -> width
Extent.width Extent meas Small vert height width
extentA
   let (height
height,nrhs
widthB) = Extent meas vert horiz height nrhs -> (height, nrhs)
forall meas vert horiz height width.
(Measure meas, C vert, C horiz) =>
Extent meas vert horiz height width -> (height, width)
Extent.dimensions Extent meas vert horiz height nrhs
extentB
   let (Int
m,Int
n) = Full meas Small vert height width -> (Int, Int)
forall meas vert horiz height width.
(Measure meas, C vert, C horiz, C height, C width) =>
Full meas vert horiz height width -> (Int, Int)
Layout.dimensions Full meas Small vert height width
shapeA
   let lda :: Int
lda = Int
m
   let nrhs :: Int
nrhs = nrhs -> Int
forall sh. C sh => sh -> Int
Shape.size nrhs
widthB
   let ldb :: Int
ldb = height -> Int
forall sh. C sh => sh -> Int
Shape.size height
height
   let ldx :: Int
ldx = width -> Int
forall sh. C sh => sh -> Int
Shape.size width
widthA
   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
mPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
m
      Ptr CInt
nPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
      Ptr CInt
nrhsPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
nrhs
      (Ptr CChar
transPtr,Ptr a
aPtr) <- Order -> Int -> ForeignPtr a -> ContT () IO (Ptr CChar, Ptr a)
forall a r.
Floating a =>
Order -> Int -> ForeignPtr a -> ContT r IO (Ptr CChar, Ptr a)
adjointA Order
orderA (Int
mInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
n) ForeignPtr a
a
      Ptr CInt
ldaPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
lda
      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
ldx
      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
$ Order -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
forall a.
Floating a =>
Order -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
copyToSubColumnMajor Order
orderB Int
ldb Int
nrhs Ptr a
bPtr Int
ldx Ptr a
xPtr
      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
$ [Char]
-> [Char] -> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
forall a.
Floating a =>
[Char]
-> [Char] -> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
withAutoWorkspaceInfo [Char]
rankMsg [Char]
"gels" ((Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ())
-> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$
         Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> IO ()
forall a.
Floating a =>
Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> IO ()
LapackGen.gels Ptr CChar
transPtr
            Ptr CInt
mPtr Ptr CInt
nPtr Ptr CInt
nrhsPtr Ptr a
aPtr Ptr CInt
ldaPtr Ptr a
xPtr Ptr CInt
ldxPtr


adjointA ::
   Class.Floating a =>
   Order -> Int -> ForeignPtr a -> ContT r IO (Ptr CChar, Ptr a)
adjointA :: forall a r.
Floating a =>
Order -> Int -> ForeignPtr a -> ContT r IO (Ptr CChar, Ptr a)
adjointA Order
order Int
size ForeignPtr a
a = do
   Ptr a
aPtr <- Int -> ForeignPtr a -> ContT r IO (Ptr a)
forall a r. Storable a => Int -> ForeignPtr a -> ContT r IO (Ptr a)
copyToTemp Int
size ForeignPtr a
a
   Char
trans <-
      case Order
order of
         Order
RowMajor -> do
            Ptr CInt
sizePtr <- Int -> FortranIO r (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
size
            Ptr CInt
incPtr <- Int -> FortranIO r (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1
            IO () -> ContT r IO ()
forall a. IO a -> ContT r IO a
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> ContT r IO ()) -> IO () -> ContT r IO ()
forall a b. (a -> b) -> a -> b
$ Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
forall a. Floating a => Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
lacgv Ptr CInt
sizePtr Ptr a
aPtr Ptr CInt
incPtr
            Char -> ContT r IO Char
forall a. a -> ContT r IO a
forall (m :: * -> *) a. Monad m => a -> m a
return (Char -> ContT r IO Char) -> Char -> ContT r IO Char
forall a b. (a -> b) -> a -> b
$ ForeignPtr a -> Char
forall a (f :: * -> *). Floating a => f a -> Char
HH.invChar ForeignPtr a
a
         Order
ColumnMajor -> Char -> ContT r IO Char
forall a. a -> ContT r IO a
forall (m :: * -> *) a. Monad m => a -> m a
return Char
'N'
   Ptr CChar
transPtr <- Char -> FortranIO r (Ptr CChar)
forall r. Char -> FortranIO r (Ptr CChar)
Call.char Char
trans
   (Ptr CChar, Ptr a) -> ContT r IO (Ptr CChar, Ptr a)
forall a. a -> ContT r IO a
forall (m :: * -> *) a. Monad m => a -> m a
return (Ptr CChar
transPtr, Ptr a
aPtr)


leastSquaresMinimumNormRCond ::
   (Extent.Measure meas, Extent.C vert, Extent.C horiz,
    Shape.C height, Eq height, Shape.C width, Shape.C nrhs, Class.Floating a) =>
   RealOf a ->
   Full meas horiz vert height width a ->
   Full meas vert horiz height nrhs a ->
   (Int, Full meas vert horiz width nrhs a)
leastSquaresMinimumNormRCond :: forall meas vert horiz height width nrhs a.
(Measure meas, C vert, C horiz, C height, Eq height, C width,
 C nrhs, Floating a) =>
RealOf a
-> Full meas horiz vert height width a
-> Full meas vert horiz height nrhs a
-> (Int, Full meas vert horiz width nrhs a)
leastSquaresMinimumNormRCond RealOf a
rcond
      (Array (Layout.Full Order
orderA Extent meas horiz vert height width
extentA) ForeignPtr a
a)
      (Array (Layout.Full Order
orderB Extent meas vert horiz height nrhs
extentB) ForeignPtr a
b) =
   case Extent meas vert horiz width height
-> Extent meas vert horiz height nrhs
-> Maybe (Extent meas vert horiz width nrhs)
forall meas vert horiz fuse height width.
(Measure meas, C vert, C horiz, Eq fuse) =>
Extent meas vert horiz height fuse
-> Extent meas vert horiz fuse width
-> Maybe (Extent meas vert horiz height width)
Extent.fuse (Extent meas horiz vert height width
-> Extent meas vert horiz width height
forall meas vert horiz height width.
(Measure meas, C vert, C horiz) =>
Extent meas vert horiz height width
-> Extent meas horiz vert width height
Extent.transpose Extent meas horiz vert height width
extentA) Extent meas vert horiz height nrhs
extentB of
      Maybe (Extent meas vert horiz width nrhs)
Nothing -> [Char] -> (Int, Full meas vert horiz width nrhs a)
forall a. HasCallStack => [Char] -> a
error [Char]
"leastSquaresMinimumNormRCond: height shapes mismatch"
      Just Extent meas vert horiz width nrhs
extent ->
         let widthA :: width
widthA = Extent meas horiz vert height width -> width
forall meas vert horiz height width.
(Measure meas, C vert, C horiz) =>
Extent meas vert horiz height width -> width
Extent.width Extent meas horiz vert height width
extentA
             (height
height,nrhs
widthB) = Extent meas vert horiz height nrhs -> (height, nrhs)
forall meas vert horiz height width.
(Measure meas, C vert, C horiz) =>
Extent meas vert horiz height width -> (height, width)
Extent.dimensions Extent meas vert horiz height nrhs
extentB
             shapeX :: Full meas vert horiz width nrhs
shapeX = Order
-> Extent meas vert horiz width nrhs
-> Full meas vert horiz width nrhs
forall meas vert horiz height width.
Order
-> Extent meas vert horiz height width
-> Full meas vert horiz height width
Layout.Full Order
ColumnMajor Extent meas vert horiz width nrhs
extent
             m :: Int
m = height -> Int
forall sh. C sh => sh -> Int
Shape.size height
height
             n :: Int
n = width -> Int
forall sh. C sh => sh -> Int
Shape.size width
widthA
             nrhs :: Int
nrhs = nrhs -> Int
forall sh. C sh => sh -> Int
Shape.size nrhs
widthB
         in  if Int
m Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0
                then (Int
0, Full meas vert horiz width nrhs
-> Full meas vert horiz width nrhs a
forall sh a. (C sh, Floating a) => sh -> Vector sh a
Vector.zero Full meas vert horiz width nrhs
shapeX)
                else
                  if Int
nrhs Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0
                     then
                        ((Int, Array (General width ()) a) -> Int
forall a b. (a, b) -> a
fst ((Int, Array (General width ()) a) -> Int)
-> (Int, Array (General width ()) a) -> Int
forall a b. (a -> b) -> a -> b
$ IO (Int, Array (General width ()) a)
-> (Int, Array (General width ()) a)
forall a. IO a -> a
unsafePerformIO (IO (Int, Array (General width ()) a)
 -> (Int, Array (General width ()) a))
-> IO (Int, Array (General width ()) a)
-> (Int, Array (General width ()) a)
forall a b. (a -> b) -> a -> b
$
                         case height -> Vector height a
forall sh a. (C sh, Floating a) => sh -> Vector sh a
Vector.zero height
height of
                           Array height
_ ForeignPtr a
b1 ->
                              RealOf a
-> General width ()
-> Order
-> ForeignPtr a
-> Order
-> ForeignPtr a
-> Int
-> Int
-> Int
-> IO (Int, Array (General width ()) a)
forall sh a.
(C sh, Floating a) =>
RealOf a
-> sh
-> Order
-> ForeignPtr a
-> Order
-> ForeignPtr a
-> Int
-> Int
-> Int
-> IO (Int, Array sh a)
leastSquaresMinimumNormIO RealOf a
rcond
                                 (Order -> width -> () -> General width ()
forall height width.
Order -> height -> width -> General height width
Layout.general Order
ColumnMajor width
widthA ())
                                 Order
orderA ForeignPtr a
a Order
orderB ForeignPtr a
b1 Int
m Int
n Int
1,
                         Full meas vert horiz width nrhs
-> Full meas vert horiz width nrhs a
forall sh a. (C sh, Floating a) => sh -> Vector sh a
Vector.zero Full meas vert horiz width nrhs
shapeX)
                     else
                        IO (Int, Full meas vert horiz width nrhs a)
-> (Int, Full meas vert horiz width nrhs a)
forall a. IO a -> a
unsafePerformIO (IO (Int, Full meas vert horiz width nrhs a)
 -> (Int, Full meas vert horiz width nrhs a))
-> IO (Int, Full meas vert horiz width nrhs a)
-> (Int, Full meas vert horiz width nrhs a)
forall a b. (a -> b) -> a -> b
$
                        RealOf a
-> Full meas vert horiz width nrhs
-> Order
-> ForeignPtr a
-> Order
-> ForeignPtr a
-> Int
-> Int
-> Int
-> IO (Int, Full meas vert horiz width nrhs a)
forall sh a.
(C sh, Floating a) =>
RealOf a
-> sh
-> Order
-> ForeignPtr a
-> Order
-> ForeignPtr a
-> Int
-> Int
-> Int
-> IO (Int, Array sh a)
leastSquaresMinimumNormIO RealOf a
rcond Full meas vert horiz width nrhs
shapeX
                           Order
orderA ForeignPtr a
a Order
orderB ForeignPtr a
b Int
m Int
n Int
nrhs

leastSquaresMinimumNormIO ::
   (Shape.C sh, Class.Floating a) =>
   RealOf a -> sh ->
   Order -> ForeignPtr a ->
   Order -> ForeignPtr a ->
   Int -> Int -> Int -> IO (Int, Array sh a)
leastSquaresMinimumNormIO :: forall sh a.
(C sh, Floating a) =>
RealOf a
-> sh
-> Order
-> ForeignPtr a
-> Order
-> ForeignPtr a
-> Int
-> Int
-> Int
-> IO (Int, Array sh a)
leastSquaresMinimumNormIO RealOf a
rcond sh
shapeX Order
orderA ForeignPtr a
a Order
orderB ForeignPtr a
b Int
m Int
n Int
nrhs =
   sh
-> Int
-> Int
-> Int
-> ((Ptr a, Int) -> IO Int)
-> IO (Int, Array sh a)
forall sh a rank.
(C sh, Floating a) =>
sh
-> Int
-> Int
-> Int
-> ((Ptr a, Int) -> IO rank)
-> IO (rank, Array sh a)
createHigherArray sh
shapeX Int
m Int
n Int
nrhs (((Ptr a, Int) -> IO Int) -> IO (Int, Array sh a))
-> ((Ptr a, Int) -> IO Int) -> IO (Int, Array sh a)
forall a b. (a -> b) -> a -> b
$ \(Ptr a
tmpPtr,Int
ldtmp) -> do

   let lda :: Int
lda = Int
m
   ContT Int IO Int -> IO Int
forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT (ContT Int IO Int -> IO Int) -> ContT Int IO Int -> IO Int
forall a b. (a -> b) -> a -> b
$ do
      Ptr a
aPtr <- Order -> Int -> Int -> ForeignPtr a -> ContT Int IO (Ptr a)
forall a r.
Floating a =>
Order -> Int -> Int -> ForeignPtr a -> ContT r IO (Ptr a)
copyToColumnMajorTemp Order
orderA Int
m Int
n ForeignPtr a
a
      Ptr CInt
ldaPtr <- Int -> FortranIO Int (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
lda
      Ptr CInt
ldtmpPtr <- Int -> FortranIO Int (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
ldtmp
      Ptr a
bPtr <- ((Ptr a -> IO Int) -> IO Int) -> ContT Int IO (Ptr a)
forall {k} (r :: k) (m :: k -> *) a.
((a -> m r) -> m r) -> ContT r m a
ContT (((Ptr a -> IO Int) -> IO Int) -> ContT Int IO (Ptr a))
-> ((Ptr a -> IO Int) -> IO Int) -> ContT Int IO (Ptr a)
forall a b. (a -> b) -> a -> b
$ ForeignPtr a -> (Ptr a -> IO Int) -> IO Int
forall a b. ForeignPtr a -> (Ptr a -> IO b) -> IO b
withForeignPtr ForeignPtr a
b
      IO () -> ContT Int IO ()
forall a. IO a -> ContT Int IO a
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> ContT Int IO ()) -> IO () -> ContT Int IO ()
forall a b. (a -> b) -> a -> b
$ Order -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
forall a.
Floating a =>
Order -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
copyToSubColumnMajor Order
orderB Int
m Int
nrhs Ptr a
bPtr Int
ldtmp Ptr a
tmpPtr
      Ptr CInt
jpvtPtr <- Int -> FortranIO Int (Ptr CInt)
forall a r. Storable a => Int -> FortranIO r (Ptr a)
Call.allocaArray Int
n
      IO () -> ContT Int IO ()
forall a. IO a -> ContT Int IO a
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> ContT Int IO ()) -> IO () -> ContT Int IO ()
forall a b. (a -> b) -> a -> b
$ Ptr CInt -> [CInt] -> IO ()
forall a. Storable a => Ptr a -> [a] -> IO ()
pokeArray Ptr CInt
jpvtPtr (Int -> CInt -> [CInt]
forall a. Int -> a -> [a]
replicate Int
n CInt
0)
      Ptr CInt
rankPtr <- FortranIO Int (Ptr CInt)
forall a r. Storable a => FortranIO r (Ptr a)
Call.alloca
      GELSY_ Int (RealOf a) a
forall a r. Floating a => GELSY_ r (RealOf a) a
gelsy Int
m Int
n Int
nrhs Ptr a
aPtr Ptr CInt
ldaPtr Ptr a
tmpPtr Ptr CInt
ldtmpPtr Ptr CInt
jpvtPtr RealOf a
rcond Ptr CInt
rankPtr
      IO Int -> ContT Int IO Int
forall a. IO a -> ContT Int IO a
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO Int -> ContT Int IO Int) -> IO Int -> ContT Int IO Int
forall a b. (a -> b) -> a -> b
$ Ptr CInt -> IO Int
peekCInt Ptr CInt
rankPtr


type GELSY_ r ar a =
   Int -> Int -> Int -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt ->
   Ptr CInt -> ar -> Ptr CInt -> ContT r IO ()

newtype GELSY r a = GELSY {forall r a. GELSY r a -> GELSY_ r (RealOf a) a
getGELSY :: GELSY_ r (RealOf a) a}

gelsy :: (Class.Floating a) => GELSY_ r (RealOf a) a
gelsy :: forall a r. Floating a => GELSY_ r (RealOf a) a
gelsy =
   GELSY r a
-> Int
-> Int
-> Int
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> RealOf a
-> Ptr CInt
-> ContT r IO ()
forall r a. GELSY r a -> GELSY_ r (RealOf a) a
getGELSY (GELSY r a
 -> Int
 -> Int
 -> Int
 -> Ptr a
 -> Ptr CInt
 -> Ptr a
 -> Ptr CInt
 -> Ptr CInt
 -> RealOf a
 -> Ptr CInt
 -> ContT r IO ())
-> GELSY r a
-> Int
-> Int
-> Int
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> RealOf a
-> Ptr CInt
-> ContT r IO ()
forall a b. (a -> b) -> a -> b
$
   GELSY r Float
-> GELSY r Double
-> GELSY r (Complex Float)
-> GELSY r (Complex Double)
-> GELSY r 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
      (GELSY_ r (RealOf Float) Float -> GELSY r Float
forall r a. GELSY_ r (RealOf a) a -> GELSY r a
GELSY GELSY_ r Float Float
GELSY_ r (RealOf Float) Float
forall a r. Real a => GELSY_ r a a
gelsyReal)
      (GELSY_ r (RealOf Double) Double -> GELSY r Double
forall r a. GELSY_ r (RealOf a) a -> GELSY r a
GELSY GELSY_ r Double Double
GELSY_ r (RealOf Double) Double
forall a r. Real a => GELSY_ r a a
gelsyReal)
      (GELSY_ r (RealOf (Complex Float)) (Complex Float)
-> GELSY r (Complex Float)
forall r a. GELSY_ r (RealOf a) a -> GELSY r a
GELSY GELSY_ r Float (Complex Float)
GELSY_ r (RealOf (Complex Float)) (Complex Float)
forall a r. Real a => GELSY_ r a (Complex a)
gelsyComplex)
      (GELSY_ r (RealOf (Complex Double)) (Complex Double)
-> GELSY r (Complex Double)
forall r a. GELSY_ r (RealOf a) a -> GELSY r a
GELSY GELSY_ r Double (Complex Double)
GELSY_ r (RealOf (Complex Double)) (Complex Double)
forall a r. Real a => GELSY_ r a (Complex a)
gelsyComplex)

gelsyReal :: (Class.Real a) => GELSY_ r a a
gelsyReal :: forall a r. Real a => GELSY_ r a a
gelsyReal Int
m Int
n Int
nrhs Ptr a
aPtr Ptr CInt
ldaPtr Ptr a
bPtr Ptr CInt
ldbPtr Ptr CInt
jpvtPtr a
rcond Ptr CInt
rankPtr = do
   Ptr CInt
mPtr <- Int -> FortranIO r (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
m
   Ptr CInt
nPtr <- Int -> FortranIO r (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
   Ptr CInt
nrhsPtr <- Int -> FortranIO r (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
nrhs
   Ptr a
rcondPtr <- a -> FortranIO r (Ptr a)
forall a r. Real a => a -> FortranIO r (Ptr a)
Call.real a
rcond
   IO () -> ContT r IO ()
forall a. IO a -> ContT r IO a
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> ContT r IO ()) -> IO () -> ContT r IO ()
forall a b. (a -> b) -> a -> b
$ [Char]
-> [Char] -> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
forall a.
Floating a =>
[Char]
-> [Char] -> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
withAutoWorkspaceInfo [Char]
errorCodeMsg [Char]
"gelsy" ((Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ())
-> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$
      Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> IO ()
forall a.
Real a =>
Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> IO ()
LapackReal.gelsy Ptr CInt
mPtr Ptr CInt
nPtr Ptr CInt
nrhsPtr
         Ptr a
aPtr Ptr CInt
ldaPtr Ptr a
bPtr Ptr CInt
ldbPtr Ptr CInt
jpvtPtr Ptr a
rcondPtr Ptr CInt
rankPtr

gelsyComplex :: (Class.Real a) => GELSY_ r a (Complex a)
gelsyComplex :: forall a r. Real a => GELSY_ r a (Complex a)
gelsyComplex Int
m Int
n Int
nrhs Ptr (Complex a)
aPtr Ptr CInt
ldaPtr Ptr (Complex a)
bPtr Ptr CInt
ldbPtr Ptr CInt
jpvtPtr a
rcond Ptr CInt
rankPtr = do
   Ptr CInt
mPtr <- Int -> FortranIO r (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
m
   Ptr CInt
nPtr <- Int -> FortranIO r (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
   Ptr CInt
nrhsPtr <- Int -> FortranIO r (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
nrhs
   Ptr a
rcondPtr <- a -> FortranIO r (Ptr a)
forall a r. Real a => a -> FortranIO r (Ptr a)
Call.real a
rcond
   Ptr a
rworkPtr <- Int -> FortranIO r (Ptr a)
forall a r. Storable a => Int -> FortranIO r (Ptr a)
Call.allocaArray (Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
n)
   IO () -> ContT r IO ()
forall a. IO a -> ContT r IO a
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> ContT r IO ()) -> IO () -> ContT r IO ()
forall a b. (a -> b) -> a -> b
$
      [Char]
-> [Char]
-> (Ptr (Complex a) -> Ptr CInt -> Ptr CInt -> IO ())
-> IO ()
forall a.
Floating a =>
[Char]
-> [Char] -> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
withAutoWorkspaceInfo [Char]
errorCodeMsg [Char]
"gelsy" ((Ptr (Complex a) -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ())
-> (Ptr (Complex a) -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \Ptr (Complex a)
workPtr Ptr CInt
lworkPtr Ptr CInt
infoPtr ->
      Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> IO ()
forall a.
Real a =>
Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> IO ()
LapackComplex.gelsy Ptr CInt
mPtr Ptr CInt
nPtr Ptr CInt
nrhsPtr
         Ptr (Complex a)
aPtr Ptr CInt
ldaPtr Ptr (Complex a)
bPtr Ptr CInt
ldbPtr Ptr CInt
jpvtPtr Ptr a
rcondPtr Ptr CInt
rankPtr
         Ptr (Complex a)
workPtr Ptr CInt
lworkPtr Ptr a
rworkPtr Ptr CInt
infoPtr


pseudoInverseRCond ::
   (Extent.Measure meas, Extent.C vert, Extent.C horiz,
    Shape.C height, Eq height, Shape.C width, Eq width, Class.Floating a) =>
   RealOf a ->
   Full meas vert horiz height width a ->
   (Int, Full meas horiz vert width height a)
pseudoInverseRCond :: forall meas vert horiz height width a.
(Measure meas, C vert, C horiz, C height, Eq height, C width,
 Eq width, Floating a) =>
RealOf a
-> Full meas vert horiz height width a
-> (Int, Full meas horiz vert width height a)
pseudoInverseRCond RealOf a
rcond Full meas vert horiz height width a
a =
   case Full meas vert horiz height width a
-> Either (Tall height width a) (Wide height width a)
forall meas vert horiz height width a.
(Measure meas, C vert, C horiz, C height, C width) =>
Full meas vert horiz height width a
-> Either (Tall height width a) (Wide height width a)
Basic.caseTallWide Full meas vert horiz height width a
a of
      Left Tall height width a
_ ->
         (Full meas vert horiz height width a
 -> Full meas horiz vert width height a)
-> (Int, Full meas vert horiz height width a)
-> (Int, Full meas horiz vert width height a)
forall b c a. (b -> c) -> (a, b) -> (a, c)
mapSnd Full meas vert horiz height width a
-> Full meas horiz vert width height a
forall meas vert horiz height width a.
(Measure meas, C vert, C horiz) =>
Full meas vert horiz height width a
-> Full meas horiz vert width height a
Basic.transpose ((Int, Full meas vert horiz height width a)
 -> (Int, Full meas horiz vert width height a))
-> (Int, Full meas vert horiz height width a)
-> (Int, Full meas horiz vert width height a)
forall a b. (a -> b) -> a -> b
$
         RealOf a
-> Full meas horiz vert width height a
-> Full meas vert horiz width width a
-> (Int, Full meas vert horiz height width a)
forall meas vert horiz height width nrhs a.
(Measure meas, C vert, C horiz, C height, Eq height, C width,
 C nrhs, Floating a) =>
RealOf a
-> Full meas horiz vert height width a
-> Full meas vert horiz height nrhs a
-> (Int, Full meas vert horiz width nrhs a)
leastSquaresMinimumNormRCond RealOf a
rcond (Full meas vert horiz height width a
-> Full meas horiz vert width height a
forall meas vert horiz height width a.
(Measure meas, C vert, C horiz) =>
Full meas vert horiz height width a
-> Full meas horiz vert width height a
Basic.transpose Full meas vert horiz height width a
a) (Full meas vert horiz width width a
 -> (Int, Full meas vert horiz height width a))
-> Full meas vert horiz width width a
-> (Int, Full meas vert horiz height width a)
forall a b. (a -> b) -> a -> b
$
         Square width a -> Full meas vert horiz 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
Square.toFull (Square width a -> Full meas vert horiz width width a)
-> Square width a -> Full meas vert horiz width width a
forall a b. (a -> b) -> a -> b
$ width -> Square width a
forall sh a. (C sh, Floating a) => sh -> Square sh a
Square.identity (width -> Square width a) -> width -> Square width a
forall a b. (a -> b) -> a -> b
$
         Full meas vert horiz height width -> width
forall meas vert horiz height width.
(Measure meas, C vert, C horiz) =>
Full meas vert horiz height width -> width
Layout.fullWidth (Full meas vert horiz height width -> width)
-> Full meas vert horiz height width -> width
forall a b. (a -> b) -> a -> b
$ Full meas vert horiz height width a
-> Full meas vert horiz height width
forall sh a. Array sh a -> sh
Array.shape Full meas vert horiz height width a
a
      Right Wide height width a
_ ->
         RealOf a
-> Full meas vert horiz height width a
-> Full meas horiz vert height height a
-> (Int, Full meas horiz vert width height a)
forall meas vert horiz height width nrhs a.
(Measure meas, C vert, C horiz, C height, Eq height, C width,
 C nrhs, Floating a) =>
RealOf a
-> Full meas horiz vert height width a
-> Full meas vert horiz height nrhs a
-> (Int, Full meas vert horiz width nrhs a)
leastSquaresMinimumNormRCond RealOf a
rcond Full meas vert horiz height width a
a (Full meas horiz vert height height a
 -> (Int, Full meas horiz vert width height a))
-> Full meas horiz vert height height a
-> (Int, Full meas horiz vert width height a)
forall a b. (a -> b) -> a -> b
$
         Square height a -> Full meas horiz vert 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
Square.toFull (Square height a -> Full meas horiz vert height height a)
-> Square height a -> Full meas horiz vert height height a
forall a b. (a -> b) -> a -> b
$ height -> Square height a
forall sh a. (C sh, Floating a) => sh -> Square sh a
Square.identity (height -> Square height a) -> height -> Square height a
forall a b. (a -> b) -> a -> b
$
         Full meas vert horiz height width -> height
forall meas vert horiz height width.
(Measure meas, C vert, C horiz) =>
Full meas vert horiz height width -> height
Layout.fullHeight (Full meas vert horiz height width -> height)
-> Full meas vert horiz height width -> height
forall a b. (a -> b) -> a -> b
$ Full meas vert horiz height width a
-> Full meas vert horiz height width
forall sh a. Array sh a -> sh
Array.shape Full meas vert horiz height width a
a


leastSquaresConstraint ::
   (Shape.C height, Eq height,
    Shape.C width, Eq width,
    Shape.C constraints, Eq constraints, Class.Floating a) =>
   General height width a -> Vector height a ->
   Wide constraints width a -> Vector constraints a ->
   Vector width a
leastSquaresConstraint :: forall height width constraints a.
(C height, Eq height, C width, Eq width, C constraints,
 Eq constraints, Floating a) =>
General height width a
-> Vector height a
-> Wide constraints width a
-> Vector constraints a
-> Vector width a
leastSquaresConstraint
   (Array (Layout.Full Order
orderA Extent Size Big Big height width
extentA) ForeignPtr a
a) Vector height a
c
   (Array (Layout.Full Order
orderB Extent Size Small Big constraints width
extentB) ForeignPtr a
b) Vector constraints a
d =

 let sameShape :: [Char] -> a -> a -> a
sameShape [Char]
name a
shape0 a
shape1 =
      if a
shape0 a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
shape1
         then a
shape0
         else [Char] -> a
forall a. HasCallStack => [Char] -> a
error ([Char] -> a) -> [Char] -> a
forall a b. (a -> b) -> a -> b
$ [Char]
"leastSquaresConstraint: " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Char]
name [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Char]
" shapes mismatch"
     width :: width
width = [Char] -> width -> width -> width
forall {a}. Eq a => [Char] -> a -> a -> a
sameShape [Char]
"width" (Extent Size Big Big height width -> width
forall meas vert horiz height width.
(Measure meas, C vert, C horiz) =>
Extent meas vert horiz height width -> width
Extent.width Extent Size Big Big height width
extentA) (Extent Size Small Big constraints width -> width
forall meas vert horiz height width.
(Measure meas, C vert, C horiz) =>
Extent meas vert horiz height width -> width
Extent.width Extent Size Small Big constraints width
extentB)
 in
   width -> (Ptr a -> IO ()) -> Array width a
forall sh a.
(C sh, Storable a) =>
sh -> (Ptr a -> IO ()) -> Array sh a
Array.unsafeCreate width
width ((Ptr a -> IO ()) -> Array width a)
-> (Ptr a -> IO ()) -> Array width a
forall a b. (a -> b) -> a -> b
$ \Ptr a
xPtr -> do

   let height :: height
height = [Char] -> height -> height -> height
forall {a}. Eq a => [Char] -> a -> a -> a
sameShape [Char]
"height" (Extent Size Big Big height width -> height
forall meas vert horiz height width.
(Measure meas, C vert, C horiz) =>
Extent meas vert horiz height width -> height
Extent.height Extent Size Big Big height width
extentA) (Vector height a -> height
forall sh a. Array sh a -> sh
Array.shape Vector height a
c)
   let constraints :: constraints
constraints =
         [Char] -> constraints -> constraints -> constraints
forall {a}. Eq a => [Char] -> a -> a -> a
sameShape [Char]
"constraints" (Extent Size Small Big constraints width -> constraints
forall meas vert horiz height width.
(Measure meas, C vert, C horiz) =>
Extent meas vert horiz height width -> height
Extent.height Extent Size Small Big constraints width
extentB) (Vector constraints a -> constraints
forall sh a. Array sh a -> sh
Array.shape Vector constraints a
d)
   let m :: Int
m = height -> Int
forall sh. C sh => sh -> Int
Shape.size height
height
   let n :: Int
n = width -> Int
forall sh. C sh => sh -> Int
Shape.size width
width
   let p :: Int
p = constraints -> Int
forall sh. C sh => sh -> Int
Shape.size constraints
constraints
   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
mPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
m
      Ptr CInt
nPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
      Ptr CInt
pPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
p
      Ptr a
aPtr <- Order -> Int -> Int -> ForeignPtr a -> ContT () IO (Ptr a)
forall a r.
Floating a =>
Order -> Int -> Int -> ForeignPtr a -> ContT r IO (Ptr a)
copyToColumnMajorTemp Order
orderA Int
m Int
n ForeignPtr a
a
      Ptr CInt
ldaPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
m
      Ptr a
bPtr <- Order -> Int -> Int -> ForeignPtr a -> ContT () IO (Ptr a)
forall a r.
Floating a =>
Order -> Int -> Int -> ForeignPtr a -> ContT r IO (Ptr a)
copyToColumnMajorTemp Order
orderB Int
p Int
n ForeignPtr a
b
      Ptr CInt
ldbPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
p
      Ptr a
cPtr <- Int -> ForeignPtr a -> ContT () IO (Ptr a)
forall a r. Storable a => Int -> ForeignPtr a -> ContT r IO (Ptr a)
copyToTemp Int
m (Vector height a -> ForeignPtr a
forall sh a. Array sh a -> ForeignPtr a
Array.buffer Vector height a
c)
      Ptr a
dPtr <- Int -> ForeignPtr a -> ContT () IO (Ptr a)
forall a r. Storable a => Int -> ForeignPtr a -> ContT r IO (Ptr a)
copyToTemp Int
p (Vector constraints a -> ForeignPtr a
forall sh a. Array sh a -> ForeignPtr a
Array.buffer Vector constraints a
d)
      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
$ [Char]
-> [Char] -> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
forall a.
Floating a =>
[Char]
-> [Char] -> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
withAutoWorkspaceInfo [Char]
rankMsg [Char]
"gglse" ((Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ())
-> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$
         Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> IO ()
forall a.
Floating a =>
Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> IO ()
LapackGen.gglse
            Ptr CInt
mPtr Ptr CInt
nPtr Ptr CInt
pPtr Ptr a
aPtr Ptr CInt
ldaPtr Ptr a
bPtr Ptr CInt
ldbPtr Ptr a
cPtr Ptr a
dPtr Ptr a
xPtr

gaussMarkovLinearModel ::
   (Shape.C height, Eq height,
    Shape.C width, Eq width,
    Shape.C opt, Eq opt, Class.Floating a) =>
   Tall height width a -> General height opt a -> Vector height a ->
   (Vector width a, Vector opt a)
gaussMarkovLinearModel :: forall height width opt a.
(C height, Eq height, C width, Eq width, C opt, Eq opt,
 Floating a) =>
Tall height width a
-> General height opt a
-> Vector height a
-> (Vector width a, Vector opt a)
gaussMarkovLinearModel
   (Array (Layout.Full Order
orderA Extent Size Big Small height width
extentA) ForeignPtr a
a)
   (Array (Layout.Full Order
orderB Extent Size Big Big height opt
extentB) ForeignPtr a
b) Vector height a
d =

   let width :: width
width = Extent Size Big Small height width -> width
forall meas vert horiz height width.
(Measure meas, C vert, C horiz) =>
Extent meas vert horiz height width -> width
Extent.width Extent Size Big Small height width
extentA in
   let opt :: opt
opt = Extent Size Big Big height opt -> opt
forall meas vert horiz height width.
(Measure meas, C vert, C horiz) =>
Extent meas vert horiz height width -> width
Extent.width Extent Size Big Big height opt
extentB in
   width
-> (Int -> Ptr a -> IO (Vector opt a))
-> (Array width a, Vector opt a)
forall sh a b.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO b) -> (Array sh a, b)
Array.unsafeCreateWithSizeAndResult width
width ((Int -> Ptr a -> IO (Vector opt a))
 -> (Array width a, Vector opt a))
-> (Int -> Ptr a -> IO (Vector opt a))
-> (Array width a, Vector opt a)
forall a b. (a -> b) -> a -> b
$ \Int
m Ptr a
xPtr -> do
   opt -> (Int -> Ptr a -> IO ()) -> IO (Vector opt a)
forall (m :: * -> *) sh a.
(PrimMonad m, C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> m (Array sh a)
ArrayIO.unsafeCreateWithSize opt
opt ((Int -> Ptr a -> IO ()) -> IO (Vector opt a))
-> (Int -> Ptr a -> IO ()) -> IO (Vector opt a)
forall a b. (a -> b) -> a -> b
$ \Int
p Ptr a
yPtr -> do

   let sameHeight :: a -> a -> a
sameHeight a
shape0 a
shape1 =
         if a
shape0 a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
shape1
            then a
shape0
            else [Char] -> a
forall a. HasCallStack => [Char] -> a
error ([Char] -> a) -> [Char] -> a
forall a b. (a -> b) -> a -> b
$ [Char]
"gaussMarkovLinearModel: height shapes mismatch"
       height :: height
height =
         height -> height -> height
forall {a}. Eq a => a -> a -> a
sameHeight (Vector height a -> height
forall sh a. Array sh a -> sh
Array.shape Vector height a
d) (height -> height) -> height -> height
forall a b. (a -> b) -> a -> b
$
         height -> height -> height
forall {a}. Eq a => a -> a -> a
sameHeight (Extent Size Big Small height width -> height
forall meas vert horiz height width.
(Measure meas, C vert, C horiz) =>
Extent meas vert horiz height width -> height
Extent.height Extent Size Big Small height width
extentA) (Extent Size Big Big height opt -> height
forall meas vert horiz height width.
(Measure meas, C vert, C horiz) =>
Extent meas vert horiz height width -> height
Extent.height Extent Size Big Big height opt
extentB)
   let n :: Int
n = height -> Int
forall sh. C sh => sh -> Int
Shape.size height
height
   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
mPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
m
      Ptr CInt
nPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
      Ptr CInt
pPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
p
      Ptr a
aPtr <- Order -> Int -> Int -> ForeignPtr a -> ContT () IO (Ptr a)
forall a r.
Floating a =>
Order -> Int -> Int -> ForeignPtr a -> ContT r IO (Ptr a)
copyToColumnMajorTemp Order
orderA Int
n Int
m ForeignPtr a
a
      Ptr CInt
ldaPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
n
      Ptr a
bPtr <- Order -> Int -> Int -> ForeignPtr a -> ContT () IO (Ptr a)
forall a r.
Floating a =>
Order -> Int -> Int -> ForeignPtr a -> ContT r IO (Ptr a)
copyToColumnMajorTemp Order
orderB Int
n Int
p ForeignPtr a
b
      Ptr CInt
ldbPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
n
      Ptr a
dPtr <- Int -> ForeignPtr a -> ContT () IO (Ptr a)
forall a r. Storable a => Int -> ForeignPtr a -> ContT r IO (Ptr a)
copyToTemp Int
n (Vector height a -> ForeignPtr a
forall sh a. Array sh a -> ForeignPtr a
Array.buffer Vector height a
d)
      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
$ [Char]
-> [Char] -> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
forall a.
Floating a =>
[Char]
-> [Char] -> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
withAutoWorkspaceInfo [Char]
rankMsg [Char]
"ggglm" ((Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ())
-> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$
         Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> IO ()
forall a.
Floating a =>
Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> IO ()
LapackGen.ggglm
            Ptr CInt
nPtr Ptr CInt
mPtr Ptr CInt
pPtr Ptr a
aPtr Ptr CInt
ldaPtr Ptr a
bPtr Ptr CInt
ldbPtr Ptr a
dPtr Ptr a
xPtr Ptr a
yPtr


determinantAbsolute ::
   (Extent.Measure meas, Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width,
    Class.Floating a) =>
   Full meas vert horiz height width a -> RealOf a
determinantAbsolute :: 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 -> RealOf a
determinantAbsolute =
   a -> RealOf a
forall a. Floating a => a -> RealOf a
absolute (a -> RealOf a)
-> (Full meas vert horiz height width a -> a)
-> Full meas vert horiz height width a
-> RealOf a
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   (Full Size Big Small height width a -> a)
-> (Wide height width a -> a)
-> Either
     (Full Size Big Small height width a) (Wide height width a)
-> a
forall a c b. (a -> c) -> (b -> c) -> Either a b -> c
either (Householder Size Big Small height width a -> a
forall meas vert height width a.
(Measure meas, C vert, C height, C width, Floating a) =>
Householder meas vert Small height width a -> a
HH.determinantR (Householder Size Big Small height width a -> a)
-> (Full Size Big Small height width a
    -> Householder Size Big Small height width a)
-> Full Size Big Small height width a
-> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Full Size Big Small height width a
-> Householder Size Big Small height width 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
-> Householder meas vert horiz height width a
HH.fromMatrix) (a -> Wide height width a -> a
forall a b. a -> b -> a
const a
forall a. Floating a => a
zero) (Either (Full Size Big Small height width a) (Wide height width a)
 -> a)
-> (Full meas vert horiz height width a
    -> Either
         (Full Size Big Small height width a) (Wide height width a))
-> Full meas vert horiz height width a
-> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   Full meas vert horiz height width a
-> Either
     (Full Size Big Small height width a) (Wide height width a)
forall meas vert horiz height width a.
(Measure meas, C vert, C horiz, C height, C width) =>
Full meas vert horiz height width a
-> Either (Tall height width a) (Wide height width a)
Basic.caseTallWide


complement ::
   (Shape.C height, Shape.C width, Class.Floating a) =>
   Tall height width a -> Tall height ShapeInt a
complement :: forall height width a.
(C height, C width, Floating a) =>
Tall height width a -> Tall height ShapeInt a
complement = Tall height width a -> Tall height ShapeInt a
forall height width a.
(C height, C width, Floating a) =>
Tall height width a -> Tall height ShapeInt a
extractComplement (Tall height width a -> Tall height ShapeInt a)
-> (Tall height width a -> Tall height width a)
-> Tall height width a
-> Tall height ShapeInt a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Tall height width a -> Tall height width 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
-> Householder meas vert horiz height width a
HH.fromMatrix

extractComplement ::
   (Shape.C height, Shape.C width, Class.Floating a) =>
   HH.Tall height width a -> Tall height ShapeInt a
extractComplement :: forall height width a.
(C height, C width, Floating a) =>
Tall height width a -> Tall height ShapeInt a
extractComplement Tall height width a
qr =
   Int
-> Full Size Big Small height ShapeInt a
-> Full Size Big Small height ShapeInt a
forall horiz height a.
(C horiz, C height, Floating a) =>
Int
-> Full Size Big horiz height ShapeInt a
-> Full Size Big horiz height ShapeInt a
Basic.dropColumns
      (width -> Int
forall sh. C sh => sh -> Int
Shape.size (width -> Int) -> width -> Int
forall a b. (a -> b) -> a -> b
$ Split Reflector Size Big Small height width -> width
forall meas vert horiz lower height width.
(Measure meas, C vert, C horiz) =>
Split lower meas vert horiz height width -> width
Layout.splitWidth (Split Reflector Size Big Small height width -> width)
-> Split Reflector Size Big Small height width -> width
forall a b. (a -> b) -> a -> b
$ Array (Split Reflector Size Big Small height width) a
-> Split Reflector Size Big Small height width
forall sh a. Array sh a -> sh
Array.shape (Array (Split Reflector Size Big Small height width) a
 -> Split Reflector Size Big Small height width)
-> Array (Split Reflector Size Big Small height width) a
-> Split Reflector Size Big Small height width
forall a b. (a -> b) -> a -> b
$ Tall height width a
-> Array (Split Reflector Size Big Small height width) a
forall xl xu lower upper meas vert horiz height width a.
Matrix Hh xl xu lower upper meas vert horiz height width a
-> SplitArray meas vert horiz height width a
HH.split_ Tall height width a
qr) (Full Size Big Small height ShapeInt a
 -> Full Size Big Small height ShapeInt a)
-> Full Size Big Small height ShapeInt a
-> Full Size Big Small height ShapeInt a
forall a b. (a -> b) -> a -> b
$
   (height -> ShapeInt)
-> Full Size Big Small height height a
-> Full Size Big Small height ShapeInt 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 (Int -> ShapeInt
shapeInt (Int -> ShapeInt) -> (height -> Int) -> height -> ShapeInt
forall b c a. (b -> c) -> (a -> b) -> a -> c
. height -> Int
forall sh. C sh => sh -> Int
Shape.size) (Full Size Big Small height height a
 -> Full Size Big Small height ShapeInt a)
-> Full Size Big Small height height a
-> Full Size Big Small height ShapeInt a
forall a b. (a -> b) -> a -> b
$ Square height a -> Full Size Big Small 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
Square.toFull (Square height a -> Full Size Big Small height height a)
-> Square height a -> Full Size Big Small height height a
forall a b. (a -> b) -> a -> b
$
   Tall height width a -> Square height a
forall meas vert horiz height width a.
(Measure meas, C vert, C horiz, C height, C width, Floating a) =>
Householder meas vert horiz height width a -> Square height a
HH.extractQ Tall height width a
qr