lapack 0.3.1 → 0.3.2

raw patch · 45 files changed

+1124/−324 lines, 45 filesdep +doctest-exitcode-stdiodep ~basedep ~comfort-arraydep ~hypernew-component:exe:lapack-debug

Dependencies added: doctest-exitcode-stdio

Dependency ranges changed: base, comfort-array, hyper, random

Files

Changes.md +14/−0
Makefile +8/−2
ReadMe.md +4/−2
bug/GGGLM.hs +19/−0
lapack.cabal +35/−10
src/Numeric/LAPACK/Example/DividedDifference.hs +116/−0
src/Numeric/LAPACK/Example/EconomicAllocation.hs +14/−0
src/Numeric/LAPACK/Linear/Plain.hs +3/−6
src/Numeric/LAPACK/Matrix.hs +13/−1
src/Numeric/LAPACK/Matrix/Array.hs +8/−0
src/Numeric/LAPACK/Matrix/Basic.hs +19/−1
src/Numeric/LAPACK/Matrix/Divide.hs +5/−5
src/Numeric/LAPACK/Matrix/Extent/Kind.hs +0/−35
src/Numeric/LAPACK/Matrix/Extent/Private.hs +127/−119
src/Numeric/LAPACK/Matrix/Hermitian.hs +8/−2
src/Numeric/LAPACK/Matrix/Hermitian/Basic.hs +2/−2
src/Numeric/LAPACK/Matrix/HermitianPositiveDefinite/Linear.hs +3/−6
src/Numeric/LAPACK/Matrix/Multiply.hs +5/−5
src/Numeric/LAPACK/Matrix/Plain/Divide.hs +3/−3
src/Numeric/LAPACK/Matrix/Plain/Multiply.hs +3/−3
src/Numeric/LAPACK/Matrix/Private.hs +1/−1
src/Numeric/LAPACK/Matrix/Shape/Private.hs +8/−0
src/Numeric/LAPACK/Matrix/Square/Basic.hs +1/−7
src/Numeric/LAPACK/Matrix/Square/Linear.hs +3/−5
src/Numeric/LAPACK/Matrix/Symmetric.hs +7/−1
src/Numeric/LAPACK/Matrix/Symmetric/Basic.hs +2/−2
src/Numeric/LAPACK/Matrix/Type.hs +17/−0
src/Numeric/LAPACK/Orthogonal.hs +109/−1
src/Numeric/LAPACK/Orthogonal/Plain.hs +101/−17
src/Numeric/LAPACK/Orthogonal/Private.hs +2/−2
src/Numeric/LAPACK/Permutation/Private.hs +30/−0
src/Numeric/LAPACK/Private.hs +10/−0
src/Numeric/LAPACK/Singular/Plain.hs +2/−5
test-module.list +3/−0
test/DocTest/Main.hs +14/−0
test/DocTest/Numeric/LAPACK/Example/DividedDifference.hs +47/−0
test/DocTest/Numeric/LAPACK/Example/EconomicAllocation.hs +23/−0
test/DocTest/Numeric/LAPACK/Permutation/Private.hs +43/−0
test/Main.hs +30/−6
test/Test/Example.hs +0/−21
test/Test/Generator.hs +18/−0
test/Test/Logic.hs +25/−0
test/Test/Orthogonal.hs +214/−3
test/Test/Permutation.hs +1/−47
test/Test/Shape.hs +4/−4

Changes.md view

@@ -1,5 +1,19 @@ # Change log for the `lapack` package +## 0.3.2++ * `Orthogonal`: `project`, `affineKernelFromSpan`, `affineSpanFromKernel`,+   `leastSquaresConstraint`, `gaussMarkovLinearModel`++ * `Symmetric.fromHermitian`, `Hermitian.fromSymmetric`++ * `instance Monoid Matrix`, especially `mempty`+   for matrices with static shapes.++ * `Extent.Dimensions`: turn from type family to data family++ * Start using `doctest-extract` for simple tests+ ## 0.3.1   * `Matrix.Symmetric`:

Makefile view

@@ -1,10 +1,16 @@ all:	ReadMe.html Changes.html -run-test:+run-test:	update-test 	runhaskell Setup configure --user --enable-tests -fbuildExamples 	runhaskell Setup build 	runhaskell Setup haddock 	./dist/build/lapack-test/lapack-test +run-test-cabal:	update-test+	cabal test --show-details=streaming++update-test:+	doctest-extract -i src/ -o test/ --module-prefix DocTest --library-main=Main $$(cat test-module.list)+ %.html: %.md-	pandoc --output=$@ $<+	pandoc --standalone --metadata pagetitle="$*" --output=$@ $<

ReadMe.md view

@@ -99,7 +99,7 @@ Vectors are `Storable.Array`s from the [comfort-array](https://hackage.haskell.org/package/comfort-array) package. An array can have a fancy shape-like a shape defined by an Haskell enumeration type,+like a shape defined by an enumeration type, the shape of two appended arrays, a shape that is compatible to a Haskell container type, a rectangular or triangular shape.@@ -157,7 +157,7 @@ ~~~~  The relations are defined using two type tags-in order to support matrix transposition flawlessly.+in order to support matrix transposition without hassle. Using `Small Small` for square matrices and `Big Big` for general matrices appears to be arbitrary, but is chosen such that altering `Small` to `Big`@@ -180,6 +180,8 @@ that always have a unit diagonal by construction. The property of a unit diagonal is preserved by some operations and enables some optimizations by LAPACK.+E.g. solving with a unit triangular matrix does not require division+and thus cannot fail due to division by zero. `NonUnit` is a bit of a misnomer. A `NonUnit` matrix can still have a unit diagonal, but in general it has not and no optimizations will take place.

+ bug/GGGLM.hs view

@@ -0,0 +1,19 @@+module Main where++import qualified Numeric.LAPACK.Orthogonal as Ortho+import qualified Numeric.LAPACK.Matrix as Matrix+import qualified Numeric.LAPACK.Vector as Vector+import Numeric.LAPACK.Matrix (shapeInt)++import Data.Array.Comfort.Shape ((:+:)((:+:)))+++main :: IO ()+main = do+   let ab = Matrix.fromList (shapeInt 0) (shapeInt 0 :+: shapeInt 3) ([]::[Double])+   let d = Vector.fromList (shapeInt 0) []+   let xy = Vector.fromList (shapeInt 0 :+: shapeInt 3) [-3,1,2]+   let a = Matrix.tallFromGeneral $ Matrix.takeLeft ab+   let b = Matrix.takeRight ab+   print $ Ortho.gaussMarkovLinearModel a b d+   print $ Ortho.project (Matrix.wideFromGeneral ab) d xy

lapack.cabal view

@@ -1,10 +1,11 @@+Cabal-Version:    2.2 Name:             lapack-Version:          0.3.1-License:          BSD3+Version:          0.3.2+License:          BSD-3-Clause License-File:     LICENSE Author:           Henning Thielemann <haskell@henning-thielemann.de> Maintainer:       Henning Thielemann <haskell@henning-thielemann.de>-Homepage:         http://hub.darcs.net/thielema/lapack/+Homepage:         https://hub.darcs.net/thielema/lapack/ Category:         Data Structures Synopsis:         Numerical Linear Algebra using LAPACK Description:@@ -17,7 +18,8 @@   .   * Based on @comfort-array@:     Allows to precisely express one-column or one-row matrices,-    as well as dense, square, triangular, banded, symmetric and block matrices.+    as well as dense, square, triangular, banded, symmetric, Hermitian,+    banded Hermitian, blocked and LU or QR decomposed matrices.   .   * Support all data types that are supported by LAPACK,     i.e. Float, Double, Complex Float, Complex Double@@ -40,22 +42,26 @@   .   See also: @hmatrix@. Tested-With:      GHC==7.4.2, GHC==7.8.4, GHC==8.2.2-Cabal-Version:    >=1.8 Build-Type:       Simple Extra-Source-Files:   ReadMe.md   Changes.md   Makefile+  test-module.list  Source-Repository this-  Tag:         0.3.1+  Tag:         0.3.2   Type:        darcs-  Location:    http://hub.darcs.net/thielema/lapack/+  Location:    https://hub.darcs.net/thielema/lapack/  Source-Repository head   Type:        darcs-  Location:    http://hub.darcs.net/thielema/lapack/+  Location:    https://hub.darcs.net/thielema/lapack/ +Flag debug+  Description: Compile programs for demonstrating bugs in the LAPACK implementation.+  Default: False+ Library   Build-Depends:     lapack-ffi >=0.0.1 && <0.1,@@ -77,6 +83,7 @@     utility-ht >=0.0.10 && <0.1,     base >=4.5 && <5 +  Default-Language: Haskell98   GHC-Options:      -Wall   Hs-Source-Dirs:   src   Exposed-Modules:@@ -107,6 +114,7 @@     Numeric.LAPACK.Format     Numeric.LAPACK.Output     Numeric.LAPACK.Example.EconomicAllocation+    Numeric.LAPACK.Example.DividedDifference   Other-Modules:     Numeric.LAPACK.Singular.Plain     Numeric.LAPACK.Orthogonal.Plain@@ -136,7 +144,6 @@     Numeric.LAPACK.Matrix.BandedHermitianPositiveDefinite.Linear     Numeric.LAPACK.Matrix.Shape.Private     Numeric.LAPACK.Matrix.Extent.Private-    Numeric.LAPACK.Matrix.Extent.Kind     Numeric.LAPACK.Matrix.Divide     Numeric.LAPACK.Matrix.Multiply     Numeric.LAPACK.Matrix.Class@@ -171,6 +178,7 @@     data-ref >=0.0.1 && <0.1,     unique-logic-tf >=0.5.1 && <0.6,     random >=1.1 && <1.2,+    doctest-exitcode-stdio >=0.0 && <0.1,     quickcheck-transformer >=0.3 && <0.4,     QuickCheck >=2.5 && <3,     ChasingBottoms >=1.2.2 && <1.4,@@ -181,12 +189,16 @@     utility-ht,     base +  Default-Language: Haskell98   GHC-Options:      -Wall   Hs-Source-Dirs:   test   Main-Is:          Main.hs   Other-Modules:+    DocTest.Numeric.LAPACK.Example.EconomicAllocation+    DocTest.Numeric.LAPACK.Example.DividedDifference+    DocTest.Numeric.LAPACK.Permutation.Private+    DocTest.Main     Test.Shape-    Test.Example     Test.Indexed     Test.Divide     Test.Multiply@@ -208,3 +220,16 @@     Test.Logic     Test.Format     Test.Utility++Executable lapack-debug+  Default-Language: Haskell98+  GHC-Options:      -Wall+  Hs-Source-Dirs:   bug+  Main-Is:          GGGLM.hs+  If flag(debug)+    Build-Depends:+      lapack,+      comfort-array,+      base+  Else+    Buildable: False

+ src/Numeric/LAPACK/Example/DividedDifference.hs view

@@ -0,0 +1,116 @@+{- |+This module demonstrates triangular matrices.++It verifies that the divided difference scheme+nicely fits into a triangular matrix,+where function addition is mapped to matrix addition+and function multiplication is mapped to matrix multiplication.++<http://en.wikipedia.org/wiki/Divided_difference>+-}+module Numeric.LAPACK.Example.DividedDifference where++import qualified Numeric.LAPACK.Matrix.Triangular as Triangular+import qualified Numeric.LAPACK.Matrix.Square as Square+import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape+import qualified Numeric.LAPACK.Matrix as Matrix+import qualified Numeric.LAPACK.Vector as Vector+import Numeric.LAPACK.Matrix (ShapeInt, (#+#), (#-#))+import Numeric.LAPACK.Vector (Vector, (|+|), (|-|))+import Numeric.LAPACK.Format ((##))++import qualified Data.Array.Comfort.Shape as Shape+import qualified Data.Array.Comfort.Storable as Array++import qualified Data.List as List+import Data.Semigroup ((<>))+++{- $setup+>>> import qualified Test.Utility as Util+>>> import Test.Utility (approxArray)+>>>+>>> import qualified Numeric.LAPACK.Vector as Vector+>>> import Numeric.LAPACK.Example.DividedDifference (dividedDifferencesMatrix)+>>> import Numeric.LAPACK.Matrix (ShapeInt, (#+#))+>>> import Numeric.LAPACK.Vector ((|+|))+>>>+>>> import qualified Data.Array.Comfort.Storable as Array+>>>+>>> import qualified Test.QuickCheck as QC+>>>+>>> import Control.Monad (liftM2)+>>> import Data.Tuple.HT (mapPair)+>>> import Data.Semigroup ((<>))+>>>+>>> type Vector = Vector.Vector ShapeInt Float+>>>+>>> genDD :: QC.Gen (Vector, (Vector, Vector))+>>> genDD = do+>>>    (ys0,ys1) <-+>>>       fmap (mapPair (Vector.autoFromList, Vector.autoFromList) .+>>>             unzip . take 10) $+>>>       QC.listOf $ liftM2 (,) (Util.genElement 10) (Util.genElement 10)+>>>    xs <- Util.genDistinct 10 10 $ Array.shape ys0+>>>    return (xs,(ys0,ys1))+-}+++size :: Vector ShapeInt a -> Int+size = Shape.zeroBasedSize . Array.shape++subSlices :: Int -> Vector ShapeInt Float -> Vector ShapeInt Float+subSlices k xs = Vector.drop k xs |-| Vector.take (size xs - k) xs++parameterDifferences :: Vector ShapeInt Float -> [Vector ShapeInt Float]+parameterDifferences xs = map (flip subSlices xs) [1 .. size xs - 1]++dividedDifferences ::+   Vector ShapeInt Float -> Vector ShapeInt Float -> [Vector ShapeInt Float]+dividedDifferences xs ys =+   scanl+      (\ddys dxs -> Vector.divide (subSlices 1 ddys) dxs)+      ys+      (parameterDifferences xs)++{- |+prop> QC.forAll genDD $ \(xs, (ys0,ys1)) -> approxArray (dividedDifferencesMatrix xs (ys0|+|ys1)) (dividedDifferencesMatrix xs ys0 #+# dividedDifferencesMatrix xs ys1)+prop> QC.forAll genDD $ \(xs, (ys0,ys1)) -> approxArray (dividedDifferencesMatrix xs (Vector.mul ys0 ys1)) (dividedDifferencesMatrix xs ys0 <> dividedDifferencesMatrix xs ys1)+-}+dividedDifferencesMatrix ::+   Vector ShapeInt Float -> Vector ShapeInt Float ->+   Triangular.Upper ShapeInt Float+dividedDifferencesMatrix xs ys =+   Triangular.upperFromList MatrixShape.RowMajor (Array.shape xs) $+   concat $ List.transpose $ map Vector.toList $ dividedDifferences xs ys+++parameterDifferencesMatrix ::+   Vector ShapeInt Float -> Triangular.Upper ShapeInt Float+parameterDifferencesMatrix xs =+   let ones = Vector.one $ Array.shape xs+       tp = Matrix.tensorProduct MatrixShape.RowMajor+   in Triangular.takeUpper $ Square.fromGeneral $ tp ones xs #-# tp xs ones+++main :: IO ()+main = do+   let xs  = Vector.autoFromList [0,1,4,9,16]+   let ys0 = Vector.autoFromList [3,1,4,1,5]+   let ys1 = Vector.autoFromList [2,7,1,8,1]++   mapM_ (## "%.4f") $ parameterDifferences xs+   parameterDifferencesMatrix xs ## "%.4f"++   let ddys0 = dividedDifferencesMatrix xs ys0+   let ddys1 = dividedDifferencesMatrix xs ys1+   ddys0 ## "%.4f"+   ddys1 ## "%.4f"+   putStrLn ""++   dividedDifferencesMatrix xs (ys0|+|ys1) ## "%.4f"+   ddys0 #+# ddys1 ## "%.4f"+   putStrLn ""++   dividedDifferencesMatrix xs (Vector.mul ys0 ys1) ## "%.4f"+   ddys0 <> ddys1 ## "%.4f"

src/Numeric/LAPACK/Example/EconomicAllocation.hs view

@@ -17,6 +17,17 @@ import Data.Array.Comfort.Shape ((:+:)((:+:)))  +{- $setup+>>> import Numeric.LAPACK.Example.EconomicAllocation+>>> import Test.Utility (approxVector)+>>>+>>> import qualified Numeric.LAPACK.Vector as Vector+>>> import Numeric.LAPACK.Vector ((+++))+>>>+>>> import qualified Data.Array.Comfort.Storable as Array+-}++ type ZeroInt2 = ShapeInt:+:ShapeInt type Vector sh = Vector.Vector sh Double type Matrix height width = Matrix.General height width Double@@ -75,6 +86,9 @@    in x |-| p #*| (k #-# Square.identityFrom k) #\| y  +{- |+prop> let result = iterated expenses0 balances0 in approxVector result $ compensated expenses0 balances0 +++ Vector.zero (Array.shape $ Vector.takeRight result)+-} main :: IO () main = do    iterated expenses0 balances0 ## "%10.2f"

src/Numeric/LAPACK/Linear/Plain.hs view

@@ -55,7 +55,8 @@ import Numeric.LAPACK.Matrix.Private (Full) import Numeric.LAPACK.Linear.Private (solver, withInfo) import Numeric.LAPACK.Vector (Vector)-import Numeric.LAPACK.Private (copyBlock, copyTransposed, copyToColumnMajor)+import Numeric.LAPACK.Private+         (copyBlock, copyTransposed, copyToColumnMajor, copyToColumnMajorTemp)  import qualified Numeric.LAPACK.FFI.Generic as LapackGen import qualified Numeric.BLAS.FFI.Generic as BlasGen@@ -166,11 +167,7 @@             Transposed -> 'T'       aPtr <-          case orderLU of-            RowMajor -> do-               aPtr <- ContT $ withForeignPtr lu-               atmpPtr <- Call.allocaArray (n*n)-               liftIO $ copyToColumnMajor orderLU n n aPtr atmpPtr-               return atmpPtr+            RowMajor -> copyToColumnMajorTemp orderLU n n lu             ColumnMajor -> ContT $ withForeignPtr lu       ldaPtr <- Call.leadingDim lda       ipivPtr <- fmap Perm.deconsElementPtr $ ContT $ withForeignPtr ipiv

src/Numeric/LAPACK/Matrix.hs view

@@ -50,6 +50,7 @@    Basic.OrderBias, leftBias, rightBias, contiguousBias,    (|||), beside,    (===), above,+   stack,     (|*-),    tensorProduct,@@ -605,6 +606,17 @@ contiguousBias = Basic.ContiguousBias  +stack ::+   (Extent.C vert, Extent.C horiz,+    Shape.C heightA, Eq heightA, Shape.C heightB, Eq heightB,+    Shape.C widthA, Eq widthA, Shape.C widthB, Eq widthB, Class.Floating a) =>+   Full vert horiz heightA widthA a -> General heightA widthB a ->+   General heightB widthA a -> Full vert horiz heightB widthB a ->+   Full vert horiz (heightA:+:heightB) (widthA:+:widthB) a+stack = ArrMatrix.lift4 Basic.stack+++ rowSums ::    (Extent.C vert, Extent.C horiz,     Shape.C height, Shape.C width, Class.Floating a) =>@@ -745,7 +757,7 @@ multiplySquare ::    (Multiply.MultiplySquare typ,     Type.HeightOf typ ~ height, Eq height, Shape.C width,-    Extent.C horiz, Extent.C vert, Class.Floating a) =>+    Extent.C vert, Extent.C horiz, Class.Floating a) =>    Mod.Transposition -> Type.Matrix typ a ->    Full vert horiz height width a -> Full vert horiz height width a multiplySquare = Multiply.transposableSquare

src/Numeric/LAPACK/Matrix/Array.hs view

@@ -1,5 +1,6 @@ {-# LANGUAGE TypeFamilies #-} {-# LANGUAGE EmptyDataDecls #-}+{-# LANGUAGE UndecidableInstances #-} module Numeric.LAPACK.Matrix.Array (    Matrix(Array),    ArrayMatrix,@@ -178,6 +179,13 @@  instance (Multiply.MultiplySame sh) => Type.MultiplySame (Array sh) where    multiplySame = lift2 Multiply.same++-- requires UndecidableInstances+instance+   (Plain.SquareShape sh, Box.WidthOf sh ~ width, Shape.Static width) =>+      Type.StaticIdentity (Array sh) where+   staticIdentity =+      lift0 $ Plain.identityOrder MatrixShape.RowMajor Shape.static   zero ::

src/Numeric/LAPACK/Matrix/Basic.hs view

@@ -11,7 +11,7 @@          (Order(RowMajor, ColumnMajor), transposeFromOrder, flipOrder) import Numeric.LAPACK.Matrix.Modifier (Conjugation(NonConjugated)) import Numeric.LAPACK.Matrix.Private-         (Full, Tall, Wide, General, ShapeInt, revealOrder)+         (Full, Tall, Wide, General, fromFull, ShapeInt, revealOrder) import Numeric.LAPACK.Vector (Vector) import Numeric.LAPACK.Scalar (RealOf, zero, one) import Numeric.LAPACK.Shape.Private (Unchecked(Unchecked))@@ -347,6 +347,24 @@    Full Extent.Big horizC (heightA:+:heightB) width a above orderBias appendMode a b =    transpose $ beside orderBias appendMode (transpose a) (transpose b)++stack ::+   (Extent.C vert, Extent.C horiz,+    Shape.C heightA, Eq heightA, Shape.C heightB, Eq heightB,+    Shape.C widthA, Eq widthA, Shape.C widthB, Eq widthB, Class.Floating a) =>+   Full vert horiz heightA widthA a -> General heightA widthB a ->+   General heightB widthA a -> Full vert horiz heightB widthB a ->+   Full vert horiz (heightA:+:heightB) (widthA:+:widthB) a+stack a b c d =+   Array.mapShape+      (\(MatrixShape.Full order _) ->+         MatrixShape.Full order $+         Extent.stack+            (MatrixShape.fullExtent $ Array.shape a)+            (MatrixShape.fullExtent $ Array.shape d)) $+   above RightBias Extent.appendAny+      (beside RightBias Extent.appendAny (fromFull a) b)+      (beside RightBias Extent.appendAny c (fromFull d))   liftRowMajor ::

src/Numeric/LAPACK/Matrix/Divide.hs view

@@ -33,7 +33,7 @@    {-# MINIMAL solve | solveLeft,solveRight #-}    solve ::       (Type.HeightOf typ ~ height, Eq height, Shape.C width,-       Extent.C horiz, Extent.C vert, Class.Floating a) =>+       Extent.C vert, Extent.C horiz, Class.Floating a) =>       Transposition -> Matrix typ a ->       Full vert horiz height width a -> Full vert horiz height width a    solve NonTransposed a b = solveRight a b@@ -41,14 +41,14 @@     solveRight ::       (Type.HeightOf typ ~ height, Eq height, Shape.C width,-       Extent.C horiz, Extent.C vert, Class.Floating a) =>+       Extent.C vert, Extent.C horiz, Class.Floating a) =>       Matrix typ a ->       Full vert horiz height width a -> Full vert horiz height width a    solveRight = solve NonTransposed     solveLeft ::       (Type.WidthOf typ ~ width, Eq width, Shape.C height,-       Extent.C horiz, Extent.C vert, Class.Floating a) =>+       Extent.C vert, Extent.C horiz, Class.Floating a) =>       Full vert horiz height width a ->       Matrix typ a -> Full vert horiz height width a    solveLeft = Basic.swapMultiply $ solve Transposed@@ -61,14 +61,14 @@  (#\##) ::    (Solve typ, Type.HeightOf typ ~ height, Eq height, Shape.C width,-    Extent.C horiz, Extent.C vert, Class.Floating a) =>+    Extent.C vert, Extent.C horiz, Class.Floating a) =>    Matrix typ a ->    Full vert horiz height width a -> Full vert horiz height width a (#\##) = solveRight  (##/#) ::    (Solve typ, Type.WidthOf typ ~ width, Eq width, Shape.C height,-    Extent.C horiz, Extent.C vert, Class.Floating a) =>+    Extent.C vert, Extent.C horiz, Class.Floating a) =>    Full vert horiz height width a ->    Matrix typ a -> Full vert horiz height width a (##/#) = solveLeft

− src/Numeric/LAPACK/Matrix/Extent/Kind.hs

@@ -1,35 +0,0 @@-{-# LANGUAGE TypeFamilies #-}-{-# LANGUAGE GADTs #-}-module Numeric.LAPACK.Matrix.Extent.Kind where---data General height width =-   General {-      generalHeight :: height,-      generalWidth :: width-   } deriving (Eq, Show)--data Tall height width =-   Tall {-      tallHeight :: height,-      tallWidth :: width-   } deriving (Eq, Show)--data Wide height width =-   Wide {-      wideHeight :: height,-      wideWidth :: width-   } deriving (Eq, Show)--data Square height width =-   (height ~ width) =>-   Square {-      squareSize :: height-   }--instance (Eq height, Eq width) => Eq (Square height width) where-   Square a == Square b  =  a==b--instance (Show height, Show width) => Show (Square height width) where-   showsPrec p (Square s) =-      showParen (p>10) (showString "Square " . showsPrec 11 s)

src/Numeric/LAPACK/Matrix/Extent/Private.hs view

@@ -4,7 +4,6 @@ {-# LANGUAGE GADTs #-} module Numeric.LAPACK.Matrix.Extent.Private where -import qualified Numeric.LAPACK.Matrix.Extent.Kind as EK import Numeric.LAPACK.Shape.Private (Unchecked(deconsUnchecked)) import Numeric.LAPACK.Wrapper (Flip(Flip, getFlip)) @@ -13,15 +12,10 @@ import Control.DeepSeq (NFData, rnf)  import Data.Maybe.HT (toMaybe)-import Data.Tuple.HT (swap) import Data.Eq.HT (equating)  -data Extent vertical horizontal height width =-   Extent {-      extentDir :: (vertical,horizontal),-      extentDim :: Dimensions vertical horizontal height width-   }+data family Extent vertical horizontal :: * -> * -> *  instance    (C vertical, C horizontal, NFData height, NFData width) =>@@ -29,10 +23,10 @@    rnf =       getAccessor $       switchTagPair-         (Accessor $ \(Extent o (EK.Square s)) -> rnf (o,s))-         (Accessor $ \(Extent o (EK.Wide h w)) -> rnf (o,(h,w)))-         (Accessor $ \(Extent o (EK.Tall h w)) -> rnf (o,(h,w)))-         (Accessor $ \(Extent o (EK.General h w)) -> rnf (o,(h,w)))+         (Accessor $ \(Square s) -> rnf s)+         (Accessor $ \(Wide h w) -> rnf (h,w))+         (Accessor $ \(Tall h w) -> rnf (h,w))+         (Accessor $ \(General h w) -> rnf (h,w))   data Big = Big deriving (Eq,Show)@@ -47,25 +41,42 @@ type Square sh = Extent Small Small sh sh  -type family Dimensions vertical horizontal :: * -> * -> *+data instance Extent Big Big height width =+   General {+      generalHeight :: height,+      generalWidth :: width+   } -type instance Dimensions Big Big = EK.General-type instance Dimensions Big Small = EK.Tall-type instance Dimensions Small Big = EK.Wide-type instance Dimensions Small Small = EK.Square+data instance Extent Big Small height width =+   Tall {+      tallHeight :: height,+      tallWidth :: width+   } +data instance Extent Small Big height width =+   Wide {+      wideHeight :: height,+      wideWidth :: width+   } +data instance Extent Small Small height width =+   (height ~ width) =>+   Square {+      squareSize :: height+   }++ general :: height -> width -> General height width-general h w = Extent (Big,Big) $ EK.General h w+general = General  tall :: height -> width -> Tall height width-tall h w = Extent (Big,Small) $ EK.Tall h w+tall = Tall  wide :: height -> width -> Wide height width-wide h w = Extent (Small,Big) $ EK.Wide h w+wide = Wide  square :: sh -> Square sh-square sh = Extent (Small,Small) $ EK.Square sh+square = Square   newtype Map vertA horizA vertB horizB height width =@@ -106,10 +117,10 @@ caseTallWide ge =    getCaseTallWide $    switchTagPair-      (CaseTallWide $ \(Extent _ (EK.Square sh)) -> Left $ tall sh sh)+      (CaseTallWide $ \(Square sh) -> Left $ tall sh sh)       (CaseTallWide Right)       (CaseTallWide Left)-      (CaseTallWide $ \(Extent _ (EK.General h w)) ->+      (CaseTallWide $ \(General h w) ->          if ge h w             then Left $ tall h w             else Right $ wide h w)@@ -138,8 +149,8 @@ generalizeTall =    getGenTall $    switchTagPair-      (GenTall id) (GenTall $ \(Extent _ (EK.Square s)) -> wide s s)-      (GenTall id) (GenTall $ \(Extent _ (EK.Tall h w)) -> general h w)+      (GenTall id) (GenTall $ \(Square s) -> wide s s)+      (GenTall id) (GenTall $ \(Tall h w) -> general h w)  newtype GenWide height width vert horiz =    GenWide {@@ -154,8 +165,8 @@    switchTagPair       (GenWide id)       (GenWide id)-      (GenWide $ \(Extent _ (EK.Square s)) -> tall s s)-      (GenWide $ \(Extent _ (EK.Wide h w)) -> general h w)+      (GenWide $ \(Square s) -> tall s s)+      (GenWide $ \(Wide h w) -> general h w)   newtype GenToTall height width vert horiz =@@ -169,8 +180,8 @@ genToTall =    getGenToTall $    switchTagPair-      (GenToTall $ \(Extent _ (EK.Square s)) -> tall s s)-      (GenToTall $ \(Extent _ (EK.Wide h w)) -> general h w)+      (GenToTall $ \(Square s) -> tall s s)+      (GenToTall $ \(Wide h w) -> general h w)       (GenToTall id)       (GenToTall id) @@ -186,16 +197,12 @@ genToWide =    getGenToWide $    switchTagPair-      (GenToWide $ \(Extent _ (EK.Square s)) -> wide s s)+      (GenToWide $ \(Square s) -> wide s s)       (GenToWide id)-      (GenToWide $ \(Extent _ (EK.Tall h w)) -> general h w)+      (GenToWide $ \(Tall h w) -> general h w)       (GenToWide id)  -squareSize :: Square sh -> sh-squareSize (Extent (Small,Small) (EK.Square sh)) = sh-- newtype Accessor a height width vert horiz =    Accessor {getAccessor :: Extent vert horiz height width -> a} @@ -203,19 +210,19 @@ height =    getAccessor $    switchTagPair-      (Accessor (\(Extent _ (EK.Square s)) -> s))-      (Accessor (EK.wideHeight . extentDim))-      (Accessor (EK.tallHeight . extentDim))-      (Accessor (EK.generalHeight . extentDim))+      (Accessor squareSize)+      (Accessor wideHeight)+      (Accessor tallHeight)+      (Accessor generalHeight)  width :: (C vert, C horiz) => Extent vert horiz height width -> width width =    getAccessor $    switchTagPair-      (Accessor (\(Extent _ (EK.Square s)) -> s))-      (Accessor (EK.wideWidth . extentDim))-      (Accessor (EK.tallWidth . extentDim))-      (Accessor (EK.generalWidth . extentDim))+      (Accessor (\(Square s) -> s))+      (Accessor wideWidth)+      (Accessor tallWidth)+      (Accessor generalWidth)   dimensions ::@@ -232,7 +239,7 @@  fromSquareLiberal :: (C vert, C horiz) =>    Extent Small Small height width -> Extent vert horiz height width-fromSquareLiberal x@(Extent _ (EK.Square _)) = genSquare $ height x+fromSquareLiberal (Square s) = genSquare s  squareFromGeneral ::    (C vert, C horiz, Eq size) =>@@ -258,10 +265,10 @@ transpose =    getTranspose $    switchTagPair-      (Transpose $ \(Extent o (EK.Square s)) -> Extent o (EK.Square s))-      (Transpose $ \(Extent o (EK.Wide h w)) -> Extent (swap o) (EK.Tall w h))-      (Transpose $ \(Extent o (EK.Tall h w)) -> Extent (swap o) (EK.Wide w h))-      (Transpose $ \(Extent o (EK.General h w)) -> Extent o (EK.General w h))+      (Transpose $ \(Square s) -> Square s)+      (Transpose $ \(Wide h w) -> Tall w h)+      (Transpose $ \(Tall h w) -> Wide w h)+      (Transpose $ \(General h w) -> General w h)   newtype Equal height width vert horiz =@@ -277,10 +284,11 @@    (==) =       getEqual $       switchTagPair-         (Equal $ equating extentDim)-         (Equal $ equating extentDim)-         (Equal $ equating extentDim)-         (Equal $ equating extentDim)+         (Equal $ \(Square a) (Square b) -> a==b)+         (Equal $ \a b -> equating wideHeight a b && equating wideWidth a b)+         (Equal $ \a b -> equating tallHeight a b && equating tallWidth a b)+         (Equal $ \a b ->+            equating generalHeight a b && equating generalWidth a b)   instance@@ -324,8 +332,8 @@ widen w =    getWiden $    switchTag-      (Widen (\(Extent o x) -> Extent o (x{EK.wideWidth = w})))-      (Widen (\(Extent o x) -> Extent o (x{EK.generalWidth = w})))+      (Widen (\x -> x{wideWidth = w}))+      (Widen (\x -> x{generalWidth = w}))  reduceWideHeight ::    (C vert) =>@@ -333,8 +341,8 @@ reduceWideHeight h =    getWiden $    switchTag-      (Widen (\(Extent o x) -> Extent o (x{EK.wideHeight = h})))-      (Widen (\(Extent o x) -> Extent o (x{EK.generalHeight = h})))+      (Widen (\x -> x{wideHeight = h}))+      (Widen (\x -> x{generalHeight = h}))   newtype Adapt heightA widthA heightB widthB vert horiz =@@ -351,10 +359,10 @@ reduceConsistent h w =    getAdapt $    switchTagPair-      (Adapt $ \(Extent o (EK.Square _)) -> Extent o (EK.Square h))-      (Adapt $ \(Extent o (EK.Wide _ _)) -> Extent o (EK.Wide h w))-      (Adapt $ \(Extent o (EK.Tall _ _)) -> Extent o (EK.Tall h w))-      (Adapt $ \(Extent o (EK.General _ _)) -> Extent o (EK.General h w))+      (Adapt $ \(Square _) -> Square h)+      (Adapt $ \(Wide _ _) -> Wide h w)+      (Adapt $ \(Tall _ _) -> Tall h w)+      (Adapt $ \(General _ _) -> General h w)   class (C vert, C horiz) => GeneralTallWide vert horiz where@@ -371,9 +379,9 @@ mapHeight f =    getAdapt $    switchTagGTW-      (Adapt $ \(Extent o (EK.Wide h w)) -> Extent o (EK.Wide (f h) w))-      (Adapt $ \(Extent o (EK.Tall h w)) -> Extent o (EK.Tall (f h) w))-      (Adapt $ \(Extent o (EK.General h w)) -> Extent o (EK.General (f h) w))+      (Adapt $ \(Wide h w) -> Wide (f h) w)+      (Adapt $ \(Tall h w) -> Tall (f h) w)+      (Adapt $ \(General h w) -> General (f h) w)  mapWidth ::    (GeneralTallWide vert horiz) =>@@ -382,12 +390,12 @@ mapWidth f =    getAdapt $    switchTagGTW-      (Adapt $ \(Extent o (EK.Wide h w)) -> Extent o (EK.Wide h (f w)))-      (Adapt $ \(Extent o (EK.Tall h w)) -> Extent o (EK.Tall h (f w)))-      (Adapt $ \(Extent o (EK.General h w)) -> Extent o (EK.General h (f w)))+      (Adapt $ \(Wide h w) -> Wide h (f w))+      (Adapt $ \(Tall h w) -> Tall h (f w))+      (Adapt $ \(General h w) -> General h (f w))  mapSquareSize :: (shA -> shB) -> Square shA -> Square shB-mapSquareSize f (Extent o (EK.Square s)) = Extent o (EK.Square (f s))+mapSquareSize f (Square s) = Square (f s)   mapWrap ::@@ -399,11 +407,10 @@ mapWrap fh fw =    getAdapt $    switchTagPair-      (Adapt $ \(Extent o (EK.Square h)) -> Extent o (EK.Square (fh h)))-      (Adapt $ \(Extent o (EK.Wide h w)) -> Extent o (EK.Wide (fh h) (fw w)))-      (Adapt $ \(Extent o (EK.Tall h w)) -> Extent o (EK.Tall (fh h) (fw w)))-      (Adapt $ \(Extent o (EK.General h w)) ->-                  Extent o (EK.General (fh h) (fw w)))+      (Adapt $ \(Square h) -> Square (fh h))+      (Adapt $ \(Wide h w) -> Wide (fh h) (fw w))+      (Adapt $ \(Tall h w) -> Tall (fh h) (fw w))+      (Adapt $ \(General h w) -> General (fh h) (fw w))  {- only admissible since GHC-7.8 mapUnwrap ::@@ -415,11 +422,10 @@ mapUnwrap fh fw =    getAdapt $    switchTagPair-      (Adapt $ \(Extent o (EK.Square h)) -> Extent o (EK.Square (fh h)))-      (Adapt $ \(Extent o (EK.Wide h w)) -> Extent o (EK.Wide (fh h) (fw w)))-      (Adapt $ \(Extent o (EK.Tall h w)) -> Extent o (EK.Tall (fh h) (fw w)))-      (Adapt $ \(Extent o (EK.General h w)) ->-                  Extent o (EK.General (fh h) (fw w)))+      (Adapt $ \(Square h) -> Square (fh h))+      (Adapt $ \(Wide h w) -> Wide (fh h) (fw w))+      (Adapt $ \(Tall h w) -> Tall (fh h) (fw w))+      (Adapt $ \(General h w) -> General (fh h) (fw w)) -}  recheck ::@@ -429,14 +435,11 @@ recheck =    getAdapt $    switchTagPair-      (Adapt $ \(Extent o (EK.Square h)) ->-         Extent o (EK.Square (deconsUnchecked h)))-      (Adapt $ \(Extent o (EK.Wide h w)) ->-         Extent o (EK.Wide (deconsUnchecked h) (deconsUnchecked w)))-      (Adapt $ \(Extent o (EK.Tall h w)) ->-         Extent o (EK.Tall (deconsUnchecked h) (deconsUnchecked w)))-      (Adapt $ \(Extent o (EK.General h w)) ->-         Extent o (EK.General (deconsUnchecked h) (deconsUnchecked w)))+      (Adapt $ \(Square h) -> Square (deconsUnchecked h))+      (Adapt $ \(Wide h w) -> Wide (deconsUnchecked h) (deconsUnchecked w))+      (Adapt $ \(Tall h w) -> Tall (deconsUnchecked h) (deconsUnchecked w))+      (Adapt $ \(General h w) ->+                              General (deconsUnchecked h) (deconsUnchecked w))   @@ -456,48 +459,18 @@ fuse =    getFuse $    switchTagPair-      (Fuse $-       \(Extent o (EK.Square s0)) (Extent _ (EK.Square s1)) ->-         toMaybe (s0==s1) $ Extent o (EK.Square s0))-      (Fuse $-       \(Extent o (EK.Wide h f0)) (Extent _ (EK.Wide f1 w)) ->-         toMaybe (f0==f1) $ Extent o (EK.Wide h w))-      (Fuse $-       \(Extent o (EK.Tall h f0)) (Extent _ (EK.Tall f1 w)) ->-         toMaybe (f0==f1) $ Extent o (EK.Tall h w))-      (Fuse $-       \(Extent o (EK.General h f0)) (Extent _ (EK.General f1 w)) ->-         toMaybe (f0==f1) $ Extent o (EK.General h w))+      (Fuse $ \(Square s0) (Square s1) -> toMaybe (s0==s1) $ Square s0)+      (Fuse $ \(Wide h f0) (Wide f1 w) -> toMaybe (f0==f1) $ Wide h w)+      (Fuse $ \(Tall h f0) (Tall f1 w) -> toMaybe (f0==f1) $ Tall h w)+      (Fuse $ \(General h f0) (General f1 w) -> toMaybe (f0==f1) $ General h w)  -newtype Kronecker heightA widthA heightB widthB vert horiz =-   Kronecker {-      getKronecker ::-         Extent vert horiz heightA widthA ->-         Extent vert horiz heightB widthB ->-         Extent vert horiz (heightA,heightB) (widthA,widthB)-   }- kronecker ::    (C vert, C horiz) =>    Extent vert horiz heightA widthA ->    Extent vert horiz heightB widthB ->    Extent vert horiz (heightA,heightB) (widthA,widthB)-kronecker =-   getKronecker $-   switchTagPair-      (Kronecker $-       \(Extent o (EK.Square s0)) (Extent _ (EK.Square s1)) ->-         Extent o (EK.Square (s0,s1)))-      (Kronecker $-       \(Extent o (EK.Wide h0 w0)) (Extent _ (EK.Wide h1 w1)) ->-         Extent o (EK.Wide (h0,h1) (w0,w1)))-      (Kronecker $-       \(Extent o (EK.Tall h0 w0)) (Extent _ (EK.Tall h1 w1)) ->-         Extent o (EK.Tall (h0,h1) (w0,w1)))-      (Kronecker $-       \(Extent o (EK.General h0 w0)) (Extent _ (EK.General h1 w1)) ->-         Extent o (EK.General (h0,h1) (w0,w1)))+kronecker = stackGen (,) (,)   @@ -561,6 +534,41 @@    AppendMode vertA vertB (Append vertA vertB) height widthA widthB appendAny =    getAppendAny $ switchTag (AppendAny appendLeftAux) (AppendAny appendRight)+++stack ::+   (C vert, C horiz) =>+   Extent vert horiz heightA widthA ->+   Extent vert horiz heightB widthB ->+   Extent vert horiz (heightA:+:heightB) (widthA:+:widthB)+stack = stackGen (:+:) (:+:)++newtype Stack f heightA widthA heightB widthB vert horiz =+   Stack {+      getStack ::+         Extent vert horiz heightA widthA ->+         Extent vert horiz heightB widthB ->+         Extent vert horiz (f heightA heightB) (f widthA widthB)+   }++stackGen ::+   (C vert, C horiz) =>+   (heightA -> heightB -> f heightA heightB) ->+   (widthA -> widthB -> f widthA widthB) ->+   Extent vert horiz heightA widthA ->+   Extent vert horiz heightB widthB ->+   Extent vert horiz (f heightA heightB) (f widthA widthB)+stackGen fh fw =+   getStack $+   switchTagPair+      (Stack $ \(Square sa) (Square sb) ->+         Square (fh sa sb))+      (Stack $ \(Wide ha wa) (Wide hb wb) ->+         Wide (fh ha hb) (fw wa wb))+      (Stack $ \(Tall ha wa) (Tall hb wb) ->+         Tall (fh ha hb) (fw wa wb))+      (Stack $ \(General ha wa) (General hb wb) ->+         General (fh ha hb) (fw wa wb))

src/Numeric/LAPACK/Matrix/Hermitian.hs view

@@ -25,7 +25,7 @@    sumRank1, sumRank1NonEmpty,    sumRank2, sumRank2NonEmpty, -   toSquare,+   toSquare, fromSymmetric,    gramian,            gramianAdjoint,    congruenceDiagonal, congruenceDiagonalAdjoint,    congruence,         congruenceAdjoint,@@ -46,7 +46,7 @@ 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 (Hermitian)+import Numeric.LAPACK.Matrix.Array.Triangular (Hermitian, Symmetric) import Numeric.LAPACK.Matrix.Array (Full, General, Square) import Numeric.LAPACK.Matrix.Shape.Private (Order) import Numeric.LAPACK.Matrix.Modifier (Transposition(NonTransposed, Transposed))@@ -186,9 +186,15 @@ sumRank2NonEmpty order (NonEmpty.Cons xy xys) =    sumRank2 order (Array.shape $ fst $ snd xy) (xy:xys) + toSquare ::    (Shape.C sh, Class.Floating a) => Hermitian sh a -> Square sh a toSquare = ArrMatrix.lift1 Basic.toSquare++fromSymmetric :: (Shape.C sh, Class.Real a) => Symmetric sh a -> Hermitian sh a+fromSymmetric =+   ArrMatrix.lift1 $ Array.mapShape MatrixShape.hermitianFromSymmetric+  {- | gramian A = A^H * A

src/Numeric/LAPACK/Matrix/Hermitian/Basic.hs view

@@ -428,7 +428,7 @@    \bPtr -> gramianIO order a bPtr $ gramianParameters order extent  gramianParameters ::-   (Extent.C horiz, Extent.C vert, Shape.C height, Shape.C width) =>+   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width) =>    Order ->    Extent.Extent vert horiz height width ->    ((Int, Int), (Char, Char, Int))@@ -450,7 +450,7 @@    \bPtr -> gramianIO order a bPtr $ gramianAdjointParameters order extent  gramianAdjointParameters ::-   (Extent.C horiz, Extent.C vert, Shape.C height, Shape.C width) =>+   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width) =>    Order ->    Extent.Extent vert horiz height width ->    ((Int, Int), (Char, Char, Int))

src/Numeric/LAPACK/Matrix/HermitianPositiveDefinite/Linear.hs view

@@ -61,8 +61,7 @@             LapackGen.pptrs uploPtr nPtr nrhsPtr apPtr xPtr ldxPtr  -inverse ::-   (Shape.C sh, Class.Floating a) => Hermitian sh a -> Hermitian sh a+inverse :: (Shape.C sh, Class.Floating a) => Hermitian sh a -> Hermitian sh a inverse    (Array shape@(MatrixShape.Hermitian order sh) a) =       Array.unsafeCreateWithSize shape $ \triSize bPtr -> do@@ -75,8 +74,7 @@          withInfo definiteMsg "pptrf" $ LapackGen.pptrf uploPtr nPtr bPtr          withInfo rankMsg "pptri" $ LapackGen.pptri uploPtr nPtr bPtr -decompose ::-   (Shape.C sh, Class.Floating a) => Hermitian sh a -> Upper sh a+decompose :: (Shape.C sh, Class.Floating a) => Hermitian sh a -> Upper sh a decompose    (Array (MatrixShape.Hermitian order sh) a) =       Array.unsafeCreateWithSize@@ -91,8 +89,7 @@          withInfo definiteMsg "pptrf" $ LapackGen.pptrf uploPtr nPtr bPtr  -determinant ::-   (Shape.C sh, Class.Floating a) => Hermitian sh a -> RealOf a+determinant :: (Shape.C sh, Class.Floating a) => Hermitian sh a -> RealOf a determinant =    getDeterminant $    Class.switchFloating

src/Numeric/LAPACK/Matrix/Multiply.hs view

@@ -75,7 +75,7 @@    {-# MINIMAL transposableSquare | fullSquare,squareFull #-}    transposableSquare ::       (Type.HeightOf typ ~ height, Eq height, Shape.C width,-       Extent.C horiz, Extent.C vert, Class.Floating a) =>+       Extent.C vert, Extent.C horiz, Class.Floating a) =>       Transposition -> Matrix typ a ->       Full vert horiz height width a -> Full vert horiz height width a    transposableSquare NonTransposed a b = squareFull a b@@ -84,14 +84,14 @@     squareFull ::       (Type.HeightOf typ ~ height, Eq height, Shape.C width,-       Extent.C horiz, Extent.C vert, Class.Floating a) =>+       Extent.C vert, Extent.C horiz, Class.Floating a) =>       Matrix typ a ->       Full vert horiz height width a -> Full vert horiz height width a    squareFull = transposableSquare NonTransposed     fullSquare ::       (Type.WidthOf typ ~ width, Eq width, Shape.C height,-       Extent.C horiz, Extent.C vert, Class.Floating a) =>+       Extent.C vert, Extent.C horiz, Class.Floating a) =>       Full vert horiz height width a ->       Matrix typ a -> Full vert horiz height width a    fullSquare b a =@@ -102,14 +102,14 @@  (#*##) ::    (MultiplySquare typ, Type.HeightOf typ ~ height, Eq height, Shape.C width,-    Extent.C horiz, Extent.C vert, Class.Floating a) =>+    Extent.C vert, Extent.C horiz, Class.Floating a) =>    Matrix typ a ->    Full vert horiz height width a -> Full vert horiz height width a (#*##) = squareFull  (##*#) ::    (MultiplySquare typ, Type.WidthOf typ ~ width, Eq width, Shape.C height,-    Extent.C horiz, Extent.C vert, Class.Floating a) =>+    Extent.C vert, Extent.C horiz, Class.Floating a) =>    Full vert horiz height width a ->    Matrix typ a -> Full vert horiz height width a (##*#) = fullSquare

src/Numeric/LAPACK/Matrix/Plain/Divide.hs view

@@ -50,7 +50,7 @@    {-# MINIMAL solve | solveLeft,solveRight #-}    solve ::       (Class.Floating a, Box.HeightOf shape ~ height, Eq height,-       Extent.C horiz, Extent.C vert, Shape.C width) =>+       Extent.C vert, Extent.C horiz, Shape.C width) =>       Transposition -> Array shape a ->       Full vert horiz height width a -> Full vert horiz height width a    solve NonTransposed a b = solveRight a b@@ -58,14 +58,14 @@     solveRight ::       (Class.Floating a, Box.HeightOf shape ~ height, Eq height,-       Extent.C horiz, Extent.C vert, Shape.C width) =>+       Extent.C vert, Extent.C horiz, Shape.C width) =>       Array shape a ->       Full vert horiz height width a -> Full vert horiz height width a    solveRight = solve NonTransposed     solveLeft ::       (Class.Floating a, Box.HeightOf shape ~ width, Eq width,-       Extent.C horiz, Extent.C vert, Shape.C height) =>+       Extent.C vert, Extent.C horiz, Shape.C height) =>       Full vert horiz height width a ->       Array shape a ->       Full vert horiz height width a

src/Numeric/LAPACK/Matrix/Plain/Multiply.hs view

@@ -104,7 +104,7 @@    {-# MINIMAL transposableSquare | fullSquare,squareFull #-}    transposableSquare ::       (Box.HeightOf shape ~ height, Eq height, Shape.C width,-       Extent.C horiz, Extent.C vert, Class.Floating a) =>+       Extent.C vert, Extent.C horiz, Class.Floating a) =>       Transposition -> Array shape a ->       Full vert horiz height width a -> Full vert horiz height width a    transposableSquare NonTransposed a b = squareFull a b@@ -112,14 +112,14 @@     squareFull ::       (Box.HeightOf shape ~ height, Eq height, Shape.C width,-       Extent.C horiz, Extent.C vert, Class.Floating a) =>+       Extent.C vert, Extent.C horiz, Class.Floating a) =>       Array shape a ->       Full vert horiz height width a -> Full vert horiz height width a    squareFull = transposableSquare NonTransposed     fullSquare ::       (Box.WidthOf shape ~ width, Eq width, Shape.C height,-       Extent.C horiz, Extent.C vert, Class.Floating a) =>+       Extent.C vert, Extent.C horiz, Class.Floating a) =>       Full vert horiz height width a ->       Array shape a -> Full vert horiz height width a    fullSquare = swapMultiply $ transposableSquare Transposed

src/Numeric/LAPACK/Matrix/Private.hs view

@@ -69,7 +69,7 @@   revealOrder ::-   (Extent.C horiz, Extent.C vert) =>+   (Extent.C vert, Extent.C horiz) =>    Full vert horiz height width a ->    Either (Array (height,width) a) (Array (width,height) a) revealOrder (Array (MatrixShape.Full order extent) a) =

src/Numeric/LAPACK/Matrix/Shape/Private.hs view

@@ -363,6 +363,10 @@    height = hermitianSize    width = hermitianSize +hermitianFromSymmetric :: Symmetric size -> Hermitian size+hermitianFromSymmetric (Triangular _diag _uplo order size) =+   Hermitian order size+ uploFromOrder :: Order -> Char uploFromOrder RowMajor = 'L' uploFromOrder ColumnMajor = 'U'@@ -413,6 +417,10 @@    type WidthOf (Triangular lo diag up size) = size    height = triangularSize    width = triangularSize++symmetricFromHermitian :: Hermitian size -> Symmetric size+symmetricFromHermitian (Hermitian order size) =+   Triangular NonUnit autoUplo order size   data Unit = Unit deriving (Eq, Show)

src/Numeric/LAPACK/Matrix/Square/Basic.hs view

@@ -196,13 +196,7 @@    Square sizeA a -> Full vert horiz sizeA sizeB a ->    Full horiz vert sizeB sizeA a -> Square sizeB a ->    Square (sizeA:+:sizeB) a-stack a b c d =-   fromGeneral $-   Basic.above Basic.RightBias Extent.appendAny-      (Basic.beside Basic.RightBias Extent.appendAny-         (Matrix.fromFull a) (Matrix.fromFull b))-      (Basic.beside Basic.RightBias Extent.appendAny-         (Matrix.fromFull c) (Matrix.fromFull d))+stack a b c d = Basic.stack a (Matrix.fromFull b) (Matrix.fromFull c) d   square :: (Shape.C sh, Class.Floating a) => Square sh a -> Square sh a

src/Numeric/LAPACK/Matrix/Square/Linear.hs view

@@ -14,7 +14,7 @@ import Numeric.LAPACK.Matrix.Shape.Private (transposeFromOrder) import Numeric.LAPACK.Matrix.Private (Full, Square, argSquare) import Numeric.LAPACK.Private-         (withAutoWorkspaceInfo, copyBlock, copyToTemp, copyToColumnMajor)+         (withAutoWorkspaceInfo, copyBlock, copyToTemp, copyToColumnMajorTemp)  import qualified Numeric.LAPACK.FFI.Generic as LapackGen import qualified Numeric.Netlib.Utility as Call@@ -55,14 +55,12 @@ _solve =    argSquare $ \orderA shA a ->    solver "Square.solve" shA $ \n nPtr nrhsPtr xPtr ldxPtr -> do-      aPtr <- ContT $ withForeignPtr a-      atmpPtr <- Call.allocaArray (n*n)+      aPtr <- copyToColumnMajorTemp orderA n n a       ldaPtr <- Call.leadingDim n       ipivPtr <- Call.allocaArray n       liftIO $ do-         copyToColumnMajor orderA n n aPtr atmpPtr          withInfo "gesv" $-            LapackGen.gesv nPtr nrhsPtr atmpPtr ldaPtr ipivPtr xPtr ldxPtr+            LapackGen.gesv nPtr nrhsPtr aPtr ldaPtr ipivPtr xPtr ldxPtr   inverse :: (Shape.C sh, Class.Floating a) => Square sh a -> Square sh a

src/Numeric/LAPACK/Matrix/Symmetric.hs view

@@ -13,6 +13,7 @@    split,     toSquare,+   fromHermitian,     gramian,            gramianTransposed,    congruenceDiagonal, congruenceDiagonalTransposed,@@ -26,7 +27,7 @@ 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.Triangular (Symmetric, Hermitian) import Numeric.LAPACK.Matrix.Array (Full, General, Square) import Numeric.LAPACK.Matrix.Shape.Private (Order) import Numeric.LAPACK.Matrix.Private (ShapeInt)@@ -35,6 +36,7 @@  import qualified Numeric.Netlib.Class as Class +import qualified Data.Array.Comfort.Storable as Array import qualified Data.Array.Comfort.Shape as Shape import Data.Array.Comfort.Shape ((:+:)) @@ -60,6 +62,10 @@  toSquare :: (Shape.C sh, Class.Floating a) => Symmetric sh a -> Square sh a toSquare = Triangular.toSquare++fromHermitian :: (Shape.C sh, Class.Real a) => Hermitian sh a -> Symmetric sh a+fromHermitian =+   ArrMatrix.lift1 $ Array.mapShape MatrixShape.symmetricFromHermitian   identity :: (Shape.C sh, Class.Floating a) => Order -> sh -> Symmetric sh a

src/Numeric/LAPACK/Matrix/Symmetric/Basic.hs view

@@ -51,7 +51,7 @@    \bPtr -> gramianIO order a bPtr $ gramianParameters order extent  gramianParameters ::-   (Extent.C horiz, Extent.C vert, Shape.C height, Shape.C width) =>+   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width) =>    Order ->    Extent.Extent vert horiz height width ->    ((Int, Int), (Char, Char, Int))@@ -73,7 +73,7 @@    \bPtr -> gramianIO order a bPtr $ gramianTransposedParameters order extent  gramianTransposedParameters ::-   (Extent.C horiz, Extent.C vert, Shape.C height, Shape.C width) =>+   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width) =>    Order ->    Extent.Extent vert horiz height width ->    ((Int, Int), (Char, Char, Int))

src/Numeric/LAPACK/Matrix/Type.hs view

@@ -16,6 +16,7 @@  import qualified Data.Array.Comfort.Shape as Shape +import Data.Monoid (Monoid, mempty, mappend) import Data.Semigroup (Semigroup, (<>))  @@ -80,6 +81,22 @@ instance (Shape.C sh, Eq sh) => MultiplySame (Perm.Permutation sh) where    multiplySame (Permutation a) (Permutation b) =       Permutation $ Perm.multiply b a+++instance+   (MultiplySame typ, StaticIdentity typ, Class.Floating a) =>+      Monoid (Matrix typ a) where+   mappend = (<>)+   mempty = staticIdentity++class StaticIdentity typ where+   staticIdentity :: (Class.Floating a) => Matrix typ a++instance (Shape.Static shape) => StaticIdentity (Scale shape) where+   staticIdentity = Scale Shape.static 1++instance (Shape.Static sh) => StaticIdentity (Perm.Permutation sh) where+   staticIdentity = Permutation $ Perm.identity Shape.static   scaleWithCheck :: (Eq shape) =>

src/Numeric/LAPACK/Orthogonal.hs view

@@ -5,9 +5,15 @@    leastSquaresMinimumNormRCond,    pseudoInverseRCond, +   project,+   leastSquaresConstraint,+   gaussMarkovLinearModel,+    determinant,    determinantAbsolute,    complement,+   affineSpanFromKernel,+   affineKernelFromSpan,     householder,    householderTall,@@ -16,12 +22,17 @@ import qualified Numeric.LAPACK.Orthogonal.Householder as HH import qualified Numeric.LAPACK.Orthogonal.Plain as Plain +import qualified Numeric.LAPACK.Matrix.Multiply as Multiply import qualified Numeric.LAPACK.Matrix.Triangular as Triangular+import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix import qualified Numeric.LAPACK.Matrix.Basic as Basic+import qualified Numeric.LAPACK.Matrix.Type as Matrix import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent-import Numeric.LAPACK.Matrix.Array (Full, Tall, Square)+import qualified Numeric.LAPACK.Vector as Vector+import Numeric.LAPACK.Matrix.Array (Full, General, Tall, Wide, Square) import Numeric.LAPACK.Matrix.Private (ShapeInt)+import Numeric.LAPACK.Vector (Vector, (|-|)) import Numeric.LAPACK.Scalar (RealOf)  import qualified Numeric.Netlib.Class as Class@@ -94,6 +105,65 @@    Basic.uncheck . ArrMatrix.toVector  +{- |+@project b d x@ projects @x@ on the plane described by @B*x = d@.++@b@ must have full rank.+-}+project ::+   (Shape.C height, Eq height, Shape.C width, Eq width, Class.Floating a) =>+   Wide height width a -> Vector height a ->+   Vector width a -> Vector width a+project b d x =+   x+   |-|+   ArrMatrix.unliftColumn MatrixShape.ColumnMajor+      (minimumNorm b) (Multiply.matrixVector b x |-| d)+++{- |+@leastSquaresConstraint a c b d@ computes @x@+with minimal @|| c - A*x ||_2@ and constraint @B*x = d@.++@b@ must be wide and @a===b@ must be tall+and both matrices must have full rank.+-}+leastSquaresConstraint ::+   (Shape.C height, Eq height,+    Shape.C width, Eq width,+    Shape.C constraints, Eq constraints, Class.Floating a) =>+   General height width a -> Vector height a ->+   Wide constraints width a -> Vector constraints a ->+   Vector width a+leastSquaresConstraint a c b d =+   Plain.leastSquaresConstraint+      (ArrMatrix.toVector a) c+      (ArrMatrix.toVector b) d++{- |+@gaussMarkovLinearModel a b d@ computes @(x,y)@+with minimal @|| y ||_2@ and constraint @d = A*x + B*y@.++@a@ must be tall and @a|||b@ must be wide+and both matrices must have full rank.+-}+gaussMarkovLinearModel ::+   (Shape.C height, Eq height,+    Shape.C width, Eq width,+    Shape.C opt, Eq opt, Class.Floating a) =>+   Tall height width a -> General height opt a -> Vector height a ->+   (Vector width a, Vector opt a)+gaussMarkovLinearModel a b d =+   {-+   Fortran-LAPACK and OpenBLAS would leave Y uninitalized+   instead of setting Y to zeros.+   -}+   if Shape.size (Matrix.height a) == 0+      then (Vector.zero (Matrix.width a), Vector.zero (Matrix.width b))+      else Plain.gaussMarkovLinearModel+               (ArrMatrix.toVector a) (ArrMatrix.toVector b) d++ {- @(q,r) = householder a@ means that @q@ is unitary and @r@ is upper triangular and @a = multiply q r@.@@ -145,3 +215,41 @@    (Shape.C height, Shape.C width, Class.Floating a) =>    Tall height width a -> Tall height ShapeInt a complement = ArrMatrix.lift1 Plain.complement+++{- |+> affineSpanFromKernel a b == (c,d)++Means:+An affine subspace is given implicitly by {x : a#*|x == b}.+Convert this into an explicit representation {c#*|y|+|d : y}.+Matrix @a@ must have full rank,+otherwise the explicit representation will miss dimensions+and we cannot easily determine the origin @d@ as a minimum norm solution.+-}+affineSpanFromKernel ::+   (Shape.C width, Eq width, Shape.C height, Eq height, Class.Floating a) =>+   Wide height width a -> Vector height a ->+   (Tall width ShapeInt a, Vector width a)+affineSpanFromKernel a b =+   let qr = HH.fromMatrix $ ArrMatrix.lift1 Basic.adjoint a+   in (ArrMatrix.lift0 $ Plain.extractComplement qr,+       ArrMatrix.unliftColumn MatrixShape.ColumnMajor (HH.minimumNorm qr) b)++{- |+This conversion is somehow inverse to 'affineSpanFromKernel'.+However, it is not precisely inverse in either direction.+This is because both 'affineSpanFromKernel' and 'affineKernelFromSpan'+accept non-orthogonal matrices but always return orthogonal ones.++In @affineKernelFromSpan c d@,+matrix @c@ should have full rank,+otherwise the implicit representation will miss dimensions.+-}+affineKernelFromSpan ::+   (Shape.C width, Eq width, Shape.C height, Eq height, Class.Floating a) =>+   Tall height width a -> Vector height a ->+   (Wide ShapeInt height a, Vector ShapeInt a)+affineKernelFromSpan c d =+   let a = Basic.adjoint $ Plain.complement $ ArrMatrix.toVector c+   in (ArrMatrix.lift0 a, Basic.multiplyVector a d)

src/Numeric/LAPACK/Orthogonal/Plain.hs view

@@ -5,8 +5,12 @@    leastSquaresMinimumNormRCond,    pseudoInverseRCond, +   leastSquaresConstraint,+   gaussMarkovLinearModel,+    determinantAbsolute,    complement,+   extractComplement,    ) where  import qualified Numeric.LAPACK.Orthogonal.Private as HH@@ -17,11 +21,14 @@ import qualified Numeric.LAPACK.Matrix.Basic as Basic import qualified Numeric.LAPACK.Vector as Vector import Numeric.LAPACK.Matrix.Shape.Private (Order(RowMajor,ColumnMajor))-import Numeric.LAPACK.Matrix.Private (Full, Tall, ShapeInt, shapeInt)+import Numeric.LAPACK.Matrix.Private+         (Full, General, Tall, Wide, ShapeInt, shapeInt)+import Numeric.LAPACK.Vector (Vector) import Numeric.LAPACK.Scalar (RealOf, zero, absolute) import Numeric.LAPACK.Private          (lacgv, peekCInt,-          copySubMatrix, copyToTemp, copyToColumnMajor, copyToSubColumnMajor,+          copySubMatrix, copyToTemp,+          copyToColumnMajorTemp, copyToSubColumnMajor,           withAutoWorkspaceInfo, rankMsg, errorCodeMsg, createHigherArray)  import qualified Numeric.LAPACK.FFI.Generic as LapackGen@@ -30,6 +37,7 @@ import qualified Numeric.Netlib.Utility as Call import qualified Numeric.Netlib.Class as Class +import qualified Data.Array.Comfort.Storable.Unchecked.Monadic as ArrayIO import qualified Data.Array.Comfort.Storable.Unchecked as Array import qualified Data.Array.Comfort.Shape as Shape import Data.Array.Comfort.Storable.Unchecked (Array(Array))@@ -56,7 +64,7 @@    Full vert horiz width nrhs a leastSquares    (Array shapeA@(MatrixShape.Full orderA extentA) a)-   (Array shapeB@(MatrixShape.Full orderB extentB) b) =+   (Array        (MatrixShape.Full orderB extentB) b) =   case Extent.fuse (Extent.generalizeWide $ Extent.transpose extentA) extentB of   Nothing -> error "leastSquares: height shapes mismatch"@@ -76,15 +84,12 @@       nrhsPtr <- Call.cint nrhs       (transPtr,aPtr) <- adjointA orderA (m*n) a       ldaPtr <- Call.leadingDim lda-      bPtr <- ContT $ withForeignPtr b       ldbPtr <- Call.leadingDim ldb-      let bSize = Shape.size shapeB-      btmpPtr <- Call.allocaArray bSize-      liftIO $ copyToColumnMajor orderB ldb nrhs bPtr btmpPtr+      bPtr <- copyToColumnMajorTemp orderB ldb nrhs b       liftIO $ withAutoWorkspaceInfo rankMsg "gels" $          LapackGen.gels transPtr-            mPtr nPtr nrhsPtr aPtr ldaPtr btmpPtr ldbPtr-      liftIO $ copySubMatrix ldx nrhs ldb btmpPtr ldx xPtr+            mPtr nPtr nrhsPtr aPtr ldaPtr bPtr ldbPtr+      liftIO $ copySubMatrix ldx nrhs ldb bPtr ldx xPtr  minimumNorm ::    (Extent.C vert, Extent.C horiz,@@ -184,12 +189,9 @@ leastSquaresMinimumNormIO rcond shapeX orderA a orderB b m n nrhs =    createHigherArray shapeX m n nrhs $ \(tmpPtr,ldtmp) -> do -   let aSize = m*n    let lda = m    evalContT $ do-      aPtr <- ContT $ withForeignPtr a-      atmpPtr <- Call.allocaArray aSize-      liftIO $ copyToColumnMajor orderA m n aPtr atmpPtr+      aPtr <- copyToColumnMajorTemp orderA m n a       ldaPtr <- Call.leadingDim lda       ldtmpPtr <- Call.leadingDim ldtmp       bPtr <- ContT $ withForeignPtr b@@ -197,7 +199,7 @@       jpvtPtr <- Call.allocaArray n       liftIO $ pokeArray jpvtPtr (replicate n 0)       rankPtr <- Call.alloca-      gelsy m n nrhs atmpPtr ldaPtr tmpPtr ldtmpPtr jpvtPtr rcond rankPtr+      gelsy m n nrhs aPtr ldaPtr tmpPtr ldtmpPtr jpvtPtr rcond rankPtr       liftIO $ peekCInt rankPtr  @@ -259,6 +261,82 @@          MatrixShape.fullHeight $ Array.shape a  +leastSquaresConstraint ::+   (Shape.C height, Eq height,+    Shape.C width, Eq width,+    Shape.C constraints, Eq constraints, Class.Floating a) =>+   General height width a -> Vector height a ->+   Wide constraints width a -> Vector constraints a ->+   Vector width a+leastSquaresConstraint+   (Array (MatrixShape.Full orderA extentA) a) c+   (Array (MatrixShape.Full orderB extentB) b) d =++ let sameShape name shape0 shape1 =+      if shape0 == shape1+         then shape0+         else error $ "leastSquaresConstraint: " ++ name ++ " shapes mismatch"+     width = sameShape "width" (Extent.width extentA) (Extent.width extentB)+ in+   Array.unsafeCreate width $ \xPtr -> do++   let height = sameShape "height" (Extent.height extentA) (Array.shape c)+   let constraints =+         sameShape "constraints" (Extent.height extentB) (Array.shape d)+   let m = Shape.size height+   let n = Shape.size width+   let p = Shape.size constraints+   evalContT $ do+      mPtr <- Call.cint m+      nPtr <- Call.cint n+      pPtr <- Call.cint p+      aPtr <- copyToColumnMajorTemp orderA m n a+      ldaPtr <- Call.leadingDim m+      bPtr <- copyToColumnMajorTemp orderB p n b+      ldbPtr <- Call.leadingDim p+      cPtr <- copyToTemp m (Array.buffer c)+      dPtr <- copyToTemp p (Array.buffer d)+      liftIO $ withAutoWorkspaceInfo rankMsg "gglse" $+         LapackGen.gglse+            mPtr nPtr pPtr aPtr ldaPtr bPtr ldbPtr cPtr dPtr xPtr++gaussMarkovLinearModel ::+   (Shape.C height, Eq height,+    Shape.C width, Eq width,+    Shape.C opt, Eq opt, Class.Floating a) =>+   Tall height width a -> General height opt a -> Vector height a ->+   (Vector width a, Vector opt a)+gaussMarkovLinearModel+   (Array (MatrixShape.Full orderA extentA) a)+   (Array (MatrixShape.Full orderB extentB) b) d =++   let width = Extent.width extentA in+   let opt = Extent.width extentB in+   Array.unsafeCreateWithSizeAndResult width $ \m xPtr -> do+   ArrayIO.unsafeCreateWithSize opt $ \p yPtr -> do++   let sameHeight shape0 shape1 =+         if shape0 == shape1+            then shape0+            else error $ "gaussMarkovLinearModel: height shapes mismatch"+       height =+         sameHeight (Array.shape d) $+         sameHeight (Extent.height extentA) (Extent.height extentB)+   let n = Shape.size height+   evalContT $ do+      mPtr <- Call.cint m+      nPtr <- Call.cint n+      pPtr <- Call.cint p+      aPtr <- copyToColumnMajorTemp orderA n m a+      ldaPtr <- Call.leadingDim n+      bPtr <- copyToColumnMajorTemp orderB n p b+      ldbPtr <- Call.leadingDim n+      dPtr <- copyToTemp n (Array.buffer d)+      liftIO $ withAutoWorkspaceInfo rankMsg "ggglm" $+         LapackGen.ggglm+            nPtr mPtr pPtr aPtr ldaPtr bPtr ldbPtr dPtr xPtr yPtr++ determinantAbsolute ::    (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width,     Class.Floating a) =>@@ -272,7 +350,13 @@ complement ::    (Shape.C height, Shape.C width, Class.Floating a) =>    Tall height width a -> Tall height ShapeInt a-complement a =-   Basic.dropColumns (Shape.size $ MatrixShape.fullWidth $ Array.shape a) $+complement = extractComplement . HH.fromMatrix++extractComplement ::+   (Shape.C height, Shape.C width, Class.Floating a) =>+   HH.Tall height width a -> Tall height ShapeInt a+extractComplement qr =+   Basic.dropColumns+      (Shape.size $ MatrixShape.splitWidth $ Array.shape $ HH.split_ qr) $    Basic.mapWidth (shapeInt . Shape.size) $ Square.toFull $-   HH.extractQ $ HH.fromMatrix a+   HH.extractQ qr

src/Numeric/LAPACK/Orthogonal/Private.hs view

@@ -474,9 +474,9 @@    (vert ~ Extent.Small, horiz ~ Extent.Small,     Shape.C height, height ~ width) =>       Divide.Solve (Hh vert horiz height width) where-   solveRight = ArrMatrix.lift1 . leastSquares . mapExtent Extent.generalizeWide+   solveRight = ArrMatrix.lift1 . leastSquares . mapExtent Extent.fromSquare    solveLeft =       flip $ \a -> ArrMatrix.lift1 $          Basic.adjoint .-         minimumNorm (mapExtent Extent.generalizeWide a) .+         minimumNorm (mapExtent Extent.fromSquare a) .          Basic.adjoint

src/Numeric/LAPACK/Permutation/Private.hs view

@@ -44,6 +44,27 @@ import Prelude hiding (odd)  +{- $setup+>>> import qualified Test.QuickCheck as QC+>>> import Test.Permutation (genPerm, genPivots)+>>>+>>> import qualified Numeric.LAPACK.Permutation as Perm+>>> import Numeric.LAPACK.Permutation (Permutation, Inversion(NonInverted), determinant, multiply, transpose)+>>> import Numeric.LAPACK.Matrix (ShapeInt)+>>>+>>> import qualified Data.Array.Comfort.Storable as Array+>>> import Data.Eq.HT (equating)+>>> import Data.Semigroup ((<>))+>>>+>>> import Control.Applicative (liftA2)+>>>+>>> genPerm2 :: QC.Gen (Permutation ShapeInt, Permutation ShapeInt)+>>> genPerm2 = do+>>>    nat <- QC.arbitrary+>>>    liftA2 (,) (genPerm nat) (genPerm nat)+-}++ newtype Permutation sh = Permutation (Vector (Shape sh) (Element sh))    deriving (Show) @@ -62,6 +83,9 @@ identity :: (Shape.C sh) => sh -> Permutation sh identity shape = Permutation $ CheckedArray.sample (Shape shape) id +{- |+prop> QC.forAll QC.arbitraryBoundedEnum $ \inv -> QC.forAll (QC.arbitrary >>= genPivots) $ \xs -> Array.toList xs == Array.toList (Perm.toPivots inv (Perm.fromPivots inv xs))+-} fromPivots ::    (Shape.C sh) =>    Inversion -> Vector (Shape sh) (Element sh) -> Permutation sh@@ -140,6 +164,9 @@ {- We could also count the cycles of even number. This might be a little faster. -}+{- |+prop> QC.forAll genPerm2 $ \(p0,p1) -> determinant (multiply p0 p1) == determinant p0 <> determinant p1+-} determinant :: (Shape.C sh) => Permutation sh -> Sign determinant =    (\oddp -> if oddp then Negative else Positive) .@@ -166,6 +193,9 @@ dropEven xs = xs  +{- |+prop> QC.forAll genPerm2 $ \(p0,p1) -> equating (Array.toList . Perm.toPivots NonInverted) (transpose $ multiply p0 p1) (multiply (transpose p1) (transpose p0))+-} transpose :: (Shape.C sh) => Permutation sh -> Permutation sh transpose (Permutation perm) =    Permutation $ runST (MutArray.unsafeFreeze =<< transposeToMutable perm)

src/Numeric/LAPACK/Private.hs view

@@ -195,6 +195,16 @@            else copySubMatrix m n m aPtr ldb bPtr  +copyToColumnMajorTemp ::+   (Class.Floating a) =>+   Order -> Int -> Int -> ForeignPtr a -> ContT r IO (Ptr a)+copyToColumnMajorTemp order m n fptr = do+   ptr <- ContT $ withForeignPtr fptr+   tmpPtr <- Call.allocaArray (m*n)+   liftIO $ copyToColumnMajor order m n ptr tmpPtr+   return tmpPtr++ pointerSeq :: (Storable a) => Int -> Ptr a -> [Ptr a] pointerSeq k ptr = iterate (flip advancePtr k) ptr

src/Numeric/LAPACK/Singular/Plain.hs view

@@ -30,7 +30,7 @@ import Numeric.LAPACK.Scalar (RealOf, zero) import Numeric.LAPACK.Private          (withAutoWorkspace, peekCInt, createHigherArray,-          copyToTemp, copyToColumnMajor, copyToSubColumnMajor)+          copyToTemp, copyToColumnMajorTemp, copyToSubColumnMajor)  import qualified Numeric.LAPACK.FFI.Complex as LapackComplex import qualified Numeric.LAPACK.FFI.Real as LapackReal@@ -325,15 +325,12 @@    createHigherArray shapeX m n nrhs $ \(tmpPtr,ldtmp) -> do     let mn = min m n-   let aSize = m*n    let lda = m    evalContT $ do       mPtr <- Call.cint m       nPtr <- Call.cint n       nrhsPtr <- Call.cint nrhs-      aPtr <- Call.allocaArray aSize-      liftIO $ withForeignPtr a $ \asrcPtr ->-         copyToColumnMajor orderA m n asrcPtr aPtr+      aPtr <- copyToColumnMajorTemp orderA m n a       ldaPtr <- Call.leadingDim lda       ldtmpPtr <- Call.leadingDim ldtmp       liftIO $ withForeignPtr b $ \bPtr ->

+ test-module.list view

@@ -0,0 +1,3 @@+Numeric.LAPACK.Permutation.Private+Numeric.LAPACK.Example.DividedDifference+Numeric.LAPACK.Example.EconomicAllocation

+ test/DocTest/Main.hs view

@@ -0,0 +1,14 @@+-- Do not edit! Automatically created with doctest-extract.+module DocTest.Main where++import qualified DocTest.Numeric.LAPACK.Permutation.Private+import qualified DocTest.Numeric.LAPACK.Example.DividedDifference+import qualified DocTest.Numeric.LAPACK.Example.EconomicAllocation++import qualified Test.DocTest.Driver as DocTest++main :: DocTest.T ()+main = do+    DocTest.Numeric.LAPACK.Permutation.Private.test+    DocTest.Numeric.LAPACK.Example.DividedDifference.test+    DocTest.Numeric.LAPACK.Example.EconomicAllocation.test

+ test/DocTest/Numeric/LAPACK/Example/DividedDifference.hs view

@@ -0,0 +1,47 @@+-- Do not edit! Automatically created with doctest-extract from src/Numeric/LAPACK/Example/DividedDifference.hs+{-# LINE 29 "src/Numeric/LAPACK/Example/DividedDifference.hs" #-}++module DocTest.Numeric.LAPACK.Example.DividedDifference where++import qualified Test.DocTest.Driver as DocTest++{-# LINE 30 "src/Numeric/LAPACK/Example/DividedDifference.hs" #-}+import     qualified Test.Utility as Util+import     Test.Utility (approxArray)++import     qualified Numeric.LAPACK.Vector as Vector+import     Numeric.LAPACK.Example.DividedDifference (dividedDifferencesMatrix)+import     Numeric.LAPACK.Matrix (ShapeInt, (#+#))+import     Numeric.LAPACK.Vector ((|+|))++import     qualified Data.Array.Comfort.Storable as Array++import     qualified Test.QuickCheck as QC++import     Control.Monad (liftM2)+import     Data.Tuple.HT (mapPair)+import     Data.Semigroup ((<>))++type     Vector = Vector.Vector ShapeInt Float++genDD     :: QC.Gen (Vector, (Vector, Vector))+genDD     = do+       (ys0,ys1) <-+          fmap (mapPair (Vector.autoFromList, Vector.autoFromList) .+                unzip . take 10) $+          QC.listOf $ liftM2 (,) (Util.genElement 10) (Util.genElement 10)+       xs <- Util.genDistinct 10 10 $ Array.shape ys0+       return (xs,(ys0,ys1))++test :: DocTest.T ()+test = do+ DocTest.printPrefix "Numeric.LAPACK.Example.DividedDifference:77: "+{-# LINE 77 "src/Numeric/LAPACK/Example/DividedDifference.hs" #-}+ DocTest.property+{-# LINE 77 "src/Numeric/LAPACK/Example/DividedDifference.hs" #-}+     (QC.forAll genDD $ \(xs, (ys0,ys1)) -> approxArray (dividedDifferencesMatrix xs (ys0|+|ys1)) (dividedDifferencesMatrix xs ys0 #+# dividedDifferencesMatrix xs ys1))+ DocTest.printPrefix "Numeric.LAPACK.Example.DividedDifference:78: "+{-# LINE 78 "src/Numeric/LAPACK/Example/DividedDifference.hs" #-}+ DocTest.property+{-# LINE 78 "src/Numeric/LAPACK/Example/DividedDifference.hs" #-}+     (QC.forAll genDD $ \(xs, (ys0,ys1)) -> approxArray (dividedDifferencesMatrix xs (Vector.mul ys0 ys1)) (dividedDifferencesMatrix xs ys0 <> dividedDifferencesMatrix xs ys1))

+ test/DocTest/Numeric/LAPACK/Example/EconomicAllocation.hs view

@@ -0,0 +1,23 @@+-- Do not edit! Automatically created with doctest-extract from src/Numeric/LAPACK/Example/EconomicAllocation.hs+{-# LINE 20 "src/Numeric/LAPACK/Example/EconomicAllocation.hs" #-}++module DocTest.Numeric.LAPACK.Example.EconomicAllocation where++import qualified Test.DocTest.Driver as DocTest++{-# LINE 21 "src/Numeric/LAPACK/Example/EconomicAllocation.hs" #-}+import     Numeric.LAPACK.Example.EconomicAllocation+import     Test.Utility (approxVector)++import     qualified Numeric.LAPACK.Vector as Vector+import     Numeric.LAPACK.Vector ((+++))++import     qualified Data.Array.Comfort.Storable as Array++test :: DocTest.T ()+test = do+ DocTest.printPrefix "Numeric.LAPACK.Example.EconomicAllocation:90: "+{-# LINE 90 "src/Numeric/LAPACK/Example/EconomicAllocation.hs" #-}+ DocTest.property+{-# LINE 90 "src/Numeric/LAPACK/Example/EconomicAllocation.hs" #-}+     (let result = iterated expenses0 balances0 in approxVector result $ compensated expenses0 balances0 +++ Vector.zero (Array.shape $ Vector.takeRight result))

+ test/DocTest/Numeric/LAPACK/Permutation/Private.hs view

@@ -0,0 +1,43 @@+-- Do not edit! Automatically created with doctest-extract from src/Numeric/LAPACK/Permutation/Private.hs+{-# LINE 47 "src/Numeric/LAPACK/Permutation/Private.hs" #-}++module DocTest.Numeric.LAPACK.Permutation.Private where++import qualified Test.DocTest.Driver as DocTest++{-# LINE 48 "src/Numeric/LAPACK/Permutation/Private.hs" #-}+import     qualified Test.QuickCheck as QC+import     Test.Permutation (genPerm, genPivots)++import     qualified Numeric.LAPACK.Permutation as Perm+import     Numeric.LAPACK.Permutation (Permutation, Inversion(NonInverted), determinant, multiply, transpose)+import     Numeric.LAPACK.Matrix (ShapeInt)++import     qualified Data.Array.Comfort.Storable as Array+import     Data.Eq.HT (equating)+import     Data.Semigroup ((<>))++import     Control.Applicative (liftA2)++genPerm2     :: QC.Gen (Permutation ShapeInt, Permutation ShapeInt)+genPerm2     = do+       nat <- QC.arbitrary+       liftA2 (,) (genPerm nat) (genPerm nat)++test :: DocTest.T ()+test = do+ DocTest.printPrefix "Numeric.LAPACK.Permutation.Private:87: "+{-# LINE 87 "src/Numeric/LAPACK/Permutation/Private.hs" #-}+ DocTest.property+{-# LINE 87 "src/Numeric/LAPACK/Permutation/Private.hs" #-}+     (QC.forAll QC.arbitraryBoundedEnum $ \inv -> QC.forAll (QC.arbitrary >>= genPivots) $ \xs -> Array.toList xs == Array.toList (Perm.toPivots inv (Perm.fromPivots inv xs)))+ DocTest.printPrefix "Numeric.LAPACK.Permutation.Private:168: "+{-# LINE 168 "src/Numeric/LAPACK/Permutation/Private.hs" #-}+ DocTest.property+{-# LINE 168 "src/Numeric/LAPACK/Permutation/Private.hs" #-}+     (QC.forAll genPerm2 $ \(p0,p1) -> determinant (multiply p0 p1) == determinant p0 <> determinant p1)+ DocTest.printPrefix "Numeric.LAPACK.Permutation.Private:197: "+{-# LINE 197 "src/Numeric/LAPACK/Permutation/Private.hs" #-}+ DocTest.property+{-# LINE 197 "src/Numeric/LAPACK/Permutation/Private.hs" #-}+     (QC.forAll genPerm2 $ \(p0,p1) -> equating (Array.toList . Perm.toPivots NonInverted) (transpose $ multiply p0 p1) (multiply (transpose p1) (transpose p0)))

test/Main.hs view

@@ -14,11 +14,12 @@ import qualified Test.Singular as Singular import qualified Test.Shape as Shape import qualified Test.Permutation as Permutation-import qualified Test.Example as Example+import qualified DocTest.Main as DocTestMain import Test.Format () import Test.Utility (Tagged(Tagged), prefix)  import qualified Test.QuickCheck as QC+import qualified Test.DocTest.Driver as DocTest  import Numeric.LAPACK.Scalar (RealOf) @@ -26,6 +27,15 @@  import Type.Base.Proxy (Proxy(Proxy)) +import System.Exit (exitFailure)++import Text.Printf (printf)++import qualified Control.Monad.Trans.Writer.Strict as MW+import qualified Control.Monad.Trans.Reader as MR+import Control.Monad.IO.Class (liftIO)+import Control.Monad (when, void)+ import qualified Data.List as List import Data.Complex (Complex) import Data.Tuple.HT (mapSnd)@@ -68,14 +78,28 @@ simpleTests :: [(String, QC.Property)] simpleTests =    prefix "Shape" Shape.tests ++-   prefix "Permutation" Permutation.tests ++-   prefix "Example" Example.tests ++    [] +++run :: MW.WriterT DocTest.Count IO () -> IO ()+run act = do+   count <- MW.execWriterT act+   putStrLn ""+   void $ printf "Total: %d\n" $ DocTest.numTotal count+   void $ printf "Failures: %d\n" $ DocTest.numFailures count+   when (DocTest.numFailures count > 0) exitFailure+ main :: IO () main =-   mapM_ (\(name,act) -> putStr (name ++ ": ") >> act) $+   run $+   (>> MR.runReaderT DocTestMain.main QC.stdArgs) $ -   map (mapSnd (QC.quickCheckWith (QC.stdArgs {QC.maxSuccess=200}))) tests+   mapM_+      (\(name,(args,act)) -> do+         liftIO (putStr (name ++ ": "))+         MR.runReaderT (DocTest.property act) args) $++   map (mapSnd ((,) (QC.stdArgs {QC.maxSuccess=200}))) tests    ++-   map (mapSnd QC.quickCheck) simpleTests+   map (mapSnd ((,) QC.stdArgs)) simpleTests

− test/Test/Example.hs

@@ -1,21 +0,0 @@-module Test.Example where--import qualified Test.Utility as Util--import qualified Numeric.LAPACK.Example.EconomicAllocation as Eco-import qualified Numeric.LAPACK.Vector as Vector--import qualified Data.Array.Comfort.Storable as Array--import qualified Test.QuickCheck as QC---tests :: [(String, QC.Property)]-tests =-   ("economicAllocation",-      QC.property $-         let iterated = Eco.iterated Eco.expenses0 Eco.balances0-         in Util.approxVector iterated-               (Vector.append (Eco.compensated Eco.expenses0 Eco.balances0) $-                Vector.zero $ Array.shape $ Vector.takeRight iterated)) :-   []

test/Test/Generator.hs view

@@ -429,6 +429,24 @@       Util.genArrayExtraDiag maxElem shape          (const $ fromReal <$> Util.genReal maxElem) +lscStack ::+   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>+   Matrix (ShapeInt:+:ShapeInt) ShapeInt a+      (Matrix.Tall (ShapeInt:+:ShapeInt) ShapeInt a)+lscStack =+   condition+      (Util.fullRankTall . Matrix.transpose .+       Matrix.wideFromGeneral . Matrix.takeTop . Matrix.fromFull) $+   condition Util.fullRankTall $+   mapGen Util.genArray $+   shapeFromDims (uncurry . MatrixShape.tall) $+   Base $ do+      height <- newVariable+      width <- newVariableWith $ Logic.between height+      return+         ((height,width),+          liftA2 (,) (Logic.query height) (Logic.query width))+  {- There cannot be a pure/point function.

test/Test/Logic.hs view

@@ -87,6 +87,31 @@    assignmentM (liftZeroBased $ \a -> choose (\ maxk -> (a,maxk))) va vb    assignmentM (liftZeroBased $ \b -> choose (\_maxk -> (0,b))) vb va +{-+We cannot split this into something like++a <!= c <> c <!= a+b++because the solver might choose in this order++a+b=6:+:4+c=5++and then it is left with a>c from which it cannot recover.+-}+between :: Variable s (ShapeInt:+:ShapeInt) -> Variable s ShapeInt -> System s+between vab vc  =  AppMn.Cons $ do+   assignmentM+      (\(Shape.ZeroBased a :+: Shape.ZeroBased b) ->+         Shape.ZeroBased <$> choose (\_maxk -> (a,a+b)))+      vab vc+   assignmentM+      (\(Shape.ZeroBased c) ->+         liftA2 (:+:)+            (Shape.ZeroBased <$> choose (\_maxk -> (0,c)))+            (Shape.ZeroBased <$> choose (\ maxk -> (c,maxk))))+      vc vab+  class (Shape.C dim) => Dim dim where chooseDim :: M s dim instance (i ~ Int) => Dim (Shape.ZeroBased i) where

test/Test/Orthogonal.hs view

@@ -1,13 +1,14 @@ {-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-} {-# LANGUAGE FlexibleContexts #-} module Test.Orthogonal (testsVar) where  import qualified Test.Divide as Divide import qualified Test.Generator as Gen import qualified Test.Utility as Util-import Test.Generator ((<#*#>), (<#\#>))+import Test.Generator ((<#*#>), (<#\#>), (<-*#>), (<#*|>), (<|=|>)) import Test.Utility-         (approx, approxReal, approxArrayTol, approxMatrix,+         (approx, approxReal, approxArrayTol, approxMatrix, approxVectorTol,           Tagged, isIdentity, isUnitary, maybeConjugate)  import qualified Numeric.LAPACK.Orthogonal.Householder as HH@@ -20,13 +21,16 @@ import qualified Numeric.LAPACK.Matrix as Matrix import qualified Numeric.LAPACK.Vector as Vector import Numeric.LAPACK.Matrix.Square (Square)-import Numeric.LAPACK.Matrix (General, ShapeInt, (#*#), (##*#), (#*##), (#\##))+import Numeric.LAPACK.Matrix+         (General, ShapeInt, (#*#), (##*#), (#*##), (#\##), (#*|))+import Numeric.LAPACK.Vector (Vector, (|+|), (|-|)) import Numeric.LAPACK.Scalar (RealOf, absolute, selectReal)  import qualified Numeric.Netlib.Class as Class  import qualified Data.Array.Comfort.Storable as Array import qualified Data.Array.Comfort.Shape as Shape+import Data.Array.Comfort.Shape ((:+:))  import Control.Applicative (liftA2, (<$>)) import Data.Semigroup ((<>))@@ -240,6 +244,166 @@ complementOrthogonal = isUnitary (selectReal 1e-3 1e-7) . Ortho.complement  +affineSpanFromKernel ::+   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>+   (Matrix.Wide ShapeInt ShapeInt a, Vector ShapeInt a) -> Bool+affineSpanFromKernel (a, by) =+   let b = Vector.take (Shape.size $ Matrix.height a) by+       y = Vector.drop (Shape.size $ Matrix.height a) by+       (c,d) = Ortho.affineSpanFromKernel a b+   in approxVectorTol+         (selectReal 1e-3 1e-7)+         b+         (a#*|(c#*|y|+|d))++affineKernelFromSpan ::+   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>+   (Matrix.Tall ShapeInt ShapeInt a, Vector ShapeInt a, Vector ShapeInt a) ->+   Bool+affineKernelFromSpan (c,y,d) =+   let (a,b) = Ortho.affineKernelFromSpan c d+   in approxVectorTol+         (selectReal 1e-3 1e-7)+         b+         (a#*|(c#*|y|+|d))+++projectHit ::+   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>+   (Matrix.Wide ShapeInt ShapeInt a,+    Vector ShapeInt a,+    Vector ShapeInt a) ->+   Bool+projectHit (b,x,d) =+   approxVectorTol+      (selectReal 1e-3 1e-9)+      d+      (b #*| Ortho.project b d x)+++leastSquaresNoConstraint ::+   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>+   (Matrix.Tall ShapeInt ShapeInt a, Vector ShapeInt a) -> Bool+leastSquaresNoConstraint (a, b) =+   approxVectorTol+      (selectReal 0.1 1e-7)+      (ArrMatrix.unliftColumn MatrixShape.ColumnMajor (Ortho.leastSquares a) b)+      (Ortho.leastSquaresConstraint+         (Matrix.fromFull a) b+         (Matrix.zero $+          MatrixShape.wide MatrixShape.ColumnMajor Shape.Zero (Matrix.width a))+         (Vector.zero Shape.Zero))++leastSquaresConstraintUnique ::+   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>+   (Vector ShapeInt a,+    Matrix.General ShapeInt ShapeInt a, Matrix.Square ShapeInt a,+    Vector ShapeInt a) ->+   Bool+leastSquaresConstraintUnique (c,a,b,d) =+   approxVectorTol+      (selectReal 0.1 1e-7)+      d+      (b #*| Ortho.leastSquaresConstraint a c (Matrix.generalizeTall b) d)++splitLSCStack ::+   (Shape.C height, Shape.C constraints, Shape.C width, Class.Floating a) =>+   Matrix.Tall (constraints:+:height) width a ->+   Vector (constraints:+:height) a ->+   ((Matrix.General height width a, Vector height a),+    (Matrix.Wide constraints width a, Vector constraints a))+splitLSCStack baTall dc =+   let ba = Matrix.fromFull baTall+       b = Matrix.wideFromGeneral $ Matrix.takeTop ba+       a = Matrix.takeBottom ba+       (d,c) = (Vector.takeLeft dc, Vector.takeRight dc)+   in ((a,c),(b,d))++leastSquaresConstraintAdmissible ::+   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>+   (Matrix.Tall (ShapeInt:+:ShapeInt) ShapeInt a,+    Vector (ShapeInt:+:ShapeInt) a) ->+   Bool+leastSquaresConstraintAdmissible (ba,dc) =+   let ((a,c),(b,d)) = splitLSCStack ba dc+   in approxVectorTol+         (selectReal 0.1 1e-7)+         d+         (b #*| Ortho.leastSquaresConstraint a c b d)++leastSquaresConstraintMinimal ::+   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>+   (Matrix.Tall (ShapeInt:+:ShapeInt) ShapeInt a,+    Vector ShapeInt a,+    Vector (ShapeInt:+:ShapeInt) a) ->+   Bool+leastSquaresConstraintMinimal (ba,x,dc) =+   let ((a,c),(b,d)) = splitLSCStack ba dc+   in Vector.norm2 (c |-| a #*| Ortho.leastSquaresConstraint a c b d)+      <=+      Vector.norm2 (c |-| a #*| Ortho.project b d x) + selectReal 1e-1 1e-10+++gaussMarkovLinearModelMinimumNorm ::+   (Class.Floating a, RealOf a ~ ar, Class.Real ar, Eq a) =>+   (Matrix.Wide ShapeInt ShapeInt a, Vector ShapeInt a) -> Bool+gaussMarkovLinearModelMinimumNorm (b, d) =+   let (x,y) =+         Ortho.gaussMarkovLinearModel+            (Matrix.zero $+             MatrixShape.tall MatrixShape.ColumnMajor+               (Matrix.height b) Shape.Zero)+            (Matrix.fromFull b) d+   in x == Vector.zero Shape.Zero+      &&+      approxVectorTol+         (selectReal 0.1 1e-7)+         y+         (ArrMatrix.unliftColumn MatrixShape.ColumnMajor+            (Ortho.minimumNorm b) d)++gaussMarkovLinearModelUnique ::+   (Class.Floating a, RealOf a ~ ar, Class.Real ar, Eq a) =>+   (Matrix.Square ShapeInt a, Vector ShapeInt a) -> Bool+gaussMarkovLinearModelUnique (a,d) =+   let (x,y) =+         Ortho.gaussMarkovLinearModel+            (Matrix.generalizeWide a)+            (Matrix.zero $+             MatrixShape.general MatrixShape.ColumnMajor+               (Matrix.height a) Shape.Zero)+            d+   in y == Vector.zero Shape.Zero+      &&+      approxVectorTol (selectReal 0.1 1e-7) d (a #*| x)++gaussMarkovLinearModelAdmissible ::+   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>+   (Matrix.Wide ShapeInt (ShapeInt:+:ShapeInt) a, Vector ShapeInt a) -> Bool+gaussMarkovLinearModelAdmissible (abWide,d) =+   let ab = Matrix.fromFull abWide+       a = Matrix.tallFromGeneral $ Matrix.takeLeft ab+       b = Matrix.takeRight ab+   in approxVectorTol+         (selectReal 0.1 1e-7)+         d+         (ab #*| uncurry Vector.append (Ortho.gaussMarkovLinearModel a b d))++gaussMarkovLinearModelMinimal ::+   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>+   (Matrix.Wide ShapeInt (ShapeInt:+:ShapeInt) a,+    Vector ShapeInt a,+    Vector (ShapeInt:+:ShapeInt) a) -> Bool+gaussMarkovLinearModelMinimal (abWide,d,xy) =+   let ab = Matrix.fromFull abWide+       a = Matrix.tallFromGeneral $ Matrix.takeLeft ab+       b = Matrix.takeRight ab+   in Vector.norm2 (snd $ Ortho.gaussMarkovLinearModel a b d)+      <=+      Vector.norm2 (Vector.takeRight $ Ortho.project abWide d xy)+         + selectReal 1e-3 1e-10++ householderReconstruction ::    (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>    Matrix.General ShapeInt ShapeInt a -> Bool@@ -392,6 +556,53 @@       checkForAll Gen.tall complementBiorthogonal) :    ("complementOrthogonal",       checkForAll Gen.tall complementOrthogonal) :++   ("affineSpanFromKernel",+      checkForAll+         ((,) <$> Gen.fullRankWide <#*|> Gen.vector) affineSpanFromKernel) :+   ("affineKernelFromSpan",+      checkForAll+         ((,,) <$> Gen.tall <#*|> Gen.vector <|=|> Gen.vector)+         affineKernelFromSpan) :++   ("projectHit",+      checkForAll+         ((,,) <$> Gen.fullRankWide <#*|> Gen.vector <|=|> Gen.vector)+         projectHit) :+   ("leastSquaresNoConstraint",+      checkForAll+         ((,) <$> Gen.transpose Gen.fullRankTall <#*|> Gen.vector)+         leastSquaresNoConstraint) :+   ("leastSquaresConstraintUnique",+      checkForAll+         ((,,,) <$>+            Gen.vector <-*#> Gen.matrix <-*#> Gen.invertible <|=|> Gen.vector)+         leastSquaresConstraintUnique) :+   ("leastSquaresConstraintAdmissible",+      checkForAll+         ((,) <$> Gen.transpose Gen.lscStack <#*|> Gen.vector)+         leastSquaresConstraintAdmissible) :+   ("leastSquaresConstraintMinimal",+      checkForAll+         ((,,) <$> Gen.lscStack <#*|> Gen.vector <|=|> Gen.vector)+         leastSquaresConstraintMinimal) :+   ("gaussMarkovLinearModelMinimumNorm",+      checkForAll+         ((,) <$> Gen.transpose Gen.fullRankWide <#*|> Gen.vector)+         gaussMarkovLinearModelMinimumNorm) :+   ("gaussMarkovLinearModelUnique",+      checkForAll+         ((,) <$> Gen.invertible <#*|> Gen.vector)+         gaussMarkovLinearModelUnique) :+   ("gaussMarkovLinearModelAdmissible",+      checkForAll+         ((,) <$> fmap Matrix.transpose Gen.lscStack <#*|> Gen.vector)+         gaussMarkovLinearModelAdmissible) :+   ("gaussMarkovLinearModelMinimal",+      checkForAll+         ((,,) <$> fmap Matrix.transpose Gen.lscStack+               <#*|> Gen.vector <|=|> Gen.vector)+         gaussMarkovLinearModelMinimal) :     ("triangularLeastSquares",       checkForAll genFullRankTallRHS triangularLeastSquares) :

test/Test/Permutation.hs view

@@ -20,13 +20,10 @@  import qualified Numeric.Netlib.Class as Class -import qualified Data.Array.Comfort.Storable as Array import qualified Data.Array.Comfort.Shape as Shape  import Control.Monad (forM)-import Control.Applicative (liftA2, (<$>))--import Data.Semigroup ((<>))+import Control.Applicative ((<$>))  import qualified Test.QuickCheck as QC @@ -44,27 +41,6 @@ genPerm = fmap (Perm.fromPivots NonInverted) . genPivots  -permutationPivots :: Inversion -> Pivots -> Bool-permutationPivots inv xs =-   Array.toList (Perm.toPivots inv (Perm.fromPivots inv xs))-   ==-   Array.toList xs--determinantMultiply :: (Permutation ShapeInt, Permutation ShapeInt) -> Bool-determinantMultiply (p0,p1) =-   Perm.determinant (Perm.multiply p0 p1)-   ==-   Perm.determinant p0 <> Perm.determinant p1--transposeMultiply :: (Permutation ShapeInt, Permutation ShapeInt) -> Bool-transposeMultiply (p0,p1) =-   (Array.toList $ Perm.toPivots NonInverted $-    Perm.transpose (Perm.multiply p0 p1))-   ==-   (Array.toList $ Perm.toPivots NonInverted $-    Perm.multiply (Perm.transpose p1) (Perm.transpose p0))-- genPermutation :: (Dim sh) => Gen.Matrix sh sh a (Permutation sh) genPermutation =    flip Gen.mapGen Gen.squareDim $ \_maxElem sh ->@@ -117,28 +93,6 @@    (Class.Floating a, Eq a) => Matrix (Permutation ShapeInt) a -> Bool determinantNumber p =    PermMatrix.determinant p == Square.determinant (PermMatrix.toMatrix p)---tests :: [(String, QC.Property)]-tests =-   ("permutationPivots",-      QC.property $-         QC.forAll QC.arbitraryBoundedEnum $ \inv ->-         QC.forAll (QC.arbitrary >>= genPivots) $ \pivot ->-            permutationPivots inv pivot) :-   ("determinantMultiply",-      QC.property $-         QC.forAll-            (do nat <- QC.arbitrary-                liftA2 (,) (genPerm nat) (genPerm nat))-            determinantMultiply) :-   ("transposeMultiply",-      QC.property $-         QC.forAll-            (do nat <- QC.arbitrary-                liftA2 (,) (genPerm nat) (genPerm nat))-            transposeMultiply) :-   []

test/Test/Shape.hs view

@@ -99,19 +99,19 @@    Banded (MatrixShape.Banded sub super vert horiz height width)  instance-   (Extent.C horiz, Extent.C vert,+   (Extent.C vert, Extent.C horiz,     Show height, Show width, Shape.C height, Shape.C width) =>       Show (Banded vert horiz height width) where    showsPrec p (Banded sh) = showsPrec p sh  instance-   (Extent.C horiz, Extent.C vert, Shape.C height, Shape.C width) =>+   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width) =>       Shape.C (Banded vert horiz height width) where    size (Banded sh) = Shape.size sh    uncheckedSize (Banded sh) = Shape.uncheckedSize sh  instance-   (Extent.C horiz, Extent.C vert,+   (Extent.C vert, Extent.C horiz,     Shape.Indexed height, Shape.Indexed width) =>       Shape.Indexed (Banded vert horiz height width) where    type Index (Banded vert horiz height width) =@@ -125,7 +125,7 @@    uncheckedSizeOffset (Banded sh) = Shape.uncheckedSizeOffset sh  instance-   (Extent.C horiz, Extent.C vert,+   (Extent.C vert, Extent.C horiz,     Shape.InvIndexed height, Shape.InvIndexed width) =>       Shape.InvIndexed (Banded vert horiz height width) where