packages feed

fpnla (empty) → 0.1

raw patch · 9 files changed

+659/−0 lines, 9 filesdep +basesetup-changed

Dependencies added: base

Files

+ LICENSE view
@@ -0,0 +1,12 @@+Copyright (c) 2014, Mauro Blanco, Pablo Perdomo+All rights reserved.++Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:++1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.++2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.++3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ fpnla.cabal view
@@ -0,0 +1,41 @@+name:           fpnla+version:        0.1+cabal-version:  >=1.2+build-type:     Simple+stability:      Experimental+author:         Pablo Perdomo, Mauro Blanco+synopsis:       A library for NLA operations+description:+    .+    This package defines a framework for linear algebra operations, allowing:+    .+    - Definition of multiple instances of BLAS and LAPACK operations.+    .+    - Definition of multiple representations of vectors and matrices.+    .+    - Arbitrary combination of strategies and structure representations.+    .+    - Type-safe manipulation of context information associated to each strategy.+    .+    - Definition of specialized strategies for a given representation.+    .+    The framework is based on BLAS (<http://www.netlib.org/blas/>) and LAPACK (<http://www.netlib.org/lapack/>) linear algebra libraries, as these are well known libraries in the area. For this reason the design of the framework is strongly oriented to these libraries. Anyway it is possible to easily define operations not considered in these libraries and still maintaining the above properties.+    .+tested-with:    GHC==7.6.3+maintainer:     Pablo Perdomo <pperdomo@fing.edu.uy>, Mauro Blanco <mblanco@fing.edu.uy>+category:       Math+license:        BSD3+license-file:   LICENSE++library+  hs-source-dirs:  src+  build-depends:   base >= 4 && < 5+  ghc-options:     -Wall -fno-warn-orphans+  exposed-modules: +                   FPNLA.Matrix,+                   FPNLA.Matrix.Utils,+                   FPNLA.Operations.BLAS,+                   FPNLA.Operations.LAPACK,+                   FPNLA.Operations.Parameters+  other-modules:   FPNLA.Utils+
+ src/FPNLA/Matrix.hs view
@@ -0,0 +1,227 @@+{-# LANGUAGE FlexibleInstances     #-}+{-# LANGUAGE MultiParamTypeClasses #-}++-- | This module defines classes used to handle Matrices and Vectors in a generic way.+module FPNLA.Matrix (+    -- * Vector+    Vector(..),+    -- * Matrix+    Matrix(..),+    -- * MatrixVectror+    MatrixVector(..),+    -- * Miscellaneous+    cantRows_m,+    cantCols_m,+    size_m,+    asColumn_vm,+    diagonalBlock,+    verticalBlock+) where++import FPNLA.Utils (iif, splitSlice)++-- | This class represents a vector structure+-- For any vector data structure used in this framework it must have an instance of this class.+-- A minimal instance of this class must provide an implementation for 'generate_v', 'elem_v' and 'length_v'.+-- The index of the first element stored in the vector must be 0.+-- This is a multi-param type class because usually libraries of Matrix and Vector imposes restrictions over the type of the elements (for example, to have a "Storable" instance).+class Vector v e where+    -- *Constructors. This functions allows to create a new vector.+    +    -- | Creates a new vector from the length of the vector and a function that returns the value for each index.+    generate_v :: Int -- ^ The length of the vector+        -> (Int -> e) -- ^ A function that returns the value for each index+        -> v e+    -- | Creates a new vector from a list, th vector will contain all (and only) the elements of the list, and in the same order.+    fromList_v :: [e] -- ^ A list that contains the elements of the vector. The first element of the list will be the first element of the vector and so on.+        -> v e+    -- | Join multiple vectors in one single vector.+    concat_v :: [v e] -- ^ A list containing the vectors to be joined. If the list is @[<1>, <2,3>, <4,5,6>]@ then the vector will be @<1,2,3,4,5,6>@.+        -> v e+    +    -- * Selectors. This functions returns information from a vector.+    -- | Returns the element from a given index. The index must be between 0 and the length of the vector.+    elem_v :: Int -> v e -> e+    -- | Returns the length of the vector.+    length_v :: v e -> Int+    +    -- * High-order. This functions are the classics /fold/, /map/ and /zipWith/ but are in the class to allow to get an specific implementation from a specific data structure in order to achieve better performance.+    -- | The common /foldr/ function but specific for this instance.+    foldr_v :: (e -> b -> b) -> b -> v e -> b+    -- | The common /map/ function but specific for this instance.+    map_v :: (e -> e) -> v e -> v e+    -- | The common /zipWith/ function but specific for this instance.+    zipWith_v :: (e -> e -> e) -> v e -> v e -> v e++    -- Default implementations:+    fromList_v l = generate_v (length l) (\i ->  l !! i)+    concat_v [] = generate_v 0 undefined+    concat_v [v] = v+    concat_v (v1:v2:vs) = concat_v (join:vs)+        where l1 = length_v v1+              l2 = length_v v2+              gen i | i < l1 = elem_v i v1+                    | otherwise = elem_v (i-l1) v2+              join = generate_v (l1 + l2) gen+    foldr_v f z v = go 0+        where len = length_v v+              go i | i == len   = z+                   | otherwise  = f (elem_v i v) (go (i+1))+    map_v f v = generate_v (length_v v) (\i -> f (elem_v i v))+    zipWith_v f v1 v2 = generate_v (length_v v1) (\i -> f (elem_v i v1) (elem_v i v2))+++-- | This class represents any Matrix structure+-- For any matrix data structure used in this framework it must have an instance of this class.+-- A minimal instance of this class must provide an implementation for 'generate_m', 'dim_m' and 'elem_m'.+-- The index of the first element stored in the vector must be (0,0).+-- This is a multi-param type class because usually libraries of Matrix and Vector imposes restrictions over the type of the elements (for example, to have an "Storable" instance).+class Matrix m e where+    -- * Constructors. This functions allows to create a new matrix.+    -- | Creates a new matrix from the dimension (length of rows and columns) of the matrix and a function that returns the value for each index.+    generate_m :: Int -- ^ The length of the rows+        -> Int -- ^ A function that returns the value for each index+        -> (Int -> Int -> e) -- ^ +        -> m e+    -- | Creates a new matrix from the dimension (length of rows and columns) of and a list.+    fromList_m :: Int -> Int -> [e] -> m e+    -- | Given a matrix, creates a new matrix that will be the transpose of the original.+    transpose_m :: m e -> m e++    -- * Selectors. This functions returns information from a vector.+    -- | Returns the dimension of a given matrix. The first element in the pair is the number of rows and the second the number of columns.+    dim_m :: m e -> (Int, Int)+    -- | Returns the element that is in the position @(i,j)@.+    elem_m :: Int -- ^ The row component of the position (@i@).+        -> Int -- ^ The column component of the position (@j@).+        -> m e -- ^ The matrix.+        -> e++    -- * High-order. This functions are the classics /fold/, /map/ and /zipWith/ but are in the class to allow to get an specific implementation from a specific data structure in order to achieve better performance.+    -- | The common /map/ function but specific for this instance.+    map_m :: (e -> e) -> m e -> m e+    -- | The common /zipWith/ function but specific for this instance.+    zipWith_m :: (e -> e -> e) -> m e -> m e -> m e++    -- * Blocks. This functions are constructors too, but this splits or allow to get a sub-matrix from a given matrix.+    -- | Given an initial position @(i,j)@, a dimension @(r,c)@ and a matrix @m@, creates a new matrix that is a sub matrix of @m@ containing the elements in positions @[(p1, p2) | p1 <- [i..i+r], p2 <- [j..j+c]]@+    -- For example, if @m@ is:+    -- @1 2 3+    -- 4 5 6+    -- 7 8 9@+    -- and @(i,j)@ and @(r,c)@ are @(0,1)@ and @(2,3)@ respectively, then the result will be a new matrix containing:+    -- @2 3+    -- 5 6+    -- 8 9@+    subMatrix_m :: Int -- ^ @i@+        -> Int -- ^ @j@+        -> Int -- ^ @r@+        -> Int -- ^ @c@+        -> m e -- ^ @m@+        -> m e+    -- | Constructs a new matrix from a list of list where the elements are sub-matrices.+    -- This operations is the analogue of 'concat_v' operation from the 'Vector' class.+    -- The \'borders\' of the sub-matrices will be joined.+    -- The dimensions of the sum-matrices must be compatibles.+    fromBlocks_m :: [[m e]] -> m e+    -- | Split a matrix @m@ in blocks of size @(r,c)@ (the dimension of the blocks in the right and the bottom may be smaller).+    toBlocks_m :: Int -- ^ @r@+        -> Int -- ^ @c@+        -> m e -- ^ @m@+        -> [[m e]]++    -- Default implementations:+    fromList_m m n l = generate_m m n (\i j ->  ls !! i !! j)+        where ls = splitSlice n l+    transpose_m m = generate_m (cantCols_m m) (cantRows_m m) (\i j -> elem_m j i m)+    map_m f m = generate_m (cantRows_m m) (cantCols_m m) (\i j -> f (elem_m i j m))+    zipWith_m f m1 m2 = generate_m (cantRows_m m1) (cantCols_m m1) (\i j -> f (elem_m i j m1) (elem_m i j m2))+    subMatrix_m posI posJ cantFilas cantCols m = generate_m cantFilas cantCols (\i j -> elem_m (i + posI) (j + posJ) m)+    fromBlocks_m [] = fromList_m 0 0 []+    fromBlocks_m [[]] = fromList_m 0 0 []+    fromBlocks_m bss = +        generate_m cr cc $ \i j -> elem_m (iToB i) (jToB j) $ bss !! iToBSS i !! jToBSS j+        where bssm = length bss+              bssn = length $ head bss+              (br, bc) = dim_m . head . head $ bss+              lbr = cantRows_m . head . last $ bss+              lbc = cantCols_m . last . head $ bss+              ilimit = (bssm-1) * br+              jlimit = (bssn-1) * bc+              cr = ilimit + lbr+              cc = jlimit + lbc+              iToBSS i = iif (i < ilimit) (i `div` br) (bssm-1)+              jToBSS j = iif (j < jlimit) (j `div` bc) (bssn-1)+              iToB i = iif (i < ilimit) (i `mod` br) ((i-ilimit) `mod` lbr)+              jToB j = iif (j < jlimit) (j `mod` bc) ((j-jlimit) `mod` lbc)+    toBlocks_m brs bcs m+        | brs == 0 || bcs == 0 = []+        | otherwise = map (rowBlocks 0) [i * brs | i <- [0 .. (rs `div` brs + iif (rs `mod` brs > 0) 1 0) - 1]]+        where+            (rs, cs) = dim_m m+            getBound s b c = iif (s + b > c) (c - s) b+            rowBlocks sJ sI+                    | sJ >= cs = []+                    | otherwise = subMatrix_m sI sJ (getBound sI brs rs) (getBound sJ bcs cs) m : rowBlocks (sJ + bcs) sI+++-- | This class allows us to leave separated the two concepts of 'Vector' and 'Matrix'.+-- This class establishes a link between a vector and a matrix structure.+-- We provide a default instance for any 'Vector' and 'Matrix' instances so any vector could be used with any matrix whenever the type of the elements are the same.+class (Vector v e, Matrix m e) => MatrixVector m v e where+    -- | Extracts the @i@-th row of the matrix @m@ and converts it in the vector structure.+    row_vm :: Int -> m e -> v e+    -- | Extracts the @j@-th column of the matrix @m@ and converts it in the vector structure.+    col_vm :: Int -> m e -> v e+    -- | Joins a list of column vectors in a matrix structure.+    -- The length of every vector in the list must be the same.+    fromCols_vm :: [v e] -> m e+    -- | Returns a list containing every column of the matrix converted in the vector structure.+    toCols_vm :: m e -> [v e]++    -- Default implementations:+    row_vm i m = generate_v (cantCols_m m) (\j -> elem_m i j m)+    col_vm j m = generate_v (cantRows_m m) (\i -> elem_m i j m)+    fromCols_vm css | cc == 0 = fromList_m 0 0 []+                  | otherwise = generate_m cr cc (\i j -> elem_v i $ css !! j)+        where cc = length css+              cr = length_v $ head css+    toCols_vm m = map (`col_vm` m) [0 .. cantCols_m m - 1]+++instance (Vector v e, Matrix m e) => MatrixVector m v e++-- | Returns the number of rows of the matrix.+cantRows_m :: (Matrix m e) => m e -> Int+cantRows_m = fst . dim_m++-- | Returns the number of columns of the matrix.+cantCols_m :: (Matrix m e) => m e -> Int+cantCols_m = snd . dim_m++-- | Returns the size of the matrix (number of rows x number of columns).+size_m :: (Matrix m e) => m e -> Int+size_m = uncurry (*) . dim_m++-- | Converts a vector in a matrix of one column.+asColumn_vm :: MatrixVector m v e => v e -> m e+asColumn_vm = fromCols_vm . (:[])++-- | Splits a matrix @m@ into blocks of size @(r,c)@ and returns the block on the position @(i,i)@.+diagonalBlock :: (Matrix m e) =>  (Int, Int) -- ^ @(r,c)@+    -> Int -- ^ @i@+    -> m e -- ^ @m@+    -> m e+diagonalBlock (blockDimM, blockDimN) i m =+    subMatrix_m (blockDimM*i) (blockDimN*i) (min blockDimM (cantRows_m m - (blockDimM * i))) (min blockDimN (cantCols_m m - (blockDimN * i))) m++++-- | Splits a matrix @m@ into blocks of size @(s, cantCols_m m)@ and returns the block on the position @i@.+verticalBlock :: (Matrix m e) =>  Int -- ^ @s@+    -> Int -- ^ @i@+    -> m e -- ^ @m@+    -> m e+verticalBlock blockDim i m =+    subMatrix_m (blockDim*i) 0 (min blockDim (cantRows_m m - (blockDim * i))) (cantCols_m m) m+
+ src/FPNLA/Matrix/Utils.hs view
@@ -0,0 +1,15 @@+-- | This module defines commonly useful functions that are related specifically with vectors and matrices.+module FPNLA.Matrix.Utils where+++import FPNLA.Matrix (Matrix(..))+++-- | Prints a matrix to the standard output.+-- This operation requires the elements of the matrix to have an instance of 'Show' but does not requires a 'Show' instance for the matrix data type.+print_m :: (Show e, Matrix m e) => m e -> IO ()+print_m mi = for 0 0+    where (m,n) = dim_m mi+          for i j | i >= m = return ()+                  | j < n = (putStr . show $ elem_m i j mi) >> putStr " " >> for i (j+1)+                  | j >= n = putStrLn "" >> for (i+1) 0
+ src/FPNLA/Operations/BLAS.hs view
@@ -0,0 +1,69 @@+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE ScopedTypeVariables   #-}+{-# LANGUAGE TypeFamilies          #-}++-- | This module defines all the BLAS (Basic Linear Algebra Subprograms) operations supported by the framework.+-- See <http://www.netlib.org/blas/> for more information about BLAS and <http://www.ugcs.caltech.edu/~srbecker/blasqr_betterFonts.pdf> for a quick description of all BLAS operation signatures and behaviour.+module FPNLA.Operations.BLAS (+    -- *  Level One+    -- | Vector-Vector operations+    DOT(..),+    -- *  Level Two+    -- | Matrix-Vector operations+    GEMV(..),+    -- *  Level Three+    -- | Matrix-Matrix operations+    SYRK(..),+    GEMM(..),+    TRSM(..)++) where++import FPNLA.Matrix                (MatrixVector, Vector)+import FPNLA.Operations.Parameters (Elt (), ResM (), ResV (), ResS(),+                                    StratCtx (), TransType (),+                                    TriangType (), UnitType ())++-- | Defines the signature of the level-1 BLAS /dot/ operation in the framework.+class (Elt e, Vector v e) => DOT s v e where+    dot :: StratCtx s -- ^ The context of the operation+        -> v e -- ^ A vector /x/+        -> v e -- ^ A vector /y/+        -> ResS s e -- ^ The scalar product between /x/ and /y/++-- | Defines the signature of the level-2 BLAS /gemv/ operation in the framework.+class (Elt e, MatrixVector m v e) => GEMV s m v e where+    gemv :: StratCtx s -- ^ The context of the operation+        -> TransType (m e) -- ^ A matrix /A/+        -> v e -- ^ A vector /x/+        -> e -- ^ A scalar /alpha/+        -> e -- ^ A scalar /beta/ +        -> v e  -- ^ A vecor /y/+        -> ResV s v e -- ^ @alpha * A * x + beta * y@++-- | Defines the signature of the level-3 BLAS /syrk/ operation in the framework.+class (Elt e, MatrixVector m v e) => SYRK s m v e where+    syrk :: StratCtx s -- ^ The context of the operation+        -> e -- ^ A scalar /alpha/+        -> TransType (m e) -- ^ A matrix /A/+        -> e -- ^ A scalar /beta/ +        -> TriangType (m e) -- ^ A triangular matrix /C/+        -> ResM s v m e -- ^ @alpha * A * A' + beta * C@ where @A'@ is the conjugate transposed of /A/.++-- | Defines the signature of the level-3 BLAS /gemm/ operation in the framework.+class (Elt e, MatrixVector m v e) => GEMM s m v e where+    gemm :: StratCtx s -- ^ The context of the operation+        -> TransType (m e) -- ^ A matrix /A/+        -> TransType (m e) -- ^ A matrix /B/+        -> e -- ^ A scalar /alpha/+        -> e -- ^ A scalar /beta/ +        -> m e -- ^ A matrix /C/+        -> ResM s v m e -- ^ @alpha * A * B + beta * C@++-- | Defines the signature of the level-3 BLAS /trsm/ operation in the framework.+class (Elt e, MatrixVector m v e) => TRSM s m v e where+    trsm :: StratCtx s -- ^ The context of the operation+        -> e -- ^ A scalar /alpha/+        -> TransType (TriangType (UnitType (m e))) -- ^ A triangular matrix /A/+        -> m e -- ^ A matrix /B/+        -> ResM s v m e -- ^ @alpha * A_inv * B@ where @A_inv@ is the inverse of /A/.
+ src/FPNLA/Operations/LAPACK.hs view
@@ -0,0 +1,25 @@+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE ScopedTypeVariables   #-}+{-# LANGUAGE TypeFamilies          #-}++-- | This module defines all the LAPACK (Linear Algebra PACKage) operations supported by the framework.+-- See <http://www.netlib.org/lapack/> for more information about LAPACK.+module FPNLA.Operations.LAPACK (++    POTRF(..),++) where++import FPNLA.Matrix                (MatrixVector)+import FPNLA.Operations.Parameters (Elt (), ResM (), StratCtx (), TriangType ())++++-- Po: Symmetric matrix or Hermitian matrix positive definite+-- Trs: Triangular Factorization++-- | Defines the signature of the LAPACK /potrf/ operation in the framework.+-- This operation takes a symmetric (or hermitian) positive definite (SPD) matrix (flagged with TriangType) and computes the Cholesky factorization of the matrix.+-- The Cholesky decomposition of an SPD matrix /M/ is a lower triangular matrix /L/ where /M = L L*/ being /L*/ the conjugate transpose of /L/.+class (Elt e, MatrixVector m v e) => POTRF s m v e where+    potrf :: StratCtx s -> TriangType (m e) -> ResM s v m e
+ src/FPNLA/Operations/Parameters.hs view
@@ -0,0 +1,176 @@+{-# LANGUAGE TypeFamilies #-}++module FPNLA.Operations.Parameters(+    -- * Elements+    Elt(..),+    -- * Strategies and contexts+    StratCtx(), +    +    -- * Result type+    -- | In BLAS it's common that operations in higher levels use operations in the lower levels, so, an operation in level three that by its signature manipulates matrices only, internally uses level two operations that manipulates vectors. In order to avoid the /show . read/ problem, the type of the vector (or any other internal data type) must appear in the signature of an operation.+    -- To solve the problem we use phantom types to pass the internally used types to the Haskell type system.+    ResM(),+    ResV(),+    ResS(),+    blasResultM,+    blasResultV,+    blasResultS,+    getResultDataM,+    getResultDataV,+    getResultDataS,+    +    -- * Miscellaneous+    TransType(..),+    UnitType(..),+    TriangType(..),+    unTransT, +    unUnitT, +    unTriangT, +    elemTrans_m,+    dimTrans_m,+    elemSymm,+    dimTriang,+    elemUnit_m,+    dimUnit_m,+    elemTransUnit_m,+    dimTransUnit_m,+    transTrans_m+    +) where++import FPNLA.Matrix(MatrixVector(..), Vector(..), Matrix(..))+import Data.Complex (Complex, conjugate)+import Data.Tuple (swap)+++-- | This class represents the elements that can be used in the BLAS operations.+-- The elements in BLAS are real or complex numbers, so we provide default instances for the Haskell 'Double', 'Float' and 'Complex' types.+class (Eq e, Floating e) => Elt e where+    -- | Returns the conjugate of a number. For real numbers it's the identity function and for complex numbers it's the common 'Complex.conjugate' function.+    getConjugate :: e -> e+    getConjugate = id++instance Elt Double+instance Elt Float+instance (RealFloat e) => Elt (Complex e) where+    getConjugate = conjugate++-- | This type family is used to represent the /context/ of an operation.+-- A particular implementation is a combination of an algorithm and a parallelism technique, and we call it a /strategy/. A particular strategy may need particular information to execute. For example, an operation that computes the matrix-matrix multiplication by splitting the matrices in blocks must require the size of the blocks.+-- With this context we allows to pass any additional information that the operation needs to execute as parameters, but maintaining a common signature.+-- The /s/ type parameter is the strategy so, there must exist a Haskell data type to represent a particular strategy.+type family StratCtx s :: *+++-- | The 'ResM' data type is used as result of level three BLAS operations and returns a matrix /m/ of elements /e/ and contains the strategy /s/ and vector /v/ as phantom types.+data ResM s (v :: * -> *) m e = ResM { unResM :: m e } deriving (Show)+-- | The 'ResV' data type is used as result of level two BLAS operations and returns a vector /v/ of elements /e/ and contains the strategy /s/ as phantom types.+data ResV s v e = ResV { unResV :: v e } deriving (Show)+-- | The 'ResS' data type is used as result of level one BLAS operations and returns an scalar /e/ and contains the strategy /s/ as phantom types.+data ResS s e = ResS { unResS :: e } deriving (Show)+++-- | Wrap a matrix into a 'ResM'.+blasResultM :: (MatrixVector m v e) => m e -> ResM s v m e+blasResultM = ResM+-- | Unwrap a matrix from a 'ResM'.+getResultDataM :: (MatrixVector m v e) => ResM s v m e -> m e+getResultDataM = unResM++-- | Wrap a vector into a 'ResV'.+blasResultV :: (Vector v e) => v e -> ResV s v e+blasResultV = ResV+-- | Unwrap a vector from a 'ResV'.+getResultDataV :: (Vector v e) => ResV s v e -> v e+getResultDataV = unResV++-- | Wrap a scalar into a 'ResS'.+blasResultS :: e -> ResS s e+blasResultS = ResS+-- | Unwrap a scalar from a 'ResS'.+getResultDataS :: ResS s e -> e+getResultDataS = unResS+++-- | Indicates if a matrix must be considered as normal, transposed or transposed conjugated.+-- This is part of the common flags in the BLAS operation signatures and it's useful to work with a transposed matrix without really computing the transposed matrix.+data TransType m = Trans m | NoTrans m | ConjTrans m deriving (Eq, Show)+-- | Indicates if a matrix must be considered as unitary or not. An unitary matrix is a matrix that contains ones in the diagonal.+-- This is part of the common flags in the BLAS operation signatures.+data UnitType m = Unit m | NoUnit m deriving (Eq, Show)+-- | Indicates that a matrix is symmetric and with which triangular part of the matrix the operation is going to work ('Upper' or 'Lower').+-- The operation only will see the indicated part of the matrix and should not try to access the other part.+-- This is part of the common flags in the BLAS operation signatures.+data TriangType m = Lower m | Upper m deriving (Eq, Show)++-- | Given a data type flagged by a TransType, returns a pair containing the TransType constructor and the data type.+unTransT :: TransType a -> (b -> TransType b, a)+unTransT (Trans a) = (Trans, a)+unTransT (NoTrans a) = (NoTrans, a)+unTransT (ConjTrans a) = (ConjTrans, a)++-- | Given a data type flagged by a UnitType, returns a pair containing the UnitType constructor and the data type.+unUnitT :: UnitType a -> (b -> UnitType b, a)+unUnitT (Unit a) = (Unit, a)+unUnitT (NoUnit a) = (NoUnit, a)++-- | Given a data type flagged by a TriangType, returns a pair containing the TriangType constructor and the data type.+unTriangT :: TriangType a -> (b -> TriangType b, a)+unTriangT (Lower a) = (Lower, a)+unTriangT (Upper a) = (Upper, a)+++-- | Given an /i,j/ position and a TransType flagged matrix, returns the element in that position without computing the transpose.+elemTrans_m :: (Elt e, Matrix m e) => Int -> Int -> TransType (m e) -> e+elemTrans_m i j (NoTrans m) = elem_m i j m+elemTrans_m i j (Trans m) = elem_m j i m+elemTrans_m i j (ConjTrans m) = getConjugate $ elem_m j i m+++-- | Given a TransType flagged matrix, returns the dimension of the matrix without computing the transpose.+dimTrans_m :: (Matrix m e) => TransType (m e) -> (Int, Int)+dimTrans_m (NoTrans m) = dim_m m+dimTrans_m (ConjTrans m) = swap $ dim_m m+dimTrans_m (Trans m) = swap $ dim_m m+++-- | Given an /i,j/ position and a TransType flagged matrix, returns the element in that position only accessing the part indicated by the TransType.+elemSymm :: (Elt e, Matrix m e) => Int -> Int -> TriangType (m e) -> e+elemSymm i j (Upper m)+    | i > j = elem_m j i m+    | otherwise = elem_m i j m+elemSymm i j (Lower m)+    | i > j = elem_m i j m+    | otherwise = elem_m j i m++-- | Given a TransType flagged matrix, returns the dimension of the matrix.+dimTriang :: (Matrix m e) => TriangType (m e) -> (Int, Int)+dimTriang = dim_m . snd . unTriangT++-- | Given an /i,j/ position and a UnitType flagged matrix, returns the element in that position. If the matrix is flagged as Unit and /i == j/ (the element is in the diagonal) returns one.+elemUnit_m :: (Elt e, Matrix m e) => Int -> Int -> UnitType (m e) -> e+elemUnit_m i j (Unit m)+    | i == j = 1+    | otherwise = elem_m i j m+elemUnit_m i j (NoUnit m) = elem_m i j m++-- | Given a UnitType flagged matrix, returns the dimension of the matrix.+dimUnit_m :: (Matrix m e) => UnitType (m e) -> (Int, Int)+dimUnit_m (Unit m) = dim_m m+dimUnit_m (NoUnit m) = dim_m m++-- | Given an /i,j/ position and a TransType-UnitType flagged matrix, returns the element in that position without computing the transpose.+elemTransUnit_m :: (Elt e, Matrix m e) => Int -> Int -> TransType (UnitType (m e)) -> e+elemTransUnit_m i j (NoTrans pmA) = elemUnit_m i j pmA+elemTransUnit_m i j (Trans pmA) = elemUnit_m j i pmA+elemTransUnit_m i j (ConjTrans pmA) = getConjugate $ elemUnit_m j i pmA++-- | Given a TransType-UnitType flagged matrix, returns the dimension of the matrix.+dimTransUnit_m :: Matrix m e => TransType (UnitType (m e)) -> (Int, Int)+dimTransUnit_m = dimUnit_m . snd . unTransT++-- | Given a TransType flagged matrix, computes and returns its transpose.+transTrans_m :: (Elt e, Matrix m e) => TransType (m e) -> m e+transTrans_m (NoTrans m) = m+transTrans_m (ConjTrans m) = map_m getConjugate $ transpose_m m+transTrans_m (Trans m) = transpose_m m
+ src/FPNLA/Utils.hs view
@@ -0,0 +1,92 @@+{-# LANGUAGE FlexibleInstances    #-}+{-# LANGUAGE ScopedTypeVariables  #-}+{-# LANGUAGE UndecidableInstances #-}++-- | This module defines commonly useful functions that are not specifically related with vectors and matrices.+module FPNLA.Utils (+    iif,+    mapPair,+    splitSlice,++    Truncable(..),++    prevEnum,+    nextEnum,+    enumerate++) where++import Data.Complex (Complex ((:+)))++-- | Takes two functions @f@ and @g@, and a pair, and applies @f@ to the first element and @g@ to the second element of the pair.+mapPair :: (a -> b) -- ^ @f@+    -> (c -> d) -- ^ @g@+    -> (a, c) -- ^ Input pair @(a,b)@+    -> (b, d) -- ^ Output pair @(f a, g b)@+mapPair f g (a, c) = (f a,  g c)++-- | It's equivalent to @if cond then e1 else e2@, but defined as a function.+iif :: Bool -- ^ @cond@+    -> a -- ^ @e1@+    -> a -- ^ @e2@+    -> a -- ^ @if cond then e1 else e2@+iif cond expTrue expFalse = if cond then expTrue else expFalse+++-- | Splits a list into pieces of size @n@.+-- The last piece could be smaller if the length of the list is not multiple of @n@.+splitSlice :: Int -- ^ @n@+    -> [a] -- ^ The input list+    -> [[a]] -- ^ The input list divided in pieces of size @n@ +splitSlice 0 _ = error "n must be greater than 0"+splitSlice _ [] = []+splitSlice n ls = sl : splitSlice n sls+    where (sl, sls) = splitAt n ls++-- | This class is used to compare two float point numbers.+-- With float point numbers, two different computations that theoretically returns the same value in fact returns different values. Even if the difference is so small that is ignorable, if we test using equality the test fail.+-- For this, instead to check for equality we define that two float point numbers are the same if the distance between them is less than some epsilon value.+class Truncable a where+    -- | Check if the distance between two float point values is less than some epsilon value.+    inEpsilonRange :: Real e => e -- ^ Epsilon+        -> a -- ^ A float point value+        -> a -- ^ Another float point value+        -> Bool++-- | We provide a default instance of 'Truncable' for any data type that has 'Floating' and 'RealFrac' instances.+instance (Floating n, RealFrac n) => Truncable n where+    inEpsilonRange epsilon n1 n2 = expN1 == expN2 && abs(mantN1 - mantN2) <= eps+        where+            eps = fromRational $ toRational epsilon+            shift n d1 = n * (10 ** fromIntegral d1)+            getExponent r+                    | r == 0 = 0+                    | r < 0 = getExponent (abs r)+                    | r > 10 = 1 + getExponent (r / 10)+                    | r < 1 = (-1) + getExponent (r * 10)+                    | otherwise = 0+            expN1 :: Integer = getExponent n1+            expN2 :: Integer = getExponent n2+            mantN1 = shift n1 (-expN1)+            mantN2 = shift n2 (-expN2)++-- | We define a specific instance of 'Truncable' for the 'Complex' data type.+instance (RealFrac n, Truncable n) => Truncable (Complex n) where+    inEpsilonRange epsilon (r1:+i1) (r2:+i2) = inEpsilonRange epsilon r1 r2 && inEpsilonRange epsilon i1 i2+++-- | Returns the predecessor of a enumerated value.+-- The operation is cyclic so, the predecessor of the first value is the last value.+prevEnum :: (Eq a, Bounded a, Enum a) => a -> a+prevEnum a | a == minBound = maxBound+           | otherwise = pred a++-- | Returns the successor of a enumerated value.+-- The operation is cyclic so, the successor of the last value is the first value.+nextEnum :: (Eq a, Bounded a, Enum a) => a -> a+nextEnum a | a == maxBound = minBound+           | otherwise = succ a++-- | Returls a list containing all the values of an enumerated.+enumerate :: forall a. (Enum a, Bounded a) => [a]+enumerate = map toEnum [fromEnum (minBound :: a) .. fromEnum (maxBound :: a)]