packages feed

lapack (empty) → 0.0

raw patch · 10 files changed

+1933/−0 lines, 10 filesdep +basedep +blas-ffidep +comfort-arraysetup-changed

Dependencies added: base, blas-ffi, comfort-array, lapack-ffi, netlib-ffi, non-empty, transformers, utility-ht

Files

+ LICENSE view
@@ -0,0 +1,27 @@+Copyright (c) Henning Thielemann 2018++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions+are met:+1. Redistributions of source code must retain the above copyright+   notice, this list of conditions and the following disclaimer.+2. Redistributions in binary form must reproduce the above copyright+   notice, this list of conditions and the following disclaimer in the+   documentation and/or other materials provided with the distribution.+3. Neither the name of the author nor the names of his contributors+   may be used to endorse or promote products derived from this software+   without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE+ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS+OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)+HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT+LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY+OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF+SUCH DAMAGE.
+ Setup.lhs view
@@ -0,0 +1,3 @@+#! /usr/bin/env runhaskell+> import Distribution.Simple+> main = defaultMain
+ lapack.cabal view
@@ -0,0 +1,66 @@+Name:             lapack+Version:          0.0+License:          BSD3+License-File:     LICENSE+Author:           Henning Thielemann <haskell@henning-thielemann.de>+Maintainer:       Henning Thielemann <haskell@henning-thielemann.de>+Homepage:         http://hub.darcs.net/thielema/lapack/+Category:         Data Structures+Synopsis:         Numerical Linear Algebra using LAPACK+Description:+  This is a high-level interface to LAPACK.+  .+  Features:+  .+  * Based on @comfort-array@:+    Allows to precisely express one-column or one-row matrices,+    as well as dense, square, triangular, banded and symmetric matrices.+  .+  * Support all data types that are supported by LAPACK,+    i.e. Float, Double, Complex Float, Complex Double+  .+  * No need for c2hs, hsc, Template Haskell or C helper functions+  .+  * Dependency only on BLAS and LAPACK, no GSL+  .+  * Separate formatting operator @(##)@:+    Works better for tuples of matrices and vectors than 'show'.+    'Show' is used for code one-liners+    that can be copied back into Haskell modules.+  .+  See also: @hmatrix@.+Tested-With:      GHC==7.4.2, GHC==7.8.4, GHC==8.2.2+Cabal-Version:    >=1.6+Build-Type:       Simple++Source-Repository this+  Tag:         0.0+  Type:        darcs+  Location:    http://hub.darcs.net/thielema/lapack/++Source-Repository head+  Type:        darcs+  Location:    http://hub.darcs.net/thielema/lapack/++Library+  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,+    transformers >=0.3 && <0.6,+    non-empty >=0.3 && <0.4,+    utility-ht >=0.0.10 && <0.1,+    base >=4.5 && <5++  GHC-Options:      -Wall+  Hs-Source-Dirs:   src+  Exposed-Modules:+    Numeric.LAPACK.Matrix+    Numeric.LAPACK.Matrix.Shape+    Numeric.LAPACK.Vector+    Numeric.LAPACK.LinearSystem+  Other-Modules:+    Numeric.LAPACK.Matrix.Shape.Private+    Numeric.LAPACK.Private+    Numeric.LAPACK.Format
+ src/Numeric/LAPACK/Format.hs view
@@ -0,0 +1,125 @@+module Numeric.LAPACK.Format where++import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape+import Numeric.LAPACK.Matrix.Shape.Private (Order(RowMajor, ColumnMajor))++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 (Array)++import Foreign.Storable (Storable)++import Text.Printf (PrintfArg, printf)++import qualified Data.List.HT as ListHT+import Data.Complex (Complex((:+)))+++infix 0 ##++(##) :: (Format a) => a -> String -> IO ()+a ## fmt = putStr $ unlines $ format fmt a+++class Format a where+   format :: String -> a -> [String]++instance Format Int where+   format _fmt a = [show a]++instance Format Float where+   format fmt a = [printf fmt a]++instance Format Double where+   format fmt a = [printf fmt a]++instance (PrintfArg a) => Format (Complex a) where+   format fmt a = [printfComplex fmt a]++instance (Format a, Format b) => Format (a,b) where+   format fmt (a,b) = format fmt a ++ [""] ++ format fmt b++instance (Format a, Format b, Format c) => Format (a,b,c) where+   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+   format = formatArray+++class (Shape.C sh) => FormatArray sh where+   formatArray ::+      (Storable a, 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]++instance (Integral i) => FormatArray (Shape.OneBased i) where+   formatArray fmt m = [unwords $ map (printfFloating fmt) $ Array.toList m]++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]+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++instance+   (Shape.C height, Shape.C width) =>+      FormatArray (MatrixShape.Householder height width) where+   formatArray = formatHouseholder++formatHouseholder ::+   (Shape.C height, Shape.C width, Storable a, 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  zipWith+         (\row xs ->+            concat $+            zipWith (\col cell -> (if row==col then '|' else ' '):cell) [0..] $+            zipWith (ListHT.padLeft ' ') strWidths xs)+         [(0::Int)..] xss++formatRows ::+   (Class.Floating a, Shape.C height, Shape.C width) =>+   String -> Order -> (height, width) -> [a] -> [[String]]+formatRows fmt order (height,width) =+   (case order of+      RowMajor -> ListHT.sliceVertical (Shape.size width)+      ColumnMajor -> ListHT.sliceHorizontal (Shape.size height)) .+   map (printfFloating fmt)++columnWidths :: [[[a]]] -> [Int]+columnWidths xss =+   case map (map length) xss of+      [] -> []+      w:ws -> foldl (zipWith max) w ws+++newtype Printf a = Printf {runPrintf :: String -> a -> String}++printfFloating :: (Class.Floating a) => String -> a -> String+printfFloating =+   runPrintf $+   Class.switchFloating+      (Printf printf)+      (Printf printf)+      (Printf printfComplex)+      (Printf printfComplex)++printfComplex :: (PrintfArg a) => String -> Complex a -> String+printfComplex fmt (r:+i) = printf (fmt ++ "+i" ++ fmt) r i
+ src/Numeric/LAPACK/LinearSystem.hs view
@@ -0,0 +1,478 @@+{-# 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
@@ -0,0 +1,596 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}+module Numeric.LAPACK.Matrix (+   General,+   (Format.##),+   Format.Format,+   Format.FormatArray,+   ZeroInt, zeroInt,+   transpose,+   fromList,+   identity,+   diagonal, getDiagonal,+   fromRows,    fromRowsWithSize,+   fromColumns, fromColumnsWithSize,+   singleRow,   singleColumn,+   flattenRow,  flattenColumn,+   pickRow, pickColumn,+   takeRows, takeColumns,+   dropRows, dropColumns,+   reverseRows, reverseColumns,+   fromRowMajor, toRowMajor, flatten,+   (|||),+   (===),++   rowSums, columnSums,+   scaleRows, scaleColumns,+   multiply,+   multiplyVector,++   trace,+   ) 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.Vector as Vector+import Numeric.LAPACK.Matrix.Shape.Private+         (Order(RowMajor, ColumnMajor), charFromOrder)+import Numeric.LAPACK.Private+         (zero, one, pointerSeq, copyTransposed, copySubMatrix, copyBlock)++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 Data.Array.Comfort.Shape ((:+:)((:+:)))++import Foreign.Marshal.Array (copyArray, advancePtr, pokeArray)+import Foreign.ForeignPtr (ForeignPtr, withForeignPtr)+import Foreign.Ptr (Ptr)+import Foreign.Storable (Storable, poke, peek)++import System.IO.Unsafe (unsafePerformIO)++import Control.Monad.Trans.Cont (ContT(ContT), evalContT)+import Control.Monad.IO.Class (liftIO)++import qualified Data.NonEmpty as NonEmpty+import Data.Foldable (forM_)+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+++type ZeroInt = Shape.ZeroBased Int++zeroInt :: Int -> ZeroInt+zeroInt = Shape.ZeroBased+++fromList ::+   (Shape.C height, Shape.C width, Storable a) =>+   height -> width -> [a] -> General height width a+fromList height width =+   Array.fromList (MatrixShape.General RowMajor height width)+++type Vector = Array+++identity, _identity ::+   (Shape.C sh, Storable a, 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++diagonal ::+   (Shape.C sh, Storable a, 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++getDiagonal ::+   (Shape.C sh, Eq sh, Storable a, 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+++singleRow :: Vector width a -> General () width a+singleRow (Array sh fptr) =+   Array (MatrixShape.General RowMajor () sh) fptr++singleColumn :: Vector width a -> General width () a+singleColumn (Array sh fptr) =+   Array (MatrixShape.General ColumnMajor sh ()) fptr++flattenRow :: General () width a -> Vector width a+flattenRow (Array (MatrixShape.General _ () sh) fptr) = Array sh fptr++flattenColumn :: General width () a -> Vector width a+flattenColumn (Array (MatrixShape.General _ sh ()) fptr) = Array sh fptr+++fromRows ::+   (Shape.C width, Eq width, Storable a) =>+   NonEmpty.T [] (Vector width a) -> General ZeroInt width a+fromRows (NonEmpty.Cons row rows) =+   fromRowsWithSize (Array.shape row) (row:rows)++fromRowsWithSize ::+   (Shape.C width, Eq width, Storable a) =>+   width -> [Vector width a] -> General ZeroInt width a+fromRowsWithSize width rows =+   Array.unsafeCreate+      (MatrixShape.General RowMajor (zeroInt $ length rows) width)+      (gather width rows)++fromColumns ::+   (Shape.C height, Eq height, Storable a) =>+   NonEmpty.T [] (Vector height a) -> General height ZeroInt a+fromColumns (NonEmpty.Cons column columns) =+   fromColumnsWithSize (Array.shape column) (column:columns)++fromColumnsWithSize ::+   (Shape.C height, Eq height, Storable a) =>+   height -> [Vector height a] -> General height ZeroInt a+fromColumnsWithSize height columns =+   Array.unsafeCreate+      (MatrixShape.General ColumnMajor height (zeroInt $ length columns))+      (gather height columns)++gather ::+   (Shape.C width, Eq width, Storable a) =>+   width -> [Array width a] -> Ptr a -> IO ()+gather width rows dstPtr =+   let widthSize = Shape.size width+   in forM_ (zip (pointerSeq widthSize dstPtr) rows) $+         \(dstRowPtr, Array.Array rowWidth srcFPtr) ->+         withForeignPtr srcFPtr $ \srcPtr -> do+            Call.assert+               "Matrix.fromRows/fromColumns: non-matching vector size"+               (width == rowWidth)+            copyArray dstRowPtr srcPtr widthSize+++pickRow ::+   (Shape.C height, Shape.C width, Shape.Index height ~ ix,+    Storable a, Class.Floating a) =>+   General height width a -> ix -> Vector width a+pickRow (Array (MatrixShape.General order height width) x) ix =+   case order of+      RowMajor -> pickConsecutive height width x ix+      ColumnMajor -> pickScattered width height x ix++pickColumn ::+   (Shape.C height, Shape.C width, Shape.Index width ~ ix,+    Storable a, Class.Floating a) =>+   General height width a -> ix -> Vector height a+pickColumn (Array (MatrixShape.General order height width) x) ix =+   case order of+      RowMajor -> pickScattered height width x ix+      ColumnMajor -> pickConsecutive width height x ix++pickConsecutive ::+   (Shape.C height, Shape.C width, Shape.Index height ~ ix,+    Storable 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+      let offset = Shape.offset height ix+      nPtr <- Call.cint n+      xPtr <- ContT $ withForeignPtr x+      incxPtr <- Call.cint 1+      incyPtr <- Call.cint 1+      liftIO $+         BlasGen.copy nPtr (advancePtr xPtr (n*offset)) incxPtr yPtr incyPtr++pickScattered ::+   (Shape.C height, Shape.C width, Shape.Index width ~ ix,+    Storable 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+      let offset = Shape.offset width ix+      nPtr <- Call.cint n+      xPtr <- ContT $ withForeignPtr x+      incxPtr <- Call.cint $ Shape.size width+      incyPtr <- Call.cint 1+      liftIO $+         BlasGen.copy nPtr (advancePtr xPtr offset) incxPtr yPtr incyPtr+++takeRows, dropRows ::+   (Shape.C width, Storable a, 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 ->+      withForeignPtr a $ \aPtr ->+      case order of+         RowMajor -> copyBlock (heightB*n) aPtr bPtr+         ColumnMajor -> copySubMatrix heightB n heightA aPtr heightB bPtr++dropRows k0+      (Array (MatrixShape.General order (Shape.ZeroBased heightA) width) a) =+   let k = min k0 heightA+       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 ->+      withForeignPtr a $ \aPtr ->+      case order of+         RowMajor -> copyBlock (heightB*n) (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) =>+   Int -> General height ZeroInt a -> General height ZeroInt a+takeColumns k = transpose . takeRows k . transpose+dropColumns k = transpose . dropRows k . transpose+++-- alternative: laswp+reverseRows ::+   (Shape.C width, Storable a, 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+      let n = Shape.size height+      let m = Shape.size width+      fwdPtr <- Call.bool True+      nPtr <- Call.cint n+      mPtr <- Call.cint m+      kPtr <- Call.allocaArray n+      aPtr <- ContT $ withForeignPtr a+      liftIO $ do+         copyBlock (n*m) 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) =>+   General height ZeroInt a -> General height ZeroInt a+reverseColumns = transpose . reverseRows . transpose+++fromRowMajor ::+   (Shape.C height, Shape.C width, Storable a, 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) =>+   General height width a -> Array (height,width) a+toRowMajor (Array (MatrixShape.General order height width) x) =+   let shape = (height, width)+   in case order of+         RowMajor -> Array shape x+         ColumnMajor -> Array.unsafeCreate shape $ \yPtr -> evalContT $ do+            let n = Shape.size width+            let m = Shape.size height+            nPtr <- Call.cint n+            xPtr <- ContT $ withForeignPtr x+            incxPtr <- Call.cint m+            incyPtr <- Call.cint 1+            liftIO $ sequence_ $ take m $+               zipWith+                  (\xkPtr ykPtr ->+                     BlasGen.copy nPtr xkPtr incxPtr ykPtr incyPtr)+                  (pointerSeq 1 xPtr)+                  (pointerSeq n yPtr)++flatten ::+   (Shape.C height, Shape.C width, Storable a, 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+++infixl 3 |||+infixl 2 ===++(|||) ::+   (Shape.C height, Eq height, Shape.C widtha, Shape.C widthb,+    Storable a, Class.Floating a) =>+   General height widtha a ->+   General height widthb a ->+   General height (widtha:+:widthb) a+(|||)+      (Array (MatrixShape.General orderA heightA widthA) a)+      (Array (MatrixShape.General orderB heightB widthB) b) =+   if heightA /= heightB+      then error "(|||): mismatching heights"+      else+         case (orderA,orderB) of+            (RowMajor,RowMajor) ->+               Array.unsafeCreate+                  (MatrixShape.General RowMajor heightA (widthA:+:widthB)) $+               \cPtr -> evalContT $ do+                  let n = Shape.size heightA+                  let ma = Shape.size widthA+                  let mb = Shape.size widthB+                  let m = ma+mb+                  maPtr <- Call.cint ma+                  mbPtr <- Call.cint mb+                  aPtr <- ContT $ withForeignPtr a+                  bPtr <- ContT $ withForeignPtr b+                  incxPtr <- Call.cint 1+                  incyPtr <- Call.cint 1+                  liftIO $+                     sequence_ $ take n $+                     zipWith3+                        (\akPtr bkPtr ckPtr -> do+                           BlasGen.copy maPtr akPtr incxPtr ckPtr incyPtr+                           BlasGen.copy mbPtr bkPtr incxPtr+                              (ckPtr `advancePtr` ma) incyPtr)+                        (pointerSeq ma aPtr)+                        (pointerSeq mb bPtr)+                        (pointerSeq m cPtr)+            (RowMajor,ColumnMajor) ->+               Array.unsafeCreate+                  (MatrixShape.General ColumnMajor heightA (widthA:+:widthB)) $+               \cPtr -> evalContT $ do+                  let n = Shape.size heightA+                  let ma = Shape.size widthA+                  let mb = Shape.size widthB+                  aPtr <- ContT $ withForeignPtr a+                  bPtr <- ContT $ withForeignPtr b+                  liftIO $ do+                     copyTransposed n ma aPtr n cPtr+                     copyBlock (n*mb) bPtr (advancePtr cPtr (n*ma))+            (ColumnMajor,RowMajor) ->+               Array.unsafeCreate+                  (MatrixShape.General ColumnMajor heightA (widthA:+:widthB)) $+               \cPtr -> evalContT $ do+                  let n = Shape.size heightA+                  let ma = Shape.size widthA+                  let mb = Shape.size widthB+                  let volA = n*ma+                  aPtr <- ContT $ withForeignPtr a+                  bPtr <- ContT $ withForeignPtr b+                  liftIO $ do+                     copyBlock volA aPtr cPtr+                     copyTransposed n mb bPtr n (advancePtr cPtr volA)+            (ColumnMajor,ColumnMajor) ->+               Array.unsafeCreate+                  (MatrixShape.General ColumnMajor heightA (widthA:+:widthB)) $+               \cPtr -> evalContT $ do+                  let n = Shape.size heightA+                  let na = n * Shape.size widthA+                  let nb = n * Shape.size widthB+                  naPtr <- Call.cint na+                  nbPtr <- Call.cint nb+                  aPtr <- ContT $ withForeignPtr a+                  bPtr <- ContT $ withForeignPtr b+                  incxPtr <- Call.cint 1+                  incyPtr <- Call.cint 1+                  liftIO $ do+                     BlasGen.copy naPtr aPtr incxPtr cPtr incyPtr+                     BlasGen.copy nbPtr bPtr incxPtr+                        (cPtr `advancePtr` na) incyPtr++(===) ::+   (Shape.C width, Eq width, Shape.C heighta, Shape.C heightb,+    Storable a, Class.Floating a) =>+   General heighta width a ->+   General heightb width a ->+   General (heighta:+:heightb) width a+(===) a b = transpose (transpose a ||| transpose b)+++rowSums ::+   (Shape.C height, Shape.C width, Storable a, 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) =>+   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) =>+   Vector height a -> General height width a -> General height width a+scaleRows+   (Array heightX x) (Array shape@(MatrixShape.General order height width) a) =+      Array.unsafeCreate shape $ \bPtr -> do+   Call.assert "scaleRows: sizes mismatch" (heightX == height)+   case order of+      RowMajor -> evalContT $ do+         let m = Shape.size height+         let n = Shape.size width+         alphaPtr <- Call.alloca+         nPtr <- Call.cint n+         xPtr <- ContT $ withForeignPtr x+         aPtr <- ContT $ withForeignPtr a+         incaPtr <- Call.cint 1+         incbPtr <- Call.cint 1+         liftIO $ sequence_ $ take m $+            zipWith3+               (\xkPtr akPtr bkPtr -> do+                  poke alphaPtr =<< peek xkPtr+                  BlasGen.copy nPtr akPtr incaPtr bkPtr incbPtr+                  BlasGen.scal nPtr alphaPtr bkPtr incbPtr)+               (pointerSeq 1 xPtr)+               (pointerSeq n aPtr)+               (pointerSeq n bPtr)+      ColumnMajor -> evalContT $ do+         let m = Shape.size width+         let n = Shape.size height+         transPtr <- Call.char 'N'+         nPtr <- Call.cint n+         klPtr <- Call.cint 0+         kuPtr <- Call.cint 0+         alphaPtr <- Call.number one+         xPtr <- ContT $ withForeignPtr x+         ldxPtr <- Call.cint 1+         aPtr <- ContT $ withForeignPtr a+         incaPtr <- Call.cint 1+         betaPtr <- Call.number zero+         incbPtr <- Call.cint 1+         liftIO $ sequence_ $ take m $+            zipWith+               (\akPtr bkPtr ->+                  BlasGen.gbmv transPtr+                     nPtr nPtr klPtr kuPtr alphaPtr xPtr ldxPtr+                     akPtr incaPtr betaPtr bkPtr incbPtr)+               (pointerSeq n aPtr)+               (pointerSeq n bPtr)++scaleColumns ::+   (Shape.C height, Shape.C width, Eq width, Storable a, 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/Shape.hs view
@@ -0,0 +1,5 @@+module Numeric.LAPACK.Matrix.Shape (+   General,+   ) where++import Numeric.LAPACK.Matrix.Shape.Private
+ src/Numeric/LAPACK/Matrix/Shape/Private.hs view
@@ -0,0 +1,127 @@+{-# LANGUAGE TypeFamilies #-}+module Numeric.LAPACK.Matrix.Shape.Private where++import qualified Data.Array.Comfort.Shape as Shape++import Data.Tuple.HT (swap)+++data Order = RowMajor | ColumnMajor+   deriving (Eq, Show)++flipOrder :: Order -> Order+flipOrder RowMajor = ColumnMajor+flipOrder ColumnMajor = RowMajor++charFromOrder :: Order -> Char+charFromOrder RowMajor = 'T'+charFromOrder ColumnMajor = 'N'+++data General height width =+   General {+      generalOrder :: Order,+      generalHeight :: height,+      generalWidth :: width+   } deriving (Eq, Show)++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)++   offset (General RowMajor height width) =+      Shape.offset (height,width)+   offset (General ColumnMajor height width) =+      Shape.offset (width,height) . swap+   uncheckedOffset (General RowMajor height width) =+      Shape.uncheckedOffset (height,width)+   uncheckedOffset (General ColumnMajor height width) =+      Shape.uncheckedOffset (width,height) . swap++   sizeOffset (General RowMajor height width) =+      Shape.sizeOffset (height,width)+   sizeOffset (General ColumnMajor height width) =+      Shape.sizeOffset (width,height) . swap+   uncheckedSizeOffset (General RowMajor height width) =+      Shape.uncheckedSizeOffset (height,width)+   uncheckedSizeOffset (General ColumnMajor height width) =+      Shape.uncheckedSizeOffset (width,height) . swap++   inBounds (General _ height width) = Shape.inBounds (height,width)+   size (General _ height width) = Shape.size (height,width)+   uncheckedSize (General _ height width) = Shape.uncheckedSize (height,width)+++transpose :: General height width -> General width height+transpose (General order height width) = General (flipOrder order) width height++dimensions ::+   (Shape.C height, Shape.C width) => General height width -> (Int, Int)+dimensions (General order height width) =+   case order of+      RowMajor -> (Shape.size width, Shape.size height)+      ColumnMajor -> (Shape.size height, Shape.size width)+++data Householder height width =+   Householder {+      householderOrder :: Order,+      householderHeight :: height,+      householderWidth :: width+   } deriving (Eq, Show)++data Reflector = Reflector deriving (Eq)+data Triangular = Triangular deriving (Eq)++householderPart ::+   (Shape.C height, Shape.C width) =>+   Householder height width ->+   (Shape.Index height, Shape.Index width) -> Either Reflector Triangular+householderPart (Householder _ height width) (r,c) =+   if Shape.offset height r > Shape.offset width c+     then Left Reflector+     else Right Triangular++instance+   (Shape.C height, Shape.C width) =>+      Shape.C (Householder height width) where++   type Index (Householder height width) =+            (Either Reflector Triangular,+             (Shape.Index height, Shape.Index width))++   indices sh@(Householder _ height width) =+      map (\ix -> (householderPart sh ix, ix)) $+      Shape.indices (height,width)++   offset sh@(Householder order height width) (part,ix) =+      if part == householderPart sh ix+        then+            case order of+               RowMajor -> Shape.offset (height,width) ix+               ColumnMajor -> Shape.offset (width,height) (swap ix)+        else error "Shape.Householder.offset: wrong matrix part"+   uncheckedOffset (Householder RowMajor height width) =+      Shape.uncheckedOffset (height,width) . snd+   uncheckedOffset (Householder ColumnMajor height width) =+      Shape.uncheckedOffset (width,height) . swap . snd++   sizeOffset sh@(Householder order height width) (part,ix) =+      if part == householderPart sh ix+        then+            case order of+               RowMajor -> Shape.sizeOffset (height,width) ix+               ColumnMajor -> Shape.sizeOffset (width,height) (swap ix)+        else error "Shape.Householder.sizeOffset: wrong matrix part"+   uncheckedSizeOffset (Householder RowMajor height width) =+      Shape.uncheckedSizeOffset (height,width) . snd+   uncheckedSizeOffset (Householder ColumnMajor height width) =+      Shape.uncheckedSizeOffset (width,height) . swap . snd++   size (Householder _ height width) = Shape.size (height,width)+   uncheckedSize (Householder _ height width) =+      Shape.uncheckedSize (height,width)+   inBounds sh@(Householder _ height width) (part,ix) =+      Shape.inBounds (height,width) ix+      &&+      part == householderPart sh ix
+ src/Numeric/LAPACK/Private.hs view
@@ -0,0 +1,137 @@+{-# LANGUAGE TypeFamilies #-}+module Numeric.LAPACK.Private where++import qualified Numeric.LAPACK.FFI.Generic as LapackGen+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.Ptr (Ptr)+import Foreign.Storable (Storable, peek)++import Control.Monad.Trans.Cont (evalContT)+import Control.Monad.IO.Class (liftIO)++import Data.Functor.Identity (Identity(Identity, runIdentity))++import Data.Complex (Complex)++import Prelude hiding (sum)+++type family RealOf x++type instance RealOf Float = Float+type instance RealOf Double = Double+type instance RealOf (Complex Float) = Float+type instance RealOf (Complex Double) = Double+++zero, one, minusOne :: Class.Floating a => a+zero =+   runIdentity $+   Class.switchFloating (Identity 0) (Identity 0) (Identity 0) (Identity 0)+one =+   runIdentity $+   Class.switchFloating (Identity 1) (Identity 1) (Identity 1) (Identity 1)+minusOne =+   runIdentity $+   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+   srcPtr <- Call.number a+   incxPtr <- Call.cint 0+   incyPtr <- Call.cint 1+   liftIO $ BlasGen.copy nPtr srcPtr incxPtr dstPtr incyPtr+++copyBlock :: (Class.Floating a) => Int -> Ptr a -> Ptr a -> IO ()+copyBlock n srcPtr dstPtr = evalContT $ do+   nPtr <- Call.cint n+   incxPtr <- Call.cint 1+   incyPtr <- Call.cint 1+   liftIO $ BlasGen.copy nPtr srcPtr incxPtr dstPtr incyPtr++{- |+In ColumnMajor:+Copy a m-by-n-matrix with lda>=m and ldb>=m.+-}+copySubMatrix ::+   (Storable 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'+   mPtr <- Call.cint m+   nPtr <- Call.cint n+   ldaPtr <- Call.cint lda+   ldbPtr <- Call.cint ldb+   liftIO $ LapackGen.lacpy uploPtr mPtr nPtr aPtr ldaPtr bPtr ldbPtr++copyTransposed ::+   (Storable a, Class.Floating a) =>+   Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()+copyTransposed n m aPtr ldb bPtr = evalContT $ do+   nPtr <- Call.cint n+   incaPtr <- Call.cint m+   incbPtr <- Call.cint 1+   liftIO $ sequence_ $ take m $+      zipWith+         (\akPtr bkPtr -> BlasGen.copy nPtr akPtr incaPtr bkPtr incbPtr)+         (pointerSeq 1 aPtr)+         (pointerSeq ldb bPtr)++pointerSeq :: (Storable a) => Int -> Ptr a -> [Ptr a]+pointerSeq k ptr = iterate (flip advancePtr k) ptr+++newtype Sum a = Sum {runSum :: Int -> Ptr a -> Int -> IO a}++sum :: Class.Floating a => Int -> Ptr a -> Int -> IO a+sum =+   runSum $+   Class.switchFloating+      (Sum sumReal)+      (Sum sumReal)+      (Sum sumComplex)+      (Sum sumComplex)++sumReal :: Class.Real a => Int -> Ptr a -> Int -> IO a+sumReal n xPtr incx =+   evalContT $ do+      nPtr <- Call.cint n+      incxPtr <- Call.cint incx+      yPtr <- Call.real oneReal+      incyPtr <- Call.cint 0+      liftIO $ BlasReal.dot nPtr xPtr incxPtr yPtr incyPtr++sumComplex :: Class.Real a => Int -> Ptr (Complex a) -> Int -> IO (Complex a)+sumComplex n xPtr incx =+   evalContT $ do+      transPtr <- Call.char 'N'+      mPtr <- Call.cint 1+      nPtr <- Call.cint n+      alphaPtr <- Call.number one+      onePtr <- Call.number one+      zeroincPtr <- Call.cint 0+      aPtr <- Call.allocaArray n+      ldaPtr <- Call.cint 1+      incxPtr <- Call.cint incx+      betaPtr <- Call.number zero+      yPtr <- Call.alloca+      incyPtr <- Call.cint 1+      liftIO $ BlasGen.copy nPtr onePtr zeroincPtr aPtr incyPtr+      liftIO $+         BlasGen.gemv+            transPtr mPtr nPtr alphaPtr aPtr ldaPtr+            xPtr incxPtr betaPtr yPtr incyPtr+      liftIO $ peek yPtr
+ src/Numeric/LAPACK/Vector.hs view
@@ -0,0 +1,369 @@+module Numeric.LAPACK.Vector (+   Vector,+   fromList,+   constant,+   dot,+   sum,+   absSum,+   norm1,+   norm2,+   argAbsMaximum,+   scale,+   add, sub ,+   mac,+   mul,+   outer,+   conjugate,+   random, RandomDistribution(..),+   ) where++import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape+import qualified Numeric.LAPACK.Private as Private+import Numeric.LAPACK.Private (RealOf, zero, one, minusOne, fill)++import qualified Numeric.LAPACK.FFI.Generic as LapackGen+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 Foreign.ForeignPtr (withForeignPtr)+import Foreign.Ptr (Ptr)+import Foreign.Storable (Storable, peek, peekElemOff, pokeElemOff)+import Foreign.C.Types (CInt)++import System.IO.Unsafe (unsafePerformIO)++import Control.Monad.Trans.Cont (ContT(ContT), evalContT)+import Control.Monad.IO.Class (liftIO)+import Control.Applicative (Const(Const,getConst), (<$>))++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 Data.Complex (Complex)+import Data.Word (Word64)+import Data.Bits (shiftR, (.&.))++import Prelude hiding (sum)+++type Vector = Array+++fromList :: (Shape.C sh, Storable a) => sh -> [a] -> Vector sh a+fromList = Array.fromList+++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)+++newtype Dot sh a = Dot {runDot :: Vector sh a -> Vector sh a -> a}++dot ::+   (Shape.C sh, Eq sh, Class.Floating a) =>+   Vector sh a -> Vector sh a -> a+dot =+   runDot $+   Class.switchFloating+      (Dot dotReal)+      (Dot dotReal)+      (Dot dotComplex)+      (Dot dotComplex)++dotReal ::+   (Shape.C sh, Eq sh, Class.Real a) =>+   Vector sh a -> Vector sh a -> a+dotReal (Array shX x) (Array shY y) = unsafePerformIO $ do+   Call.assert "dot: shapes mismatch" (shX == shY)+   evalContT $ do+      nPtr <- Call.cint $ Shape.size shX+      sxPtr <- ContT $ withForeignPtr x+      incxPtr <- Call.cint 1+      syPtr <- ContT $ withForeignPtr y+      incyPtr <- Call.cint 1+      liftIO $ BlasReal.dot nPtr sxPtr incxPtr syPtr incyPtr++{-+We cannot use 'cdot' because Haskell's FFI+does not support Complex numbers as return values.+-}+dotComplex ::+   (Shape.C sh, Eq sh, Class.Real a) =>+   Vector sh (Complex a) -> Vector sh (Complex a) -> Complex a+dotComplex (Array shX x) (Array shY y) = unsafePerformIO $ do+   Call.assert "dot: shapes mismatch" (shX == shY)+   evalContT $ do+      transPtr <- Call.char 'N'+      mPtr <- Call.cint 1+      nPtr <- Call.cint $ Shape.size shX+      alphaPtr <- Call.number one+      xPtr <- ContT $ withForeignPtr x+      ldxPtr <- Call.cint 1+      yPtr <- ContT $ withForeignPtr y+      incyPtr <- Call.cint 1+      betaPtr <- Call.number zero+      zPtr <- Call.alloca+      inczPtr <- Call.cint 1+      liftIO $+         BlasGen.gemv+            transPtr mPtr nPtr alphaPtr xPtr ldxPtr+            yPtr incyPtr betaPtr zPtr inczPtr+      liftIO $ peek zPtr++sum :: (Shape.C sh, Class.Floating a) => Vector sh a -> a+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 (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++{- |+Sum of the absolute values of real numbers or components of complex numbers.+For real numbers it is equivalent to 'norm1'.+-}+absSum :: (Shape.C sh, Class.Floating a) => Vector sh a -> RealOf a+absSum (Array sh x) = unsafePerformIO $+   evalContT $ do+      nPtr <- Call.cint $ Shape.size sh+      sxPtr <- ContT $ withForeignPtr x+      incxPtr <- Call.cint 1+      liftIO $ asum nPtr sxPtr incxPtr++asum :: Class.Floating a => Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a)+asum =+   getNorm $+   Class.switchFloating+      (Norm BlasReal.asum) (Norm BlasReal.asum)+      (Norm BlasComplex.casum) (Norm BlasComplex.casum)+++{- |+Euclidean norm of a vector or Frobenius norm of a matrix.+-}+norm2 :: (Shape.C sh, Class.Floating a) => Vector sh a -> RealOf a+norm2 (Array sh x) = unsafePerformIO $+   evalContT $ do+      nPtr <- Call.cint $ Shape.size sh+      sxPtr <- ContT $ withForeignPtr x+      incxPtr <- Call.cint 1+      liftIO $ nrm2 nPtr sxPtr incxPtr++nrm2 :: Class.Floating a => Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a)+nrm2 =+   getNorm $+   Class.switchFloating+      (Norm BlasReal.nrm2) (Norm BlasReal.nrm2)+      (Norm BlasComplex.cnrm2) (Norm BlasComplex.cnrm2)++newtype Norm a =+   Norm {getNorm :: Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a)}++{- |+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) =>+   Vector sh a -> (Shape.Index sh, a)+argAbsMaximum (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)+++scale, _scale ::+   (Shape.C sh, Storable a, Class.Floating a) =>+   a -> Vector sh a -> Vector sh a+scale alpha (Array sh x) = Array.unsafeCreate sh $ \syPtr -> do+   evalContT $ do+      alphaPtr <- Call.number alpha+      nPtr <- Call.cint $ Shape.size sh+      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+   let m = 1+   let k = 1+   let n = Shape.size sh+   evalContT $ do+      transaPtr <- Call.char 'N'+      transbPtr <- Call.char 'N'+      mPtr <- Call.cint m+      kPtr <- Call.cint k+      nPtr <- Call.cint n+      alphaPtr <- Call.number one+      aPtr <- Call.number a+      ldaPtr <- Call.cint m+      bPtr <- ContT $ withForeignPtr b+      ldbPtr <- Call.cint k+      betaPtr <- Call.number zero+      ldcPtr <- Call.cint m+      liftIO $+         BlasGen.gemm+            transaPtr transbPtr mPtr nPtr kPtr alphaPtr+            aPtr ldaPtr bPtr ldbPtr betaPtr cPtr ldcPtr++add, sub ::+   (Shape.C sh, Eq sh, Storable a, 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) =>+   a -> Vector sh a -> Vector sh a -> Vector sh a+mac alpha (Array shX x) (Array shY y) = Array.unsafeCreate shX $ \szPtr -> do+   Call.assert "mac: shapes mismatch" (shX == shY)+   evalContT $ do+      nPtr <- Call.cint $ Shape.size shX+      saPtr <- Call.number alpha+      sxPtr <- ContT $ withForeignPtr x+      incxPtr <- Call.cint 1+      syPtr <- ContT $ withForeignPtr y+      incyPtr <- Call.cint 1+      inczPtr <- Call.cint 1+      liftIO $ BlasGen.copy nPtr syPtr incyPtr szPtr inczPtr+      liftIO $ BlasGen.axpy nPtr saPtr sxPtr incxPtr szPtr inczPtr++mul ::+   (Shape.C sh, Eq sh, Storable a, Class.Floating a) =>+   Vector sh a -> Vector sh a -> Vector sh a+mul (Array shA a) (Array shX x) = Array.unsafeCreate shX $ \yPtr -> do+   Call.assert "mul: shapes mismatch" (shA == shX)+   let n = Shape.size shX+   evalContT $ do+      transPtr <- Call.char 'N'+      mPtr <- Call.cint n+      nPtr <- Call.cint n+      klPtr <- Call.cint 0+      kuPtr <- Call.cint 0+      alphaPtr <- Call.number one+      aPtr <- ContT $ withForeignPtr a+      ldaPtr <- Call.cint 1+      xPtr <- ContT $ withForeignPtr x+      incxPtr <- Call.cint 1+      betaPtr <- Call.number zero+      incyPtr <- Call.cint 1+      liftIO $+         BlasGen.gbmv transPtr+            mPtr nPtr klPtr kuPtr alphaPtr aPtr ldaPtr+            xPtr incxPtr betaPtr yPtr incyPtr+++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+outer (Array shX x) (Array shY y) =+   Array.unsafeCreate (MatrixShape.General MatrixShape.ColumnMajor shX shY) $+      \cPtr -> do+   let m = Shape.size shX+   let n = Shape.size shY+   evalContT $ do+      transaPtr <- Call.char 'N'+      transbPtr <- Call.char 'N'+      mPtr <- Call.cint m+      nPtr <- Call.cint n+      kPtr <- Call.cint 1+      alphaPtr <- Call.number one+      aPtr <- ContT $ withForeignPtr x+      ldaPtr <- Call.cint m+      bPtr <- ContT $ withForeignPtr y+      ldbPtr <- Call.cint 1+      betaPtr <- Call.number zero+      ldcPtr <- Call.cint m+      liftIO $+         BlasGen.gemm+            transaPtr transbPtr mPtr nPtr kPtr alphaPtr+            aPtr ldaPtr bPtr ldbPtr betaPtr cPtr ldcPtr+++newtype Conjugate sh a = Conjugate {getConjugate :: Vector sh a -> Vector sh a}++conjugate ::+   (Shape.C sh, Storable a, Class.Floating a) =>+   Vector sh a -> Vector sh a+conjugate =+   getConjugate $+   Class.switchFloating+      (Conjugate id)+      (Conjugate id)+      (Conjugate complexConjugate)+      (Conjugate complexConjugate)++complexConjugate ::+   (Shape.C sh, Storable a, Class.Real a) =>+   Vector sh (Complex a) -> Vector sh (Complex a)+complexConjugate (Array sh x) = Array.unsafeCreate sh $ \syPtr ->+   evalContT $ do+      nPtr <- Call.cint $ Shape.size sh+      sxPtr <- ContT $ withForeignPtr x+      incxPtr <- Call.cint 1+      incyPtr <- Call.cint 1+      liftIO $ do+         BlasGen.copy nPtr sxPtr incxPtr syPtr incyPtr+         LapackComplex.lacgv nPtr syPtr incyPtr+++data RandomDistribution =+     UniformBox01+   | UniformBoxPM1+   | Normal+   | UniformDisc+   | UniformCircle+   deriving (Eq, Ord, Show, Enum)++{-+@random distribution shape seed@++Only the least significant 47 bits of @seed@ are used.+-}+random ::+   (Shape.C sh, Storable a, Class.Floating a) =>+   RandomDistribution -> sh -> Word64 -> Vector sh a+random dist sh seed = Array.unsafeCreate sh $ \xPtr ->+   evalContT $ do+      nPtr <- Call.cint $ Shape.size sh+      distPtr <-+         Call.cint $+         case (getConst $ isComplexInFunctor xPtr, dist) of+            (_, UniformBox01) -> 1+            (_, UniformBoxPM1) -> 2+            (_, Normal) -> 3+            (True, UniformDisc) -> 4+            (True, UniformCircle) -> 5+            (False, UniformDisc) -> 2+            (False, UniformCircle) ->+               error+                  "Vector.random: UniformCircle not supported for real numbers"+      iseedPtr <- Call.allocaArray 4+      liftIO $ do+         pokeElemOff iseedPtr 0 $ fromIntegral ((seed `shiftR` 35) .&. 0xFFF)+         pokeElemOff iseedPtr 1 $ fromIntegral ((seed `shiftR` 23) .&. 0xFFF)+         pokeElemOff iseedPtr 2 $ fromIntegral ((seed `shiftR` 11) .&. 0xFFF)+         pokeElemOff iseedPtr 3 $ fromIntegral ((seed.&.0x7FF)*2+1)+         LapackGen.larnv distPtr iseedPtr nPtr xPtr++isComplexInFunctor :: (Class.Floating a) => f a -> Const Bool a+isComplexInFunctor _ =+   Class.switchFloating (Const False) (Const False) (Const True) (Const True)