lapack 0.3.0.1 → 0.3.1
raw patch · 11 files changed
+659/−44 lines, 11 files
Files
- Changes.md +6/−0
- lapack.cabal +5/−2
- src/Numeric/LAPACK/Matrix/Hermitian/Basic.hs +4/−26
- src/Numeric/LAPACK/Matrix/Symmetric.hs +178/−0
- src/Numeric/LAPACK/Matrix/Symmetric/Basic.hs +209/−0
- src/Numeric/LAPACK/Matrix/Triangular/Basic.hs +1/−13
- src/Numeric/LAPACK/Matrix/Triangular/Private.hs +14/−0
- src/Numeric/LAPACK/Split.hs +34/−1
- test/Main.hs +2/−0
- test/Test/Symmetric.hs +205/−0
- test/Test/Triangular.hs +1/−2
Changes.md view
@@ -1,5 +1,11 @@ # Change log for the `lapack` package +## 0.3.1++ * `Matrix.Symmetric`:+ You can now import many functions for symmetric matrices from this module.+ This is more natural than importing them from `Triangular`.+ ## 0.3 * Matrix data family
lapack.cabal view
@@ -1,5 +1,5 @@ Name: lapack-Version: 0.3.0.1+Version: 0.3.1 License: BSD3 License-File: LICENSE Author: Henning Thielemann <haskell@henning-thielemann.de>@@ -48,7 +48,7 @@ Makefile Source-Repository this- Tag: 0.3.0.1+ Tag: 0.3.1 Type: darcs Location: http://hub.darcs.net/thielema/lapack/ @@ -89,6 +89,7 @@ Numeric.LAPACK.Matrix.Square Numeric.LAPACK.Matrix.Hermitian Numeric.LAPACK.Matrix.HermitianPositiveDefinite+ Numeric.LAPACK.Matrix.Symmetric Numeric.LAPACK.Matrix.Triangular Numeric.LAPACK.Matrix.Banded Numeric.LAPACK.Matrix.BandedHermitian@@ -126,6 +127,7 @@ Numeric.LAPACK.Matrix.Hermitian.Linear Numeric.LAPACK.Matrix.Hermitian.Eigen Numeric.LAPACK.Matrix.HermitianPositiveDefinite.Linear+ Numeric.LAPACK.Matrix.Symmetric.Basic Numeric.LAPACK.Matrix.Symmetric.Private Numeric.LAPACK.Matrix.Banded.Basic Numeric.LAPACK.Matrix.Banded.Linear@@ -194,6 +196,7 @@ Test.Matrix Test.Square Test.Triangular+ Test.Symmetric Test.Hermitian Test.Orthogonal Test.LowerUpper
src/Numeric/LAPACK/Matrix/Hermitian/Basic.hs view
@@ -545,7 +545,8 @@ congruence b = skipCheckCongruence Basic.mapWidth $ \a -> scaledAnticommutator one a $- Split.tallMultiplyR NonTransposed (takeHalf b) a+ Split.tallMultiplyR NonTransposed+ (Split.takeHalf MatrixShape.hermitianOrder b) a congruenceAdjoint :: (Shape.C height, Shape.C width, Eq width, Class.Floating a) =>@@ -553,31 +554,8 @@ congruenceAdjoint = flip $ \b -> skipCheckCongruence Basic.mapHeight $ \a -> scaledAnticommutatorAdjoint one a $- Basic.swapMultiply (Split.tallMultiplyR Transposed) a (takeHalf b)---data Corrupt = Corrupt- deriving (Eq)--{- |-> let b = takeHalf a-> ==>-> isTriangular b && a == addAdjoint b--}-takeHalf ::- (Shape.C sh, Class.Floating a) =>- Hermitian sh a -> Split.Square Corrupt sh a-takeHalf (Array (MatrixShape.Hermitian order sh) a) =- Array.unsafeCreate (MatrixShape.Split Corrupt order (Extent.square sh)) $- \bPtr -> evalContT $ do- let n = Shape.size sh- aPtr <- ContT $ withForeignPtr a- nPtr <- Call.cint n- alphaPtr <- Call.number 0.5- incxPtr <- Call.cint (n+1)- liftIO $ do- unpack order n aPtr bPtr- BlasGen.scal nPtr alphaPtr bPtr incxPtr+ Basic.swapMultiply (Split.tallMultiplyR Transposed)+ a (Split.takeHalf MatrixShape.hermitianOrder b) scaledAnticommutator ::
+ src/Numeric/LAPACK/Matrix/Symmetric.hs view
@@ -0,0 +1,178 @@+{-# LANGUAGE TypeOperators #-}+module Numeric.LAPACK.Matrix.Symmetric (+ Symmetric,+ size,+ fromList, autoFromList,+ identity,+ diagonal,+ takeDiagonal,+ transpose,+ adjoint,++ stack, (#%%%#),+ split,++ toSquare,++ gramian, gramianTransposed,+ congruenceDiagonal, congruenceDiagonalTransposed,+ congruence, congruenceTransposed,+ anticommutator, anticommutatorTransposed,+ ) where++import qualified Numeric.LAPACK.Matrix.Symmetric.Basic as Basic+import qualified Numeric.LAPACK.Matrix.Triangular as Triangular++import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix+import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape+import qualified Numeric.LAPACK.Matrix.Extent as Extent+import Numeric.LAPACK.Matrix.Array.Triangular (Symmetric)+import Numeric.LAPACK.Matrix.Array (Full, General, Square)+import Numeric.LAPACK.Matrix.Shape.Private (Order)+import Numeric.LAPACK.Matrix.Private (ShapeInt)+import Numeric.LAPACK.Vector (Vector)+import Numeric.LAPACK.Scalar (one)++import qualified Numeric.Netlib.Class as Class++import qualified Data.Array.Comfort.Shape as Shape+import Data.Array.Comfort.Shape ((:+:))++import Foreign.Storable (Storable)+++size :: Symmetric sh a -> sh+size = MatrixShape.triangularSize . ArrMatrix.shape++transpose :: Symmetric sh a -> Symmetric sh a+transpose = Triangular.transpose++adjoint :: (Shape.C sh, Class.Floating a) => Symmetric sh a -> Symmetric sh a+adjoint = Triangular.adjoint+++fromList :: (Shape.C sh, Storable a) => Order -> sh -> [a] -> Symmetric sh a+fromList = Triangular.symmetricFromList++autoFromList :: (Storable a) => Order -> [a] -> Symmetric ShapeInt a+autoFromList = Triangular.autoSymmetricFromList+++toSquare :: (Shape.C sh, Class.Floating a) => Symmetric sh a -> Square sh a+toSquare = Triangular.toSquare+++identity :: (Shape.C sh, Class.Floating a) => Order -> sh -> Symmetric sh a+identity order = Triangular.relaxUnitDiagonal . Triangular.identity order++diagonal ::+ (Shape.C sh, Class.Floating a) => Order -> Vector sh a -> Symmetric sh a+diagonal = Triangular.diagonal++takeDiagonal :: (Shape.C sh, Class.Floating a) => Symmetric sh a -> Vector sh a+takeDiagonal = Triangular.takeDiagonal+++stack ::+ (Shape.C sh0, Eq sh0, Shape.C sh1, Eq sh1, Class.Floating a) =>+ Symmetric sh0 a ->+ General sh0 sh1 a ->+ Symmetric sh1 a ->+ Symmetric (sh0:+:sh1) a+stack = Triangular.stackSymmetric++infixr 2 #%%%#++(#%%%#) ::+ (Shape.C sh0, Eq sh0, Shape.C sh1, Eq sh1, Class.Floating a) =>+ (Symmetric sh0 a, General sh0 sh1 a) ->+ Symmetric sh1 a ->+ Symmetric (sh0:+:sh1) a+(#%%%#) = uncurry stack+++split ::+ (Shape.C sh0, Eq sh0, Shape.C sh1, Eq sh1, Class.Floating a) =>+ Symmetric (sh0:+:sh1) a ->+ (Symmetric sh0 a, General sh0 sh1 a, Symmetric sh1 a)+split = Triangular.splitSymmetric++++{- |+gramian A = A^T * A+-}+gramian ::+ (Shape.C height, Shape.C width, Class.Floating a) =>+ General height width a -> Symmetric width a+gramian = ArrMatrix.lift1 Basic.gramian++{- |+gramianTransposed A = A * A^T = gramian (A^T)+-}+gramianTransposed ::+ (Shape.C height, Shape.C width, Class.Floating a) =>+ General height width a -> Symmetric height a+gramianTransposed = ArrMatrix.lift1 Basic.gramianTransposed++{- |+congruenceDiagonal D A = A^T * D * A+-}+congruenceDiagonal ::+ (Shape.C height, Eq height, Shape.C width, Class.Floating a) =>+ Vector height a -> General height width a -> Symmetric width a+congruenceDiagonal = ArrMatrix.lift1 . Basic.congruenceDiagonal++{- |+congruenceDiagonalTransposed A D = A * D * A^T+-}+congruenceDiagonalTransposed ::+ (Shape.C height, Shape.C width, Eq width, Class.Floating a) =>+ General height width a -> Vector width a -> Symmetric height a+congruenceDiagonalTransposed a =+ ArrMatrix.lift0 . Basic.congruenceDiagonalTransposed (ArrMatrix.toVector a)++{- |+congruence B A = A^T * B * A+-}+congruence ::+ (Shape.C height, Eq height, Shape.C width, Class.Floating a) =>+ Symmetric height a -> General height width a -> Symmetric width a+congruence = ArrMatrix.lift2 Basic.congruence++{- |+congruenceTransposed B A = A * B * A^T+-}+congruenceTransposed ::+ (Shape.C height, Shape.C width, Eq width, Class.Floating a) =>+ General height width a -> Symmetric width a -> Symmetric height a+congruenceTransposed = ArrMatrix.lift2 Basic.congruenceTransposed+++{- |+anticommutator A B = A^T * B + B^T * A++Not exactly a matrix anticommutator,+thus I like to call it Symmetric anticommutator.+-}+anticommutator ::+ (Extent.C vert, Extent.C horiz,+ Shape.C height, Eq height, Shape.C width, Eq width,+ Class.Floating a) =>+ Full vert horiz height width a ->+ Full vert horiz height width a -> Symmetric width a+anticommutator = ArrMatrix.lift2 $ Basic.scaledAnticommutator one++{- |+anticommutatorTransposed A B+ = A * B^T + B * A^T+ = anticommutator (transpose A) (transpose B)+-}+anticommutatorTransposed ::+ (Extent.C vert, Extent.C horiz,+ Shape.C height, Eq height, Shape.C width, Eq width,+ Class.Floating a) =>+ Full vert horiz height width a ->+ Full vert horiz height width a -> Symmetric height a+anticommutatorTransposed =+ ArrMatrix.lift2 $ Basic.scaledAnticommutatorTransposed one
+ src/Numeric/LAPACK/Matrix/Symmetric/Basic.hs view
@@ -0,0 +1,209 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}+module Numeric.LAPACK.Matrix.Symmetric.Basic (+ Symmetric,++ gramian, gramianTransposed,+ congruenceDiagonal, congruenceDiagonalTransposed,+ congruence, congruenceTransposed,+ scaledAnticommutator, scaledAnticommutatorTransposed,+ ) where++import qualified Numeric.LAPACK.Matrix.Private as Matrix+import qualified Numeric.LAPACK.Matrix.Basic as Basic+import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape+import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent+import qualified Numeric.LAPACK.Split as Split+import Numeric.LAPACK.Matrix.Triangular.Private (pack, recheck)+import Numeric.LAPACK.Matrix.Shape.Private+ (Order(RowMajor,ColumnMajor), NonUnit(NonUnit))+import Numeric.LAPACK.Matrix.Modifier+ (Transposition(NonTransposed, Transposed))+import Numeric.LAPACK.Matrix.Private (Full)+import Numeric.LAPACK.Vector (Vector)+import Numeric.LAPACK.Scalar (zero, one)+import Numeric.LAPACK.Shape.Private (Unchecked(Unchecked))++import qualified Numeric.BLAS.FFI.Generic as BlasGen+import qualified Numeric.Netlib.Utility as Call+import qualified Numeric.Netlib.Class as Class++import qualified Data.Array.Comfort.Storable.Unchecked as Array+import qualified Data.Array.Comfort.Shape as Shape+import Data.Array.Comfort.Storable.Unchecked (Array(Array))++import Foreign.ForeignPtr (ForeignPtr, withForeignPtr)+import Foreign.Ptr (Ptr)++import Control.Monad.Trans.Cont (ContT(ContT), evalContT)+import Control.Monad.IO.Class (liftIO)+++type Symmetric sh = Array (MatrixShape.FlexSymmetric NonUnit sh)+++-- cf. Hermitian.Basic+gramian ::+ (Shape.C height, Shape.C width, Class.Floating a) =>+ Matrix.General height width a -> Symmetric width a+gramian (Array (MatrixShape.Full order extent) a) =+ Array.unsafeCreate (symmetricShape order $ Extent.width extent) $+ \bPtr -> gramianIO order a bPtr $ gramianParameters order extent++gramianParameters ::+ (Extent.C horiz, Extent.C vert, Shape.C height, Shape.C width) =>+ Order ->+ Extent.Extent vert horiz height width ->+ ((Int, Int), (Char, Char, Int))+gramianParameters order extent =+ let (height, width) = Extent.dimensions extent+ n = Shape.size width+ k = Shape.size height+ in ((n,k),+ case order of+ ColumnMajor -> ('U', 'T', k)+ RowMajor -> ('L', 'N', n))+++gramianTransposed ::+ (Shape.C height, Shape.C width, Class.Floating a) =>+ Matrix.General height width a -> Symmetric height a+gramianTransposed (Array (MatrixShape.Full order extent) a) =+ Array.unsafeCreate (symmetricShape order $ Extent.height extent) $+ \bPtr -> gramianIO order a bPtr $ gramianTransposedParameters order extent++gramianTransposedParameters ::+ (Extent.C horiz, Extent.C vert, Shape.C height, Shape.C width) =>+ Order ->+ Extent.Extent vert horiz height width ->+ ((Int, Int), (Char, Char, Int))+gramianTransposedParameters order extent =+ let (height, width) = Extent.dimensions extent+ n = Shape.size height+ k = Shape.size width+ in ((n,k),+ case order of+ ColumnMajor -> ('U', 'N', n)+ RowMajor -> ('L', 'T', k))++gramianIO ::+ (Class.Floating a) =>+ Order ->+ ForeignPtr a -> Ptr a ->+ ((Int, Int), (Char, Char, Int)) -> IO ()+gramianIO order a bPtr ((n,k), (uplo,trans,lda)) =+ evalContT $ do+ uploPtr <- Call.char uplo+ transPtr <- Call.char trans+ nPtr <- Call.cint n+ kPtr <- Call.cint k+ alphaPtr <- Call.number one+ aPtr <- ContT $ withForeignPtr a+ ldaPtr <- Call.leadingDim lda+ betaPtr <- Call.number zero+ cPtr <- Call.allocaArray (n*n)+ ldcPtr <- Call.leadingDim n+ liftIO $ do+ BlasGen.syrk uploPtr transPtr+ nPtr kPtr alphaPtr aPtr ldaPtr betaPtr cPtr ldcPtr+ pack order n cPtr bPtr++skipCheckCongruence ::+ ((sh -> Unchecked sh) -> matrix0 -> matrix1) ->+ (matrix1 -> Symmetric (Unchecked sh) a) -> matrix0 -> Symmetric sh a+skipCheckCongruence mapSize f a =+ recheck $ f $ mapSize Unchecked a+++congruenceDiagonal ::+ (Shape.C height, Eq height, Shape.C width, Class.Floating a) =>+ Vector height a -> Matrix.General height width a -> Symmetric width a+congruenceDiagonal d =+ skipCheckCongruence Basic.mapWidth $ \a ->+ scaledAnticommutator 0.5 a $ Basic.scaleRows d a++congruenceDiagonalTransposed ::+ (Shape.C height, Shape.C width, Eq width, Class.Floating a) =>+ Matrix.General height width a -> Vector width a -> Symmetric height a+congruenceDiagonalTransposed =+ flip $ \d -> skipCheckCongruence Basic.mapHeight $ \a ->+ scaledAnticommutatorTransposed 0.5 a $ Basic.scaleColumns d a+++congruence ::+ (Shape.C height, Eq height, Shape.C width, Class.Floating a) =>+ Symmetric height a -> Matrix.General height width a -> Symmetric width a+congruence b =+ skipCheckCongruence Basic.mapWidth $ \a ->+ scaledAnticommutator one a $+ Split.tallMultiplyR NonTransposed+ (Split.takeHalf MatrixShape.triangularOrder b) a++congruenceTransposed ::+ (Shape.C height, Shape.C width, Eq width, Class.Floating a) =>+ Matrix.General height width a -> Symmetric width a -> Symmetric height a+congruenceTransposed =+ flip $ \b -> skipCheckCongruence Basic.mapHeight $ \a ->+ scaledAnticommutatorTransposed one a $+ Basic.swapMultiply (Split.tallMultiplyR Transposed)+ a (Split.takeHalf MatrixShape.triangularOrder b)+++scaledAnticommutator ::+ (Extent.C vert, Extent.C horiz,+ Shape.C height, Eq height, Shape.C width, Eq width, Class.Floating a) =>+ a ->+ Full vert horiz height width a ->+ Full vert horiz height width a -> Symmetric width a+scaledAnticommutator alpha arr (Array (MatrixShape.Full order extentB) b) = do+ let (Array (MatrixShape.Full _ extentA) a) = Basic.forceOrder order arr+ Array.unsafeCreate (symmetricShape order $ Extent.width extentB) $+ \cpPtr -> do+ Call.assert "Symmetric.anticommutator: extents mismatch"+ (extentA==extentB)+ scaledAnticommutatorIO alpha order a b cpPtr $+ gramianParameters order extentB++scaledAnticommutatorTransposed ::+ (Extent.C vert, Extent.C horiz,+ Shape.C height, Eq height, Shape.C width, Eq width, Class.Floating a) =>+ a ->+ Full vert horiz height width a ->+ Full vert horiz height width a -> Symmetric height a+scaledAnticommutatorTransposed+ alpha arr (Array (MatrixShape.Full order extentB) b) = do+ let (Array (MatrixShape.Full _ extentA) a) = Basic.forceOrder order arr+ Array.unsafeCreate (symmetricShape order $ Extent.height extentB) $+ \cpPtr -> do+ Call.assert "Symmetric.anticommutatorTransposed: extents mismatch"+ (extentA==extentB)+ scaledAnticommutatorIO alpha order a b cpPtr $+ gramianTransposedParameters order extentB++scaledAnticommutatorIO ::+ (Class.Floating a) =>+ a ->+ Order -> ForeignPtr a -> ForeignPtr a -> Ptr a ->+ ((Int, Int), (Char, Char, Int)) -> IO ()+scaledAnticommutatorIO alpha order a b cpPtr ((n,k), (uplo,trans,lda)) =+ evalContT $ do+ uploPtr <- Call.char uplo+ transPtr <- Call.char trans+ nPtr <- Call.cint n+ kPtr <- Call.cint k+ alphaPtr <- Call.number alpha+ aPtr <- ContT $ withForeignPtr a+ ldaPtr <- Call.leadingDim lda+ bPtr <- ContT $ withForeignPtr b+ let ldbPtr = ldaPtr+ betaPtr <- Call.number zero+ cPtr <- Call.allocaArray (n*n)+ ldcPtr <- Call.leadingDim n+ liftIO $ do+ BlasGen.syr2k uploPtr transPtr nPtr kPtr alphaPtr+ aPtr ldaPtr bPtr ldbPtr betaPtr cPtr ldcPtr+ pack order n cPtr cpPtr+++symmetricShape :: Order -> size -> MatrixShape.Symmetric size+symmetricShape = MatrixShape.Triangular NonUnit MatrixShape.autoUplo
src/Numeric/LAPACK/Matrix/Triangular/Basic.hs view
@@ -52,7 +52,7 @@ import qualified Numeric.LAPACK.Vector as Vector import Numeric.LAPACK.Matrix.Triangular.Private (Triangular, FlexDiagonal, diagonalPointers, diagonalPointerPairs,- pack, packRect, unpack, unpackZero, unpackToTemp)+ pack, packRect, unpack, unpackZero, unpackToTemp, uncheck, recheck) import Numeric.LAPACK.Matrix.Shape.Private (Order(RowMajor,ColumnMajor), flipOrder, transposeFromOrder, uploFromOrder, uploOrder,@@ -587,18 +587,6 @@ Array.unsafeCreate shape $ Symmetric.square NonConjugated order (Shape.size sh) a --uncheck :: Triangular lo diag up sh a -> Triangular lo diag up (Unchecked sh) a-uncheck =- Array.mapShape $- \(MatrixShape.Triangular diag uplo order sh) ->- MatrixShape.Triangular diag uplo order (Unchecked sh)--recheck :: Triangular lo diag up (Unchecked sh) a -> Triangular lo diag up sh a-recheck =- Array.mapShape $- \(MatrixShape.Triangular diag uplo order (Unchecked sh)) ->- MatrixShape.Triangular diag uplo order sh {- Requires frequent unpacking and packing of triangles.
src/Numeric/LAPACK/Matrix/Triangular/Private.hs view
@@ -13,6 +13,7 @@ import Numeric.LAPACK.Matrix.Modifier (Conjugation(NonConjugated)) import Numeric.LAPACK.Matrix.Private (Full) import Numeric.LAPACK.Scalar (zero)+import Numeric.LAPACK.Shape.Private (Unchecked(Unchecked)) import Numeric.LAPACK.Private (pointerSeq, copyBlock, copyCondConjugateToTemp, pokeCInt, fill, withInfo, errorCodeMsg)@@ -138,6 +139,19 @@ zPtr <- Call.number z liftIO $ LapackGen.laset uploPtr nPtr nPtr zPtr zPtr aPtr nPtr +++uncheck :: Triangular lo diag up sh a -> Triangular lo diag up (Unchecked sh) a+uncheck =+ Array.mapShape $+ \(MatrixShape.Triangular diag uplo order sh) ->+ MatrixShape.Triangular diag uplo order (Unchecked sh)++recheck :: Triangular lo diag up (Unchecked sh) a -> Triangular lo diag up sh a+recheck =+ Array.mapShape $+ \(MatrixShape.Triangular diag uplo order (Unchecked sh)) ->+ MatrixShape.Triangular diag uplo order sh stack ::
src/Numeric/LAPACK/Split.hs view
@@ -1,12 +1,14 @@+{-# LANGUAGE TypeFamilies #-} module Numeric.LAPACK.Split where import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape+import qualified Numeric.LAPACK.Matrix.Shape.Box as Box import qualified Numeric.LAPACK.Matrix.Triangular.Private as TriPriv import qualified Numeric.LAPACK.Matrix.Triangular.Basic as Tri import qualified Numeric.LAPACK.Matrix.Private as Matrix import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent import qualified Numeric.LAPACK.Private as Private-import Numeric.LAPACK.Matrix.Triangular.Private (diagonalPointers)+import Numeric.LAPACK.Matrix.Triangular.Private (diagonalPointers, unpack) import Numeric.LAPACK.Matrix.Triangular.Basic (UnitLower, Upper) import Numeric.LAPACK.Matrix.Shape.Private (Order(RowMajor, ColumnMajor), transposeFromOrder,@@ -238,3 +240,34 @@ withInfo "trtrs" $ LapackGen.trtrs uploPtr transPtr diagPtr nPtr nrhsPtr aPtr ldaPtr xPtr ldxPtr+++data Corrupt = Corrupt+ deriving (Eq)++{-+We could use Plain.Class.ShapeOrder+but this would currently cause an import cycle.+-}+{- |+> let b = takeHalf a+> ==>+> isTriangular b && a == addTransposed b+-}+takeHalf ::+ (Box.Box symShape, Box.HeightOf symShape ~ sh, Shape.C sh,+ Class.Floating a) =>+ (symShape -> Order) -> Array symShape a -> Square Corrupt sh a+takeHalf shapeOrder (Array symShape a) =+ let sh = Box.height symShape+ order = shapeOrder symShape+ in Array.unsafeCreate (MatrixShape.Split Corrupt order (Extent.square sh)) $+ \bPtr -> evalContT $ do+ let n = Shape.size sh+ aPtr <- ContT $ withForeignPtr a+ nPtr <- Call.cint n+ alphaPtr <- Call.number 0.5+ incxPtr <- Call.cint (n+1)+ liftIO $ do+ unpack order n aPtr bPtr+ BlasGen.scal nPtr alphaPtr bPtr incxPtr
test/Main.hs view
@@ -5,6 +5,7 @@ import qualified Test.Matrix as Matrix import qualified Test.Square as Square import qualified Test.Triangular as Triangular+import qualified Test.Symmetric as Symmetric import qualified Test.Hermitian as Hermitian import qualified Test.Banded as Banded import qualified Test.BandedHermitian as BandedHermitian@@ -39,6 +40,7 @@ prefix "Matrix" Matrix.testsVar ++ prefix "Square" Square.testsVar ++ prefix "Triangular" Triangular.testsVar +++ prefix "Symmetric" Symmetric.testsVar ++ prefix "Hermitian" Hermitian.testsVar ++ prefix "Banded" Banded.testsVar ++ prefix "BandedHermitian" BandedHermitian.testsVar ++
+ test/Test/Symmetric.hs view
@@ -0,0 +1,205 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}+module Test.Symmetric (testsVar) where++import qualified Test.Generator as Gen+import qualified Test.Utility as Util+import Test.Generator ((<-*#>), (<#*|>), (<#*#>), (<#=#>))+import Test.Utility+ (approxArray, approxMatrix, equalArray, Tagged, genOrder, (!===))++import qualified Numeric.LAPACK.Matrix.Symmetric as Symmetric+import qualified Numeric.LAPACK.Matrix.Triangular as Triangular+import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix+import qualified Numeric.LAPACK.Matrix as Matrix+import qualified Numeric.LAPACK.Vector as Vector+import Numeric.LAPACK.Matrix.Symmetric (Symmetric)+import Numeric.LAPACK.Matrix.Shape (Order)+import Numeric.LAPACK.Matrix (General, ShapeInt, (#+#), (|||))+import Numeric.LAPACK.Vector (Vector)+import Numeric.LAPACK.Scalar (RealOf)++import qualified Numeric.Netlib.Class as Class++import qualified Data.Array.Comfort.Shape as Shape+import Data.Array.Comfort.Shape ((:+:))++import Control.Applicative ((<$>))++import Data.Semigroup ((<>))+import Data.Tuple.HT (uncurry3)++import qualified Test.QuickCheck as QC++++generalFromSymmetric ::+ (Shape.C sh, Class.Floating a) => Symmetric sh a -> General sh sh a+generalFromSymmetric = Matrix.fromFull . Triangular.toSquare++stack ::+ (Class.Floating a) =>+ (Symmetric ShapeInt a, General ShapeInt ShapeInt a, Symmetric ShapeInt a) ->+ Bool+stack (a,b,c) =+ let abc = generalFromSymmetric $ Symmetric.stack a b c+ in equalArray abc $+ (Matrix.fromFull (Triangular.toSquare a) ||| b+ !===+ Matrix.transpose b ||| Matrix.fromFull (Triangular.toSquare c))++split :: (Class.Floating a) => Symmetric (ShapeInt:+:ShapeInt) a -> Bool+split abc = equalArray abc $ uncurry3 Symmetric.stack $ Symmetric.split abc+++gramian ::+ (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>+ General ShapeInt ShapeInt a -> Bool+gramian x =+ approxArray+ (generalFromSymmetric $ Symmetric.gramian x)+ (Matrix.transpose x <> x)++gramianTransposed ::+ (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>+ General ShapeInt ShapeInt a -> Bool+gramianTransposed x =+ approxArray+ (generalFromSymmetric $ Symmetric.gramianTransposed x)+ (Matrix.adaptOrder x $ x <> Matrix.transpose x)++gramianNonTransposed ::+ (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>+ General ShapeInt ShapeInt a -> Bool+gramianNonTransposed x =+ approxArray+ (Matrix.forceOrder (ArrMatrix.shapeOrder $ ArrMatrix.shape x) $+ Symmetric.gramian $ Matrix.transpose x)+ (Symmetric.gramianTransposed x)++congruenceDiagonal ::+ (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>+ (Vector ShapeInt a, General ShapeInt ShapeInt a) -> Bool+congruenceDiagonal (d,a) =+ approxArray+ (generalFromSymmetric $ Symmetric.congruenceDiagonal d a)+ (Matrix.transpose a <> Matrix.scaleRows d a)++congruenceDiagonalTransposed ::+ (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>+ (General ShapeInt ShapeInt a, Vector ShapeInt a) -> Bool+congruenceDiagonalTransposed (a,d) =+ approxMatrix 1e-5+ (generalFromSymmetric $ Symmetric.congruenceDiagonalTransposed a d)+ (Matrix.scaleColumns d a <> Matrix.transpose a)++congruenceDiagonalGramian ::+ (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>+ General ShapeInt ShapeInt a -> Bool+congruenceDiagonalGramian a =+ approxArray+ (Symmetric.congruenceDiagonal (Vector.one $ Matrix.height a) a)+ (Symmetric.gramian a)++congruence ::+ (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>+ (Symmetric ShapeInt a, General ShapeInt ShapeInt a) -> Bool+congruence (b,a) =+ approxArray+ (generalFromSymmetric $ Symmetric.congruence b a)+ (Matrix.transpose a <> generalFromSymmetric b <> a)++congruenceTransposed ::+ (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>+ (General ShapeInt ShapeInt a, Symmetric ShapeInt a) -> Bool+congruenceTransposed (a,b) =+ approxMatrix 1e-5+ (generalFromSymmetric $ Symmetric.congruenceTransposed a b)+ (a <> generalFromSymmetric b <> Matrix.transpose a)++congruenceCongruenceDiagonal ::+ (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>+ Order -> (Vector ShapeInt a, General ShapeInt ShapeInt a) -> Bool+congruenceCongruenceDiagonal order (d,a) =+ approxArray+ (Symmetric.congruenceDiagonal d a)+ (Symmetric.congruence (Triangular.diagonal order d) a)++anticommutator ::+ (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>+ (General ShapeInt ShapeInt a, General ShapeInt ShapeInt a) -> Bool+anticommutator (a,b) =+ approxArray+ (generalFromSymmetric $ Symmetric.anticommutator a b)+ ((Matrix.transpose b <> a) #+# (Matrix.transpose a <> b))++anticommutatorCommutative ::+ (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>+ (General ShapeInt ShapeInt a, General ShapeInt ShapeInt a) -> Bool+anticommutatorCommutative (a,b) =+ approxMatrix 1e-5+ (Symmetric.anticommutator a b)+ (Symmetric.anticommutator b a)++anticommutatorTransposed ::+ (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>+ (General ShapeInt ShapeInt a, General ShapeInt ShapeInt a) -> Bool+anticommutatorTransposed (a,b) =+ approxArray+ (Matrix.forceOrder (ArrMatrix.shapeOrder $ ArrMatrix.shape b) $+ Symmetric.anticommutator (Matrix.transpose a) (Matrix.transpose b))+ (Symmetric.anticommutatorTransposed a b)++++checkForAll ::+ (Show a, QC.Testable test) =>+ Gen.T dim tag a -> (a -> test) -> Tagged tag QC.Property+checkForAll gen = Util.checkForAll (Gen.run gen 3 5)++checkForAllExtra ::+ (Show a, Show b, QC.Testable test) =>+ QC.Gen a -> Gen.T dim tag b ->+ (a -> b -> test) -> Tagged tag QC.Property+checkForAllExtra = Gen.withExtra checkForAll+++testsVar ::+ (Show a, Class.Floating a, Eq a, RealOf a ~ ar, Class.Real ar) =>+ [(String, Tagged a QC.Property)]+testsVar =+ ("stack",+ checkForAll (Gen.stack3 Gen.triangular Gen.matrix Gen.triangular) stack) :+ ("split",+ checkForAll Gen.triangular split) :++ ("gramian",+ checkForAll Gen.matrix gramian) :+ ("gramianTransposed",+ checkForAll Gen.matrix gramianTransposed) :+ ("gramianNonTransposed",+ checkForAll Gen.matrix gramianNonTransposed) :+ ("congruenceDiagonal",+ checkForAll ((,) <$> Gen.vector <-*#> Gen.matrix) congruenceDiagonal) :+ ("congruence",+ checkForAll ((,) <$> Gen.triangular <#*#> Gen.matrix) congruence) :+ ("congruenceDiagonalTransposed",+ checkForAll ((,) <$> Gen.matrix <#*|> Gen.vector)+ congruenceDiagonalTransposed) :+ ("congruenceDiagonalGramian",+ checkForAll Gen.matrix congruenceDiagonalGramian) :+ ("congruenceTransposed",+ checkForAll ((,) <$> Gen.matrix <#*#> Gen.triangular)+ congruenceTransposed) :+ ("congruenceCongruenceDiagonal",+ checkForAllExtra genOrder+ ((,) <$> Gen.vector <-*#> Gen.matrix) congruenceCongruenceDiagonal) :+ ("anticommutator",+ checkForAll ((,) <$> Gen.matrix <#=#> Gen.matrix) anticommutator) :+ ("anticommutatorCommutative",+ checkForAll ((,) <$> Gen.matrix <#=#> Gen.matrix)+ anticommutatorCommutative) :+ ("anticommutatorTransposed",+ checkForAll ((,) <$> Gen.matrix <#=#> Gen.matrix)+ anticommutatorTransposed) :+ []
test/Test/Triangular.hs view
@@ -116,8 +116,7 @@ Triangular.FlexUpper diag ShapeInt a) -> Bool stackUpper (a,b,c) =- let abc = Matrix.fromFull $- Triangular.toSquare $ Triangular.stackUpper a b c+ let abc = expandTriangle $ Triangular.stackUpper a b c in equalArray abc $ (expandTriangle a ||| b !===