lapack 0.2.1 → 0.2.2

raw patch · 29 files changed

+1086/−216 lines, 29 filesdep +monoid-transformerPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies added: monoid-transformer

API changes (from Hackage documentation)

Files

• lapack.cabal +6/−2
• src/Numeric/LAPACK/Matrix.hs +149/−10
• src/Numeric/LAPACK/Matrix/Banded/Basic.hs +2/−1
• src/Numeric/LAPACK/Matrix/BandedHermitian/Basic.hs +2/−1
• src/Numeric/LAPACK/Matrix/Basic.hs +46/−0
• src/Numeric/LAPACK/Matrix/Divide.hs +3/−3
• src/Numeric/LAPACK/Matrix/Extent/Private.hs +35/−3
• src/Numeric/LAPACK/Matrix/Hermitian/Basic.hs +57/−2
• src/Numeric/LAPACK/Matrix/Indexed.hs +86/−0
• src/Numeric/LAPACK/Matrix/Multiply.hs +8/−3
• src/Numeric/LAPACK/Matrix/Private.hs +5/−8
• src/Numeric/LAPACK/Matrix/Shape.hs +6/−2
• src/Numeric/LAPACK/Matrix/Shape/Box.hs +18/−0
• src/Numeric/LAPACK/Matrix/Shape/Private.hs +44/−18
• src/Numeric/LAPACK/Matrix/Square/Basic.hs +2/−1
• src/Numeric/LAPACK/Matrix/Triangular/Basic.hs +85/−7
• src/Numeric/LAPACK/Matrix/Triangular/Private.hs +75/−0
• src/Numeric/LAPACK/Orthogonal.hs +5/−12
• src/Numeric/LAPACK/Vector.hs +28/−2
• test/Test/Banded.hs +10/−0
• test/Test/BandedHermitian.hs +10/−3
• test/Test/Generator.hs +91/−68
• test/Test/Hermitian.hs +19/−2
• test/Test/Indexed.hs +53/−0
• test/Test/Logic.hs +14/−11
• test/Test/Matrix.hs +62/−44
• test/Test/Triangular.hs +88/−8
• test/Test/Utility.hs +43/−2
• test/Test/Vector.hs +34/−3

@@ -1,5 +1,5 @@ Name:             lapack-Version:          0.2.1+Version:          0.2.2 License:          BSD3 License-File:     LICENSE Author:           Henning Thielemann <haskell@henning-thielemann.de>@@ -45,7 +45,7 @@   default:     False  Source-Repository this-  Tag:         0.2.1+  Tag:         0.2.2   Type:        darcs   Location:    http://hub.darcs.net/thielema/lapack/ @@ -75,6 +75,7 @@   Exposed-Modules:     Numeric.LAPACK.Matrix     Numeric.LAPACK.Matrix.Extent+    Numeric.LAPACK.Matrix.Shape.Box     Numeric.LAPACK.Matrix.Shape     Numeric.LAPACK.Matrix.Square     Numeric.LAPACK.Matrix.Hermitian@@ -120,6 +121,7 @@     Numeric.LAPACK.Matrix.Extent.Kind     Numeric.LAPACK.Matrix.Multiply     Numeric.LAPACK.Matrix.Divide+    Numeric.LAPACK.Matrix.Indexed     Numeric.LAPACK.Matrix.Basic     Numeric.LAPACK.Matrix.Private     Numeric.LAPACK.Private@@ -139,6 +141,7 @@     QuickCheck >=2.5 && <3,     ChasingBottoms >=1.2.2 && <1.4,     transformers,+    monoid-transformer >=0.0.3 && <0.1,     semigroups >=0.18.3 && <0.19,     non-empty >=0.3.1,     utility-ht,@@ -149,6 +152,7 @@   Main-Is:          Main.hs   Other-Modules:     Test.Shape+    Test.Indexed     Test.Permutation     Test.Vector     Test.Matrix

@@ -11,12 +11,15 @@    fromList,    mapExtent, fromFull,    generalizeTall, generalizeWide,+   mapHeight, mapWidth,    identity,    diagonal,    fromRowsNonEmpty,    fromRowArray,    fromRows,    fromColumnsNonEmpty, fromColumnArray, fromColumns,    Basic.singleRow,   Basic.singleColumn,    Basic.flattenRow,  Basic.flattenColumn,+   Basic.liftRow,     Basic.liftColumn,+   Basic.unliftRow,   Basic.unliftColumn,    toRows, toColumns,    toRowArray, toColumnArray,    takeRow, takeColumn,@@ -24,6 +27,8 @@    dropRows, dropColumns, dropEqually,    takeTopRows, takeBottomRows,    takeLeftColumns, takeRightColumns,+   takeRowArray, takeColumnArray,+   swapRows, swapColumns,    reverseRows, reverseColumns,    fromRowMajor, toRowMajor, flatten,    forceOrder, adaptOrder,@@ -46,6 +51,7 @@    Multiply, (<#>),    MultiplyLeft, (<#),    MultiplyRight, (#>),+   Indexed, (#!),     Solve, solve, solveVector,    Inverse, inverse,@@ -62,6 +68,7 @@ import Numeric.LAPACK.Matrix.Multiply          (Multiply((<#>)), MultiplyLeft((<#)), MultiplyRight((#>)),           multiplyVector, multiply, multiplyVectorUnchecked)+import Numeric.LAPACK.Matrix.Indexed (Indexed((#!))) import Numeric.LAPACK.Matrix.Divide          (Solve(solve), solveVector, Inverse(inverse)) import Numeric.LAPACK.Matrix.Basic@@ -81,6 +88,7 @@  import qualified Data.Array.Comfort.Boxed as BoxedArray import qualified Data.Array.Comfort.Storable.Unchecked as Array+import qualified Data.Array.Comfort.Storable as CheckedArray import qualified Data.Array.Comfort.Shape as Shape import Data.Array.Comfort.Storable.Unchecked (Array(Array)) import Data.Array.Comfort.Shape ((:+:)((:+:)))@@ -94,10 +102,13 @@  import Control.Monad.Trans.Cont (ContT(ContT), evalContT) import Control.Monad.IO.Class (liftIO)+import Control.Monad (when, mfilter) +import qualified Data.NonEmpty.Mixed as NonEmptyM import qualified Data.NonEmpty as NonEmpty import Data.Complex (Complex) import Data.Foldable (forM_)+import Data.Maybe (listToMaybe) import Data.Bool.HT (if')  @@ -138,7 +149,7 @@    (Shape.C height, Shape.C width, Storable a) =>    height -> width -> [a] -> General height width a fromList height width =-   Array.fromList (MatrixShape.general RowMajor height width)+   CheckedArray.fromList (MatrixShape.general RowMajor height width)   identity ::@@ -152,6 +163,40 @@ diagonal = Square.toGeneral . Square.diagonal  +{- |+The number of rows must be maintained by the height mapping function.+-}+mapHeight ::+   (Shape.C heightA, Shape.C heightB,+    Extent.GeneralTallWide vert horiz,+    Extent.GeneralTallWide horiz vert) =>+   (heightA -> heightB) ->+   Full vert horiz heightA width a ->+   Full vert horiz heightB width a+mapHeight f = Basic.mapHeight $ withSizeCheck "mapHeight" f++{- |+The number of columns must be maintained by the width mapping function.+-}+mapWidth ::+   (Shape.C widthA, Shape.C widthB,+    Extent.GeneralTallWide vert horiz,+    Extent.GeneralTallWide horiz vert) =>+   (widthA -> widthB) ->+   Full vert horiz height widthA a ->+   Full vert horiz height widthB a+mapWidth f = Basic.mapWidth $ withSizeCheck "mapWidth" f++withSizeCheck ::+   (Shape.C sha, Shape.C shb) =>+   String -> (sha -> shb) -> sha -> shb+withSizeCheck name f sizeA =+   let sizeB = f sizeA+   in if Shape.size sizeA == Shape.size sizeB+         then sizeB+         else error $ name ++ ": sizes mismatch"++ fromRowsNonEmpty ::    (Shape.C width, Eq width, Storable a) =>    NonEmpty.T [] (Vector width a) -> General ZeroInt width a@@ -162,7 +207,7 @@    (Shape.C height, Shape.C width, Eq width, Storable a) =>    width -> BoxedArray.Array height (Vector width a) -> General height width a fromRowArray width rows =-   Array.reshape (MatrixShape.general RowMajor (BoxedArray.shape rows) width) $+   Basic.mapHeight (const $ BoxedArray.shape rows) $    fromRows width $ BoxedArray.toList rows  fromRows ::@@ -183,8 +228,7 @@    (Shape.C height, Eq height, Shape.C width, Storable a) =>    height -> BoxedArray.Array width (Vector height a) -> General height width a fromColumnArray height columns =-   Array.reshape-      (MatrixShape.general ColumnMajor height (BoxedArray.shape columns)) $+   Basic.mapWidth (const $ BoxedArray.shape columns) $    fromColumns height $ BoxedArray.toList columns  fromColumns ::@@ -223,19 +267,21 @@  toRowArray ::    (Extent.C vert, Extent.C horiz,-    Shape.Indexed height, Shape.C width, Class.Floating a) =>+    Shape.C height, Shape.C width, Class.Floating a) =>    Full vert horiz height width a -> BoxedArray.Array height (Vector width a) toRowArray a =-   let height = Matrix.height a-   in BoxedArray.fromList height $ map (takeRow a) $ Shape.indices height+   let ad = Basic.mapHeight Shape.Deferred $ fromFull a+   in BoxedArray.fromList (Matrix.height a) $+      map (takeRow ad) $ Shape.indices (Matrix.height ad)  toColumnArray ::    (Extent.C vert, Extent.C horiz,-    Shape.C height, Shape.Indexed width, Class.Floating a) =>+    Shape.C height, Shape.C width, Class.Floating a) =>    Full vert horiz height width a -> BoxedArray.Array width (Vector height a) toColumnArray a =-   let width = Matrix.width a-   in BoxedArray.fromList width $ map (takeColumn a) $ Shape.indices width+   let ad = Basic.mapWidth Shape.Deferred $ fromFull a+   in BoxedArray.fromList (Matrix.width a) $+      map (takeColumn ad) $ Shape.indices (Matrix.width ad)   takeRow ::@@ -429,6 +475,40 @@                heightA (advancePtr aPtr (left*heightA+top)) heightB bPtr  +swapRows ::+   (Extent.C vert, Extent.C horiz,+    Shape.Indexed height, Shape.C width, Class.Floating a) =>+   Shape.Index height -> Shape.Index height ->+   Full vert horiz height width a -> Full vert horiz height width a+swapRows i j (Array shape@(MatrixShape.Full order extent) a) =+   Array.unsafeCreateWithSize shape $ \blockSize bPtr -> evalContT $ do+      let (height,width) = Extent.dimensions extent+      let m = Shape.size height+      let n = Shape.size width+      nPtr <- Call.cint n+      aPtr <- ContT $ withForeignPtr a+      let offsetI = Shape.offset height i+      let offsetJ = Shape.offset height j+      let (incVert,incHoriz) =+            case order of+               RowMajor -> (n,1)+               ColumnMajor -> (1,m)+      incPtr <- Call.cint incHoriz+      liftIO $ do+         copyBlock blockSize aPtr bPtr+         when (offsetI/=offsetJ) $+            BlasGen.swap nPtr+               (advancePtr bPtr (incVert*offsetI)) incPtr+               (advancePtr bPtr (incVert*offsetJ)) incPtr++swapColumns ::+   (Extent.C vert, Extent.C horiz,+    Shape.C height, Shape.Indexed width, Class.Floating a) =>+   Shape.Index width -> Shape.Index width ->+   Full vert horiz height width a -> Full vert horiz height width a+swapColumns i j = transpose . swapRows i j . transpose++ -- alternative: laswp reverseRows ::    (Extent.C vert, Extent.C horiz, Shape.C width, Class.Floating a) =>@@ -454,6 +534,65 @@    (Extent.C vert, Extent.C horiz, Shape.C height, Class.Floating a) =>    Full vert horiz height ZeroInt a -> Full vert horiz height ZeroInt a reverseColumns = transpose . reverseRows . transpose+++{-+The function is optimized for blocks of consecutive rows.+For scattered rows in column major order+the function has quite ugly memory access patterns.+-}+takeRowArray ::+   (Shape.Indexed height, Shape.C width, Shape.C sh, Class.Floating a) =>+   BoxedArray.Array sh (Shape.Index height) ->+   General height width a -> General sh width a+takeRowArray ixs (Array (MatrixShape.Full order extent) a) =+   let (heightA, width) = Extent.dimensions extent+       heightB = BoxedArray.shape ixs+       offsets = map (Shape.offset heightA) $ BoxedArray.toList ixs+       startBlocks blocks = zip (scanl (+) 0 $ map fst blocks) blocks+       ma = Shape.size heightA+       mb = Shape.size heightB+       n = Shape.size width+   in Array.unsafeCreate (MatrixShape.general order heightB width) $ \bPtr ->+      withForeignPtr a $ \aPtr ->+      case order of+         RowMajor -> do+            forM_ (startBlocks $ chopRowBlocks offsets) $+               \(dest, (numRows, (start,step))) ->+                  copySubMatrix n numRows+                     (step*n) (advancePtr aPtr (start*n))+                     n (advancePtr bPtr (dest*n))+         ColumnMajor -> do+            forM_ (startBlocks $ chopColumnBlocks offsets) $+               \(dest, (numRows, start)) ->+                  copySubMatrix numRows n+                     ma (advancePtr aPtr start)+                     mb (advancePtr bPtr dest)++chopRowBlocks :: (Integral i) => [i] -> [(Int,(i,i))]+chopRowBlocks =+   let go [] = []+       go is@(i0:is0) =+         case mfilter (i0<) $ listToMaybe is0 of+            Nothing -> (1,(i0,0)) : go is0+            Just i1 ->+               let (consecutive,remainder) =+                     span (uncurry (==)) $ zip [i0,i1..] is+               in (length consecutive, (i0,i1-i0)) : go (map snd remainder)+   in go++chopColumnBlocks :: (Integral i) => [i] -> [(Int,i)]+chopColumnBlocks =+   map (\is -> (length $ NonEmpty.flatten is, NonEmpty.head is)) .+   NonEmptyM.groupBy (\i j -> i+1 == j)+++takeColumnArray ::+   (Shape.C height, Shape.Indexed width, Shape.C sh, Class.Floating a) =>+   BoxedArray.Array sh (Shape.Index width) ->+   General height width a -> General height sh a+takeColumnArray ixs = transpose . takeRowArray ixs . transpose+   fromRowMajor ::

@@ -51,6 +51,7 @@ import Type.Base.Proxy (Proxy(Proxy))  import qualified Data.Array.Comfort.Storable.Unchecked as Array+import qualified Data.Array.Comfort.Storable as CheckedArray import qualified Data.Array.Comfort.Shape as Shape import Data.Array.Comfort.Storable.Unchecked (Array(Array)) @@ -120,7 +121,7 @@    Extent.Extent vert horiz height width -> [a] ->    Banded sub super vert horiz height width a fromListGen offDiag order extent =-   Array.fromList (MatrixShape.Banded offDiag order extent)+   CheckedArray.fromList (MatrixShape.Banded offDiag order extent)   mapExtent ::

@@ -50,6 +50,7 @@ import Type.Base.Proxy (Proxy(Proxy))  import qualified Data.Array.Comfort.Storable.Unchecked as Array+import qualified Data.Array.Comfort.Storable as CheckedArray import qualified Data.Array.Comfort.Shape as Shape import Data.Array.Comfort.Storable.Unchecked (Array(Array)) @@ -80,7 +81,7 @@    UnaryProxy offDiag -> Order -> size -> [a] ->    BandedHermitian offDiag size a fromList numOff order size =-   Array.fromList (MatrixShape.BandedHermitian numOff order size)+   CheckedArray.fromList (MatrixShape.BandedHermitian numOff order size)  identity ::    (Shape.C sh, Class.Floating a) => sh -> Diagonal sh a

@@ -31,7 +31,29 @@    Full vert horiz height width a -> Full horiz vert width height a transpose = Array.mapShape MatrixShape.transpose +mapHeight ::+   (Extent.GeneralTallWide vert horiz,+    Extent.GeneralTallWide horiz vert) =>+   (heightA -> heightB) ->+   Full vert horiz heightA width a ->+   Full vert horiz heightB width a+mapHeight f =+   Array.mapShape+      (\(MatrixShape.Full order extent) ->+         MatrixShape.Full order $ Extent.mapHeight f extent) +mapWidth ::+   (Extent.GeneralTallWide vert horiz,+    Extent.GeneralTallWide horiz vert) =>+   (widthA -> widthB) ->+   Full vert horiz height widthA a ->+   Full vert horiz height widthB a+mapWidth f =+   Array.mapShape+      (\(MatrixShape.Full order extent) ->+         MatrixShape.Full order $ Extent.mapWidth f extent)++ singleRow :: Order -> Vector width a -> General () width a singleRow order = Array.mapShape (MatrixShape.general order ()) @@ -43,6 +65,30 @@  flattenColumn :: General height () a -> Vector height a flattenColumn = Array.mapShape MatrixShape.fullHeight++liftRow ::+   Order ->+   (Vector height0 a -> Vector height1 b) ->+   General () height0 a -> General () height1 b+liftRow order f = singleRow order . f . flattenRow++liftColumn ::+   Order ->+   (Vector height0 a -> Vector height1 b) ->+   General height0 () a -> General height1 () b+liftColumn order f = singleColumn order . f . flattenColumn++unliftRow ::+   Order ->+   (General () height0 a -> General () height1 b) ->+   Vector height0 a -> Vector height1 b+unliftRow order f = flattenRow . f . singleRow order++unliftColumn ::+   Order ->+   (General height0 () a -> General height1 () b) ->+   Vector height0 a -> Vector height1 b+unliftColumn order f = flattenColumn . f . singleColumn order   forceRowMajor ::

@@ -16,8 +16,8 @@  import qualified Numeric.LAPACK.Matrix.Basic as Basic import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape+import qualified Numeric.LAPACK.Matrix.Shape.Box as Box import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent-import Numeric.LAPACK.Matrix.Shape.Private (HeightOf) import Numeric.LAPACK.Matrix.Extent.Private (Small) import Numeric.LAPACK.Matrix.Private (Full) import Numeric.LAPACK.Vector (Vector)@@ -32,7 +32,7 @@  class (Shape.C shape) => Solve shape where    solve ::-      (Class.Floating a, HeightOf shape ~ height, Eq height,+      (Class.Floating a, Box.HeightOf shape ~ height, Eq height,        Extent.C horiz, Extent.C vert, Shape.C nrhs) =>       Array shape a ->       Full vert horiz height nrhs a -> Full vert horiz height nrhs a@@ -41,7 +41,7 @@    inverse :: (Class.Floating a) => Array shape a -> Array shape a  solveVector ::-   (Solve shape, HeightOf shape ~ height, Eq height, Class.Floating a) =>+   (Solve shape, Box.HeightOf shape ~ height, Eq height, Class.Floating a) =>    Array shape a -> Vector height a -> Vector height a solveVector m =    Basic.flattenColumn . solve m . Basic.singleColumn MatrixShape.ColumnMajor

@@ -1,4 +1,5 @@ {-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE MultiParamTypeClasses #-} {-# LANGUAGE GADTs #-} module Numeric.LAPACK.Matrix.Extent.Private where @@ -332,11 +333,11 @@       (Widen (\(Extent o x) -> Extent o (x{EK.generalHeight = h})))  -newtype Adapt height width vert horiz =+newtype Adapt heightA widthA heightB widthB vert horiz =    Adapt {       getAdapt ::-         Extent vert horiz height width ->-         Extent vert horiz height width+         Extent vert horiz heightA widthA ->+         Extent vert horiz heightB widthB    }  reduceConsistent ::@@ -350,6 +351,37 @@       (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))+++class (C vert, C horiz) => GeneralTallWide vert horiz where+   switchTagGTW :: f Small Big -> f Big Small -> f Big Big -> f vert horiz++instance GeneralTallWide Small Big where switchTagGTW f _ _ = f+instance GeneralTallWide Big Small where switchTagGTW _ f _ = f+instance GeneralTallWide Big Big where switchTagGTW _ _ f = f++mapHeight ::+   (GeneralTallWide vert horiz) =>+   (heightA -> heightB) ->+   Extent vert horiz heightA width -> Extent vert horiz heightB width+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))++mapWidth ::+   (GeneralTallWide vert horiz) =>+   (widthA -> widthB) ->+   Extent vert horiz height widthA -> Extent vert horiz height widthB+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)))+   newtype Fuse height fuse width vert horiz =

@@ -1,4 +1,5 @@ {-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-} module Numeric.LAPACK.Matrix.Hermitian.Basic (    Hermitian,    Transposition(..),@@ -7,6 +8,8 @@    identity,    diagonal,    takeDiagonal,+   forceOrder,+   stack,     multiplyVector,    square,@@ -21,11 +24,13 @@    ) where  import qualified Numeric.LAPACK.Matrix.Symmetric.Private as Symmetric+import qualified Numeric.LAPACK.Matrix.Triangular.Private as Triangular import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent+import qualified Numeric.LAPACK.Matrix.Basic as Basic import Numeric.LAPACK.Matrix.Hermitian.Private (Diagonal(..), TakeDiagonal(..)) import Numeric.LAPACK.Matrix.Triangular.Private-         (forPointers, pack, unpack, unpackToTemp,+         (forPointers, pack, packRect, unpack, unpackToTemp,           diagonalPointers, diagonalPointerPairs,           rowMajorPointers, columnMajorPointers) import Numeric.LAPACK.Matrix.Shape.Private@@ -47,8 +52,10 @@ import qualified Numeric.Netlib.Class as Class  import qualified Data.Array.Comfort.Storable.Unchecked as Array+import qualified Data.Array.Comfort.Storable as CheckedArray import qualified Data.Array.Comfort.Shape as Shape import Data.Array.Comfort.Storable.Unchecked (Array(Array))+import Data.Array.Comfort.Shape ((:+:))  import Foreign.C.Types (CInt, CChar) import Foreign.ForeignPtr (withForeignPtr)@@ -68,7 +75,7 @@  fromList :: (Shape.C sh, Storable a) => Order -> sh -> [a] -> Hermitian sh a fromList order sh =-   Array.fromList (MatrixShape.Hermitian order sh)+   CheckedArray.fromList (MatrixShape.Hermitian order sh)  autoFromList :: (Storable a) => Order -> [a] -> Hermitian ZeroInt a autoFromList order xs =@@ -123,6 +130,54 @@    withForeignPtr a $ \aPtr ->       forM_ (diagonalPointerPairs order n xPtr aPtr) $          \(dstPtr,srcPtr) -> poke dstPtr . realPart =<< peek srcPtr+++{-+This is not maximally efficient.+It fills up a whole square.+This wastes memory but enables more regular memory access patterns.+Additionally, it fills unused parts of the square with mirrored values.+-}+forceOrder ::+   (Shape.C sh, Class.Floating a) =>+   Order -> Hermitian sh a -> Hermitian sh a+forceOrder newOrder a =+   if MatrixShape.hermitianOrder (Array.shape a) == newOrder+      then a+      else fromUpperPart $ Basic.forceOrder newOrder $ toSquare a++fromUpperPart ::+   (Extent.C vert, Shape.C height, Shape.C width, Class.Floating a) =>+   Full vert Extent.Small height width a -> Hermitian width a+fromUpperPart (Array (MatrixShape.Full order extent) a) =+   let (height,width) = Extent.dimensions extent+       m = Shape.size height+       n = Shape.size width+       k = case order of RowMajor -> n; ColumnMajor -> m+   in Array.unsafeCreate (MatrixShape.Hermitian order width) $ \bPtr ->+      withForeignPtr a $ \aPtr -> packRect order n k aPtr bPtr+++{- |+> toSquare (stack a b c)+>+> =+>+> toSquare a ||| b+> ===+> adjoint b ||| toSquare c++It holds @order (stack a b c) = order b@.+The function is most efficient when the order of all blocks match.+-}+stack ::+   (Shape.C sh0, Eq sh0, Shape.C sh1, Eq sh1, Class.Floating a) =>+   Hermitian sh0 a -> General sh0 sh1 a -> Hermitian sh1 a ->+   Hermitian (sh0:+:sh1) a+stack a b c =+   let order = MatrixShape.fullOrder $ Array.shape b+   in Triangular.stack "Hermitian" (MatrixShape.Hermitian order)+         (forceOrder order a) b (forceOrder order c)   multiplyVector ::

@@ -0,0 +1,86 @@+{-# LANGUAGE TypeFamilies #-}+module Numeric.LAPACK.Matrix.Indexed where++import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape+import qualified Numeric.LAPACK.Matrix.Shape.Box as Box+import qualified Numeric.LAPACK.Matrix.Extent as Extent+import Numeric.LAPACK.Scalar (conjugate, zero)++import qualified Numeric.Netlib.Class as Class++import qualified Type.Data.Num.Unary as Unary++import qualified Data.Array.Comfort.Storable.Unchecked as UArray+import qualified Data.Array.Comfort.Shape as Shape+import Data.Array.Comfort.Storable (Array, (!))++import Data.Tuple.HT (swap)+import Data.Bool.HT (if')+++infixl 9 #!++class (Box.Box sh) => Indexed sh where+   (#!) ::+      (Class.Floating a) =>+      Array sh a ->+      (Shape.Index (Box.HeightOf sh), Shape.Index (Box.WidthOf sh)) -> a++instance+   (Extent.C vert, Extent.C horiz, Shape.Indexed height, Shape.Indexed width) =>+      Indexed (MatrixShape.Full vert horiz height width) where+   (#!) = (!)++instance (Shape.Indexed size) => Indexed (MatrixShape.Hermitian size) where+   arr#!ij =+      if Shape.inBounds (UArray.shape arr) ij+         then arr UArray.! ij+         else conjugate $ arr ! swap ij++instance+   (MatrixShape.Content lo, MatrixShape.TriDiag diag, MatrixShape.Content up,+    Shape.Indexed size) =>+      Indexed (MatrixShape.Triangular lo diag up size) where+   arr#!ij =+      let sh = UArray.shape arr+      in if Shape.inBounds sh ij+            then arr UArray.! ij+            else+               MatrixShape.caseDiagUpLoSym (MatrixShape.triangularUplo sh)+                  (checkedZero "Diagonal" sh ij)+                  (checkedZero "UpperTriangular" sh ij)+                  (checkedZero "LowerTriangular" sh ij)+                  (arr ! swap ij)++instance+   (Unary.Natural sub, Unary.Natural super,+    Extent.C vert, Extent.C horiz, Shape.Indexed height, Shape.Indexed width) =>+      Indexed (MatrixShape.Banded sub super vert horiz height width) where+   arr#!ij =+      let boxIx = uncurry MatrixShape.InsideBox+          sh = UArray.shape arr+      in if Shape.inBounds sh $ boxIx ij+            then arr UArray.! boxIx ij+            else checkedZero "Banded" sh ij++instance+   (Unary.Natural off, Shape.Indexed size) =>+      Indexed (MatrixShape.BandedHermitian off size) where+   arr#!ij =+      let boxIx = uncurry MatrixShape.InsideBox+          sh = UArray.shape arr+      in if' (Shape.inBounds sh $ boxIx ij)+            (arr UArray.! boxIx ij) $+         if' (Shape.inBounds sh $ boxIx $ swap ij)+            (conjugate $ arr UArray.! boxIx (swap ij))+         (checkedZero "BandedHermitian" sh ij)++checkedZero ::+   (Box.Box shape, Class.Floating a,+    Box.HeightOf shape ~ height, Shape.Indexed height,+    Box.WidthOf shape ~ width, Shape.Indexed width) =>+   String -> shape -> (Shape.Index height, Shape.Index width) -> a+checkedZero name sh ij =+   if Shape.inBounds (Box.height sh, Box.width sh) ij+      then zero+      else error $ "Matrix.Indexed." ++ name ++ ": index out of range"

@@ -7,6 +7,7 @@ module Numeric.LAPACK.Matrix.Multiply where  import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape+import qualified Numeric.LAPACK.Matrix.Shape.Box as Box import qualified Numeric.LAPACK.Matrix.Extent.Private as ExtentPriv import qualified Numeric.LAPACK.Matrix.Extent as Extent import qualified Numeric.LAPACK.Matrix.BandedHermitian.Basic as BandedHermitian@@ -15,7 +16,7 @@ import qualified Numeric.LAPACK.Matrix.Hermitian.Basic as Hermitian import qualified Numeric.LAPACK.Private as Private import Numeric.LAPACK.Matrix.Shape.Private-         (HeightOf, WidthOf, Empty, Filled, Unit, NonUnit,+         (Empty, Filled, Unit, NonUnit,           Order(RowMajor,ColumnMajor), flipOrder, transposeFromOrder) import Numeric.LAPACK.Matrix.Extent.Private (Small) import Numeric.LAPACK.Matrix.Triangular.Basic (Triangular)@@ -139,12 +140,16 @@ class (Shape.C shape) => MultiplyRight shape where    (#>) ::       (Class.Floating a) =>-      Array shape a -> Vector (WidthOf shape) a -> Vector (HeightOf shape) a+      Array shape a ->+      Vector (Box.WidthOf shape) a ->+      Vector (Box.HeightOf shape) a  class (Shape.C shape) => MultiplyLeft shape where    (<#) ::       (Class.Floating a) =>-      Vector (HeightOf shape) a -> Array shape a -> Vector (WidthOf shape) a+      Vector (Box.HeightOf shape) a ->+      Array shape a ->+      Vector (Box.WidthOf shape) a   instance

@@ -1,6 +1,7 @@ module Numeric.LAPACK.Matrix.Private where  import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape+import qualified Numeric.LAPACK.Matrix.Shape.Box as Box import qualified Numeric.LAPACK.Matrix.Extent as Extent import Numeric.LAPACK.Matrix.Shape.Private (Order, flipOrder) @@ -61,15 +62,11 @@ generalizeWide = mapExtent Extent.generalizeWide  -height ::-   (Extent.C vert, Extent.C horiz) =>-   Full vert horiz height width a -> height-height = MatrixShape.fullHeight . Array.shape+height :: (Box.Box shape) => Array shape a -> Box.HeightOf shape+height = Box.height . Array.shape -width ::-   (Extent.C vert, Extent.C horiz) =>-   Full vert horiz height width a -> width-width = MatrixShape.fullWidth . Array.shape+width :: (Box.Box shape) => Array shape a -> Box.WidthOf shape+width = Box.width . Array.shape   data Transposition = NonTransposed | Transposed

@@ -33,6 +33,9 @@    autoDiag,    autoUplo,    DiagUpLo,+   switchDiagUpLoSym,+   TriDiag,+   switchTriDiag,    Unit(Unit),    NonUnit(NonUnit), @@ -48,15 +51,16 @@    bandedFromFull,    UnaryProxy,    addOffDiagonals,-   TriDiag,-   switchTriDiag,    Content,     BandedHermitian(..),    bandedHermitian,++   Box.Box, Box.HeightOf, Box.WidthOf, Box.height, Box.width,    ) where  import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent+import qualified Numeric.LAPACK.Matrix.Shape.Box as Box import Numeric.LAPACK.Matrix.Shape.Private  

@@ -0,0 +1,18 @@+{-# LANGUAGE TypeFamilies #-}+module Numeric.LAPACK.Matrix.Shape.Box where++import qualified Data.Array.Comfort.Shape as Shape+++class (Shape.C shape) => Box shape where+   type HeightOf shape+   type WidthOf shape+   height :: shape -> HeightOf shape+   width :: shape -> WidthOf shape++indices ::+   (Box shape,+    HeightOf shape ~ height, Shape.Indexed height,+    WidthOf shape ~ width, Shape.Indexed width) =>+   shape -> [(Shape.Index height, Shape.Index width)]+indices sh = Shape.indices (height sh, width sh)

@@ -4,6 +4,7 @@ {-# LANGUAGE MultiParamTypeClasses #-} module Numeric.LAPACK.Matrix.Shape.Private where +import qualified Numeric.LAPACK.Matrix.Shape.Box as Box import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent import Numeric.LAPACK.Matrix.Extent.Private (Extent) import Numeric.LAPACK.Wrapper (Flip(Flip, getFlip))@@ -54,18 +55,19 @@    in (side,(m,n))  -type family HeightOf shape-type family WidthOf shape-- data Full vert horiz height width =    Full {       fullOrder :: Order,       fullExtent :: Extent vert horiz height width    } deriving (Eq, Show) -type instance HeightOf (Full vert horiz height width) = height-type instance WidthOf (Full vert horiz height width) = width+instance+   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width) =>+      Box.Box (Full vert horiz height width) where+   type HeightOf (Full vert horiz height width) = height+   type WidthOf (Full vert horiz height width) = width+   height = fullHeight+   width = fullWidth  instance    (Extent.C vert, Extent.C horiz, NFData height, NFData width) =>@@ -238,9 +240,13 @@    either (Left . Split lowerPart order) (Right . Split lowerPart order) $    Extent.caseTallWide (\h w -> Shape.size h >= Shape.size w) extent --type instance HeightOf (Split lower vert horiz height width) = height-type instance WidthOf (Split lower vert horiz height width) = width+instance+   (Eq lower, Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width) =>+      Box.Box (Split lower vert horiz height width) where+   type HeightOf (Split lower vert horiz height width) = height+   type WidthOf (Split lower vert horiz height width) = width+   height = splitHeight+   width = splitWidth  data Reflector = Reflector deriving (Eq, Show) data Triangle = Triangle deriving (Eq, Show)@@ -273,7 +279,8 @@       Shape.uncheckedSize (Extent.dimensions extent)  instance-   (Eq lower, Extent.C vert, Extent.C horiz, Shape.Indexed height, Shape.Indexed width) =>+   (Eq lower, Extent.C vert, Extent.C horiz,+    Shape.Indexed height, Shape.Indexed width) =>       Shape.Indexed (Split lower vert horiz height width) where     type Index (Split lower vert horiz height width) =@@ -340,8 +347,11 @@       hermitianSize :: size    } deriving (Eq, Show) -type instance HeightOf (Hermitian size) = size-type instance WidthOf (Hermitian size) = size+instance (Shape.C size) => Box.Box (Hermitian size) where+   type HeightOf (Hermitian size) = size+   type WidthOf (Hermitian size) = size+   height = hermitianSize+   width = hermitianSize  uploFromOrder :: Order -> Char uploFromOrder RowMajor = 'L'@@ -386,8 +396,13 @@       triangularSize :: size    } deriving (Eq, Show) -type instance HeightOf (Triangular lo diag up size) = size-type instance WidthOf (Triangular lo diag up size) = size+instance+   (Content lo, TriDiag diag, Content up, Shape.C size) =>+      Box.Box (Triangular lo diag up size) where+   type HeightOf (Triangular lo diag up size) = size+   type WidthOf (Triangular lo diag up size) = size+   height = triangularSize+   width = triangularSize   data Unit = Unit deriving (Eq, Show)@@ -709,8 +724,14 @@ bandedMapExtent f (Banded numDiag order extent) =    Banded numDiag order $ Extent.apply f extent -type instance HeightOf (Banded sub super vert horiz height width) = height-type instance WidthOf (Banded sub super vert horiz height width) = width+instance+   (Unary.Natural sub, Unary.Natural super, Extent.C vert, Extent.C horiz,+    Shape.C height, Shape.C width) =>+      Box.Box (Banded sub super vert horiz height width) where+   type HeightOf (Banded sub super vert horiz height width) = height+   type WidthOf (Banded sub super vert horiz height width) = width+   height = bandedHeight+   width = bandedWidth  bandedBreadth ::    (Unary.Natural sub, Unary.Natural super) =>@@ -929,8 +950,13 @@       bandedHermitianSize :: size    } deriving (Eq, Show) -type instance HeightOf (BandedHermitian off size) = size-type instance WidthOf (BandedHermitian off size) = size+instance+   (Unary.Natural off, Shape.C size) =>+      Box.Box (BandedHermitian off size) where+   type HeightOf (BandedHermitian off size) = size+   type WidthOf (BandedHermitian off size) = size+   height = bandedHermitianSize+   width = bandedHermitianSize  instance (Unary.Natural off, NFData size) =>       NFData (BandedHermitian off size) where

@@ -47,6 +47,7 @@ import qualified Numeric.Netlib.Class as Class  import qualified Data.Array.Comfort.Storable.Unchecked as Array+import qualified Data.Array.Comfort.Storable as CheckedArray import qualified Data.Array.Comfort.Shape as Shape import Data.Array.Comfort.Storable.Unchecked (Array(Array)) @@ -85,7 +86,7 @@  fromList :: (Shape.C sh, Storable a) => sh -> [a] -> Square sh a fromList sh =-   Array.fromList (MatrixShape.square RowMajor sh)+   CheckedArray.fromList (MatrixShape.square RowMajor sh)  autoFromList :: (Storable a) => [a] -> Square ZeroInt a autoFromList xs =

@@ -1,10 +1,12 @@ {-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-} {-# LANGUAGE ConstraintKinds #-} module Numeric.LAPACK.Matrix.Triangular.Basic (    Triangular, MatrixShape.UpLo,-   Upper, UnitUpper,-   Lower, UnitLower,-   Symmetric, Diagonal,+   Diagonal,+   Upper, FlexUpper, UnitUpper,+   Lower, FlexLower, UnitLower,+   Symmetric, FlexSymmetric,    fromList, autoFromList,    lowerFromList, autoLowerFromList,    upperFromList, autoUpperFromList,@@ -18,10 +20,13 @@    takeDiagonal,    transpose,    adjoint,+   stackLower,+   stackUpper,+   stackSymmetric,     toSquare,-   takeUpper,    takeLower,+   takeUpper,     fromLowerRowMajor, toLowerRowMajor,    fromUpperRowMajor, toUpperRowMajor,@@ -38,6 +43,7 @@  import qualified Numeric.LAPACK.Matrix.Symmetric.Private as Symmetric import qualified Numeric.LAPACK.Matrix.Triangular.Private as Tri+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@@ -62,8 +68,10 @@ import qualified Numeric.Netlib.Class as Class  import qualified Data.Array.Comfort.Storable.Unchecked as Array+import qualified Data.Array.Comfort.Storable as CheckedArray import qualified Data.Array.Comfort.Shape as Shape import Data.Array.Comfort.Storable.Unchecked (Array(Array))+import Data.Array.Comfort.Shape ((:+:))  import Foreign.C.Types (CChar, CInt) import Foreign.ForeignPtr (withForeignPtr)@@ -85,8 +93,8 @@ type FlexUpper diag sh = Array (MatrixShape.UpperTriangular diag sh) type FlexSymmetric diag sh = Array (MatrixShape.FlexSymmetric diag sh) -type UnitLower sh = Array (MatrixShape.LowerTriangular Unit sh)-type UnitUpper sh = Array (MatrixShape.UpperTriangular Unit sh)+type UnitLower sh = FlexLower Unit sh+type UnitUpper sh = FlexUpper Unit sh  transpose ::    (MatrixShape.Content lo, MatrixShape.Content up,@@ -106,7 +114,8 @@    (MatrixShape.Content lo, MatrixShape.Content up, Shape.C sh, Storable a) =>    Order -> sh -> [a] -> Triangular lo NonUnit up sh a fromList order sh =-   Array.fromList (MatrixShape.Triangular NonUnit MatrixShape.autoUplo order sh)+   CheckedArray.fromList+      (MatrixShape.Triangular NonUnit MatrixShape.autoUplo order sh)  lowerFromList :: (Shape.C sh, Storable a) => Order -> sh -> [a] -> Lower sh a lowerFromList = fromList@@ -372,6 +381,75 @@    (MatrixShape.TriDiag diag) =>    Triangular lo diag up sh a -> Triangular lo NonUnit up sh a strictNonUnitDiagonal = Array.mapShape MatrixShape.strictNonUnitDiagonal+++{-+It does not make much sense to put+'stackLower', 'stackUpper', 'stackSymmetric' in one function+because in 'stackLower' and 'stackUpper'+the height and width of matrix b is swapped.+-}+stackLower ::+   (MatrixShape.TriDiag diag,+    Shape.C sh0, Eq sh0, Shape.C sh1, Eq sh1, Class.Floating a) =>+   FlexLower diag sh0 a ->+   Matrix.General sh1 sh0 a ->+   FlexLower diag sh1 a ->+   FlexLower diag (sh0:+:sh1) a+stackLower a b c =+   transpose $+   stackAux "LowerTriangular" (transpose a) (Basic.transpose b) (transpose c)++{- |+For upper triangular matrices it holds:++> toSquare (stack a b c)+>+> =+>+> toSquare a ||| b+> ===+> constant 0 ||| toSquare c++For the other triangular flavors it holds accordingly++It holds @order (stack a b c) = order b@.+The function is most efficient when the order of all blocks match.+-}+stackUpper ::+   (MatrixShape.TriDiag diag,+    Shape.C sh0, Eq sh0, Shape.C sh1, Eq sh1, Class.Floating a) =>+   FlexUpper diag sh0 a ->+   Matrix.General sh0 sh1 a ->+   FlexUpper diag sh1 a ->+   FlexUpper diag (sh0:+:sh1) a+stackUpper = stackAux "UpperTriangular"++stackSymmetric ::+   (MatrixShape.TriDiag diag,+    Shape.C sh0, Eq sh0, Shape.C sh1, Eq sh1, Class.Floating a) =>+   FlexSymmetric diag sh0 a ->+   Matrix.General sh0 sh1 a ->+   FlexSymmetric diag sh1 a ->+   FlexSymmetric diag (sh0:+:sh1) a+stackSymmetric = stackAux "Symmetric"++stackAux ::+   (MatrixShape.Content lo, MatrixShape.TriDiag diag,+    Shape.C sh0, Eq sh0, Shape.C sh1, Eq sh1, Class.Floating a) =>+   String ->+   Triangular lo diag MatrixShape.Filled sh0 a ->+   Matrix.General sh0 sh1 a ->+   Triangular lo diag MatrixShape.Filled sh1 a ->+   Triangular lo diag MatrixShape.Filled (sh0:+:sh1) a+stackAux name a b c =+   let order = MatrixShape.fullOrder $ Array.shape b+   in Tri.stack name+         (\sh ->+            (Array.shape a) {+               MatrixShape.triangularOrder = order,+               MatrixShape.triangularSize = sh})+         (forceOrder order a) b (forceOrder order c)   multiplyVector ::

@@ -1,7 +1,10 @@ {-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-} module Numeric.LAPACK.Matrix.Triangular.Private where +import qualified Numeric.LAPACK.Matrix.Private as Matrix import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape+import qualified Numeric.LAPACK.Matrix.Shape.Box as Box import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent import Numeric.LAPACK.Matrix.Shape.Private          (Order(RowMajor,ColumnMajor), flipOrder, uploFromOrder,@@ -19,6 +22,7 @@ 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 Data.Array.Comfort.Shape ((:+:)((:+:)))  import Foreign.Marshal.Alloc (alloca) import Foreign.Marshal.Array (advancePtr)@@ -128,6 +132,77 @@    nPtr <- Call.cint n    zPtr <- Call.number z    liftIO $ LapackGen.laset uploPtr nPtr nPtr zPtr zPtr aPtr nPtr++++stack ::+   (Box.Box sh0, Box.HeightOf sh0 ~ height, Shape.C height, Eq height,+    Box.Box sh1, Box.WidthOf sh1 ~ width, Shape.C width, Eq width,+    Shape.C sh2, Class.Floating a) =>+   String -> (height:+:width -> sh2) ->+   Array sh0 a -> Matrix.General height width a -> Array sh1 a -> Array sh2 a+stack name consShape+      (Array sha a) (Array (MatrixShape.Full order extent) b) (Array shc c) =+   let (height,width) = Extent.dimensions extent+   in Array.unsafeCreate (consShape (height :+: width)) $ \xPtr ->+   withForeignPtr b $ \bPtr -> do+      Call.assert (name++".stack: height shapes mismatch") $+         height == Box.height sha+      Call.assert (name++".stack: width shapes mismatch") $+         width == Box.width shc+      let m = Shape.size height+      let n = Shape.size width+      copyTriangleA order m n a xPtr+      copyRectangle order m n bPtr xPtr+      copyTriangleC order m n c xPtr++copyTriangleA ::+   (Class.Floating a) =>+   Order -> Int -> Int -> ForeignPtr a -> Ptr a -> IO ()+copyTriangleA order m n a xPtr =+   withForeignPtr a $ \aPtr ->+      case order of+         ColumnMajor -> copyBlock (Shape.triangleSize m) aPtr xPtr+         RowMajor ->+            forM_ (zip (iterate pred m) $+                   zip (diagonalPointers order m aPtr)+                       (diagonalPointers order (m+n) xPtr)) $+               \(k,(aiPtr,xiPtr)) -> copyBlock k aiPtr xiPtr++copyTriangleC ::+   (Class.Floating a) =>+   Order -> Int -> Int -> ForeignPtr a -> Ptr a -> IO ()+copyTriangleC order m n c xPtr =+   withForeignPtr c $ \cPtr ->+      case order of+         RowMajor ->+            let triSize = Shape.triangleSize n+            in copyBlock triSize cPtr+                  (advancePtr xPtr $ Shape.triangleSize (m+n) - triSize)+         ColumnMajor ->+            forM_ (zip (iterate succ 0) $+                   zip (diagonalPointers order n cPtr)+                       (drop m $ diagonalPointers order (m+n) xPtr)) $+               \(k,(aiPtr,xiPtr)) ->+                  copyBlock (k+1)+                     (advancePtr aiPtr (-k)) (advancePtr xiPtr (-k))++copyRectangle ::+   (Class.Floating a) => Order -> Int -> Int -> Ptr a -> Ptr a -> IO ()+copyRectangle order m n bPtr xPtr =+   case order of+      RowMajor ->+         forM_ (take m $ zip (iterate pred m) $+                zip (pointerSeq n bPtr) (diagonalPointers order (m+n) xPtr)) $+            \(k,(biPtr,xiPtr)) ->+               copyBlock n biPtr (advancePtr xiPtr k)+      ColumnMajor ->+         forM_ (take n $ zip (iterate succ m) $+                zip (pointerSeq m bPtr)+                    (drop m $ diagonalPointers order (m+n) xPtr)) $+            \(k,(biPtr,xiPtr)) ->+               copyBlock m biPtr (advancePtr xiPtr (-k))+   type Triangular lo diag up sh = Array (MatrixShape.Triangular lo diag up sh)

@@ -16,7 +16,7 @@  import qualified Numeric.LAPACK.Matrix.Square.Basic as Square import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent-import qualified Numeric.LAPACK.Matrix.Extent.Kind as EK+import qualified Numeric.LAPACK.Matrix.Basic as Basic import Numeric.LAPACK.Matrix.Square.Basic (Square) import Numeric.LAPACK.Matrix.Private (Full, Tall, ZeroInt, zeroInt) import Numeric.LAPACK.Matrix (transpose, dropColumns)@@ -160,9 +160,9 @@   {- |-If @x = leastSquaresMinimumNormRCond rcond a b@+If @(rank,x) = leastSquaresMinimumNormRCond rcond a b@ then @x@ is the vector with minimum @Vector.norm2 x@-that minimizes @Vector.norm2 (multiply a x `sub` b)@.+that minimizes @Vector.norm2 (a#>x `sub` b)@.  Matrix @a@ can have any rank but you must specify the reciprocal condition of the rank-truncated matrix.@@ -334,12 +334,5 @@    Tall height width a -> Tall height ZeroInt a complement a =    dropColumns (Shape.size $ MatrixShape.fullWidth $ Array.shape a) $-   Array.mapShape zeroIntWidth $ Square.toFull $-   HH.extractQ $ HH.fromMatrix $ Matrix.fromFull a--zeroIntWidth ::-   (Shape.C width) =>-   MatrixShape.Tall height width -> MatrixShape.Tall height ZeroInt-zeroIntWidth (MatrixShape.Full order (Extent.Extent o (EK.Tall height width))) =-   MatrixShape.Full order-      (Extent.Extent o (EK.Tall height (zeroInt $ Shape.size width)))+   Basic.mapWidth (zeroInt . Shape.size) $ Square.toFull $+   HH.extractQ $ HH.fromMatrix a

@@ -9,6 +9,8 @@    autoFromList,    append, take, drop,    takeLeft, takeRight,+   swap,+   singleton,    constant,    unit,    dot, inner,@@ -56,16 +58,20 @@ import Foreign.Marshal.Array (copyArray, advancePtr) import Foreign.ForeignPtr (withForeignPtr) import Foreign.Ptr (Ptr, castPtr)-import Foreign.Storable (Storable, peek, peekElemOff, pokeElemOff)+import Foreign.Storable (Storable, peek, poke, peekElemOff, pokeElemOff) import Foreign.C.Types (CInt)  import System.IO.Unsafe (unsafePerformIO)  import Control.Monad.Trans.Cont (ContT(ContT), evalContT) import Control.Monad.IO.Class (liftIO)+import Control.Monad.ST (runST) import Control.Applicative (Const(Const,getConst), liftA3, (<$>)) +import qualified Data.Array.Comfort.Storable.Mutable.Unchecked as UMutArray+import qualified Data.Array.Comfort.Storable.Mutable as MutArray import qualified Data.Array.Comfort.Storable.Unchecked as Array+import qualified Data.Array.Comfort.Storable as CheckedArray import qualified Data.Array.Comfort.Shape as Shape import Data.Array.Comfort.Storable.Unchecked (Array(Array)) import Data.Array.Comfort.Shape ((:+:)((:+:)))@@ -85,12 +91,20 @@ toList = Array.toList  fromList :: (Shape.C sh, Storable a) => sh -> [a] -> Vector sh a-fromList = Array.fromList+fromList = CheckedArray.fromList  autoFromList :: (Storable a) => [a] -> Vector (Shape.ZeroBased Int) a autoFromList = Array.vectorFromList  +{- |+> singleton = constant ()++However, singleton does not need 'Class.Floating' constraint.+-}+singleton :: (Storable a) => a -> Vector () a+singleton a = Array.unsafeCreate () $ flip poke a+ constant :: (Shape.C sh, Class.Floating a) => sh -> a -> Vector sh a constant sh a = Array.unsafeCreateWithSize sh $ fill a @@ -144,6 +158,18 @@    Array.unsafeCreateWithSize sh1 $ \k yPtr ->    withForeignPtr x $ \xPtr ->       copyArray yPtr (advancePtr xPtr (Shape.size sh0)) k++swap ::+   (Shape.Indexed sh, Storable a) =>+   Shape.Index sh -> Shape.Index sh -> Vector sh a -> Vector sh a+swap i j x =+   runST (do+      y <- MutArray.thaw x+      xi <- MutArray.read y i+      xj <- MutArray.read y j+      MutArray.write y i xj+      MutArray.write y j xi+      UMutArray.unsafeFreeze y)   newtype Dot sh a = Dot {runDot :: Vector sh a -> Vector sh a -> a}

@@ -4,6 +4,7 @@ {-# LANGUAGE GADTs #-} module Test.Banded (testsVar) where +import qualified Test.Indexed as Indexed import qualified Test.Generator as Gen import qualified Test.Utility as Util import Test.Banded.Utility@@ -15,6 +16,7 @@           genOrder, genArray, Tagged, equalListWith)  import qualified Numeric.LAPACK.Matrix.Banded as Banded+import qualified Numeric.LAPACK.Matrix.Shape.Box as Box import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape import qualified Numeric.LAPACK.Matrix.Square as Square import qualified Numeric.LAPACK.Matrix as Matrix@@ -30,6 +32,7 @@ import qualified Type.Data.Num.Unary as Unary import Type.Data.Num.Unary (unary, (:+:)) +import qualified Data.Array.Comfort.Storable as Array import qualified Data.Array.Comfort.Shape as Shape  import Foreign.Storable (Storable)@@ -230,6 +233,13 @@    (Show a, Class.Floating a, Eq a, RealOf a ~ ar, Class.Real ar) =>    [(String, Tagged a QC.Property)] testsVar =+   ("index",+      checkForAll+         (Indexed.genMatrixIndexGen+            (\(Banded arr) -> Box.indices $ Array.shape arr)+            genBanded)+         (\(mix, Banded arr) -> Indexed.unitDot (mix, arr))) :+    ("multiplyFullIdentity",       checkForAll genBanded multiplyFullIdentity) :    ("multiplyFullAny",

@@ -5,6 +5,7 @@ {-# LANGUAGE GADTs #-} module Test.BandedHermitian (testsVar) where +import qualified Test.Indexed as Indexed import qualified Test.Generator as Gen import qualified Test.Utility as Util import Test.Banded.Utility@@ -19,6 +20,7 @@ import qualified Numeric.LAPACK.Matrix.BandedHermitian as BandedHermitian import qualified Numeric.LAPACK.Matrix.Banded as Banded import qualified Numeric.LAPACK.Matrix.Hermitian as Hermitian+import qualified Numeric.LAPACK.Matrix.Shape.Box as Box import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape import qualified Numeric.LAPACK.Matrix.Square as Square import qualified Numeric.LAPACK.Matrix as Matrix@@ -42,7 +44,6 @@ import Control.Applicative (liftA2, (<$>))  import qualified Data.List.HT as ListHT-import Data.Traversable (for) import Data.Tuple.HT (mapSnd)  import qualified Test.QuickCheck as QC@@ -81,8 +82,8 @@    k <- QC.choose (0, toInteger maxDim)    Unary.reifyNatural k $ \numOff -> do       let shape = MatrixShape.bandedHermitian (unary numOff) order size-      BandedHermitian . Array.fromList shape <$>-         (for (Shape.indices shape) $ \ix ->+      BandedHermitian <$>+         (Util.genArrayIndexed shape $ \ix ->             let real =                   case ix of                      MatrixShape.InsideBox r c -> r==c@@ -359,6 +360,12 @@     Class.Floating a, Eq a, RealOf a ~ ar, Class.Real ar) =>    [(String, Tagged a QC.Property)] testsVar =+   ("index",+      checkForAll+         (Indexed.genMatrixIndexGen+            (\(BandedHermitian arr) -> Box.indices $ Array.shape arr)+            genBandedHermitian)+         (\(mix, BandedHermitian arr) -> Indexed.unitDot (mix, arr))) :    ("convertToFull",       checkForAll genBandedHermitian convertToFull) :    ("takeDiagonal",

@@ -6,7 +6,7 @@  import qualified Test.Logic as Logic import qualified Test.Utility as Util-import Test.Logic (Dim, MatchMode(DontForceMatch,ForceMatch))+import Test.Logic (Dim, MatchMode(DontForceMatch,ForceMatch), (=!=), (<!=)) import Test.Utility (Match)  import qualified UniqueLogic.ST.TF.System.Simple as Sys@@ -22,17 +22,20 @@  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.Storable (Array) import Data.Array.Comfort.Shape ((:+:))  import qualified Control.Monad.Trans.RWS as MRWS+import qualified Control.Monad.Trans.Class as MT import qualified Control.Applicative.HT as AppHT+import Control.Monad.Trans.RWS (RWST, evalRWST) import Control.Applicative (liftA2, (<*>), (<$>)) -import Data.Traversable (for)-import Data.Tuple.HT (swap, mapFst3, mapThd3)+import qualified Data.Ref as Ref+import Data.Semigroup ((<>))+import Data.Monoid (Monoid, mempty)+import Data.Tuple.HT (swap, mapFst, mapSnd)  import qualified Test.QuickCheck as QC @@ -44,7 +47,7 @@ Dimensions will be choosen arbitrarily from the range @(0,maxDim)@. Elements are choosen from the range @(-maxElem,maxElem)@. -I moved to the 's' tag to within the 'Cons' constructor+I moved the 's' tag to within the 'Cons' constructor and furthermore defined 'TaggedVariables' to strip the 's' tag from the Variables in 'dim'. This way, we can easily define 'checkForAll' in the test modules.@@ -54,8 +57,8 @@ newtype T tag dim array =    Cons {       decons :: forall s.-         Integer ->-         Logic.M s (TaggedVariables s dim, Logic.System s, Logic.M s array)+         RWST Integer (Logic.System s) () (Logic.M s)+            (TaggedVariables s dim, Logic.M s array)    }  data Variable dim@@ -67,8 +70,11 @@                   (TaggedVariables s a, TaggedVariables s b)  instance Functor (T tag dim) where-   fmap f (Cons gen) = Cons $ \mapElem -> mapThd3 (fmap f) <$> gen mapElem+   fmap f (Cons gen) = Cons $ mapSnd (fmap f) <$> gen +newVariable :: (Ref.C m, Monoid w) => RWST r w () m (Sys.Variable m a)+newVariable = MT.lift Sys.globalVariable+ run ::    T tag dim array -> Integer -> Int ->    Util.TaggedGen tag (array, Match)@@ -76,8 +82,8 @@    Util.Tagged $       QC.elements [DontForceMatch, ForceMatch] >>=       Logic.runSTInGen-         (do (_dim, sys, queries) <- decons gen maxElem-             Sys.solve sys+         (do ((_dim, queries), sys) <- evalRWST (decons gen) maxElem ()+             Logic.solve sys              queries)          maxDim @@ -92,16 +98,18 @@    (Integer -> a -> QC.Gen b) ->    T tag dim a -> T tag dim b mapGen f (Cons gen) =-   Cons $ \maxElem ->-      mapThd3 (Logic.liftGen . f maxElem =<<) <$> gen maxElem+   Cons $ do+      maxElem <- MRWS.ask+      mapSnd (Logic.liftGen . f maxElem =<<) <$> gen  mapGenDim ::    (Integer -> Int -> a -> QC.Gen b) ->    T tag dim a -> T tag dim b mapGenDim f (Cons gen) =-   Cons $ \maxElem -> do-      (maxDim, _matchMode) <- Logic.M MRWS.ask-      mapThd3 (Logic.liftGen . f maxElem maxDim =<<) <$> gen maxElem+   Cons $ do+      maxElem <- MRWS.ask+      (maxDim, _matchMode) <- MT.lift $ Logic.M MRWS.ask+      mapSnd (Logic.liftGen . f maxElem maxDim =<<) <$> gen   combine ::@@ -112,26 +120,28 @@    T tag dimA a ->    T tag dimB b combine combineDim (Cons genF) (Cons genA) =-   Cons $ \maxElem -> do-      (dimF,sysF,f) <- genF maxElem-      (dimA,sysA,a) <- genA maxElem+   Cons $ do+      (dimF,f) <- genF+      (dimA,a) <- genA       let (dimB, constraint) = combineDim dimF dimA-      return (dimB, sysF >> sysA >> constraint, f <*> a)+      MRWS.tell constraint+      return (dimB, f <*> a)   type Scalar tag = T tag ()  scalar :: (Class.Floating a) => Scalar a a scalar =-   Cons $ \ maxElem ->-      return ((), return (), Logic.liftGen $ Util.genElement maxElem)+   Cons $ do+      maxElem <- MRWS.ask+      return ((), Logic.liftGen $ Util.genElement maxElem)  (<.*.>) ::    (Dim size, Eq size) =>    Vector tag size (a -> b) ->    Vector tag size a ->    Scalar tag b-(<.*.>) = combine (\dimF dimA -> ((), Logic.ruleEqualDim dimF dimA))+(<.*.>) = combine (\dimF dimA -> ((), dimF=!=dimA))   queryZeroInt :: Logic.Variable s Int -> Logic.M s ZeroInt@@ -141,9 +151,9 @@  vectorDim :: (Class.Floating a) => Vector a Int ZeroInt vectorDim =-   Cons $ \ _maxElem -> do-      dim <- Sys.globalVariable-      return (dim, return (), queryZeroInt dim)+   Cons $ do+      dim <- newVariable+      return (dim, queryZeroInt dim)  vector :: (Class.Floating a) => Vector a Int (Vector.Vector ZeroInt a) vector = mapGen Util.genArray vectorDim@@ -158,24 +168,21 @@    Vector tag height (a -> b) ->    Matrix tag height width a ->    Vector tag width b-(<.*|>) =-   combine (\dim (height,width) -> (width, Logic.ruleEqualDim dim height))+(<.*|>) = combine (\dim (height,width) -> (width, dim=!=height))  (<|*.>) ::    (Dim width, Eq width) =>    Matrix tag height width (a -> b) ->    Vector tag width a ->    Vector tag height b-(<|*.>) =-   combine (\(height,width) dim -> (height, Logic.ruleEqualDim width dim))+(<|*.>) = combine (\(height,width) dim -> (height, width=!=dim))  (<.=.>) ::    (Dim size, Eq size) =>    Vector tag size (a -> b) ->    Vector tag size a ->    Vector tag size b-(<.=.>) =-   combine (\sizeF sizeA -> (sizeF, Logic.ruleEqualDim sizeF sizeA))+(<.=.>) = combine (\sizeF sizeA -> (sizeF, sizeF=!=sizeA))   type Matrix tag height width = T tag (Variable height, Variable width)@@ -183,9 +190,9 @@ matrixDims ::    (Class.Floating a) => Matrix a Int Int (ZeroInt, ZeroInt) matrixDims =-   Cons $ \ _maxElem -> do-      dims <- liftA2 (,) Sys.globalVariable Sys.globalVariable-      return (dims, return (), AppHT.mapPair (queryZeroInt,queryZeroInt) dims)+   Cons $ do+      dims <- liftA2 (,) newVariable newVariable+      return (dims, AppHT.mapPair (queryZeroInt,queryZeroInt) dims)  matrix ::    (Class.Floating a) => Matrix a Int Int (Matrix.General ZeroInt ZeroInt a)@@ -197,9 +204,9 @@  squareDim :: (Class.Floating a) => Matrix a Int Int ZeroInt squareDim =-   Cons $ \ _maxElem -> do-      dim <- Sys.globalVariable-      return ((dim,dim), return (), queryZeroInt dim)+   Cons $ do+      dim <- newVariable+      return ((dim,dim), queryZeroInt dim)  squareShaped ::    (Shape.C sh, Class.Floating a) =>@@ -278,20 +285,17 @@    runGenTriangularDiag $    MatrixShape.switchTriDiag       (GenTriangularDiag $ \shape ->-         Array.fromList shape <$>-            (for (Shape.indices shape) $ \(r,c) ->-               if r==c-                  then return one-                  else Util.genElement maxElem))+         Util.genArrayExtraDiag maxElem shape (const $ return one))       (GenTriangularDiag $ Util.genArray maxElem)   tallDims :: (Class.Floating a) => Matrix a Int Int (ZeroInt, ZeroInt) tallDims =-   Cons $ \ _maxElem -> do-      height <- Sys.globalVariable-      width  <- Sys.globalVariable-      return ((height,width), Logic.ruleLessOrEqual width height,+   Cons $ do+      height <- newVariable+      width  <- newVariable+      MRWS.tell $  width <!= height+      return ((height,width),               liftA2 (,) (queryZeroInt height) (queryZeroInt width))  tall ::@@ -336,11 +340,8 @@    flip mapGen squareDim $ \maxElem size -> do       order <- Util.genOrder       let shape = MatrixShape.hermitian order size-      (Array.fromList shape <$>-         (for (Shape.indices shape) $ \(r,c) ->-            if r==c-               then fromReal <$> Util.genReal maxElem-               else Util.genElement maxElem))+      (Util.genArrayExtraDiag maxElem shape+          (const $ fromReal <$> Util.genReal maxElem))          `QC.suchThat` cond  @@ -353,14 +354,12 @@    Matrix tag fuse width a ->    Matrix tag height width b (<|*|>) =-   combine-      (\(height,fuseF) (fuseA,width) ->-         ((height,width), Logic.ruleEqualDim fuseF fuseA))+   combine (\(height,fuseF) (fuseA,width) -> ((height,width), fuseF=!=fuseA))  transpose ::    Matrix tag height width a ->    Matrix tag width height a-transpose (Cons gen) = Cons $ \maxElem -> mapFst3 swap <$> gen maxElem+transpose (Cons gen) = Cons $ mapFst swap <$> gen  (<|\|>) ::    (Dim height, Eq height) =>@@ -373,7 +372,7 @@    Vector tag height (a -> b) ->    Vector tag width a ->    Matrix tag height width b-(<***>) = combine (\height width -> ((height,width), return ()))+(<***>) = combine (\height width -> ((height,width), mempty))   (<|=|>) ::@@ -384,25 +383,30 @@    Matrix tag height width b (<|=|>) =    combine $ \(heightF,widthF) (heightA,widthA) ->-      ((heightF,widthF),-       Logic.ruleEqualDim heightF heightA >> Logic.ruleEqualDim widthF widthA)+      ((heightF,widthF), heightF=!=heightA <> widthF=!=widthA) ++(!+!) ::+   Logic.Variable s dimA ->+   Logic.Variable s dimB ->+   RWST r (Logic.System s) () (Logic.M s) (Logic.Variable s (dimA :+: dimB))+a!+!b = do+   c <- newVariable+   MRWS.tell $ (a Logic.!+! b) c+   return c+ (<===>) ::    (Dim width, Eq width) =>    Matrix tag heightA width (a -> b) ->    Matrix tag heightB width a ->    Matrix tag (heightA:+:heightB) width b (<===>) (Cons genF) (Cons genA) =-   Cons $ \maxElem -> do-      heightB <- Sys.globalVariable-      ((heightF,widthF),sysF,f) <- genF maxElem-      ((heightA,widthA),sysA,a) <- genA maxElem-      return-         ((heightB,widthF),-          do sysF >> sysA-             Logic.ruleEqualDim widthF widthA-             Logic.ruleAppendDim heightF heightA heightB,-          f <*> a)+   Cons $ do+      ((heightF,widthF),f) <- genF+      ((heightA,widthA),a) <- genA+      MRWS.tell $ widthF=!=widthA+      heightB <- heightF!+!heightA+      return ((heightB,widthF), f <*> a)  (<|||>) ::    (Dim height, Eq height) =>@@ -410,6 +414,25 @@    Matrix tag height widthB a ->    Matrix tag height (widthA:+:widthB) b (<|||>) f a = transpose $ transpose f <===> transpose a++++stack3 ::+   (Dim heightA, Eq heightA) =>+   (Dim widthB, Eq widthB) =>+   Matrix tag heightA widthA a ->+   Matrix tag heightA widthB b ->+   Matrix tag heightB widthB c ->+   Matrix tag (heightA:+:heightB) (widthA:+:widthB) (a,b,c)+stack3 (Cons genF) (Cons genA) (Cons genB) =+   Cons $ do+      ((heightF,widthF),f) <- genF+      ((heightA,widthA),a) <- genA+      ((heightB,widthB),b) <- genB+      MRWS.tell $  heightF=!=heightA  <>  widthA=!=widthB+      heightC <- heightF!+!heightB+      widthC <- widthF!+!widthB+      return ((heightC,widthC), (,,) <$> f <*> a <*> b)   infixl 4 <.*.>, <.*|>, <|*.>, <|*|>, <|\|>, <***>, <.=.>, <|=|>, <===>, <|||>

@@ -1,12 +1,13 @@ {-# LANGUAGE TypeFamilies #-} module Test.Hermitian (testsVar) where +import qualified Test.Indexed as Indexed import qualified Test.Generator as Gen import qualified Test.Utility as Util import Test.Generator ((<.*|>), (<|*.>), (<|*|>), (<|\|>)) import Test.Utility          (approx, approxReal, approxArray, approxArrayTol, approxMatrix,-          Tagged, genOrder)+          equalArray, Tagged, genOrder)  import qualified Numeric.LAPACK.Orthogonal.Householder as HH import qualified Numeric.LAPACK.Matrix.HermitianPositiveDefinite as HermitianPD@@ -19,7 +20,8 @@ import Numeric.LAPACK.Matrix.Hermitian (Hermitian) import Numeric.LAPACK.Matrix.Square (Square) import Numeric.LAPACK.Matrix.Shape (Order)-import Numeric.LAPACK.Matrix (General, ZeroInt, zeroInt, (<#), (<#>), (#>))+import Numeric.LAPACK.Matrix+         (General, ZeroInt, zeroInt, (<#), (<#>), (#>), (|||), (===)) import Numeric.LAPACK.Vector (Vector) import Numeric.LAPACK.Scalar (RealOf, fromReal, selectReal) @@ -36,6 +38,17 @@ import qualified Test.QuickCheck as QC  +stack ::+   (Class.Floating a) =>+   (Hermitian ZeroInt a, General ZeroInt ZeroInt a, Hermitian ZeroInt a) ->+   Bool+stack (a,b,c) =+   let abc = Matrix.fromFull $ Hermitian.toSquare $ Hermitian.stack a b c+   in equalArray abc $ Matrix.adaptOrder abc $+         (Matrix.fromFull (Hermitian.toSquare a) ||| b+          ===+          Matrix.adjoint b ||| Matrix.fromFull (Hermitian.toSquare c))+ covariance ::    (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>    General ZeroInt ZeroInt a -> Bool@@ -316,6 +329,10 @@    (Show a, Show ar, Class.Floating a, Eq a, RealOf a ~ ar, Class.Real ar) =>    [(String, Tagged a QC.Property)] testsVar =+   ("index",+      checkForAll (Indexed.genMatrixIndex Gen.hermitian) Indexed.unitDot) :+   ("stack",+      checkForAll (Gen.stack3 Gen.hermitian Gen.matrix Gen.hermitian) stack) :    ("covariance",       checkForAll Gen.matrix covariance) :    ("outer",

@@ -0,0 +1,53 @@+{-# LANGUAGE TypeFamilies #-}+module Test.Indexed where++import qualified Test.Generator as Gen+import Test.Utility (maybeProperty)++import qualified Numeric.LAPACK.Matrix.Shape.Box as Box+import qualified Numeric.LAPACK.Matrix as Matrix+import qualified Numeric.LAPACK.Vector as Vector+import Numeric.LAPACK.Matrix ((#!), (#>))++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.Storable (Array)++import qualified Data.Traversable as Trav+import Data.Maybe.HT (toMaybe)++import qualified Test.QuickCheck as QC+++genMatrixIndexGen ::+   (Class.Floating a) =>+   (array -> [ix]) ->+   Gen.Matrix a Int Int array ->+   Gen.Matrix a Int Int (Maybe ix, array)+genMatrixIndexGen indices gen =+   flip Gen.mapGen gen $ \_maxElem m -> do+      let set = indices m+      ij <- Trav.mapM QC.elements $ toMaybe (not $ null set) set+      return (ij,m)++genMatrixIndex ::+   (Shape.Indexed shape, Class.Floating a) =>+   Gen.Matrix a Int Int (Array shape a) ->+   Gen.Matrix a Int Int (Maybe (Shape.Index shape), Array shape a)+genMatrixIndex = genMatrixIndexGen (Shape.indices . Array.shape)++unitDot ::+   (Matrix.Indexed shape, Matrix.MultiplyRight shape,+    Box.HeightOf shape ~ height, Shape.Indexed height, Eq height,+    Box.WidthOf shape ~ width, Shape.Indexed width,+    Class.Floating a, Eq a) =>+   (Maybe (Shape.Index height, Shape.Index width), Array shape a) -> QC.Property+unitDot (mij,m) =+   maybeProperty $+   flip fmap mij $ \(i,j) ->+      m#!(i,j) ==+      Vector.dot+         (Vector.unit (Matrix.height m) i)+         (m #> Vector.unit (Matrix.width m) j)

@@ -19,6 +19,9 @@ import Control.Monad.ST (ST, runST) import Control.Applicative (Applicative, liftA2, (<$>)) +import qualified Data.Monoid.Applicative as AppMn+import Data.Semigroup ((<>))+ import qualified Test.QuickCheck as QC import qualified Test.QuickCheck.GenT as GenT import Test.QuickCheck.Gen (Gen(MkGen))@@ -48,7 +51,7 @@   type Variable s a = Sys.Variable (M s) a-type System s = Sys.T (M s) ()+type System s = AppMn.T (Sys.T (M s)) ()   example :: Int -> MatchMode -> QC.Gen ([Int], Match)@@ -60,18 +63,18 @@       d <- Sys.globalVariable       e <- Sys.globalVariable       f <- Sys.globalVariable-      Sys.solve $ do-         ruleLessOrEqual a b-         ruleEqualDim d f-         ruleEqualDim c d+      solve $  a<!=b  <>  d=!=f  <>  c=!=d       mapM query [a,b,c,d,e,f])  +solve :: System s -> M s ()+solve = Sys.solve . AppMn.run+ choose :: (Random a) => (Int -> (a,a)) -> M s a choose f = liftGen . QC.choose . f . fst =<< M MRWS.ask -ruleLessOrEqual :: Variable s Int -> Variable s Int -> System s-ruleLessOrEqual va vb = do+(<!=) :: Variable s Int -> Variable s Int -> System s+va <!= vb  =  AppMn.Cons $ do    assignmentM (\a -> choose (\ maxk -> (a,maxk))) va vb    assignmentM (\b -> choose (\_maxk -> (0,b))) vb va @@ -81,9 +84,9 @@ instance (Dim dimA, Dim dimB) => Dim (dimA:+:dimB) where    chooseDim = liftA2 (:+:) chooseDim chooseDim -ruleEqualDim ::+(=!=) ::    (Dim dim, Eq dim) => Variable s dim -> Variable s dim -> System s-ruleEqualDim va vb = do+va =!= vb  =  AppMn.Cons $ do    let equalM =          assignmentM $ \x -> do             matchMode <- M $ MRWS.asks snd@@ -101,9 +104,9 @@ assignmentM f vx vy = Sys.runApplyM (f <$> Sys.arg vx) vy  -ruleAppendDim ::+(!+!) ::    Variable s dimA -> Variable s dimB -> Variable s (dimA:+:dimB) -> System s-ruleAppendDim va vb vab = do+(!+!) va vb vab = AppMn.Cons $ do    Sys.assignment3 (:+:) va vb vab    Sys.assignment2 (\(a:+:_) -> a) vab va    Sys.assignment2 (\(_:+:b) -> b) vab vb

@@ -2,13 +2,14 @@ {-# LANGUAGE TypeOperators #-} module Test.Matrix (testsVar) where +import qualified Test.Indexed as Indexed import qualified Test.Generator as Gen import qualified Test.Utility as Util import Test.Generator          ((<|*|>), (<|*.>), (<.*.>), (<***>), (<|=|>), (<|||>), (<===>)) import Test.Utility-         (approx, approxArray, approxMatrix,-          genOrder, Tagged(Tagged), TaggedGen)+         (equalArray, approx, approxArray, approxMatrix,+          maybeProperty, genOrder, Tagged(Tagged), TaggedGen)  import qualified Numeric.LAPACK.Matrix.Triangular as Triangular import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape@@ -22,6 +23,7 @@  import qualified Numeric.Netlib.Class as Class +import qualified Data.Array.Comfort.Boxed as BoxedArray import qualified Data.Array.Comfort.Storable as Array import qualified Data.Array.Comfort.Shape as Shape import Data.Array.Comfort.Storable (Array)@@ -29,8 +31,9 @@  import Control.Applicative (liftA2, (<$>)) +import Data.Traversable (forM)+import Data.Maybe.HT (toMaybe) import Data.Tuple.HT (mapPair, swap)-import Data.Eq.HT (equating)  import qualified Test.QuickCheck as QC @@ -39,28 +42,7 @@    (Shape.C shape, Class.Floating a) => shape -> QC.Gen (Array shape a) genArray = Util.genArray 10 -equalArray ::-   (Shape.C shape, Eq shape, Class.Floating a) =>-   Array shape a -> Array shape a -> Bool-equalArray x y =-   if Array.shape x == Array.shape y-     then equalArrayBody x y-     else error "equalArray: shapes mismatch" -equalArrayBody ::-   (Shape.C shape, Class.Floating a) =>-   Array shape a -> Array shape a -> Bool-equalArrayBody =-   getEqualArray $-   Class.switchFloating-      (EqualArray $ equating Array.toList)-      (EqualArray $ equating Array.toList)-      (EqualArray $ equating Array.toList)-      (EqualArray $ equating Array.toList)--newtype EqualArray f a = EqualArray {getEqualArray :: f a -> f a -> Bool}-- dotProduct ::    (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>    (Vector ZeroInt a, Vector ZeroInt a) -> Bool@@ -221,35 +203,48 @@    equalArray x $ Matrix.reverseColumns (Matrix.reverseColumns x)  -mapHeight ::-   (heightA -> heightB) ->-   MatrixShape.General heightA width ->-   MatrixShape.General heightB width-mapHeight f shape =-   MatrixShape.general-      (MatrixShape.fullOrder shape)-      (f $ MatrixShape.fullHeight shape)-      (MatrixShape.fullWidth shape)+-- cf. Vector.genSwapVector+genSwapVector ::+   (Class.Floating a) =>+   Gen.Vector a Int+      (Maybe (Int,Int), (General ZeroInt ZeroInt a, Vector ZeroInt a))+genSwapVector =+   flip Gen.mapGen ((,) <$> Gen.matrix <|*.> Gen.vector) $ \_maxElem (m,x) -> do+      let set = Shape.indices $ Array.shape x+      ij <-+         forM (toMaybe (not $ null set) set) $ \s ->+            liftA2 (,) (QC.elements s) (QC.elements s)+      return (ij,(m,x)) -mapWidth ::-   (widthA -> widthB) ->-   MatrixShape.General height widthA ->-   MatrixShape.General height widthB-mapWidth f shape =-   MatrixShape.general-      (MatrixShape.fullOrder shape)-      (MatrixShape.fullHeight shape)-      (f $ MatrixShape.fullWidth shape)+swapColumns ::+   (Class.Floating a) =>+   (Maybe (Int,Int), (General ZeroInt ZeroInt a, Vector ZeroInt a)) ->+   QC.Property+swapColumns (mij,(m,x)) =+   let mx = m#>x+   in seq (Array.shape mx) $+      maybeProperty $+      flip fmap mij $ \(i,j) ->+         equalArray mx (Matrix.swapColumns i j m #> Vector.swap i j x) +_swapColumns ::+   (Class.Floating a) =>+   ((Int,Int), (General ZeroInt ZeroInt a, Vector ZeroInt a)) -> Bool+_swapColumns ((i,j),(m,x)) =+   equalArray+      (m#>x)+      (Matrix.swapColumns i j m #> Vector.swap i j x)++ zeroIntHeight ::    (Shape.C height, Shape.C width) =>    General height width a -> General ZeroInt width a-zeroIntHeight = Array.mapShape (mapHeight (zeroInt . Shape.size))+zeroIntHeight = Matrix.mapHeight (zeroInt . Shape.size)  zeroIntWidth ::    (Shape.C height, Shape.C width) =>    General height width a -> General height ZeroInt a-zeroIntWidth = Array.mapShape (mapWidth (zeroInt . Shape.size))+zeroIntWidth = Matrix.mapWidth (zeroInt . Shape.size)  reverseRowsStack ::    (Class.Floating a) =>@@ -312,6 +307,15 @@       (Matrix.dropRows k (Matrix.dropColumns k x))  +takeRowArray ::+   (Class.Floating a) => [Int] -> General ZeroInt ZeroInt a -> Bool+takeRowArray ixs x =+   equalArray+      (Matrix.adaptOrder x $+       Matrix.fromRows (Matrix.width x) $ map (Matrix.takeRow x) ixs)+      (Matrix.takeRowArray (BoxedArray.vectorFromList ixs) x)++ stackSplitRows ::    (Class.Floating a) => Int -> General ZeroInt ZeroInt a -> Bool stackSplitRows k x =@@ -433,6 +437,8 @@    (Show a, Class.Floating a, Eq a, RealOf a ~ ar, Class.Real ar) =>    [(String, Tagged a QC.Property)] testsVar =+   ("index",+      checkForAll (Indexed.genMatrixIndex Gen.matrix) Indexed.unitDot) :    ("dotProduct",       checkForAll ((,) <$> Gen.vector <.*.> Gen.vector) dotProduct) :    ("innerDot",@@ -497,6 +503,18 @@       checkForAllExtra (QC.choose (0,5)) Gen.matrix takeEqually) :    ("dropEqually",       checkForAllExtra (QC.choose (0,5)) Gen.matrix dropEqually) :+   ("takeRowArray",+      checkForAll+         (Gen.mapGen+            (\_maxElem x -> do+               let set = Shape.indices $ Matrix.height x+               ixs <-+                  if null set then return [] else QC.listOf $ QC.elements set+               return (ixs,x))+            Gen.matrix)+         (uncurry takeRowArray)) :+   ("swapColumns",+      checkForAll genSwapVector swapColumns) :    ("stackSplitRows",       checkForAllExtra (QC.choose (0,5)) Gen.matrix stackSplitRows) :    ("stackSplitColumns",

@@ -3,10 +3,12 @@ {-# LANGUAGE Rank2Types #-} module Test.Triangular (testsVar) where +import qualified Test.Indexed as Indexed import qualified Test.Generator as Gen import qualified Test.Utility as Util import Test.Generator ((<.*|>), (<|*.>), (<|*|>), (<|\|>), (<|=|>))-import Test.Utility (approx, approxArray, approxArrayTol, approxMatrix, Tagged)+import Test.Utility+         (approx, approxArray, approxArrayTol, approxMatrix, equalArray, Tagged)  import qualified Numeric.LAPACK.Matrix.Triangular as Triangular import qualified Numeric.LAPACK.Matrix.Square as Square@@ -14,7 +16,7 @@ import qualified Numeric.LAPACK.Matrix as Matrix import qualified Numeric.LAPACK.Vector as Vector import Numeric.LAPACK.Matrix.Triangular (Triangular)-import Numeric.LAPACK.Matrix (General, ZeroInt, (<#), (<#>), (#>))+import Numeric.LAPACK.Matrix (General, ZeroInt, (<#), (<#>), (#>), (|||), (===)) import Numeric.LAPACK.Vector (Vector) import Numeric.LAPACK.Scalar (RealOf, selectReal, absolute) @@ -26,7 +28,6 @@  import Control.Applicative ((<$>)) -import Data.Traversable (for) import Data.Tuple.HT (mapFst)  import qualified Test.QuickCheck as QC@@ -72,6 +73,58 @@       (Vector.sub (a#>x) (b#>x))  +expandTriangle ::+   (MatrixShape.Content lo, MatrixShape.TriDiag diag, MatrixShape.Content up,+    Shape.C sh, Class.Floating a) =>+   Triangular lo diag up sh a -> General sh sh a+expandTriangle = Matrix.fromFull . Triangular.toSquare++transposedZero ::+   (Class.Floating a, Shape.C width, Shape.C height) =>+   General height width a -> General width height a+transposedZero a = Vector.constant (Array.shape (Matrix.transpose a)) 0++stackLower ::+   (MatrixShape.TriDiag diag, Class.Floating a) =>+   (Triangular.FlexLower diag ZeroInt a,+    General ZeroInt ZeroInt a,+    Triangular.FlexLower diag ZeroInt a) ->+   Bool+stackLower (a,b,c) =+   let abc = expandTriangle $ Triangular.stackLower a b c+   in equalArray abc $ Matrix.adaptOrder abc $+         (expandTriangle a ||| transposedZero b+          ===+          b ||| expandTriangle c)++stackUpper ::+   (MatrixShape.TriDiag diag, Class.Floating a) =>+   (Triangular.FlexUpper diag ZeroInt a,+    General ZeroInt ZeroInt a,+    Triangular.FlexUpper diag ZeroInt a) ->+   Bool+stackUpper (a,b,c) =+   let abc = Matrix.fromFull $+         Triangular.toSquare $ Triangular.stackUpper a b c+   in equalArray abc $ Matrix.adaptOrder abc $+         (expandTriangle a ||| b+          ===+          transposedZero b ||| expandTriangle c)++stackSymmetric ::+   (MatrixShape.TriDiag diag, Class.Floating a) =>+   (Triangular.FlexSymmetric diag ZeroInt a,+    General ZeroInt ZeroInt a,+    Triangular.FlexSymmetric diag ZeroInt a) ->+   Bool+stackSymmetric (a,b,c) =+   let abc = expandTriangle $ Triangular.stackSymmetric a b c+   in equalArray abc $ Matrix.adaptOrder abc $+         (expandTriangle a ||| b+          ===+          Matrix.transpose b ||| expandTriangle c)++ multiplyIdentityVector ::    (MatrixShape.Content lo, MatrixShape.Content up, MatrixShape.TriDiag diag,     Class.Floating a, RealOf a ~ ar, Class.Real ar) =>@@ -260,11 +313,7 @@       let shape =             MatrixShape.Triangular                MatrixShape.NonUnit MatrixShape.autoUplo order size-      Array.fromList shape <$>-         (for (Shape.indices shape) $ \(r,c) ->-            if r==c-               then return (d!r)-               else Util.genElement maxElem)+      Util.genArrayExtraDiag maxElem shape (\r -> return (d!r))  eigensystem ::    (MatrixShape.DiagUpLo lo up, Eq lo, Eq up,@@ -377,10 +426,41 @@       (Triangular.forceNonUnitDiagonal <$> gen)  +checkFlexDiag ::+   String ->+   (forall diag.+    (MatrixShape.TriDiag diag, Eq diag, Show diag) =>+    GenTriangular lo diag up a -> Tagged a QC.Property) ->+   (forall diag.+    (MatrixShape.TriDiag diag, Eq diag, Show diag) =>+    GenTriangular lo diag up a) ->+   [(String, Tagged a QC.Property)]+checkFlexDiag name checker gen =+   (name++".Unit", checker (Triangular.forceUnitDiagonal <$> gen)) :+   (name++".NonUnit", checker (Triangular.forceNonUnitDiagonal <$> gen)) :+   []++ testsVar ::    (Show a, Show ar, Class.Floating a, Eq a, RealOf a ~ ar, Class.Real ar) =>    [(String, Tagged a QC.Property)] testsVar =+   checkAny "index"+      (\gen -> checkForAll (Indexed.genMatrixIndex gen) Indexed.unitDot)+      Gen.triangular ++++   checkFlexDiag "stackLower"+      (\gen ->+         checkForAll (Gen.stack3 gen (Matrix.transpose <$> Gen.matrix) gen)+         stackLower)+      Gen.triangular +++   checkFlexDiag "stackUpper"+      (\gen -> checkForAll (Gen.stack3 gen Gen.matrix gen) stackUpper)+      Gen.triangular +++   checkFlexDiag "stackSymmetric"+      (\gen -> checkForAll (Gen.stack3 gen Gen.matrix gen) stackSymmetric)+      Gen.triangular +++    checkAny "forceOrder"       (\gen ->          checkForAllExtra Util.genOrder

@@ -25,8 +25,10 @@ import qualified Data.List.HT as ListHT import qualified Data.Complex as Complex import Data.Complex (Complex((:+)))+import Data.Traversable (traverse) import Data.Monoid (Monoid(mempty,mappend)) import Data.Semigroup (Semigroup((<>)))+import Data.Eq.HT (equating)  import qualified Test.QuickCheck as QC import Test.ChasingBottoms.IsBottom (isBottom)@@ -45,6 +47,28 @@       (map Just ys ++ repeat Nothing)  +equalArray ::+   (Shape.C shape, Eq shape, Class.Floating a) =>+   Array shape a -> Array shape a -> Bool+equalArray x y =+   if Array.shape x == Array.shape y+     then equalArrayBody x y+     else error "equalArray: shapes mismatch"++equalArrayBody ::+   (Shape.C shape, Class.Floating a) =>+   Array shape a -> Array shape a -> Bool+equalArrayBody =+   getEqualArray $+   Class.switchFloating+      (EqualArray $ equating Array.toList)+      (EqualArray $ equating Array.toList)+      (EqualArray $ equating Array.toList)+      (EqualArray $ equating Array.toList)++newtype EqualArray f a = EqualArray {getEqualArray :: f a -> f a -> Bool}++ approx ::    (Class.Floating a, RealOf a ~ ar, Class.Real ar) => ar -> a -> a -> Bool approx tol x y = absolute (x-y) <= tol@@ -87,6 +111,10 @@       (Matrix.toRowMajor $ Matrix.fromFull y)  +maybeProperty :: (QC.Testable a) => Maybe a -> QC.Property+maybeProperty = maybe (QC.property QC.Discard) QC.property++ genReal :: (Class.Real a) => Integer -> QC.Gen a genReal n = fromInteger <$> QC.choose (-n,n) @@ -100,8 +128,21 @@ genArray ::    (Shape.C shape, Class.Floating a) =>    Integer -> shape -> QC.Gen (Array shape a)-genArray n shape =-   Array.fromList shape <$> replicateM (Shape.size shape) (genElement n)+genArray maxElem shape =+   Array.fromList shape <$> replicateM (Shape.size shape) (genElement maxElem)++genArrayIndexed ::+   (Shape.Indexed shape, Class.Floating a) =>+   shape -> (Shape.Index shape -> QC.Gen a) -> QC.Gen (Array shape a)+genArrayIndexed shape f =+   Array.fromList shape <$> traverse f (Shape.indices shape)++genArrayExtraDiag ::+   (Shape.Indexed shape, Shape.Index shape ~ (i,i), Eq i, Class.Floating a) =>+   Integer -> shape -> (i -> QC.Gen a) -> QC.Gen (Array shape a)+genArrayExtraDiag maxElem shape diag =+   genArrayIndexed shape $+      \(r,c) -> if r==c then diag r else genElement maxElem   select :: [a] -> QC.Gen (a, [a])

@@ -26,11 +26,14 @@ import Test.ChasingBottoms.IsBottom (isBottom)  +singleton :: (Class.Floating a) => a -> Bool+singleton x = Util.equalArray (Vector.singleton x) (Vector.constant () x)+ appendTakeDrop ::    (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>    Int -> Vector ZeroInt a -> Bool appendTakeDrop n x =-   Util.approxArray x $+   Util.equalArray x $    Array.mapShape (zeroInt . Shape.size)       (Vector.append (Vector.take n x) (Vector.drop n x)) @@ -39,11 +42,33 @@    (Vector ZeroInt a, Vector ZeroInt a) -> Bool takeLeftRightAppend (x,y) =    let xy = Vector.append x y-   in Util.approxArray x (Vector.takeLeft xy)+   in Util.equalArray x (Vector.takeLeft xy)       &&-      Util.approxArray y (Vector.takeRight xy)+      Util.equalArray y (Vector.takeRight xy)  +genSwapVector :: (Class.Floating a) => TaggedGen a ((Int,Int), Vector ZeroInt a)+genSwapVector =+   flip fmap genVector $ \gen -> do+      x <- gen `QC.suchThat` (not . null . Vector.toList)+      let set = Shape.indices $ Array.shape x+      i <- QC.elements set+      j <- QC.elements set+      return ((i,j),x)++swapInverse ::+   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>+   ((Int,Int), Vector ZeroInt a) -> Bool+swapInverse ((i,j),x) =+   Util.equalArray x $ Vector.swap i j $ Vector.swap i j x++swapCommutative ::+   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>+   ((Int,Int), Vector ZeroInt a) -> Bool+swapCommutative ((i,j),x) =+   Util.equalArray (Vector.swap i j x) (Vector.swap j i x)++ normInf ::    (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>    Vector ZeroInt a -> Bool@@ -111,12 +136,18 @@     Class.Floating a, Eq a, RealOf a ~ ar, Class.Real ar, RealOf ar ~ ar) =>    [(String, Tagged a QC.Property)] testsVar =+   ("singleton",+      checkForAll Gen.scalar singleton) :    ("appendTakeDrop",       Gen.withExtra checkForAll          (QC.getNonNegative <$> QC.arbitrary) Gen.vector appendTakeDrop) :    ("takeLeftRightAppend",       Util.checkForAllPlain          (liftA2 (liftA2 (,)) genVector genVector) takeLeftRightAppend) :+   ("swapInverse",+      Util.checkForAllPlain genSwapVector swapInverse) :+   ("swapCommutative",+      Util.checkForAllPlain genSwapVector swapCommutative) :    ("normInf",       checkForAll Gen.vector normInf) :    ("normInf1",