packages feed

fpnla-examples (empty) → 0.1

raw patch · 46 files changed

+2903/−0 lines, 46 filesdep +HUnitdep +QuickCheckdep +acceleratesetup-changed

Dependencies added: HUnit, QuickCheck, accelerate, array, base, data-default, deepseq, fpnla, fpnla-examples, hmatrix, linear-algebra-cblas, monad-par, parallel, random, repa, tagged, test-framework, test-framework-hunit, test-framework-quickcheck2, time, vector

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-examples.cabal view
@@ -0,0 +1,151 @@+name:           fpnla-examples+version:        0.1+cabal-version:  >=1.8+build-type:     Simple+stability:      Experimental+author:         Pablo Perdomo, Mauro Blanco+synopsis:       Example implementations for FPNLA library+description:+    .+    This package provides many example implementations of the operations defined in package @fpnla@ using various data structures and algorithms.+    .+    IMPORTANT: Many algorithms are only partially implemented+    .+    The structures currently used are:+    .+    - @hmatrix@ <http://hackage.haskell.org/package/hmatrix>+    .+    - @vector@ <http://hackage.haskell.org/package/vector>+    .+    - @array@ <http://hackage.haskell.org/package/array>+    .+    - @repa arrays@ <http://hackage.haskell.org/package/repa>+    .+    - @accelerate arrays@ <http://hackage.haskell.org/package/accelerate>+    .+    - @array@ <http://hackage.haskell.org/package/array>+    .+    - Haskell lists <http://hackage.haskell.org/package/base-4.6.0.1/docs/Data-List.html>+    .+    Operations are parallelized using the @parallel@ framework (<http://hackage.haskell.org/package/parallel>), @Monad-Par@ (<http://hackage.haskell.org/package/monad-par>), @Repa@ (<http://hackage.haskell.org/package/repa>) and @Accelerate@ (<http://hackage.haskell.org/package/accelerate>). +    .+tested-with:    GHC==7.6.3+maintainer:     Pablo Perdomo <pperdomo@fing.edu.uy>, Mauro Blanco <mblanco@fing.edu.uy>+license-file:   LICENSE+license:        BSD3+category:       Math, Parallelism++library+  hs-source-dirs:   src+  build-depends:    +                   base >= 4 && < 5,+                   deepseq >= 1.3.0.1,+                   accelerate >= 0.13.0.5,+                   fpnla,+                   array >= 0.4.0.1,+                   parallel >= 3.2.0.3,+                   hmatrix >= 0.14.1.0,+                   repa >= 3.2.3.3,+                   vector >= 0.10.0.1,+                   linear-algebra-cblas,+                   monad-par >= 0.3.4.2+  ghc-options:      -Wall -fno-warn-orphans -fno-warn-name-shadowing+  exposed-modules:  +                   FPNLA.Matrix.Instances.AccMatrix,+                   FPNLA.Matrix.Instances.ArrayMatrix,+                   FPNLA.Matrix.Instances.HMatrix,+                   FPNLA.Matrix.Instances.LLMatrixByCols,+                   FPNLA.Matrix.Instances.LLMatrixByRows,+                   FPNLA.Matrix.Instances.RepaMatrix,+                   FPNLA.Matrix.Instances.VectorMatrix,+                   FPNLA.Operations.BLAS.Strategies,+                   FPNLA.Operations.BLAS.Strategies.DOT.DefSeq,+                   FPNLA.Operations.BLAS.Strategies.DataTypes,+                   FPNLA.Operations.BLAS.Strategies.GEMM.Accelerate.DefPar,+                   FPNLA.Operations.BLAS.Strategies.GEMM.BlocksSeq,+                   FPNLA.Operations.BLAS.Strategies.GEMM.CBindSeq,+                   FPNLA.Operations.BLAS.Strategies.GEMM.DefSeq,+                   FPNLA.Operations.BLAS.Strategies.GEMM.MonadPar.BlocksPar,+                   FPNLA.Operations.BLAS.Strategies.GEMM.MonadPar.DefPar,+                   FPNLA.Operations.BLAS.Strategies.GEMM.MonadPar.StrassenMPar,+                   FPNLA.Operations.BLAS.Strategies.GEMM.MonadPar.StrassenPar,+                   FPNLA.Operations.BLAS.Strategies.GEMM.Strategies.BlocksPar,+                   FPNLA.Operations.BLAS.Strategies.GEMM.Strategies.DefPar,+                   FPNLA.Operations.BLAS.Strategies.GEMM.Strategies.StrassenPar,+                   FPNLA.Operations.BLAS.Strategies.GEMV,+                   FPNLA.Operations.BLAS.Strategies.SYRK.DefSeq,+                   FPNLA.Operations.BLAS.Strategies.SYRK.MonadPar.DefPar,+                   FPNLA.Operations.BLAS.Strategies.SYRK.Strategies.DefPar,+                   FPNLA.Operations.BLAS.Strategies.TRSM.BlocksByCols,+                   FPNLA.Operations.BLAS.Strategies.TRSM.BlocksByRows,+                   FPNLA.Operations.BLAS.Strategies.TRSM.CBindSeq,+                   FPNLA.Operations.BLAS.Strategies.TRSM.DefSeq,+                   FPNLA.Operations.BLAS.Strategies.TRSM.MonadPar.ColsPar,+                   FPNLA.Operations.BLAS.Strategies.TRSM.Strategies.ColsPar,+                   FPNLA.Operations.LAPACK.Strategies,+                   FPNLA.Operations.LAPACK.Strategies.DataTypes,+                   FPNLA.Operations.LAPACK.Strategies.POTRF.BlocksSeq,+                   FPNLA.Operations.LAPACK.Strategies.POTRF.DefSeq,+                   FPNLA.Operations.LAPACK.Strategies.POTRF.HMatrixBindSeq,+                   FPNLA.Operations.LAPACK.Strategies.POTRF.Repa.BlocksPar,+                   FPNLA.Operations.LAPACK.Strategies.POTRF.Repa.DefPar,+                   FPNLA.Operations.Utils,+                   FPNLA.Utils++test-suite test-fpnla-blas+  type:            exitcode-stdio-1.0+  ghc-options:     -Wall -rtsopts -fno-warn-orphans+  build-depends:   base >= 4, +                   HUnit >= 1.2 && < 2, +                   QuickCheck >= 2.4, +                   test-framework >= 0.4.1, +                   test-framework-quickcheck2, +                   test-framework-hunit,+                   data-default,+                   tagged,+                   time,+                   random,+                   fpnla,+                   fpnla-examples,  +                   deepseq,+                   accelerate,+                   array,+                   parallel,+                   hmatrix,+                   repa,+                   vector,+                   linear-algebra-cblas,+                   monad-par+  hs-source-dirs:  src-test+  main-is:         TestBLAS.hs+  extra-lib-dirs:  /usr/lib+  other-modules:   TestBase++test-suite test-fpnla-lapack+  type:            exitcode-stdio-1.0+  ghc-options:     -Wall -rtsopts -fno-warn-orphans+  build-depends:   base >= 4, +                   HUnit >= 1.2 && < 2, +                   QuickCheck >= 2.4, +                   test-framework >= 0.4.1, +                   test-framework-quickcheck2, +                   test-framework-hunit,+                   data-default,+                   tagged,+                   time,+                   random,+                   fpnla,+                   fpnla-examples,  +                   deepseq,+                   accelerate,+                   array,+                   parallel,+                   hmatrix,+                   repa,+                   vector,+                   linear-algebra-cblas,+                   monad-par+  hs-source-dirs:  src-test+  main-is:         TestLAPACK.hs+  extra-lib-dirs:  /usr/lib+  other-modules:   TestBase
+ src-test/TestBLAS.hs view
@@ -0,0 +1,134 @@+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE TypeSynonymInstances #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE OverlappingInstances #-}++module Main(main) where++import Test.QuickCheck (Arbitrary,quickCheckResult)+import Data.Tagged (Tagged(Tagged, unTagged))+import Data.Default (Default(def))++import qualified Data.Packed.Matrix as HM (Matrix)+import qualified Data.Packed.Vector as HV (Vector)++import TestBase (Truncable, GEMMTestData(..), StrassenTestData(..), TRSMTestData(..), Tag2(), eqMatrix, handleResult)++import FPNLA.Matrix.Instances.HMatrix ()+import FPNLA.Matrix.Instances.AccMatrix (AccMatrix, AccVector)+import Data.Array.Repa.Repr.Unboxed ()+import FPNLA.Operations.BLAS (GEMM(gemm), TRSM(trsm))+import FPNLA.Operations.Parameters (StratCtx, ResM, getResultDataM)+import FPNLA.Operations.BLAS.Strategies (CBindSeq, DefSeq, BlocksSeq, BlocksByRows, BlocksByCols, DefPar_MP, BlocksPar_MP, StrassenPar_MP, StrassenMPar_MP, ColsPar_MP, DefPar_ST, DefPar_Acc, BlocksPar_ST, StrassenPar_ST, ColsPar_ST)++import qualified Data.Array.Accelerate as A (+    IsNum, Elt)+++-- Versión general de parar realizar pruebas de GEMM:+propGEMMStratIsOK :: forall st m v e. (Show e, Eq e, Truncable e, Arbitrary e, GEMM st m v e, Default (Tagged st (StratCtx st)), GEMM CBindSeq m v e) => GEMMTestData m e -> Tagged (Tag2 st v) Bool+propGEMMStratIsOK (GEMMTestData pmA pmB alpha beta mC) = Tagged $ eqMatrix r1 r2+    where+        r1 = getResultDataM (gemm (unTagged (def :: Tagged st (StratCtx st))) pmA pmB alpha beta mC :: ResM st v m e)+        r2 = getResultDataM (gemm (unTagged (def :: Tagged CBindSeq (StratCtx CBindSeq))) pmA pmB alpha beta mC :: ResM CBindSeq v m e)++propGEMMAccStratIsOK :: forall st e. (A.Elt e, A.IsNum e, Eq e, Truncable e, Arbitrary e, GEMM st AccMatrix AccVector e, Default (Tagged st (StratCtx st)), GEMM CBindSeq AccMatrix AccVector e) => GEMMTestData AccMatrix e -> Tagged (Tag2 st AccVector) Bool+propGEMMAccStratIsOK (GEMMTestData pmA pmB alpha beta mC) = Tagged $ eqMatrix r1 r2+    where+        r1 = getResultDataM (gemm (unTagged (def :: Tagged st (StratCtx st))) pmA pmB alpha beta mC :: ResM st AccVector AccMatrix e)+        r2 = getResultDataM (gemm (unTagged (def :: Tagged CBindSeq (StratCtx CBindSeq))) pmA pmB alpha beta mC :: ResM CBindSeq AccVector AccMatrix e)+++propStrassenStratIsOK :: forall st m v e. (Show e, Eq e, Truncable e, Arbitrary e, GEMM st m v e, Default (Tagged st (StratCtx st)), GEMM CBindSeq m v e) => StrassenTestData m e -> Tagged (Tag2 st v) Bool+propStrassenStratIsOK (StrassenTestData pmA pmB alpha beta mC) = Tagged $ eqMatrix r1 r2+    where+        r1 = getResultDataM (gemm (unTagged (def :: Tagged st (StratCtx st))) pmA pmB alpha beta mC :: ResM st v m e)+        r2 = getResultDataM (gemm (unTagged (def :: Tagged CBindSeq (StratCtx CBindSeq))) pmA pmB alpha beta mC :: ResM CBindSeq v m e)++propTRSMStratIsOK :: forall st m v e. (Show e, Eq e, Truncable e, Arbitrary e, TRSM st m v e, Default (Tagged st (StratCtx st)), TRSM CBindSeq m v e) => TRSMTestData m e -> Tagged (Tag2 st v) Bool+propTRSMStratIsOK (TRSMTestData alpha mA mB) = Tagged $ eqMatrix r1 r2+    where+        r1 = getResultDataM (trsm (unTagged (def :: Tagged st (StratCtx st))) alpha mA mB :: ResM st v m e)+        r2 = getResultDataM (trsm (unTagged (def :: Tagged CBindSeq (StratCtx CBindSeq))) alpha mA mB :: ResM CBindSeq v m e)++--+type TMatrix = {-VMatrix-} {--VMatrix--} HM.Matrix {-LLR.LLMatrix-} {-RepaMatrix-}+type TVector = {-VVector-} {--A.Array Int--} HV.Vector {-LLR.LVector-} {-RepaVector-}+type TElem = Double++propDefSeqStratIsOK_Inst :: GEMMTestData TMatrix TElem -> Tagged (Tag2 DefSeq TVector) Bool+propDefSeqStratIsOK_Inst = propGEMMStratIsOK+propBlocksSeqStratIsOK_Inst :: GEMMTestData TMatrix TElem -> Tagged (Tag2 (BlocksSeq DefSeq) TVector) Bool+propBlocksSeqStratIsOK_Inst = propGEMMStratIsOK+propTRSMStratIsOK_Inst :: TRSMTestData TMatrix TElem -> Tagged (Tag2 DefSeq TVector) Bool+propTRSMStratIsOK_Inst = propTRSMStratIsOK+propTRSMBlockByColsIsOK_Inst :: TRSMTestData TMatrix TElem -> Tagged (Tag2 (BlocksByCols DefSeq DefSeq) TVector) Bool+propTRSMBlockByColsIsOK_Inst = propTRSMStratIsOK+propTRSMBlockByRowsIsOK_Inst :: TRSMTestData TMatrix TElem -> Tagged (Tag2 (BlocksByRows DefSeq DefSeq) TVector) Bool+propTRSMBlockByRowsIsOK_Inst = propTRSMStratIsOK++propDefParMPStratIsOK_Inst :: GEMMTestData TMatrix TElem -> Tagged (Tag2 DefPar_MP TVector) Bool+propDefParMPStratIsOK_Inst = propGEMMStratIsOK+propBlocksParMPStratIsOK_Inst :: GEMMTestData TMatrix TElem -> Tagged (Tag2 (BlocksPar_MP DefSeq) TVector) Bool+propBlocksParMPStratIsOK_Inst = propGEMMStratIsOK+propStrassenMPStratIsOK_Inst :: StrassenTestData TMatrix TElem -> Tagged (Tag2 (StrassenPar_MP DefSeq) TVector) Bool+propStrassenMPStratIsOK_Inst = propStrassenStratIsOK+propStrassenMMPStratIsOK_Inst :: StrassenTestData TMatrix TElem -> Tagged (Tag2 (StrassenMPar_MP DefSeq) TVector) Bool+propStrassenMMPStratIsOK_Inst = propStrassenStratIsOK+propTRSMMPStratIsOK_Inst :: TRSMTestData TMatrix TElem -> Tagged (Tag2 (ColsPar_MP DefSeq) TVector) Bool+propTRSMMPStratIsOK_Inst = propTRSMStratIsOK++propDefParSTStratIsOK_Inst :: GEMMTestData TMatrix TElem -> Tagged (Tag2 DefPar_ST TVector) Bool+propDefParSTStratIsOK_Inst = propGEMMStratIsOK+propBlocksParSTStratIsOK_Inst :: GEMMTestData TMatrix TElem -> Tagged (Tag2 (BlocksPar_ST DefSeq) TVector) Bool+propBlocksParSTStratIsOK_Inst = propGEMMStratIsOK+propStrassenSTStratIsOK_Inst :: StrassenTestData TMatrix TElem -> Tagged (Tag2 (StrassenPar_ST DefSeq) TVector) Bool+propStrassenSTStratIsOK_Inst = propStrassenStratIsOK+propTRSMSTStratIsOK_Inst :: TRSMTestData TMatrix TElem -> Tagged (Tag2 (ColsPar_ST DefSeq) TVector) Bool+propTRSMSTStratIsOK_Inst = propTRSMStratIsOK++propDefParAccStratIsOK_Inst :: GEMMTestData AccMatrix TElem -> Tagged (Tag2 DefPar_Acc AccVector) Bool+propDefParAccStratIsOK_Inst = propGEMMAccStratIsOK++main :: IO ()+main = do+    putStrLn "execGEMMQC"+    res1 <- quickCheckResult propDefSeqStratIsOK_Inst+    res2 <- quickCheckResult propBlocksSeqStratIsOK_Inst+    res3 <- quickCheckResult propDefParMPStratIsOK_Inst+    res4 <- quickCheckResult propBlocksParMPStratIsOK_Inst+    res5 <- quickCheckResult propDefParSTStratIsOK_Inst+    res6 <- quickCheckResult propBlocksParSTStratIsOK_Inst+    putStrLn "execGEMMAccQC"+    --res7 <- quickCheckResult propDefParAccStratIsOK_Inst+    putStrLn "execStrassenQC"+    res8 <- quickCheckResult propStrassenMPStratIsOK_Inst+    res9 <- quickCheckResult propStrassenMMPStratIsOK_Inst+    res10 <- quickCheckResult propStrassenSTStratIsOK_Inst+    putStrLn "execTRSMQC"+    res11 <- quickCheckResult propTRSMStratIsOK_Inst+    res12 <- quickCheckResult propTRSMBlockByColsIsOK_Inst+    res13 <- quickCheckResult propTRSMBlockByRowsIsOK_Inst+    res14 <- quickCheckResult propTRSMMPStratIsOK_Inst+    res15 <- quickCheckResult propTRSMSTStratIsOK_Inst+    handleResult res1+    handleResult res2+    handleResult res3+    handleResult res4+    handleResult res5+    handleResult res6+    --handleResult res7+    handleResult res8+    handleResult res9+    handleResult res10+    handleResult res11+    handleResult res12+    handleResult res13+    handleResult res14+    handleResult res15+    +    +
+ src-test/TestBase.hs view
@@ -0,0 +1,239 @@+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE TypeSynonymInstances #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE OverlappingInstances #-}+{-# LANGUAGE EmptyDataDecls #-}++module TestBase(+    GEMMTestData(..), +    StrassenTestData(..), +    TRSMTestData(..), +    PotrfTestData(..), +    Tag2(), +    eqMatrix, +    handleResult,+    Truncable(..)+    +) where++import           Data.Complex (Complex ((:+)))+import Test.QuickCheck (Result(..), Arbitrary(arbitrary), Testable(property), Gen, choose, vectorOf, elements)+import Data.Tagged (Tagged(Tagged, unTagged))+import Data.Default (Default(def))+import System.Random (Random)+import Debug.Trace (trace)+import Control.Applicative (liftA)++import FPNLA.Matrix+       (Matrix(fromList_m, dim_m, elem_m),+        Matrix(generate_m, transpose_m, zipWith_m))+import FPNLA.Matrix.Instances.HMatrix ()+import FPNLA.Operations.Parameters (Elt, TriangType(..), TransType(..), UnitType(..), StratCtx)+import FPNLA.Operations.LAPACK.Strategies (HMatrixBindSeq, CholLLVSeq, CholLLVPar_Repa, CholLLVBlocksSeq, newNullContext, newSqrBlockContext)+import FPNLA.Operations.BLAS.Strategies (DefPar_ST, CBindSeq, DefSeq, BlocksSeq, BlocksByRows, BlocksByCols, DefPar_MP, BlocksPar_MP, StrassenPar_MP, StrassenMPar_MP, ColsPar_MP, DefPar_ST, DefPar_Acc, BlocksPar_ST, StrassenPar_ST, ColsPar_ST, newBlockContext, newStrassenContext)+import FPNLA.Utils (iif)++-- Links interesantes:+-- http://www.haskell.org/haskellwiki/Scoped_type_variables+-- http://www.haskell.org/ghc/docs/7.4.1/html/users_guide/other-type-extensions.html seccion 7.11.7+-- http://hackage.haskell.org/trac/ghc/ticket/3485++blockSize :: Int+blockSize = 2++-- Valor de epsilon para hacer las comparaciones:+epsilon :: Double+epsilon = 1e-5++class Truncable a where+    inEpsilonRange :: Real e => e -> a -> a -> Bool++instance (Floating n, RealFrac n) => Truncable n where+    inEpsilonRange eps n1 n2 = expN1 == expN2 && abs(mantN1 - mantN2) <= epsFractional+        where+            epsFractional = fromRational $ toRational eps+            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)++--+instance (RealFrac n, Truncable n) => Truncable (Complex n) where+    inEpsilonRange eps (r1:+i1) (r2:+i2) = inEpsilonRange eps r1 r2 && inEpsilonRange eps i1 i2+++-- Compara una matriz elemento a elemento truncando los dígitos. Si falla imprime los índices del primer elemento que falló:+eqMatrix :: (Show e, Eq e, Truncable e, Matrix m e) => m e -> m e -> Bool+eqMatrix m1 m2 = dimm1 == dimm2 && eqElems 0 0+    where+        dimm1 = dim_m m1+        dimm2 = dim_m m2+        (m,n) = dimm1+        getElem = {--truncDigits truncValue $--} elem_m+        eqElems i j+          | j >= n = eqElems (i+1) 0+          | i >= m = True+          -- | getElem i j m1 /= getElem i j m2 = trace ("(i,j) = (" ++ show i ++ "," ++ show j ++ ")") False+          | not $ inEpsilonRange epsilon (getElem i j m1) (getElem i j m2) = trace ("(i,j) = (" ++ show i ++ "," ++ show j ++ ")" ++ " / m1 = " ++ show (getElem i j m1) ++ " / m2 = " ++ show (getElem i j m2)) False+          | otherwise =  eqElems i (j+1)++arbSDPMatrix :: (Elt e, Arbitrary e, Random e, Matrix m e) => Int -> Gen (m e)+arbSDPMatrix n = do+    la <- vectorOf (n*n) $ choose (0, 1)+    ld <- vectorOf n $ choose (100, 200)+    let a = fromList_m n n la+        d = generate_m n n (\i j -> iif (i == j) (ld !! i) 0)+        s = zipWith_m (+) a (transpose_m a)+        r = zipWith_m (+) s d+    return r+++handleResult :: Monad m => Result -> m ()+handleResult (Success{}) = return ()+handleResult (GaveUp{}) = fail "GaveUp"+handleResult (Failure{}) = fail "Failure"+handleResult (NoExpectedFailure{}) = fail "NoExpectedFailure"++--+instance Arbitrary (TransType a) where+    arbitrary = elements [NoTrans undefined, Trans undefined, ConjTrans undefined]+--+instance Arbitrary (a -> TransType a) where+    arbitrary = elements [NoTrans, Trans, ConjTrans]+instance Arbitrary (a -> TriangType a) where+    arbitrary = elements [Lower, Upper]+instance Arbitrary (a -> UnitType a) where+    arbitrary = elements [Unit, NoUnit]++--+arbitraryIntBounded :: Gen Int+arbitraryIntBounded = choose (2, 10)++arbitraryMatrix :: (Arbitrary e, Matrix m e) => Int -> Int -> Gen (m e)+arbitraryMatrix m n = do+    l <- vectorOf (m*n) arbitrary+    return $ fromList_m m n l++arbitraryTransMatrix :: (Arbitrary e, Matrix m e) => TransType a -> Int -> Int -> Gen (TransType (m e))+arbitraryTransMatrix (NoTrans _) m n = liftA NoTrans $ arbitraryMatrix m n+arbitraryTransMatrix (Trans _) m n = liftA Trans $ arbitraryMatrix n m+arbitraryTransMatrix (ConjTrans _) m n = liftA ConjTrans $ arbitraryMatrix n m+++data GEMMTestData m e = GEMMTestData (TransType (m e)) (TransType (m e)) e e (m e) deriving (Show)+data StrassenTestData m e = StrassenTestData (TransType (m e)) (TransType (m e)) e e (m e) deriving (Show)+data TRSMTestData m e = TRSMTestData e (TransType (TriangType (UnitType (m e)))) (m e) deriving (Show)++instance (Matrix m e, Arbitrary e) => Arbitrary (GEMMTestData m e) where+    arbitrary = do+        m <- arbitraryIntBounded+        n <- arbitraryIntBounded+        k <- arbitraryIntBounded+        ttmA <- arbitrary+        ttmB <- arbitrary+        pmA <- arbitraryTransMatrix ttmA m k+        pmB <- arbitraryTransMatrix ttmB k n+        alpha <- arbitrary+        beta <- arbitrary+        mC <- arbitraryMatrix m n+        return $ GEMMTestData pmA pmB alpha beta mC++instance (Matrix m e, Arbitrary e) => Arbitrary (StrassenTestData m e) where+    arbitrary = do+        e <- choose (2, 4) :: Gen Int+        let n = truncate $ (2 :: Double) ** fromIntegral e+        ttmA <- arbitrary+        ttmB <- arbitrary+        alpha <- arbitrary+        mA <- arbitraryTransMatrix ttmA n n+        mB <- arbitraryTransMatrix ttmB n n+        beta <- arbitrary+        mC <- arbitraryMatrix n n+        return $ StrassenTestData mA mB alpha beta mC+--+instance (Matrix m e, Arbitrary e) => Arbitrary (TRSMTestData m e) where+    arbitrary = do+        m <- arbitraryIntBounded+        n <- arbitraryIntBounded+        trans <- arbitrary :: Gen (a -> TransType a)+        triang <- arbitrary :: Gen (a -> TriangType a)+        unit <- arbitrary :: Gen (a -> UnitType a)+        alpha <- arbitrary+        mA <- arbitraryMatrix m m+        mB <- arbitraryMatrix m n+        return $ TRSMTestData alpha (trans . triang . unit $ mA) mB++data PotrfTestData m e = PotrfTestData (TriangType (m e)) deriving (Show)++instance (Matrix m e, Elt e, Arbitrary e, Random e) => Arbitrary (PotrfTestData m e) where+    arbitrary = do+        _ <- arbitrary :: Gen (a -> TriangType a)+        n <- choose (3, 5)+        m <- arbSDPMatrix (blockSize ^ (n :: Int))+        return $ PotrfTestData $ Lower m++-- Cualquier propiedad tageada es una propiedad:+instance (Testable p) => Testable (Tagged t p) where+    property p = property $ unTagged p+    +data Tag2 (s :: *) (v :: * -> *)++instance (StratCtx HMatrixBindSeq ~ ctx_Seq) => Default (Tagged HMatrixBindSeq ctx_Seq) where+    def = Tagged newNullContext++instance (StratCtx dots ~ ctx_dots, Default (Tagged dots ctx_dots), StratCtx gs ~ ctx_gs, Default (Tagged gs ctx_gs), StratCtx (CholLLVSeq dots gs) ~ ctx_bp) => Default (Tagged (CholLLVSeq dots gs) ctx_bp) where+    def = Tagged (unTagged (def :: Tagged dots ctx_dots), unTagged (def :: Tagged gs ctx_gs))++instance (StratCtx dots ~ ctx_dots, Default (Tagged dots ctx_dots), StratCtx gs ~ ctx_gs, Default (Tagged gs ctx_gs), StratCtx (CholLLVPar_Repa dots gs) ~ ctx_bp) => Default (Tagged (CholLLVPar_Repa dots gs) ctx_bp) where+    def = Tagged (unTagged (def :: Tagged dots ctx_dots), unTagged (def :: Tagged gs ctx_gs))++instance (StratCtx ps ~ ctx_ps, Default (Tagged ps ctx_ps), StratCtx ss ~ ctx_ss, Default (Tagged ss ctx_ss), StratCtx gs ~ ctx_gs, Default (Tagged gs ctx_gs), StratCtx ts ~ ctx_ts, Default (Tagged ts ctx_ts), StratCtx (CholLLVBlocksSeq ps ss gs ts) ~ ctx_cs) => Default (Tagged (CholLLVBlocksSeq ps ss gs ts) ctx_cs) where+    def = Tagged (newSqrBlockContext blockSize, unTagged (def :: Tagged ps ctx_ps), unTagged (def :: Tagged ss ctx_ss), unTagged (def :: Tagged gs ctx_gs), unTagged (def :: Tagged ts ctx_ts))+++--+instance (StratCtx DefSeq ~ ctx_DefSeq) => Default (Tagged DefSeq ctx_DefSeq) where+    def = Tagged newNullContext+instance (StratCtx CBindSeq ~ ctx_DefSeq) => Default (Tagged CBindSeq ctx_DefSeq) where+    def = Tagged newNullContext+instance (StratCtx gs ~ ctx_gs, Default (Tagged gs ctx_gs), StratCtx (BlocksSeq gs) ~ ctx_bp) => Default (Tagged (BlocksSeq gs) ctx_bp) where+    def = Tagged (newBlockContext 2 2, unTagged (def :: Tagged gs ctx_gs))+instance (StratCtx gs ~ ctx_gs, Default (Tagged gs ctx_gs), StratCtx ts ~ ctx_ts, Default (Tagged ts ctx_ts), StratCtx (BlocksByCols gs ts) ~ ctx_bbcp) => Default (Tagged (BlocksByCols gs ts) ctx_bbcp) where+    def = Tagged (newSqrBlockContext 2, unTagged (def :: Tagged gs ctx_gs), unTagged (def :: Tagged ts ctx_ts))+instance (StratCtx gs ~ ctx_gs, Default (Tagged gs ctx_gs), StratCtx ts ~ ctx_ts, Default (Tagged ts ctx_ts), StratCtx (BlocksByRows gs ts) ~ ctx_bbcp) => Default (Tagged (BlocksByRows gs ts) ctx_bbcp) where+    def = Tagged (newSqrBlockContext 2, unTagged (def :: Tagged gs ctx_gs), unTagged (def :: Tagged ts ctx_ts))++instance (StratCtx DefPar_MP ~ ctx_DefPar) => Default (Tagged DefPar_MP ctx_DefPar) where+    def = Tagged newNullContext+instance (StratCtx gs ~ ctx_gs, Default (Tagged gs ctx_gs), StratCtx (BlocksPar_MP gs) ~ ctx_bp) => Default (Tagged (BlocksPar_MP gs) ctx_bp) where+    def = Tagged (newBlockContext 2 2, unTagged (def :: Tagged gs ctx_gs))+instance (StratCtx ts ~ ctx_ts, Default (Tagged ts ctx_ts), StratCtx (ColsPar_MP ts) ~ ctx_cp) => Default (Tagged (ColsPar_MP ts) ctx_cp) where+    def = Tagged $ unTagged (def :: Tagged ts ctx_ts)+instance (StratCtx gs ~ ctx_gs, Default (Tagged gs ctx_gs), StratCtx (StrassenPar_MP gs) ~ ctx_bp) => Default (Tagged (StrassenPar_MP gs) ctx_bp) where+   def = Tagged (newStrassenContext 2, unTagged (def :: Tagged gs ctx_gs))+instance (StratCtx gs ~ ctx_gs, Default (Tagged gs ctx_gs), StratCtx (StrassenMPar_MP gs) ~ ctx_bp) => Default (Tagged (StrassenMPar_MP gs) ctx_bp) where+   def = Tagged (newStrassenContext 2, unTagged (def :: Tagged gs ctx_gs))++instance (StratCtx DefPar_ST ~ ctx_DefPar) => Default (Tagged DefPar_ST ctx_DefPar) where+    def = Tagged newNullContext+instance (StratCtx gs ~ ctx_gs, Default (Tagged gs ctx_gs), StratCtx (BlocksPar_ST gs) ~ ctx_bp) => Default (Tagged (BlocksPar_ST gs) ctx_bp) where+    def = Tagged (newBlockContext 2 2, unTagged (def :: Tagged gs ctx_gs))+instance (StratCtx ts ~ ctx_ts, Default (Tagged ts ctx_ts), StratCtx (ColsPar_ST ts) ~ ctx_cp) => Default (Tagged (ColsPar_ST ts) ctx_cp) where+    def = Tagged $ unTagged (def :: Tagged ts ctx_ts)+instance (StratCtx gs ~ ctx_gs, Default (Tagged gs ctx_gs), StratCtx (StrassenPar_ST gs) ~ ctx_bp) => Default (Tagged (StrassenPar_ST gs) ctx_bp) where+   def = Tagged (newStrassenContext 2, unTagged (def :: Tagged gs ctx_gs))++instance (StratCtx DefPar_Acc ~ ctx_DefPar) => Default (Tagged DefPar_Acc ctx_DefPar) where+    def = Tagged newNullContext++
+ src-test/TestLAPACK.hs view
@@ -0,0 +1,72 @@+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE TypeSynonymInstances #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE OverlappingInstances #-}++module Main (++    main,++) where++import Test.QuickCheck (Arbitrary, quickCheckResult)+import Data.Tagged (Tagged(Tagged, unTagged))+import Data.Default (Default(def))+import System.Random (Random)++--import qualified Data.Array as A (Array)+--import qualified Data.Array.Repa as R (D, DIM1, DIM2, Array)+import qualified Data.Packed.Matrix as HM (Matrix)+import qualified Data.Packed.Vector as HV (Vector)++import TestBase (Truncable, PotrfTestData(..), Tag2(), eqMatrix, handleResult)++import FPNLA.Matrix.Instances.RepaMatrix(RepaMatrix, RepaVector)+import FPNLA.Matrix.Instances.HMatrix ()+import FPNLA.Operations.Parameters (ResM, getResultDataM, StratCtx)+import FPNLA.Operations.LAPACK.Strategies (HMatrixBindSeq, CholLLVSeq, CholLLVPar_Repa, CholLLVBlocksSeq)+import FPNLA.Operations.BLAS.Strategies (DefSeq, DefPar_ST)+import FPNLA.Operations.LAPACK (POTRF(potrf))++-- Propiedades:+propPotrfIsOK :: forall st m v e. (Show e, Eq e, Truncable e, Arbitrary e, Random e, POTRF st m v e, Default (Tagged st (StratCtx st)), POTRF HMatrixBindSeq m v e) => PotrfTestData m e -> Tagged (Tag2 st v) Bool+propPotrfIsOK (PotrfTestData pmA) = Tagged $ eqMatrix r1 r2+    where+        r1 = getResultDataM (potrf (unTagged (def :: Tagged st (StratCtx st))) pmA :: ResM st v m e)+        r2 = getResultDataM (potrf (unTagged (def :: Tagged HMatrixBindSeq (StratCtx HMatrixBindSeq))) pmA :: ResM HMatrixBindSeq v m e)++type TMatrix = {-VMatrix-} {--VMatrix--} HM.Matrix {-LLR.LLMatrix-} {-RepaMatrix-}+type TVector = {-VVector-} {--A.Array Int--} HV.Vector {-LLR.LVector-} {-RepaVector-}+--type TElem = Complex Double+type TElem = Double++propTest :: PotrfTestData HM.Matrix Double -> Tagged (Tag2 HMatrixBindSeq HV.Vector) Bool+propTest = propPotrfIsOK++propPotrfDefSeqStratIsOK_Inst :: PotrfTestData TMatrix TElem -> Tagged (Tag2 (CholLLVSeq DefSeq DefSeq) TVector) Bool+propPotrfDefSeqStratIsOK_Inst = propPotrfIsOK++propPotrfRepaDefSeqStratIsOK_Inst :: PotrfTestData RepaMatrix TElem -> Tagged (Tag2 (CholLLVPar_Repa DefSeq DefSeq) RepaVector) Bool+propPotrfRepaDefSeqStratIsOK_Inst = propPotrfIsOK++propPotrfCholLLVBlocksSeqStratIsOK_Inst :: PotrfTestData TMatrix TElem -> Tagged (Tag2 (CholLLVBlocksSeq DefPar_ST DefSeq DefSeq (CholLLVSeq DefSeq DefSeq)) TVector) Bool+propPotrfCholLLVBlocksSeqStratIsOK_Inst = propPotrfIsOK++execPotrfQC :: IO [()]+execPotrfQC = do+    putStrLn "execPotrfQC"+    res1 <- quickCheckResult propTest+    res2 <- quickCheckResult propPotrfDefSeqStratIsOK_Inst+    res3 <- quickCheckResult propPotrfRepaDefSeqStratIsOK_Inst+    res4 <- quickCheckResult propPotrfCholLLVBlocksSeqStratIsOK_Inst+    mapM handleResult [res1, res2, res3, res4]++execQuickCheck :: IO [()]+execQuickCheck = execPotrfQC++main :: IO [()]+main = execQuickCheck+
+ src/FPNLA/Matrix/Instances/AccMatrix.hs view
@@ -0,0 +1,68 @@+{-# LANGUAGE FlexibleInstances     #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE OverlappingInstances  #-}+{-# LANGUAGE TypeOperators         #-}++module FPNLA.Matrix.Instances.AccMatrix (+    AccMatrix(..),+    AccVector(..),+) where++import           FPNLA.Matrix                      (Matrix (generate_m, fromList_m, transpose_m, dim_m, elem_m, map_m, zipWith_m, subMatrix_m, fromBlocks_m, toBlocks_m), MatrixVector (row_vm, col_vm, fromCols_vm, toCols_vm), Vector (generate_v, fromList_v, concat_v, elem_v, length_v, foldr_v, map_v, zipWith_v))++import           Data.Array.Accelerate             as A ((:.) (..), Acc (),+                                                         Array (), DIM1 (),+                                                         DIM2 (), Elt (),+                                                         IsNum (), Z (Z),+                                                         arrayShape, foldAll,+                                                         fromList, indexArray,+                                                         use)+import           Data.Array.Accelerate.Interpreter as R (run)++import           Control.DeepSeq                   (NFData (rnf))++import           Debug.Trace                       (trace)++newtype AccVector e = AccVector (Acc (Array DIM1 e))+newtype AccMatrix e = AccMatrix (Acc (Array DIM2 e))++instance (NFData e, Elt e, IsNum e) => (NFData (AccMatrix e)) where+    rnf (AccMatrix m) = rnf $ seq (R.run $ foldAll (\_ _ -> 0) 0 m) ()++instance (NFData e, Elt e, IsNum e) => (NFData (AccVector e)) where+    rnf (AccVector v) = rnf $ seq (R.run $ foldAll (\_ _ -> 0) 0 v) ()++instance (Elt e) => Show (AccMatrix e) where+    show (AccMatrix m) = show m++instance (Elt e) => Show (AccVector e) where+    show (AccVector m) = show m++instance (Elt e) => Vector AccVector e where+    generate_v size gen = trace "generate_v" undefined+    fromList_v l = AccVector . use $ fromList (Z:.(length l)) l+    concat_v l = trace "concat_v" undefined+    elem_v pos v = trace "elem_v" undefined+    length_v v = trace "length_v" undefined+    foldr_v cons zero v = trace "foldr_v" undefined+    map_v = trace "map_v" undefined+    zipWith_v = trace "zipWith_v" undefined++instance (Elt e) => Matrix AccMatrix e where+    generate_m rows cols gen = trace "generate_m" undefined+    fromList_m m n l = AccMatrix . use $ fromList (Z:.m:.n) l+    transpose_m = trace "transpose_m" undefined+    dim_m (AccMatrix m) = (rows, cols)+        where (Z :. rows :. cols) = arrayShape (R.run m)+    elem_m i j (AccMatrix m) = A.indexArray (R.run m) (Z:.i:.j)+    map_m = trace "map_m" undefined+    zipWith_m = trace "zipWith_m" undefined+    subMatrix_m posI posJ cantRows cantCols = trace "subMatrix_m" undefined+    fromBlocks_m = trace "fromBlocks_m" undefined+    toBlocks_m = trace "toBlocks_m" undefined++instance (Elt e) => MatrixVector AccMatrix AccVector e  where+    row_vm pos m = trace "row_vm" undefined+    col_vm pos m = trace "col_vm" undefined+    fromCols_vm = trace "fromCols_vm" undefined+    toCols_vm = trace "toCols_vm" undefined
+ src/FPNLA/Matrix/Instances/ArrayMatrix.hs view
@@ -0,0 +1,48 @@+{-# LANGUAGE FlexibleInstances     #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE OverlappingInstances  #-}++module FPNLA.Matrix.Instances.ArrayMatrix (+    ArrayMatrix(ArrayMatrix),+) where++import           Control.DeepSeq (NFData (rnf))+import           Data.Array      (Array, bounds, elems, listArray, (!))+import qualified Data.Foldable   (foldr)+import           FPNLA.Matrix    (Matrix (fromList_m, dim_m, elem_m, generate_m), Vector (fromList_v, elem_v, generate_v, concat_v, foldr_v, length_v))+import           FPNLA.Utils     (mapPair)++-- Las instancias de NFData para Array estan definidas en Control.DeepSeq+instance (NFData arr) => NFData (ArrayMatrix arr) where+   rnf (ArrayMatrix arr) = rnf arr+++instance Vector (Array Int) e where+    -- Los indices siempre comienzan en 0!+    generate_v n f = listArray (0, n-1) [f x | x <- [0..]]+    fromList_v l = listArray (0, length l -1) l+    concat_v arrs = fromList_v $ concatMap elems arrs+    elem_v = flip (!)+    length_v = (+1) . snd . bounds+    foldr_v = Data.Foldable.foldr+++instance Show e => Matrix (Array (Int,Int)) e where+    generate_m rs cs gen = fromList_m rs cs [gen i j | i <- [0..rs-1], j <- [0..cs-1]]+    fromList_m rs cs = listArray ((0,0), (rs-1, cs-1))+    dim_m = mapPair (+1) (+1) . snd . bounds+    elem_m i j = flip (!) (i,j)+++newtype ArrayMatrix e = ArrayMatrix (Array (Int,Int) e) deriving (Show)+wrapIn :: Array (Int,Int) e -> ArrayMatrix e+wrapIn = ArrayMatrix+wrapOut :: ArrayMatrix e -> Array (Int,Int) e+wrapOut (ArrayMatrix arr) = arr++instance Show e => Matrix (ArrayMatrix) e where+    elem_m i j = flip (!) (i,j) . wrapOut+    dim_m = mapPair (+1) (+1) . snd . bounds . wrapOut+    fromList_m rs cs ls = wrapIn $ listArray ((0,0), (rs-1, cs-1)) ls+    generate_m rs cs gen = wrapIn $ fromList_m rs cs [gen i j | i <- [0..rs-1], j <- [0..cs-1]]+
+ src/FPNLA/Matrix/Instances/HMatrix.hs view
@@ -0,0 +1,57 @@+{-# LANGUAGE FlexibleInstances     #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE OverlappingInstances  #-}++module FPNLA.Matrix.Instances.HMatrix () where++import           Control.DeepSeq             (NFData (rnf))+import           Control.Parallel.Strategies (rseq, withStrategy)+import qualified Data.Packed.Matrix          as HM (Element, Matrix,+                                                    buildMatrix, cols,+                                                    fromBlocks, fromColumns,+                                                    mapMatrix, rows, subMatrix,+                                                    toBlocksEvery, toColumns,+                                                    toRows, trans, (><), (@@>))+import qualified Data.Packed.Vector          as HV (Vector, buildVector, dim,+                                                    foldVector, fromList, join,+                                                    mapVector, zipVectorWith,+                                                    (@>))+import           FPNLA.Matrix                (Matrix (generate_m, fromList_m, transpose_m, dim_m, elem_m, subMatrix_m, fromBlocks_m, toBlocks_m, map_m), MatrixVector (fromCols_vm, toCols_vm, row_vm, col_vm), Vector (generate_v, fromList_v, concat_v, elem_v, length_v, foldr_v, map_v, zipWith_v))+{-+instance (NFData e) => NFData (HV.Vector e) where+    -- Asumo que Vector es estricto+    -- http://haskell.1045720.n5.nabble.com/NFData-instance-for-Numeric-LinearAlgebra-Matrix-td4265725.html+    rnf v = (withStrategy rseq v) `seq` ()+-}+instance (NFData e) => NFData (HM.Matrix e) where+    -- Matrix es estricto+    -- http://haskell.1045720.n5.nabble.com/NFData-instance-for-Numeric-LinearAlgebra-Matrix-td4265725.html+    rnf m = withStrategy rseq m `seq` ()++instance (HM.Element e) => Vector HV.Vector e where+    generate_v = HV.buildVector+    fromList_v = HV.fromList+    concat_v = HV.join+    elem_v pos v = v HV.@> pos+    length_v = HV.dim+    foldr_v = HV.foldVector+    map_v = HV.mapVector+    zipWith_v = HV.zipVectorWith++instance (HM.Element e) => Matrix HM.Matrix e where+    generate_m rows cols gen = HM.buildMatrix rows cols (uncurry gen)+    fromList_m = (HM.><)+    transpose_m = HM.trans+    dim_m m = (HM.rows m, HM.cols m)+    elem_m i j m = m HM.@@> (i, j)+    map_m = HM.mapMatrix+    --zipWith_m+    subMatrix_m posI posJ cantRows cantCols = HM.subMatrix (posI, posJ) (cantRows, cantCols)+    fromBlocks_m = HM.fromBlocks+    toBlocks_m = HM.toBlocksEvery++instance (HM.Element e) => MatrixVector HM.Matrix HV.Vector e where+    row_vm pos m = HM.toRows m !! pos+    col_vm pos m = HM.toColumns m !! pos+    fromCols_vm = HM.fromColumns+    toCols_vm = HM.toColumns
+ src/FPNLA/Matrix/Instances/LLMatrixByCols.hs view
@@ -0,0 +1,70 @@+{-# LANGUAGE FlexibleInstances     #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE OverlappingInstances  #-}+{-# LANGUAGE TypeFamilies          #-}++module FPNLA.Matrix.Instances.LLMatrixByCols (+    LVector(LVector),+    LLMatrix(LLMatrix),+) where++import           Control.DeepSeq             (NFData (rnf))+import           Control.Parallel.Strategies (rdeepseq, withStrategy)+import           Data.Foldable               (foldr')+import           Data.List                   (transpose)+import           FPNLA.Matrix                (Matrix (generate_m, transpose_m, dim_m, elem_m, map_m, zipWith_m, subMatrix_m, fromBlocks_m), MatrixVector (fromCols_vm, toCols_vm, row_vm, col_vm), Vector (generate_v, fromList_v, concat_v, elem_v, length_v, foldr_v, map_v, zipWith_v),+                                              cantCols_m)++data LVector e = LVector [e]+data LLMatrix e = LLMatrix [[e]]++getLVectorInner :: LVector t -> [t]+getLVectorInner (LVector l) = l+getLLMatrixInner :: LLMatrix t -> [[t]]+getLLMatrixInner (LLMatrix ll) = ll+newLVector :: [e] -> LVector e+newLVector = LVector+newLLMatrix :: [[e]] -> LLMatrix e+newLLMatrix = LLMatrix++instance Show e => Show (LLMatrix e) where+    show = show . getLLMatrixInner . transpose_m++instance (NFData e) => NFData (LVector e) where+    rnf v =+        withStrategy rdeepseq (getLVectorInner v) `seq` ()++instance (NFData e) => NFData (LLMatrix e) where+    rnf m =+        withStrategy rdeepseq (getLLMatrixInner m) `seq` ()++instance Vector LVector e where+    generate_v size gen = newLVector [gen x | x <- [0..(size - 1)]]+    fromList_v = newLVector+    concat_v = newLVector . concatMap getLVectorInner+    elem_v pos v = getLVectorInner v !! pos+    length_v = length . getLVectorInner+    foldr_v cons zero = foldr' cons zero . getLVectorInner+    map_v f = newLVector . map f . getLVectorInner+    zipWith_v f v1 v2 = newLVector $ zipWith f (getLVectorInner v1) (getLVectorInner v2)++instance Matrix LLMatrix e where+    generate_m rows cols gen = newLLMatrix [[gen r c | r <- [0..(rows - 1)]] | c <- [0..(cols - 1)]]+    --fromList_m m n l = se puede mejorar+    transpose_m = newLLMatrix . transpose . getLLMatrixInner+    dim_m m = (length $ head (getLLMatrixInner m), length $ getLLMatrixInner m)+    elem_m i j m = getLLMatrixInner m !! j !! i+    map_m f m = newLLMatrix . map (map f) $ getLLMatrixInner m+    zipWith_m f m1 m2 = newLLMatrix $ zipWith (zipWith f) (getLLMatrixInner m1) (getLLMatrixInner m2)+    subMatrix_m posI posJ cantRows cantCols = newLLMatrix . map (take cantRows . drop posI) . take cantCols . drop posJ . getLLMatrixInner+    fromBlocks_m = expandVert . map expandHoriz+        where+            expandVert lm = newLLMatrix [concatMap (getLVectorInner . col_vm c) lm | c <- [0 .. (cantCols_m (head lm) - 1)]]+            expandHoriz = newLLMatrix . concatMap getLLMatrixInner+    --toBlocks_m m n m1 =++instance MatrixVector LLMatrix LVector e where+    row_vm pos m = newLVector $ map (!! pos) (getLLMatrixInner m)+    col_vm pos m = newLVector $ getLLMatrixInner m !! pos+    fromCols_vm = newLLMatrix . map getLVectorInner+    toCols_vm m = map newLVector $ getLLMatrixInner m
+ src/FPNLA/Matrix/Instances/LLMatrixByRows.hs view
@@ -0,0 +1,74 @@+{-# LANGUAGE FlexibleInstances     #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE OverlappingInstances  #-}+{-# LANGUAGE TypeFamilies          #-}++module FPNLA.Matrix.Instances.LLMatrixByRows (+    LVector(LVector),+    LLMatrix(LLMatrix),+) where++import           Control.DeepSeq             (NFData (rnf))+import           Control.Parallel.Strategies (rdeepseq, withStrategy)+import           Data.Foldable               (foldr')+import           Data.List                   (transpose)+import           FPNLA.Matrix                (Matrix (generate_m, fromList_m, transpose_m, dim_m, elem_m, map_m, zipWith_m, subMatrix_m, fromBlocks_m), MatrixVector (fromCols_vm, toCols_vm, row_vm, col_vm), Vector (generate_v, fromList_v, concat_v, elem_v, length_v, foldr_v, map_v, zipWith_v),+                                              cantCols_m, cantRows_m)+import           FPNLA.Utils                 (splitSlice)++data LVector e = LVector [e]+data LLMatrix e = LLMatrix [[e]]++getLVectorInner :: LVector t -> [t]+getLVectorInner (LVector l) = l+getLLMatrixInner :: LLMatrix t -> [[t]]+getLLMatrixInner (LLMatrix ll) = ll+newLVector :: [e] -> LVector e+newLVector = LVector+newLLMatrix :: [[e]] -> LLMatrix e+newLLMatrix = LLMatrix++instance Show e => Show (LLMatrix e) where+    show = show . getLLMatrixInner++instance (NFData e) => NFData (LVector e) where+    rnf v =+        withStrategy rdeepseq (getLVectorInner v) `seq` ()++instance (NFData e) => NFData (LLMatrix e) where+    rnf m =+        withStrategy rdeepseq (getLLMatrixInner m) `seq` ()++instance Vector LVector e where+    generate_v size gen = newLVector [gen x | x <- [0..(size - 1)]]+    fromList_v = newLVector+    concat_v = newLVector . concatMap getLVectorInner+    elem_v pos v = getLVectorInner v !! pos+    length_v = length . getLVectorInner+    foldr_v cons zero = foldr' cons zero . getLVectorInner+    map_v f = newLVector . map f . getLVectorInner+    zipWith_v f v1 v2 = newLVector $ zipWith f (getLVectorInner v1) (getLVectorInner v2)++instance Matrix LLMatrix e where+    generate_m rows cols gen = newLLMatrix [[gen r c | c <- [0..(cols - 1)]] | r <- [0..(rows - 1)]]+    fromList_m m n l = newLLMatrix . splitSlice n $ take (m*n) l+    transpose_m = newLLMatrix . transpose . getLLMatrixInner+    dim_m m = (length $ getLLMatrixInner m, length $ head (getLLMatrixInner m))+    elem_m i j m = getLLMatrixInner m !! i !! j+    map_m f m = newLLMatrix . map (map f) $ getLLMatrixInner m+    zipWith_m f m1 m2 = newLLMatrix $ zipWith (zipWith f) (getLLMatrixInner m1) (getLLMatrixInner m2)+    subMatrix_m posI posJ cantRows cantCols = newLLMatrix . map (take cantCols . drop posJ) . take cantRows . drop posI . getLLMatrixInner+    fromBlocks_m = expandVert . map expandHoriz+        where+            expandHoriz lm = newLLMatrix [concatMap (getLVectorInner . row_vm r) lm | r <- [0 .. (cantRows_m (head lm) - 1)]]+            expandVert = newLLMatrix . concatMap getLLMatrixInner+    --toBlocks_m m n m1 =++instance MatrixVector LLMatrix LVector e where+    row_vm pos m = newLVector $ getLLMatrixInner m !! pos+    col_vm pos m = newLVector $ map (!! pos) (getLLMatrixInner m)+    fromCols_vm lv = newLLMatrix [[getLVectorInner (lv !! c) !! r | c <- [0 .. (cols - 1)]]| r <- [0 .. (rows - 1)]]+        where+            rows = length_v (head lv)+            cols = length lv+    toCols_vm m = map (`col_vm` m) [0 .. (cantCols_m m - 1)]
+ src/FPNLA/Matrix/Instances/RepaMatrix.hs view
@@ -0,0 +1,62 @@+{-# LANGUAGE FlexibleInstances     #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE OverlappingInstances  #-}++module FPNLA.Matrix.Instances.RepaMatrix (RepaMatrix, RepaVector) where++import           Control.DeepSeq              (NFData (rnf))+import           FPNLA.Matrix                 (Matrix (generate_m, fromList_m, transpose_m, dim_m, elem_m, map_m, zipWith_m, subMatrix_m, fromBlocks_m), MatrixVector (row_vm, col_vm, fromCols_vm), Vector (generate_v, fromList_v, concat_v, elem_v, length_v, foldr_v, map_v, zipWith_v))++import           Data.Array.Repa              ((:.) (..), All (..), Any (..),+                                               Array, D, DIM1, DIM2, Shape,+                                               Source, Z (..), append,+                                               computeUnboxedS, deepSeqArray,+                                               delay, extend, extent, extract,+                                               fromFunction, fromFunction,+                                               fromListUnboxed, index, map,+                                               size, slice, toList, transpose,+                                               zipWith)+import           Data.Array.Repa.Eval         ()+import           Data.Array.Repa.Repr.Unboxed (Unbox ())+import           Data.Foldable                (foldr')+import           Prelude                      (Int, Show (..), foldr1, length,+                                               ($), (.))+import qualified Prelude                      as P hiding (Show)++type RepaVector = Array D DIM1+type RepaMatrix = Array D DIM2++instance (Shape sh, Show sh, Show e, Unbox e) => Show (Array D sh e) where+    show = show . computeUnboxedS++instance (Source r e, Shape sh, NFData e) => NFData (Array r sh e) where+    rnf m = deepSeqArray m ()++instance (Unbox e) => Vector RepaVector e where+    generate_v size gen = fromFunction (Z:.size) (\(Z:.pos) -> gen pos)+    fromList_v l = delay $ fromListUnboxed (Z :. length l) l+    concat_v = foldr1 append+    elem_v pos v = index v (Z:.pos)+    length_v v = size $ extent v+    --foldr_v cons zero v = foldAllS cons zero v -- foldAllS no sirve por su tipo+    foldr_v cons zero v = foldr' cons zero (toList v) -- es lo que hay...+    map_v = map+    zipWith_v = zipWith++instance (Unbox e) => Matrix RepaMatrix e where+    generate_m rows cols gen = fromFunction (Z:.rows:.cols) (\(Z:.i:.j) -> gen i j)+    fromList_m m n l = delay $ fromListUnboxed (Z:.m:.n) l+    transpose_m = transpose+    dim_m = (\(Z:.rows:.cols) -> (rows, cols)) . extent+    elem_m i j m = index m (Z:.i:.j)+    map_m = map+    zipWith_m = zipWith+    subMatrix_m posI posJ cantRows cantCols = extract (Z:.posI:.posJ) (Z:.cantRows:.cantCols)+    fromBlocks_m = transpose . foldr1 append . P.map (transpose . foldr1 append)+    --toBlocks_m = -- Queda el default usado submatrices (extract)++instance (Unbox e) => MatrixVector RepaMatrix RepaVector e  where+    row_vm pos m = slice m (Any:.pos:.All)+    col_vm pos m = slice m (Any:.pos)+    fromCols_vm = foldr1 append . P.map (transpose . extend (Any:.(1::Int):.All))+    --toCols_vm = -- Queda el default usado submatrices (extract)
+ src/FPNLA/Matrix/Instances/VectorMatrix.hs view
@@ -0,0 +1,72 @@+{-# LANGUAGE FlexibleInstances     #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE OverlappingInstances  #-}++module FPNLA.Matrix.Instances.VectorMatrix (+    VVector,+    VMatrix(VMatrix),+) where++import           Control.DeepSeq             (NFData, rnf)+import           Control.Parallel.Strategies (rdeepseq, withStrategy)+import qualified Data.Foldable               (foldr)+import qualified Data.Vector                 as DV (Vector, concat, drop,+                                                    fromList, generate, head,+                                                    length, map, take, toList,+                                                    zipWith, (!), (++))+import           FPNLA.Matrix                (Matrix (..), MatrixVector (..),+                                              Vector (..), cantCols_m)+import           FPNLA.Utils                 (iif)+{-+instance (NFData a) => NFData (DV.Vector a) where+  rnf v = evalVect (DV.length v)+    where evalVect i+            | i > 0 = withStrategy rdeepseq (v DV.! (i-1)) `seq` evalVect (i-1)+            | otherwise = ()+-}+instance (NFData a) => NFData (VMatrix a) where+    rnf (VMatrix v) = withStrategy rdeepseq v `seq` ()++--+type VVector = DV.Vector+instance Vector (DV.Vector) e where+    generate_v = DV.generate+    fromList_v = DV.fromList+    concat_v = DV.concat+    elem_v = flip (DV.!)+    length_v = DV.length+    foldr_v = Data.Foldable.foldr+    map_v = DV.map+    zipWith_v = DV.zipWith++--+newtype VMatrix e = VMatrix (DV.Vector (DV.Vector e)) deriving (Show)+wrapIn :: DV.Vector (DV.Vector e) -> VMatrix e+wrapIn = VMatrix+wrapOut :: VMatrix e -> DV.Vector (DV.Vector e)+wrapOut (VMatrix vv) = vv++instance Matrix VMatrix e where+    generate_m cr cc gen = wrapIn $ DV.generate cc (\j -> DV.generate cr (`gen` j))+    --fromList_m =+    --transpose_m =+    -- Se guarda por columnas!+    dim_m m = (iif (cc > 0) cr 0, cc)+        where cc = DV.length (wrapOut m)+              cr = DV.length $ DV.head (wrapOut m)+    elem_m i j m = wrapOut m DV.! j DV.! i+    map_m f m = wrapIn $ DV.map (DV.map f) (wrapOut m)+    zipWith_m f m1 m2 = wrapIn $ DV.zipWith (DV.zipWith f) (wrapOut m1) (wrapOut m2)+    subMatrix_m posI posJ cantRows cantCols =+        wrapIn . DV.map (DV.take cantRows . DV.drop posI) . DV.take cantCols . DV.drop posJ . wrapOut+    fromBlocks_m = expandVert . map expandHoriz+        where+            expandVert lm = wrapIn . DV.generate (cantCols_m (head lm)) $ (\c -> foldr1 (DV.++) $ map (col_vm c) lm)+            expandHoriz = wrapIn . foldr1 (DV.++) . map wrapOut+    --toBlocks_m++instance MatrixVector VMatrix DV.Vector e where+    row_vm i m = DV.map (DV.! i) (wrapOut m)+    col_vm j m = wrapOut m DV.! j+    fromCols_vm = wrapIn . DV.fromList+    toCols_vm = DV.toList . wrapOut
+ src/FPNLA/Operations/BLAS/Strategies.hs view
@@ -0,0 +1,39 @@+module FPNLA.Operations.BLAS.Strategies (+    module FPNLA.Operations.BLAS.Strategies.DataTypes+) where++import FPNLA.Operations.BLAS.Strategies.DataTypes++-- Solo importan las instancias+import FPNLA.Operations.BLAS.Strategies.GEMM.DefSeq ()+import FPNLA.Operations.BLAS.Strategies.GEMM.CBindSeq ()+import FPNLA.Operations.BLAS.Strategies.GEMM.BlocksSeq ()++import FPNLA.Operations.BLAS.Strategies.GEMM.MonadPar.DefPar ()+import FPNLA.Operations.BLAS.Strategies.GEMM.MonadPar.StrassenPar ()+import FPNLA.Operations.BLAS.Strategies.GEMM.MonadPar.StrassenMPar ()+import FPNLA.Operations.BLAS.Strategies.GEMM.MonadPar.BlocksPar ()++import FPNLA.Operations.BLAS.Strategies.GEMM.Strategies.DefPar ()+import FPNLA.Operations.BLAS.Strategies.GEMM.Strategies.StrassenPar ()+import FPNLA.Operations.BLAS.Strategies.GEMM.Strategies.BlocksPar ()++import FPNLA.Operations.BLAS.Strategies.GEMM.Accelerate.DefPar ()++import FPNLA.Operations.BLAS.Strategies.GEMV ()++import FPNLA.Operations.BLAS.Strategies.DOT.DefSeq ()++import FPNLA.Operations.BLAS.Strategies.SYRK.DefSeq ()+import FPNLA.Operations.BLAS.Strategies.SYRK.Strategies.DefPar ()+import FPNLA.Operations.BLAS.Strategies.SYRK.MonadPar.DefPar ()+++import FPNLA.Operations.BLAS.Strategies.TRSM.DefSeq ()+import FPNLA.Operations.BLAS.Strategies.TRSM.BlocksByCols ()+import FPNLA.Operations.BLAS.Strategies.TRSM.BlocksByRows ()+import FPNLA.Operations.BLAS.Strategies.TRSM.CBindSeq ()++import FPNLA.Operations.BLAS.Strategies.TRSM.MonadPar.ColsPar ()++import FPNLA.Operations.BLAS.Strategies.TRSM.Strategies.ColsPar ()
+ src/FPNLA/Operations/BLAS/Strategies/DOT/DefSeq.hs view
@@ -0,0 +1,14 @@+{-# LANGUAGE FlexibleInstances     #-}+{-# LANGUAGE MultiParamTypeClasses #-}++module FPNLA.Operations.BLAS.Strategies.DOT.DefSeq (++) where++import           FPNLA.Matrix                (Vector, foldr_v, zipWith_v)+import           FPNLA.Operations.BLAS       (DOT (dot))+import           FPNLA.Operations.Parameters (Elt, blasResultS)+import           FPNLA.Operations.BLAS.Strategies.DataTypes (DefSeq)++instance (Elt e, Vector v e) => DOT DefSeq v e where+    dot _ pvA pvB = blasResultS . foldr_v (+) 0 $ zipWith_v (*) pvA pvB
+ src/FPNLA/Operations/BLAS/Strategies/DataTypes.hs view
@@ -0,0 +1,102 @@+{-# LANGUAGE EmptyDataDecls        #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE TypeFamilies          #-}++module FPNLA.Operations.BLAS.Strategies.DataTypes (+    DefSeq,+    CBindSeq,+    BlocksByCols,+    BlocksByRows,+    BlocksSeq,++    DefPar_ST,+    ColsPar_ST,+    StrassenPar_ST,+    BlocksPar_ST,++    DefPar_MP,+    ColsPar_MP,+    StrassenPar_MP,+    StrassenMPar_MP,+    BlocksPar_MP,++    DefPar_Acc,++    NullContext,+    BlockContext(getBlockDim),+    SqrBlockContext(getSqrBlockDim),+    StrassenContext(getStrassenLimit),++    newNullContext,+    newBlockContext,+    newSqrBlockContext,+    newStrassenContext,+) where++import           FPNLA.Operations.Parameters (StratCtx)++-- Contextos:+data NullContext = NullContext deriving (Show)+data BlockContext = BlockContext {+        getBlockDim :: (Int, Int)+    } deriving (Show)+data SqrBlockContext = SqrBlockContext {+        getSqrBlockDim :: Int+    } deriving (Show)+data StrassenContext = StrassenContext {+        getStrassenExponentialLimit :: Int,+        getStrassenLimit            :: Int+    } deriving (Show)++newNullContext :: NullContext+newNullContext = NullContext++newBlockContext :: Int -> Int -> BlockContext+newBlockContext m n = BlockContext { getBlockDim = (m,n) }++newSqrBlockContext :: Int -> SqrBlockContext+newSqrBlockContext d = SqrBlockContext { getSqrBlockDim = d }++newStrassenContext :: Int -> StrassenContext+newStrassenContext e = StrassenContext { getStrassenExponentialLimit = e, getStrassenLimit = 2 ^ e }++-- Estrategias:+data DefSeq+data CBindSeq+data BlocksSeq gs+data BlocksByCols gs ts+data BlocksByRows gs ts++-- Strategies+data DefPar_ST+data ColsPar_ST s+data StrassenPar_ST gs+data BlocksPar_ST gs++-- Monad-Par+data DefPar_MP+data ColsPar_MP s+data StrassenPar_MP gs+data StrassenMPar_MP gs+data BlocksPar_MP gs++data DefPar_Acc++type instance StratCtx DefSeq = NullContext+type instance StratCtx CBindSeq = NullContext+type instance StratCtx (BlocksByCols gs ts) = (SqrBlockContext, StratCtx gs, StratCtx ts)+type instance StratCtx (BlocksByRows gs ts) = (SqrBlockContext, StratCtx gs, StratCtx ts)+type instance StratCtx (BlocksSeq gs) = (BlockContext, StratCtx gs)++type instance StratCtx DefPar_ST = NullContext+type instance StratCtx (ColsPar_ST s) = StratCtx s+type instance StratCtx (StrassenPar_ST gs)= (StrassenContext, StratCtx gs)+type instance StratCtx (BlocksPar_ST gs) = (BlockContext, StratCtx gs)++type instance StratCtx DefPar_MP = NullContext+type instance StratCtx (ColsPar_MP s) = StratCtx s+type instance StratCtx (StrassenPar_MP gs)= (StrassenContext, StratCtx gs)+type instance StratCtx (StrassenMPar_MP gs)= (StrassenContext, StratCtx gs)+type instance StratCtx (BlocksPar_MP gs) = (BlockContext, StratCtx gs)++type instance StratCtx DefPar_Acc = NullContext
+ src/FPNLA/Operations/BLAS/Strategies/GEMM/Accelerate/DefPar.hs view
@@ -0,0 +1,64 @@+{-# LANGUAGE FlexibleInstances     #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE ScopedTypeVariables   #-}+{-# LANGUAGE TypeOperators         #-}++module FPNLA.Operations.BLAS.Strategies.GEMM.Accelerate.DefPar () where++import           FPNLA.Matrix.Instances.AccMatrix           (AccMatrix (..),+                                                             AccVector (..))+import           FPNLA.Operations.BLAS                      (GEMM (gemm))+import           FPNLA.Operations.BLAS.Strategies.DataTypes (DefPar_Acc)+import           FPNLA.Operations.Parameters                (Elt,+                                                             TransType (..),+                                                             blasResultM)++import           Data.Array.Accelerate                      ((:.) (..), Acc,+                                                             All (All), Array,+                                                             DIM2, Exp, IsNum,+                                                             Z (Z),+                                                             constant,+                                                             shape)+import qualified Data.Array.Accelerate                      as A (Elt, fold,+                                                                  lift, map,+                                                                  replicate,+                                                                  transpose,+                                                                  unlift,+                                                                  zipWith)++instance (Elt e, A.Elt e, IsNum e) => GEMM DefPar_Acc AccMatrix AccVector e where+    gemm _ pmA pmB alpha beta (AccMatrix mC) =+        blasResultM (AccMatrix mC')++        where+            alpha_expr = constant alpha+            beta_expr = constant beta+            mC' = A.map (beta_expr*) . A.zipWith (+) mC $ matMul mA mB+            mA = A.map (alpha_expr*) $ unAccTrans pmA+            mB = unAccTrans pmB+            unAccTrans pm = case pm of (NoTrans (AccMatrix m)) -> m+                                       (Trans (AccMatrix m)) -> A.transpose m+                                       (ConjTrans (AccMatrix m)) -> A.transpose m --no hay complejos en A.Elt+++matMul :: (A.Elt e, IsNum e) => Acc (Array DIM2 e) -> Acc (Array DIM2 e) -> Acc (Array DIM2 e)+matMul arr brr+  = A.fold (+) 0 $ A.zipWith (*) arrRepl brrRepl+  where+    Z :. rowsA :. _     = A.unlift (shape arr) :: Z :. Exp Int :. Exp Int+    Z :. _     :. colsB = A.unlift (shape brr) :: Z :. Exp Int :. Exp Int++    arrRepl             = A.replicate (A.lift $ Z :. All   :. colsB :. All) arr+    brrRepl             = A.replicate (A.lift $ Z :. rowsA :. All   :. All) (A.transpose brr)+{-+matMul2 :: (A.Elt e, IsNum e) => Acc (Array DIM2 e) -> Acc (Array DIM2 e) -> Acc (Array DIM2 e)+matMul2 arr brr =+    let+        (Z:.rowsA:.colsA) = A.unlift (shape arr) :: Z :. Exp Int :. Exp Int+        (Z:.rowsB:.colsB) = A.unlift (shape brr) :: Z :. Exp Int :. Exp Int+        -- Transpongo mB:+        mB'_acc = backpermute (A.lift $ Z:.colsB:.rowsB) (\e -> A.uncurry index2 $ A.lift (A.snd $ unindex2 e, A.fst $ unindex2 e)) brr+        repB = A.replicate (A.lift $ Z:.rowsA:.All:.All) mB'_acc+        repA = A.replicate (A.lift $ Z:.All:.colsB:.All) arr+    in A.fold1 (+) $ A.zipWith (*) repA repB+-}
+ src/FPNLA/Operations/BLAS/Strategies/GEMM/BlocksSeq.hs view
@@ -0,0 +1,33 @@+{-# LANGUAGE FlexibleContexts      #-}+{-# LANGUAGE FlexibleInstances     #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE ScopedTypeVariables   #-}++module FPNLA.Operations.BLAS.Strategies.GEMM.BlocksSeq () where++import           FPNLA.Matrix                               (fromBlocks_m,+                                                             map_m, toBlocks_m,+                                                             zipWith_m)+import           FPNLA.Operations.BLAS                      (GEMM (gemm))+import           FPNLA.Operations.BLAS.Strategies.DataTypes (BlockContext (getBlockDim),+                                                             BlocksSeq)+import           FPNLA.Operations.Parameters                (ResM, TransType (NoTrans),+                                                             blasResultM,+                                                             getResultDataM,+                                                             transTrans_m)++instance (GEMM s m v e) => GEMM (BlocksSeq s) m v e where+    gemm (bctx, sctx) tmA tmB alpha beta mC = blasResultM . fromBlocks_m $ generate rA cB (\i j -> add_m (matMultIJ i j (bmC !! i !! j)) (map_m (beta *) (bmC !! i !! j)))+        where+            (r, c) = getBlockDim bctx+            add_m = zipWith_m (+)+            generate = generateBlocks+            generateBlocks m n f = [[f i j | j <- [0 .. (n - 1)]] | i <- [0 .. (m - 1)]]+            matMultIJ i j bmCIJ = foldr1 add_m  [ callGemm (bmA !! i !! k) (bmB !! k !! j) bmCIJ |  k <- [0 .. (cA - 1)]]+            callGemm mA mB mC = getResultDataM (gemm sctx (NoTrans mA) (NoTrans mB) alpha 0 mC :: ResM s v m e)+            bmA = toBlocks_m r c $ transTrans_m tmA+            bmB = toBlocks_m r c $ transTrans_m tmB+            bmC = toBlocks_m r c mC+            rA = length bmA+            cA = length (head bmA)+            cB = length (head bmB)
+ src/FPNLA/Operations/BLAS/Strategies/GEMM/CBindSeq.hs view
@@ -0,0 +1,55 @@+{-# LANGUAGE FlexibleContexts      #-}+{-# LANGUAGE FlexibleInstances     #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE ScopedTypeVariables   #-}++module FPNLA.Operations.BLAS.Strategies.GEMM.CBindSeq () where++import           FPNLA.Matrix                               (Matrix,+                                                             MatrixVector,+                                                             cantRows_m)+import           FPNLA.Operations.BLAS                      (GEMM (gemm))+import           FPNLA.Operations.BLAS.Strategies.DataTypes (CBindSeq)+import           FPNLA.Operations.Parameters                (Elt, blasResultM,+                                                             dimTrans_m,+                                                             unTransT)+import           FPNLA.Operations.Utils                     (matrixToPtr,+                                                             ptrToMatrix,+                                                             transToForeign,+                                                             unsafePerformIO)++import           Control.DeepSeq                            (NFData, rnf)+import           Control.Exception                          (evaluate)+import           Foreign.BLAS                               as L3 (BLAS3, Trans,+                                                                   gemm)+import           Foreign.Marshal.Alloc                      (free)++instance  (BLAS3 e, Elt e, NFData (m e), MatrixVector m v e) => GEMM CBindSeq m v e where+    gemm _ pmA pmB alfa beta mC =+        let ttA = transToForeign pmA+            ttB = transToForeign pmB+            mA = snd $ unTransT pmA+            mB = snd $ unTransT pmB+            (ldA, ldB, ldC) = (cantRows_m mA, cantRows_m mB, cantRows_m mC)+            (m, k, n) = (fst $ dimTrans_m pmA, snd $ dimTrans_m pmA, snd $ dimTrans_m pmB)+        in+            blasResultM $ ioGemm ttA ttB m n k alfa mA mB beta mC ldA ldB ldC++ioGemm :: (Matrix m a, Matrix m1 a,+             Matrix m2 a, Matrix m3 a, Elt a, BLAS3 a,+             NFData (m3 a)) =>+            Trans -> Trans -> Int -> Int -> Int+            -> a -> m a -> m1 a -> a -> m2 a +            -> Int -> Int -> Int -> m3 a+ioGemm ttA ttB m n k alfa mA mB beta mC ldA ldB ldC= unsafePerformIO $ do+        mAptr <- matrixToPtr mA+        mBptr <- matrixToPtr mB+        mCptr <- matrixToPtr mC+        --L3.gemm  :: Trans -> Trans -> Int -> Int -> Int -> a -> Ptr a -> Int -> Ptr a -> Int -> a -> Ptr a -> Int -> IO ()+        L3.gemm ttA ttB m n k alfa mAptr ldA mBptr ldB beta mCptr ldC+        res <- ptrToMatrix m n mCptr+        evaluate $ rnf res+        free mAptr+        free mBptr+        free mCptr+        return res
+ src/FPNLA/Operations/BLAS/Strategies/GEMM/DefSeq.hs view
@@ -0,0 +1,24 @@+{-# LANGUAGE FlexibleInstances     #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE ScopedTypeVariables   #-}++module FPNLA.Operations.BLAS.Strategies.GEMM.DefSeq () where++import           FPNLA.Matrix                               (MatrixVector,+                                                             elem_m, foldr_v,+                                                             generate_m,+                                                             generate_v)+import           FPNLA.Operations.BLAS                      (GEMM (gemm))+import           FPNLA.Operations.BLAS.Strategies.DataTypes (DefSeq)+import           FPNLA.Operations.Parameters                (Elt, blasResultM,+                                                             dimTrans_m,+                                                             elemTrans_m)++instance  (Elt e, MatrixVector m v e) => GEMM DefSeq m v e where+    gemm _ pmA pmB alpha beta mC+        | p /= p' = error "gemm: incompatible ranges"+        | otherwise = blasResultM . generate_m m n $ (\i j -> (alpha * matMultIJ i j) + beta * elem_m i j mC)+        where+            (m, p) = dimTrans_m pmA+            (p',n) = dimTrans_m pmB+            matMultIJ i j = foldr_v (+) 0 (generate_v p (\k -> (*) (elemTrans_m i k pmA) (elemTrans_m k j pmB)) :: v e)
+ src/FPNLA/Operations/BLAS/Strategies/GEMM/MonadPar/BlocksPar.hs view
@@ -0,0 +1,37 @@+{-# LANGUAGE FlexibleContexts      #-}+{-# LANGUAGE FlexibleInstances     #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE ScopedTypeVariables   #-}++module FPNLA.Operations.BLAS.Strategies.GEMM.MonadPar.BlocksPar () where++import           FPNLA.Matrix                               (fromBlocks_m,+                                                             map_m, toBlocks_m,+                                                             zipWith_m)+import           FPNLA.Operations.BLAS                      (GEMM (gemm))+import           FPNLA.Operations.BLAS.Strategies.DataTypes (BlockContext (getBlockDim),+                                                             BlocksPar_MP)+import           FPNLA.Operations.Parameters                (ResM, TransType (NoTrans),+                                                             blasResultM,+                                                             getResultDataM,+                                                             transTrans_m)++import           Control.DeepSeq                            (NFData)+import           Control.Monad.Par                          as MP (parMap,+                                                                   runPar)++instance (NFData (m e), GEMM s m v e) => GEMM (BlocksPar_MP s) m v e where+    gemm (bctx, sctx) tmA tmB alpha beta mC = blasResultM . fromBlocks_m $ generatePar rA cB (\i j -> add_m (matMultIJ i j (bmC !! i !! j)) (map_m (beta *) (bmC !! i !! j)))+        where+            (r, c) = getBlockDim bctx+            add_m = zipWith_m (+)+            generatePar m n f = runPar . parMap id $ generateBlocks m n f+            generateBlocks m n f = [[f i j | j <- [0 .. (n - 1)]] | i <- [0 .. (m - 1)]]+            matMultIJ i j bmCIJ = foldr1 add_m  [callGemm (bmA !! i !! k) (bmB !! k !! j) bmCIJ |  k <- [0 .. (cA - 1)]]+            callGemm mA mB mC = getResultDataM (gemm sctx (NoTrans mA) (NoTrans mB) alpha 0 mC :: ResM s v m e)+            bmA = toBlocks_m r c $ transTrans_m tmA+            bmB = toBlocks_m r c $ transTrans_m tmB+            bmC = toBlocks_m r c mC+            rA = length bmA+            cA = length (head bmA)+            cB = length (head bmB)
+ src/FPNLA/Operations/BLAS/Strategies/GEMM/MonadPar/DefPar.hs view
@@ -0,0 +1,29 @@+{-# LANGUAGE FlexibleContexts      #-}+{-# LANGUAGE FlexibleInstances     #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE ScopedTypeVariables   #-}++module FPNLA.Operations.BLAS.Strategies.GEMM.MonadPar.DefPar () where++import           Control.DeepSeq                            (NFData)+import           Control.Monad.Par                          as MP (parMap,+                                                                   runPar)+import           FPNLA.Matrix                               (MatrixVector (fromCols_vm),+                                                             elem_m, foldr_v,+                                                             generate_v)+import           FPNLA.Operations.BLAS                      (GEMM (gemm))+import           FPNLA.Operations.BLAS.Strategies.DataTypes (DefPar_MP)+import           FPNLA.Operations.Parameters                (Elt, blasResultM,+                                                             dimTrans_m,+                                                             elemTrans_m)++instance (NFData (v e), Elt e, MatrixVector m v e) => GEMM DefPar_MP m v e where+    gemm _ pmA pmB alpha beta mC+        | p /= p' = error "gemm: incompatible ranges"+        | otherwise = blasResultM $ generatePar_m m n (\i j -> (alpha * matMultIJ i j) + beta * elem_m i j mC)+        where+            (m, p) = dimTrans_m pmA+            (p',n) = dimTrans_m pmB+            matMultIJ i j = foldr_v (+) 0 (generate_v p (\k -> (*) (elemTrans_m i k pmA) (elemTrans_m k j pmB)) :: v e)+            generatePar_m m n gen = fromCols_vm . MP.runPar . MP.parMap (\j -> generate_v m (`gen` j) :: v e) $ [0 .. (n - 1)]+
+ src/FPNLA/Operations/BLAS/Strategies/GEMM/MonadPar/StrassenMPar.hs view
@@ -0,0 +1,80 @@+{-# LANGUAGE FlexibleContexts      #-}+{-# LANGUAGE FlexibleInstances     #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE ScopedTypeVariables   #-}++module FPNLA.Operations.BLAS.Strategies.GEMM.MonadPar.StrassenMPar () where++import           Control.DeepSeq                            (NFData)+import           Control.Monad.Par                          as MP (Par, parMap,+                                                                   runPar)+import           FPNLA.Matrix                               (cantCols_m,+                                                             cantRows_m,+                                                             fromBlocks_m,+                                                             generate_m,+                                                             toBlocks_m,+                                                             zipWith_m)+import           FPNLA.Operations.BLAS                      (GEMM (gemm))+import           FPNLA.Operations.BLAS.Strategies.DataTypes (StrassenMPar_MP,+                                                             getStrassenLimit)+import           FPNLA.Operations.Parameters                (Elt, ResM, TransType (NoTrans),+                                                             blasResultM,+                                                             getResultDataM,+                                                             transTrans_m)++instance (NFData (m e), Elt e, GEMM s m v e) => GEMM (StrassenMPar_MP s) m v e where+    gemm (ctx, sctx) pmA pmB alpha beta mC =+        blasResultM . MP.runPar $ strassen mA mB alpha beta mC+        where++            callGemm mA mB alpha beta mC = getResultDataM (gemm sctx (NoTrans mA) (NoTrans mB) alpha beta mC :: ResM s v m e)++            mA = transTrans_m pmA+            mB = transTrans_m pmB++            strassenLimit = getStrassenLimit ctx++            strassen :: m e -> m e -> e -> e -> m e -> Par (m e)+            strassen mA mB alpha beta mC =+                if cantRows_m mA <= strassenLimit+                then return $ callGemm mA mB alpha beta mC+                else do+                    m <- toM (quadPartition mA, quadPartition mB, quadPartition mC)+                    c <- toC m+                    return $ quadJoin c++            quadPartition :: m e -> (m e, m e, m e, m e)+            quadPartition m =+                let+                blocks = toBlocks_m (div (cantRows_m m) 2) (div (cantCols_m m) 2) m in+                (head (head blocks), head blocks !! 1, head (blocks !! 1), blocks !! 1 !! 1)++            fillWithZeros :: m e -> m e+            fillWithZeros m = generate_m (cantRows_m m) (cantCols_m m) (\_ _ -> 0)++            toM :: ((m e, m e, m e, m e), (m e, m e, m e, m e), (m e, m e, m e, m e)) -> Par (m e, m e, m e, m e, m e, m e, m e)+            toM ((a11, a12, a21, a22), (b11, b12, b21, b22), (c11, c12, c21, c22)) = do+                x <- sequence [strassen (zipWith_m (+) a11 a22) (zipWith_m (+) b11 b22) alpha 0 (fillWithZeros c11),+                    strassen (zipWith_m (+) a21 a22) b11 alpha 0 (fillWithZeros c11),+                    strassen a11 (zipWith_m (-) b12 b22) alpha beta c12,+                    strassen a22 (zipWith_m (-) b21 b11) alpha beta c21,+                    strassen (zipWith_m (+) a11 a12) b22 alpha 0 (fillWithZeros c11),+                    strassen (zipWith_m (-) a21 a11) (zipWith_m (+) b11 b12) alpha beta (zipWith_m (-) c22 c12),+                    strassen (zipWith_m (-) a12 a22) (zipWith_m (+) b21 b22) alpha beta (zipWith_m (-) c11 c21)]+                y <- MP.parMap id x+                return $ (\[a, b, c, d, e, f, g] -> (a, b, c, d, e, f, g)) y++            toC :: (m e, m e, m e, m e, m e, m e, m e) -> Par (m e, m e, m e, m e)+            toC (m1, m2, m3, m4, m5, m6, m7) = do+                x <- MP.parMap id+                    [zipWith_m (+) m1 . zipWith_m (+) m4 . zipWith_m (-) m7 $ m5,+                    zipWith_m (+) m3 m5,+                    zipWith_m (+) m2 m4,+                    zipWith_m (+) m1 . zipWith_m (+) m3 . zipWith_m (-) m6 $ m2]+                return $ (\[a, b, c, d] -> (a, b, c, d)) x++            quadJoin :: (m e, m e, m e, m e) -> m e+            quadJoin (c11, c12, c21, c22) =+                fromBlocks_m [+                    [c11, c12],+                    [c21, c22]]
+ src/FPNLA/Operations/BLAS/Strategies/GEMM/MonadPar/StrassenPar.hs view
@@ -0,0 +1,76 @@+{-# LANGUAGE FlexibleContexts      #-}+{-# LANGUAGE FlexibleInstances     #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE ScopedTypeVariables   #-}++module FPNLA.Operations.BLAS.Strategies.GEMM.MonadPar.StrassenPar () where++import           Control.DeepSeq                            (NFData)+import           Control.Monad.Par                          as MP (parMap,+                                                                   runPar)+import           FPNLA.Matrix                               (cantCols_m,+                                                             cantRows_m,+                                                             fromBlocks_m,+                                                             fromBlocks_m,+                                                             generate_m,+                                                             toBlocks_m,+                                                             zipWith_m)+import           FPNLA.Operations.BLAS                      (GEMM (gemm))+import           FPNLA.Operations.BLAS.Strategies.DataTypes (StrassenPar_MP,+                                                             getStrassenLimit)+import           FPNLA.Operations.Parameters                (Elt, ResM, TransType (NoTrans),+                                                             blasResultM,+                                                             getResultDataM,+                                                             transTrans_m)++instance (NFData (m e), Elt e, GEMM s m v e) => GEMM (StrassenPar_MP s) m v e where+    gemm (ctx, sctx) pmA pmB alpha beta mC =+        blasResultM $ strassen mA mB alpha beta mC+        where++            callGemm mA mB alpha beta mC = getResultDataM (gemm sctx (NoTrans mA) (NoTrans mB) alpha beta mC :: ResM s v m e)++            mA = transTrans_m pmA+            mB = transTrans_m pmB++            strassenLimit = getStrassenLimit ctx++            strassen :: m e -> m e -> e -> e -> m e -> m e+            strassen mA mB alpha beta mC =+                if cantRows_m mA <= strassenLimit+                then callGemm mA mB alpha beta mC+                else quadJoin . toC . toM $ (quadPartition mA, quadPartition mB, quadPartition mC)++            quadPartition :: m e -> (m e, m e, m e, m e)+            quadPartition m =+                let+                blocks = toBlocks_m (div (cantRows_m m) 2) (div (cantCols_m m) 2) m in+                (head (head blocks), head blocks !! 1, head (blocks !! 1), blocks !! 1 !! 1)++            fillWithZeros :: m e -> m e+            fillWithZeros m = generate_m (cantRows_m m) (cantCols_m m) (\_ _ -> 0)++            toM :: ((m e, m e, m e, m e), (m e, m e, m e, m e), (m e, m e, m e, m e)) -> (m e, m e, m e, m e, m e, m e, m e)+            toM ((a11, a12, a21, a22), (b11, b12, b21, b22), (c11, c12, c21, c22)) =+                (\[a, b, c, d, e, f, g] -> (a, b, c, d, e, f, g)) . runPar . MP.parMap id $+                    [strassen (zipWith_m (+) a11 a22) (zipWith_m (+) b11 b22) alpha 0 (fillWithZeros c11),+                    strassen (zipWith_m (+) a21 a22) b11 alpha 0 (fillWithZeros c11),+                    strassen a11 (zipWith_m (-) b12 b22) alpha beta c12,+                    strassen a22 (zipWith_m (-) b21 b11) alpha beta c21,+                    strassen (zipWith_m (+) a11 a12) b22 alpha 0 (fillWithZeros c11),+                    strassen (zipWith_m (-) a21 a11) (zipWith_m (+) b11 b12) alpha beta (zipWith_m (-) c22 c12),+                    strassen (zipWith_m (-) a12 a22) (zipWith_m (+) b21 b22) alpha beta (zipWith_m (-) c11 c21)]++            toC :: (m e, m e, m e, m e, m e, m e, m e) -> (m e, m e, m e, m e)+            toC (m1, m2, m3, m4, m5, m6, m7) =+                (\[a, b, c, d] -> (a, b, c, d)) . runPar . MP.parMap id $+                    [zipWith_m (+) m1 . zipWith_m (+) m4 . zipWith_m (-) m7 $ m5,+                    zipWith_m (+) m3 m5,+                    zipWith_m (+) m2 m4,+                    zipWith_m (+) m1 . zipWith_m (+) m3 . zipWith_m (-) m6 $ m2]++            quadJoin :: (m e, m e, m e, m e) -> m e+            quadJoin (c11, c12, c21, c22) =+                fromBlocks_m [+                    [c11, c12],+                    [c21, c22]]
+ src/FPNLA/Operations/BLAS/Strategies/GEMM/Strategies/BlocksPar.hs view
@@ -0,0 +1,38 @@+{-# LANGUAGE FlexibleContexts      #-}+{-# LANGUAGE FlexibleInstances     #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE ScopedTypeVariables   #-}++module FPNLA.Operations.BLAS.Strategies.GEMM.Strategies.BlocksPar () where++import           FPNLA.Matrix                               (fromBlocks_m,+                                                             map_m, toBlocks_m,+                                                             zipWith_m)+import           FPNLA.Operations.BLAS                      (GEMM (gemm))+import           FPNLA.Operations.BLAS.Strategies.DataTypes (BlockContext (getBlockDim),+                                                             BlocksPar_ST)+import           FPNLA.Operations.Parameters                (ResM, TransType (NoTrans),+                                                             blasResultM,+                                                             getResultDataM,+                                                             transTrans_m)++import           Control.DeepSeq                            (NFData)+import           Control.Parallel.Strategies                (parTraversable,+                                                             rdeepseq,+                                                             withStrategy)++instance (NFData (m e), GEMM s m v e) => GEMM (BlocksPar_ST s) m v e where+    gemm (bctx, sctx) tmA tmB alpha beta mC = blasResultM . fromBlocks_m $ generatePar rA cB (\i j -> add_m (matMultIJ i j (bmC !! i !! j)) (map_m (beta *) (bmC !! i !! j)))+        where+            (r, c) = getBlockDim bctx+            add_m = zipWith_m (+)+            generatePar m n f = withStrategy (parTraversable rdeepseq) $ generateBlocks m n f+            generateBlocks m n f = [[f i j | j <- [0 .. (n - 1)]] | i <- [0 .. (m - 1)]]+            matMultIJ i j bmCIJ = foldr1 add_m  [callGemm (bmA !! i !! k) (bmB !! k !! j) bmCIJ |  k <- [0 .. (cA - 1)]]+            callGemm mA mB mC = getResultDataM (gemm sctx (NoTrans mA) (NoTrans mB) alpha 0 mC :: ResM s v m e)+            bmA = toBlocks_m r c $ transTrans_m tmA+            bmB = toBlocks_m r c $ transTrans_m tmB+            bmC = toBlocks_m r c mC+            rA = length bmA+            cA = length (head bmA)+            cB = length (head bmB)
+ src/FPNLA/Operations/BLAS/Strategies/GEMM/Strategies/DefPar.hs view
@@ -0,0 +1,28 @@+{-# LANGUAGE FlexibleContexts      #-}+{-# LANGUAGE FlexibleInstances     #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE ScopedTypeVariables   #-}++module FPNLA.Operations.BLAS.Strategies.GEMM.Strategies.DefPar () where++import           Control.DeepSeq                            (NFData)+import           Control.Parallel.Strategies                (parMap, rdeepseq)+import           FPNLA.Matrix                               (MatrixVector (fromCols_vm),+                                                             elem_m, foldr_v,+                                                             generate_v)+import           FPNLA.Operations.BLAS                      (GEMM (gemm))+import           FPNLA.Operations.BLAS.Strategies.DataTypes (DefPar_ST)+import           FPNLA.Operations.Parameters                (Elt, blasResultM,+                                                             dimTrans_m,+                                                             elemTrans_m)++instance (NFData (v e), Elt e, MatrixVector m v e) => GEMM DefPar_ST m v e where+    gemm _ pmA pmB alpha beta mC+        | p /= p' = error "gemm: incompatible ranges"+        | otherwise = blasResultM $ generatePar_m m n (\i j -> (alpha * matMultIJ i j) + beta * elem_m i j mC)+        where+            (m, p) = dimTrans_m pmA+            (p',n) = dimTrans_m pmB+            matMultIJ i j = foldr_v (+) 0 (generate_v p (\k -> (*) (elemTrans_m i k pmA) (elemTrans_m k j pmB)) :: v e)+            generatePar_m m n gen = fromCols_vm . parMap rdeepseq (\j -> generate_v m (`gen` j) :: v e) $ [0 .. (n - 1)]+
+ src/FPNLA/Operations/BLAS/Strategies/GEMM/Strategies/StrassenPar.hs view
@@ -0,0 +1,77 @@+{-# LANGUAGE FlexibleContexts      #-}+{-# LANGUAGE FlexibleInstances     #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE ScopedTypeVariables   #-}++module FPNLA.Operations.BLAS.Strategies.GEMM.Strategies.StrassenPar () where++import           Control.DeepSeq                            (NFData)+import           Control.Parallel.Strategies                (parTraversable,+                                                             rdeepseq,+                                                             withStrategy)+import           FPNLA.Matrix                               (cantCols_m,+                                                             cantRows_m,+                                                             fromBlocks_m,+                                                             fromBlocks_m,+                                                             generate_m,+                                                             toBlocks_m,+                                                             zipWith_m)+import           FPNLA.Operations.BLAS                      (GEMM (gemm))+import           FPNLA.Operations.BLAS.Strategies.DataTypes (StrassenPar_ST,+                                                             getStrassenLimit)+import           FPNLA.Operations.Parameters                (Elt, ResM, TransType (NoTrans),+                                                             blasResultM,+                                                             getResultDataM,+                                                             transTrans_m)++instance (NFData (m e), Elt e, GEMM s m v e) => GEMM (StrassenPar_ST s) m v e where+    gemm (ctx, sctx) pmA pmB alpha beta mC =+        blasResultM $ strassen mA mB alpha beta mC+        where++            callGemm mA mB alpha beta mC = getResultDataM (gemm sctx (NoTrans mA) (NoTrans mB) alpha beta mC :: ResM s v m e)++            mA = transTrans_m pmA+            mB = transTrans_m pmB++            strassenLimit = getStrassenLimit ctx++            strassen :: m e -> m e -> e -> e -> m e -> m e+            strassen mA mB alpha beta mC =+                if cantRows_m mA <= strassenLimit+                then callGemm mA mB alpha beta mC+                else quadJoin . toC . toM $ (quadPartition mA, quadPartition mB, quadPartition mC)++            quadPartition :: m e -> (m e, m e, m e, m e)+            quadPartition m =+                let+                blocks = toBlocks_m (div (cantRows_m m) 2) (div (cantCols_m m) 2) m in+                (head (head blocks), head blocks !! 1, head (blocks !! 1), blocks !! 1 !! 1)++            fillWithZeros :: m e -> m e+            fillWithZeros m = generate_m (cantRows_m m) (cantCols_m m) (\_ _ -> 0)++            toM :: ((m e, m e, m e, m e), (m e, m e, m e, m e), (m e, m e, m e, m e)) -> (m e, m e, m e, m e, m e, m e, m e)+            toM ((a11, a12, a21, a22), (b11, b12, b21, b22), (c11, c12, c21, c22)) =+                (\[a, b, c, d, e, f, g] -> (a, b, c, d, e, f, g)) . withStrategy (parTraversable rdeepseq) $+                    [strassen (zipWith_m (+) a11 a22) (zipWith_m (+) b11 b22) alpha 0 (fillWithZeros c11),+                    strassen (zipWith_m (+) a21 a22) b11 alpha 0 (fillWithZeros c11),+                    strassen a11 (zipWith_m (-) b12 b22) alpha beta c12,+                    strassen a22 (zipWith_m (-) b21 b11) alpha beta c21,+                    strassen (zipWith_m (+) a11 a12) b22 alpha 0 (fillWithZeros c11),+                    strassen (zipWith_m (-) a21 a11) (zipWith_m (+) b11 b12) alpha beta (zipWith_m (-) c22 c12),+                    strassen (zipWith_m (-) a12 a22) (zipWith_m (+) b21 b22) alpha beta (zipWith_m (-) c11 c21)]++            toC :: (m e, m e, m e, m e, m e, m e, m e) -> (m e, m e, m e, m e)+            toC (m1, m2, m3, m4, m5, m6, m7) =+                (\[a, b, c, d] -> (a, b, c, d)) . withStrategy (parTraversable rdeepseq) $+                    [zipWith_m (+) m1 . zipWith_m (+) m4 . zipWith_m (-) m7 $ m5,+                    zipWith_m (+) m3 m5,+                    zipWith_m (+) m2 m4,+                    zipWith_m (+) m1 . zipWith_m (+) m3 . zipWith_m (-) m6 $ m2]++            quadJoin :: (m e, m e, m e, m e) -> m e+            quadJoin (c11, c12, c21, c22) =+                fromBlocks_m [+                    [c11, c12],+                    [c21, c22]]
+ src/FPNLA/Operations/BLAS/Strategies/GEMV.hs view
@@ -0,0 +1,21 @@+{-# LANGUAGE FlexibleInstances     #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE ScopedTypeVariables   #-}+{-# LANGUAGE UndecidableInstances  #-}++module FPNLA.Operations.BLAS.Strategies.GEMV (++) where++import           FPNLA.Matrix                (asColumn_vm, toCols_vm)+import           FPNLA.Operations.BLAS       (GEMM (gemm), GEMV (gemv))+import           FPNLA.Operations.Parameters (ResM, TransType (..), blasResultV,+                                              getResultDataM)++instance (GEMM s m v e) => GEMV s m v e where+    gemv strat tmA vB alpha beta vC =+        blasResultV . head . toCols_vm. getResultDataM $+            call_gemm tmA pmB alpha beta pmC+        where pmB = NoTrans $ asColumn_vm vB+              pmC = asColumn_vm vC+              call_gemm mA mB alpha beta mC = gemm strat mA mB alpha beta mC :: ResM s v m e
+ src/FPNLA/Operations/BLAS/Strategies/SYRK/DefSeq.hs view
@@ -0,0 +1,33 @@+{-# LANGUAGE FlexibleInstances     #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE ScopedTypeVariables   #-}++module FPNLA.Operations.BLAS.Strategies.SYRK.DefSeq () where+++import           FPNLA.Matrix                               (MatrixVector,+                                                             foldr_v,+                                                             generate_m,+                                                             generate_v)+import           FPNLA.Operations.BLAS                      (SYRK (syrk))+import           FPNLA.Operations.BLAS.Strategies.DataTypes (DefSeq)+import           FPNLA.Operations.Parameters                (Elt,+                                                             TransType (..),+                                                             blasResultM,+                                                             dimTrans_m,+                                                             dimTriang,+                                                             elemSymm,+                                                             elemTrans_m)++instance  (Elt e, MatrixVector m v e) => SYRK DefSeq m v e where+    syrk _ alpha pmA beta pmB+        | p /= p' = error "syrk: incompatible ranges"+        | otherwise = blasResultM $ generate_m p p (\i j -> (alpha * pmAMultIJ i j) + beta * elemSymm i j pmB)+        where+            (p, p') = dimTriang pmB+            matMultIJ i j tmA tmB = foldr_v (+) 0 (generate_v (snd $ dimTrans_m tmA) (\k -> (*) (elemTrans_m i k tmA) (elemTrans_m k j tmB)) :: v e)+            pmAMultIJ i j =+                case pmA of+                    (NoTrans mA) -> matMultIJ i j pmA (Trans mA)+                    (Trans mA) -> matMultIJ i j (Trans mA) pmA+                    (ConjTrans mA) -> matMultIJ i j (Trans mA) pmA
+ src/FPNLA/Operations/BLAS/Strategies/SYRK/MonadPar/DefPar.hs view
@@ -0,0 +1,38 @@+{-# LANGUAGE FlexibleContexts      #-}+{-# LANGUAGE FlexibleInstances     #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE ScopedTypeVariables   #-}++module FPNLA.Operations.BLAS.Strategies.SYRK.MonadPar.DefPar () where+++import           Control.DeepSeq                            (NFData)+import           Control.Monad.Par                          as MP (parMap,+                                                                   runPar)+import           FPNLA.Matrix                               (MatrixVector,+                                                             foldr_v,+                                                             fromCols_vm,+                                                             generate_v)+import           FPNLA.Operations.BLAS                      (SYRK (syrk))+import           FPNLA.Operations.BLAS.Strategies.DataTypes (DefPar_MP)+import           FPNLA.Operations.Parameters                (Elt,+                                                             TransType (..),+                                                             blasResultM,+                                                             dimTrans_m,+                                                             dimTriang,+                                                             elemSymm,+                                                             elemTrans_m)++instance  (NFData (v e), Elt e, MatrixVector m v e) => SYRK DefPar_MP m v e where+    syrk _ alpha pmA beta pmB+        | p /= p' = error "syrk: incompatible ranges"+        | otherwise = blasResultM $ generatePar_m p p (\i j -> (alpha * pmAMultIJ i j) + beta * elemSymm i j pmB)+        where+            (p, p') = dimTriang pmB+            matMultIJ i j tmA tmB = foldr_v (+) 0 (generate_v (snd $ dimTrans_m tmA) (\k -> (*) (elemTrans_m i k tmA) (elemTrans_m k j tmB)) :: v e)+            pmAMultIJ i j =+                case pmA of+                    (NoTrans mA) -> matMultIJ i j pmA (Trans mA)+                    (Trans mA) -> matMultIJ i j (Trans mA) pmA+                    (ConjTrans mA) -> matMultIJ i j (Trans mA) pmA+            generatePar_m m n gen = fromCols_vm . MP.runPar . MP.parMap (\j -> generate_v m (`gen` j) :: v e) $ [0 .. (n - 1)]
+ src/FPNLA/Operations/BLAS/Strategies/SYRK/Strategies/DefPar.hs view
@@ -0,0 +1,37 @@+{-# LANGUAGE FlexibleContexts      #-}+{-# LANGUAGE FlexibleInstances     #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE ScopedTypeVariables   #-}++module FPNLA.Operations.BLAS.Strategies.SYRK.Strategies.DefPar () where+++import           Control.DeepSeq                            (NFData)+import           Control.Parallel.Strategies                (parMap, rdeepseq)+import           FPNLA.Matrix                               (MatrixVector,+                                                             foldr_v,+                                                             fromCols_vm,+                                                             generate_v)+import           FPNLA.Operations.BLAS                      (SYRK (syrk))+import           FPNLA.Operations.BLAS.Strategies.DataTypes (DefPar_ST)+import           FPNLA.Operations.Parameters                (Elt,+                                                             TransType (..),+                                                             blasResultM,+                                                             dimTrans_m,+                                                             dimTriang,+                                                             elemSymm,+                                                             elemTrans_m)++instance  (NFData (v e), Elt e, MatrixVector m v e) => SYRK DefPar_ST m v e where+    syrk _ alpha pmA beta pmB+        | p /= p' = error "syrk: incompatible ranges"+        | otherwise = blasResultM $ generatePar_m p p (\i j -> (alpha * pmAMultIJ i j) + beta * elemSymm i j pmB)+        where+            (p, p') = dimTriang pmB+            matMultIJ i j tmA tmB = foldr_v (+) 0 (generate_v (snd $ dimTrans_m tmA) (\k -> (*) (elemTrans_m i k tmA) (elemTrans_m k j tmB)) :: v e)+            pmAMultIJ i j =+                case pmA of+                    (NoTrans mA) -> matMultIJ i j pmA (Trans mA)+                    (Trans mA) -> matMultIJ i j (Trans mA) pmA+                    (ConjTrans mA) -> matMultIJ i j (Trans mA) pmA+            generatePar_m m n gen = fromCols_vm . parMap rdeepseq (\j -> generate_v m (`gen` j) :: v e) $ [0 .. (n - 1)]
+ src/FPNLA/Operations/BLAS/Strategies/TRSM/BlocksByCols.hs view
@@ -0,0 +1,116 @@+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE ScopedTypeVariables #-}++module FPNLA.Operations.BLAS.Strategies.TRSM.BlocksByCols () where++import           FPNLA.Matrix                               (cantCols_m,+                                                             cantRows_m,+                                                             diagonalBlock,+                                                             fromBlocks_m,+                                                             generate_m,+                                                             subMatrix_m,+                                                             verticalBlock)+import           FPNLA.Operations.BLAS                      (GEMM (gemm),+                                                             TRSM (trsm))+import           FPNLA.Operations.BLAS.Strategies.DataTypes (BlocksByCols, SqrBlockContext (getSqrBlockDim))+import           FPNLA.Operations.Parameters                (ResM, TransType (NoTrans, Trans, ConjTrans), TriangType (Upper, Lower),+                                                             blasResultM,+                                                             getResultDataM,+                                                             transTrans_m,+                                                             unTransT,+                                                             unTriangT, unUnitT)++instance (GEMM gs m v e, TRSM ts m v e) => TRSM (BlocksByCols gs ts) m v e where+    --trsm :: StratCtx s -> e -> TransType (TriangType (UnitType (m e))) -> m e -> Res s v m e+    trsm _ 0 _ mB = blasResultM $ generate_m (cantRows_m mB) (cantCols_m mB) (\_ _ -> 0)+    trsm (ctx, gctx, tctx) alpha tratmA mB =+        blasResultM $ solveTrsmAux alpha (tritAux (transTrans_m (trat mA))) mB []+        where++            (trat, tritmA) = unTransT tratmA++            (_, utmA) = unTriangT tritmA+            (ut, mA) = unUnitT utmA++            tritAux = case (tritmA, tratmA) of+                          (Lower _, NoTrans _) -> Lower+                          (Upper _, Trans _) -> Lower+                          (Upper _, ConjTrans _) -> Lower+                          _ -> Upper++            blockDim = getSqrBlockDim ctx+            callTrsmBase trit ut alpha mA mB = getResultDataM (trsm tctx alpha (NoTrans (trit (ut mA))) mB :: ResM ts v m e)+            callGemmBase alpha mA mB beta mC = getResultDataM (gemm gctx (NoTrans mA) (NoTrans mB) alpha beta mC :: ResM gs v m e)+            solveTrsmAux alpha (tritmA@(Lower _)) mB listX =+                if fin then+                    fromBlocks_m . map (:[]) . reverse $ listX'+                else+                    solveTrsmAux alpha (trit mA) mB' listX'+                where+                    dim_mA = cantRows_m mA+                    (trit, mA) = unTriangT tritmA+                    listX' = mX_i:listX+                    i = length listX+                    fin = blockDim*(i + 1) >= dim_mA+                    tailVerticalBlock m = subMatrix_m blockDim 0 (cantRows_m m - blockDim) (cantCols_m m) m+                    (mX_i, mB')+                        | null listX =+                            let+                                mA_00 = diagonalBlock (getSqrBlockDim ctx, getSqrBlockDim ctx) 0 mA+                                mB_0 = verticalBlock (getSqrBlockDim ctx) 0 mB+                                mB_inf = tailVerticalBlock mB+                            in+                                (callTrsmBase trit ut alpha mA_00 mB_0, mB_inf)+                        | otherwise =+                            let+                                mX = head listX+                                mA_next = subMatrix_m+                                    (blockDim*i)+                                    (blockDim*(i - 1))+                                    (dim_mA - (blockDim*i))+                                    (min blockDim (dim_mA - (blockDim*(i - 1))))+                                    mA+                                mB' = callGemmBase (negate (1/alpha)) mA_next mX 1 mB+                                mB'_inf = tailVerticalBlock mB'+                                mA_ii = diagonalBlock (getSqrBlockDim ctx, getSqrBlockDim ctx) i mA+                                mB_i' = verticalBlock (getSqrBlockDim ctx) 0 mB'+                            in+                                (callTrsmBase trit ut alpha mA_ii mB_i', mB'_inf)+            solveTrsmAux alpha (tritmA@(Upper _)) mB listX =+                if fin then+                    fromBlocks_m . map (:[])  $ listX'+                else+                    solveTrsmAux alpha (trit mA) mB' listX'+                where+                    dim_mA = cantRows_m mA+                    (trit, mA) = unTriangT tritmA+                    mCantBlocksMN m = ceiling ((fromIntegral (cantRows_m m) / fromIntegral blockDim) :: Double)+                    invertBlockIndex i m = mCantBlocksMN m - i -1+                    listX' = mX_i:listX+                    i = length listX+                    fin = blockDim*(i + 1) >= dim_mA+                    torsoVerticalBlock m = subMatrix_m 0 0 (blockDim*(mCantBlocksMN m - 1)) (cantCols_m m) m+                    (mX_i, mB')+                        | null listX =+                            let+                                mA_00 = diagonalBlock (getSqrBlockDim ctx, getSqrBlockDim ctx) (invertBlockIndex 0 mA) mA+                                mB_0 = verticalBlock (getSqrBlockDim ctx) (invertBlockIndex 0 mB) mB+                                mB_sup = torsoVerticalBlock mB+                            in+                                (callTrsmBase trit ut alpha mA_00 mB_0, mB_sup)+                        | otherwise =+                            let+                                mX = head listX+                                mA_next = subMatrix_m+                                    0+                                    (blockDim * invertBlockIndex (i - 1) mA)+                                    (blockDim*(mCantBlocksMN mA - i))+                                    (min blockDim (dim_mA - (blockDim * invertBlockIndex (i - 1) mA))) mA+                                mB' = callGemmBase (negate (1/alpha)) mA_next mX 1 mB+                                mB'_sup = torsoVerticalBlock mB'+                                mA_ii = diagonalBlock (getSqrBlockDim ctx, getSqrBlockDim ctx) (invertBlockIndex i mA) mA+                                mB_i' = verticalBlock (getSqrBlockDim ctx) (invertBlockIndex 0 mB') mB'+                            in+                                (callTrsmBase trit ut alpha mA_ii mB_i', mB'_sup)
+ src/FPNLA/Operations/BLAS/Strategies/TRSM/BlocksByRows.hs view
@@ -0,0 +1,89 @@+{-# LANGUAGE FlexibleContexts      #-}+{-# LANGUAGE FlexibleInstances     #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE ScopedTypeVariables   #-}++module FPNLA.Operations.BLAS.Strategies.TRSM.BlocksByRows () where++import           FPNLA.Matrix                               (cantCols_m,+                                                             cantRows_m,+                                                             diagonalBlock,+                                                             fromBlocks_m,+                                                             generate_m,+                                                             subMatrix_m,+                                                             verticalBlock)+import           FPNLA.Operations.BLAS                      (GEMM (gemm),+                                                             TRSM (trsm))+import           FPNLA.Operations.BLAS.Strategies.DataTypes (BlocksByRows, SqrBlockContext (getSqrBlockDim))+import           FPNLA.Operations.Parameters                (ResM, TransType (NoTrans, Trans, ConjTrans), TriangType (Upper, Lower),+                                                             blasResultM,+                                                             getResultDataM,+                                                             transTrans_m,+                                                             unTransT,+                                                             unTriangT, unUnitT)++instance (GEMM gs m v e, TRSM ts m v e) => TRSM (BlocksByRows gs ts) m v e where+    --trsm :: StratCtx s -> e -> TransType (TriangType (UnitType (m e))) -> m e -> Res s v m e+    trsm _ 0 _ mB = blasResultM $ generate_m (cantRows_m mB) (cantCols_m mB) (\_ _ -> 0)+    trsm (ctx, gctx, tctx) alpha tratmA mB = blasResultM $ solveTrsmAux alpha (tritAux (transTrans_m (trat mA))) mB undefined 0+        where++            (trat, tritmA) = unTransT tratmA+            (_, utmA) = unTriangT tritmA+            (ut, mA) = unUnitT utmA++            tritAux = case (tritmA, tratmA) of+                          (Lower _, NoTrans _) -> Lower+                          (Upper _, Trans _) -> Lower+                          (Upper _, ConjTrans _) -> Lower+                          _ -> Upper++            blockDim = getSqrBlockDim ctx+            callTrsmBase trit ut alpha mA mB = getResultDataM (trsm tctx alpha (NoTrans (trit (ut mA))) mB :: ResM ts v m e)+            callGemmBase alpha mA mB beta mC = getResultDataM (gemm gctx (NoTrans mA) (NoTrans mB) alpha beta mC :: ResM gs v m e)++            solveTrsmAux alpha tritmA mB mX i =+                if fin then+                    mX+                else+                    mX'+                where+                    (trit, mA) = unTriangT tritmA+                    dim_mA = cantRows_m mA+                    mCantBlocksMN = ceiling ((fromIntegral dim_mA / fromIntegral blockDim) :: Double)+                    invertBlockIndex i = case tritmA of+                        (Upper _) -> mCantBlocksMN - i - 1+                        (Lower _) -> i+                    fin = (i*blockDim) >= dim_mA+                    mX'+                        | i == 0 =+                            let+                                mA_00 = diagonalBlock (getSqrBlockDim ctx, getSqrBlockDim ctx) (invertBlockIndex 0) mA+                                mB_0 = verticalBlock (getSqrBlockDim ctx) (invertBlockIndex 0) mB+                                mX_0 = callTrsmBase trit ut alpha mA_00 mB_0+                            in+                                solveTrsmAux alpha (trit mA) mB mX_0 (i + 1)+                        | otherwise =+                            let+                                mA_ii = diagonalBlock (getSqrBlockDim ctx, getSqrBlockDim ctx) (invertBlockIndex i) mA+                                mB_i = verticalBlock (getSqrBlockDim ctx) (invertBlockIndex i) mB+                                mA_sub = case tritmA of+                                    (Upper _) -> subMatrix_m+                                        (blockDim * invertBlockIndex i)+                                        (blockDim*(invertBlockIndex i + 1))+                                        (min blockDim (dim_mA - blockDim * invertBlockIndex i))+                                        (min (blockDim*i) (dim_mA - blockDim*(invertBlockIndex i + 1)))+                                        mA+                                    (Lower _) -> subMatrix_m+                                        (blockDim*i)+                                        0+                                        (min blockDim (dim_mA - blockDim*i))+                                        (blockDim*i)+                                        mA+                                mB'_i = callGemmBase (negate (1/alpha)) mA_sub mX 1 mB_i+                                mX_i = callTrsmBase trit ut alpha mA_ii mB'_i+                                mX_join = case tritmA of+                                    (Upper _) -> fromBlocks_m [[mX_i], [mX]]+                                    (Lower _) -> fromBlocks_m [[mX], [mX_i]]+                            in+                                solveTrsmAux alpha (trit mA) mB mX_join (i + 1)
+ src/FPNLA/Operations/BLAS/Strategies/TRSM/CBindSeq.hs view
@@ -0,0 +1,58 @@+{-# LANGUAGE FlexibleContexts      #-}+{-# LANGUAGE FlexibleInstances     #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE ScopedTypeVariables   #-}++module FPNLA.Operations.BLAS.Strategies.TRSM.CBindSeq () where++import           FPNLA.Matrix                               (Matrix, MatrixVector,+                                                             cantCols_m,+                                                             cantRows_m)+import           FPNLA.Operations.BLAS                      (TRSM(..))+import           FPNLA.Operations.BLAS.Strategies.DataTypes (CBindSeq)+import           FPNLA.Operations.Parameters                (Elt, blasResultM,+                                                             unTransT,+                                                             unTriangT, unUnitT)+import           FPNLA.Operations.Utils                     (matrixToPtr,+                                                             ptrToMatrix,+                                                             transToForeign,+                                                             triangToForeign,+                                                             unitToForeign,+                                                             unsafePerformIO)+++import           Control.DeepSeq                            (NFData (rnf))+import           Control.Exception                          (evaluate)+import           Foreign.BLAS                               as L3 (BLAS3, Side (LeftSide),+                                                                   trsm, Uplo, Trans, Diag)+import           Foreign.Marshal.Alloc                      (free)++instance (BLAS3 e, Elt e, NFData (m e), MatrixVector m v e) => TRSM CBindSeq m v e where+    trsm _ alpha pmA mB =+        let+            (_, pmA') = unTransT pmA+            (_, pmA'') = unTriangT pmA'+            (_, mA) = unUnitT pmA''+            trans = transToForeign pmA+            uplo = triangToForeign pmA'+            diag = unitToForeign pmA''+            (m, n) = (cantRows_m mA, cantCols_m mB)+            (ldA, ldB) = (cantRows_m mA, cantRows_m mB)+        in+            blasResultM $ ioTrsm uplo trans diag m n alpha mA mB ldA ldB++ioTrsm :: (Matrix m a, Matrix m1 a,+           Matrix m2 a, Elt a, BLAS3 a, NFData (m2 a)) =>+           Uplo -> Trans -> Diag -> Int -> Int+           -> a -> m a -> m1 a+           -> Int -> Int -> m2 a+ioTrsm uplo trans diag m n alpha mA mB ldA ldB = unsafePerformIO $ do+        mAptr <- matrixToPtr mA+        mBptr <- matrixToPtr mB+        --L3.trsm  :: Side -> Uplo -> Trans -> Diag -> Int -> Int -> a -> Ptr a -> Int -> Ptr a -> Int -> IO ()+        L3.trsm LeftSide uplo trans diag m n alpha mAptr ldA mBptr ldB+        res <- ptrToMatrix m n mBptr+        evaluate $ rnf res+        free mAptr+        free mBptr+        return res
+ src/FPNLA/Operations/BLAS/Strategies/TRSM/DefSeq.hs view
@@ -0,0 +1,74 @@+{-# LANGUAGE FlexibleInstances     #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE ScopedTypeVariables   #-}++module FPNLA.Operations.BLAS.Strategies.TRSM.DefSeq () where+++import           FPNLA.Matrix                               (MatrixVector,+                                                             dim_m, elem_v,+                                                             fromCols_vm,+                                                             fromList_v,+                                                             toCols_vm)+import           FPNLA.Operations.BLAS                      (TRSM (trsm))+import           FPNLA.Operations.BLAS.Strategies.DataTypes (DefSeq)+import           FPNLA.Operations.Parameters                (Elt, TransType (NoTrans, Trans, ConjTrans), TriangType (Upper, Lower),+                                                             UnitType,+                                                             blasResultM,+                                                             dimTransUnit_m,+                                                             dimUnit_m,+                                                             elemTransUnit_m,+                                                             unTransT,+                                                             unTriangT)++forwardElim :: (Elt e, MatrixVector m v e) => e -> TransType (UnitType (m e)) -> v e -> v e+forwardElim alpha pmA vB = fromList_v res+    where+        cr = fst . dimTransUnit_m $ pmA -- Cantidad de ecuaciones+        res = solve cr -- Usada para la evaluacion perezosa+        solve 0 = []+        solve i = calcX i : solve (i-1)+        calcX i -- Calcula el elemento X_i+            | v_ind == 0  = 0 -- Devuelve 0 si el valor del vector vale 0.+            | otherwise   = (v_ind - sumValues 0 res) / elemTransUnit_m ind ind pmA+            where+                  ind = cr - i -- Indice real+                  ind_v = ind+                  eqr = cr - i -- Cantidad de ecuaciones resueltas+                  v_ind = alpha * elem_v ind_v vB+                  sumValues j (x:xs)+                      | j < eqr = x * elemTransUnit_m ind j pmA + sumValues (j+1) xs+                      | otherwise = 0+--+backwardElim :: (Elt e, MatrixVector m v e) => e -> TransType (UnitType (m e)) -> v e -> v e+backwardElim alpha pmA vB = fromList_v $ reverse res+    where+        cr = fst . dimTransUnit_m $ pmA -- Cantidad de ecuaciones+        res = solve cr -- Usada para la evaluación perezosa+        solve 0 = []+        solve i = calcX i : solve (i-1)+        calcX i -- Calcula el elemento X_[cr - i]+            | v_ind == 0  = 0 -- Devuelve 0 si el valor del vector vale 0.+            | otherwise   = (v_ind - sumValues (cr-1) res) / elemTransUnit_m ind ind pmA+            where+                  ind = i - 1 -- Indice real+                  eqr = cr - i -- Cantidad de ecuaciones resueltas+                  v_ind = alpha * elem_v ind vB+                  sumValues j (x:xs)+                      | j > cr - 1 - eqr = x * elemTransUnit_m ind j pmA + sumValues (j-1) xs+                      | otherwise = 0+--+instance (Elt e, MatrixVector m v e) => TRSM DefSeq m v e where+    trsm _ alpha pmA mB+        | rowsmA /= colsmA = error "trsm: matrix A must be squared"+        | colsmA /= rowsmB = error "trsm: incompatible ranges"+        | otherwise = blasResultM . fromCols_vm . map (selectStrat pmA) $ (toCols_vm mB :: [v e])+        where+            (rowsmA, colsmA) = dimUnit_m . snd . unTriangT . snd . unTransT $ pmA+            rowsmB = fst . dim_m $ mB+            selectStrat (NoTrans (Lower umA)) = forwardElim alpha (NoTrans umA)+            selectStrat (Trans (Lower umA)) = backwardElim alpha (Trans umA)+            selectStrat (ConjTrans (Lower umA)) = backwardElim alpha (ConjTrans umA)+            selectStrat (NoTrans (Upper umA)) = backwardElim alpha (NoTrans umA)+            selectStrat (Trans (Upper umA)) = forwardElim alpha (Trans umA)+            selectStrat (ConjTrans (Upper umA)) = forwardElim alpha (ConjTrans umA)
+ src/FPNLA/Operations/BLAS/Strategies/TRSM/MonadPar/ColsPar.hs view
@@ -0,0 +1,32 @@+{-# LANGUAGE FlexibleContexts      #-}+{-# LANGUAGE FlexibleInstances     #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE ScopedTypeVariables   #-}++module FPNLA.Operations.BLAS.Strategies.TRSM.MonadPar.ColsPar () where++import           Control.DeepSeq                            (NFData)+import           Control.Monad.Par                          as MP (parMap,+                                                                   runPar)+import           FPNLA.Matrix                               (asColumn_vm, dim_m,+                                                             fromBlocks_m,+                                                             toCols_vm)+import           FPNLA.Operations.BLAS                      (TRSM (trsm))+import           FPNLA.Operations.BLAS.Strategies.DataTypes (ColsPar_MP)+import           FPNLA.Operations.Parameters                (ResM, blasResultM,+                                                             dimUnit_m,+                                                             getResultDataM,+                                                             unTransT,+                                                             unTriangT)++-- Optimizacion de la version secuencial resolviendo las columnas en paralelo:+instance (NFData (m e), TRSM s m v e) => TRSM (ColsPar_MP s) m v e where+    trsm ctx alpha pmA mB+        | rowsmA /= colsmA = error "trsm: matrix A must be squared"+        | colsmA /= rowsmB = error "trsm: incompatible ranges"+        | otherwise = blasResultM . fromBlocks_m . (:[]) . MP.runPar . MP.parMap (getResultDataM . callTrsmBase) $ toCols_vm mB+        where+            (rowsmA, colsmA) = dimUnit_m . snd . unTriangT . snd . unTransT $ pmA+            rowsmB = fst . dim_m $ mB+            callTrsmBase :: v e -> ResM s v m e+            callTrsmBase vB = trsm ctx alpha pmA $ asColumn_vm vB
+ src/FPNLA/Operations/BLAS/Strategies/TRSM/Strategies/ColsPar.hs view
@@ -0,0 +1,31 @@+{-# LANGUAGE FlexibleContexts      #-}+{-# LANGUAGE FlexibleInstances     #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE ScopedTypeVariables   #-}++module FPNLA.Operations.BLAS.Strategies.TRSM.Strategies.ColsPar () where++import           Control.DeepSeq                            (NFData)+import           Control.Parallel.Strategies                (parMap, rdeepseq)+import           FPNLA.Matrix                               (asColumn_vm, dim_m,+                                                             fromBlocks_m,+                                                             toCols_vm)+import           FPNLA.Operations.BLAS                      (TRSM (trsm))+import           FPNLA.Operations.BLAS.Strategies.DataTypes (ColsPar_ST)+import           FPNLA.Operations.Parameters                (ResM, blasResultM,+                                                             dimUnit_m,+                                                             getResultDataM,+                                                             unTransT,+                                                             unTriangT)++-- Optimizacion de la version secuencial resolviendo las columnas en paralelo:+instance (NFData (m e), TRSM s m v e) => TRSM (ColsPar_ST s) m v e where+    trsm ctx alpha pmA mB+        | rowsmA /= colsmA = error "trsm: matrix A must be squared"+        | colsmA /= rowsmB = error "trsm: incompatible ranges"+        | otherwise = blasResultM . fromBlocks_m . (:[]) . parMap rdeepseq (getResultDataM . callTrsmBase) $ toCols_vm mB+        where+            (rowsmA, colsmA) = dimUnit_m . snd . unTriangT . snd . unTransT $ pmA+            rowsmB = fst . dim_m $ mB+            callTrsmBase :: v e -> ResM s v m e+            callTrsmBase vB = trsm ctx alpha pmA $ asColumn_vm vB
+ src/FPNLA/Operations/LAPACK/Strategies.hs view
@@ -0,0 +1,12 @@+module FPNLA.Operations.LAPACK.Strategies (+    module FPNLA.Operations.LAPACK.Strategies.DataTypes+) where++import FPNLA.Operations.LAPACK.Strategies.DataTypes++import FPNLA.Operations.LAPACK.Strategies.POTRF.DefSeq ()+import FPNLA.Operations.LAPACK.Strategies.POTRF.Repa.DefPar ()+import FPNLA.Operations.LAPACK.Strategies.POTRF.BlocksSeq ()+import FPNLA.Operations.LAPACK.Strategies.POTRF.Repa.BlocksPar ()+import FPNLA.Operations.LAPACK.Strategies.POTRF.HMatrixBindSeq ()+
+ src/FPNLA/Operations/LAPACK/Strategies/DataTypes.hs view
@@ -0,0 +1,45 @@+{-# LANGUAGE EmptyDataDecls        #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE TypeFamilies          #-}++module FPNLA.Operations.LAPACK.Strategies.DataTypes (++    HMatrixBindSeq,+    CholLLVBlocksSeq,+    CholLLVBlocksPar_Repa,+    CholLLVSeq,+    CholLLVPar_Repa,++    NullContext,++    SqrBlockContext,+    getSqrBlockDim,++    newNullContext,+    newSqrBlockContext,++) where+++import           FPNLA.Operations.Parameters                (StratCtx)++-- TODO:Ver como sacar esto!!+import           FPNLA.Operations.BLAS.Strategies.DataTypes++data HMatrixBindSeq++data CholLLVSeq dots gemvs+data CholLLVPar_Repa dots gemvs++data CholLLVBlocksSeq potrfs syrks gemms trsms+data CholLLVBlocksPar_Repa potrfs syrks gemms trsms++type instance StratCtx HMatrixBindSeq = NullContext++type instance StratCtx (CholLLVSeq dots gemvs) = (StratCtx dots, StratCtx gemvs)+type instance StratCtx (CholLLVPar_Repa dots gemvs) = (StratCtx dots, StratCtx gemvs)++type instance StratCtx (CholLLVBlocksSeq syrks gemms trsms potrfs) =+    (SqrBlockContext, StratCtx syrks, StratCtx gemms, StratCtx trsms, StratCtx potrfs)+type instance StratCtx (CholLLVBlocksPar_Repa syrks gemms trsms potrfs) =+    (SqrBlockContext, StratCtx syrks, StratCtx gemms, StratCtx trsms, StratCtx potrfs)
+ src/FPNLA/Operations/LAPACK/Strategies/POTRF/BlocksSeq.hs view
@@ -0,0 +1,86 @@+{-# LANGUAGE FlexibleInstances     #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE ScopedTypeVariables   #-}++module FPNLA.Operations.LAPACK.Strategies.POTRF.BlocksSeq (++) where+++import           FPNLA.Matrix                                 (Matrix,+                                                               MatrixVector,+                                                               cantCols_m,+                                                               cantRows_m,+                                                               elem_m,+                                                               fromCols_vm,+                                                               generate_m,+                                                               subMatrix_m,+                                                               toCols_vm,+                                                               transpose_m)+import           FPNLA.Operations.BLAS                        (GEMM (gemm),+                                                               SYRK (syrk),+                                                               TRSM (trsm))+import           FPNLA.Operations.LAPACK                      (POTRF (potrf))+import           FPNLA.Operations.LAPACK.Strategies.DataTypes (CholLLVBlocksSeq,+                                                               SqrBlockContext,+                                                               getSqrBlockDim)+import           FPNLA.Operations.Parameters                  (Elt, ResM,+                                                               TransType (..),+                                                               TriangType (..),+                                                               UnitType (..),+                                                               blasResultM,+                                                               getResultDataM)+import           FPNLA.Utils                                  (iif)++--import Debug.Trace++--trace' s a = trace (s ++ ": " ++ (show a)) a++instance (Elt e, MatrixVector m v e, POTRF potrfs m v e, SYRK syrks m v e,+         GEMM gemms m v e, TRSM trsms m v e) =>+         POTRF (CholLLVBlocksSeq syrks gemms trsms potrfs) m v e where+    potrf ctx (Lower mA)+        =  blasResultM $ chol_blk_l ctx 0 mA (generate_m 0 0 undefined)+        where+            chol_blk_l ctx@(block_ctx, syrk_ctx, gemm_ctx, trsm_ctx, potrf_ctx) k mA mAL+                | k == 0           =  let  mA11'  = call_chol_unb (Lower mA11)+                                           mA21'  = transpose_m $ call_trsm 1 (NoTrans . Lower $ NoUnit mA11') (transpose_m mA21)+                                           mAx1   = iif (k == cantBlocks - 1) mA11' $ concatByCol_m mA11' mA21'+                                      in chol_blk_l ctx (k + 1) mA mAx1+                | k == cantBlocks  =  mAL+                | otherwise        =  chol_blk_l ctx (k + 1) mA mAL'++                where+                    mA_dim      = cantCols_m mA+                    mA10        = subMatrix_m (block * k) 0 block (block*k) mAL+                    mA11        = subMatrix_m (block * k) (block * k) block block mA+                    mA11'       = call_syrk (-1) (NoTrans mA10) 1 (Lower mA11)+                    mA11''      = call_chol_unb (Lower mA11')+                    mA20        = subMatrix_m ((k + 1) * block) 0 (mA_dim - (k + 1)*block) (k*block) mAL+                    mA21        = subMatrix_m ((k + 1) * block) (k * block) (mA_dim - (k + 1)*block) block mA+                    mA21'       = call_gemm  (NoTrans mA20) (Trans mA10) (-1) 1 mA21+                    mA21''      = transpose_m $ call_trsm 1 (NoTrans . Lower $ NoUnit mA11'') (transpose_m mA21')++                    mAx1        = add_zeros k block . iif (k == cantBlocks - 1) mA11'' $ concatByCol_m mA11'' mA21''++                    mAL'        = fromCols_vm $ (toCols_vm mAL :: [v e]) ++ (toCols_vm mAx1 :: [v e])++                    call_syrk m1 alpha beta m2 = getResultDataM (syrk syrk_ctx m1 alpha beta m2 :: ResM syrks v m e)+                    call_chol_unb m = getResultDataM (potrf potrf_ctx m :: ResM potrfs v m e)+                    call_gemm m1 m2 alpha beta m3 = getResultDataM (gemm gemm_ctx m1 m2 alpha beta m3 :: ResM gemms v m e)+                    call_trsm alpha m1 m2 = getResultDataM (trsm trsm_ctx alpha m1 m2 :: ResM trsms v m e)++                    add_zeros :: (Matrix m e) => Int -> Int -> m e -> m e+                    add_zeros k block = concatByCol_m (generate_m (k*block) block (\_ _ -> 0) :: m e)++                    block       = getSqrBlockDim block_ctx+                    cantBlocks  = mA_dim `div` block++concatByCol_m :: (Matrix m e) => m e -> m e -> m e+concatByCol_m m1 m2 = let rows_m1 = cantRows_m m1+                          cols_m1 = cantCols_m m1+                          rows_m2 = cantRows_m m2+                      in generate_m (rows_m1 + rows_m2) cols_m1+                                    (\i j -> iif (i >= rows_m1) (elem_m (i - rows_m1) j m2)+                                                                (elem_m i j m1))+
+ src/FPNLA/Operations/LAPACK/Strategies/POTRF/DefSeq.hs view
@@ -0,0 +1,119 @@+{-# LANGUAGE FlexibleInstances     #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE ScopedTypeVariables   #-}+--{-# LANGUAGE FlexibleContexts #-} --TODO sacar al terminar de debugear junto con el Show++module FPNLA.Operations.LAPACK.Strategies.POTRF.DefSeq (++) where+++import           FPNLA.Matrix                                 (MatrixVector,+                                                               cantCols_m,+                                                               concat_v, elem_m,++                                                               fromCols_vm,+                                                               fromList_v,+                                                               generate_v,+                                                               map_v,+                                                               subMatrix_m,+                                                               toCols_vm)+import           FPNLA.Operations.BLAS                        (DOT (dot),+                                                               GEMV (gemv))+import           FPNLA.Operations.LAPACK                      (POTRF (potrf))+import           FPNLA.Operations.LAPACK.Strategies.DataTypes (CholLLVSeq)+import           FPNLA.Operations.Parameters                  (Elt, ResV, ResS,+                                                               TransType (..),+                                                               TriangType (..),+                                                               blasResultM,+                                                               getConjugate,+                                                               getResultDataV,+                                                               getResultDataS,+                                                               unTriangT)+import           FPNLA.Utils                                  (iif)+--import Debug.Trace++{-+    Yo entiendo que la implementacion de lapack es el "The Cholesky Banachiewicz+    and Cholesky Crout algorithms", que seria lo mismo al "Left-looking variant".+-}+++instance (Elt e, MatrixVector m v e, DOT dots v e, GEMV gemvs m v e) =>+    POTRF (CholLLVSeq dots gemvs) m v e where+    potrf (dot_ctx, gemv_ctx) pmA =+        case pmA of+            (Lower _) -> blasResultM $ potrf_l 0 undefined+            (Upper _) -> undefined+        where (_, mA) = unTriangT pmA+              mA_dim = cantCols_m mA+              call_gemv m v1 alpha beta v2 = getResultDataV (gemv gemv_ctx (NoTrans m) v1 alpha beta v2 :: ResV gemvs v e)+              potrf_l :: Int -> m e -> m e+              potrf_l (k@0) _  = let mAkk = sqrt $ elem_m k k mA+                                     col_mAk :: v e = generate_v (mA_dim - k - 1) (\i -> elem_m (i + k + 1) k mA)+                                     col_mLk :: v e = iif (k == mA_dim - 1) (fromList_v [mAkk]) $+                                                      concat_v [fromList_v [mAkk], map_v (/mAkk) col_mAk]+                                 in potrf_l (k + 1) $ fromCols_vm [col_mLk]+              potrf_l k mL+                  | k == mA_dim = mL+                  | otherwise = let row_mLk :: v e = generate_v k (\j -> elem_m k j mL)+                                    row_mLk_conj = map_v getConjugate row_mLk+                                    mAkk = sqrt $ elem_m k k mA - getResultDataS (dot dot_ctx row_mLk row_mLk_conj :: ResS dots e)+                                    col_mAk :: v e = generate_v (mA_dim - k - 1) (\i -> elem_m (i + k + 1) k mA)+                                    sub_mL = subMatrix_m (k + 1) 0 (mA_dim - k - 1) k mL+                                    col_mAk' = call_gemv sub_mL row_mLk (-1) 1 col_mAk+                                    col_mLk :: v e = add_zeros . iif (k == mA_dim - 1) (fromList_v [mAkk]) $+                                                     concat_v [fromList_v [mAkk], map_v (/mAkk) col_mAk']+                                    add_zeros v = concat_v [generate_v k (const 0) , v]+                                in potrf_l (k + 1) $ fromCols_vm $ toCols_vm mL ++ [col_mLk]+++++{-+import Lapack (Res, ResV, Elt, getConjugate, LapackPotrfOp(potrf), BlasDotOp(dot), BlasGemvOp(gemv),+              blasResult, getResultDataV, TriangType(..), unTriangT, TransType(..))+import Matrix (MatrixVector, cantCols_m, generate_v, generate_m, toCols_vm, fromCols_vm,+              subMatrix_m, elem_m, elem_v, fromList_v, concat_v, map_v)+import Lapack.Strategies.DataTypes (CholLLVSeq)+import Utils (iif)++--import Debug.Trace++--trace' s a = trace (s ++ ": " ++ (show a)) a++instance (Elt e, MatrixVector m v e, BlasDotOp v e, BlasGemvOp gemvs m v e) =>+    LapackPotrfOp (CholLLVSeq gemvs) m v e where+    potrf gemv_ctx (Lower mA)+        =  blasResult $ chol_l gemv_ctx 0 mA (generate_m 0 0 undefined)+        where+            chol_l gemv_ctx k mA mAL+                | k == 0  =  let eA11'  = sqrt $ eA11+                                 vA11'  = fromList_v [eA11']+                                 vA21'  = map_v (/eA11') mA21+                                 vAx1   = iif (k == mA_dim - 1) vA11' (concat_v [vA11', vA21'])+                                 mAL'   = fromCols_vm [vAx1]+                             in chol_l gemv_ctx (k + 1) mA mAL'+                | k == mA_dim  =  mAL+                | otherwise    =  chol_l gemv_ctx (k + 1) mA mAL'++                where+                    vA10    = generate_v k (\j -> elem_m k j mAL)+                    vA10_c  = map_v getConjugate vA10+                    eA11    = elem_m k k mA+                    eA11'   = sqrt $ eA11 - (dot vA10 vA10_c)+                    vA11'   = fromList_v [eA11']+                    mA21    = generate_v (mA_dim - k - 1) (\i -> elem_m (i + k + 1) k mA)+                    mA20    = subMatrix_m (k + 1) 0 (mA_dim - k - 1) k mAL+                    vA21'   = call_gemv (NoTrans mA20) vA10 (-1) 1 mA21++                    vAx1  = add_zeros . iif (k == mA_dim - 1) vA11' $+                                     concat_v [vA11', map_v (/eA11') vA21']++                    mAL'  = fromCols_vm $ (toCols_vm mAL) ++ [vAx1]++                    call_gemv m v1 alpha beta v2 = getResultDataV $ ((gemv gemv_ctx m v1 alpha beta v2) :: ResV gemvs v e)++                    add_zeros v  = concat_v [generate_v k (\i -> 0) , v]+                    mA_dim       = cantCols_m mA+-}
+ src/FPNLA/Operations/LAPACK/Strategies/POTRF/HMatrixBindSeq.hs view
@@ -0,0 +1,28 @@+{-# LANGUAGE FlexibleContexts      #-}+{-# LANGUAGE FlexibleInstances     #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE ScopedTypeVariables   #-}++module FPNLA.Operations.LAPACK.Strategies.POTRF.HMatrixBindSeq (++) where++import qualified Data.Packed.Matrix                           as HM+import           FPNLA.Matrix                                 (MatrixVector,+                                                               cantCols_m,+                                                               cantRows_m,+                                                               elem_m,+                                                               generate_m,+                                                               transpose_m)+import           FPNLA.Operations.LAPACK                      (POTRF (potrf))+import           FPNLA.Operations.LAPACK.Strategies.DataTypes (HMatrixBindSeq)+import           FPNLA.Operations.Parameters                  (Elt, blasResultM,+                                                               unTriangT)+import           Numeric.LinearAlgebra.Algorithms++instance (Field e, Elt e, MatrixVector m v e) => POTRF HMatrixBindSeq m v e where+    potrf _ pmA = blasResultM . fromHMatrix . chol . toHMatrix $ mA+        where mA = snd $ unTriangT pmA+              fromHMatrix m = transpose_m $ generate_m (HM.rows m) (HM.cols m) (\i j -> m HM.@@> (i, j))+              toHMatrix m = HM.buildMatrix (cantRows_m m) (cantCols_m m) (\(i, j)-> elem_m i j m)+
+ src/FPNLA/Operations/LAPACK/Strategies/POTRF/Repa/BlocksPar.hs view
@@ -0,0 +1,132 @@+{-# LANGUAGE BangPatterns          #-}+{-# LANGUAGE FlexibleContexts      #-}+{-# LANGUAGE FlexibleInstances     #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE ScopedTypeVariables   #-}+--{-# LANGUAGE FlexibleContexts #-} --TODO sacar al terminar de debugear junto con el Show++module FPNLA.Operations.LAPACK.Strategies.POTRF.Repa.BlocksPar (++) where+++import           FPNLA.Matrix                                 (Matrix,+                                                               MatrixVector,+                                                               cantCols_m,+                                                               cantRows_m,+                                                               elem_m,+                                                               fromCols_vm,+                                                               generate_m,+                                                               subMatrix_m,+                                                               toCols_vm,+                                                               transpose_m)+import           FPNLA.Matrix.Instances.RepaMatrix            (RepaMatrix,+                                                               RepaVector)+import           FPNLA.Operations.BLAS                        (GEMM (gemm),+                                                               SYRK (syrk),+                                                               TRSM (trsm))+import           FPNLA.Operations.LAPACK                      (POTRF (potrf))+import           FPNLA.Operations.LAPACK.Strategies.DataTypes (CholLLVBlocksPar_Repa,+                                                               SqrBlockContext,+                                                               getSqrBlockDim)+import           FPNLA.Operations.Parameters                  (Elt, ResM,+                                                               TransType (..),+                                                               TriangType (..),+                                                               UnitType (..),+                                                               blasResultM,+                                                               getResultDataM)+import           FPNLA.Utils                                  (iif)++import           Data.Array.Repa                              (Array, D, Shape,+                                                               U, deepSeqArray,+                                                               delay)+import           Data.Array.Repa.Eval                         (suspendedComputeP)+import           Data.Array.Repa.Repr.Unboxed                 (Unbox)++import           Control.DeepSeq                              (NFData)++--import Debug.Trace++--trace' s a = trace (s ++ ": " ++ (show a)) a++instance (Elt e, Unbox e, NFData e,+         MatrixVector RepaMatrix RepaVector e,+         POTRF potrfs RepaMatrix RepaVector e,+         SYRK syrks RepaMatrix RepaVector e,+         GEMM gemms RepaMatrix RepaVector e,+         TRSM trsms RepaMatrix RepaVector e) =>+         POTRF (CholLLVBlocksPar_Repa syrks gemms trsms potrfs) RepaMatrix RepaVector e where+    potrf ctx (Lower mA)+        =  blasResultM $ chol_blk_l ctx 0 mA (generate_m 0 0 undefined)+        where+            chol_blk_l ctx@(block_ctx, syrk_ctx, gemm_ctx, trsm_ctx, potrf_ctx) k mA mAL+                | k == 0           =  let  mA11'  = delay . computeUnboxedP $+                                                    call_chol_unb (Lower mA11)+                                           mA21'  = delay . computeUnboxedP .+                                                    deepSeqArray mA11' .+                                                    deepSeqArray mA21 $+                                                    transpose_m $ call_trsm 1 (NoTrans . Lower $ NoUnit mA11') (transpose_m mA21)+                                           mAx1   = delay . computeUnboxedP .+                                                    deepSeqArray mA21 $+                                                    iif (k == cantBlocks - 1) mA11' $ concatByCol_m mA11' mA21'+                                      in chol_blk_l ctx (k + 1) mA mAx1+                | k == cantBlocks  =  mAL+                | otherwise        =  chol_blk_l ctx (k + 1) mA mAL'++                where++                    {-# INLINE mA_dim #-}+                    mA_dim      = cantCols_m mA+                    {-# INLINE block #-}+                    block       = getSqrBlockDim block_ctx+                    {-# INLINE cantBlocks #-}+                    cantBlocks  = mA_dim `div` block++                    mA10        = delay . computeUnboxedP .+                                  deepSeqArray mAL $+                                  subMatrix_m (block * k) 0 block (block*k) mAL+                    mA11        = delay . computeUnboxedP $+                                  subMatrix_m (block * k) (block * k) block block mA+                    mA11'       = delay . computeUnboxedP .+                                  deepSeqArray mA10+                                  deepSeqArray mA11 $+                                  call_syrk (-1) (NoTrans mA10) 1 (Lower mA11)+                    mA11''      = delay . computeUnboxedP .+                                  deepSeqArray mA11' $+                                  call_chol_unb (Lower mA11')+                    mA20        = delay . computeUnboxedP $+                                  subMatrix_m ((k + 1) * block) 0 (mA_dim - (k + 1)*block) (k*block) mAL+                    mA21        = delay . computeUnboxedP $+                                  subMatrix_m ((k + 1) * block) (k * block) (mA_dim - (k + 1)*block) block mA+                    mA21'       = delay . computeUnboxedP .+                                  deepSeqArray mA20 .+                                  deepSeqArray mA21 $+                                  call_gemm  (NoTrans mA20) (Trans mA10) (-1) 1 mA21+                    mA21''      = delay . computeUnboxedP .+                                  transpose_m $ call_trsm 1 (NoTrans . Lower $ NoUnit mA11'') (transpose_m mA21')++                    mAx1        = delay . computeUnboxedP .+                                  add_zeros k block . iif (k == cantBlocks - 1) mA11'' $ concatByCol_m mA11'' mA21''++                    mAL'        = delay . computeUnboxedP .+                                  fromCols_vm $ (toCols_vm mAL :: [RepaVector e]) ++ (toCols_vm mAx1 :: [RepaVector e])++                    call_syrk !m1 !alpha !beta !m2 = getResultDataM $ (syrk syrk_ctx m1 alpha beta m2 :: ResM syrks RepaVector RepaMatrix e)+                    call_chol_unb !m = getResultDataM $ (potrf potrf_ctx m :: ResM potrfs RepaVector RepaMatrix e)+                    call_gemm !m1 !m2 !alpha !beta !m3 = getResultDataM $ (gemm gemm_ctx m1 m2 alpha beta m3 :: ResM gemms RepaVector RepaMatrix e)+                    call_trsm !alpha !m1 !m2 = getResultDataM $ (trsm trsm_ctx alpha m1 m2 :: ResM trsms RepaVector RepaMatrix e)++                    add_zeros :: (Matrix RepaMatrix e) => Int -> Int -> RepaMatrix e -> RepaMatrix e+                    add_zeros !k !block !v = concatByCol_m (generate_m (k*block) (block) (\_ _ -> 0) :: RepaMatrix e) v++concatByCol_m :: (Unbox e, Matrix RepaMatrix e) => RepaMatrix e -> RepaMatrix e -> RepaMatrix e+concatByCol_m !m1 !m2 = let rows_m1 = cantRows_m m1+                            cols_m1 = cantCols_m m1+                            rows_m2 = cantRows_m m2+                            cols_m2 = cantCols_m m2+                        in generate_m (rows_m1 + rows_m2) cols_m1+                                      (\i j -> iif (i >= rows_m1) (elem_m (i - rows_m1) j m2)+                                                                  (elem_m i j m1))++computeUnboxedP :: (Unbox e, Shape sh) => Array D sh e -> Array U sh e+computeUnboxedP = suspendedComputeP
+ src/FPNLA/Operations/LAPACK/Strategies/POTRF/Repa/DefPar.hs view
@@ -0,0 +1,113 @@+{-# LANGUAGE BangPatterns          #-}+{-# LANGUAGE FlexibleContexts      #-}+{-# LANGUAGE FlexibleInstances     #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE ScopedTypeVariables   #-}+--{-# LANGUAGE FlexibleContexts #-} --TODO sacar al terminar de debugear junto con el Show++module FPNLA.Operations.LAPACK.Strategies.POTRF.Repa.DefPar (++) where++import           FPNLA.Matrix                                 (MatrixVector,+                                                               cantCols_m,+                                                               concat_v, elem_m,++                                                               fromCols_vm,+                                                               fromList_v,+                                                               generate_m,+                                                               generate_v,+                                                               map_v,+                                                               subMatrix_m,+                                                               toCols_vm)+import           FPNLA.Matrix.Instances.RepaMatrix            (RepaMatrix,+                                                               RepaVector)+import           FPNLA.Operations.BLAS                        (DOT (dot),+                                                               GEMV (gemv))+import           FPNLA.Operations.LAPACK                      (POTRF (potrf))+import           FPNLA.Operations.LAPACK.Strategies.DataTypes (CholLLVPar_Repa)+import           FPNLA.Operations.Parameters                  (Elt, ResS, ResV,+                                                               TransType (..),+                                                               TriangType (..),+                                                               blasResultM,+                                                               getConjugate,+                                                               getResultDataS,+                                                               getResultDataV)+import           FPNLA.Utils                                  (iif)++import           Data.Array.Repa                              (Array, D, Shape,+                                                               U, deepSeqArray,+                                                               delay)+import           Data.Array.Repa.Eval                         (suspendedComputeP)+import           Data.Array.Repa.Repr.Unboxed                 (Unbox)++import           Control.DeepSeq                              (NFData, deepseq)++--import Debug.Trace++--trace' s a = trace (s ++ ": " ++ (show a)) a++instance (Elt e, Unbox e, NFData e,+          MatrixVector RepaMatrix RepaVector e,+          DOT dots RepaVector e,+          GEMV gemvs RepaMatrix RepaVector e) =>+          POTRF (CholLLVPar_Repa dots gemvs) RepaMatrix RepaVector e where+    potrf (dot_ctx, gemv_ctx) (Lower mA)+        =  blasResultM $ chol_l gemv_ctx 0 mA (generate_m 0 0 undefined)+        where+            chol_l gemv_ctx k mA mAL+                | k == 0  =  let {-# INLINE eA11' #-}+                                 eA11'  = sqrt eA11+                                 vA11'  = fromList_v [eA11']+                                 vA21'  =+                                          map_v (/eA11') mA21+                                 vAx1   =+                                          iif (k == mA_dim - 1) vA11' (concat_v [vA11', vA21'])+                                 mAL'   =+                                          fromCols_vm [vAx1]+                             in chol_l gemv_ctx (k + 1) mA mAL'+                | k == mA_dim  =  mAL+                | otherwise    =  chol_l gemv_ctx (k + 1) mA mAL'++                where+                    vA10    =  delay . computeUnboxedP .+                               deepSeqArray mAL $+                               generate_v k (\j -> elem_m k j mAL)+                    vA10_c  =+                               map_v getConjugate vA10+                    {-# INLINE vA10_p #-}+                    vA10_p  =  deepSeqArray vA10_c $+                               getResultDataS (dot dot_ctx vA10 vA10_c :: ResS dots e)+                    {-# INLINE eA11 #-}+                    eA11    =  elem_m k k mA+                    {-# INLINE eA11' #-}+                    eA11'   =  sqrt $ eA11 - vA10_p+                    vA11'   =  delay . computeUnboxedP .+                               deepseq eA11' $+                               fromList_v [eA11']+                    mA21    =  delay . computeUnboxedP $+                               generate_v (mA_dim - k - 1) (\i -> elem_m (i + k + 1) k mA)+                    mA20    =  delay . computeUnboxedP $+                               subMatrix_m (k + 1) 0 (mA_dim - k - 1) k mAL+                    vA21'   =  delay . computeUnboxedP .+                               deepSeqArray mA20 .+                               deepSeqArray mA21 $+                               call_gemv (NoTrans mA20) vA10 (-1) 1 mA21++                    vAx1  = delay . computeUnboxedP .+                            deepSeqArray vA11' .+                            deepSeqArray vA21' .+                            add_zeros $ iif (k == mA_dim - 1) vA11' (concat_v [vA11', map_v (/eA11') vA21'])++                    mAL'  = delay . computeUnboxedP .+                            deepSeqArray vAx1 .+                            fromCols_vm $ toCols_vm mAL ++ [vAx1]++                    call_gemv !m !v1 !alpha !beta !v2 = getResultDataV (gemv gemv_ctx m v1 alpha beta v2 :: ResV gemvs RepaVector e)++                    add_zeros !v  = concat_v [generate_v k (const 0) , v]+                    {-# INLINE mA_dim #-}+                    mA_dim       = cantCols_m mA++computeUnboxedP :: (Unbox e, Shape sh) => Array D sh e -> Array U sh e+computeUnboxedP = suspendedComputeP
+ src/FPNLA/Operations/Utils.hs view
@@ -0,0 +1,43 @@++module FPNLA.Operations.Utils (+    matrixToPtr,+    ptrToMatrix,+    transToForeign,+    triangToForeign,+    unitToForeign,+    unsafePerformIO,+) where++import FPNLA.Matrix (Matrix(dim_m, elem_m, generate_m))+import FPNLA.Operations.Parameters (UnitType(NoUnit, Unit), TransType(NoTrans, Trans, ConjTrans), TriangType(Upper, Lower), Elt)++import Foreign.Storable (Storable, peekElemOff)+import Foreign.Marshal.Array (newArray)+import Foreign.Ptr (Ptr)+import System.IO.Unsafe (unsafePerformIO)++import qualified Foreign.BLAS as L3 (Trans(ConjTrans, NoTrans, Trans), Uplo(Lower, Upper), Diag(Unit, NonUnit))++-- ----------------------------------------------------------------------------+-- CBind:+-- ----------------------------------------------------------------------------+matrixToPtr :: (Storable e, Elt e, Matrix m e) => m e -> IO (Ptr e)+matrixToPtr m = do+    let (rs, cs) = dim_m m+    newArray [elem_m i j m | j <- [0..cs-1], i <- [0..rs-1]]++ptrToMatrix :: (Storable e, Elt e, Matrix m e) => Int -> Int -> Ptr e -> IO (m e)+ptrToMatrix m n ptr = return $ generate_m m n (\i j -> unsafePerformIO $ peekElemOff ptr $ j*m+i)++transToForeign :: TransType t -> L3.Trans+transToForeign (NoTrans _)    = L3.NoTrans+transToForeign (Trans _)      = L3.Trans+transToForeign (ConjTrans _)  = L3.ConjTrans++triangToForeign :: TriangType t -> L3.Uplo+triangToForeign (Lower _) = L3.Lower+triangToForeign (Upper _) = L3.Upper++unitToForeign :: UnitType t -> L3.Diag+unitToForeign (Unit _) = L3.Unit+unitToForeign (NoUnit _) = L3.NonUnit
+ src/FPNLA/Utils.hs view
@@ -0,0 +1,39 @@+{-# LANGUAGE FlexibleInstances    #-}+{-# LANGUAGE IncoherentInstances  #-}+{-# LANGUAGE ScopedTypeVariables  #-}+{-# LANGUAGE UndecidableInstances #-}++module FPNLA.Utils (+    iif,+    mapPair,+    splitSlice,++    prevEnum,+    nextEnum,+    enumerate,+) where++mapPair :: (a -> b) -> (c -> d)-> (a, c) -> (b, d)+mapPair f g (a, c) = (f a,  g c)++-- inline if+iif :: Bool -> a -> a -> a+iif cond expTrue expFalse = if cond then expTrue else expFalse++-- Parte la lista parametro en trozos de tamaño 'n'. El último puede ser mas pequeño si no existe k tal que length ls = k * n.+splitSlice :: Int -> [a] -> [[a]]+splitSlice 0 _ = error "El tamaño de n debe ser > 0"+splitSlice _ [] = []+splitSlice n ls = sl : splitSlice n sls+    where (sl, sls) = splitAt n ls++prevEnum :: (Eq a, Bounded a, Enum a) => a -> a+prevEnum a | a == minBound = maxBound+           | otherwise = pred a+nextEnum :: (Eq a, Bounded a, Enum a) => a -> a+nextEnum a | a == maxBound = minBound+           | otherwise = succ a++-- Devuelve una lista con los elementos de cualquier tipo que sea instancia de Enum y Bounded.+enumerate :: forall a. (Enum a, Bounded a) => [a]+enumerate = map toEnum [fromEnum (minBound :: a) .. fromEnum (maxBound :: a)]