packages feed

lapack 0.0 → 0.1

raw patch · 22 files changed

+3980/−735 lines, 22 filesdep ~comfort-arraydep ~netlib-ffi

Dependency ranges changed: comfort-array, netlib-ffi

Files

lapack.cabal view
@@ -1,5 +1,5 @@ Name:             lapack-Version:          0.0+Version:          0.1 License:          BSD3 License-File:     LICENSE Author:           Henning Thielemann <haskell@henning-thielemann.de>@@ -9,6 +9,9 @@ Synopsis:         Numerical Linear Algebra using LAPACK Description:   This is a high-level interface to LAPACK.+  It provides solvers for simultaneous linear equations,+  linear least-squares problems, eigenvalue and singular value problems+  for matrices with certain kinds of structures.   .   Features:   .@@ -34,7 +37,7 @@ Build-Type:       Simple  Source-Repository this-  Tag:         0.0+  Tag:         0.1   Type:        darcs   Location:    http://hub.darcs.net/thielema/lapack/ @@ -46,8 +49,8 @@   Build-Depends:     lapack-ffi >=0.0.1 && <0.1,     blas-ffi >=0.0 && <0.1,-    netlib-ffi >=0.0.1 && <0.1,-    comfort-array >=0.0 && <0.1,+    netlib-ffi >=0.1 && <0.2,+    comfort-array >=0.0.1 && <0.1,     transformers >=0.3 && <0.6,     non-empty >=0.3 && <0.4,     utility-ht >=0.0.10 && <0.1,@@ -58,9 +61,23 @@   Exposed-Modules:     Numeric.LAPACK.Matrix     Numeric.LAPACK.Matrix.Shape+    Numeric.LAPACK.Matrix.Square+    Numeric.LAPACK.Matrix.Hermitian+    Numeric.LAPACK.Matrix.Triangular     Numeric.LAPACK.Vector-    Numeric.LAPACK.LinearSystem+    Numeric.LAPACK.Orthogonal+    Numeric.LAPACK.Linear.General+    Numeric.LAPACK.Linear.HermitianPositiveDefinite+    Numeric.LAPACK.Linear.Hermitian+    Numeric.LAPACK.Linear.Triangular+    Numeric.LAPACK.Eigen.General+    Numeric.LAPACK.Eigen.Hermitian+    Numeric.LAPACK.Eigen.Triangular+    Numeric.LAPACK.Singular   Other-Modules:+    Numeric.LAPACK.Matrix.Triangular.Private     Numeric.LAPACK.Matrix.Shape.Private+    Numeric.LAPACK.Matrix.Multiply+    Numeric.LAPACK.Matrix.Private     Numeric.LAPACK.Private     Numeric.LAPACK.Format
+ src/Numeric/LAPACK/Eigen/General.hs view
@@ -0,0 +1,319 @@+{-# LANGUAGE TypeFamilies #-}+module Numeric.LAPACK.Eigen.General (+   values,+   schur,+   decompose,+   ComplexOf,+   ) where++import Numeric.LAPACK.Matrix.Square (Square)++import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape+import Numeric.LAPACK.Matrix.Shape.Private (Order(RowMajor,ColumnMajor))+import Numeric.LAPACK.Vector (Vector)+import Numeric.LAPACK.Private+         (ComplexOf, RealOf, zero, withAutoWorkspaceInfo,+          copyToTemp, copyToColumnMajor, allocArray)++import qualified Numeric.LAPACK.FFI.Complex as LapackComplex+import qualified Numeric.LAPACK.FFI.Real as LapackReal+import qualified Numeric.BLAS.FFI.Complex as BlasComplex+import qualified Numeric.BLAS.FFI.Real as BlasReal+import qualified Numeric.Netlib.Utility as Call+import qualified Numeric.Netlib.Class as Class++import qualified Data.Array.Comfort.Storable.Internal as Array+import qualified Data.Array.Comfort.Shape as Shape+import Data.Array.Comfort.Storable.Internal (Array(Array))++import System.IO.Unsafe (unsafePerformIO)++import Foreign.Marshal.Array (advancePtr, peekArray)+import Foreign.C.Types (CInt, CChar)+import Foreign.ForeignPtr (withForeignPtr)+import Foreign.Ptr (Ptr, nullPtr, nullFunPtr, castPtr)+import Foreign.Storable (Storable)++import Control.Monad.Trans.Cont (ContT(ContT), evalContT)+import Control.Monad.IO.Class (liftIO)++import Data.Complex (Complex)+++values ::+   (Shape.C sh, Class.Floating a) =>+   Square sh a -> Vector sh (ComplexOf a)+values =+   getValues $+   Class.switchFloating+      (Values valuesAux) (Values valuesAux)+      (Values valuesAux) (Values valuesAux)++type Values_ sh a = Square sh a -> Vector sh (ComplexOf a)++newtype Values sh a = Values {getValues :: Values_ sh a}++valuesAux ::+   (Shape.C sh, Class.Floating a, RealOf a ~ ar, Storable ar) =>+   Values_ sh a+valuesAux (Array (MatrixShape.Square _order size) a) =+      Array.unsafeCreateWithSize size $ \n wPtr -> do+   let lda = n+   evalContT $ do+      jobvsPtr <- Call.char 'N'+      sortPtr <- Call.char 'N'+      aPtr <- copyToTemp (n*n) a+      ldaPtr <- Call.cint lda+      sdimPtr <- Call.alloca+      let vsPtr = nullPtr+      ldvsPtr <- Call.cint n+      let bworkPtr = nullPtr+      liftIO $ withAutoWorkspaceInfo "gees" $ \workPtr lworkPtr infoPtr ->+         gees+            jobvsPtr sortPtr n aPtr ldaPtr sdimPtr+            wPtr vsPtr ldvsPtr workPtr lworkPtr bworkPtr infoPtr+++{- |+If @(q,r) = schur a@, then @a = q \<#\> r \<#\> adjoint q@,+where @q@ is unitary (orthogonal)+and @r@ is a right-upper triangular matrix for complex @a@+and a 1x1-or-2x2-block upper triangular matrix for real @a@.+With @getDiagonal r@ you get all eigenvalues of @a@ if @a@ is complex+and the real parts of the eigenvalues if @a@ is real.+Complex conjugated eigenvalues of a real matrix @a@+are encoded as 2x2 blocks along the diagonal.+-}+schur ::+   (Shape.C sh, Class.Floating a) =>+   Square sh a -> (Square sh a, Square sh a)+schur =+   getSchur $+   Class.switchFloating+      (Schur schurAux) (Schur schurAux)+      (Schur schurAux) (Schur schurAux)++type Schur_ sh a = Square sh a -> (Square sh a, Square sh a)++newtype Schur sh a = Schur {getSchur :: Schur_ sh a}++schurAux ::+   (Shape.C sh, Class.Floating a, RealOf a ~ ar, Storable ar) =>+   Schur_ sh a+schurAux (Array (MatrixShape.Square order size) a) = unsafePerformIO $ do+   let n = Shape.size size+   let lda = n+   let sh = MatrixShape.Square ColumnMajor size+   evalContT $ do+      jobvsPtr <- Call.char 'V'+      sortPtr <- Call.char 'N'+      aPtr <- ContT $ withForeignPtr a+      (s,sPtr) <- allocArray sh+      liftIO $ copyToColumnMajor order n n aPtr sPtr+      ldaPtr <- Call.cint lda+      sdimPtr <- Call.alloca+      wPtr <- Call.allocaArray n+      (vs,vsPtr) <- allocArray sh+      ldvsPtr <- Call.cint n+      let bworkPtr = nullPtr+      liftIO $ withAutoWorkspaceInfo "gees" $ \workPtr lworkPtr infoPtr ->+         gees+            jobvsPtr sortPtr n sPtr ldaPtr sdimPtr+            wPtr vsPtr ldvsPtr workPtr lworkPtr bworkPtr infoPtr+      return (vs, s)++++type GEES_ ar a =+   Ptr CChar -> Ptr CChar -> Int -> Ptr a -> Ptr CInt ->+   Ptr CInt -> Ptr (Complex ar) -> Ptr a -> Ptr CInt ->+   Ptr a -> Ptr CInt -> Ptr Bool -> Ptr CInt -> IO ()++newtype GEES a = GEES {getGEES :: GEES_ (RealOf a) a}++gees :: Class.Floating a => GEES_ (RealOf a) a+gees =+   getGEES $+   Class.switchFloating+      (GEES geesReal) (GEES geesReal) (GEES geesComplex) (GEES geesComplex)++geesReal :: Class.Real a => GEES_ a a+geesReal+      jobvsPtr sortPtr n aPtr ldaPtr sdimPtr+      wPtr vsPtr ldvsPtr workPtr lworkPtr bworkPtr infoPtr =+   evalContT $ do+      let selectPtr = nullFunPtr+      nPtr <- Call.cint n+      wrPtr <- Call.allocaArray n+      wiPtr <- Call.allocaArray n+      liftIO $+         LapackReal.gees+            jobvsPtr sortPtr selectPtr nPtr aPtr ldaPtr sdimPtr+            wrPtr wiPtr vsPtr ldvsPtr workPtr lworkPtr bworkPtr infoPtr+      liftIO $ zipComplex n wrPtr wiPtr wPtr++geesComplex :: Class.Real a => GEES_ a (Complex a)+geesComplex+      jobvsPtr sortPtr n aPtr ldaPtr sdimPtr+      wPtr vsPtr ldvsPtr workPtr lworkPtr bworkPtr infoPtr =+   evalContT $ do+      let selectPtr = nullFunPtr+      nPtr <- Call.cint n+      rworkPtr <- Call.allocaArray n+      liftIO $+         LapackComplex.gees+            jobvsPtr sortPtr selectPtr nPtr aPtr ldaPtr sdimPtr+            wPtr vsPtr ldvsPtr workPtr lworkPtr rworkPtr bworkPtr infoPtr++++{- |+@(vr,d,vl) = Eigen.decompose a@++Counterintuitively, @vr@ contains the right eigenvectors+and @vl@ contains the left eigenvectors as columns.+The idea is to provide a decomposition of @a@.+If @a@ is diagonalizable, then @vr@ and @vl@ are almost inverse to each other.+More precisely, @adjoint vl \<#\> vr@ is a diagonal matrix.+This is because all eigenvectors are normalized to Euclidean norm 1.+With the following scaling, the decomposition becomes perfect:++> let scal = Array.map recip $ getDiagonal $ adjoint vl <#> vr+> a == vr <#> diagonal d <#> diagonal scal <#> adjoint vl++If @a@ is non-diagonalizable then some columns of @vr@ and @vl@ are left zero+and the above property does not hold.+-}+decompose ::+   (Shape.C sh, Class.Floating a) =>+   Square sh a ->+   (Square sh (ComplexOf a),+    Vector sh (ComplexOf a),+    Square sh (ComplexOf a))+decompose =+   getDecompose $+   Class.switchFloating+      (Decompose decomposeReal)+      (Decompose decomposeReal)+      (Decompose decomposeComplex)+      (Decompose decomposeComplex)++newtype Decompose sh a =+   Decompose {+      getDecompose ::+         Square sh a ->+         (Square sh (ComplexOf a),+          Vector sh (ComplexOf a),+          Square sh (ComplexOf a))+   }++decomposeReal ::+   (Shape.C sh, Class.Real a) =>+   Square sh a ->+   (Square sh (Complex a), Vector sh (Complex a), Square sh (Complex a))+decomposeReal (Array (MatrixShape.Square order size) a) =+      unsafePerformIO $ do+   let n = Shape.size size+   let lda = n+   evalContT $ do+      jobvlPtr <- Call.char 'V'+      jobvrPtr <- Call.char 'V'+      nPtr <- Call.cint n+      aPtr <- copyToTemp (n*n) a+      ldaPtr <- Call.cint lda+      wrPtr <- Call.allocaArray n+      wiPtr <- Call.allocaArray n+      vlPtr <- Call.allocaArray (n*n)+      ldvlPtr <- Call.cint n+      vrPtr <- Call.allocaArray (n*n)+      ldvrPtr <- Call.cint n+      liftIO $ withAutoWorkspaceInfo "geev" $+         LapackReal.geev+            jobvlPtr jobvrPtr nPtr aPtr ldaPtr+            wrPtr wiPtr vlPtr ldvlPtr vrPtr ldvrPtr+      (w,wPtr) <- allocArray size+      liftIO $ zipComplex n wrPtr wiPtr wPtr+      let sh = MatrixShape.Square ColumnMajor size+      (vlc,vlcPtr) <- allocArray sh+      (vrc,vrcPtr) <- allocArray sh+      liftIO $ eigenvectorsToComplex n wiPtr vlPtr vlcPtr+      liftIO $ eigenvectorsToComplex n wiPtr vrPtr vrcPtr+      return $+         case order of+            RowMajor -> (vlc, w, vrc)+            ColumnMajor -> (vrc, w, vlc)++eigenvectorsToComplex ::+   (Eq a, Class.Real a) =>+   Int -> Ptr a -> Ptr a -> Ptr (Complex a) -> IO ()+eigenvectorsToComplex n wiPtr vPtr vcPtr = evalContT $ do+   nPtr <- Call.cint n+   zeroPtr <- Call.real zero+   inc0Ptr <- Call.cint 0+   inc1Ptr <- Call.cint 1+   inc2Ptr <- Call.cint 2+   liftIO $ do+      let go _ _ [] = return ()+          go xPtr yPtr (False:wi) = do+            let yrPtr = castPtr yPtr+            let yiPtr = advancePtr yrPtr 1+            BlasReal.copy nPtr xPtr    inc1Ptr yrPtr inc2Ptr+            BlasReal.copy nPtr zeroPtr inc0Ptr yiPtr inc2Ptr+            go (advancePtr xPtr n) (advancePtr yPtr n) wi+          go xPtr yPtr (True:True:wi) = do+            let xrPtr = xPtr+            let xiPtr = advancePtr xPtr n+            let yrPtr = castPtr yPtr+            let yiPtr = advancePtr yrPtr 1+            let y1Ptr = advancePtr yPtr n+            BlasReal.copy nPtr xrPtr inc1Ptr yrPtr inc2Ptr+            BlasReal.copy nPtr xiPtr inc1Ptr yiPtr inc2Ptr+            BlasComplex.copy nPtr yPtr inc1Ptr y1Ptr inc1Ptr+            LapackComplex.lacgv nPtr y1Ptr inc1Ptr+            go (advancePtr xPtr (2*n)) (advancePtr yPtr (2*n)) wi+          go _xPtr _yPtr wi =+            error $ "eigenvectorToComplex: invalid non-real pattern " ++ show wi+      go vPtr vcPtr . map (zero/=) =<< peekArray n wiPtr++decomposeComplex ::+   (Shape.C sh, Class.Real a) =>+   Square sh (Complex a) ->+   (Square sh (Complex a), Vector sh (Complex a), Square sh (Complex a))+decomposeComplex (Array (MatrixShape.Square order size) a) =+      unsafePerformIO $ do+   let n = Shape.size size+   let lda = n+   evalContT $ do+      jobvlPtr <- Call.char 'V'+      jobvrPtr <- Call.char 'V'+      nPtr <- Call.cint n+      aPtr <- copyToTemp (n*n) a+      ldaPtr <- Call.cint lda+      (w,wPtr) <- allocArray size+      let sh = MatrixShape.Square ColumnMajor size+      (vl,vlPtr) <- allocArray sh+      ldvlPtr <- Call.cint n+      (vr,vrPtr) <- allocArray sh+      ldvrPtr <- Call.cint n+      rworkPtr <- Call.allocaArray (2*n)+      liftIO $ withAutoWorkspaceInfo "geev" $ \workPtr lworkPtr infoPtr ->+         LapackComplex.geev+            jobvlPtr jobvrPtr nPtr aPtr ldaPtr+            wPtr vlPtr ldvlPtr vrPtr ldvrPtr+            workPtr lworkPtr rworkPtr infoPtr+      return $+         case order of+            RowMajor -> (vl, w, vr)+            ColumnMajor -> (vr, w, vl)+++zipComplex ::+   (Class.Real a) => Int -> Ptr a -> Ptr a -> Ptr (Complex a) -> IO ()+zipComplex n vr vi vc =+   evalContT $ do+      nPtr <- Call.cint n+      incxPtr <- Call.cint 1+      incyPtr <- Call.cint 2+      let yPtr = castPtr vc+      liftIO $ BlasReal.copy nPtr vr incxPtr yPtr incyPtr+      liftIO $ BlasReal.copy nPtr vi incxPtr (advancePtr yPtr 1) incyPtr
+ src/Numeric/LAPACK/Eigen/Hermitian.hs view
@@ -0,0 +1,142 @@+{-# LANGUAGE TypeFamilies #-}+module Numeric.LAPACK.Eigen.Hermitian (+   values,+   decompose,+   ) where++import Numeric.LAPACK.Matrix.Hermitian (Hermitian)+import Numeric.LAPACK.Matrix.Square (Square)++import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape+import Numeric.LAPACK.Matrix.Shape.Private+         (Order(ColumnMajor), uploFromOrder, triangleSize)+import Numeric.LAPACK.Vector (Vector)+import Numeric.LAPACK.Private (RealOf, copyToTemp, allocArray)++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.Internal as Array+import qualified Data.Array.Comfort.Shape as Shape+import Data.Array.Comfort.Storable.Internal (Array(Array))++import System.IO.Unsafe (unsafePerformIO)++import Foreign.Marshal.Alloc (alloca)+import Foreign.C.Types (CInt, CChar)+import Foreign.Ptr (Ptr, nullPtr)+import Foreign.Storable (Storable, peek)++import Control.Monad.Trans.Cont (evalContT)+import Control.Monad.IO.Class (liftIO)+import Control.Applicative ((<$>))++import Text.Printf (printf)++import Data.Complex (Complex)+++values ::+   (Shape.C sh, Class.Floating a) =>+   Hermitian sh a -> Vector sh (RealOf a)+values =+   getValues $+   Class.switchFloating+      (Values valuesAux) (Values valuesAux)+      (Values valuesAux) (Values valuesAux)++newtype Values sh a =+   Values {getValues :: Hermitian sh a -> Vector sh (RealOf a)}++valuesAux ::+   (Shape.C sh, Class.Floating a, RealOf a ~ ar, Storable ar) =>+   Hermitian sh a -> Vector sh ar+valuesAux (Array (MatrixShape.Hermitian order size) a) =+      Array.unsafeCreateWithSize size $ \n wPtr -> do+   evalContT $ do+      jobzPtr <- Call.char 'N'+      uploPtr <- Call.char $ uploFromOrder order+      aPtr <- copyToTemp (triangleSize n) a+      let zPtr = nullPtr+      ldzPtr <- Call.cint (max 1 n)+      liftIO $ withInfo "hpev" $+         hpev jobzPtr uploPtr n aPtr wPtr zPtr ldzPtr+++{- |+For symmetric eigenvalue problems, @Eigen.decompose@ and @schur@ coincide.+-}+decompose ::+   (Shape.C sh, Class.Floating a) =>+   Hermitian sh a -> (Square sh a, Vector sh (RealOf a))+decompose =+   getDecompose $+   Class.switchFloating+      (Decompose decomposeAux) (Decompose decomposeAux)+      (Decompose decomposeAux) (Decompose decomposeAux)++type Decompose_ sh a = Hermitian sh a -> (Square sh a, Vector sh (RealOf a))++newtype Decompose sh a = Decompose {getDecompose :: Decompose_ sh a}++decomposeAux ::+   (Shape.C sh, Class.Floating a, RealOf a ~ ar, Storable ar) =>+   Decompose_ sh a+decomposeAux (Array (MatrixShape.Hermitian order size) a) = unsafePerformIO $ do+   let n = Shape.size size+   let shZ = MatrixShape.Square ColumnMajor size+   evalContT $ do+      jobzPtr <- Call.char 'V'+      uploPtr <- Call.char $ uploFromOrder order+      aPtr <- copyToTemp (triangleSize n) a+      (w,wPtr) <- allocArray size+      (z,zPtr) <- allocArray shZ+      ldzPtr <- Call.cint n+      liftIO $ withInfo "hpev" $+         hpev jobzPtr uploPtr n aPtr wPtr zPtr ldzPtr+      return (z, w)+++withInfo :: String -> (Ptr CInt -> IO ()) -> IO ()+withInfo name computation = alloca $ \infoPtr -> do+   computation infoPtr+   info <- fromIntegral <$> peek infoPtr+   case compare info (0::Int) of+      EQ -> return ()+      LT -> error $ printf "%s: illegal value in %d-th argument" name (-info)+      GT -> error $+         printf "%s: %d off-diagonal elements not converging" name info+++type HPEV_ ar a =+   Ptr CChar -> Ptr CChar -> Int -> Ptr a -> Ptr ar ->+   Ptr a -> Ptr CInt -> Ptr CInt -> IO ()++newtype HPEV a = HPEV {getHPEV :: HPEV_ (RealOf a) a}++hpev :: Class.Floating a => HPEV_ (RealOf a) a+hpev =+   getHPEV $+   Class.switchFloating+      (HPEV spevReal) (HPEV spevReal) (HPEV hpevComplex) (HPEV hpevComplex)++spevReal :: Class.Real a => HPEV_ a a+spevReal jobzPtr uploPtr n apPtr wPtr zPtr ldzPtr infoPtr =+   evalContT $ do+      nPtr <- Call.cint n+      workPtr <- Call.allocaArray (3*n)+      liftIO $+         LapackReal.spev+            jobzPtr uploPtr nPtr apPtr wPtr zPtr ldzPtr workPtr infoPtr++hpevComplex :: Class.Real a => HPEV_ a (Complex a)+hpevComplex jobzPtr uploPtr n apPtr wPtr zPtr ldzPtr infoPtr =+   evalContT $ do+      nPtr <- Call.cint n+      workPtr <- Call.allocaArray (max 1 (2*n-1))+      rworkPtr <- Call.allocaArray (max 1 (3*n-2))+      liftIO $+         LapackComplex.hpev+            jobzPtr uploPtr nPtr apPtr wPtr zPtr ldzPtr workPtr rworkPtr infoPtr
+ src/Numeric/LAPACK/Eigen/Triangular.hs view
@@ -0,0 +1,176 @@+{-# LANGUAGE TypeFamilies #-}+module Numeric.LAPACK.Eigen.Triangular (+   values,+   decompose,+   ) where++import qualified Numeric.LAPACK.Matrix.Triangular as Triangular+import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape+import Numeric.LAPACK.Matrix.Triangular.Private+         (unpackZero, pack, unpackToTemp, fillTriangle,+          forPointers, rowMajorPointers)+import Numeric.LAPACK.Matrix.Triangular (Triangular)+import Numeric.LAPACK.Matrix.Shape.Private+         (Order(ColumnMajor,RowMajor), caseUplo, uploOrder, triangleSize)+import Numeric.LAPACK.Vector (Vector)+import Numeric.LAPACK.Private (zero, lacgv, allocArray)++import qualified Numeric.LAPACK.FFI.Complex as LapackComplex+import qualified Numeric.LAPACK.FFI.Real as LapackReal+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.Shape as Shape+import Data.Array.Comfort.Storable.Internal (Array(Array))++import System.IO.Unsafe (unsafePerformIO)++import Foreign.Marshal.Alloc (alloca)+import Foreign.C.Types (CInt, CChar)+import Foreign.Ptr (Ptr, nullPtr)+import Foreign.Storable (peek)++import Control.Monad.Trans.Cont (evalContT)+import Control.Monad.IO.Class (liftIO)+import Control.Applicative ((<$>))++import Text.Printf (printf)++import Data.Complex (Complex)+++values ::+   (MatrixShape.Uplo uplo, Shape.C sh, Class.Floating a) =>+   Triangular uplo sh a -> Vector sh a+values = Triangular.getDiagonal+++{- |+@(vr,d,vlAdj) = TriEigen.decompose a@++Counterintuitively, @vr@ contains the right eigenvectors as columns+and @vlAdj@ contains the left conjugated eigenvectors as rows.+The idea is to provide a decomposition of @a@.+If @a@ is diagonalizable, then @vr@ and @vlAdj@+are almost inverse to each other.+More precisely, @vlAdj \<#\> vr@ is a diagonal matrix.+This is because the eigenvectors are not normalized.+With the following scaling, the decomposition becomes perfect:++> let scal = Array.map recip $ getDiagonal $ vlAdj <#> vr+> a == vr <#> diagonal d <#> diagonal scal <#> vlAdj++If @a@ is non-diagonalizable+then some columns of @vr@ and corresponding rows of @vlAdj@ are left zero+and the above property does not hold.+-}+decompose ::+   (MatrixShape.Uplo uplo, Shape.C sh, Class.Floating a) =>+   Triangular uplo sh a ->+   (Triangular uplo sh a, Vector sh a, Triangular uplo sh a)+decompose a =+   let (vr,vl) = decomposePlain a+   in  (vr, values a, vl)++decomposePlain ::+   (MatrixShape.Uplo uplo, Shape.C sh, Class.Floating a) =>+   Triangular uplo sh a -> (Triangular uplo sh a, Triangular uplo sh a)+decomposePlain (Array (MatrixShape.Triangular uplo order sh) a) =+      unsafePerformIO $ do+   let n = Shape.size sh+   let n2 = n*n+   let triSize = triangleSize n+   evalContT $ do+      sidePtr <- Call.char 'B'+      howManyPtr <- Call.char 'A'+      let selectPtr = nullPtr+      let unpk =+            case uploOrder uplo order of+               ColumnMajor -> unpackZero ColumnMajor+               RowMajor -> unpackZeroRowMajor+      aPtr <- unpackToTemp unpk n a+      ldaPtr <- Call.cint n+      vlPtr <- Call.allocaArray n2+      vrPtr <- Call.allocaArray n2+      mmPtr <- Call.cint n+      mPtr <- Call.alloca+      liftIO $ withInfo "trevc" $+         trevc sidePtr howManyPtr selectPtr n+            aPtr ldaPtr vlPtr ldaPtr vrPtr ldaPtr mmPtr mPtr+      (vl,vlpPtr) <-+         allocArray $+         MatrixShape.Triangular uplo (uploOrder uplo RowMajor) sh+      (vr,vrpPtr) <-+         allocArray $+         MatrixShape.Triangular uplo (uploOrder uplo ColumnMajor) sh+      sizePtr <- Call.cint triSize+      incPtr <- Call.cint 1+      liftIO $ do+         pack ColumnMajor n vrPtr vrpPtr+         pack RowMajor n vlPtr vlpPtr+         lacgv sizePtr vlpPtr incPtr+      return $ caseUplo uplo (vl,vr) (vr,vl)+++unpackZeroRowMajor :: Class.Floating a => Int -> Ptr a -> Ptr a -> IO ()+unpackZeroRowMajor n packedPtr fullPtr = do+   fillTriangle zero RowMajor n fullPtr+   unpackRowMajor n packedPtr fullPtr++unpackRowMajor :: Class.Floating a => Int -> Ptr a -> Ptr a -> IO ()+unpackRowMajor n packedPtr fullPtr = evalContT $ do+   incxPtr <- Call.cint 1+   incyPtr <- Call.cint n+   liftIO $+      forPointers (rowMajorPointers n fullPtr packedPtr) $+            \nPtr (dstPtr,srcPtr) ->+         BlasGen.copy nPtr srcPtr incxPtr dstPtr incyPtr+++withInfo :: String -> (Ptr CInt -> IO ()) -> IO ()+withInfo name computation = alloca $ \infoPtr -> do+   computation infoPtr+   info <- fromIntegral <$> peek infoPtr+   case compare info (0::Int) of+      EQ -> return ()+      LT -> error $ printf "%s: illegal value in %d-th argument" name (-info)+      GT -> error $+         printf "%s: %d off-diagonal elements not converging" name info+++type TREVC_ a =+   Ptr CChar -> Ptr CChar -> Ptr Bool ->+   Int -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt ->+   Ptr CInt -> Ptr CInt -> Ptr CInt -> IO ()++newtype TREVC a = TREVC {getTREVC :: TREVC_ a}++trevc :: Class.Floating a => TREVC_ a+trevc =+   getTREVC $+   Class.switchFloating+      (TREVC trevcReal) (TREVC trevcReal)+      (TREVC trevcComplex) (TREVC trevcComplex)++trevcReal :: Class.Real a => TREVC_ a+trevcReal sidePtr howmnyPtr selectPtr n+      tPtr ldtPtr vlPtr ldvlPtr vrPtr ldvrPtr mmPtr mPtr infoPtr =+   evalContT $ do+      nPtr <- Call.cint n+      workPtr <- Call.allocaArray (3*n)+      liftIO $+         LapackReal.trevc sidePtr howmnyPtr selectPtr nPtr+            tPtr ldtPtr vlPtr ldvlPtr vrPtr ldvrPtr mmPtr mPtr workPtr infoPtr++trevcComplex :: Class.Real a => TREVC_ (Complex a)+trevcComplex sidePtr howmnyPtr selectPtr n+      tPtr ldtPtr vlPtr ldvlPtr vrPtr ldvrPtr mmPtr mPtr infoPtr =+   evalContT $ do+      nPtr <- Call.cint n+      workPtr <- Call.allocaArray (2*n)+      rworkPtr <- Call.allocaArray n+      liftIO $+         LapackComplex.trevc sidePtr howmnyPtr selectPtr nPtr+            tPtr ldtPtr vlPtr ldvlPtr vrPtr ldvrPtr mmPtr mPtr+            workPtr rworkPtr infoPtr
src/Numeric/LAPACK/Format.hs view
@@ -1,7 +1,13 @@-module Numeric.LAPACK.Format where+module Numeric.LAPACK.Format (+   (##),+   Format(format),+   FormatArray(formatArray),+   ) where  import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape+import qualified Numeric.LAPACK.Matrix.Square as Square import Numeric.LAPACK.Matrix.Shape.Private (Order(RowMajor, ColumnMajor))+import Numeric.LAPACK.Matrix.Private (General)  import qualified Numeric.Netlib.Class as Class @@ -9,11 +15,12 @@ import qualified Data.Array.Comfort.Shape as Shape import Data.Array.Comfort.Storable (Array) -import Foreign.Storable (Storable)- import Text.Printf (PrintfArg, printf)  import qualified Data.List.HT as ListHT+import qualified Data.Complex as Complex+import Data.Monoid (Endo(Endo,appEndo))+import Data.List (mapAccumL, transpose) import Data.Complex (Complex((:+)))  @@ -35,7 +42,7 @@ instance Format Double where    format fmt a = [printf fmt a] -instance (PrintfArg a) => Format (Complex a) where+instance (PrintfArg a, Class.Real a) => Format (Complex a) where    format fmt a = [printfComplex fmt a]  instance (Format a, Format b) => Format (a,b) where@@ -45,15 +52,12 @@    format fmt (a,b,c) =       format fmt a ++ [""] ++ format fmt b ++ [""] ++ format fmt c -instance-   (FormatArray sh, Class.Floating a, Storable a) =>-      Format (Array sh a) where+instance (FormatArray sh, Class.Floating a) => Format (Array sh a) where    format = formatArray   class (Shape.C sh) => FormatArray sh where-   formatArray ::-      (Storable a, Class.Floating a) => String -> Array sh a -> [String]+   formatArray :: (Class.Floating a) => String -> Array sh a -> [String]  instance (Integral i) => FormatArray (Shape.ZeroBased i) where    formatArray fmt m = [unwords $ map (printfFloating fmt) $ Array.toList m]@@ -61,19 +65,20 @@ instance (Integral i) => FormatArray (Shape.OneBased i) where    formatArray fmt m = [unwords $ map (printfFloating fmt) $ Array.toList m] +instance (Shape.C sh) => FormatArray (MatrixShape.Square sh) where+   formatArray fmt = formatGeneral fmt . Square.toGeneral+ instance    (Shape.C height, Shape.C width) =>       FormatArray (MatrixShape.General height width) where    formatArray = formatGeneral  formatGeneral ::-   (Shape.C height, Shape.C width, Storable a, Class.Floating a) =>-   String -> Array (MatrixShape.General height width) a -> [String]+   (Shape.C height, Shape.C width, Class.Floating a) =>+   String -> General height width a -> [String] formatGeneral fmt m =    let MatrixShape.General order height width = Array.shape m-       xss = formatRows fmt order (height,width) $ Array.toList m-       strWidths = columnWidths xss-   in  map (unwords . zipWith (ListHT.padLeft ' ') strWidths) xss+   in  formatAligned $ formatRows fmt order (height,width) $ Array.toList m  instance    (Shape.C height, Shape.C width) =>@@ -81,12 +86,94 @@    formatArray = formatHouseholder  formatHouseholder ::-   (Shape.C height, Shape.C width, Storable a, Class.Floating a) =>+   (Shape.C height, Shape.C width, Class.Floating a) =>    String -> Array (MatrixShape.Householder height width) a -> [String] formatHouseholder fmt m =    let MatrixShape.Householder order height width = Array.shape m-       xss = formatRows fmt order (height,width) $ Array.toList m-       strWidths = columnWidths xss+   in formatSeparateTriangle $+      formatRows fmt order (height,width) $ Array.toList m++instance (Shape.C size) => FormatArray (MatrixShape.Hermitian size) where+   formatArray = formatHermitian++formatHermitian ::+   (Shape.C size, Class.Floating a) =>+   String -> Array (MatrixShape.Hermitian size) a -> [String]+formatHermitian fmt m =+   let MatrixShape.Hermitian order size = Array.shape m+   in  formatSeparateTriangle $+       map (map (printfFloating fmt)) $+       complementTriangle order (Shape.size size) $ Array.toList m++complementTriangle :: (Class.Floating a) => Order -> Int -> [a] -> [[a]]+complementTriangle order n xs =+   let mergeTriangles lower upper =+         zipWith (++) (map (map conjugate . init) lower) upper+   in case order of+         RowMajor ->+            let tri = slice (take n $ iterate pred n) xs+                trans = reverse $ transpose $ map reverse tri+            in  mergeTriangles trans tri+         ColumnMajor ->+            let tri = slice (take n [1..]) xs+            in  mergeTriangles tri (transpose tri)++conjugate :: (Class.Floating a) => a -> a+conjugate =+   appEndo $+   Class.switchFloating+      (Endo id) (Endo id) (Endo Complex.conjugate) (Endo Complex.conjugate)++instance+   (MatrixShape.Uplo uplo, Shape.C size) =>+      FormatArray (MatrixShape.Triangular uplo size) where+   formatArray = formatTriangular++formatTriangular ::+   (MatrixShape.Uplo uplo, Shape.C size, Class.Floating a) =>+   String -> Array (MatrixShape.Triangular uplo size) a -> [String]+formatTriangular fmt m =+   let MatrixShape.Triangular uplo order size = Array.shape m+   in  formatAligned $+       MatrixShape.caseUplo uplo+         padLowerTriangle padUpperTriangle order (Shape.size size) $+       map (printfFloating fmt) $ Array.toList m++padUpperTriangle :: Order -> Int -> [String] -> [[String]]+padUpperTriangle order n xs =+   case order of+      RowMajor ->+         zipWith (++) (iterate ("":) []) (slice (take n $ iterate pred n) xs)+      ColumnMajor ->+         transpose $+         zipWith (++)+            (slice (take n [1..]) xs)+            (reverse $ take n $ iterate ("":) [])++padLowerTriangle :: Order -> Int -> [String] -> [[String]]+padLowerTriangle order n xs =+   case order of+      RowMajor ->+         map (take n) $ map (++ repeat "") $ slice (take n [1..]) xs+      ColumnMajor ->+         transpose $+         zipWith (++) (iterate ("":) []) (slice (take n $ iterate pred n) xs)++_padLowerTriangle :: Order -> Int -> [a] -> [[a]]+_padLowerTriangle order n xs =+   case order of+      RowMajor -> slice (take n [1..]) xs+      ColumnMajor ->+         foldr (\(y:ys) zs -> [y] : zipWith (:) ys zs) []+            (slice (take n $ iterate pred n) xs)++slice :: [Int] -> [a] -> [[a]]+slice ns xs =+   snd $ mapAccumL (\ys n -> let (vs,ws) = splitAt n ys in (ws,vs)) xs ns++formatSeparateTriangle :: [[String]] -> [String]+formatSeparateTriangle xss =+   let strWidths = columnWidths xss    in  zipWith          (\row xs ->             concat $@@ -103,6 +190,11 @@       ColumnMajor -> ListHT.sliceHorizontal (Shape.size height)) .    map (printfFloating fmt) +formatAligned :: [[String]] -> [String]+formatAligned xss =+   let strWidths = columnWidths xss+   in  map (unwords . zipWith (ListHT.padLeft ' ') strWidths) xss+ columnWidths :: [[[a]]] -> [Int] columnWidths xss =    case map (map length) xss of@@ -121,5 +213,8 @@       (Printf printfComplex)       (Printf printfComplex) -printfComplex :: (PrintfArg a) => String -> Complex a -> String-printfComplex fmt (r:+i) = printf (fmt ++ "+i" ++ fmt) r i+printfComplex :: (PrintfArg a, Class.Real a) => String -> Complex a -> String+printfComplex fmt (r:+i) =+   if i<0 || isNegativeZero i+     then printf (fmt ++ "-i" ++ fmt) r (-i)+     else printf (fmt ++ "+i" ++ fmt) r i
+ src/Numeric/LAPACK/Linear/General.hs view
@@ -0,0 +1,88 @@+{-# LANGUAGE TypeFamilies #-}+module Numeric.LAPACK.Linear.General (+   solve,+   inverse,+   ) where++import Numeric.LAPACK.Matrix.Square (Square)+import Numeric.LAPACK.Matrix (General)++import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape+import Numeric.LAPACK.Matrix.Shape.Private (Order(ColumnMajor))+import Numeric.LAPACK.Private (withAutoWorkspace, copyBlock, copyToColumnMajor)++import qualified Numeric.LAPACK.FFI.Generic as LapackGen+import qualified Numeric.Netlib.Utility as Call+import qualified Numeric.Netlib.Class as Class++import qualified Data.Array.Comfort.Storable.Internal as Array+import qualified Data.Array.Comfort.Shape as Shape+import Data.Array.Comfort.Storable.Internal (Array(Array))++import Foreign.Marshal.Alloc (alloca)+import Foreign.C.Types (CInt)+import Foreign.ForeignPtr (withForeignPtr)+import Foreign.Ptr (Ptr)+import Foreign.Storable (peek)++import Control.Monad.Trans.Cont (ContT(ContT), evalContT)+import Control.Monad.IO.Class (liftIO)+import Control.Applicative ((<$>))++import Text.Printf (printf)+++solve ::+   (Shape.C sh, Eq sh, Shape.C nrhs, Class.Floating a) =>+   Square sh a -> General sh nrhs a -> General sh nrhs a+solve+   (Array (MatrixShape.Square orderA shA) a)+   (Array (MatrixShape.General orderB heightB widthB) b) =+      Array.unsafeCreate (MatrixShape.General ColumnMajor heightB widthB) $+         \xPtr -> do+   Call.assert "Square.solve: height shapes mismatch"+      (shA == heightB)+   let n = Shape.size heightB+   let nrhs = Shape.size widthB+   let ldb = n+   evalContT $ do+      nPtr <- Call.cint n+      nrhsPtr <- Call.cint nrhs+      aPtr <- ContT $ withForeignPtr a+      atmpPtr <- Call.allocaArray (n*n)+      ldaPtr <- Call.cint ldb+      ipivPtr <- Call.allocaArray n+      bPtr <- ContT $ withForeignPtr b+      ldbPtr <- Call.cint ldb+      liftIO $ do+         copyToColumnMajor orderA n n aPtr atmpPtr+         copyToColumnMajor orderB n nrhs bPtr xPtr+         withInfo "gesv" $+            LapackGen.gesv nPtr nrhsPtr atmpPtr ldaPtr ipivPtr xPtr ldbPtr+++inverse :: (Shape.C sh, Class.Floating a) => Square sh a -> Square sh a+inverse (Array shape@(MatrixShape.Square _order sh) a) =+      Array.unsafeCreateWithSize shape $ \blockSize bPtr -> do+   let n = Shape.size sh+   evalContT $ do+      nPtr <- Call.cint n+      aPtr <- ContT $ withForeignPtr a+      ldbPtr <- Call.cint n+      ipivPtr <- Call.allocaArray n+      liftIO $ do+         copyBlock blockSize aPtr bPtr+         withInfo "getrf" $ LapackGen.getrf nPtr nPtr bPtr ldbPtr ipivPtr+         withInfo "getri" $ \infoPtr ->+            withAutoWorkspace $ \workPtr lworkPtr ->+               LapackGen.getri nPtr bPtr ldbPtr ipivPtr workPtr lworkPtr infoPtr+++withInfo :: String -> (Ptr CInt -> IO ()) -> IO ()+withInfo name computation = alloca $ \infoPtr -> do+   computation infoPtr+   info <- fromIntegral <$> peek infoPtr+   case compare info (0::Int) of+      EQ -> return ()+      LT -> error $ printf "%s: illegal value in %d-th argument" name (-info)+      GT -> error $ printf "%s: %d-th diagonal value is zero" name info
+ src/Numeric/LAPACK/Linear/Hermitian.hs view
@@ -0,0 +1,87 @@+{-# LANGUAGE TypeFamilies #-}+module Numeric.LAPACK.Linear.Hermitian (+   solve,+   inverse,+   ) where++import Numeric.LAPACK.Matrix.Hermitian (Hermitian)+import Numeric.LAPACK.Matrix (General)++import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape+import Numeric.LAPACK.Matrix.Triangular.Private (copyTriangleToTemp)+import Numeric.LAPACK.Matrix.Shape.Private (Order(ColumnMajor), uploFromOrder)+import Numeric.LAPACK.Private (copyBlock, copyToColumnMajor)++import qualified Numeric.LAPACK.FFI.Generic as LapackGen+import qualified Numeric.Netlib.Utility as Call+import qualified Numeric.Netlib.Class as Class++import qualified Data.Array.Comfort.Storable.Internal as Array+import qualified Data.Array.Comfort.Shape as Shape+import Data.Array.Comfort.Storable.Internal (Array(Array))++import Foreign.Marshal.Alloc (alloca)+import Foreign.C.Types (CInt)+import Foreign.ForeignPtr (withForeignPtr)+import Foreign.Ptr (Ptr)+import Foreign.Storable (peek)++import Control.Monad.Trans.Cont (ContT(ContT), evalContT)+import Control.Monad.IO.Class (liftIO)+import Control.Applicative ((<$>))++import Text.Printf (printf)+++solve ::+   (Shape.C sh, Eq sh, Shape.C nrhs, Class.Floating a) =>+   Hermitian sh a -> General sh nrhs a -> General sh nrhs a+solve+   (Array (MatrixShape.Hermitian orderA shA) a)+   (Array (MatrixShape.General orderB heightB widthB) b) =+      Array.unsafeCreate (MatrixShape.General ColumnMajor heightB widthB) $+         \xPtr -> do+   Call.assert "Hermitian.solve: height shapes mismatch"+      (shA == heightB)+   let n = Shape.size heightB+   let nrhs = Shape.size widthB+   let ldb = n+   evalContT $ do+      uploPtr <- Call.char $ uploFromOrder orderA+      nPtr <- Call.cint n+      nrhsPtr <- Call.cint nrhs+      apPtr <- copyTriangleToTemp orderA n a+      ipivPtr <- Call.allocaArray n+      bPtr <- ContT $ withForeignPtr b+      ldbPtr <- Call.cint ldb+      liftIO $ do+         copyToColumnMajor orderB n nrhs bPtr xPtr+         withInfo "hpsv" $+            LapackGen.hpsv uploPtr nPtr nrhsPtr apPtr ipivPtr xPtr ldbPtr+++inverse ::+   (Shape.C sh, Class.Floating a) => Hermitian sh a -> Hermitian sh a+inverse (Array shape@(MatrixShape.Hermitian order sh) a) =+      Array.unsafeCreateWithSize shape $ \triSize bPtr -> do+   let n = Shape.size sh+   evalContT $ do+      uploPtr <- Call.char $ uploFromOrder order+      nPtr <- Call.cint n+      aPtr <- ContT $ withForeignPtr a+      ipivPtr <- Call.allocaArray n+      workPtr <- Call.allocaArray n+      liftIO $ do+         copyBlock triSize aPtr bPtr+         withInfo "hptrf" $ LapackGen.hptrf uploPtr nPtr bPtr ipivPtr+         withInfo "hptri" $ LapackGen.hptri uploPtr nPtr bPtr ipivPtr workPtr+++withInfo :: String -> (Ptr CInt -> IO ()) -> IO ()+withInfo name computation = alloca $ \infoPtr -> do+   computation infoPtr+   info <- fromIntegral <$> peek infoPtr+   case compare info (0::Int) of+      EQ -> return ()+      LT -> error $ printf "%s: illegal value in %d-th argument" name (-info)+      GT -> error $ printf "%s: %d-th diagonal value is zero" name info
+ src/Numeric/LAPACK/Linear/HermitianPositiveDefinite.hs view
@@ -0,0 +1,105 @@+{-# LANGUAGE TypeFamilies #-}+module Numeric.LAPACK.Linear.HermitianPositiveDefinite (+   solve,+   inverse,+   decompose,+   ) where++import Numeric.LAPACK.Matrix.Hermitian (Hermitian)+import Numeric.LAPACK.Matrix.Triangular (Upper)+import Numeric.LAPACK.Matrix (General)++import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape+import Numeric.LAPACK.Matrix.Triangular.Private (copyTriangleToTemp)+import Numeric.LAPACK.Matrix.Shape.Private (Order(ColumnMajor), uploFromOrder)+import Numeric.LAPACK.Private (copyBlock, copyToColumnMajor)++import qualified Numeric.LAPACK.FFI.Generic as LapackGen+import qualified Numeric.Netlib.Utility as Call+import qualified Numeric.Netlib.Class as Class++import qualified Data.Array.Comfort.Storable.Internal as Array+import qualified Data.Array.Comfort.Shape as Shape+import Data.Array.Comfort.Storable.Internal (Array(Array))++import Foreign.Marshal.Alloc (alloca)+import Foreign.C.Types (CInt)+import Foreign.ForeignPtr (withForeignPtr)+import Foreign.Ptr (Ptr)+import Foreign.Storable (peek)++import Control.Monad.Trans.Cont (ContT(ContT), evalContT)+import Control.Monad.IO.Class (liftIO)+import Control.Applicative ((<$>))++import Text.Printf (printf)+++solve ::+   (Shape.C sh, Eq sh, Shape.C nrhs, Class.Floating a) =>+   Hermitian sh a -> General sh nrhs a -> General sh nrhs a+solve+   (Array (MatrixShape.Hermitian orderA shA) a)+   (Array (MatrixShape.General orderB heightB widthB) b) =+      Array.unsafeCreate (MatrixShape.General ColumnMajor heightB widthB) $+         \xPtr -> do+   Call.assert "Hermitian.solve: height shapes mismatch"+      (shA == heightB)+   let n = Shape.size heightB+   let nrhs = Shape.size widthB+   let ldb = n+   evalContT $ do+      uploPtr <- Call.char $ uploFromOrder orderA+      nPtr <- Call.cint n+      nrhsPtr <- Call.cint nrhs+      apPtr <- copyTriangleToTemp orderA n a+      bPtr <- ContT $ withForeignPtr b+      ldbPtr <- Call.cint ldb+      liftIO $ do+         copyToColumnMajor orderB n nrhs bPtr xPtr+         withInfo "ppsv" $+            LapackGen.ppsv uploPtr nPtr nrhsPtr apPtr xPtr ldbPtr+++inverse ::+   (Shape.C sh, Class.Floating a) => Hermitian sh a -> Hermitian sh a+inverse+   (Array shape@(MatrixShape.Hermitian order sh) a) =+      Array.unsafeCreateWithSize shape $ \triSize bPtr -> do+   evalContT $ do+      uploPtr <- Call.char $ uploFromOrder order+      nPtr <- Call.cint $ Shape.size sh+      aPtr <- ContT $ withForeignPtr a+      liftIO $ do+         copyBlock triSize aPtr bPtr+         withInfo "pptrf" $ LapackGen.pptrf uploPtr nPtr bPtr+         withInfo "pptri" $ LapackGen.pptri uploPtr nPtr bPtr++{- |+Cholesky decomposition+-}+decompose ::+   (Shape.C sh, Class.Floating a) => Hermitian sh a -> Upper sh a+decompose+   (Array (MatrixShape.Hermitian order sh) a) =+      Array.unsafeCreateWithSize+         (MatrixShape.Triangular MatrixShape.Upper order sh) $+            \triSize bPtr -> do+   evalContT $ do+      uploPtr <- Call.char $ uploFromOrder order+      nPtr <- Call.cint $ Shape.size sh+      aPtr <- ContT $ withForeignPtr a+      liftIO $ do+         copyBlock triSize aPtr bPtr+         withInfo "pptrf" $ LapackGen.pptrf uploPtr nPtr bPtr+++withInfo :: String -> (Ptr CInt -> IO ()) -> IO ()+withInfo name computation = alloca $ \infoPtr -> do+   computation infoPtr+   info <- fromIntegral <$> peek infoPtr+   case compare info (0::Int) of+      EQ -> return ()+      LT -> error $ printf "%s: illegal value in %d-th argument" name (-info)+      GT -> error $+         printf "%s: minor of order %d not positive definite" name info
+ src/Numeric/LAPACK/Linear/Triangular.hs view
@@ -0,0 +1,88 @@+{-# LANGUAGE TypeFamilies #-}+module Numeric.LAPACK.Linear.Triangular (+   solve,+   inverse,+   ) where++import Numeric.LAPACK.Matrix.Triangular (Triangular)+import Numeric.LAPACK.Matrix (General)++import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape+import Numeric.LAPACK.Matrix.Shape.Private+         (Order(ColumnMajor),+          transposeFromOrder, uploFromOrder, uploOrder, triangleSize)+import Numeric.LAPACK.Private+         (copyBlock, copyToTemp, copyToColumnMajor)++import qualified Numeric.LAPACK.FFI.Generic as LapackGen+import qualified Numeric.Netlib.Utility as Call+import qualified Numeric.Netlib.Class as Class++import qualified Data.Array.Comfort.Storable.Internal as Array+import qualified Data.Array.Comfort.Shape as Shape+import Data.Array.Comfort.Storable.Internal (Array(Array))++import Foreign.Marshal.Alloc (alloca)+import Foreign.C.Types (CInt)+import Foreign.ForeignPtr (withForeignPtr)+import Foreign.Ptr (Ptr)+import Foreign.Storable (peek)++import Control.Monad.Trans.Cont (ContT(ContT), evalContT)+import Control.Monad.IO.Class (liftIO)+import Control.Applicative ((<$>))++import Text.Printf (printf)+++solve ::+   (MatrixShape.Uplo uplo, Shape.C sh, Eq sh, Shape.C nrhs, Class.Floating a) =>+   Triangular uplo sh a -> General sh nrhs a -> General sh nrhs a+solve+   (Array (MatrixShape.Triangular uplo orderA shA) a)+   (Array (MatrixShape.General orderB heightB widthB) b) =+      Array.unsafeCreate (MatrixShape.General ColumnMajor heightB widthB) $+         \xPtr -> do+   Call.assert "Triangular.solve: height shapes mismatch" (shA == heightB)+   let n = Shape.size heightB+   let nrhs = Shape.size widthB+   let ldb = n+   evalContT $ do+      uploPtr <- Call.char $ uploFromOrder $ uploOrder uplo orderA+      transPtr <- Call.char $ transposeFromOrder orderA+      diagPtr <- Call.char 'N'+      nPtr <- Call.cint n+      nrhsPtr <- Call.cint nrhs+      apPtr <- copyToTemp (triangleSize n) a+      bPtr <- ContT $ withForeignPtr b+      ldbPtr <- Call.cint ldb+      liftIO $ do+         copyToColumnMajor orderB n nrhs bPtr xPtr+         withInfo "tptrs" $+            LapackGen.tptrs uploPtr transPtr diagPtr+               nPtr nrhsPtr apPtr xPtr ldbPtr+++inverse ::+   (MatrixShape.Uplo uplo, Shape.C sh, Class.Floating a) =>+   Triangular uplo sh a -> Triangular uplo sh a+inverse (Array shape@(MatrixShape.Triangular uplo order sh) a) =+      Array.unsafeCreateWithSize shape $ \triSize bPtr -> do+   evalContT $ do+      uploPtr <- Call.char $ uploFromOrder $ uploOrder uplo order+      diagPtr <- Call.char 'N'+      nPtr <- Call.cint $ Shape.size sh+      aPtr <- ContT $ withForeignPtr a+      liftIO $ do+         copyBlock triSize aPtr bPtr+         withInfo "tptri" $ LapackGen.tptri uploPtr diagPtr nPtr bPtr+++withInfo :: String -> (Ptr CInt -> IO ()) -> IO ()+withInfo name computation = alloca $ \infoPtr -> do+   computation infoPtr+   info <- fromIntegral <$> peek infoPtr+   case compare info (0::Int) of+      EQ -> return ()+      LT -> error $ printf "%s: illegal value in %d-th argument" name (-info)+      GT -> error $ printf "%s: %d-th diagonal element zero" name info
− src/Numeric/LAPACK/LinearSystem.hs 
@@ -1,478 +0,0 @@-{-# LANGUAGE TypeFamilies #-}-module Numeric.LAPACK.LinearSystem (-   leastSquares,-   minimumNorm,-   leastSquaresMinimumNorm,-   pseudoInverseRCond,--   Householder,-   householder,-   householderDecompose,-   householderDeterminant,-   determinant,-   householderExtractQ,-   householderExtractR,-   orthogonalComplement,-   ) where--import Numeric.LAPACK.Matrix-         (General, ZeroInt, zeroInt, transpose, identity, dropColumns)--import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape-import Numeric.LAPACK.Matrix.Shape.Private-         (Order(RowMajor, ColumnMajor), charFromOrder)-import Numeric.LAPACK.Vector (Vector)-import Numeric.LAPACK.Private-         (RealOf, zero, one, fill, pointerSeq,-          copyTransposed, copySubMatrix, copyBlock)--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.Internal as Array-import qualified Data.Array.Comfort.Shape as Shape-import Data.Array.Comfort.Storable.Internal (Array(Array))--import System.IO.Unsafe (unsafePerformIO)--import Foreign.Marshal.Array (allocaArray, advancePtr)-import Foreign.Marshal.Alloc (alloca)-import Foreign.C.Types (CInt)-import Foreign.ForeignPtr (withForeignPtr, mallocForeignPtrArray)-import Foreign.Ptr (Ptr)-import Foreign.Storable (Storable, poke, peek)--import Text.Printf (printf)--import Control.Monad.Trans.Cont (ContT(ContT), evalContT)-import Control.Monad.IO.Class (liftIO)-import Control.Monad (when, foldM)-import Control.Applicative ((<$>))--import qualified Data.Complex as Complex-import Data.Complex (Complex)-import Data.Tuple.HT (mapSnd)---{- |-If @x = leastSquares a b@-then @x@ minimizes @Vector.norm2 (multiply a x `sub` b)@.--Precondition: @a@ must have full rank and @height a >= width a@.--}-leastSquares ::-   (Shape.C height, Eq height, Shape.C width, Shape.C nrhs,-    Storable a, Class.Floating a) =>-   General height width a -> General height nrhs a -> General width nrhs a-leastSquares-   (Array shapeA@(MatrixShape.General orderA heightA widthA) a)-   (Array        (MatrixShape.General orderB heightB widthB) b) =-      Array.unsafeCreate (MatrixShape.General ColumnMajor widthA widthB) $-         \xPtr -> do-   Call.assert "leastSquares: height shapes mismatch" (heightA == heightB)-   Call.assert "leastSquares: height of 'a' must be at least the width"-      (Shape.size heightA >= Shape.size widthA)-   let (m,n) = MatrixShape.dimensions shapeA-   let lda = m-   let nrhs = Shape.size widthB-   let ldb = Shape.size heightB-   let ldx = Shape.size widthA-   evalContT $ do-      transPtr <- Call.char $ charFromOrder orderA-      mPtr <- Call.cint m-      nPtr <- Call.cint n-      nrhsPtr <- Call.cint nrhs-      aPtr <- ContT $ withForeignPtr a-      ldaPtr <- Call.cint lda-      let aSize = Shape.size (heightA,widthA)-      atmpPtr <- Call.allocaArray aSize-      liftIO $ copyBlock aSize aPtr atmpPtr-      bPtr <- ContT $ withForeignPtr b-      ldbPtr <- Call.cint ldb-      let bSize = Shape.size (heightB,widthB)-      btmpPtr <- Call.allocaArray bSize-      liftIO $ copyToColumnMajor orderB ldb nrhs bPtr btmpPtr-      liftIO $ withAutoWorkspaceInfo "gels" $-         LapackGen.gels transPtr-            mPtr nPtr nrhsPtr atmpPtr ldaPtr btmpPtr ldbPtr-      liftIO $ copySubMatrix ldx nrhs ldb btmpPtr ldx xPtr--{- |-The vector @x@ with @x = minimumNorm a b@-is the vector with minimal @Vector.norm2 x@-that satisfies @multiply a x == b@.--Precondition: @a@ must have full rank and @height a <= width a@.--}-minimumNorm ::-   (Shape.C height, Eq height, Shape.C width, Shape.C nrhs,-    Storable a, Class.Floating a) =>-   General height width a -> General height nrhs a -> General width nrhs a-minimumNorm-   (Array shapeA@(MatrixShape.General orderA heightA widthA) a)-   (Array        (MatrixShape.General orderB heightB widthB) b) =-      Array.unsafeCreate (MatrixShape.General ColumnMajor widthA widthB) $-         \xPtr -> do-   Call.assert "minimumNorm: height shapes mismatch" (heightA == heightB)-   Call.assert "minimumNorm: width of 'a' must be at least the height"-      (Shape.size widthA >= Shape.size heightA)-   let (m,n) = MatrixShape.dimensions shapeA-   let lda = m-   let nrhs = Shape.size widthB-   let ldb = Shape.size heightB-   let ldx = Shape.size widthA-   evalContT $ do-      transPtr <- Call.char $ charFromOrder orderA-      mPtr <- Call.cint m-      nPtr <- Call.cint n-      nrhsPtr <- Call.cint nrhs-      aPtr <- ContT $ withForeignPtr a-      ldaPtr <- Call.cint lda-      let aSize = Shape.size (heightA,widthA)-      atmpPtr <- Call.allocaArray aSize-      liftIO $ copyBlock aSize aPtr atmpPtr-      bPtr <- ContT $ withForeignPtr b-      ldxPtr <- Call.cint ldx-      liftIO $ copyToSubColumnMajor orderB ldb nrhs bPtr ldx xPtr-      liftIO $ withAutoWorkspaceInfo "gels" $-         LapackGen.gels transPtr-            mPtr nPtr nrhsPtr atmpPtr ldaPtr xPtr ldxPtr--{- |-If @x = leastSquaresMinimumNorm a b@-then @x@ is the vector with minimum @Vector.norm2 x@-that minimizes @Vector.norm2 (multiply a x `sub` b)@.--Matrix @a@ can have any rank-but you must specify the reciprocal condition of the rank-truncated matrix.--}-leastSquaresMinimumNorm ::-   (Shape.C height, Eq height, Shape.C width, Shape.C nrhs,-    Storable a, Class.Floating a) =>-   RealOf a ->-   General height width a -> General height nrhs a ->-   (Int, General width nrhs a)-leastSquaresMinimumNorm rcond-   (Array (MatrixShape.General orderA heightA widthA) a)-   (Array (MatrixShape.General orderB heightB widthB) b) =-      unsafePerformIO $ do-   Call.assert "minimumNorm: height shapes mismatch" (heightA == heightB)-   let shapeX = MatrixShape.General ColumnMajor widthA widthB-   let m = Shape.size heightA-   let n = Shape.size widthA-   let nrhs = Shape.size widthB-   let aSize = m*n-   let lda = m-   let ldtmp = max m n-   let tmpSize = ldtmp*nrhs-   evalContT $ do-      aPtr <- ContT $ withForeignPtr a-      atmpPtr <- Call.allocaArray aSize-      liftIO $ copyToColumnMajor orderA m n aPtr atmpPtr-      ldaPtr <- Call.cint lda-      bPtr <- ContT $ withForeignPtr b-      let needTmp = m>n-      x <- liftIO $ mallocForeignPtrArray $ Shape.size shapeX-      tmpPtr <--         ContT $ if needTmp then allocaArray tmpSize else withForeignPtr x-      ldtmpPtr <- Call.cint ldtmp-      liftIO $ copyToSubColumnMajor orderB m nrhs bPtr ldtmp tmpPtr-      jpvtPtr <- Call.allocaArray n-      rankPtr <- Call.alloca-      gelsy m n nrhs atmpPtr ldaPtr tmpPtr ldtmpPtr jpvtPtr rcond rankPtr-      when needTmp $ liftIO $-         withForeignPtr x $ copySubMatrix n nrhs ldtmp tmpPtr n-      rank <- liftIO $ fromIntegral <$> peek rankPtr-      return (rank, Array shapeX x)---newtype GELSY r a =-   GELSY {-      getGELSY ::-         Int -> Int -> Int -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt ->-         Ptr CInt -> RealOf a -> Ptr CInt -> ContT r IO ()-   }--gelsy ::-   (Class.Floating a) =>-   Int -> Int -> Int ->-   Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt ->-   Ptr CInt -> RealOf a -> Ptr CInt -> ContT r IO ()-gelsy =-   getGELSY $-   Class.switchFloating-      (GELSY gelsyReal)-      (GELSY gelsyReal)-      (GELSY gelsyComplex)-      (GELSY gelsyComplex)--gelsyReal ::-   (Class.Real a, Class.Floating a) =>-   Int -> Int -> Int ->-   Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt ->-   Ptr CInt -> a -> Ptr CInt -> ContT r IO ()-gelsyReal m n nrhs aPtr ldaPtr bPtr ldbPtr jpvtPtr rcond rankPtr = do-   mPtr <- Call.cint m-   nPtr <- Call.cint n-   nrhsPtr <- Call.cint nrhs-   rcondPtr <- Call.real rcond-   liftIO $ withAutoWorkspaceInfo "gelsy" $-      LapackReal.gelsy mPtr nPtr nrhsPtr-         aPtr ldaPtr bPtr ldbPtr jpvtPtr rcondPtr rankPtr--gelsyComplex ::-   (Class.Real a) =>-   Int -> Int -> Int ->-   Ptr (Complex a) -> Ptr CInt -> Ptr (Complex a) -> Ptr CInt ->-   Ptr CInt -> a -> Ptr CInt -> ContT r IO ()-gelsyComplex m n nrhs aPtr ldaPtr bPtr ldbPtr jpvtPtr rcond rankPtr = do-   mPtr <- Call.cint m-   nPtr <- Call.cint n-   nrhsPtr <- Call.cint nrhs-   rcondPtr <- Call.real rcond-   rworkPtr <- Call.allocaArray (2*n)-   liftIO $ withAutoWorkspaceInfo "gelsy" $ \workPtr lworkPtr infoPtr ->-      LapackComplex.gelsy mPtr nPtr nrhsPtr-         aPtr ldaPtr bPtr ldbPtr jpvtPtr rcondPtr rankPtr-         workPtr lworkPtr rworkPtr infoPtr---pseudoInverseRCond ::-   (Shape.C height, Eq height, Shape.C width, Eq width,-    Storable a, Class.Floating a) =>-   RealOf a -> General height width a -> (Int, General width height a)-pseudoInverseRCond rcond a =-   let (MatrixShape.General _ height width) = Array.shape a-   in if Shape.size height < Shape.size width-         then leastSquaresMinimumNorm rcond a $ identity height-         else mapSnd transpose $-              leastSquaresMinimumNorm rcond (transpose a) $-              identity width---type Householder height width = Array (MatrixShape.Householder height width)--{--@(q,r) = householder a@-means that @q@ is unitary and @r@ is upper triangular and @a = multiply q r@.--}-householder ::-   (Shape.C height, Shape.C width, Eq width, Storable a, Class.Floating a) =>-   General height width a ->-   (General height height a, General height width a)-householder a =-   let hh = householderDecompose a-   in  (householderExtractQ hh, householderExtractR $ snd hh)--householderDecompose ::-   (Shape.C height, Shape.C width, Storable a, Class.Floating a) =>-   General height width a -> (Vector width a, Householder height width a)-householderDecompose (Array (MatrixShape.General order height width) a) =-   unsafePerformIO $ do--   let (m,n) =-         case order of-            RowMajor -> (Shape.size width, Shape.size height)-            ColumnMajor -> (Shape.size height, Shape.size width)-   let lda = m-   let mn = min m n-   evalContT $ do-      mPtr <- Call.cint m-      nPtr <- Call.cint n-      aPtr <- ContT $ withForeignPtr a-      ldaPtr <- Call.cint lda-      qr <- liftIO $ mallocForeignPtrArray (m*n)-      qrPtr <- ContT $ withForeignPtr qr-      liftIO $ copyBlock (m*n) aPtr qrPtr-      tau <- liftIO $ mallocForeignPtrArray n-      tauPtr <- ContT $ withForeignPtr tau-      liftIO $ fill zero (n-mn) (advancePtr tauPtr mn)-      liftIO $-         case order of-            RowMajor ->-               withAutoWorkspaceInfo "gelqf" $-                  LapackGen.gelqf mPtr nPtr qrPtr ldaPtr tauPtr-            ColumnMajor ->-               withAutoWorkspaceInfo "geqrf" $-                  LapackGen.geqrf mPtr nPtr qrPtr ldaPtr tauPtr-      return (Array width tau,-              Array (MatrixShape.Householder order height width) qr)--householderDeterminant ::-   (Shape.C height, Shape.C width, Storable a, Class.Floating a) =>-   Householder height width a -> a-householderDeterminant-      (Array (MatrixShape.Householder order height width) a) =-   let m = Shape.size height-       n = Shape.size width-       k = case order of RowMajor -> n; ColumnMajor -> m-   in unsafePerformIO $-      withForeignPtr a $ \aPtr ->-         foldM (\x ptr -> do y <- peek ptr; return $! mul x y) one $-         take (min m n) $ pointerSeq (k+1) aPtr--newtype Mul a = Mul {getMul :: a -> a -> a}--mul :: (Class.Floating a) => a -> a -> a-mul = getMul $ Class.switchFloating (Mul (*)) (Mul (*)) (Mul (*)) (Mul (*))---{-|-Generalized determinant - works also for non-square matrices.-In contrast to the square root of the Gramian determinant-it has the proper sign.--}-determinant ::-   (Shape.C height, Shape.C width, Eq a, Storable a, Class.Floating a) =>-   General height width a -> a-determinant a =-   let (tau,hh) = householderDecompose a-   in  foldl (\x _ -> neg x)-         (householderDeterminant hh)-         (takeWhile (/=zero) $ Array.toList tau)--newtype Neg a = Neg {getNeg :: a -> a}--neg :: (Class.Floating a) => a -> a-neg =-   getNeg $-   Class.switchFloating (Neg negate) (Neg negate) (Neg negate) (Neg negate)---householderExtractQ ::-   (Shape.C height, Shape.C width, Eq width, Storable a, Class.Floating a) =>-   (Vector width a, Householder height width a) -> General height height a-householderExtractQ-   (Array widthTau tau,-    Array (MatrixShape.Householder order height width) qr) =--   Array.unsafeCreate (MatrixShape.General order height height) $ \qPtr -> do--   Call.assert "householderExtractQ: width shapes mismatch" (widthTau == width)--   let m = Shape.size height-   let k = min m $ Shape.size width-   let lda = m-   evalContT $ do-      mPtr <- Call.cint m-      kPtr <- Call.cint k-      qrPtr <- ContT $ withForeignPtr qr-      ldaPtr <- Call.cint lda-      tauPtr <- ContT $ withForeignPtr tau-      liftIO $-         case order of-            RowMajor -> do-               copySubMatrix k m k qrPtr lda qPtr-               withAutoWorkspaceInfo "unglq" $-                  LapackGen.unglq mPtr mPtr kPtr qPtr ldaPtr tauPtr-            ColumnMajor -> do-               copyBlock (m*k) qrPtr qPtr-               withAutoWorkspaceInfo "ungqr" $-                  LapackGen.ungqr mPtr mPtr kPtr qPtr ldaPtr tauPtr--householderExtractR ::-   (Shape.C height, Shape.C width, Eq width, Storable a, Class.Floating a) =>-   Householder height width a -> General height width a-householderExtractR-      (Array (MatrixShape.Householder order height width) qr) =--   Array.unsafeCreate (MatrixShape.General order height width) $-      \rPtr -> do--   let (uplo, (m,n)) =-         case order of-            RowMajor -> ('L', (Shape.size width, Shape.size height))-            ColumnMajor -> ('U', (Shape.size height, Shape.size width))-   fill zero (m*n) rPtr-   evalContT $ do-      uploPtr <- Call.char uplo-      mPtr <- Call.cint m-      nPtr <- Call.cint n-      qrPtr <- ContT $ withForeignPtr qr-      ldqrPtr <- Call.cint m-      ldrPtr <- Call.cint m-      liftIO $ LapackGen.lacpy uploPtr mPtr nPtr qrPtr ldqrPtr rPtr ldrPtr--{- |-For an m-by-n-matrix @a@ with m>=n-this function computes an m-by-(m-n)-matrix @b@-such that @Matrix.multiply (transpose b) a@ is a zero matrix.-The function does not try to compensate a rank deficiency of @a@.-That is, @a|||b@ has full rank if and only if @a@ has full rank.--For full-rank matrices you might also call this @kernel@ or @nullspace@.--}-orthogonalComplement ::-   (Shape.C height, Shape.C width, Eq width, Storable a, Class.Floating a) =>-   General height width a -> General height ZeroInt a-orthogonalComplement a =-   dropColumns (Shape.size $ MatrixShape.generalWidth $ Array.shape a) $-   Array.mapShape zeroIntWidth $ householderExtractQ $ householderDecompose a--zeroIntWidth ::-   (Shape.C width) =>-   MatrixShape.General height width -> MatrixShape.General height ZeroInt-zeroIntWidth (MatrixShape.General order height width) =-   MatrixShape.General order height (zeroInt $ Shape.size width)----withAutoWorkspaceInfo ::-   (Storable a, Class.Floating a) =>-   String -> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()-withAutoWorkspaceInfo name computation = evalContT $ do-   infoPtr <- Call.alloca-   liftIO $ withAutoWorkspace $ \workPtr lworkPtr ->-      computation workPtr lworkPtr infoPtr-   info <- liftIO $ fromIntegral <$> peek infoPtr-   case compare info (0::Int) of-      EQ -> return ()-      LT -> error $ printf "%s: illegal value in %d-th argument" name (-info)-      GT -> error $ printf "%s: deficient rank %d" name info--withAutoWorkspace ::-   (Storable a, Class.Floating a) =>-   (Ptr a -> Ptr CInt -> IO ()) -> IO ()-withAutoWorkspace computation = evalContT $ do-   lworkPtr <- Call.cint (-1)-   lwork <- liftIO $ alloca $ \workPtr -> do-      computation workPtr lworkPtr-      ceilingSize <$> peek workPtr-   workPtr <- Call.allocaArray lwork-   liftIO $ poke lworkPtr $ fromIntegral lwork-   liftIO $ computation workPtr lworkPtr---copyToColumnMajor ::-   (Storable a, Class.Floating a) =>-   Order -> Int -> Int -> Ptr a -> Ptr a -> IO ()-copyToColumnMajor order m n aPtr bPtr =-   case order of-      RowMajor -> copyTransposed m n aPtr m bPtr-      ColumnMajor -> copyBlock (m*n) aPtr bPtr--copyToSubColumnMajor ::-   (Storable a, Class.Floating a) =>-   Order -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()-copyToSubColumnMajor order m n aPtr ldb bPtr =-   case order of-      RowMajor -> copyTransposed m n aPtr ldb bPtr-      ColumnMajor ->-         if m==ldb-           then copyBlock (m*n) aPtr bPtr-           else copySubMatrix m n m aPtr ldb bPtr---newtype FuncArg b a = FuncArg {runFuncArg :: a -> b}--ceilingSize :: (Class.Floating a) => a -> Int-ceilingSize =-   runFuncArg $-   Class.switchFloating-      (FuncArg ceiling)-      (FuncArg ceiling)-      (FuncArg $ ceiling . Complex.realPart)-      (FuncArg $ ceiling . Complex.realPart)
src/Numeric/LAPACK/Matrix.hs view
@@ -2,11 +2,12 @@ {-# LANGUAGE TypeOperators #-} module Numeric.LAPACK.Matrix (    General,-   (Format.##),-   Format.Format,-   Format.FormatArray,+   (##),+   Format,+   FormatArray,    ZeroInt, zeroInt,-   transpose,+   transpose, adjoint,+   fromScalar, toScalar,    fromList,    identity,    diagonal, getDiagonal,@@ -27,15 +28,21 @@    multiply,    multiplyVector, -   trace,+   Multiply, (<#>),+   MultiplyLeft, (<#),+   MultiplyRight, (#>),    ) where  import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape-import qualified Numeric.LAPACK.Private as Private-import qualified Numeric.LAPACK.Format as Format+import qualified Numeric.LAPACK.Matrix.Square as Square import qualified Numeric.LAPACK.Vector as Vector-import Numeric.LAPACK.Matrix.Shape.Private-         (Order(RowMajor, ColumnMajor), charFromOrder)+import Numeric.LAPACK.Format (Format, FormatArray, (##))+import Numeric.LAPACK.Matrix.Shape.Private (Order(RowMajor, ColumnMajor))+import Numeric.LAPACK.Matrix.Multiply+         (Multiply((<#>)), MultiplyLeft((<#)), MultiplyRight((#>)),+          transpose, multiplyVector, multiply, multiplyVectorUnchecked)+import Numeric.LAPACK.Matrix.Private (General, ZeroInt, zeroInt)+import Numeric.LAPACK.Vector (Vector) import Numeric.LAPACK.Private          (zero, one, pointerSeq, copyTransposed, copySubMatrix, copyBlock) @@ -64,18 +71,21 @@ import Data.Bool.HT (if')  -type General height width = Array (MatrixShape.General height width)---transpose :: General height width a -> General width height a-transpose = Array.mapShape MatrixShape.transpose+{- |+conjugate transpose+-}+adjoint ::+   (Shape.C height, Shape.C width, Class.Floating a) =>+   General height width a -> General width height a+adjoint = transpose . Vector.conjugate  -type ZeroInt = Shape.ZeroBased Int--zeroInt :: Int -> ZeroInt-zeroInt = Shape.ZeroBased+fromScalar :: (Storable a) => a -> General () () a+fromScalar = Square.toGeneral . Square.fromScalar +toScalar :: (Storable a) => General () () a -> a+toScalar (Array (MatrixShape.General _ () ()) a) =+   unsafePerformIO $ withForeignPtr a peek  fromList ::    (Shape.C height, Shape.C width, Storable a) =>@@ -84,70 +94,20 @@    Array.fromList (MatrixShape.General RowMajor height width)  -type Vector = Array---identity, _identity ::-   (Shape.C sh, Storable a, Class.Floating a) =>+identity ::+   (Shape.C sh, Class.Floating a) =>    sh -> General sh sh a-identity sh =-   Array.unsafeCreate (MatrixShape.General ColumnMajor sh sh) $ \aPtr ->-   evalContT $ do-      uploPtr <- Call.char 'A'-      nPtr <- Call.cint $ Shape.size sh-      alphaPtr <- Call.number zero-      betaPtr <- Call.number one-      liftIO $ LapackGen.laset uploPtr nPtr nPtr alphaPtr betaPtr aPtr nPtr--_identity sh =-   Array.unsafeCreate (MatrixShape.General ColumnMajor sh sh) $ \yPtr ->-   evalContT $ do-      nPtr <- Call.alloca-      xPtr <- Call.number zero-      incxPtr <- Call.cint 0-      incyPtr <- Call.cint 1-      liftIO $ do-         let n = fromIntegral $ Shape.size sh-         poke nPtr $ n*n-         BlasGen.copy nPtr xPtr incxPtr yPtr incyPtr-         poke nPtr n-         poke xPtr one-         poke incyPtr (n+1)-         BlasGen.copy nPtr xPtr incxPtr yPtr incyPtr+identity = Square.toGeneral . Square.identity  diagonal ::-   (Shape.C sh, Storable a, Class.Floating a) =>+   (Shape.C sh, Class.Floating a) =>    Vector sh a -> General sh sh a-diagonal (Array sh x) =-   Array.unsafeCreate (MatrixShape.General ColumnMajor sh sh) $ \yPtr ->-   evalContT $ do-      nPtr <- Call.alloca-      xPtr <- ContT $ withForeignPtr x-      zPtr <- Call.number zero-      incxPtr <- Call.cint 1-      incyPtr <- Call.cint 1-      inczPtr <- Call.cint 0-      liftIO $ do-         let n = fromIntegral $ Shape.size sh-         poke nPtr $ n*n-         BlasGen.copy nPtr zPtr inczPtr yPtr incyPtr-         poke nPtr n-         poke incyPtr (n+1)-         BlasGen.copy nPtr xPtr incxPtr yPtr incyPtr+diagonal = Square.toGeneral . Square.diagonal  getDiagonal ::-   (Shape.C sh, Eq sh, Storable a, Class.Floating a) =>+   (Shape.C sh, Eq sh, Class.Floating a) =>    General sh sh a -> Vector sh a-getDiagonal (Array (MatrixShape.General _ height width) x) =-      Array.unsafeCreate height $ \yPtr -> do-   Call.assert "getDiagonal: non-square matrix" (height==width)-   evalContT $ do-      let n = Shape.size height-      nPtr <- Call.cint n-      xPtr <- ContT $ withForeignPtr x-      incxPtr <- Call.cint (n+1)-      incyPtr <- Call.cint 1-      liftIO $ BlasGen.copy nPtr xPtr incxPtr yPtr incyPtr+getDiagonal = Square.getDiagonal . Square.fromGeneral   singleRow :: Vector width a -> General () width a@@ -209,7 +169,7 @@  pickRow ::    (Shape.C height, Shape.C width, Shape.Index height ~ ix,-    Storable a, Class.Floating a) =>+    Class.Floating a) =>    General height width a -> ix -> Vector width a pickRow (Array (MatrixShape.General order height width) x) ix =    case order of@@ -218,7 +178,7 @@  pickColumn ::    (Shape.C height, Shape.C width, Shape.Index width ~ ix,-    Storable a, Class.Floating a) =>+    Class.Floating a) =>    General height width a -> ix -> Vector height a pickColumn (Array (MatrixShape.General order height width) x) ix =    case order of@@ -227,11 +187,10 @@  pickConsecutive ::    (Shape.C height, Shape.C width, Shape.Index height ~ ix,-    Storable a, Class.Floating a) =>+    Class.Floating a) =>    height -> width -> ForeignPtr a -> ix -> Vector width a pickConsecutive height width x ix =-   Array.unsafeCreate width $ \yPtr -> evalContT $ do-      let n = Shape.size width+   Array.unsafeCreateWithSize width $ \n yPtr -> evalContT $ do       let offset = Shape.offset height ix       nPtr <- Call.cint n       xPtr <- ContT $ withForeignPtr x@@ -242,11 +201,10 @@  pickScattered ::    (Shape.C height, Shape.C width, Shape.Index width ~ ix,-    Storable a, Class.Floating a) =>+    Class.Floating a) =>    height -> width -> ForeignPtr a -> ix -> Vector height a pickScattered height width x ix =-   Array.unsafeCreate height $ \yPtr -> evalContT $ do-      let n = Shape.size height+   Array.unsafeCreateWithSize height $ \n yPtr -> evalContT $ do       let offset = Shape.offset width ix       nPtr <- Call.cint n       xPtr <- ContT $ withForeignPtr x@@ -257,18 +215,19 @@   takeRows, dropRows ::-   (Shape.C width, Storable a, Class.Floating a) =>+   (Shape.C width, Class.Floating a) =>    Int -> General ZeroInt width a -> General ZeroInt width a takeRows k       (Array (MatrixShape.General order (Shape.ZeroBased heightA) width) a) =    let heightB = min k heightA        n = Shape.size width    in if' (k<0) (error "take: negative number") $-      Array.unsafeCreate-         (MatrixShape.General order (Shape.ZeroBased heightB) width) $ \bPtr ->+      Array.unsafeCreateWithSize+         (MatrixShape.General order (Shape.ZeroBased heightB) width) $+            \blockSize bPtr ->       withForeignPtr a $ \aPtr ->       case order of-         RowMajor -> copyBlock (heightB*n) aPtr bPtr+         RowMajor -> copyBlock blockSize aPtr bPtr          ColumnMajor -> copySubMatrix heightB n heightA aPtr heightB bPtr  dropRows k0@@ -277,17 +236,18 @@        heightB = heightA - k        n = Shape.size width    in if' (k<0) (error "take: negative number") $-      Array.unsafeCreate-         (MatrixShape.General order (Shape.ZeroBased heightB) width) $ \bPtr ->+      Array.unsafeCreateWithSize+         (MatrixShape.General order (Shape.ZeroBased heightB) width) $+            \blockSize bPtr ->       withForeignPtr a $ \aPtr ->       case order of-         RowMajor -> copyBlock (heightB*n) (advancePtr aPtr (k*n)) bPtr+         RowMajor -> copyBlock blockSize (advancePtr aPtr (k*n)) bPtr          ColumnMajor ->             copySubMatrix heightB n heightA (advancePtr aPtr k) heightB bPtr   takeColumns, dropColumns ::-   (Shape.C height, Storable a, Class.Floating a) =>+   (Shape.C height, Class.Floating a) =>    Int -> General height ZeroInt a -> General height ZeroInt a takeColumns k = transpose . takeRows k . transpose dropColumns k = transpose . dropRows k . transpose@@ -295,10 +255,10 @@  -- alternative: laswp reverseRows ::-   (Shape.C width, Storable a, Class.Floating a) =>+   (Shape.C width, Class.Floating a) =>    General ZeroInt width a -> General ZeroInt width a reverseRows (Array shape@(MatrixShape.General order height width) a) =-   Array.unsafeCreate shape $ \bPtr -> evalContT $ do+   Array.unsafeCreateWithSize shape $ \blockSize bPtr -> evalContT $ do       let n = Shape.size height       let m = Shape.size width       fwdPtr <- Call.bool True@@ -307,26 +267,26 @@       kPtr <- Call.allocaArray n       aPtr <- ContT $ withForeignPtr a       liftIO $ do-         copyBlock (n*m) aPtr bPtr+         copyBlock blockSize aPtr bPtr          pokeArray kPtr $ take n $ iterate (subtract 1) $ fromIntegral n          case order of             RowMajor -> LapackGen.lapmt fwdPtr mPtr nPtr bPtr mPtr kPtr             ColumnMajor -> LapackGen.lapmr fwdPtr nPtr mPtr bPtr nPtr kPtr  reverseColumns ::-   (Shape.C height, Storable a, Class.Floating a) =>+   (Shape.C height, Class.Floating a) =>    General height ZeroInt a -> General height ZeroInt a reverseColumns = transpose . reverseRows . transpose   fromRowMajor ::-   (Shape.C height, Shape.C width, Storable a, Class.Floating a) =>+   (Shape.C height, Shape.C width, Class.Floating a) =>    Array (height,width) a -> General height width a fromRowMajor (Array (height,width) x) =    Array (MatrixShape.General RowMajor height width) x  toRowMajor ::-   (Shape.C height, Shape.C width, Storable a, Class.Floating a) =>+   (Shape.C height, Shape.C width, Class.Floating a) =>    General height width a -> Array (height,width) a toRowMajor (Array (MatrixShape.General order height width) x) =    let shape = (height, width)@@ -347,12 +307,11 @@                   (pointerSeq n yPtr)  flatten ::-   (Shape.C height, Shape.C width, Storable a, Class.Floating a) =>+   (Shape.C height, Shape.C width, Class.Floating a) =>    General height width a -> Vector ZeroInt a flatten x =    case toRowMajor x of-      Array (height,width) fptr ->-         Array (zeroInt $ Shape.size height * Shape.size width) fptr+      Array shape fptr -> Array (zeroInt $ Shape.size shape) fptr   infixl 3 |||@@ -360,7 +319,7 @@  (|||) ::    (Shape.C height, Eq height, Shape.C widtha, Shape.C widthb,-    Storable a, Class.Floating a) =>+    Class.Floating a) =>    General height widtha a ->    General height widthb a ->    General height (widtha:+:widthb) a@@ -440,7 +399,7 @@  (===) ::    (Shape.C width, Eq width, Shape.C heighta, Shape.C heightb,-    Storable a, Class.Floating a) =>+    Class.Floating a) =>    General heighta width a ->    General heightb width a ->    General (heighta:+:heightb) width a@@ -448,93 +407,22 @@   rowSums ::-   (Shape.C height, Shape.C width, Storable a, Class.Floating a) =>+   (Shape.C height, Shape.C width, Class.Floating a) =>    General height width a -> Vector height a rowSums m =    let MatrixShape.General _ _ width = Array.shape m    in  multiplyVectorUnchecked m (Vector.constant width one)  columnSums ::-   (Shape.C height, Shape.C width, Storable a, Class.Floating a) =>+   (Shape.C height, Shape.C width, Class.Floating a) =>    General height width a -> Vector width a columnSums m =    let MatrixShape.General _ height _ = Array.shape m    in  multiplyVectorUnchecked (transpose m) (Vector.constant height one) -multiplyVector ::-   (Shape.C height, Shape.C width, Eq width,-    Storable a, Class.Floating a) =>-   General height width a -> Vector width a -> Vector height a-multiplyVector a x =-   let MatrixShape.General _order _height width = Array.shape a-   in if width == Array.shape x-         then multiplyVectorUnchecked a x-         else error "multiplyVector: width shapes mismatch" -multiplyVectorUnchecked ::-   (Shape.C height, Shape.C width,-    Storable a, Class.Floating a) =>-   General height width a -> Vector width a -> Vector height a-multiplyVectorUnchecked-   (Array shape@(MatrixShape.General order height _width) a) (Array _ x) =-      Array.unsafeCreate height $ \yPtr -> do-   let (m,n) = MatrixShape.dimensions shape-   let lda = m-   evalContT $ do-      transPtr <- Call.char $ charFromOrder order-      mPtr <- Call.cint m-      nPtr <- Call.cint n-      alphaPtr <- Call.number one-      aPtr <- ContT $ withForeignPtr a-      ldaPtr <- Call.cint lda-      xPtr <- ContT $ withForeignPtr x-      incxPtr <- Call.cint 1-      betaPtr <- Call.number zero-      incyPtr <- Call.cint 1-      liftIO $-         BlasGen.gemv-            transPtr mPtr nPtr alphaPtr aPtr ldaPtr-            xPtr incxPtr betaPtr yPtr incyPtr--multiply ::-   (Shape.C height,-    Shape.C fuse, Eq fuse,-    Shape.C width,-    Storable a, Class.Floating a) =>-   General height fuse a -> General fuse width a -> General height width a-multiply-   (Array (MatrixShape.General orderA height fuseA) a)-   (Array (MatrixShape.General orderB fuseB width) b) =-      Array.unsafeCreate (MatrixShape.General ColumnMajor height width) $-         \cPtr -> do-   Call.assert "multiply: fuse shapes mismatch" (fuseA == fuseB)-   let m = Shape.size height-   let n = Shape.size width-   let k = Shape.size fuseA-   let lda = case orderA of RowMajor -> k; ColumnMajor -> m-   let ldb = case orderB of RowMajor -> n; ColumnMajor -> k-   let ldc = m-   evalContT $ do-      transaPtr <- Call.char $ charFromOrder orderA-      transbPtr <- Call.char $ charFromOrder orderB-      mPtr <- Call.cint m-      nPtr <- Call.cint n-      kPtr <- Call.cint k-      alphaPtr <- Call.number one-      aPtr <- ContT $ withForeignPtr a-      ldaPtr <- Call.cint lda-      bPtr <- ContT $ withForeignPtr b-      ldbPtr <- Call.cint ldb-      betaPtr <- Call.number zero-      ldcPtr <- Call.cint ldc-      liftIO $-         BlasGen.gemm-            transaPtr transbPtr mPtr nPtr kPtr alphaPtr aPtr ldaPtr-            bPtr ldbPtr betaPtr cPtr ldcPtr-- scaleRows ::-   (Shape.C height, Eq height, Shape.C width, Storable a, Class.Floating a) =>+   (Shape.C height, Eq height, Shape.C width, Class.Floating a) =>    Vector height a -> General height width a -> General height width a scaleRows    (Array heightX x) (Array shape@(MatrixShape.General order height width) a) =@@ -583,14 +471,6 @@                (pointerSeq n bPtr)  scaleColumns ::-   (Shape.C height, Shape.C width, Eq width, Storable a, Class.Floating a) =>+   (Shape.C height, Shape.C width, Eq width, Class.Floating a) =>    Vector width a -> General height width a -> General height width a scaleColumns x = transpose . scaleRows x . transpose----trace :: (Shape.C sh, Eq sh, Class.Floating a) => General sh sh a -> a-trace (Array (MatrixShape.General _ height width) x) = unsafePerformIO $ do-   Call.assert "trace: non-square matrix" (height==width)-   let n = Shape.size height-   withForeignPtr x $ \xPtr -> Private.sum n xPtr (n+1)
+ src/Numeric/LAPACK/Matrix/Hermitian.hs view
@@ -0,0 +1,540 @@+{-# LANGUAGE TypeFamilies #-}+module Numeric.LAPACK.Matrix.Hermitian (+   Hermitian,+   fromList,+   autoFromList,+   identity,+   diagonal,+   getDiagonal,++   multiplyVector,+   square,+   multiplySquareLeft,+   multiplyGeneralLeft,+   multiplySquareRight,+   multiplyGeneralRight,+   outer,+   sumRank1,+   sumRank2,++   toSquare,+   covariance,+   addTransposed,+   ) where++import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape+import Numeric.LAPACK.Matrix.Triangular.Private+         (forPointers, pack, unpack, unpackToTemp,+          diagonalPointers, rowMajorPointers, columnMajorPointers)+import Numeric.LAPACK.Matrix.Shape.Private+         (Order(RowMajor,ColumnMajor), flipOrder, uploFromOrder)+import Numeric.LAPACK.Matrix.Square (Square)+import Numeric.LAPACK.Matrix.Private (General, ZeroInt, zeroInt)+import Numeric.LAPACK.Vector (Vector)+import Numeric.LAPACK.Private+         (RealOf, fill, zero, one, lacgv, fromReal, realPart, copyToTemp)++import qualified Numeric.LAPACK.FFI.Complex as LapackComplex+import qualified Numeric.BLAS.FFI.Generic as BlasGen+import qualified Numeric.BLAS.FFI.Complex as BlasComplex+import qualified Numeric.BLAS.FFI.Real as BlasReal+import qualified Numeric.Netlib.Utility as Call+import qualified Numeric.Netlib.Class as Class++import qualified Data.Array.Comfort.Storable.Internal as Array+import qualified Data.Array.Comfort.Shape as Shape+import Data.Array.Comfort.Storable.Internal (Array(Array))++import Foreign.C.Types (CInt, CChar)+import Foreign.ForeignPtr (ForeignPtr, withForeignPtr)+import Foreign.Ptr (Ptr)+import Foreign.Storable (Storable, poke, peek)++import Control.Monad.Trans.Cont (ContT(ContT), evalContT)+import Control.Monad.IO.Class (liftIO)+import Control.Monad (when)++import qualified Data.NonEmpty as NonEmpty+import Data.Foldable (forM_)+import Data.Complex (Complex)+++type Hermitian sh = Array (MatrixShape.Hermitian sh)+++fromList :: (Shape.C sh, Storable a) => Order -> sh -> [a] -> Hermitian sh a+fromList order sh =+   Array.fromList (MatrixShape.Hermitian order sh)++autoFromList :: (Storable a) => Order -> [a] -> Hermitian ZeroInt a+autoFromList order xs =+   fromList order+      (zeroInt $ MatrixShape.triangleExtent "Hermitian.autoFromList" $+       length xs)+      xs+++identity :: (Shape.C sh, Class.Floating a) => Order -> sh -> Hermitian sh a+identity order sh =+   Array.unsafeCreateWithSize (MatrixShape.Hermitian order sh) $+      \triSize aPtr -> do+   fill zero triSize aPtr+   mapM_ (flip poke one . snd) $+      diagonalPointers order (Shape.size sh) aPtr aPtr++diagonal ::+   (Shape.C sh, Class.Floating a) =>+   Order -> Vector sh (RealOf a) -> Hermitian sh a+diagonal order =+   runDiagonal $+   Class.switchFloating+      (Diagonal $ diagonalAux order) (Diagonal $ diagonalAux order)+      (Diagonal $ diagonalAux order) (Diagonal $ diagonalAux order)++newtype Diagonal sh a =+   Diagonal {runDiagonal :: Vector sh (RealOf a) -> Hermitian sh a}++diagonalAux ::+   (Shape.C sh, Class.Floating a, RealOf a ~ ar, Storable ar) =>+   Order -> Vector sh ar -> Hermitian sh a+diagonalAux order (Array sh x) =+   Array.unsafeCreateWithSize (MatrixShape.Hermitian order sh) $+      \triSize aPtr -> do+   fill zero triSize aPtr+   withForeignPtr x $ \xPtr ->+      forM_ (diagonalPointers order (Shape.size sh) xPtr aPtr) $+         \(srcPtr,dstPtr) -> poke dstPtr . fromReal =<< peek srcPtr+++getDiagonal ::+   (Shape.C sh, Class.Floating a) =>+   Hermitian sh a -> Vector sh (RealOf a)+getDiagonal =+   runGetDiagonal $+   Class.switchFloating+      (GetDiagonal $ getDiagonalAux) (GetDiagonal $ getDiagonalAux)+      (GetDiagonal $ getDiagonalAux) (GetDiagonal $ getDiagonalAux)++newtype GetDiagonal sh a =+   GetDiagonal {runGetDiagonal :: Hermitian sh a -> Vector sh (RealOf a)}++getDiagonalAux ::+   (Shape.C sh, Class.Floating a, RealOf a ~ ar, Storable ar) =>+   Hermitian sh a -> Vector sh ar+getDiagonalAux (Array (MatrixShape.Hermitian order sh) a) =+   Array.unsafeCreateWithSize sh $ \n xPtr ->+   withForeignPtr a $ \aPtr ->+      forM_ (diagonalPointers order n xPtr aPtr) $+         \(dstPtr,srcPtr) -> poke dstPtr . realPart =<< peek srcPtr+++multiplyVector ::+   (Shape.C sh, Eq sh, Class.Floating a) =>+   Hermitian sh a -> Vector sh a -> Vector sh a+multiplyVector (Array (MatrixShape.Hermitian order shA) a) (Array shX x) =+      Array.unsafeCreateWithSize shX $ \n yPtr -> do+   Call.assert "Hermitian.multiplyVector: width shapes mismatch" (shA == shX)+   evalContT $ do+      uploPtr <- Call.char $ uploFromOrder order+      nPtr <- Call.cint n+      alphaPtr <- Call.number one+      aPtr <- ContT $ withForeignPtr a+      xPtr <- ContT $ withForeignPtr x+      incxPtr <- Call.cint 1+      betaPtr <- Call.number zero+      incyPtr <- Call.cint 1+      liftIO $+         BlasGen.hpmv+            uploPtr nPtr alphaPtr aPtr xPtr incxPtr betaPtr yPtr incyPtr+++square ::+   (Shape.C sh, Eq sh, Class.Floating a) =>+   Hermitian sh a -> Hermitian sh a+square+   (Array shape@(MatrixShape.Hermitian order sh) a) =+      Array.unsafeCreate shape $ \cpPtr -> do+   let n = Shape.size sh+   evalContT $ do+      sidePtr <- Call.char 'L'+      uploPtr <- Call.char 'U'+      nPtr <- Call.cint n+      let ldPtr = nPtr+      bPtr <- unpackToTemp (unpackFull order) n a+      cPtr <- Call.allocaArray (n*n)+      alphaPtr <- Call.number one+      betaPtr <- Call.number zero+      liftIO $ do+         BlasGen.hemm sidePtr uploPtr+            nPtr nPtr alphaPtr bPtr ldPtr+            bPtr ldPtr betaPtr cPtr ldPtr+         pack order n cPtr cpPtr+++multiplySquareLeft ::+   (Shape.C sh, Eq sh, Class.Floating a) =>+   Square sh a -> Hermitian sh a -> Square sh a+multiplySquareLeft+   (Array shapeB@(MatrixShape.Square orderB shB) b)+   (Array        (MatrixShape.Hermitian orderA shA) a) =+      Array.unsafeCreate shapeB $ \cPtr -> do+   Call.assert "Hermitian.multiplySquareLeft: shapes mismatch" (shA == shB)+   let n = Shape.size shB+   multiplyAux True orderA n a (flipOrder orderB) n b cPtr++multiplyGeneralLeft ::+   (Shape.C height, Shape.C width, Eq width, Class.Floating a) =>+   General height width a -> Hermitian width a -> General height width a+multiplyGeneralLeft+   (Array shapeB@(MatrixShape.General orderB height width) b)+   (Array        (MatrixShape.Hermitian orderA shA) a) =+      Array.unsafeCreate shapeB $ \cPtr -> do+   Call.assert "Hermitian.multiplyGeneralLeft: shapes mismatch" (shA == width)+   multiplyAux True+      orderA (Shape.size width) a (flipOrder orderB) (Shape.size height) b cPtr++multiplySquareRight ::+   (Shape.C sh, Eq sh, Class.Floating a) =>+   Hermitian sh a -> Square sh a -> Square sh a+multiplySquareRight+   (Array        (MatrixShape.Hermitian orderA shA) a)+   (Array shapeB@(MatrixShape.Square orderB shB) b) =+      Array.unsafeCreate shapeB $ \cPtr -> do+   Call.assert "Hermitian.multiplySquareRight: shapes mismatch" (shA == shB)+   let n = Shape.size shB+   multiplyAux False orderA n a orderB n b cPtr++multiplyGeneralRight ::+   (Shape.C height, Eq height, Shape.C width, Class.Floating a) =>+   Hermitian height a -> General height width a -> General height width a+multiplyGeneralRight+   (Array        (MatrixShape.Hermitian orderA shA) a)+   (Array shapeB@(MatrixShape.General orderB height width) b) =+      Array.unsafeCreate shapeB $ \cPtr -> do+   Call.assert "Hermitian.multiplyGeneralRight: shapes mismatch" (shA == height)+   multiplyAux False+      orderA (Shape.size height) a orderB (Shape.size width) b cPtr++multiplyAux ::+   Class.Floating a =>+   Bool ->+   Order -> Int -> ForeignPtr a ->+   Order -> Int -> ForeignPtr a -> Ptr a -> IO ()+multiplyAux extraConjugate orderA m0 a orderB n0 b cPtr = do+   let size = m0*m0+   evalContT $ do+      let (side,(m,n)) =+            case orderB of+               ColumnMajor -> ('L',(m0,n0))+               RowMajor -> ('R',(n0,m0))+      sidePtr <- Call.char side+      uploPtr <- Call.char $ uploFromOrder orderA+      mPtr <- Call.cint m+      nPtr <- Call.cint n+      alphaPtr <- Call.number one+      aPtr <- unpackToTemp (unpack orderA) m0 a+      ldaPtr <- Call.cint m0+      incaPtr <- Call.cint 1+      sizePtr <- Call.cint size+      bPtr <- ContT $ withForeignPtr b+      ldbPtr <- Call.cint m+      betaPtr <- Call.number zero+      ldcPtr <- Call.cint m+      liftIO $ do+         when ((orderA/=orderB) /= extraConjugate) $+            lacgv sizePtr aPtr incaPtr+         BlasGen.hemm sidePtr uploPtr+            mPtr nPtr alphaPtr aPtr ldaPtr+            bPtr ldbPtr betaPtr cPtr ldcPtr+++outer :: (Shape.C sh, Class.Floating a) => Vector sh a -> Hermitian sh a+outer =+   getMap $+   Class.switchFloating+      (Map outerAux) (Map outerAux)+      (Map outerAux) (Map outerAux)++outerAux ::+   (Shape.C sh, Class.Floating a, RealOf a ~ ar, Class.Real ar) =>+   Vector sh a -> Hermitian sh a+outerAux (Array sh x) =+   Array.unsafeCreateWithSize (MatrixShape.Hermitian ColumnMajor sh) $+      \triSize aPtr -> do+   let n = Shape.size sh+   evalContT $ do+      uploPtr <- Call.char $ uploFromOrder ColumnMajor+      nPtr <- Call.cint n+      alphaPtr <- Call.real one+      xPtr <- ContT $ withForeignPtr x+      incxPtr <- Call.cint 1+      liftIO $ fill zero triSize aPtr+      liftIO $ hpr uploPtr nPtr alphaPtr xPtr incxPtr aPtr+++sumRank1 ::+   (Shape.C sh, Eq sh, Class.Floating a) =>+   NonEmpty.T [] (RealOf a, Vector sh a) -> Hermitian sh a+sumRank1 =+   getSumRank1 $+   Class.switchFloating+      (SumRank1 sumRank1Aux) (SumRank1 sumRank1Aux)+      (SumRank1 sumRank1Aux) (SumRank1 sumRank1Aux)++type SumRank1_ sh a = NonEmpty.T [] (RealOf a, Vector sh a) -> Hermitian sh a++newtype SumRank1 sh a = SumRank1 {getSumRank1 :: SumRank1_ sh a}++sumRank1Aux ::+   (Shape.C sh, Eq sh, Class.Floating a, RealOf a ~ ar, Class.Real ar) =>+   SumRank1_ sh a+sumRank1Aux xs@(NonEmpty.Cons (_, Array sh _) _) =+   Array.unsafeCreateWithSize (MatrixShape.Hermitian ColumnMajor sh) $+      \triSize aPtr -> do+   let n = Shape.size sh+   evalContT $ do+      uploPtr <- Call.char $ uploFromOrder ColumnMajor+      nPtr <- Call.cint n+      alphaPtr <- Call.alloca+      incxPtr <- Call.cint 1+      liftIO $ do+         fill zero triSize aPtr+         forM_ xs $ \(alpha, Array shX x) ->+            withForeignPtr x $ \xPtr -> do+               Call.assert+                  "Hermitian.sumRank1: non-matching vector size" (sh==shX)+               poke alphaPtr alpha+               hpr uploPtr nPtr alphaPtr xPtr incxPtr aPtr+++type HPR_ a =+   Ptr CChar -> Ptr CInt ->+   Ptr (RealOf a) -> Ptr a -> Ptr CInt -> Ptr a -> IO ()++newtype HPR a = HPR {getHPR :: HPR_ a}++hpr :: Class.Floating a => HPR_ a+hpr =+   getHPR $+   Class.switchFloating+      (HPR BlasReal.spr) (HPR BlasReal.spr)+      (HPR BlasComplex.hpr) (HPR BlasComplex.hpr)+++sumRank2 ::+   (Shape.C sh, Eq sh, Class.Floating a) =>+   NonEmpty.T [] (a, (Vector sh a, Vector sh a)) -> Hermitian sh a+sumRank2 xys@(NonEmpty.Cons (_, (Array sh _, _)) _) =+   Array.unsafeCreateWithSize (MatrixShape.Hermitian ColumnMajor sh) $+      \triSize aPtr -> do+   let n = Shape.size sh+   evalContT $ do+      uploPtr <- Call.char $ uploFromOrder ColumnMajor+      nPtr <- Call.cint n+      alphaPtr <- Call.alloca+      incPtr <- Call.cint 1+      liftIO $ do+         fill zero triSize aPtr+         forM_ xys $ \(alpha, (Array shX x, Array shY y)) ->+            withForeignPtr x $ \xPtr ->+            withForeignPtr y $ \yPtr -> do+               Call.assert+                  "Hermitian.sumRank2: non-matching x vector size" (sh==shX)+               Call.assert+                  "Hermitian.sumRank2: non-matching y vector size" (sh==shY)+               poke alphaPtr alpha+               BlasGen.hpr2 uploPtr nPtr alphaPtr xPtr incPtr yPtr incPtr aPtr+++{-+It is not strictly necessary to keep the 'order'.+It would be neither more complicated nor less efficient+to change the order via the conversion.+-}+toSquare, _toSquare ::+   (Shape.C sh, Class.Floating a) => Hermitian sh a -> Square sh a+_toSquare (Array (MatrixShape.Hermitian order sh) a) =+      Array.unsafeCreate (MatrixShape.Square order sh) $ \bPtr ->+   evalContT $ do+      let n = Shape.size sh+      aPtr <- ContT $ withForeignPtr a+      conjPtr <- conjugateToTemp (MatrixShape.triangleSize n) a+      liftIO $ do+         unpack (flipOrder order) n conjPtr bPtr -- wrong+         unpack order n aPtr bPtr++toSquare (Array (MatrixShape.Hermitian order sh) a) =+      Array.unsafeCreate (MatrixShape.Square order sh) $ \bPtr ->+   withForeignPtr a $ \aPtr ->+      unpackFull order (Shape.size sh) aPtr bPtr+++{- |+Make a temporary copy only for complex matrices.+-}+conjugateToTemp ::+   (Class.Floating a) => Int -> ForeignPtr a -> ContT r IO (Ptr a)+conjugateToTemp n =+   runCopyToTemp $+   Class.switchFloating+      (CopyToTemp $ ContT . withForeignPtr)+      (CopyToTemp $ ContT . withForeignPtr)+      (CopyToTemp $ complexConjugateToTemp n)+      (CopyToTemp $ complexConjugateToTemp n)++newtype CopyToTemp r a =+   CopyToTemp {runCopyToTemp :: ForeignPtr a -> ContT r IO (Ptr a)}++complexConjugateToTemp ::+   Class.Real a =>+   Int -> ForeignPtr (Complex a) -> ContT r IO (Ptr (Complex a))+complexConjugateToTemp n x = do+   nPtr <- Call.cint n+   xPtr <- copyToTemp n x+   incxPtr <- Call.cint 1+   liftIO $ LapackComplex.lacgv nPtr xPtr incxPtr+   return xPtr+++{- |+A^H * A+-}+covariance ::+   (Shape.C height, Shape.C width, Class.Floating a) =>+   General height width a -> Hermitian width a+covariance =+   getMap $+   Class.switchFloating+      (Map covarianceAux) (Map covarianceAux)+      (Map covarianceAux) (Map covarianceAux)++newtype Map f g a = Map {getMap :: f a -> g a}++covarianceAux ::+   (Shape.C height, Shape.C width,+    Class.Floating a, RealOf a ~ ar, Class.Real ar) =>+   General height width a -> Hermitian width a+covarianceAux (Array (MatrixShape.General order height width) a) =+      Array.unsafeCreate (MatrixShape.Hermitian order width) $ \bPtr -> do+   let n = Shape.size width+   let k = Shape.size height+   evalContT $ do+      nPtr <- Call.cint n+      kPtr <- Call.cint k+      alphaPtr <- Call.number one+      aPtr <- ContT $ withForeignPtr a+      betaPtr <- Call.number zero+      cPtr <- Call.allocaArray (n*n)+      ldcPtr <- Call.cint n++      case order of+         ColumnMajor -> do+            uploPtr <- Call.char 'U'+            transPtr <- Call.char 'C'+            ldaPtr <- Call.cint k+            liftIO $ do+               herk uploPtr transPtr+                  nPtr kPtr alphaPtr aPtr ldaPtr betaPtr cPtr ldcPtr+               pack ColumnMajor n cPtr bPtr++         RowMajor -> do+            uploPtr <- Call.char 'L'+            transPtr <- Call.char 'N'+            ldaPtr <- Call.cint n+            liftIO $ do+               herk uploPtr transPtr+                  nPtr kPtr alphaPtr aPtr ldaPtr betaPtr cPtr ldcPtr+               pack RowMajor n cPtr bPtr+++type HERK_ a =+   Ptr CChar -> Ptr CChar -> Ptr CInt -> Ptr CInt -> Ptr (RealOf a) -> Ptr a ->+   Ptr CInt -> Ptr (RealOf a) -> Ptr a -> Ptr CInt -> IO ()++newtype HERK a = HERK {getHERK :: HERK_ a}++herk :: Class.Floating a => HERK_ a+herk =+   getHERK $+   Class.switchFloating+      (HERK BlasReal.syrk)+      (HERK BlasReal.syrk)+      (HERK BlasComplex.herk)+      (HERK BlasComplex.herk)+++{- |+A^H + A+-}+addTransposed, _addTransposed ::+   (Shape.C sh, Class.Floating a) => Square sh a -> Hermitian sh a+_addTransposed (Array (MatrixShape.Square order sh) a) =+      Array.unsafeCreateWithSize (MatrixShape.Hermitian order sh) $ \bSize bPtr -> do+   let n = Shape.size sh+   evalContT $ do+      alphaPtr <- Call.number one+      incxPtr <- Call.cint 1+      aPtr <- ContT $ withForeignPtr a+      sizePtr <- Call.cint bSize+      conjPtr <- Call.allocaArray bSize+      liftIO $ do+         pack order n aPtr bPtr+         pack (flipOrder order) n aPtr conjPtr -- wrong+         lacgv sizePtr conjPtr incxPtr+         BlasGen.axpy sizePtr alphaPtr conjPtr incxPtr bPtr incxPtr++addTransposed (Array (MatrixShape.Square order sh) a) =+      Array.unsafeCreate (MatrixShape.Hermitian order sh) $ \bPtr -> do+   let n = Shape.size sh+   evalContT $ do+      alphaPtr <- Call.number one+      incxPtr <- Call.cint 1+      incnPtr <- Call.cint n+      aPtr <- ContT $ withForeignPtr a+      liftIO $ case order of+         RowMajor ->+            forPointers (rowMajorPointers n aPtr bPtr) $+               \nPtr (srcPtr,dstPtr) -> do+            BlasGen.copy nPtr srcPtr incnPtr dstPtr incxPtr+            lacgv nPtr dstPtr incxPtr+            BlasGen.axpy nPtr alphaPtr srcPtr incxPtr dstPtr incxPtr+         ColumnMajor ->+            forPointers (columnMajorPointers n aPtr bPtr) $+               \nPtr ((srcRowPtr,srcColumnPtr),dstPtr) -> do+            BlasGen.copy nPtr srcRowPtr incnPtr dstPtr incxPtr+            lacgv nPtr dstPtr incxPtr+            BlasGen.axpy nPtr alphaPtr srcColumnPtr incxPtr dstPtr incxPtr+++unpackFull :: Class.Floating a => Order -> Int -> Ptr a -> Ptr a -> IO ()+unpackFull order n packedPtr fullPtr = evalContT $ do+   incxPtr <- Call.cint 1+   incyPtr <- Call.cint n+   liftIO $ case order of+      RowMajor ->+         forPointers (rowMajorPointers n fullPtr packedPtr) $+               \nPtr (dstPtr,srcPtr) -> do+            BlasGen.copy nPtr srcPtr incxPtr dstPtr incyPtr+            lacgv nPtr dstPtr incyPtr+            BlasGen.copy nPtr srcPtr incxPtr dstPtr incxPtr+      ColumnMajor ->+         forPointers (columnMajorPointers n fullPtr packedPtr) $+               \nPtr ((dstRowPtr,dstColumnPtr),srcPtr) -> do+            BlasGen.copy nPtr srcPtr incxPtr dstRowPtr incyPtr+            lacgv nPtr dstRowPtr incyPtr+            BlasGen.copy nPtr srcPtr incxPtr dstColumnPtr incxPtr++_pack :: Class.Floating a => Order -> Int -> Ptr a -> Ptr a -> IO ()+_pack order n fullPtr packedPtr =+   evalContT $ do+      incxPtr <- Call.cint 1+      liftIO $+         case order of+            ColumnMajor ->+               forPointers (columnMajorPointers n fullPtr packedPtr) $+                  \nPtr ((_,srcPtr),dstPtr) ->+                     BlasGen.copy nPtr srcPtr incxPtr dstPtr incxPtr+            RowMajor ->+               forPointers (rowMajorPointers n fullPtr packedPtr) $+                  \nPtr (srcPtr,dstPtr) ->+                     BlasGen.copy nPtr srcPtr incxPtr dstPtr incxPtr
+ src/Numeric/LAPACK/Matrix/Multiply.hs view
@@ -0,0 +1,328 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE MultiParamTypeClasses #-}+module Numeric.LAPACK.Matrix.Multiply where++import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape+import qualified Numeric.LAPACK.Matrix.Triangular as Triangular+import qualified Numeric.LAPACK.Matrix.Hermitian as Hermitian+import qualified Numeric.LAPACK.Matrix.Square as Square+import qualified Numeric.LAPACK.Vector as Vector+import qualified Numeric.LAPACK.Private as Private+import Numeric.LAPACK.Matrix.Shape.Private+         (HeightOf, WidthOf, Order(ColumnMajor), transposeFromOrder)+import Numeric.LAPACK.Matrix.Triangular (Triangular)+import Numeric.LAPACK.Matrix.Private (General)+import Numeric.LAPACK.Vector (Vector)+import Numeric.LAPACK.Private (zero, one)++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.Internal as Array+import qualified Data.Array.Comfort.Shape as Shape+import Data.Array.Comfort.Storable.Internal (Array(Array))++import Foreign.ForeignPtr (withForeignPtr)++import Control.Monad.Trans.Cont (ContT(ContT), evalContT)+import Control.Monad.IO.Class (liftIO)+++transpose :: General height width a -> General width height a+transpose = Array.mapShape MatrixShape.transpose++multiplyVector ::+   (Shape.C height, Shape.C width, Eq width, Class.Floating a) =>+   General height width a -> Vector width a -> Vector height a+multiplyVector a x =+   let MatrixShape.General _order _height width = Array.shape a+   in if width == Array.shape x+         then multiplyVectorUnchecked a x+         else error "multiplyVector: width shapes mismatch"++multiplyVectorUnchecked ::+   (Shape.C height, Shape.C width, Class.Floating a) =>+   General height width a -> Vector width a -> Vector height a+multiplyVectorUnchecked+   (Array shape@(MatrixShape.General order height _width) a) (Array _ x) =+      Array.unsafeCreate height $ \yPtr -> do+   let (m,n) = MatrixShape.dimensions shape+   let lda = m+   evalContT $ do+      transPtr <- Call.char $ transposeFromOrder order+      mPtr <- Call.cint m+      nPtr <- Call.cint n+      alphaPtr <- Call.number one+      aPtr <- ContT $ withForeignPtr a+      ldaPtr <- Call.cint lda+      xPtr <- ContT $ withForeignPtr x+      incxPtr <- Call.cint 1+      betaPtr <- Call.number zero+      incyPtr <- Call.cint 1+      liftIO $+         BlasGen.gemv+            transPtr mPtr nPtr alphaPtr aPtr ldaPtr+            xPtr incxPtr betaPtr yPtr incyPtr++multiply ::+   (Shape.C height,+    Shape.C fuse, Eq fuse,+    Shape.C width,+    Class.Floating a) =>+   General height fuse a -> General fuse width a -> General height width a+multiply+   (Array (MatrixShape.General orderA height fuseA) a)+   (Array (MatrixShape.General orderB fuseB width) b) =+      Array.unsafeCreate (MatrixShape.General ColumnMajor height width) $+         \cPtr -> do+   Call.assert "multiply: fuse shapes mismatch" (fuseA == fuseB)+   let m = Shape.size height+   let n = Shape.size width+   let k = Shape.size fuseA+   Private.multiplyMatrix orderA orderB m k n a b cPtr+++infixl 7 <#, <#>+infixr 7 #>++class MultiplyRight shape where+   (#>) ::+      (Class.Floating a) =>+      Array shape a -> Array (WidthOf shape) a -> Array (HeightOf shape) a++class MultiplyLeft shape where+   (<#) ::+      (Class.Floating a) =>+      Array (HeightOf shape) a -> Array shape a -> Array (WidthOf shape) a++class Multiply shapeA shapeB where+   type Multiplied shapeA shapeB+   (<#>) ::+      (Class.Floating a) =>+      Array shapeA a -> Array shapeB a -> Array (Multiplied shapeA shapeB) a+++instance+   (Eq width, Shape.C width, Shape.C height) =>+      MultiplyRight (MatrixShape.General height width) where+   (#>) = multiplyVector++instance+   (Eq height, Shape.C width, Shape.C height) =>+      MultiplyLeft (MatrixShape.General height width) where+   v <# m = multiplyVector (transpose m) v+++instance+   (Eq shape, Shape.C shape) =>+      MultiplyRight (MatrixShape.Square shape) where+   m #> v = multiplyVector (Square.toGeneral m) v++instance+   (Eq shape, Shape.C shape) =>+      MultiplyLeft (MatrixShape.Square shape) where+   v <# m = multiplyVector (transpose $ Square.toGeneral m) v+++instance+   (Eq shape, Shape.C shape) =>+      MultiplyRight (MatrixShape.Hermitian shape) where+   m #> v = Hermitian.multiplyVector m v++instance+   (Eq shape, Shape.C shape) =>+      MultiplyLeft (MatrixShape.Hermitian shape) where+   v <# m = Hermitian.multiplyVector (Vector.conjugate m) v+++instance+   (MatrixShape.Uplo uplo, Eq shape, Shape.C shape) =>+      MultiplyRight (MatrixShape.Triangular uplo shape) where+   m #> v = Triangular.multiplyVectorRight m v++instance+   (MatrixShape.Uplo uplo, Eq shape, Shape.C shape) =>+      MultiplyLeft (MatrixShape.Triangular uplo shape) where+   v <# m = Triangular.multiplyVectorLeft m v+++instance+   (Shape.C heightA, Shape.C widthA, Shape.C widthB,+    widthA ~ heightB, Eq heightB) =>+      Multiply+         (MatrixShape.General heightA widthA)+         (MatrixShape.General heightB widthB) where+   type Multiplied+         (MatrixShape.General heightA widthA)+         (MatrixShape.General heightB widthB) =+            MatrixShape.General heightA widthB+   (<#>) = multiply++instance+   (Shape.C shapeA, Shape.C widthB, shapeA ~ heightB, Eq heightB) =>+      Multiply+         (MatrixShape.Square shapeA)+         (MatrixShape.General heightB widthB) where+   type Multiplied+         (MatrixShape.Square shapeA)+         (MatrixShape.General heightB widthB) =+            MatrixShape.General heightB widthB+   a <#> b = multiply (Square.toGeneral a) b++instance+   (Shape.C heightA, Shape.C widthA, widthA ~ shapeB, Eq shapeB) =>+      Multiply+         (MatrixShape.General heightA widthA)+         (MatrixShape.Square shapeB) where+   type Multiplied+         (MatrixShape.General heightA widthA)+         (MatrixShape.Square shapeB) =+            MatrixShape.General heightA widthA+   a <#> b = multiply a (Square.toGeneral b)++instance+   (Shape.C shapeA, shapeA ~ shapeB, Eq shapeB) =>+      Multiply (MatrixShape.Square shapeA) (MatrixShape.Square shapeB) where+   type Multiplied (MatrixShape.Square shapeA) (MatrixShape.Square shapeB) =+            MatrixShape.Square shapeA+   (<#>) = Square.multiply+++instance+   (Shape.C shapeA, shapeA ~ width, Eq width, Shape.C height) =>+      Multiply+         (MatrixShape.General height width)+         (MatrixShape.Hermitian shapeA)+            where+   type Multiplied+         (MatrixShape.General height width) (MatrixShape.Hermitian shapeA) =+            MatrixShape.General height width+   (<#>) = Hermitian.multiplyGeneralLeft++instance+   (Shape.C shapeA, shapeA ~ shapeB, Eq shapeB) =>+      Multiply (MatrixShape.Square shapeB) (MatrixShape.Hermitian shapeA)+         where+   type Multiplied+         (MatrixShape.Square shapeB) (MatrixShape.Hermitian shapeA) =+            MatrixShape.Square shapeA+   (<#>) = Hermitian.multiplySquareLeft++instance+   (Shape.C shapeA, shapeA ~ height, Eq height, Shape.C width) =>+      Multiply+         (MatrixShape.Hermitian shapeA)+         (MatrixShape.General height width)+            where+   type Multiplied+         (MatrixShape.Hermitian shapeA) (MatrixShape.General height width) =+            MatrixShape.General height width+   (<#>) = Hermitian.multiplyGeneralRight++instance+   (Shape.C shapeA, shapeA ~ shapeB, Eq shapeB) =>+      Multiply (MatrixShape.Hermitian shapeA) (MatrixShape.Square shapeB)+         where+   type Multiplied+         (MatrixShape.Hermitian shapeA) (MatrixShape.Square shapeB) =+            MatrixShape.Square shapeA+   (<#>) = Hermitian.multiplySquareRight++instance+   (Shape.C shapeA, shapeA ~ shapeB, Eq shapeB) =>+      Multiply (MatrixShape.Hermitian shapeA) (MatrixShape.Hermitian shapeB)+         where+   type Multiplied+         (MatrixShape.Hermitian shapeA) (MatrixShape.Hermitian shapeB) =+            MatrixShape.Square shapeA+   a <#> b = Hermitian.multiplySquareRight a (Hermitian.toSquare b)+++++instance+   (MatrixShape.Uplo uplo,+    Shape.C shapeA, shapeA ~ width, Eq width, Shape.C height) =>+      Multiply+         (MatrixShape.General height width)+         (MatrixShape.Triangular uplo shapeA)+            where+   type Multiplied+         (MatrixShape.General height width)+         (MatrixShape.Triangular uplo shapeA) =+            MatrixShape.General height width+   (<#>) = Triangular.multiplyGeneralLeft++instance+   (MatrixShape.Uplo uplo, Shape.C shapeA, shapeA ~ shapeB, Eq shapeB) =>+      Multiply (MatrixShape.Square shapeB) (MatrixShape.Triangular uplo shapeA)+         where+   type Multiplied+         (MatrixShape.Square shapeB) (MatrixShape.Triangular uplo shapeA) =+            MatrixShape.Square shapeA+   (<#>) = Triangular.multiplySquareLeft++instance+   (MatrixShape.Uplo uplo,+    Shape.C shapeA, shapeA ~ height, Eq height, Shape.C width) =>+      Multiply+         (MatrixShape.Triangular uplo shapeA)+         (MatrixShape.General height width)+            where+   type Multiplied+         (MatrixShape.Triangular uplo shapeA)+         (MatrixShape.General height width) =+            MatrixShape.General height width+   (<#>) = Triangular.multiplyGeneralRight++instance+   (MatrixShape.Uplo uplo, Shape.C shapeA, shapeA ~ shapeB, Eq shapeB) =>+      Multiply (MatrixShape.Triangular uplo shapeA) (MatrixShape.Square shapeB)+         where+   type Multiplied+         (MatrixShape.Triangular uplo shapeA) (MatrixShape.Square shapeB) =+            MatrixShape.Square shapeA+   (<#>) = Triangular.multiplySquareRight++instance+   (Shape.C shapeA, shapeA ~ shapeB, Eq shapeB,+    MultiplyTriangular uploA uploB) =>+      Multiply+         (MatrixShape.Triangular uploA shapeA)+         (MatrixShape.Triangular uploB shapeB) where+   type Multiplied+         (MatrixShape.Triangular uploA shapeA)+         (MatrixShape.Triangular uploB shapeB) =+            MultipliedTriangular uploA uploB shapeB+   (<#>) = multiplyTriangular++class MultiplyTriangular uploA uploB where+   type MultipliedTriangular uploA uploB :: * -> *+   multiplyTriangular ::+      (Class.Floating a, Shape.C shape, Eq shape) =>+      Triangular uploA shape a ->+      Triangular uploB shape a ->+      Array (MultipliedTriangular uploA uploB shape) a++instance MultiplyTriangular MatrixShape.Lower MatrixShape.Lower where+   type MultipliedTriangular MatrixShape.Lower MatrixShape.Lower =+         MatrixShape.Triangular MatrixShape.Lower+   multiplyTriangular = Triangular.multiply++instance MultiplyTriangular MatrixShape.Upper MatrixShape.Upper where+   type MultipliedTriangular MatrixShape.Upper MatrixShape.Upper =+         MatrixShape.Triangular MatrixShape.Upper+   multiplyTriangular = Triangular.multiply++instance MultiplyTriangular MatrixShape.Lower MatrixShape.Upper where+   type MultipliedTriangular MatrixShape.Lower MatrixShape.Upper =+         MatrixShape.Square+   multiplyTriangular a b =+      Square.multiply (Triangular.toSquare a) (Triangular.toSquare b)++instance MultiplyTriangular MatrixShape.Upper MatrixShape.Lower where+   type MultipliedTriangular MatrixShape.Upper MatrixShape.Lower =+         MatrixShape.Square+   multiplyTriangular a b =+      Square.multiply (Triangular.toSquare a) (Triangular.toSquare b)
+ src/Numeric/LAPACK/Matrix/Private.hs view
@@ -0,0 +1,15 @@+module Numeric.LAPACK.Matrix.Private where++import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape+import Data.Array.Comfort.Storable (Array)++import qualified Data.Array.Comfort.Shape as Shape+++type General height width = Array (MatrixShape.General height width)+++type ZeroInt = Shape.ZeroBased Int++zeroInt :: Int -> ZeroInt+zeroInt = Shape.ZeroBased
src/Numeric/LAPACK/Matrix/Shape/Private.hs view
@@ -3,6 +3,10 @@  import qualified Data.Array.Comfort.Shape as Shape +import Control.Applicative (Const(Const, getConst))++import Data.Functor.Identity (Identity(Identity), runIdentity)+import Data.List (tails) import Data.Tuple.HT (swap)  @@ -13,11 +17,15 @@ flipOrder RowMajor = ColumnMajor flipOrder ColumnMajor = RowMajor -charFromOrder :: Order -> Char-charFromOrder RowMajor = 'T'-charFromOrder ColumnMajor = 'N'+transposeFromOrder :: Order -> Char+transposeFromOrder RowMajor = 'T'+transposeFromOrder ColumnMajor = 'N'  +type family HeightOf shape+type family WidthOf shape++ data General height width =    General {       generalOrder :: Order,@@ -25,6 +33,9 @@       generalWidth :: width    } deriving (Eq, Show) +type instance HeightOf (General height width) = height+type instance WidthOf (General height width) = width+ instance (Shape.C height, Shape.C width) => Shape.C (General height width) where    type Index (General height width) = (Shape.Index height, Shape.Index width)    indices (General _ height width) = Shape.indices (height,width)@@ -63,6 +74,48 @@       ColumnMajor -> (Shape.size height, Shape.size width)  +data Square size =+   Square {+      squareOrder :: Order,+      squareSize :: size+   } deriving (Eq, Show)++type instance HeightOf (Square size) = size+type instance WidthOf (Square size) = size++generalFromSquare :: Square size -> General size size+generalFromSquare (Square order sh) = General order sh sh++transposeSquare :: Square sh -> Square sh+transposeSquare (Square order size) = Square (flipOrder order) size++instance (Shape.C size) => Shape.C (Square size) where+   type Index (Square size) = (Shape.Index size, Shape.Index size)+   indices (Square _ size) = Shape.indices (size,size)++   offset (Square RowMajor size) =+      Shape.offset (size,size)+   offset (Square ColumnMajor size) =+      Shape.offset (size,size) . swap+   uncheckedOffset (Square RowMajor size) =+      Shape.uncheckedOffset (size,size)+   uncheckedOffset (Square ColumnMajor size) =+      Shape.uncheckedOffset (size,size) . swap++   sizeOffset (Square RowMajor size) =+      Shape.sizeOffset (size,size)+   sizeOffset (Square ColumnMajor size) =+      Shape.sizeOffset (size,size) . swap+   uncheckedSizeOffset (Square RowMajor size) =+      Shape.uncheckedSizeOffset (size,size)+   uncheckedSizeOffset (Square ColumnMajor size) =+      Shape.uncheckedSizeOffset (size,size) . swap++   inBounds (Square _ size) = Shape.inBounds (size,size)+   size (Square _ size) = Shape.size (size,size)+   uncheckedSize (Square _ size) = Shape.uncheckedSize (size,size)++ data Householder height width =    Householder {       householderOrder :: Order,@@ -70,24 +123,27 @@       householderWidth :: width    } deriving (Eq, Show) +type instance HeightOf (Householder height width) = height+type instance WidthOf (Householder height width) = width+ data Reflector = Reflector deriving (Eq)-data Triangular = Triangular deriving (Eq)+data Triangle = Triangle deriving (Eq)  householderPart ::    (Shape.C height, Shape.C width) =>    Householder height width ->-   (Shape.Index height, Shape.Index width) -> Either Reflector Triangular+   (Shape.Index height, Shape.Index width) -> Either Reflector Triangle householderPart (Householder _ height width) (r,c) =    if Shape.offset height r > Shape.offset width c      then Left Reflector-     else Right Triangular+     else Right Triangle  instance    (Shape.C height, Shape.C width) =>       Shape.C (Householder height width) where     type Index (Householder height width) =-            (Either Reflector Triangular,+            (Either Reflector Triangle,              (Shape.Index height, Shape.Index width))     indices sh@(Householder _ height width) =@@ -125,3 +181,147 @@       Shape.inBounds (height,width) ix       &&       part == householderPart sh ix+++{- |+Store the upper triangular half of a real symmetric or complex Hermitian matrix.+-}+data Hermitian size =+   Hermitian {+      hermitianOrder :: Order,+      hermitianSize :: size+   } deriving (Eq, Show)++type instance HeightOf (Hermitian size) = size+type instance WidthOf (Hermitian size) = size++uploFromOrder :: Order -> Char+uploFromOrder RowMajor = 'L'+uploFromOrder ColumnMajor = 'U'++instance (Shape.C size) => Shape.C (Hermitian size) where+   type Index (Hermitian size) = (Shape.Index size, Shape.Index size)++   indices (Hermitian _ size) =+      let ixs = Shape.indices size+      in  concat $ zipWith (\r cs -> map ((,) r) cs) ixs $ tails ixs++   uncheckedOffset sh ix =+      snd $ Shape.uncheckedSizeOffset sh ix++   sizeOffset sh ix =+      if Shape.inBounds sh ix+        then Shape.uncheckedSizeOffset sh ix+        else error "Shape.Hermitian.sizeOffset: wrong matrix part"++   uncheckedSizeOffset (Hermitian RowMajor size) (rs,cs) =+      let (s,r) = Shape.uncheckedSizeOffset size rs+          c = Shape.uncheckedOffset size cs+      in  (s, triangleSize s - triangleSize (s-r) + c-r)+   uncheckedSizeOffset (Hermitian ColumnMajor size) (rs,cs) =+      let (s,r) = Shape.uncheckedSizeOffset size rs+          c = Shape.uncheckedOffset size cs+      in  (s, triangleSize c + r)++   size (Hermitian _ size) = triangleSize $ Shape.size size+   uncheckedSize (Hermitian _ size) = triangleSize $ Shape.uncheckedSize size+   inBounds (Hermitian _ size) ix@(r,c) =+      Shape.inBounds (size,size) ix+      &&+      Shape.offset size r <= Shape.offset size c+++data Triangular uplo size =+   Triangular {+      triangularUplo :: uplo,+      triangularOrder :: Order,+      triangularSize :: size+   } deriving (Eq, Show)++type instance HeightOf (Triangular uplo size) = size+type instance WidthOf (Triangular uplo size) = size++data Lower = Lower deriving (Eq, Show)+data Upper = Upper deriving (Eq, Show)++type LowerTriangular = Triangular Lower+type UpperTriangular = Triangular Upper++triangularTransposeUp :: LowerTriangular sh -> UpperTriangular sh+triangularTransposeUp (Triangular Lower order size) =+   Triangular Upper (flipOrder order) size++triangularTransposeDown :: UpperTriangular sh -> LowerTriangular sh+triangularTransposeDown (Triangular Upper order size) =+   Triangular Lower (flipOrder order) size+++class Uplo uplo where switchUplo :: f Lower -> f Upper -> f uplo+instance Uplo Lower where switchUplo f _ = f+instance Uplo Upper where switchUplo _ f = f++autoUplo :: Uplo uplo => uplo+autoUplo = runIdentity $ switchUplo (Identity Lower) (Identity Upper)++uploOrder :: Uplo uplo => uplo -> Order -> Order+uploOrder uplo order = caseUplo uplo (flipOrder order) order++getUploConst :: uplo -> Const a uplo -> a+getUploConst _ = getConst++caseUplo :: Uplo uplo => uplo -> a -> a -> a+caseUplo uplo lo up =+   getUploConst uplo $ switchUplo (Const lo) (Const up)++instance (Uplo uplo, Shape.C size) => Shape.C (Triangular uplo size) where+   type Index (Triangular uplo size) = (Shape.Index size, Shape.Index size)++   indices (Triangular uplo _ size) =+      let ixs = Shape.indices size+          rcs = concat $ zipWith (\r cs -> map ((,) r) cs) ixs $ tails ixs+      in  caseUplo uplo (map swap rcs) rcs++   uncheckedOffset sh ix =+      snd $ Shape.uncheckedSizeOffset sh ix++   sizeOffset sh ix =+      if Shape.inBounds sh ix+        then Shape.uncheckedSizeOffset sh ix+        else error "Shape.Triangular.sizeOffset: wrong matrix part"++   uncheckedSizeOffset (Triangular uplo RowMajor size) (rs,cs) =+      let (s,r) = Shape.uncheckedSizeOffset size rs+          c = Shape.uncheckedOffset size cs+      in  (s,+           caseUplo uplo+               (triangleSize r + c)+               (triangleSize s - triangleSize (s-r) + c-r))+   uncheckedSizeOffset (Triangular uplo ColumnMajor size) (rs,cs) =+      let (s,r) = Shape.uncheckedSizeOffset size rs+          c = Shape.uncheckedOffset size cs+      in  (s,+           caseUplo uplo+               (triangleSize s - triangleSize (s-c) + r-c)+               (triangleSize c + r))++   size (Triangular _ _ size) = triangleSize $ Shape.size size+   uncheckedSize (Triangular _ _ size) = triangleSize $ Shape.uncheckedSize size+   inBounds (Triangular uplo _ size) ix@(r,c) =+      Shape.inBounds (size,size) ix+      &&+      caseUplo uplo+         (Shape.offset size r >= Shape.offset size c)+         (Shape.offset size r <= Shape.offset size c)++triangleSize :: Int -> Int+triangleSize n = div (n*(n+1)) 2++triangleRoot :: Floating a => a -> a+triangleRoot size = (sqrt (8*size+1)-1)/2++triangleExtent :: String -> Int -> Int+triangleExtent name size =+   let n = round (triangleRoot (fromIntegral size :: Double))+   in if size == triangleSize n+        then n+        else error (name ++ ": no triangular number of elements")
+ src/Numeric/LAPACK/Matrix/Square.hs view
@@ -0,0 +1,202 @@+module Numeric.LAPACK.Matrix.Square (+   Square,+   size,+   toGeneral,+   fromGeneral,+   fromScalar,+   toScalar,+   fromList,+   autoFromList,++   transpose,+   adjoint,++   identity,+   identityFrom,+   diagonal,+   getDiagonal,+   trace,++   multiply,+   square,+   power,+   ) where+++import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape+import qualified Numeric.LAPACK.Vector as Vector+import qualified Numeric.LAPACK.Private as Private+import Numeric.LAPACK.Matrix.Shape.Private (Order(RowMajor, ColumnMajor))+import Numeric.LAPACK.Matrix.Private (General, ZeroInt, zeroInt)+import Numeric.LAPACK.Vector (Vector)+import Numeric.LAPACK.Private (zero, one)++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.Internal as Array+import qualified Data.Array.Comfort.Shape as Shape+import Data.Array.Comfort.Storable.Internal (Array(Array))++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)+++type Square sh = Array (MatrixShape.Square sh)++size :: Square sh a -> sh+size = MatrixShape.squareSize . Array.shape++toGeneral :: Square sh a -> General sh sh a+toGeneral (Array sh a) = Array (MatrixShape.generalFromSquare sh) a++fromGeneral :: (Eq sh) => General sh sh a -> Square sh a+fromGeneral (Array (MatrixShape.General order height width) a) =+   if height==width+     then Array (MatrixShape.Square order height) a+     else error "Square.fromGeneral: no square shape"+++fromScalar :: (Storable a) => a -> Square () a+fromScalar a =+   Array.unsafeCreate (MatrixShape.Square RowMajor ()) $ flip poke a++toScalar :: (Storable a) => Square () a -> a+toScalar (Array (MatrixShape.Square _ ()) a) =+   unsafePerformIO $ withForeignPtr a peek++fromList :: (Shape.C sh, Storable a) => sh -> [a] -> Square sh a+fromList sh =+   Array.fromList (MatrixShape.Square RowMajor sh)++autoFromList :: (Storable a) => [a] -> Square ZeroInt a+autoFromList xs =+   let n = length xs+       m = round $ sqrt (fromIntegral n :: Double)+   in if n == m*m+        then fromList (zeroInt m) xs+        else error "Square.autoFromList: no quadratic number of elements"+++transpose :: Square sh a -> Square sh a+transpose = Array.mapShape MatrixShape.transposeSquare++{- |+conjugate transpose+-}+adjoint :: (Shape.C sh, Class.Floating a) => Square sh a -> Square sh a+adjoint = transpose . Vector.conjugate+++identity :: (Shape.C sh, Class.Floating a) => sh -> Square sh a+identity = identityOrder ColumnMajor++identityFrom :: (Shape.C sh, Class.Floating a) => Square sh a -> Square sh a+identityFrom (Array (MatrixShape.Square order sh) _) = identityOrder order sh++identityOrder, _identityOrder ::+   (Shape.C sh, Class.Floating a) => Order -> sh -> Square sh a+identityOrder order sh =+   Array.unsafeCreate (MatrixShape.Square order sh) $ \aPtr ->+   evalContT $ do+      uploPtr <- Call.char 'A'+      nPtr <- Call.cint $ Shape.size sh+      alphaPtr <- Call.number zero+      betaPtr <- Call.number one+      liftIO $ LapackGen.laset uploPtr nPtr nPtr alphaPtr betaPtr aPtr nPtr++_identityOrder order sh =+   Array.unsafeCreateWithSize (MatrixShape.Square order sh) $ \blockSize yPtr ->+   evalContT $ do+      nPtr <- Call.alloca+      xPtr <- Call.number zero+      incxPtr <- Call.cint 0+      incyPtr <- Call.cint 1+      liftIO $ do+         poke nPtr $ fromIntegral blockSize+         BlasGen.copy nPtr xPtr incxPtr yPtr incyPtr+         let n = fromIntegral $ Shape.size sh+         poke nPtr n+         poke xPtr one+         poke incyPtr (n+1)+         BlasGen.copy nPtr xPtr incxPtr yPtr incyPtr++diagonal :: (Shape.C sh, Class.Floating a) => Vector sh a -> Square sh a+diagonal (Array sh x) =+   Array.unsafeCreateWithSize (MatrixShape.Square ColumnMajor sh) $+      \blockSize yPtr ->+   evalContT $ do+      nPtr <- Call.alloca+      xPtr <- ContT $ withForeignPtr x+      zPtr <- Call.number zero+      incxPtr <- Call.cint 1+      incyPtr <- Call.cint 1+      inczPtr <- Call.cint 0+      liftIO $ do+         poke nPtr $ fromIntegral blockSize+         BlasGen.copy nPtr zPtr inczPtr yPtr incyPtr+         let n = fromIntegral $ Shape.size sh+         poke nPtr n+         poke incyPtr (n+1)+         BlasGen.copy nPtr xPtr incxPtr yPtr incyPtr++getDiagonal :: (Shape.C sh, Class.Floating a) => Square sh a -> Vector sh a+getDiagonal (Array (MatrixShape.Square _ sh) x) =+   Array.unsafeCreateWithSize sh $ \n yPtr -> evalContT $ do+      nPtr <- Call.cint n+      xPtr <- ContT $ withForeignPtr x+      incxPtr <- Call.cint (n+1)+      incyPtr <- Call.cint 1+      liftIO $ BlasGen.copy nPtr xPtr incxPtr yPtr incyPtr++trace :: (Shape.C sh, Class.Floating a) => Square sh a -> a+trace (Array (MatrixShape.Square _ sh) x) = unsafePerformIO $ do+   let n = Shape.size sh+   withForeignPtr x $ \xPtr -> Private.sum n xPtr (n+1)+++multiply ::+   (Shape.C sh, Eq sh, Class.Floating a) =>+   Square sh a -> Square sh a -> Square sh a+multiply+   (Array (MatrixShape.Square orderA shA) a)+   (Array (MatrixShape.Square orderB shB) b) =+      Array.unsafeCreate (MatrixShape.Square ColumnMajor shA) $ \cPtr -> do+   Call.assert "Square.multiply: shapes mismatch" (shA == shB)+   let n = Shape.size shA+   Private.multiplyMatrix orderA orderB n n n a b cPtr++square :: (Shape.C sh, Class.Floating a) => Square sh a -> Square sh a+square a = multiplyCommutativeUnchecked a a++power ::+   (Shape.C sh, Class.Floating a) =>+   Integer -> Square sh a -> Square sh a+power n a =+   powerAssociative multiplyCommutativeUnchecked (identityFrom a) a n++{-+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+   (Array shape@(MatrixShape.Square  order  sh) a)+   (Array       (MatrixShape.Square _order _sh) b) =+      Array.unsafeCreate shape $ \cPtr ->+   let n = Shape.size sh+       (at,bt) =+         case order of+            ColumnMajor -> (a,b)+            RowMajor -> (b,a)+   in  Private.multiplyMatrix ColumnMajor ColumnMajor n n n at bt cPtr
+ src/Numeric/LAPACK/Matrix/Triangular.hs view
@@ -0,0 +1,315 @@+{-# LANGUAGE TypeFamilies #-}+module Numeric.LAPACK.Matrix.Triangular (+   Triangular, MatrixShape.Uplo(..),+   Upper, Lower,+   fromList, autoFromList,+   lowerFromList, autoLowerFromList,+   upperFromList, autoUpperFromList,+   identity,+   diagonal,+   getDiagonal,+   transposeUp, transposeDown,+   adjointUp, adjointDown,++   toSquare,++   multiplyVectorLeft,+   multiplyVectorRight,+   square,+   multiply,+   multiplySquareLeft,+   multiplyGeneralLeft,+   multiplySquareRight,+   multiplyGeneralRight,+   ) where++import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape+import qualified Numeric.LAPACK.Vector as Vector+import Numeric.LAPACK.Matrix.Triangular.Private+         (diagonalPointers, pack, unpack, unpackZero, unpackToTemp)+import Numeric.LAPACK.Matrix.Shape.Private+         (Order(RowMajor,ColumnMajor),+          flipOrder, transposeFromOrder, uploFromOrder, uploOrder)+import Numeric.LAPACK.Matrix.Square (Square)+import Numeric.LAPACK.Matrix.Private (General, ZeroInt, zeroInt)+import Numeric.LAPACK.Vector (Vector)+import Numeric.LAPACK.Private (fill, zero, one, copyBlock)++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.Internal as Array+import qualified Data.Array.Comfort.Shape as Shape+import Data.Array.Comfort.Storable.Internal (Array(Array))++import Foreign.ForeignPtr (ForeignPtr, withForeignPtr)+import Foreign.Ptr (Ptr)+import Foreign.Storable (Storable, poke, peek)++import Control.Monad.Trans.Cont (ContT(ContT), evalContT)+import Control.Monad.IO.Class (liftIO)++import Data.Foldable (forM_)+++type Triangular uplo sh = Array (MatrixShape.Triangular uplo sh)++type Lower sh = Array (MatrixShape.LowerTriangular sh)+type Upper sh = Array (MatrixShape.UpperTriangular sh)++transposeUp :: Lower sh a -> Upper sh a+transposeUp (Array sh a) =+   Array (MatrixShape.triangularTransposeUp sh) a++transposeDown :: Upper sh a -> Lower sh a+transposeDown (Array sh a) =+   Array (MatrixShape.triangularTransposeDown sh) a++adjointUp :: (Shape.C sh, Class.Floating a) => Lower sh a -> Upper sh a+adjointUp = Vector.conjugate . transposeUp++adjointDown :: (Shape.C sh, Class.Floating a) => Upper sh a -> Lower sh a+adjointDown = Vector.conjugate . transposeDown+++fromList ::+   (MatrixShape.Uplo uplo, Shape.C sh, Storable a) =>+   Order -> sh -> [a] -> Triangular uplo sh a+fromList order sh =+   Array.fromList (MatrixShape.Triangular MatrixShape.autoUplo order sh)++lowerFromList :: (Shape.C sh, Storable a) => Order -> sh -> [a] -> Lower sh a+lowerFromList = fromList++upperFromList :: (Shape.C sh, Storable a) => Order -> sh -> [a] -> Upper sh a+upperFromList = fromList+++autoFromList ::+   (MatrixShape.Uplo uplo, Storable a) =>+   Order -> [a] -> Triangular uplo ZeroInt a+autoFromList order xs =+   fromList order+      (zeroInt $ MatrixShape.triangleExtent "Triangular.autoFromList" $+       length xs)+      xs++autoLowerFromList :: (Storable a) => Order -> [a] -> Lower ZeroInt a+autoLowerFromList = autoFromList++autoUpperFromList :: (Storable a) => Order -> [a] -> Upper ZeroInt a+autoUpperFromList = autoFromList+++toSquare ::+   (MatrixShape.Uplo uplo, Shape.C sh, Class.Floating a) =>+   Triangular uplo sh a -> Square sh a+toSquare (Array (MatrixShape.Triangular uplo order sh) a) =+   Array.unsafeCreate (MatrixShape.Square order sh) $ \bPtr ->+      withForeignPtr a $ \aPtr ->+         unpackZero (uploOrder uplo order) (Shape.size sh) aPtr bPtr+++identity ::+   (MatrixShape.Uplo uplo, Shape.C sh, Class.Floating a) =>+   Order -> sh -> Triangular uplo sh a+identity order sh =+   let (realOrder, uplo) = autoUploOrder order+   in Array.unsafeCreate (MatrixShape.Triangular uplo order sh) $ \aPtr -> do+      let n = Shape.size sh+      fill zero (MatrixShape.triangleSize n) aPtr+      forM_ (diagonalPointers realOrder n aPtr aPtr) $ flip poke one . snd++diagonal ::+   (MatrixShape.Uplo uplo, Shape.C sh, Class.Floating a) =>+   Order -> Vector sh a -> Triangular uplo sh a+diagonal order (Array sh x) =+   let (realOrder, uplo) = autoUploOrder order+   in Array.unsafeCreate (MatrixShape.Triangular uplo order sh) $ \aPtr -> do+      let n = Shape.size sh+      fill zero (MatrixShape.triangleSize n) aPtr+      withForeignPtr x $ \xPtr ->+         forM_ (diagonalPointers realOrder n xPtr aPtr) $+            \(srcPtr,dstPtr) -> poke dstPtr =<< peek srcPtr++getDiagonal ::+   (MatrixShape.Uplo uplo, Shape.C sh, Class.Floating a) =>+   Triangular uplo sh a -> Vector sh a+getDiagonal (Array (MatrixShape.Triangular uplo order sh) a) =+      Array.unsafeCreate sh $ \xPtr -> do+   withForeignPtr a $ \aPtr ->+      mapM_+         (\(dstPtr,srcPtr) -> poke dstPtr =<< peek srcPtr)+         (diagonalPointers (uploOrder uplo order) (Shape.size sh) xPtr aPtr)+++multiplyVectorLeft, multiplyVectorRight ::+   (MatrixShape.Uplo uplo, Shape.C sh, Eq sh, Class.Floating a) =>+   Triangular uplo sh a -> Vector sh a -> Vector sh a+multiplyVectorLeft = multiplyVector True+multiplyVectorRight = multiplyVector False++multiplyVector ::+   (MatrixShape.Uplo uplo, Shape.C sh, Eq sh, Class.Floating a) =>+   Bool -> Triangular uplo sh a -> Vector sh a -> Vector sh a+multiplyVector transp+   (Array (MatrixShape.Triangular uplo order shA) a) (Array shX x) =+      Array.unsafeCreate shX $ \yPtr -> do+   Call.assert "Triangular.multiplyVector: width shapes mismatch" (shA == shX)+   let n = Shape.size shA+   evalContT $ do+      uploPtr <- Call.char $ uploFromOrder $ uploOrder uplo order+      transPtr <-+         Call.char $ transposeFromOrder $+         (if transp then flipOrder else id) order+      diagPtr <- Call.char 'N'+      nPtr <- Call.cint n+      aPtr <- ContT $ withForeignPtr a+      xPtr <- ContT $ withForeignPtr x+      incyPtr <- Call.cint 1+      liftIO $ do+         copyBlock n xPtr yPtr+         BlasGen.tpmv uploPtr transPtr diagPtr nPtr aPtr yPtr incyPtr+++square ::+   (MatrixShape.Uplo uplo, Shape.C sh, Eq sh, Class.Floating a) =>+   Triangular uplo sh a -> Triangular uplo sh a+square+   (Array shape@(MatrixShape.Triangular uplo order sh) a) =+      Array.unsafeCreate shape $ \bpPtr -> do+   let n = Shape.size sh+   evalContT $ do+      sidePtr <- Call.char 'L'+      let realOrder = uploOrder uplo order+      uploPtr <- Call.char $ uploFromOrder realOrder+      transPtr <- Call.char 'N'+      diagPtr <- Call.char 'N'+      nPtr <- Call.cint n+      let ldPtr = nPtr+      aPtr <- unpackToTemp (unpack realOrder) n a+      bPtr <- unpackToTemp (unpackZero realOrder) n a+      alphaPtr <- Call.number one+      liftIO $ do+         BlasGen.trmm sidePtr uploPtr transPtr diagPtr+            nPtr nPtr alphaPtr aPtr ldPtr bPtr ldPtr+         pack realOrder n bPtr bpPtr++multiply ::+   (MatrixShape.Uplo uplo, Shape.C sh, Eq sh, Class.Floating a) =>+   Triangular uplo sh a -> Triangular uplo sh a -> Triangular uplo sh a+multiply+   (Array        (MatrixShape.Triangular uploA orderA shA) a)+   (Array shapeB@(MatrixShape.Triangular uploB orderB shB) b) =+      Array.unsafeCreate shapeB $ \cpPtr -> do+   Call.assert "Triangular.multiply: width shapes mismatch" (shA == shB)+   let n = Shape.size shA+   evalContT $ do+      let (side,trans) =+            case orderB of+               ColumnMajor -> ('L', orderA)+               RowMajor -> ('R', flipOrder orderA)+      sidePtr <- Call.char side+      let realOrderA = uploOrder uploA orderA+      let realOrderB = uploOrder uploB orderB+      uploPtr <- Call.char $ uploFromOrder realOrderA+      transPtr <- Call.char $ transposeFromOrder trans+      diagPtr <- Call.char 'N'+      nPtr <- Call.cint n+      let ldPtr = nPtr+      aPtr <- unpackToTemp (unpack realOrderA) n a+      bPtr <- unpackToTemp (unpackZero realOrderB) n b+      alphaPtr <- Call.number one+      liftIO $ do+         BlasGen.trmm sidePtr uploPtr transPtr diagPtr+            nPtr nPtr alphaPtr aPtr ldPtr bPtr ldPtr+         pack realOrderB n bPtr cpPtr+++multiplySquareLeft ::+   (MatrixShape.Uplo uplo, Shape.C sh, Eq sh, Class.Floating a) =>+   Square sh a -> Triangular uplo sh a -> Square sh a+multiplySquareLeft+   (Array shapeB@(MatrixShape.Square orderB shB) b)+   (Array        (MatrixShape.Triangular uploA orderA shA) a) =+      Array.unsafeCreate shapeB $ \cPtr -> do+   Call.assert "Triangular.multiplySquareLeft: shapes mismatch" (shA == shB)+   let n = Shape.size shB+   MatrixShape.caseUplo uploA+      (multiplyAux MatrixShape.Upper)+      (multiplyAux MatrixShape.Lower)+      (flipOrder orderA) n a (flipOrder orderB) n b cPtr++multiplyGeneralLeft ::+   (MatrixShape.Uplo uplo,+    Shape.C height, Shape.C width, Eq width, Class.Floating a) =>+   General height width a -> Triangular uplo width a -> General height width a+multiplyGeneralLeft+   (Array shapeB@(MatrixShape.General orderB height width) b)+   (Array        (MatrixShape.Triangular uploA orderA shA) a) =+      Array.unsafeCreate shapeB $ \cPtr -> do+   Call.assert "Triangular.multiplyGeneralLeft: shapes mismatch" (shA == width)+   MatrixShape.caseUplo uploA+      (multiplyAux MatrixShape.Upper)+      (multiplyAux MatrixShape.Lower)+      (flipOrder orderA) (Shape.size width) a+      (flipOrder orderB) (Shape.size height) b cPtr++multiplySquareRight ::+   (MatrixShape.Uplo uplo, Shape.C sh, Eq sh, Class.Floating a) =>+   Triangular uplo sh a -> Square sh a -> Square sh a+multiplySquareRight+   (Array        (MatrixShape.Triangular uploA orderA shA) a)+   (Array shapeB@(MatrixShape.Square orderB shB) b) =+      Array.unsafeCreate shapeB $ \cPtr -> do+   Call.assert "Triangular.multiplySquareRight: shapes mismatch" (shA == shB)+   let n = Shape.size shB+   multiplyAux uploA orderA n a orderB n b cPtr++multiplyGeneralRight ::+   (MatrixShape.Uplo uplo,+    Shape.C height, Eq height, Shape.C width, Class.Floating a) =>+   Triangular uplo height a -> General height width a -> General height width a+multiplyGeneralRight+   (Array        (MatrixShape.Triangular uploA orderA shA) a)+   (Array shapeB@(MatrixShape.General orderB height width) b) =+      Array.unsafeCreate shapeB $ \cPtr -> do+   Call.assert "Triangular.multiplyGeneralRight: shapes mismatch"+      (shA == height)+   multiplyAux+      uploA orderA (Shape.size height) a orderB (Shape.size width) b cPtr++multiplyAux ::+   (MatrixShape.Uplo uplo, Class.Floating a) =>+   uplo ->+   Order -> Int -> ForeignPtr a ->+   Order -> Int -> ForeignPtr a -> Ptr a -> IO ()+multiplyAux uploA orderA m0 a orderB n0 b cPtr =+   evalContT $ do+      let (side,trans,(m,n)) =+            case orderB of+               ColumnMajor -> ('L', orderA, (m0,n0))+               RowMajor -> ('R', flipOrder orderA, (n0,m0))+      sidePtr <- Call.char side+      let realOrderA = uploOrder uploA orderA+      uploPtr <- Call.char $ uploFromOrder realOrderA+      transPtr <- Call.char $ transposeFromOrder trans+      diagPtr <- Call.char 'N'+      mPtr <- Call.cint m+      nPtr <- Call.cint n+      alphaPtr <- Call.number one+      aPtr <- unpackToTemp (unpack realOrderA) m0 a+      ldaPtr <- Call.cint m0+      bPtr <- ContT $ withForeignPtr b+      ldbPtr <- Call.cint m+      liftIO $ do+         copyBlock (m0*n0) bPtr cPtr+         BlasGen.trmm sidePtr uploPtr transPtr diagPtr+            mPtr nPtr alphaPtr aPtr ldaPtr cPtr ldbPtr+++autoUploOrder :: MatrixShape.Uplo uplo => Order -> (Order, uplo)+autoUploOrder order =+   case MatrixShape.autoUplo of+      uplo -> (uploOrder uplo order, uplo)
+ src/Numeric/LAPACK/Matrix/Triangular/Private.hs view
@@ -0,0 +1,121 @@+module Numeric.LAPACK.Matrix.Triangular.Private where++import Numeric.LAPACK.Matrix.Shape.Private+         (Order(RowMajor,ColumnMajor), flipOrder, uploFromOrder, triangleSize)+import Numeric.LAPACK.Private (pointerSeq, copyToTemp, lacgv, fill, zero)++import qualified Numeric.LAPACK.FFI.Generic as LapackGen+import qualified Numeric.Netlib.Utility as Call+import qualified Numeric.Netlib.Class as Class++import Foreign.Marshal.Alloc (alloca)+import Foreign.Marshal.Array (advancePtr)+import Foreign.C.Types (CInt)+import Foreign.ForeignPtr (ForeignPtr, withForeignPtr)+import Foreign.Ptr (Ptr)+import Foreign.Storable (Storable, poke)++import Control.Monad.Trans.Cont (ContT(ContT), evalContT)+import Control.Monad.IO.Class (liftIO)++import Data.Foldable (forM_)+++diagonalPointers ::+   (Storable a, Storable ar) =>+   Order -> Int -> Ptr ar -> Ptr a -> [(Ptr ar, Ptr a)]+diagonalPointers order n xPtr aPtr =+   take n $ zip (pointerSeq 1 xPtr) $ scanl advancePtr aPtr $+   case order of+      RowMajor -> iterate pred n+      ColumnMajor -> iterate succ 2+++columnMajorPointers ::+   (Storable a) => Int -> Ptr a -> Ptr a -> [(Int, ((Ptr a, Ptr a), Ptr a))]+columnMajorPointers n fullPtr packedPtr =+   let ds = iterate succ 1+   in  take n $ zip ds $+       zip+         (zip (pointerSeq 1 fullPtr) (pointerSeq n fullPtr))+         (scanl advancePtr packedPtr ds)++rowMajorPointers ::+   (Storable a) => Int -> Ptr a -> Ptr a -> [(Int, (Ptr a, Ptr a))]+rowMajorPointers n fullPtr packedPtr =+   let ds = iterate pred n+   in  take n $ zip ds $+       zip (pointerSeq (n+1) fullPtr) (scanl advancePtr packedPtr ds)+++forPointers :: [(Int, a)] -> (Ptr CInt -> a -> IO ()) -> IO ()+forPointers xs act =+   alloca $ \nPtr ->+   forM_ xs $ \(d,ptrs) -> do+      poke nPtr $ fromIntegral d+      act nPtr ptrs+++copyTriangleToTemp ::+   Class.Floating a =>+   Order -> Int -> ForeignPtr a -> ContT r IO (Ptr a)+copyTriangleToTemp order n a = do+   let aSize = triangleSize n+   apPtr <- copyToTemp aSize a+   liftIO $ evalContT $ do+      aSizePtr <- Call.cint aSize+      incPtr <- Call.cint 1+      case order of+         RowMajor -> liftIO $ lacgv aSizePtr apPtr incPtr+         ColumnMajor -> return ()+   return apPtr+++unpackToTemp ::+   Storable a =>+   (Int -> Ptr a -> Ptr a -> IO ()) ->+   Int -> ForeignPtr a -> ContT r IO (Ptr a)+unpackToTemp f n a = do+   apPtr <- ContT $ withForeignPtr a+   aPtr <- Call.allocaArray (n*n)+   liftIO $ f n apPtr aPtr+   return aPtr+++unpack :: Class.Floating a => Order -> Int -> Ptr a -> Ptr a -> IO ()+unpack order n packedPtr fullPtr =+   evalContT $ do+      uploPtr <- Call.char $ uploFromOrder order+      nPtr <- Call.cint n+      ldaPtr <- Call.cint n+      liftIO $ withInfo $ LapackGen.tpttr uploPtr nPtr packedPtr fullPtr ldaPtr++pack :: Class.Floating a => Order -> Int -> Ptr a -> Ptr a -> IO ()+pack order n fullPtr packedPtr =+   evalContT $ do+      uploPtr <- Call.char $ uploFromOrder order+      nPtr <- Call.cint n+      ldaPtr <- Call.cint n+      liftIO $ withInfo $ LapackGen.trttp uploPtr nPtr fullPtr ldaPtr packedPtr+++unpackZero, _unpackZero ::+   Class.Floating a => Order -> Int -> Ptr a -> Ptr a -> IO ()+_unpackZero order n packedPtr fullPtr = do+   fill zero (n*n) fullPtr+   unpack order n packedPtr fullPtr++unpackZero order n packedPtr fullPtr = do+   fillTriangle zero (flipOrder order) n fullPtr+   unpack order n packedPtr fullPtr++fillTriangle :: Class.Floating a => a -> Order -> Int -> Ptr a -> IO ()+fillTriangle z order n aPtr = evalContT $ do+   uploPtr <- Call.char $ uploFromOrder order+   nPtr <- Call.cint n+   zPtr <- Call.number z+   liftIO $ LapackGen.laset uploPtr nPtr nPtr zPtr zPtr aPtr nPtr+++withInfo :: (Ptr CInt -> IO ()) -> IO ()+withInfo = alloca
+ src/Numeric/LAPACK/Orthogonal.hs view
@@ -0,0 +1,372 @@+{-# LANGUAGE TypeFamilies #-}+module Numeric.LAPACK.Orthogonal (+   leastSquares,+   minimumNorm,+   leastSquaresMinimumNorm,+   pseudoInverseRCond,++   Householder,+   householder,+   householderDecompose,+   householderDeterminant,+   determinant,+   householderExtractQ,+   householderExtractR,+   orthogonalComplement,+   ) where++import qualified Numeric.LAPACK.Matrix.Square as Square+import Numeric.LAPACK.Matrix+         (General, ZeroInt, zeroInt, transpose, identity, dropColumns)+import Numeric.LAPACK.Matrix.Square (Square)++import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape+import qualified Numeric.LAPACK.Private as Private+import Numeric.LAPACK.Matrix.Shape.Private+         (Order(RowMajor, ColumnMajor), transposeFromOrder)+import Numeric.LAPACK.Vector (Vector)+import Numeric.LAPACK.Private+         (RealOf, zero, fill,+          copySubMatrix, copyBlock, copyToTemp,+          copyToColumnMajor, copyToSubColumnMajor,+          withAutoWorkspaceInfo, allocArray, allocHigherArray)++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.Internal as Array+import qualified Data.Array.Comfort.Shape as Shape+import Data.Array.Comfort.Storable.Internal (Array(Array))++import System.IO.Unsafe (unsafePerformIO)++import Foreign.Marshal.Array (advancePtr)+import Foreign.C.Types (CInt)+import Foreign.ForeignPtr (withForeignPtr)+import Foreign.Ptr (Ptr)+import Foreign.Storable (peek)++import Control.Monad.Trans.Cont (ContT(ContT), evalContT)+import Control.Monad.IO.Class (liftIO)+import Control.Applicative ((<$>))++import Data.Complex (Complex)+import Data.Tuple.HT (mapSnd)+++{- |+If @x = leastSquares a b@+then @x@ minimizes @Vector.norm2 (multiply a x `sub` b)@.++Precondition: @a@ must have full rank and @height a >= width a@.+-}+leastSquares ::+   (Shape.C height, Eq height, Shape.C width, Shape.C nrhs, Class.Floating a) =>+   General height width a -> General height nrhs a -> General width nrhs a+leastSquares+   (Array shapeA@(MatrixShape.General orderA heightA widthA) a)+   (Array        (MatrixShape.General orderB heightB widthB) b) =+      Array.unsafeCreate (MatrixShape.General ColumnMajor widthA widthB) $+         \xPtr -> do+   Call.assert "leastSquares: height shapes mismatch" (heightA == heightB)+   Call.assert "leastSquares: height of 'a' must be at least the width"+      (Shape.size heightA >= Shape.size widthA)+   let (m,n) = MatrixShape.dimensions shapeA+   let lda = m+   let nrhs = Shape.size widthB+   let ldb = Shape.size heightB+   let ldx = Shape.size widthA+   evalContT $ do+      transPtr <- Call.char $ transposeFromOrder orderA+      mPtr <- Call.cint m+      nPtr <- Call.cint n+      nrhsPtr <- Call.cint nrhs+      aPtr <- copyToTemp (Shape.size shapeA) a+      ldaPtr <- Call.cint lda+      bPtr <- ContT $ withForeignPtr b+      ldbPtr <- Call.cint ldb+      let bSize = Shape.size (heightB,widthB)+      btmpPtr <- Call.allocaArray bSize+      liftIO $ copyToColumnMajor orderB ldb nrhs bPtr btmpPtr+      liftIO $ withAutoWorkspaceInfo "gels" $+         LapackGen.gels transPtr+            mPtr nPtr nrhsPtr aPtr ldaPtr btmpPtr ldbPtr+      liftIO $ copySubMatrix ldx nrhs ldb btmpPtr ldx xPtr++{- |+The vector @x@ with @x = minimumNorm a b@+is the vector with minimal @Vector.norm2 x@+that satisfies @multiply a x == b@.++Precondition: @a@ must have full rank and @height a <= width a@.+-}+minimumNorm ::+   (Shape.C height, Eq height, Shape.C width, Shape.C nrhs, Class.Floating a) =>+   General height width a -> General height nrhs a -> General width nrhs a+minimumNorm+   (Array shapeA@(MatrixShape.General orderA heightA widthA) a)+   (Array        (MatrixShape.General orderB heightB widthB) b) =+      Array.unsafeCreate (MatrixShape.General ColumnMajor widthA widthB) $+         \xPtr -> do+   Call.assert "minimumNorm: height shapes mismatch" (heightA == heightB)+   Call.assert "minimumNorm: width of 'a' must be at least the height"+      (Shape.size widthA >= Shape.size heightA)+   let (m,n) = MatrixShape.dimensions shapeA+   let lda = m+   let nrhs = Shape.size widthB+   let ldb = Shape.size heightB+   let ldx = Shape.size widthA+   evalContT $ do+      transPtr <- Call.char $ transposeFromOrder orderA+      mPtr <- Call.cint m+      nPtr <- Call.cint n+      nrhsPtr <- Call.cint nrhs+      aPtr <- copyToTemp (Shape.size shapeA) a+      ldaPtr <- Call.cint lda+      bPtr <- ContT $ withForeignPtr b+      ldxPtr <- Call.cint ldx+      liftIO $ copyToSubColumnMajor orderB ldb nrhs bPtr ldx xPtr+      liftIO $ withAutoWorkspaceInfo "gels" $+         LapackGen.gels transPtr+            mPtr nPtr nrhsPtr aPtr ldaPtr xPtr ldxPtr++{- |+If @x = leastSquaresMinimumNorm a b@+then @x@ is the vector with minimum @Vector.norm2 x@+that minimizes @Vector.norm2 (multiply a x `sub` b)@.++Matrix @a@ can have any rank+but you must specify the reciprocal condition of the rank-truncated matrix.+-}+leastSquaresMinimumNorm ::+   (Shape.C height, Eq height, Shape.C width, Shape.C nrhs, Class.Floating a) =>+   RealOf a ->+   General height width a -> General height nrhs a ->+   (Int, General width nrhs a)+leastSquaresMinimumNorm rcond+   (Array (MatrixShape.General orderA heightA widthA) a)+   (Array (MatrixShape.General orderB heightB widthB) b) =+      unsafePerformIO $ do+   Call.assert "leastSquaresMinimumNorm: height shapes mismatch"+      (heightA == heightB)+   let shapeX = MatrixShape.General ColumnMajor widthA widthB+   let m = Shape.size heightA+   let n = Shape.size widthA+   let nrhs = Shape.size widthB+   let aSize = m*n+   let lda = m+   evalContT $ do+      aPtr <- ContT $ withForeignPtr a+      atmpPtr <- Call.allocaArray aSize+      liftIO $ copyToColumnMajor orderA m n aPtr atmpPtr+      ldaPtr <- Call.cint lda+      (x,(tmpPtr,ldtmp)) <- allocHigherArray shapeX m n nrhs+      ldtmpPtr <- Call.cint ldtmp+      bPtr <- ContT $ withForeignPtr b+      liftIO $ copyToSubColumnMajor orderB m nrhs bPtr ldtmp tmpPtr+      jpvtPtr <- Call.allocaArray n+      rankPtr <- Call.alloca+      gelsy m n nrhs atmpPtr ldaPtr tmpPtr ldtmpPtr jpvtPtr rcond rankPtr+      rank <- liftIO $ fromIntegral <$> peek rankPtr+      return (rank, x)+++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 {getGELSY :: GELSY_ r (RealOf a) a}++gelsy :: (Class.Floating a) => GELSY_ r (RealOf a) a+gelsy =+   getGELSY $+   Class.switchFloating+      (GELSY gelsyReal)+      (GELSY gelsyReal)+      (GELSY gelsyComplex)+      (GELSY gelsyComplex)++gelsyReal :: (Class.Real a) => GELSY_ r a a+gelsyReal m n nrhs aPtr ldaPtr bPtr ldbPtr jpvtPtr rcond rankPtr = do+   mPtr <- Call.cint m+   nPtr <- Call.cint n+   nrhsPtr <- Call.cint nrhs+   rcondPtr <- Call.real rcond+   liftIO $ withAutoWorkspaceInfo "gelsy" $+      LapackReal.gelsy mPtr nPtr nrhsPtr+         aPtr ldaPtr bPtr ldbPtr jpvtPtr rcondPtr rankPtr++gelsyComplex :: (Class.Real a) => GELSY_ r a (Complex a)+gelsyComplex m n nrhs aPtr ldaPtr bPtr ldbPtr jpvtPtr rcond rankPtr = do+   mPtr <- Call.cint m+   nPtr <- Call.cint n+   nrhsPtr <- Call.cint nrhs+   rcondPtr <- Call.real rcond+   rworkPtr <- Call.allocaArray (2*n)+   liftIO $ withAutoWorkspaceInfo "gelsy" $ \workPtr lworkPtr infoPtr ->+      LapackComplex.gelsy mPtr nPtr nrhsPtr+         aPtr ldaPtr bPtr ldbPtr jpvtPtr rcondPtr rankPtr+         workPtr lworkPtr rworkPtr infoPtr+++pseudoInverseRCond ::+   (Shape.C height, Eq height, Shape.C width, Eq width, Class.Floating a) =>+   RealOf a -> General height width a -> (Int, General width height a)+pseudoInverseRCond rcond a =+   let (MatrixShape.General _ height width) = Array.shape a+   in if Shape.size height < Shape.size width+         then leastSquaresMinimumNorm rcond a $ identity height+         else mapSnd transpose $+              leastSquaresMinimumNorm rcond (transpose a) $+              identity width+++type Householder height width = Array (MatrixShape.Householder height width)++{-+@(q,r) = householder a@+means that @q@ is unitary and @r@ is upper triangular and @a = multiply q r@.+-}+householder ::+   (Shape.C height, Shape.C width, Eq width, Class.Floating a) =>+   General height width a ->+   (Square height a, General height width a)+householder a =+   let hh = householderDecompose a+   in  (householderExtractQ hh, householderExtractR $ snd hh)++householderDecompose ::+   (Shape.C height, Shape.C width, Class.Floating a) =>+   General height width a -> (Vector width a, Householder height width a)+householderDecompose (Array shape@(MatrixShape.General order height width) a) =+   unsafePerformIO $ do++   let (m,n) = MatrixShape.dimensions shape+   let lda = m+   let mn = min m n+   evalContT $ do+      mPtr <- Call.cint m+      nPtr <- Call.cint n+      aPtr <- ContT $ withForeignPtr a+      ldaPtr <- Call.cint lda+      (qr,qrPtr) <- allocArray $ MatrixShape.Householder order height width+      liftIO $ copyBlock (m*n) aPtr qrPtr+      (tau,tauPtr) <- allocArray width+      liftIO $ fill zero (n-mn) (advancePtr tauPtr mn)+      liftIO $+         case order of+            RowMajor ->+               withAutoWorkspaceInfo "gelqf" $+                  LapackGen.gelqf mPtr nPtr qrPtr ldaPtr tauPtr+            ColumnMajor ->+               withAutoWorkspaceInfo "geqrf" $+                  LapackGen.geqrf mPtr nPtr qrPtr ldaPtr tauPtr+      return (tau, qr)++householderDeterminant ::+   (Shape.C height, Shape.C width, Class.Floating a) =>+   Householder height width a -> a+householderDeterminant+      (Array (MatrixShape.Householder order height width) a) =+   let m = Shape.size height+       n = Shape.size width+       k = case order of RowMajor -> n; ColumnMajor -> m+   in unsafePerformIO $+      withForeignPtr a $ \aPtr ->+      Private.product (min m n) aPtr (k+1)+++{-|+Generalized determinant - works also for non-square matrices.+In contrast to the square root of the Gramian determinant+it has the proper sign.+-}+determinant ::+   (Shape.C height, Shape.C width, Eq a, Class.Floating a) =>+   General height width a -> a+determinant a =+   let (tau,hh) = householderDecompose a+   in  foldl (\x _ -> negate x)+         (householderDeterminant hh)+         (takeWhile (/=zero) $ Array.toList tau)+++householderExtractQ ::+   (Shape.C height, Shape.C width, Eq width, Class.Floating a) =>+   (Vector width a, Householder height width a) -> Square height a+householderExtractQ+   (Array widthTau tau,+    Array (MatrixShape.Householder order height width) qr) =++   Array.unsafeCreate (MatrixShape.Square order height) $ \qPtr -> do++   Call.assert "householderExtractQ: width shapes mismatch" (widthTau == width)++   let m = Shape.size height+   let k = min m $ Shape.size width+   let lda = m+   evalContT $ do+      mPtr <- Call.cint m+      kPtr <- Call.cint k+      qrPtr <- ContT $ withForeignPtr qr+      ldaPtr <- Call.cint lda+      tauPtr <- ContT $ withForeignPtr tau+      liftIO $+         case order of+            RowMajor -> do+               copySubMatrix k m k qrPtr lda qPtr+               withAutoWorkspaceInfo "unglq" $+                  LapackGen.unglq mPtr mPtr kPtr qPtr ldaPtr tauPtr+            ColumnMajor -> do+               copyBlock (m*k) qrPtr qPtr+               withAutoWorkspaceInfo "ungqr" $+                  LapackGen.ungqr mPtr mPtr kPtr qPtr ldaPtr tauPtr++householderExtractR ::+   (Shape.C height, Shape.C width, Eq width, Class.Floating a) =>+   Householder height width a -> General height width a+householderExtractR+      (Array (MatrixShape.Householder order height width) qr) =++   Array.unsafeCreate (MatrixShape.General order height width) $+      \rPtr -> do++   let (uplo, (m,n)) =+         case order of+            RowMajor -> ('L', (Shape.size width, Shape.size height))+            ColumnMajor -> ('U', (Shape.size height, Shape.size width))+   fill zero (m*n) rPtr+   evalContT $ do+      uploPtr <- Call.char uplo+      mPtr <- Call.cint m+      nPtr <- Call.cint n+      qrPtr <- ContT $ withForeignPtr qr+      ldqrPtr <- Call.cint m+      ldrPtr <- Call.cint m+      liftIO $ LapackGen.lacpy uploPtr mPtr nPtr qrPtr ldqrPtr rPtr ldrPtr++{- |+For an m-by-n-matrix @a@ with m>=n+this function computes an m-by-(m-n)-matrix @b@+such that @Matrix.multiply (transpose b) a@ is a zero matrix.+The function does not try to compensate a rank deficiency of @a@.+That is, @a|||b@ has full rank if and only if @a@ has full rank.++For full-rank matrices you might also call this @kernel@ or @nullspace@.+-}+orthogonalComplement ::+   (Shape.C height, Shape.C width, Eq width, Class.Floating a) =>+   General height width a -> General height ZeroInt a+orthogonalComplement a =+   dropColumns (Shape.size $ MatrixShape.generalWidth $ Array.shape a) $+   Array.mapShape zeroIntWidth $+   Square.toGeneral $ householderExtractQ $ householderDecompose a++zeroIntWidth ::+   (Shape.C width) =>+   MatrixShape.General height width -> MatrixShape.General height ZeroInt+zeroIntWidth (MatrixShape.General order height width) =+   MatrixShape.General order height (zeroInt $ Shape.size width)
src/Numeric/LAPACK/Private.hs view
@@ -1,23 +1,38 @@ {-# LANGUAGE TypeFamilies #-} module Numeric.LAPACK.Private where +import Numeric.LAPACK.Matrix.Shape.Private+         (Order(RowMajor, ColumnMajor), transposeFromOrder)+ 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.BLAS.FFI.Generic as BlasGen import qualified Numeric.Netlib.Utility as Call import qualified Numeric.Netlib.Class as Class  import Foreign.Marshal.Array (advancePtr)+import Foreign.Marshal.Alloc (alloca)+import Foreign.C.Types (CInt)+import Foreign.ForeignPtr (ForeignPtr, withForeignPtr, mallocForeignPtrArray) import Foreign.Ptr (Ptr)-import Foreign.Storable (Storable, peek)+import Foreign.Storable (Storable, poke, peek) -import Control.Monad.Trans.Cont (evalContT)+import Text.Printf (printf)++import Control.Monad.Trans.Cont (ContT(ContT), evalContT) import Control.Monad.IO.Class (liftIO)+import Control.Monad (foldM)+import Control.Applicative ((<$>))  import Data.Functor.Identity (Identity(Identity, runIdentity)) -import Data.Complex (Complex)+import qualified Data.Array.Comfort.Shape as Shape+import Data.Array.Comfort.Storable.Internal (Array(Array)) +import qualified Data.Complex as Complex+import Data.Complex (Complex((:+)))+ import Prelude hiding (sum)  @@ -25,10 +40,12 @@  type instance RealOf Float = Float type instance RealOf Double = Double-type instance RealOf (Complex Float) = Float-type instance RealOf (Complex Double) = Double+type instance RealOf (Complex a) = a  +type ComplexOf x = Complex (RealOf x)++ zero, one, minusOne :: Class.Floating a => a zero =    runIdentity $@@ -41,11 +58,8 @@    Class.switchFloating       (Identity (-1)) (Identity (-1)) (Identity (-1)) (Identity (-1)) -oneReal :: Class.Real a => a-oneReal = runIdentity $ Class.switchReal (Identity 1) (Identity 1)  - fill :: (Class.Floating a) => a -> Int -> Ptr a -> IO () fill a n dstPtr = evalContT $ do    nPtr <- Call.cint n@@ -62,12 +76,20 @@    incyPtr <- Call.cint 1    liftIO $ BlasGen.copy nPtr srcPtr incxPtr dstPtr incyPtr +copyToTemp :: (Class.Floating a) => Int -> ForeignPtr a -> ContT r IO (Ptr a)+copyToTemp n fptr = do+   ptr <- ContT $ withForeignPtr fptr+   tmpPtr <- Call.allocaArray n+   liftIO $ copyBlock n ptr tmpPtr+   return tmpPtr++ {- | In ColumnMajor: Copy a m-by-n-matrix with lda>=m and ldb>=m. -} copySubMatrix ::-   (Storable a, Class.Floating a) =>+   (Class.Floating a) =>    Int -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO () copySubMatrix m n lda aPtr ldb bPtr = evalContT $ do    uploPtr <- Call.char 'A'@@ -78,7 +100,7 @@    liftIO $ LapackGen.lacpy uploPtr mPtr nPtr aPtr ldaPtr bPtr ldbPtr  copyTransposed ::-   (Storable a, Class.Floating a) =>+   (Class.Floating a) =>    Int -> Int -> Ptr a -> Int -> Ptr a -> IO () copyTransposed n m aPtr ldb bPtr = evalContT $ do    nPtr <- Call.cint n@@ -90,10 +112,57 @@          (pointerSeq 1 aPtr)          (pointerSeq ldb bPtr) ++{- |+Copy a m-by-n-matrix to ColumnMajor order.+-}+copyToColumnMajor ::+   (Class.Floating a) =>+   Order -> Int -> Int -> Ptr a -> Ptr a -> IO ()+copyToColumnMajor order m n aPtr bPtr =+   case order of+      RowMajor -> copyTransposed m n aPtr m bPtr+      ColumnMajor -> copyBlock (m*n) aPtr bPtr++copyToSubColumnMajor ::+   (Class.Floating a) =>+   Order -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()+copyToSubColumnMajor order m n aPtr ldb bPtr =+   case order of+      RowMajor -> copyTransposed m n aPtr ldb bPtr+      ColumnMajor ->+         if m==ldb+           then copyBlock (m*n) aPtr bPtr+           else copySubMatrix m n m aPtr ldb bPtr++ pointerSeq :: (Storable a) => Int -> Ptr a -> [Ptr a] pointerSeq k ptr = iterate (flip advancePtr k) ptr  +allocArray :: (Shape.C sh, Storable a) => sh -> ContT r IO (Array sh a, Ptr a)+allocArray sh = do+   fptr <- liftIO $ mallocForeignPtrArray $ Shape.size sh+   ptr <- ContT $ withForeignPtr fptr+   return (Array sh fptr, ptr)+++allocHigherArray ::+   (Shape.C sh, Class.Floating a) =>+   sh -> Int -> Int -> Int -> ContT r IO (Array sh a, (Ptr a, Int))+allocHigherArray shapeX m n nrhs = do+   (x,xPtr) <- allocArray shapeX+   if m>n+      then do+         tmpPtr <- Call.allocaArray (m*nrhs)+         ContT $ \act -> do+            r <- act (x,(tmpPtr,m))+            copySubMatrix n nrhs m tmpPtr n xPtr+            return r+      else return (x,(xPtr,n))+++ newtype Sum a = Sum {runSum :: Int -> Ptr a -> Int -> IO a}  sum :: Class.Floating a => Int -> Ptr a -> Int -> IO a@@ -110,7 +179,7 @@    evalContT $ do       nPtr <- Call.cint n       incxPtr <- Call.cint incx-      yPtr <- Call.real oneReal+      yPtr <- Call.real one       incyPtr <- Call.cint 0       liftIO $ BlasReal.dot nPtr xPtr incxPtr yPtr incyPtr @@ -135,3 +204,110 @@             transPtr mPtr nPtr alphaPtr aPtr ldaPtr             xPtr incxPtr betaPtr yPtr incyPtr       liftIO $ peek yPtr+++product :: Class.Floating a => Int -> Ptr a -> Int -> IO a+product n xPtr incx =+   foldM (\x ptr -> do y <- peek ptr; return $! x*y) one $+   take n $ pointerSeq incx xPtr+++newtype LACGV a = LACGV {getLACGV :: Ptr CInt -> Ptr a -> Ptr CInt -> IO ()}++lacgv :: Class.Floating a => Ptr CInt -> Ptr a -> Ptr CInt -> IO ()+lacgv =+   getLACGV $+   Class.switchFloating+      (LACGV $ const $ const $ const $ return ())+      (LACGV $ const $ const $ const $ return ())+      (LACGV LapackComplex.lacgv)+      (LACGV LapackComplex.lacgv)+++multiplyMatrix ::+   (Class.Floating a) =>+   Order -> Order -> Int -> Int -> Int ->+   ForeignPtr a -> ForeignPtr a -> Ptr a -> IO ()+multiplyMatrix orderA orderB m k n a b cPtr = do+   let lda = case orderA of RowMajor -> k; ColumnMajor -> m+   let ldb = case orderB of RowMajor -> n; ColumnMajor -> k+   let ldc = m+   evalContT $ do+      transaPtr <- Call.char $ transposeFromOrder orderA+      transbPtr <- Call.char $ transposeFromOrder orderB+      mPtr <- Call.cint m+      nPtr <- Call.cint n+      kPtr <- Call.cint k+      alphaPtr <- Call.number one+      aPtr <- ContT $ withForeignPtr a+      ldaPtr <- Call.cint lda+      bPtr <- ContT $ withForeignPtr b+      ldbPtr <- Call.cint ldb+      betaPtr <- Call.number zero+      ldcPtr <- Call.cint ldc+      liftIO $+         BlasGen.gemm+            transaPtr transbPtr mPtr nPtr kPtr alphaPtr aPtr ldaPtr+            bPtr ldbPtr betaPtr cPtr ldcPtr++++withAutoWorkspaceInfo ::+   (Class.Floating a) =>+   String -> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()+withAutoWorkspaceInfo name computation = evalContT $ do+   infoPtr <- Call.alloca+   liftIO $ withAutoWorkspace $ \workPtr lworkPtr ->+      computation workPtr lworkPtr infoPtr+   info <- liftIO $ fromIntegral <$> peek infoPtr+   case compare info (0::Int) of+      EQ -> return ()+      LT -> error $ printf "%s: illegal value in %d-th argument" name (-info)+      GT -> error $ printf "%s: deficient rank %d" name info++withAutoWorkspace ::+   (Class.Floating a) =>+   (Ptr a -> Ptr CInt -> IO ()) -> IO ()+withAutoWorkspace computation = evalContT $ do+   lworkPtr <- Call.cint (-1)+   lwork <- liftIO $ alloca $ \workPtr -> do+      computation workPtr lworkPtr+      ceilingSize <$> peek workPtr+   workPtr <- Call.allocaArray lwork+   liftIO $ poke lworkPtr $ fromIntegral lwork+   liftIO $ computation workPtr lworkPtr+++newtype FromReal a = FromReal {getFromReal :: RealOf a -> a}++fromReal :: (Class.Floating a) => RealOf a -> a+fromReal =+   getFromReal $+   Class.switchFloating+      (FromReal id)+      (FromReal id)+      (FromReal (:+0))+      (FromReal (:+0))++newtype RealPart a = RealPart {getRealPart :: a -> RealOf a}++realPart :: (Class.Floating a) => a -> RealOf a+realPart =+   getRealPart $+   Class.switchFloating+      (RealPart id)+      (RealPart id)+      (RealPart Complex.realPart)+      (RealPart Complex.realPart)+++newtype FuncArg b a = FuncArg {runFuncArg :: a -> b}++ceilingSize :: (Class.Floating a) => a -> Int+ceilingSize =+   runFuncArg $+   Class.switchFloating+      (FuncArg ceiling)+      (FuncArg ceiling)+      (FuncArg $ ceiling . Complex.realPart)+      (FuncArg $ ceiling . Complex.realPart)
+ src/Numeric/LAPACK/Singular.hs view
@@ -0,0 +1,406 @@+{-# LANGUAGE TypeFamilies #-}+module Numeric.LAPACK.Singular (+   values,+   decompose,+   decomposeNarrow,+   decomposeSquat,+   determinantAbsolute,+   leastSquaresMinimumNormRCond,+   pseudoInverseRCond,+   RealOf,+   ) where++import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape+import qualified Numeric.LAPACK.Matrix.Square as Square+import qualified Numeric.LAPACK.Matrix as Matrix+import qualified Numeric.LAPACK.Vector as Vector+import Numeric.LAPACK.Matrix.Square (Square)+import Numeric.LAPACK.Matrix.Shape.Private (Order(RowMajor,ColumnMajor))+import Numeric.LAPACK.Matrix.Private (General, ZeroInt, zeroInt)+import Numeric.LAPACK.Vector (Vector)+import Numeric.LAPACK.Private+         (RealOf, withAutoWorkspace, fromReal, allocArray, allocHigherArray,+          copyToTemp, copyToColumnMajor, copyToSubColumnMajor)++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.Internal as Array+import qualified Data.Array.Comfort.Shape as Shape+import Data.Array.Comfort.Storable.Internal (Array(Array))++import System.IO.Unsafe (unsafePerformIO)++import Foreign.Marshal.Array (allocaArray)+import Foreign.Marshal.Alloc (alloca)+import Foreign.C.Types (CInt, CChar)+import Foreign.ForeignPtr (withForeignPtr)+import Foreign.Ptr (Ptr, nullPtr)+import Foreign.Storable (Storable, peek)++import Control.Monad.Trans.Cont (evalContT)+import Control.Monad.IO.Class (liftIO)+import Control.Applicative ((<$>))++import Text.Printf (printf)++import Data.Complex (Complex)+++values ::+   (Shape.C height, Shape.C width, Class.Floating a) =>+   General height width a -> Vector ZeroInt (RealOf a)+values =+   getValues $+   Class.switchFloating+      (Values valuesAux) (Values valuesAux)+      (Values valuesAux) (Values valuesAux)++type Values_ height width a =+   General height width a -> Vector ZeroInt (RealOf a)++newtype Values height width a = Values {getValues :: Values_ height width a}++valuesAux ::+   (Shape.C height, Shape.C width,+    Class.Floating a, RealOf a ~ ar, Storable ar) =>+   Values_ height width a+valuesAux (Array shape@(MatrixShape.General _order height width) a) =+   Array.unsafeCreateWithSize+      (zeroInt $ min (Shape.size height) (Shape.size width)) $ \mn sPtr -> do+   let (m,n) = MatrixShape.dimensions shape+   let lda = m+   evalContT $ do+      jobuPtr <- Call.char 'N'+      jobvtPtr <- Call.char 'N'+      mPtr <- Call.cint m+      nPtr <- Call.cint n+      aPtr <- copyToTemp (m*n) a+      ldaPtr <- Call.cint lda+      let uPtr = nullPtr+      let vtPtr = nullPtr+      lduPtr <- Call.cint m+      ldvtPtr <- Call.cint n+      liftIO $+         withInfo "gesvd" $ \infoPtr ->+         gesvd jobuPtr jobvtPtr mPtr nPtr+            aPtr ldaPtr sPtr uPtr lduPtr vtPtr ldvtPtr mn infoPtr+++determinantAbsolute ::+   (Shape.C height, Shape.C width, Class.Floating a) =>+   General height width a -> RealOf a+determinantAbsolute =+   getDeterminantAbsolute $+   Class.switchFloating+      (DeterminantAbsolute determinantAbsoluteAux)+      (DeterminantAbsolute determinantAbsoluteAux)+      (DeterminantAbsolute determinantAbsoluteAux)+      (DeterminantAbsolute determinantAbsoluteAux)++newtype DeterminantAbsolute f a =+   DeterminantAbsolute {+      getDeterminantAbsolute :: f a -> RealOf a+   }++determinantAbsoluteAux ::+   (Shape.C height, Shape.C width,+    Class.Floating a, RealOf a ~ ar, Class.Floating ar) =>+   General height width a -> ar+determinantAbsoluteAux = Vector.product . values+++decompose ::+   (Shape.C height, Shape.C width, Class.Floating a) =>+   General height width a ->+   (Square height a, Vector ZeroInt (RealOf a), Square width a)+decompose =+   getDecompose $+   Class.switchFloating+      (Decompose decomposeAux)+      (Decompose decomposeAux)+      (Decompose decomposeAux)+      (Decompose decomposeAux)++newtype Decompose m f v g a =+   Decompose {+      getDecompose :: m a -> (f a, v (RealOf a), g a)+   }++decomposeAux ::+   (Shape.C height, Shape.C width,+    Class.Floating a, RealOf a ~ ar, Storable ar) =>+   General height width a ->+   (Square height a, Vector ZeroInt ar, Square width a)+decomposeAux (Array (MatrixShape.General order height width) a) =+   unsafePerformIO $ evalContT $ do+      (u,uPtr0) <- allocArray (MatrixShape.Square order height)+      (vt,vtPtr0) <- allocArray (MatrixShape.Square order width)+      let ((m,n),(uPtr,vtPtr)) =+            case order of+               RowMajor ->+                  ((Shape.size width, Shape.size height), (vtPtr0,uPtr0))+               ColumnMajor ->+                  ((Shape.size height, Shape.size width), (uPtr0,vtPtr0))+      let mn = min m n+      let lda = m+      jobuPtr <- Call.char 'A'+      jobvtPtr <- Call.char 'A'+      mPtr <- Call.cint m+      nPtr <- Call.cint n+      aPtr <- copyToTemp (m*n) a+      ldaPtr <- Call.cint lda+      (s,sPtr) <- allocArray (zeroInt mn)+      lduPtr <- Call.cint m+      ldvtPtr <- Call.cint n+      liftIO $+         withInfo "gesvd" $ \infoPtr ->+         gesvd jobuPtr jobvtPtr mPtr nPtr+            aPtr ldaPtr sPtr uPtr lduPtr vtPtr ldvtPtr mn infoPtr+      return (u, s, vt)+++decomposeSquat ::+   (Shape.C height, Shape.C width, Class.Floating a) =>+   General height width a ->+   (Square height a, Vector height (RealOf a), General height width a)+decomposeSquat a =+   let (u,s,vt) = decomposeNarrow $ Matrix.transpose a+   in  (Square.transpose vt, s, Matrix.transpose u)++decomposeNarrow ::+   (Shape.C height, Shape.C width, Class.Floating a) =>+   General height width a ->+   (General height width a, Vector width (RealOf a), Square width a)+decomposeNarrow =+   getDecompose $+   Class.switchFloating+      (Decompose decomposeThin)+      (Decompose decomposeThin)+      (Decompose decomposeThin)+      (Decompose decomposeThin)++decomposeThin ::+   (Shape.C height, Shape.C width,+    Class.Floating a, RealOf a ~ ar, Storable ar) =>+   General height width a ->+   (General height width a, Vector width ar, Square width a)+decomposeThin (Array (MatrixShape.General order height width) a) =+   unsafePerformIO $ do+     Call.assert "Singular.decomposeThin: matrix is wider than high"+         (Shape.size height >= Shape.size width)+     evalContT $ do+      (u,uPtr0) <- allocArray (MatrixShape.General order height width)+      (vt,vtPtr0) <- allocArray (MatrixShape.Square order width)+      let ((m,n),(uPtr,vtPtr)) =+            case order of+               RowMajor ->+                  ((Shape.size width, Shape.size height), (vtPtr0,uPtr0))+               ColumnMajor ->+                  ((Shape.size height, Shape.size width), (uPtr0,vtPtr0))+      let mn = min m n+      let lda = m+      jobuPtr <- Call.char 'S'+      jobvtPtr <- Call.char 'S'+      mPtr <- Call.cint m+      nPtr <- Call.cint n+      aPtr <- copyToTemp (m*n) a+      ldaPtr <- Call.cint lda+      (s,sPtr) <- allocArray width+      lduPtr <- Call.cint m+      ldvtPtr <- Call.cint mn+      liftIO $+         withInfo "gesvd" $ \infoPtr ->+         gesvd jobuPtr jobvtPtr mPtr nPtr+            aPtr ldaPtr sPtr uPtr lduPtr vtPtr ldvtPtr mn infoPtr+      return (u, s, vt)+++type GESVD_ ar a =+   Ptr CChar -> Ptr CChar -> Ptr CInt -> Ptr CInt ->+   Ptr a -> Ptr CInt -> Ptr ar ->+   Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> Int -> Ptr CInt -> IO ()++newtype GESVD a = GESVD {getGESVD :: GESVD_ (RealOf a) a}++gesvd :: Class.Floating a => GESVD_ (RealOf a) a+gesvd =+   getGESVD $+   Class.switchFloating+      (GESVD gesvdReal)+      (GESVD gesvdReal)+      (GESVD gesvdComplex)+      (GESVD gesvdComplex)++gesvdReal :: (Class.Real a) => GESVD_ a a+gesvdReal jobuPtr jobvtPtr mPtr nPtr+      aPtr ldaPtr sPtr uPtr lduPtr vtPtr ldvtPtr _mn infoPtr =+   withAutoWorkspace $ \workPtr lworkPtr ->+   LapackReal.gesvd jobuPtr jobvtPtr+      mPtr nPtr aPtr ldaPtr sPtr uPtr lduPtr vtPtr ldvtPtr+      workPtr lworkPtr infoPtr++gesvdComplex :: (Class.Real a) => GESVD_ a (Complex a)+gesvdComplex jobuPtr jobvtPtr+      mPtr nPtr aPtr ldaPtr sPtr uPtr lduPtr vtPtr ldvtPtr mn infoPtr =+   allocaArray (5*mn) $ \rworkPtr ->+   withAutoWorkspace $ \workPtr lworkPtr ->+   LapackComplex.gesvd jobuPtr jobvtPtr+      mPtr nPtr aPtr ldaPtr sPtr uPtr lduPtr vtPtr ldvtPtr+      workPtr lworkPtr rworkPtr infoPtr+++leastSquaresMinimumNormRCond ::+   (Shape.C height, Eq height, Shape.C width, Shape.C nrhs, Class.Floating a) =>+   RealOf a -> General height width a -> General height nrhs a ->+   (Int, General width nrhs a)+leastSquaresMinimumNormRCond =+   getLeastSquaresMinimumNormRCond $+   Class.switchFloating+      (LeastSquaresMinimumNormRCond leastSquaresMinimumNormRCondAux)+      (LeastSquaresMinimumNormRCond leastSquaresMinimumNormRCondAux)+      (LeastSquaresMinimumNormRCond leastSquaresMinimumNormRCondAux)+      (LeastSquaresMinimumNormRCond leastSquaresMinimumNormRCondAux)++newtype LeastSquaresMinimumNormRCond f g h a =+   LeastSquaresMinimumNormRCond {+      getLeastSquaresMinimumNormRCond ::+         RealOf a -> f a -> g a -> (Int, h a)+   }++leastSquaresMinimumNormRCondAux ::+   (Shape.C height, Eq height, Shape.C width, Shape.C nrhs,+    Class.Floating a, RealOf a ~ ar, Class.Floating ar) =>+   ar -> General height width a -> General height nrhs a ->+   (Int, General width nrhs a)+leastSquaresMinimumNormRCondAux rcond+   (Array (MatrixShape.General orderA heightA widthA) a)+   (Array (MatrixShape.General orderB heightB widthB) b) =+      unsafePerformIO $ do+   Call.assert "leastSquaresMinimumNorm: height shapes mismatch"+      (heightA == heightB)+   let shapeX = MatrixShape.General ColumnMajor widthA widthB+   let m = Shape.size heightA+   let n = Shape.size widthA+   let nrhs = Shape.size widthB+   let mn = min m n+   let aSize = m*n+   let lda = m+   evalContT $ do+      mPtr <- Call.cint m+      nPtr <- Call.cint n+      nrhsPtr <- Call.cint nrhs+      aPtr <- Call.allocaArray aSize+      liftIO $ withForeignPtr a $ \asrcPtr ->+         copyToColumnMajor orderA m n asrcPtr aPtr+      ldaPtr <- Call.cint lda+      (x,(tmpPtr,ldtmp)) <- allocHigherArray shapeX m n nrhs+      ldtmpPtr <- Call.cint ldtmp+      liftIO $ withForeignPtr b $ \bPtr ->+         copyToSubColumnMajor orderB m nrhs bPtr ldtmp tmpPtr++      sPtr <- Call.allocaArray mn+      rcondPtr <- Call.number rcond+      rankPtr <- Call.alloca+      liftIO $+         withInfo "gelss" $ \infoPtr ->+         gelss mPtr nPtr nrhsPtr aPtr ldaPtr tmpPtr ldtmpPtr sPtr rcondPtr+            rankPtr mn infoPtr++      rank <- liftIO $ fromIntegral <$> peek rankPtr+      return (rank, x)+++type GELSS_ ar a =+   Ptr CInt -> Ptr CInt -> Ptr CInt ->+   Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt ->+   Ptr ar -> Ptr ar -> Ptr CInt -> Int -> Ptr CInt -> IO ()++newtype GELSS a = GELSS {getGELSS :: GELSS_ (RealOf a) a}++gelss :: Class.Floating a => GELSS_ (RealOf a) a+gelss =+   getGELSS $+   Class.switchFloating+      (GELSS gelssReal)+      (GELSS gelssReal)+      (GELSS gelssComplex)+      (GELSS gelssComplex)++gelssReal :: (Class.Real a) => GELSS_ a a+gelssReal mPtr nPtr nrhsPtr aPtr ldaPtr bPtr ldbPtr sPtr rcondPtr+      rankPtr _mn infoPtr =+   withAutoWorkspace $ \workPtr lworkPtr ->+   LapackReal.gelss+      mPtr nPtr nrhsPtr aPtr ldaPtr bPtr ldbPtr sPtr rcondPtr+      rankPtr workPtr lworkPtr infoPtr++gelssComplex :: (Class.Real a) => GELSS_ a (Complex a)+gelssComplex mPtr nPtr nrhsPtr aPtr ldaPtr bPtr ldbPtr sPtr rcondPtr+      rankPtr mn infoPtr =+   allocaArray (5*mn) $ \rworkPtr ->+   withAutoWorkspace $ \workPtr lworkPtr ->+   LapackComplex.gelss+      mPtr nPtr nrhsPtr aPtr ldaPtr bPtr ldbPtr sPtr rcondPtr+      rankPtr workPtr lworkPtr rworkPtr infoPtr+++pseudoInverseRCond ::+   (Shape.C height, Eq height, Shape.C width, Eq width, Class.Floating a) =>+   RealOf a -> General height width a -> (Int, General width height a)+pseudoInverseRCond =+   getPseudoInverseRCond $+   Class.switchFloating+      (PseudoInverseRCond pseudoInverseRCondAux)+      (PseudoInverseRCond pseudoInverseRCondAux)+      (PseudoInverseRCond pseudoInverseRCondAux)+      (PseudoInverseRCond pseudoInverseRCondAux)++newtype PseudoInverseRCond f g a =+   PseudoInverseRCond {+      getPseudoInverseRCond :: RealOf a -> f a -> (Int, g a)+   }++pseudoInverseRCondAux ::+   (Shape.C height, Eq height, Shape.C width, Eq width,+    Class.Floating a, RealOf a ~ ar, Class.Real ar) =>+   RealOf a -> General height width a -> (Int, General width height a)+pseudoInverseRCondAux rcond a =+   let (MatrixShape.General _ height width) = Array.shape a+   in if Shape.size height < Shape.size width+        then+          let (u,s,vt) = decomposeSquat a+              (rank,recipS) = recipSigma rcond s+          in  (rank,+               Matrix.multiply (Matrix.adjoint vt) $+               Matrix.scaleRows recipS $+               Square.toGeneral $ Square.adjoint u)+        else+          let (u,s,vt) = decomposeNarrow a+              (rank,recipS) = recipSigma rcond s+          in  (rank,+               Matrix.multiply (Square.toGeneral $ Square.adjoint vt) $+               Matrix.scaleRows recipS $ Matrix.adjoint u)++recipSigma ::+   (Shape.C sh, Class.Floating a, RealOf a ~ ar, Class.Real ar) =>+   ar -> Array sh ar -> (Int, Array sh a)+recipSigma rcond sigmas =+   case Array.toList sigmas of+      [] -> (0, Array.map fromReal sigmas)+      xs@(x:_) ->+         let smin = x * rcond+         in (length (takeWhile (>=smin) xs),+             Array.map (\s -> if s>=smin then fromReal (recip s) else 0) sigmas)+++withInfo :: String -> (Ptr CInt -> IO ()) -> IO ()+withInfo name computation = alloca $ \infoPtr -> do+   computation infoPtr+   info <- fromIntegral <$> peek infoPtr+   case compare info (0::Int) of+      EQ -> return ()+      LT -> error $ printf "%s: illegal value in %d-th argument" name (-info)+      GT -> error $ printf "%s: %d superdiagonals did not converge" name info
src/Numeric/LAPACK/Vector.hs view
@@ -1,6 +1,7 @@ module Numeric.LAPACK.Vector (    Vector,    fromList,+   autoFromList,    constant,    dot,    sum,@@ -8,6 +9,8 @@    norm1,    norm2,    argAbsMaximum,+   argAbs1Maximum,+   product,    scale,    add, sub ,    mac,@@ -18,6 +21,7 @@    ) where  import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape+import qualified Numeric.LAPACK.Matrix.Private as Matrix import qualified Numeric.LAPACK.Private as Private import Numeric.LAPACK.Private (RealOf, zero, one, minusOne, fill) @@ -48,7 +52,7 @@ import Data.Word (Word64) import Data.Bits (shiftR, (.&.)) -import Prelude hiding (sum)+import Prelude hiding (sum, product)   type Vector = Array@@ -57,11 +61,14 @@ fromList :: (Shape.C sh, Storable a) => sh -> [a] -> Vector sh a fromList = Array.fromList +autoFromList :: (Storable a) => [a] -> Vector (Shape.ZeroBased Int) a+autoFromList = Array.vectorFromList -constant :: (Shape.C sh, Storable a, Class.Floating a) => sh -> a -> Vector sh a-constant sh a = Array.unsafeCreate sh $ fill a (Shape.size sh) +constant :: (Shape.C sh, Class.Floating a) => sh -> a -> Vector sh a+constant sh a = Array.unsafeCreateWithSize sh $ fill a + newtype Dot sh a = Dot {runDot :: Vector sh a -> Vector sh a -> a}  dot ::@@ -119,14 +126,24 @@ sum (Array sh x) = unsafePerformIO $    withForeignPtr x $ \xPtr -> Private.sum (Shape.size sh) xPtr 1 -norm1 :: (Shape.C sh, Class.Real a) => Vector sh a -> a+norm1 :: (Shape.C sh, Class.Floating a) => Vector sh a -> RealOf a norm1 (Array sh x) = unsafePerformIO $    evalContT $ do       nPtr <- Call.cint $ Shape.size sh       sxPtr <- ContT $ withForeignPtr x       incxPtr <- Call.cint 1-      liftIO $ BlasReal.asum nPtr sxPtr incxPtr+      liftIO $ csum1 nPtr sxPtr incxPtr +csum1 :: Class.Floating a => Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a)+csum1 =+   getNorm $+   Class.switchFloating+      (Norm BlasReal.asum)+      (Norm BlasReal.asum)+      (Norm LapackComplex.sum1)+      (Norm LapackComplex.sum1)++ {- | Sum of the absolute values of real numbers or components of complex numbers. For real numbers it is equivalent to 'norm1'.@@ -170,12 +187,10 @@  {- | Returns the index and value of the element with the maximal absolute value.-The function does not strictly compare the absolute value of a complex number-but the sum of the absolute complex components. Caution: It actually returns the value of the element, not its absolute value! -} argAbsMaximum ::-   (Shape.C sh, Storable a, Class.Floating a) =>+   (Shape.C sh, Class.Floating a) =>    Vector sh a -> (Shape.Index sh, a) argAbsMaximum (Array sh x) = unsafePerformIO $    evalContT $ do@@ -183,28 +198,64 @@       sxPtr <- ContT $ withForeignPtr x       incxPtr <- Call.cint 1       liftIO $ do+         k <- fromIntegral . subtract 1 <$> absMax nPtr sxPtr incxPtr+         xmax <- peekElemOff sxPtr k+         return (Shape.indices sh !! k, xmax)++newtype ArgMaximum a =+   ArgMaximum {runArgMaximum :: Ptr CInt -> Ptr a -> Ptr CInt -> IO CInt}++absMax :: Class.Floating a => Ptr CInt -> Ptr a -> Ptr CInt -> IO CInt+absMax =+   runArgMaximum $+   Class.switchFloating+      (ArgMaximum BlasGen.iamax)+      (ArgMaximum BlasGen.iamax)+      (ArgMaximum LapackComplex.imax1)+      (ArgMaximum LapackComplex.imax1)+++{- |+Returns the index and value of the element with the maximal absolute value.+The function does not strictly compare the absolute value of a complex number+but the sum of the absolute complex components.+Caution: It actually returns the value of the element, not its absolute value!+-}+argAbs1Maximum ::+   (Shape.C sh, Class.Floating a) =>+   Vector sh a -> (Shape.Index sh, a)+argAbs1Maximum (Array sh x) = unsafePerformIO $+   evalContT $ do+      nPtr <- Call.cint $ Shape.size sh+      sxPtr <- ContT $ withForeignPtr x+      incxPtr <- Call.cint 1+      liftIO $ do          k <- fromIntegral . subtract 1 <$> BlasGen.iamax nPtr sxPtr incxPtr          xmax <- peekElemOff sxPtr k          return (Shape.indices sh !! k, xmax)  +product :: (Shape.C sh, Class.Floating a) => Vector sh a -> a+product (Array sh x) = unsafePerformIO $+   withForeignPtr x $ \xPtr -> Private.product (Shape.size sh) xPtr 1++ scale, _scale ::-   (Shape.C sh, Storable a, Class.Floating a) =>+   (Shape.C sh, Class.Floating a) =>    a -> Vector sh a -> Vector sh a-scale alpha (Array sh x) = Array.unsafeCreate sh $ \syPtr -> do+scale alpha (Array sh x) = Array.unsafeCreateWithSize sh $ \n syPtr -> do    evalContT $ do       alphaPtr <- Call.number alpha-      nPtr <- Call.cint $ Shape.size sh+      nPtr <- Call.cint n       sxPtr <- ContT $ withForeignPtr x       incxPtr <- Call.cint 1       incyPtr <- Call.cint 1       liftIO $ BlasGen.copy nPtr sxPtr incxPtr syPtr incyPtr       liftIO $ BlasGen.scal nPtr alphaPtr syPtr incyPtr -_scale a (Array sh b) = Array.unsafeCreate sh $ \cPtr -> do+_scale a (Array sh b) = Array.unsafeCreateWithSize sh $ \n cPtr -> do    let m = 1    let k = 1-   let n = Shape.size sh    evalContT $ do       transaPtr <- Call.char 'N'       transbPtr <- Call.char 'N'@@ -224,18 +275,19 @@             aPtr ldaPtr bPtr ldbPtr betaPtr cPtr ldcPtr  add, sub ::-   (Shape.C sh, Eq sh, Storable a, Class.Floating a) =>+   (Shape.C sh, Eq sh, Class.Floating a) =>    Vector sh a -> Vector sh a -> Vector sh a add = mac one sub x y = mac minusOne y x  mac ::-   (Shape.C sh, Eq sh, Storable a, Class.Floating a) =>+   (Shape.C sh, Eq sh, Class.Floating a) =>    a -> Vector sh a -> Vector sh a -> Vector sh a-mac alpha (Array shX x) (Array shY y) = Array.unsafeCreate shX $ \szPtr -> do+mac alpha (Array shX x) (Array shY y) =+      Array.unsafeCreateWithSize shX $ \n szPtr -> do    Call.assert "mac: shapes mismatch" (shX == shY)    evalContT $ do-      nPtr <- Call.cint $ Shape.size shX+      nPtr <- Call.cint n       saPtr <- Call.number alpha       sxPtr <- ContT $ withForeignPtr x       incxPtr <- Call.cint 1@@ -246,11 +298,11 @@       liftIO $ BlasGen.axpy nPtr saPtr sxPtr incxPtr szPtr inczPtr  mul ::-   (Shape.C sh, Eq sh, Storable a, Class.Floating a) =>+   (Shape.C sh, Eq sh, Class.Floating a) =>    Vector sh a -> Vector sh a -> Vector sh a-mul (Array shA a) (Array shX x) = Array.unsafeCreate shX $ \yPtr -> do+mul (Array shA a) (Array shX x) =+      Array.unsafeCreateWithSize shX $ \n yPtr -> do    Call.assert "mul: shapes mismatch" (shA == shX)-   let n = Shape.size shX    evalContT $ do       transPtr <- Call.char 'N'       mPtr <- Call.cint n@@ -271,9 +323,8 @@   outer ::-   (Shape.C shx, Eq shx, Shape.C shy, Eq shy,-    Storable a, Class.Floating a) =>-   Vector shx a -> Vector shy a -> Array (MatrixShape.General shx shy) a+   (Shape.C shx, Eq shx, Shape.C shy, Eq shy, Class.Floating a) =>+   Vector shx a -> Vector shy a -> Matrix.General shx shy a outer (Array shX x) (Array shY y) =    Array.unsafeCreate (MatrixShape.General MatrixShape.ColumnMajor shX shY) $       \cPtr -> do@@ -301,7 +352,7 @@ newtype Conjugate sh a = Conjugate {getConjugate :: Vector sh a -> Vector sh a}  conjugate ::-   (Shape.C sh, Storable a, Class.Floating a) =>+   (Shape.C sh, Class.Floating a) =>    Vector sh a -> Vector sh a conjugate =    getConjugate $@@ -312,11 +363,11 @@       (Conjugate complexConjugate)  complexConjugate ::-   (Shape.C sh, Storable a, Class.Real a) =>+   (Shape.C sh, Class.Real a) =>    Vector sh (Complex a) -> Vector sh (Complex a)-complexConjugate (Array sh x) = Array.unsafeCreate sh $ \syPtr ->+complexConjugate (Array sh x) = Array.unsafeCreateWithSize sh $ \n syPtr ->    evalContT $ do-      nPtr <- Call.cint $ Shape.size sh+      nPtr <- Call.cint n       sxPtr <- ContT $ withForeignPtr x       incxPtr <- Call.cint 1       incyPtr <- Call.cint 1@@ -339,11 +390,11 @@ Only the least significant 47 bits of @seed@ are used. -} random ::-   (Shape.C sh, Storable a, Class.Floating a) =>+   (Shape.C sh, Class.Floating a) =>    RandomDistribution -> sh -> Word64 -> Vector sh a-random dist sh seed = Array.unsafeCreate sh $ \xPtr ->+random dist sh seed = Array.unsafeCreateWithSize sh $ \n xPtr ->    evalContT $ do-      nPtr <- Call.cint $ Shape.size sh+      nPtr <- Call.cint n       distPtr <-          Call.cint $          case (getConst $ isComplexInFunctor xPtr, dist) of