text:
-----------------------------------------------------------------------------
-- |
-- Module      :  Finance.Hqfl.Pricer.BlackScholes
-- Copyright   :  (C) 2016 Mika'il Khan
-- License     :  (see the file LICENSE)
-- Maintainer  :  Mika'il Khan <co.kleisli@gmail.com>
-- Stability   :  stable
-- Portability :  portable
--
----------------------------------------------------------------------------
{-# LANGUAGE FlexibleInstances #-}

module Finance.Hqfl.Pricer.BlackScholes where

import Finance.Hqfl.Instrument
import Statistics.Distribution.Normal
import Data.Random

class BS a where
  price :: a -> Double -> Double -> Double

instance BS (Option Equity) where
  price (Option (Equity s q) m European k t) r v =
    case m of
       Call -> s * exp (-q * t) * cdf normal d1 - k * exp (-r * t) * cdf normal d2
       Put  -> k * exp (-r * t) * cdf normal (-d2) - s * exp (-q * t) * cdf normal (-d1)
    where d1 = (log (s / k) + (r - q + (v * v) / 2) * t) / (v * sqrt t)
          d2 = d1 - v * sqrt t
          normal = Normal (0 :: Double) 1

instance BS (Option StockIndex) where
  price (Option (StockIndex s q) m European k t) r v =
    case m of
       Call -> s * exp (-q * t) * cdf normal d1 - k * exp (-r * t) * cdf normal d2
       Put  -> k * exp (-r * t) * cdf normal (-d2) - s * exp (-q * t) * cdf normal (-d1)
    where d1 = (log (s / k) + (r - q + (v * v) / 2) * t) / (v * sqrt t)
          d2 = d1 - v * sqrt t
          normal = Normal (0 :: Double) 1

-- Garman & Kohlhagen (1983)
--instance BS (Option Currency) where
--  price (Option (Currency s rf) m European k t) r v =
--    case m of
--       Call -> s * exp (-rf * t) * cdf normal d1 - k * exp (-r * t) * cdf normal d2
--       Put  -> k * exp (-r * t) * cdf normal (-d2) - s * exp (-rf * t) * cdf normal (-d1)
--    where d1 = (log (s / k) + (r - rf + (v * v) / 2) * t) / (v * sqrt t)
--          d2 = d1 - v * sqrt t
--          normal = Normal (0 :: Double) 1

-- cdi :: Double -> Double -> Double -> Double -> Double -> Double -> Double -> Double
-- cdi s k r v t q h = s * exp(-q * t) * (h / s)**(2 * lambda) * (cdf normal y)
--                     - k * exp(-r * t) * (h / s)**(2 * lambda - 2) * (cdf normal (y - (v * sqrt t)))
--   where lambda = ((r - q + (v * v) / 2)) / (v * v)
--         y = (log ( (h*h) / (s * k)) / v * sqrt t) + (lambda * v * sqrt t)
--         normal = Normal (0 :: Double) 1

-- cdo :: Double -> Double -> Double -> Double -> Double -> Double -> Double -> Double
-- cdo s k r v t q h = s * (cdf normal x1) * exp (-q * t) - k * exp(-r * t) * (cdf normal (x1 - v * sqrt t)) -
--                     s * exp(-q * t) * (h / s)**(2 * lambda) * (cdf normal y1) + k * exp(-r * t) * (h/ s)**(2 * lambda - 2) * (cdf normal (y1 - v * sqrt t))
--   where lambda = ((r - q + (v * v) / 2)) / (v * v)
--         y1 = (log (h / s) / v * sqrt t) + (lambda * v * sqrt t)
--         x1 = (log ( s / h ) / (v * sqrt t)) + (lambda * v * sqrt t)
--         normal = Normal (0 :: Double) 1

-- cui :: Double -> Double -> Double -> Double -> Double -> Double -> Double -> Double
-- cui s k r v t q h = s * (cdf normal x1) * exp(-q * t) - k * exp(-r * t) * (cdf normal (x1 - v * sqrt t)) -
--       s * exp(-q *t) * (h / s)**(2 * lambda) * ((cdf normal (-y)) - (cdf normal (-y1))) +
--       k * exp(-r * t) * (h / s)**(2 * lambda - 2) * ((cdf normal (-y + v * sqrt t)) - (cdf normal (-y1 + v * sqrt t)))
--   where lambda = ((r - q + (v * v) / 2)) / (v * v)
--         y1 = (log (h / s) / v * sqrt t) + (lambda * v * sqrt t)
--         x1 = (log ( s / h ) / (v * sqrt t)) + (lambda * v * sqrt t)
--         y = (log ( (h*h) / (s * k)) / v * sqrt t) + (lambda * v * sqrt t)
--         normal = Normal (0 :: Double) 1

-- pui :: Double -> Double -> Double -> Double -> Double -> Double -> Double -> Double
-- pui s k r v t q h = -s * exp(-q * t) * (h / s)**(2 * lambda) * (cdf normal (-y)) + k * exp(-r * t) * (h / s)**(2 * lambda - 2) * (cdf normal ((-y) + (v * sqrt t)))
--   where lambda = ((r - q + (v * v) / 2)) / (v * v)
--         y = (log ( (h*h) / (s * k)) / v * sqrt t) + (lambda * v * sqrt t)
--         normal = Normal (0 :: Double) 1

-- pdi :: Double -> Double -> Double -> Double -> Double -> Double -> Double -> Double
-- pdi s k r v t q h = -s * (cdf normal (-x1)) * exp(-q * t)
--                     + k * exp(-r * t) * (cdf normal (-x1 + v * sqrt t))
--                     + s * exp(-q * t) * (h / s)**(2 * lambda) * (cdf normal y - cdf normal y1)
--                     - k * exp(-r * t) * (h / s)**(2 * lambda -2) * (cdf normal (y - v * sqrt t) - cdf normal (y1 - v * sqrt t))
--   where lambda = ((r - q + (v * v) / 2)) / (v * v)
--         y1 = (log (h / s) / v * sqrt t) + (lambda * v * sqrt t)
--         x1 = (log ( s / h ) / (v * sqrt t)) + (lambda * v * sqrt t)
--         y = (log ( (h*h) / (s * k)) / v * sqrt t) + (lambda * v * sqrt t)
--         normal = Normal (0 :: Double) 1

-- puo s k r v t q h = -s * (cdf normal (-x1)) * exp(-q * t) + k * exp(-r * t) * (cdf normal (-x1 + v * sqrt t)) +
--                     s * exp(-q * t) * (h / s)**(2 * lambda) * (cdf normal (-y1)) - k * exp(-r * t) * (h / s)**(2 * lambda -2) * (cdf normal ((-y1) + v * sqrt t))
--   where lambda = ((r - q + (v * v) / 2)) / (v * v)
--         y1 = (log (h / s) / v * sqrt t) + (lambda * v * sqrt t)
--         x1 = (log ( s / h ) / (v * sqrt t)) + (lambda * v * sqrt t)
--         y = (log ( (h*h) / (s * k)) / v * sqrt t) + (lambda * v * sqrt t)
--         normal = Normal (0 :: Double) 1

-- instance BS (Option StockIndex) where
--   blackscholesprice (Option (StockIndex s) Call European k t) r v = call s k r v t q
--   blackscholesprice (Option (StockIndex s) Put European k t) r v = put s k r v t q

-- instance BS (Option Future) where
--   blackscholesprice (Option (Future f) Call European k t) r v = call f k r v t q
--   blackscholesprice (Option (Future f) Put European k t) r v = put f k r v t q

-- TODO : factor in rebates
-- instance BS (BarrierOption Equity) where
--   blackscholesprice (BarrierOption (Equity s q) x European y h k t) r v =
--     case (x, y) of
--       (DownAndIn, Call)
-- --        | s <= h -> Left "Contract is converted to a vanilla call"
--         | h <= k ->  cdi'
--         | h >= k ->   call' - cdo'
--       (DownAndOut, Call)
-- --        | s <= h -> Left "Contract is void"
--         | h <= k ->   call' - cdi'
--         | h >= k ->  cdo'
--       (UpAndIn, Call)
-- --        | h <= s -> Left "Contract is converted to a vanilla call"
--         | h <= k ->  call'
--         | h >= k ->   cui'
--       (UpAndOut, Call)
-- --        | s <= h -> Left "Contract is void"
--         | h <= k ->   call' - cdi'
--         | h >= k ->  cdo'
--       (DownAndIn, Put)
-- --        | s <= h -> Left "Contract is converted to a vanilla call"
--         | h <  k ->  pdi'
--         | h >= k ->  put'
--       (DownAndOut, Put)
-- --        | s <= h -> Left "Contract is void"
--         | h <  k ->   put' - pdi'
--         | h >=  k ->  0
--       (UpAndIn, Put)
-- --        | h <= s -> Left "Contract is converted to a vanilla call"
--         | h <= k ->   put' - puo'
--         | h >= k ->  pui'
--       (UpAndOut, Put)
-- --        | h <= s -> Left "Contract is void"
--         | h <= k ->   puo'
--         | h >= k ->   put' - pui'
--     where
--       call' = call s k r v t q
--       cdi'  = cdi  s k r v t q h
--       cdo'  = cdo  s k r v t q h
--       cui'  = cui  s k r v t q h
--       pdi'  = pdi  s k r v t q h
--       put'  = put  s k r v t q
--       pui'  = pui  s k r v t q h
--       puo'  = puo  s k r v t q h

max_stars_repo_name: co-category/hqfl

max_issues_repo_name: cokleisli/hqfl

max_forks_repo_name: cokleisli/hqfl


text:
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE ViewPatterns     #-}
{-# LANGUAGE TypeFamilies     #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE UnicodeSyntax #-}
--
-- |
-- Module      : Main
-- Description : Measures execution time of MC loop for 1D Heisenberg
-- Hamiltonian. Try @--help@ to get usage information.
-- Copyright   : (c) Tom Westerhout, 2018
-- License     : BSD3
-- Maintainer  : t.westerhout@student.ru.nl
-- Stability   : experimental
module Main where

import           Prelude                 hiding ( map
                                                , zipWithM
                                                )

import           Debug.Trace
import qualified System.Random.MWC             as MWC
-- import           Control.Lens hiding((<.>))
import           Control.Monad.Reader    hiding ( zipWithM )
import           Control.Monad.Primitive
import qualified Data.List                     as L
import           Data.Complex
import           Data.Semigroup
import           Data.Vector.Storable           ( Vector
                                                , (!)
                                                )
import qualified Data.Vector.Storable          as V
import qualified Data.Vector.Storable.Mutable  as MV
import           System.Exit
import           System.IO               hiding ( hGetLine )
import           System.Environment             ( getArgs
                                                , getProgName
                                                )
import           Foreign.Storable
import           Data.Text                      ( Text )
import           Data.Text.IO                   ( hGetLine )
import qualified Data.Text.IO                  as T

import           Lens.Micro
import           Lens.Micro.Extras

import           PSO.Random
import           PSO.FromPython
import           NQS.Rbm (Rbm, mkRbm)
import           NQS.Rbm.Mutable -- (sampleGradients)
import           NQS.Internal.Hamiltonian
import           NQS.Internal.Types -- (ℂ, ℝ)
import           NQS.Internal.Rbm (unsafeThawRbm)
import           NQS.Internal.LAPACK

fromPyFile :: FilePath -> IO Rbm
fromPyFile name = withFile name ReadMode toRbm
  where toRight :: Either String (Vector ℂ) -> Vector ℂ
        toRight (Right x) = x
        toRight (Left x)  = error x
        toRbm h = do
          hGetLine h
          a <- trace ("a...") $ toRight <$> readVector <$> hGetLine h
          b <- trace ("b...") $ toRight <$> readVector <$> hGetLine h
          s <- hGetLine h
          T.putStrLn s
          let !w = trace ("w...") $ toRight $ readMatrix s
          return $ mkRbm a b w

randomRbm :: Int -> Int -> (ℝ, ℝ) -> (ℝ, ℝ) -> (ℝ, ℝ) -> IO Rbm
randomRbm n m (lowV, highV) (lowH, highH) (lowW, highW) = do
  g <- mkMWCGen (Just 123)
  flip runReaderT g $ do
    visible <- uniformVector n (lowV :+ lowV, highV :+ highV)
    hidden  <- uniformVector m (lowH :+ lowH, highH :+ highH)
    weights <- uniformVector (n * m) (lowW :+ lowW, highW :+ highW)
    return $ trace ("mkRbm...") (mkRbm visible hidden weights)

numberSteps :: (Int, Int, Int) -> Int
numberSteps (low, high, step) = (high - low - 1) `div` step + 1

main = do -- NQS.Internal.LAPACK.test
  let filename = "/home/tom/src/tcm-swarm/cbits/test/input/rbm_6_6_0.in"
  rbm <- unsafeThawRbm =<< fromPyFile filename :: IO (MRbm (PrimState IO))
  -- rbm <- unsafeThawRbm =<< randomRbm 100 200 (-0.1, 0.1) (-0.1, 0.1) (-0.05, 0.05)
  let n = sizeVisible rbm
      config = defaultMCConfig & steps .~ (1000, 21000 * n, n)
                               & magnetisation .~ (Just 0)
      nParams = size rbm
      nSteps = numberSteps $ config ^. steps
  print n
  hamiltonian <- heisenberg (1, 1) (Just 5.0) (zip [0..] ([1..(n - 1)] ++ [0]))
  moments <- MV.new 4
  -- f <- newDenseVector nParams
  -- grad <- newDenseMatrix (nSteps * config ^. runs) nParams
  sampleMoments config hamiltonian rbm moments -- f grad
  e <- MV.read moments 0
  print ["Hello!", show e]

ext: hs

max_stars_repo_name: twesterhout/tcm-swarm

max_issues_repo_name: twesterhout/tcm-swarm

max_forks_repo_name: twesterhout/tcm-swarm


text:
module Data.Csv.HMatrix
    (
      -- * Usage example
      -- $example

      -- * Information
      -- $info

      decodeMatrix
    , decodeMatrixWith
    , encodeMatrix
    , encodeMatrixWith
    ) where

import Data.Csv
import Numeric.LinearAlgebra.HMatrix
import qualified Data.Vector as V
import Data.ByteString.Lazy (ByteString)
import qualified Data.ByteString.Char8 as C
import Data.Char (ord)

-- $example
--
-- For decoding we have to specify, if there is a header. This can be done with the HasHeader/NoHeader type:
--
-- > >>> decodeMatrix NoHeader "1.0,2.0,3.0\r\n4.0,5.0,6.0\r\n7.0,8.0,9.0\r\n"
-- > (3><3)
-- > [ 1.0, 2.0, 3.0
-- > , 4.0, 5.0, 6.0
-- > , 7.0, 8.0, 9.0 ]
--
-- Cassava, which is used for parsing the .csv files uses overloaded string literals. If you try this in ghci, make sure to start it with the right flag:
--
-- > >>> ghci -XOverloadedStrings
--
-- Encoding a file works pretty much the same, except that we do not need to specify a header.
--
-- > >>> encodeMatrix $ matrix 3 [1,2,3,4,5,6,7,8,9]
-- > "1.0,2.0,3.0\r\n4.0,5.0,6.0\r\n7.0,8.0,9.0\r\n"

-- $info
-- If you want to help improve this library, feel free to file an issue or send a pull-request on github. Every feedback is appreciated.
-- As of now only matrices of type Double are supported.

-- | Decodes a matrix.
decodeMatrix :: HasHeader       -- ^ From Data.Csv: specify if the CSV string has a header
             -> ByteString      -- ^ The 'ByteString' containing the CSVs
             -> Matrix Double   -- ^ The parsed 'Matrix'
decodeMatrix header = decodeMatrixWith header ','

-- | Decodes a matrix from ByteString and additionally allow
-- to specify the delimter which was used.
decodeMatrixWith :: HasHeader       -- ^ From Data.Csv: specify if the CSV string has a header
                 -> Char            -- ^ The delimiter
                 -> ByteString      -- ^ The 'ByteString' containing the CSVs
                 -> Matrix Double   -- ^ The parsed 'Matrix'
decodeMatrixWith header del s =
    case decodeWith opt header s of
        Left err -> error err
        Right v  -> fromLists . V.toList . V.map V.toList $ v
    where opt = defaultDecodeOptions { decDelimiter = fromIntegral (ord del) }

rowToRecord :: [Double] -> Record
rowToRecord x = record $ map (C.pack . show) x

-- | Encodes a matrix with comma as delimiter.
encodeMatrix :: Matrix Double   -- ^ The 'Matrix' to encode
             -> ByteString      -- ^ The resulting 'ByteString'
encodeMatrix = encodeMatrixWith ','

-- | Encodes a matrix but allows to specify a delimiter.
encodeMatrixWith :: Char            -- ^ The delimiter for separating the values
                 -> Matrix Double   -- ^ The 'Matrix' to encode
                 -> ByteString      -- ^ The resulting 'ByteString'
encodeMatrixWith del m = encodeWith opt s
    where opt = defaultEncodeOptions { encDelimiter = fromIntegral (ord del) }
          s = map rowToRecord $ toLists m

max_stars_repo_name: grtlr/hmatrix-csv

max_issues_repo_name: grtlr/hmatrix-csv

max_forks_repo_name: grtlr/hmatrix-csv


text:
{-# LANGUAGE BangPatterns     #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE GADTs            #-}
{-# LANGUAGE RecordWildCards  #-}
{-# LANGUAGE OverloadedStrings #-}

module Taiji.View.Commands.Rank
    ( plotRankParser
    , viewRanks
    ) where

import           Bio.Utils.Functions    (scale)
import qualified Data.Text as T
import qualified Data.Text.IO as T
import qualified Data.HashSet as S
import Data.Function (on)
import AI.Clustering.Hierarchical hiding (normalize)
import Control.Arrow (first)
import           Data.Colour            (blend)
import           Options.Applicative
import qualified Data.Matrix            as M
import qualified Data.Vector            as V
import           Diagrams.Backend.Cairo (B)
import           Diagrams.Backend.Cairo (renderCairo)
import           Diagrams.Prelude       hiding (option, scale, value, normalize)
import           Graphics.SVGFonts      (textSVG)
import Data.Colour.Palette.BrewerSet
import           Statistics.Sample
import Statistics.Correlation (spearman)
import           Text.Printf

import           Taiji.View.Types
import           Taiji.View.Utils hiding (zip, unzip)
import qualified Taiji.View.Utils as U

plotRankParser :: Parser Command
plotRankParser = ViewRanks
    <$> strArgument
      ( metavar "INPUT"
     <> help "rank file" )
    <*> strArgument
      ( metavar "OUTPUT"
     <> help "output file" )
    <*> (optional . strOption)
      ( long "expression" )
    <*> option auto
      ( long "cv"
     <> value 0.5
     <> help "TFs with coefficient of variance less than the specified value will be removed. (default: 1)"
      )
    <*> option auto
      ( long "min"
     <> value 1e-4
     <> help "lowerBound of TF rank." )
    <*> (optional . strOption) ( long "rowNamesFilter" )
    <*> (optional . strOption) ( long "output-values" )
    <*> (optional . option (maybeReader f)) ( long "rank-range" )
    <*> (optional . strOption) ( long "groups" )
  where
    f x = let [a,b] = T.splitOn "," $ T.pack x
          in Just (read $ T.unpack a, read $ T.unpack b)

viewRanks :: Command -> IO ()
viewRanks ViewRanks{..} = do
    grp <- case colGroup of
        Nothing -> return Nothing
        Just fl -> fmap (Just . map (T.splitOn ",") . T.lines) $ T.readFile fl

    rowFilt <- case rowNamesFilter of
        Nothing -> return $ const True
        Just fl -> do
            names <- T.lines <$> T.readFile fl
            return (`elem` names)

    df <- case exprFile of
        Just fl -> do
            mat <- cbind .
                map (reorderColumns (orderByCluster fst)) .
                groupDataFrame grp .
                reorderRows (orderByCluster fst) .
                uncurry U.zip . first normalize . U.unzip .
                filterRows (const $ filtCV cv . fst . V.unzip) .
                filterRows (const $ V.any ((>=minRank) . fst)) .
                filterRows (flip $ const rowFilt) <$> readData rankFile fl
            case outputValues of
                Nothing -> return ()
                Just x -> writeTable x (T.pack . show . fst) mat
            return $ Left mat
        Nothing -> do
            mat <- cbind .
                map (reorderColumns (orderByCluster id)) .
                groupDataFrame grp .
                reorderRows (orderByCluster id) .
                normalize .
                filterRows (const $ filtCV cv) .
                filterRows (const $ V.any (>=minRank)) .
                filterRows (flip $ const rowFilt) <$> readTable rankFile
            case outputValues of
                Nothing -> return ()
                Just x -> writeTable x (T.pack . show) mat
            return $ Right mat

    let w = width dia
        h = height dia
        n = fromIntegral (either (M.cols . _dataframe_data)
            (M.cols . _dataframe_data) df) * 50
        dia = bubblePlot df $ BubblePlotOpts 15 reds rankRange
    renderCairo output (dims2D n (n*(h/w))) dia

groupDataFrame :: Maybe [[T.Text]] -> DataFrame a -> [DataFrame a]
groupDataFrame Nothing df = [df]
groupDataFrame (Just grp) df = map (df `csub`) $ grp ++ [others]
  where
    others = S.toList $
        S.fromList (colNames df) `S.difference` S.fromList (concat grp)

orderByCluster :: (a -> Double) -> ReodrderFn a
orderByCluster f xs = flatten $ hclust Ward (V.fromList xs) dist
  where
    dist = euclidean `on` V.map f . snd

data BubblePlotOpts = BubblePlotOpts
    { _radius :: Double                      -- ^ The size of bubble
    , _colours :: V.Vector (Colour Double)   -- ^ Color scheme
    , _rank_range :: Maybe (Double, Double)        -- ^ The range of ranks
    }

bubblePlot :: Either (DataFrame (Double, Double)) (DataFrame Double)
           -> BubblePlotOpts -> Diagram B
bubblePlot (Left df) opts@BubblePlotOpts{..}
    | null (rowNames df) = error "Nothing to plot"
    | otherwise = center bubbles === strutY 30 === center legend
  where
    bubbles = drawBubbles df{_dataframe_data = M.zip ranks' expr'} opts
    legend = mkLegend opts
        (V.minimum $ M.flatten expr, V.maximum $ M.flatten expr)
        (case _rank_range of
            Nothing -> (V.minimum $ M.flatten ranks, V.maximum $ M.flatten ranks)
            Just r -> r )
    expr' = M.fromVector (M.dim expr) $ linearMap (4, _radius) $ M.flatten expr
    ranks' = M.fromVector (M.dim ranks) $ case _rank_range of
        Nothing -> linearMap (0, 1) $ M.flatten ranks
        Just rng -> linearMapBounded rng (0, 1) $ M.flatten ranks
    (ranks, expr) = M.unzip $ _dataframe_data df
bubblePlot (Right df) opts@BubblePlotOpts{..}
    | null (rowNames df) = error "Nothing to plot"
    | otherwise = center bubbles === strutY 30 === center legend
  where
    bubbles = drawBubbles df{_dataframe_data = M.zip expr' ranks'} opts
    expr' = M.fromVector (M.dim ranks) $ V.replicate (uncurry (*) $ M.dim ranks) 0.7
    ranks' = M.fromVector (M.dim ranks) $ linearMap (4, _radius) $ M.flatten ranks
    ranks = _dataframe_data df
    legend = mkCircleLegend "normalized rank score"
        ( V.minimum $ M.flatten $ _dataframe_data df
        , V.maximum $ M.flatten $ _dataframe_data df) _radius

drawBubbles :: DataFrame (Double, Double)
            -> BubblePlotOpts
            -> Diagram B
drawBubbles df BubblePlotOpts{..} = vsep 1 $ (colnames:) $
    zipWith drawBubble (rowNames df) $ M.toLists $ M.zip ranks' expr'
  where
    drawBubble lab xs = alignR $ textBounded lab ||| strutX 5 ||| bb ||| cor
      where
        bb = hsep 1 $ flip map xs $ \(x,y) -> withEnvelope unitEnvelope $
            circle y # lw 0 # fc (colorMapSmooth x _colours)
        cor = let r = spearman $ V.fromList xs
              in withEnvelope unitEnvelope $
                    square _radius # lw 0 # fc (colorMapSmooth ((r+1)/2) buYlRd)
    colnames = alignR $ hsep 1 $ map
        (\x -> (alignB $ textBounded x # rotate (90 @@ deg)) <> unitEnvelope) $
        colNames df ++ [""]
    (ranks', expr') = M.unzip $ _dataframe_data df
    unitEnvelope = circle _radius # lw 0

mkLegend :: BubblePlotOpts -> (Double, Double) -> (Double, Double) -> Diagram B
mkLegend BubblePlotOpts{..} (min_expr, max_expr) (min_rank, max_rank) =
    hsep 50 [rank_legend, expr_legend, cor_legend]
  where
    expr_legend = mkCircleLegend "inverse hyperbolic sine of expression level"
        (min_expr, max_expr) _radius
    rank_legend = mkColorLegend "normalized rank score" (min_rank, max_rank) _colours
    cor_legend = vsep 2
        [ rect 100 25 # lw 0 # fillTexture (mkGradient buYlRd)
        , center $ position $ zip [0^&0, 50^&0, 100^&0] $ map textBounded ["-1", "0", "1"]
        , textBounded "Spearman's correlation" ]
{-# INLINE mkLegend #-}

mkColorLegend :: T.Text   -- ^ Title
              -> (Double, Double)  -- ^ lower and upper bound
              -> V.Vector (Colour Double)   -- ^ color scheme
              -> Diagram B
mkColorLegend title (lo, hi) colors = vsep 2 $
    [ rect 100 25 # lw 0 # fillTexture (mkGradient colors)
    , center $ position $ zip [0^&0, 50^&0, 100^&0] $
        map (textBounded . T.pack . printf "%.2f") values
    , textBounded title ]
  where
    values = [lo, lo + (hi - lo) / 2, hi]

mkCircleLegend :: T.Text   -- ^ Title
               -> (Double, Double)     -- ^ Lower and upper bound
               -> Double    -- ^ radius
               -> Diagram B
mkCircleLegend title (lo, hi) r = center
    ( hsep 1 $ map (\(x,s) -> withEnvelope unitEnvelope
    (circle s # lw 0 # fc black) === strutY 2 === textBounded (T.pack $ printf "%.2f" x)) $
    zip values $ V.toList $ linearMap (4, r) $ V.fromList values ) ===
    strutY 2 === textBounded title
  where
    values = [lo, lo + (hi - lo) / 3 .. hi]
    unitEnvelope = circle r # lw 0 :: Diagram B

mkGradient :: V.Vector (Colour Double) -> Texture Double
mkGradient cs = mkLinearGradient stops ((-50) ^& 0) (50 ^& 0) GradPad
  where
    stops = mkStops $ zipWith (\c x -> (c, x, 1)) (V.toList cs) [0, 1/(n-1) .. 1]
    n = fromIntegral $ V.length cs
{-# INLINE mkGradient #-}

-- | Map numbers to colors
colorMapSmooth :: Double -- a value from 0 to 1
               -> V.Vector (Colour Double) -> Colour Double
colorMapSmooth x colors
    | isNaN x = black
    | x <0 || x > 1 = error "input value is out of range."
    | x == 1 = V.last colors
    | n == 2 = blend (1-p) (V.head colors) $ V.last colors
    | otherwise = blend (1 - p) (colors V.! i) $ colors V.! (i+1)
  where
    (i,p) = properFraction $ x * fromIntegral (n - 1)
    n = V.length colors
{-# INLINE colorMapSmooth #-}

textBounded :: T.Text -> Diagram B
textBounded x = stroke (textSVG (T.unpack x) 15) # lw 0 # fc black
{-# INLINE textBounded #-}

linearMap :: (Double, Double) -> V.Vector Double -> V.Vector Double
linearMap (lo, hi) xs = V.map f xs
  where
    f x = lo + (x - min') / (max' - min') * (hi - lo)
    min' = V.minimum xs
    max' = V.maximum xs
{-# INLINE linearMap #-}

linearMapBounded :: (Double, Double)    -- ^ Range of input data
                 -> (Double, Double)    -- ^ Range of output
                 -> V.Vector Double
                 -> V.Vector Double
linearMapBounded (min', max') (lo, hi) xs = V.map f xs
  where
    f x | x <= min' = lo
        | x >= max' = hi
        | otherwise = lo + (x - min') / (max' - min') * (hi - lo)
{-# INLINE linearMapBounded #-}

-------------------------------------------------------------------------------
-- Colours
-------------------------------------------------------------------------------

buYlRd :: V.Vector (Colour Double)
buYlRd = V.fromList $ reverse $ brewerSet RdYlBu 9

reds :: V.Vector (Colour Double)
reds = V.fromList [white, red]

-- | Normalize the data frame.
normalize :: DataFrame Double -> DataFrame Double
normalize = mapRows f
  where
    f xs | V.length xs <= 2 = V.map (logBase 2 . (/ V.head xs)) xs
         | otherwise = scale xs

-- | Determine whether the input pass the CV cutoff
filtCV :: Double -> V.Vector Double -> Bool
filtCV cutoff xs = sqrt v / m >= cutoff
  where
    (m, v) = meanVarianceUnb xs

-- | Determine whether the input pass the fold-change cutoff
filtFC :: Double -> V.Vector Double -> Bool
filtFC cutoff xs = V.maximum xs / V.minimum xs >= cutoff

max_stars_repo_name: Taiji-pipeline/Taiji-view

max_issues_repo_name: Taiji-pipeline/Taiji-view

max_forks_repo_name: Taiji-pipeline/Taiji-view


text:
module Environment where

import           Prelude
import           Control.Applicative
import           Control.Arrow
import qualified Control.Exception as E
import           Control.Monad
import           Control.Monad.Fix
import           Control.Monad.Reader
import           Control.Monad.State
import           Control.Monad.Trans
import           Control.Monad.Writer
import           Control.Parallel
import           Data.Bits
import qualified Data.ByteString.Lazy as B
import           Data.Char
import           Data.Complex
import           Data.Function
import           Data.Functor
import           Data.List
import           Data.Map (Map)
import qualified Data.Map as Map
import           Data.Maybe
import           Data.Monoid
import           Data.Ratio
import           Data.Set (Set)
import qualified Data.Set as Set
import           Data.Typeable
import           Data.Word
import           Numeric
import           System
import           System.IO
import           System.Process
import           System.Random
import           Text.Printf
import           Text.Regex
import           Text.Show.Functions

instance Typeable a => Show (IO a) where
  show io = '<' : show (typeOf io) ++ ">"

data PutStr = PutStr String

instance Show PutStr where
  show (PutStr s) = s

max_stars_repo_name: robotlolita/nanobot

max_issues_repo_name: GrimIRC/grim-bot

max_forks_repo_name: GrimIRC/grim-bot


text:
{-# LANGUAGE ScopedTypeVariables #-}

module Carrier (CarrierOpts(..), run, parseOpts) where

import           Data.ByteString.Builder
import           Data.Complex
import           IQ
import           Options.Applicative
import           Pipes
import qualified Pipes.Prelude as P
import           Render
import           System.IO
import Control.Monad (replicateM_)

data CarrierOpts = CarrierOpts
  { carrierLength :: Int
  , outputFile :: String
  , sampleRate :: Int }

parseOpts :: Parser CarrierOpts
parseOpts = CarrierOpts
         <$> option auto
             ( long "length"
            <> short 'l'
            <> metavar "LENGTH"
            <> help "Carrier length in seconds" )
         <*> strOption
             ( long "out-file"
            <> short 'o'
            <> metavar "OUTPUT_FILE"
            <> help "Output IQ file")
         <*> option auto
             ( long "sample-rate"
            <> short 's'
            <> value 1000000
            <> showDefault
            <> metavar "SAMPLE_RATE"
            <> help "IQ Sample rate in Hz")

synthesize :: Monad m => Int -> Int -> Producer IQ m ()
synthesize sr length = replicateM_ n $ yield sample
  where n = sr * length
        sample = 1.0 :+ 0.0

run :: CarrierOpts -> IO ()
run opts = do
  putStr "Synthesizing Samples ..."
  hFlush stdout

  withBinaryFile (outputFile opts) WriteMode  $ \f ->
    runEffect $ synthesize (sampleRate opts) (carrierLength opts)
            >-> cfileSink f

  putStrLn "Done."

max_stars_repo_name: hexagonal-sun/ayeQ

max_issues_repo_name: hexagonal-sun/ayeQ

max_forks_repo_name: hexagonal-sun/ayeQ


text:
{-# LANGUAGE BangPatterns, TypeFamilies #-}

module Fractal2D ( mandelbrot
                 , juliaAnimated
                 ) where

import Control.Loop
import Control.Monad
import Data.Complex
import Data.Word
import Data.Bits
import qualified Data.Vector.Storable.Mutable as VSM

import ConcurrentSegments

-- A few simple 2D fractals, just for testing

magnitudeSq :: RealFloat a => Complex a -> a
magnitudeSq c = realPart c * realPart c + imagPart c * imagPart c

-- http://linas.org/art-gallery/escape/escape.html
-- http://en.wikipedia.org/wiki/Mandelbrot_set#Continuous_.28smooth.29_coloring
fractionalIterCnt :: Int -> Complex Float -> Float
fractionalIterCnt iter escZ = max 0 $ fromIntegral iter - (log (log $ magnitudeSq escZ)) / log 2

-- Mandelbrot Set
--
-- http://en.wikipedia.org/wiki/Mandelbrot_set#Computer_drawings
mandelbrot :: Int -> Int -> VSM.IOVector Word32 -> Bool -> IO ()
mandelbrot w h fb smooth =
  forLoop 0 (< h) (+ 1) $ \py -> forLoop 0 (< w) (+ 1) $ \px ->
    let idx          = px + py * w
        fpx          = fromIntegral px :: Float
        fpy          = fromIntegral py :: Float
        fw           = fromIntegral w  :: Float
        fh           = fromIntegral h  :: Float
        ratio        = fw / fh
        y            = (fpy / fh) * 2 - 1 -- Y axis is [-1, +1]
        xshift       = (- 2) - ((2 * ratio - 2.5) * 0.5)
        x            = (fpx / fw) * 2 * ratio + xshift -- Keep aspect and center [-1, +0.5]
        c            = x :+ y
        maxIter      = 40
        (iCnt, escZ) = go (0 :: Int) (0 :+ 0)
        go iter z | (iter == maxIter) || -- Iteration limit?
                    magnitudeSq z > 4 * 4 = (iter, z) -- Hit escape radius?
                  | otherwise = let newZ = z * z + c
                                 in if   newZ == z -- Simple 1-cycle detection
                                    then (maxIter, z)
                                    else go (iter + 1) newZ
        icCont | iCnt == maxIter = fromIntegral maxIter -- Interior in case of limit
               | otherwise       = fractionalIterCnt iCnt escZ
        toGreen v = v `unsafeShiftL` 8
     in VSM.unsafeWrite fb idx . toGreen . truncate $
          if   smooth
          then icCont            / fromIntegral maxIter * 255 :: Float
          else fromIntegral iCnt / fromIntegral maxIter * 255 :: Float

-- Julia Set, computed in parallel
--
-- http://en.wikipedia.org/wiki/Julia_set
-- http://www.relativitybook.com/CoolStuff/julia_set.html
juliaAnimated :: Int -> Int -> VSM.IOVector Word32 -> Bool -> Double -> IO ()
juliaAnimated w h fb smooth tick =
  let !fTick       = realToFrac tick :: Float
      !scaledTick  = snd (properFraction $ fTick / 17 :: (Int, Float))
      !scaledTick2 = snd (properFraction $ fTick / 61 :: (Int, Float))
      !scaledTick3 = snd (properFraction $ fTick / 71 :: (Int, Float))
      !twoPi       = scaledTick * 2 * pi
      !juliaR      = sin twoPi * max 0.7 scaledTick2
      !juliaI      = cos twoPi * max 0.7 scaledTick3
      !fw          = fromIntegral w  :: Float
      !fh          = fromIntegral h  :: Float
      !ratio       = fw / fh
      !xshift      = 1.45 * ratio
      !maxIter     = 40
      doSeg lo hi  = forLoop lo (< hi) (+ 1) $ \py -> forLoop 0 (< w) (+ 1) $ \px ->
        let idx          = px + py * w
            fpx          = fromIntegral px :: Float
            fpy          = fromIntegral py :: Float
            y            = (fpy / fh) * 2.9 - 1.45 -- Y axis is [-1.45, +1.45]
            x            = (fpx / fw) * 2.9 * ratio - xshift -- Keep aspect and center
            c            = x :+ y
            (iCnt, escZ) = go (0 :: Int) c
            go iter z | (iter == maxIter) || -- Iteration limit?
                        magnitudeSq z > 4 * 4 = (iter, z) -- Hit escape radius?
                      | otherwise = let newZ = z * z + (juliaR :+ juliaI)
                                     in if   newZ == z -- Simple 1-cycle detection
                                        then (maxIter, z)
                                        else go (iter + 1) newZ
            icCont | iCnt == maxIter = fromIntegral maxIter -- Interior in case of limit
                   | otherwise       = fractionalIterCnt iCnt escZ
            toGreen v = v `unsafeShiftL` 8
         in VSM.unsafeWrite fb idx . toGreen . truncate $
              if   smooth
              then icCont            / fromIntegral maxIter * 255 :: Float
              else fromIntegral iCnt / fromIntegral maxIter * 255 :: Float
   in void $ forSegmentsConcurrently Nothing 0 h doSeg

max_stars_repo_name: blitzcode/ray-marching-distance-fields

max_issues_repo_name: blitzcode/ray-marching-distance-fields

max_forks_repo_name: blitzcode/ray-marching-distance-fields


text:
{-# LANGUAGE FlexibleContexts #-}

module Learning
    ( ols
    , gls
    , kmeans
    , Centroids
    , Clusters
    , Clustering
    , Initializer
    , average
    , invertMap
    ) where

import Debug.Trace
import Data.Set
import Data.Map.Strict
import Numeric.LinearAlgebra hiding ((!))

--covariance :: (Field n, Ord n, Num (Vector n))
--           => Set (Vector n)
--           -> Matrix n

--covariance set = exxt - uut where
--  exxt = mean $ Prelude.map (\a -> asColumn a <> asRow a) (Data.Set.elems set)
--  uut  = asColumn mu <> asRow mu
--  mu   = average set

-- | Ordinary least squares
ols :: (Field n, Monoid (Matrix n))
    => Matrix n
    -> Vector n
    -> Vector n

ols x y = head $ toColumns $ optimiseMult [l', tr l', tr x, asColumn y] where
  xTx = mTm x
  l   = chol xTx
  l'  = inv l

-- | Generalized least squares
gls :: (Eq n, Num (Vector n), Monoid (Matrix n), Field n)
    => Matrix n
    -> Matrix n
    -> Vector n
    -> Vector n

gls c x y = t where
  x' = tr x
  ic = inv c
  r  = inv (optimiseMult [x', ic, x])
  s  = optimiseMult [r, x, ic]
  t  = s #> y

type Centroids n = Set (Vector n)

type Clusters n = Map (Vector n) (Set (Vector n))

type Clustering n = (Centroids n, Clusters n)

type Initializer m n = Int -> Set (Vector n) -> m (Centroids n)

-- | KMeans clustering algorithm
kmeans :: (Functor m, Container Vector n, Ord n, Num (Vector n), Normed (Vector n), Field n)
       => Initializer m n
       -> Int
       -> Set (Vector n)
       -> m (Clustering n)

-- KMeans runs by first getting an original clustering, then reclustering
-- until a fixed point is reached.
kmeans init k set = (untilEqual fst (recluster k) . (\centroids -> (centroids, reassign centroids set))) <$> init k set

untilEqual :: Eq a => (e -> a) -> (e -> e) -> e -> e
untilEqual g f e = if g e' == g e then e else untilEqual g f e' where
  e' = f e

-- Reclustering means taking each cluster, computing a new centroid based on
-- the contents, and then reassigning each vector to the cluster which has a
-- centroid closest to it.
recluster :: (Container Vector n, Ord n, Fractional n, Normed (Vector n), Num (Vector n), Field n)
          => Int
          -> Clustering n
          -> Clustering n
recluster k (centroids, clusters) = (centroids', reassign centroids' (Data.Set.unions $ Data.Map.Strict.elems clusters)) where
  centroids' = Data.Set.map (\centroid -> average (clusters ! centroid)) centroids

average :: (Container c n, Fractional n, Ord n, Num (c n), Linear n c)
        => Set (c n) -> c n

average set = scale (1 / (fromInteger . toInteger $ Data.Set.size set)) $ Data.Set.foldr (+) one set' where
  (one, set') = Data.Set.deleteFindMin set

-- Reassignment goes like this: given a set of vectors and a set of centroids,
-- return a map from centroid to sets of vectors which have been found to be within this
-- centroid.
reassign :: (Normed (Vector n), Ord n, Eq n, Num (Vector n), Field n)
         => Set (Vector n)
         -> Set (Vector n)
         -> Clusters n

reassign centroids points = invertMap $ fromSet (getClosest centroids) points

-- | Given a map from x to y, returns a map from y to the preimage of y in the
--   original map
invertMap :: (Ord y, Ord x) => Map x y -> Map y (Set x)
invertMap map = st where
  st = fromListWith Data.Set.union [ (e, Data.Set.singleton k)
                                   | k <- ks, e <- es, map ! k == e ]
  es = Data.Map.Strict.elems map
  ks = keys map

getClosest :: (Normed (Vector n), Num (Vector n), Field n, Eq n) => Set (Vector n) -> Vector n -> Vector n
getClosest points point = if Data.Set.null points then error "getClosest" else fst $ Data.Set.foldr
  (\candidate (c, d) -> let d' = norm_2 (candidate - point)
                        in if d > d' then (candidate, d') else (c, d))
  (candidate, norm_2 (candidate - point))
   points' where
  (candidate, points') = Data.Set.deleteFindMin points

max_stars_repo_name: SamuelSchlesinger/learning

max_issues_repo_name: SamuelSchlesinger/learning

max_forks_repo_name: SamuelSchlesinger/learning


text:
{-|
  Copyright   :  (C) 2019, Google Inc.
  License     :  BSD2 (see the file LICENSE)
  Maintainer  :  Christiaan Baaij <christiaan.baaij@gmail.com>
-}
{-# LANGUAGE CPP             #-}
{-# LANGUAGE TemplateHaskell #-}

{-# OPTIONS_GHC -Wno-orphans #-}

module Clash.Class.AutoReg.Instances where

import           Clash.Class.AutoReg.Internal
import           Clash.CPP                           (maxTupleSize)

#if MIN_VERSION_base(4,12,0)
import           Data.Complex (Complex)
import           Data.Ord (Down)
#endif

import           Data.Ratio (Ratio)

#if MIN_VERSION_base(4,12,0)
deriveAutoReg ''Complex
deriveAutoReg ''Down
#endif

deriveAutoReg ''Ratio

-- TODO: this needs NFDataX instances to maxTupleSize that in turn need
-- TODO: generic instances for large tuples, but this slows GHC down by a
-- TODO: lot. See: https://phabricator.haskell.org/D2899. Maybe we could add
-- TODO: -flarge-tuples to clash-prelude?
-- deriveAutoRegTuples [2..maxTupleSize]
deriveAutoRegTuples [2..min maxTupleSize 15]

max_stars_repo_name: acairncross/clash-compiler

max_issues_repo_name: acairncross/clash-compiler

max_forks_repo_name: acairncross/clash-compiler


text:
-- {-# LANGUAGE NoImplicitPrelude #-}
{-# LANGUAGE DataKinds, GADTs, TypeFamilies #-}
{-# LANGUAGE ScopedTypeVariables  #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FunctionalDependencies #-}

module HBLAS.BLAS.Level2Spec(main, spec) where

import Data.Complex

import Test.Hspec

import Numerical.HBLAS.MatrixTypes as Matrix
import Numerical.HBLAS.BLAS.Level2 as BLAS
import qualified Data.Vector.Storable as S
import qualified Data.Vector.Storable.Mutable as SM

main :: IO ()
main = hspec spec

spec :: Spec
spec = do
  gbmvSpec
  gemvSpec
  gerSpec
  gercSpec
  geruSpec
  hbmvSpec
  hemvSpec
  herSpec
  her2Spec
  hpmvSpec
  hprSpec
  hpr2Spec
  sbmvSpec
  spmvSpec
  sprSpec
  spr2Spec
  symvSpec
  syrSpec
  syr2Spec
  tbmvSpec
  tbsvSpec
  tpmvSpec
  tpsvSpec
  trmvSpec
  trsvSpec
  gemvBugJune2017

gemvBugJune2017 :: Spec
gemvBugJune2017  =
  --context "?gemv bug"

 {-

 -}
    describe "dgemv abstraction" $ do
      it "3*2 matrix vector product sadness" $
             let m = S.fromList [1,2,3, 4,5,6]
                 x = S.fromList [10,20,30]
             in
             do

              m' <- S.unsafeThaw m
              x' <- S.unsafeThaw x
              y' <- SM.new 2

              let mat   = MutableDenseMatrix SRow 3 2 3 m'
                  xvec  = MutableDenseVector SDirect 3 1 x'
                  yvec  = MutableDenseVector SDirect 2 1 y'

              dgemv NoTranspose 1 0 mat xvec yvec
              resList <-mutableVectorToList $ _bufferMutDenseVector yvec
              resList `shouldBe` [140, 320]

gbmvSpec :: Spec
gbmvSpec =
  context "?GBMV" $ do
    describe "SGBMV" $ do
      it "3x5 a(5x5 matrix) all 1's" $ do
        matvecTest1SGBMV
    describe "DBGMV" $ do
      it "3x5 a(10x5 matrix) all 1's with beta 1.0" $ do
        matvecTest1DGBMV
    describe "CGBMV" $ do
      it "2x4 a(4x4 matrix) all 1+i's" $ do
        matvecTest1CGBMV
     -- it "gbmv on 2x4 a(4x4 matrix) all 1+i s with conjnotranspose" $ do
     --   matvecTest2CGBMV -- conjnotranspose is invalid
    describe "ZGBMV" $ do
      it "2x10 a(5x10 matrix) all 1+i s with transpose and alpha 1.0" $ do
        matvecTest1ZGBMV
      it "2x10 a(5x10 matrix) all 1+i s with conjtranspose and alpha 1.0" $ do
        matvecTest2ZGBMV

matvecTest1SGBMV :: IO ()
matvecTest1SGBMV = do
    a <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (3, 5) (\_ -> (1.0))
    x <- Matrix.generateMutableDenseVector 5 (\_ -> (1.0))
    res <- Matrix.generateMutableDenseVector 5 (\_ -> (1.0))
    BLAS.sgbmv Matrix.NoTranspose 5 5 1 1 1.0 a x 0.0 res
    resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector res
    resList `shouldBe` [2, 3, 3, 3, 2]

matvecTest1DGBMV :: IO ()
matvecTest1DGBMV = do
    a <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (3, 5) (\_ -> (1.0))
    x <- Matrix.generateMutableDenseVector 5 (\_ -> (1.0))
    res <- Matrix.generateMutableDenseVector 10 (\_ -> (1.0))
    BLAS.dgbmv Matrix.NoTranspose 10 5 1 1 1.0 a x 1.0 res
    resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector res
    resList `shouldBe` [3, 4, 4, 4, 3, 1, 1, 1, 1, 1]
    -- not [3, 4, 4, 4, 3, 2, 1, 1, 1, 1]

matvecTest1CGBMV :: IO ()
matvecTest1CGBMV = do
    a <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (2, 4) (\(x, y) -> [1.0:+1.0, 1.0:+1.0, 1.0:+1.0, 1.0:+1.0, 1.0:+1.0, 1.0:+1.0, 1.0:+1.0, 0] !! (x * 4 + y))
    x <- Matrix.generateMutableDenseVector 4 (\_ -> (1.0:+1.0))
    res <- Matrix.generateMutableDenseVector 4 (\_ -> (1.0:+1.0))
    BLAS.cgbmv Matrix.NoTranspose 4 4 0 1 1.0 a x 0.0 res
    resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector res
    resList `shouldBe` [0:+4.0, 0:+4.0, 0:+4.0, 0:+2.0]

matvecTest2CGBMV :: IO ()
matvecTest2CGBMV = do
    a <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (2, 4) (\(x, y) -> [1.0:+1.0, 1.0:+1.0, 1.0:+1.0, 1.0:+1.0, 1.0:+1.0, 1.0:+1.0, 1.0:+1.0, 0] !! (x * 4 + y))
    x <- Matrix.generateMutableDenseVector 4 (\_ -> (1.0:+1.0))
    res <- Matrix.generateMutableDenseVector 4 (\_ -> (1.0:+1.0))
    BLAS.cgbmv Matrix.ConjNoTranspose 4 4 0 1 1.0 a x 0.0 res
    resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector res
    resList `shouldBe` [4.0:+0, 4.0:+0, 4.0:+0, 2.0:+0]

matvecTest1ZGBMV :: IO ()
matvecTest1ZGBMV = do
    a <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (2, 10) (\_ -> (1.0:+1.0))
    x <- Matrix.generateMutableDenseVector 10 (\_ -> (1.0:+1.0))
    res <- Matrix.generateMutableDenseVector 5 (\_ -> (1.0:+1.0))
    BLAS.cgbmv Matrix.Transpose 10 5 1 0 2.0 a x 0.0 res
    resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector res
    resList `shouldBe` [0:+8.0, 0:+8.0, 0:+8.0, 0:+8.0, 0:+8.0]
    -- not [0:+8.0, 0:+8.0, 0:+8.0, 0:+8.0, 0:+4.0]

matvecTest2ZGBMV :: IO ()
matvecTest2ZGBMV = do
    a <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (2, 10) (\_ -> (1.0:+1.0))
    x <- Matrix.generateMutableDenseVector 10 (\_ -> (1.0:+1.0))
    res <- Matrix.generateMutableDenseVector 5 (\_ -> (1.0:+1.0))
    BLAS.cgbmv Matrix.ConjTranspose 10 5 1 0 2.0 a x 0.0 res
    resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector res
    resList `shouldBe` [8.0:+0, 8.0:+0, 8.0:+0, 8.0:+0, 8.0:+0]
    -- not [0:+8.0, 0:+8.0, 0:+8.0, 0:+8.0, 0:+4.0]

-- notes for gbmv
-- column > rows: the tail of the supper diagonals is considered.
-- rows > column: the tail of the sub diagonals is not considered.

gemvSpec :: Spec
gemvSpec =
  context "?GEMV" $ do
    describe "SGEMV" $ do
      it "2x2 all 1's" $ do
        matmatTest1SGEMV
    describe "DGEMV" $ do
      it "2x2 all 1's" $ do
        matmatTest1DGEMV
    describe "CGEMV" $ do
      it "2x2 all 1's" $ do
        matmatTest1CGEMV
    describe "ZGEMV" $ do
      it "2x2 all 1's" $ do
        matmatTest1ZGEMV

matmatTest1SGEMV :: IO ()
matmatTest1SGEMV = do
    left  <- Matrix.generateMutableDenseMatrix (Matrix.SRow)  (2,2) (\_ -> (1.0::Float))
    right <- Matrix.generateMutableDenseVector 2 (\_ -> (1.0 :: Float))
    res  <- Matrix.generateMutableDenseVector  2 (\_ -> (0.0 :: Float))
    BLAS.sgemv Matrix.NoTranspose  1.0 1.0 left right res
    resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector res
    resList `shouldBe` [2,2]

matmatTest1DGEMV :: IO ()
matmatTest1DGEMV = do
    left  <- Matrix.generateMutableDenseMatrix (Matrix.SRow)  (2,2) (\_ -> (1.0))
    right <- Matrix.generateMutableDenseVector 2 (\_ -> (1.0 ))
    res  <- Matrix.generateMutableDenseVector 2  (\_ -> (0.0 ))
    BLAS.dgemv Matrix.NoTranspose 1.0 1.0 left right res
    resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector res
    resList `shouldBe` [2.0,2.0]

matmatTest1CGEMV :: IO ()
matmatTest1CGEMV = do
    left  <- Matrix.generateMutableDenseMatrix (Matrix.SRow)  (2,2) (\_ -> (1.0))
    right <- Matrix.generateMutableDenseVector 2 (\_ -> (1.0 ))
    res  <- Matrix.generateMutableDenseVector  2 (\_ -> (0.0 ))
    BLAS.cgemv Matrix.NoTranspose  1.0 1.0 left right res
    resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector res
    resList `shouldBe` [2.0,2.0]

matmatTest1ZGEMV :: IO ()
matmatTest1ZGEMV = do
    left  <- Matrix.generateMutableDenseMatrix (Matrix.SRow)  (2,2) (\_ -> (1.0))
    right <- Matrix.generateMutableDenseVector 2 (\_ -> (1.0 ))
    res  <- Matrix.generateMutableDenseVector 2 (\_ -> (0.0 ))
    BLAS.zgemv Matrix.NoTranspose  1.0 1.0 left right res
    resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector res
    resList `shouldBe` [2.0,2.0]

----
----
gerSpec :: Spec
gerSpec =
  context "?GER" $ do
    describe "SGER" $ do
      it "2x2 all 1's" $ do
        matmatTest1SGER
    describe "DGER" $ do
      it "2x2 all 1's" $ do
        matmatTest1DGER

matmatTest1SGER :: IO ()
matmatTest1SGER = do
  res <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (2,2) (\_ -> 1.0)
  x <- Matrix.generateMutableDenseVector 2 (\_ -> 2.0)
  y <- Matrix.generateMutableDenseVector 2 (\_ -> 3.0)
  BLAS.sger 2.0 x y res
  resList <- Matrix.mutableVectorToList $ _bufferDenMutMat res
  resList `shouldBe` [13.0,13.0,13.0,13.0]

matmatTest1DGER :: IO ()
matmatTest1DGER = do
  res <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (2,2) (\_ -> 1.0)
  x <- Matrix.generateMutableDenseVector 2 (\_ -> 2.0)
  y <- Matrix.generateMutableDenseVector 2 (\_ -> 3.0)
  BLAS.sger 2.0 x y res
  resList <- Matrix.mutableVectorToList $ _bufferDenMutMat res
  resList `shouldBe` [13.0,13.0,13.0,13.0]

gercSpec :: Spec
gercSpec =
  context "?GERC" $ do
    describe "CGERC" $ do
      it "2x3 all 1+i's" $ do
        matmatTest1CGERC
    describe "ZGERC" $ do
      it "2x3 all 1+i's" $ do
        matmatTest1ZGERC

matmatTest1CGERC :: IO ()
matmatTest1CGERC = do
  res <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (3,2) (\_ -> 1.0:+1.0)
  x <- Matrix.generateMutableDenseVector 3 (\_ -> 2.0:+3.0)
  y <- Matrix.generateMutableDenseVector 2 (\_ -> 3.0:+2.0)
  BLAS.cgerc 2.0 x y res
  resList <- Matrix.mutableVectorToList $ _bufferDenMutMat res
  resList `shouldBe` [25.0:+11.0, 25.0:+11.0, 1.0:+1.0, 25.0:+11.0, 25.0:+11.0, 1.0:+1.0]
  -- why the following is not correct...
  -- resList `shouldBe` [25.0:+11.0, 25.0:+11.0, 25.0:+11.0, 25.0:+11.0, 25.0:+11.0, 25.0:+11.0]

matmatTest1ZGERC :: IO ()
matmatTest1ZGERC = do
  res <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (3,2) (\_ -> 1.0:+1.0)
  x <- Matrix.generateMutableDenseVector 3 (\_ -> 2.0:+3.0)
  y <- Matrix.generateMutableDenseVector 2 (\_ -> 3.0:+2.0)
  BLAS.zgerc 2.0 x y res
  resList <- Matrix.mutableVectorToList $ _bufferDenMutMat res
  resList `shouldBe` [25.0:+11.0, 25.0:+11.0, 1.0:+1.0, 25.0:+11.0, 25.0:+11.0, 1.0:+1.0]
  -- why the following is not correct...
  -- resList `shouldBe` [25.0:+11.0, 25.0:+11.0, 25.0:+11.0, 25.0:+11.0, 25.0:+11.0, 25.0:+11.0]

geruSpec :: Spec
geruSpec =
  context "?GERU" $ do
    describe "CGERU" $ do
      it "2x3 all 1+i's" $ do
        matmatTest1CGERU
    describe "ZGERU" $ do
      it "2x3 all 1+i's" $ do
        matmatTest1ZGERU

matmatTest1CGERU :: IO ()
matmatTest1CGERU = do
  res <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (3,2) (\_ -> 1.0:+1.0)
  x <- Matrix.generateMutableDenseVector 3 (\_ -> 2.0:+(-3.0))
  y <- Matrix.generateMutableDenseVector 2 (\_ -> 3.0:+(-2.0))
  BLAS.cgeru 2.0 x y res
  resList <- Matrix.mutableVectorToList $ _bufferDenMutMat res
  resList `shouldBe` [1.0:+(-25.0), 1.0:+(-25.0), 1.0:+1.0, 1.0:+(-25.0), 1.0:+(-25.0), 1.0:+1.0]
  -- why the following is not correct...
  -- resList `shouldBe` [1.0:+(-25.0), 1.0:+(-25.0), 1.0:+(-25.0), 1.0:+(-25.0), 1.0:+(-25.0), 1.0:+(-25.0)]

matmatTest1ZGERU :: IO ()
matmatTest1ZGERU = do
  res <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (3,2) (\_ -> 1.0:+1.0)
  x <- Matrix.generateMutableDenseVector 3 (\_ -> 2.0:+(-3.0))
  y <- Matrix.generateMutableDenseVector 2 (\_ -> 3.0:+(-2.0))
  BLAS.zgeru 2.0 x y res
  resList <- Matrix.mutableVectorToList $ _bufferDenMutMat res
  resList `shouldBe` [1.0:+(-25.0), 1.0:+(-25.0), 1.0:+1.0, 1.0:+(-25.0), 1.0:+(-25.0), 1.0:+1.0]
  -- why the following is not correct...
  -- resList `shouldBe` [1.0:+(-25.0), 1.0:+(-25.0), 1.0:+(-25.0), 1.0:+(-25.0), 1.0:+(-25.0), 1.0:+(-25.0)]

-- [1:+0    1:+1    1:+1    0:+0]
-- [1:+(-1) 1:+0    1:+1    1:+1]
-- [1:+(-1) 1:+(-1) 1:+0    1:+1]
-- [0:+0    1:+(-1) 1:+(-1) 1:+0]
--
-- [1:+1]
-- [1:+1]
-- [1:+1]
-- [1:+1]

hbmvSpec :: Spec
hbmvSpec =
  context "?HBMV" $ do
    describe "CHBMV" $ do
      it "4*3 a(4x4 matrix) upper all 1+i's" $ do
        matvecTest1CHBMV
      it "4*3 a(4x4 matrix) upper all 1+i's (column oriented)" $ do
        matvecTest2CHBMV
    describe "ZHBMV" $ do
      it "4*3 a(4x4 matrix) lower all 1+i's" $ do
        matvecTest1ZHBMV

matvecTest1CHBMV :: IO ()
matvecTest1CHBMV = do
  a <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (4, 3) (\_ -> 1.0:+1.0)
  x <- Matrix.generateMutableDenseVector 4 (\_ -> 1.0:+1.0)
  y <- Matrix.generateMutableDenseVector 4 (\_ -> 0.0:+0.0)
  BLAS.chbmv Matrix.MatUpper 2 1.0 a x 0.0 y
  resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector y
  resList `shouldBe` [1.0:+5.0, 3.0:+5.0, 5.0:+3.0, 5.0:+1.0]

matvecTest2CHBMV :: IO ()
matvecTest2CHBMV = do
  a <- Matrix.generateMutableDenseMatrix (Matrix.SColumn) (4, 3) (\_ -> 1.0:+1.0)
  x <- Matrix.generateMutableDenseVector 4 (\_ -> 1.0:+1.0)
  y <- Matrix.generateMutableDenseVector 4 (\_ -> 0.0:+0.0)
  BLAS.chbmv Matrix.MatUpper 2 1.0 a x 0.0 y
  resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector y
  resList `shouldBe` [1.0:+5.0, 3.0:+5.0, 5.0:+3.0, 5.0:+1.0]

-- [1:+0    1:+1    1:+1    0:+0]
-- [1:+(-1) 1:+0    1:+1    1:+1]
-- [1:+(-1) 1:+(-1) 1:+0    1:+1]
-- [0:+0    1:+(-1) 1:+(-1) 1:+0]
--
-- [1:+1]
-- [1:+1]
-- [1:+1]
-- [1:+1]
matvecTest1ZHBMV :: IO ()
matvecTest1ZHBMV = do
  a <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (4, 3) (\_ -> 1.0:+(-1.0))
  x <- Matrix.generateMutableDenseVector 4 (\_ -> 1.0:+1.0)
  y <- Matrix.generateMutableDenseVector 4 (\_ -> 0.0:+0.0)
  BLAS.zhbmv Matrix.MatLower 2 1.0 a x 0.0 y
  resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector y
  resList `shouldBe` [1.0:+5.0, 3.0:+5.0, 5.0:+3.0, 5.0:+1.0]

-- [1:+0    1:+1    1:+1    1:+1]
-- [1:+(-1) 1:+0    1:+1    1:+1]
-- [1:+(-1) 1:+(-1) 1:+0    1:+1]
-- [1:+(-1) 1:+(-1) 1:+(-1) 1:+0]
--
-- [1:+1]
-- [1:+1]
-- [1:+1]
-- [1:+1]
--
hemvSpec :: Spec
hemvSpec =
  context "?HEMV" $ do
    describe "CHEMV" $ do
      it "4*3 a(4x4 matrix) upper all 1+i's" $ do
        matvecTest1CHEMV
      it "4*3 a(4x4 matrix) upper all 1+i s (column oriented)" $ do
        matvecTest2CHEMV
    describe "ZHEMV" $ do
      it "4*3 a(4x4 matrix) lower all 1+i's" $ do
        matvecTest1ZHEMV

matvecTest1CHEMV :: IO ()
matvecTest1CHEMV = do
  a <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (4, 4) (\_ -> 1.0:+1.0)
  x <- Matrix.generateMutableDenseVector 4 (\_ -> 1.0:+1.0)
  y <- Matrix.generateMutableDenseVector 4 (\_ -> 0.0:+0.0)
  BLAS.chemv Matrix.MatUpper 1.0 a x 0.0 y
  resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector y
  resList `shouldBe` [1.0:+7.0, 3.0:+5.0, 5.0:+3.0, 7.0:+1.0]

matvecTest2CHEMV :: IO ()
matvecTest2CHEMV = do
  a <- Matrix.generateMutableDenseMatrix (Matrix.SColumn) (4, 4) (\_ -> 1.0:+1.0)
  x <- Matrix.generateMutableDenseVector 4 (\_ -> 1.0:+1.0)
  y <- Matrix.generateMutableDenseVector 4 (\_ -> 0.0:+0.0)
  BLAS.chemv Matrix.MatUpper 1.0 a x 0.0 y
  resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector y
  resList `shouldBe` [1.0:+7.0, 3.0:+5.0, 5.0:+3.0, 7.0:+1.0]

matvecTest1ZHEMV :: IO ()
matvecTest1ZHEMV = do
  a <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (4, 4) (\_ -> 1.0:+(-1.0))
  x <- Matrix.generateMutableDenseVector 4 (\_ -> 1.0:+1.0)
  y <- Matrix.generateMutableDenseVector 4 (\_ -> 0.0:+0.0)
  BLAS.zhemv Matrix.MatLower 1.0 a x 0.0 y
  resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector y
  resList `shouldBe` [1.0:+7.0, 3.0:+5.0, 5.0:+3.0, 7.0:+1.0]

herSpec :: Spec
herSpec =
  context "?HER" $ do
    describe "CHER" $ do
      it "4*4 a upper all 1+i's" $ do
        matvecTest1CHER
    describe "ZHER" $ do
      it "4*4 a upper all 1+i's" $ do
        matvecTest1ZHER

matvecTest1CHER :: IO ()
matvecTest1CHER = do
  a <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (4, 4) (\_ -> 1.0:+1.0)
  x <- Matrix.generateMutableDenseVector 4 (\idx -> [1.0:+1.0, 2.0:+2.0, 3.0:+3.0, 4.0:+4.0] !! idx)
  BLAS.cher Matrix.MatUpper 1.0 x a
  resList <- Matrix.mutableVectorToList $ _bufferDenMutMat a
  resList `shouldBe` [3.0:+0.0, 5.0:+1.0,  7.0:+1.0,  9.0:+1.0,
               1.0:+1.0, 9.0:+0.0, 13.0:+1.0, 17.0:+1.0,
               1.0:+1.0, 1.0:+1.0, 19.0:+0.0, 25.0:+1.0,
               1.0:+1.0, 1.0:+1.0,  1.0:+1.0, 33.0:+0.0]

matvecTest1ZHER :: IO ()
matvecTest1ZHER = do
  a <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (4, 4) (\_ -> 1.0:+1.0)
  x <- Matrix.generateMutableDenseVector 4 (\idx -> [1.0:+1.0, 2.0:+2.0, 3.0:+3.0, 4.0:+4.0] !! idx)
  BLAS.zher Matrix.MatLower 1.0 x a
  resList <- Matrix.mutableVectorToList $ _bufferDenMutMat a
  resList `shouldBe` [3.0:+0.0,  1.0:+1.0,  1.0:+1.0,  1.0:+1.0,
               5.0:+1.0,  9.0:+0.0,  1.0:+1.0,  1.0:+1.0,
               7.0:+1.0, 13.0:+1.0, 19.0:+0.0,  1.0:+1.0,
               9.0:+1.0, 17.0:+1.0, 25.0:+1.0, 33.0:+0.0]

her2Spec :: Spec
her2Spec =
  context "?HER2" $ do
    describe "CHER2" $ do
      it "4*4 a upper all 1+i's" $ do
        matvecTest1CHER2
    describe "ZHER2" $ do
      it "4*4 a upper all 1+i's" $ do
        matvecTest1ZHER2

matvecTest1CHER2 :: IO ()
matvecTest1CHER2 = do
  a <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (4, 4) (\_ -> 1.0:+1.0)
  x <- Matrix.generateMutableDenseVector 4 (\idx -> [1.0:+1.0, 2.0:+2.0, 3.0:+3.0, 4.0:+4.0] !! idx)
  y <- Matrix.generateMutableDenseVector 4 (\idx -> [1.0:+1.0, 2.0:+2.0, 3.0:+3.0, 4.0:+4.0] !! idx)
  BLAS.cher2 Matrix.MatUpper 1.0 x y a
  resList <- Matrix.mutableVectorToList $ _bufferDenMutMat a
  resList `shouldBe` [5.0:+0.0,  9.0:+1.0, 13.0:+1.0, 17.0:+1.0,
               1.0:+1.0, 17.0:+0.0, 25.0:+1.0, 33.0:+1.0,
               1.0:+1.0,  1.0:+1.0, 37.0:+0.0, 49.0:+1.0,
               1.0:+1.0,  1.0:+1.0,  1.0:+1.0, 65.0:+0.0]

matvecTest1ZHER2 :: IO ()
matvecTest1ZHER2 = do
  a <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (4, 4) (\_ -> 1.0:+1.0)
  x <- Matrix.generateMutableDenseVector 4 (\idx -> [1.0:+1.0, 2.0:+2.0, 3.0:+3.0, 4.0:+4.0] !! idx)
  y <- Matrix.generateMutableDenseVector 4 (\idx -> [1.0:+1.0, 2.0:+2.0, 3.0:+3.0, 4.0:+4.0] !! idx)
  BLAS.zher2 Matrix.MatLower 1.0 x y a
  resList <- Matrix.mutableVectorToList $ _bufferDenMutMat a
  resList `shouldBe` [ 5.0:+0.0,  1.0:+1.0,  1.0:+1.0,  1.0:+1.0,
                9.0:+1.0, 17.0:+0.0,  1.0:+1.0,  1.0:+1.0,
               13.0:+1.0, 25.0:+1.0, 37.0:+0.0,  1.0:+1.0,
               17.0:+1.0, 33.0:+1.0, 49.0:+1.0, 65.0:+0.0]

hpmvSpec :: Spec
hpmvSpec =
  context "?HPMV" $ do
    describe "CHPMV" $ do
      it "4*4 a upper (row oriented)" $ do
        matvecTest1CHPMV
    describe "ZHPMV" $ do
      it "4*4 a upper (column oriented)" $ do
        matvecTest1ZHPMV

-- [0:+0    1:+1    2:+2    3:+3]
-- [1:+(-1) 4:+0    5:+5    6:+6]
-- [2:+(-2) 5:+(-5) 7:+0    8:+8]
-- [3:+(-3) 6:+(-6) 8:+(-8) 9:+0]
--
-- [2:+2]
-- [2:+2]
-- [2:+2]
-- [2:+2]
--
-- [ 3:+27]
-- [15:+55]
-- [45:+49]
-- [89:+21]
matvecTest1CHPMV :: IO ()
matvecTest1CHPMV = do
  a <- Matrix.generateMutableDenseVector 10 (\idx -> [0.0:+0.0, 1.0:+1.0, 2.0:+2.0, 3.0:+3.0, 4.0:+4.0, 5.0:+5.0, 6.0:+6.0, 7.0:+7.0, 8.0:+8.0, 9.0:+9.0] !! idx)
  x <- Matrix.generateMutableDenseVector 4 (\idx -> [2.0:+2.0, 2.0:+2.0, 2.0:+2.0, 2.0:+2.0] !! idx)
  y <- Matrix.generateMutableDenseVector 4 (\idx -> [3.0:+3.0, 3.0:+3.0, 3.0:+3.0, 3.0:+3.0] !! idx)
  BLAS.chpmv Matrix.SRow Matrix.MatUpper 4 1.0 a x 1.0 y
  resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector y
  resList `shouldBe` [ 3.0:+27.0, 15.0:+55.0, 45.0:+49.0, 89.0:+21.0]

-- [0:+0    1:+1    3:+3    6:+6]
-- [1:+(-1) 2:+0    4:+4    7:+7]
-- [3:+(-3) 4:+(-4) 5:+0    8:+8]
-- [6:+(-6) 7:+(-7) 8:+(-8) 9:+0]
--
-- [2:+2]
-- [2:+2]
-- [2:+2]
-- [2:+2]
--
-- [  6:+46]
-- [ 14:+54]
-- [ 44:+48]
-- [108:+24]
matvecTest1ZHPMV :: IO ()
matvecTest1ZHPMV = do
  a <- Matrix.generateMutableDenseVector 10 (\idx -> [0.0:+0.0, 1.0:+1.0, 2.0:+2.0, 3.0:+3.0, 4.0:+4.0, 5.0:+5.0, 6.0:+6.0, 7.0:+7.0, 8.0:+8.0, 9.0:+9.0] !! idx)
  x <- Matrix.generateMutableDenseVector 4 (\idx -> [2.0:+2.0, 2.0:+2.0, 2.0:+2.0, 2.0:+2.0] !! idx)
  y <- Matrix.generateMutableDenseVector 4 (\idx -> [3.0:+3.0, 3.0:+3.0, 3.0:+3.0, 3.0:+3.0] !! idx)
  BLAS.zhpmv Matrix.SColumn Matrix.MatUpper 4 1.0 a x 2.0 y
  resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector y
  resList `shouldBe` [ 6.0:+46.0, 14.0:+54.0, 44.0:+48.0, 108.0:+24.0]

-- [1:+1]
-- [2:+2]
-- [3:+3]
-- [4:+4]
--
-- [2:+0    4:+0    6:+0    8:+0 ]
-- [4:+0    8:+0    12:+0   16:+0]
-- [6:+0    12:+0   18:+0   24:+0]
-- [8:+0    16:+0   24:+0   32:+0]
--
-- [0:+0    1:+1    3:+3    6:+6]
-- [1:+(-1) 2:+0    4:+4    7:+7]
-- [3:+(-3) 4:+(-4) 5:+0    8:+8]
-- [6:+(-6) 7:+(-7) 8:+(-8) 9:+0]
--
-- [2:+0     5:+1     9:+3     14:+6]
-- [5:+(-1)  10:+0    16:+4    23:+7]
-- [9:+(-3)  16:+(-4) 23:+0    32:+8]
-- [14:+(-6) 23:+(-7) 32:+(-8) 41:+0]

hprSpec :: Spec
hprSpec =
  context "?HPR" $ do
    describe "CHPR" $ do
      it "4*4 a upper (column oriented)" $ do
        matvecTest1CHPR
    describe "ZHPR" $ do
      it "4*4 a lower (row oriented)" $ do
        matvecTest1ZHPR

matvecTest1CHPR :: IO ()
matvecTest1CHPR = do
  a <- Matrix.generateMutableDenseVector 10 (\idx -> [0.0:+0.0, 1.0:+1.0, 2.0:+2.0, 3.0:+3.0, 4.0:+4.0, 5.0:+5.0, 6.0:+6.0, 7.0:+7.0, 8.0:+8.0, 9.0:+9.0] !! idx)
  x <- Matrix.generateMutableDenseVector 4 (\idx -> [1.0:+1.0, 2.0:+2.0, 3.0:+3.0, 4.0:+4.0] !! idx)
  BLAS.chpr Matrix.SColumn Matrix.MatUpper 4 1.0 x a
  resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector a
  resList `shouldBe` [2.0:+0.0, 5.0:+1.0, 10.0:+0.0, 9.0:+3.0, 16.0:+4.0, 23.0:+0.0, 14.0:+6.0, 23.0:+7.0,  32.0:+8.0, 41.0:+0.0]

-- [1:+1]
-- [2:+2]
-- [3:+3]
-- [4:+4]
--
-- [2:+0    4:+0    6:+0    8:+0 ]
-- [4:+0    8:+0    12:+0   16:+0]
-- [6:+0    12:+0   18:+0   24:+0]
-- [8:+0    16:+0   24:+0   32:+0]
--
-- [0:+0    1:+1    3:+3    6:+6]
-- [1:+(-1) 2:+0    4:+4    7:+7]
-- [3:+(-3) 4:+(-4) 5:+0    8:+8]
-- [6:+(-6) 7:+(-7) 8:+(-8) 9:+0]
--
-- [2:+0     5:+1     9:+3     14:+6]
-- [5:+(-1)  10:+0    16:+4    23:+7]
-- [9:+(-3)  16:+(-4) 23:+0    32:+8]
-- [14:+(-6) 23:+(-7) 32:+(-8) 41:+0]
matvecTest1ZHPR :: IO ()
matvecTest1ZHPR = do
  a <- Matrix.generateMutableDenseVector 10 (\idx -> [0.0:+(-0.0), 1.0:+(-1.0), 2.0:+(-2.0), 3.0:+(-3.0), 4.0:+(-4.0), 5.0:+(-5.0), 6.0:+(-6.0), 7.0:+(-7.0), 8.0:+(-8.0), 9.0:+(-9.0)] !! idx)
  x <- Matrix.generateMutableDenseVector 4 (\idx -> [1.0:+1.0, 2.0:+2.0, 3.0:+3.0, 4.0:+4.0] !! idx)
  BLAS.zhpr Matrix.SRow Matrix.MatLower 4 1.0 x a
  resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector a
  resList `shouldBe` [2.0:+0.0, 5.0:+(-1.0), 10.0:+0.0, 9.0:+(-3.0), 16.0:+(-4.0), 23.0:+0.0, 14.0:+(-6.0), 23.0:+(-7.0), 32.0:+(-8.0), 41.0:+0.0]

hpr2Spec :: Spec
hpr2Spec =
  context "?HPR2" $ do
    describe "CHRP2" $ do
      it "4*4 a upper (column oriented)" $ do
        matvecTest1CHPR2
    describe "ZHRP2" $ do
      it "4*4 a upper (row oriented)" $ do
        matvecTest1ZHPR2

-- [12:+0 ...]
-- [.       ...]
-- [.       ...]
-- [.       ...]
--
matvecTest1CHPR2 :: IO ()
matvecTest1CHPR2 = do
  a <- Matrix.generateMutableDenseVector 10 (\_ -> 0.0:+0.0)
  x <- Matrix.generateMutableDenseVector 4 (\idx -> [1.0:+1.0, 1.0:+1.0, 1.0:+1.0, 1.0:+1.0] !! idx)
  y <- Matrix.generateMutableDenseVector 4 (\idx -> [1.0:+2.0, 1.0:+2.0, 1.0:+2.0, 1.0:+2.0] !! idx)
  BLAS.chpr2 Matrix.SColumn Matrix.MatUpper 4 2.0 x y a
  resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector a
  resList `shouldBe` [12.0:+0.0, 12.0:+0.0, 12.0:+0.0, 12.0:+0.0, 12.0:+0.0, 12.0:+0.0, 12.0:+0.0, 12.0:+0.0, 12.0:+0.0, 12.0:+0.0]

matvecTest1ZHPR2 :: IO ()
matvecTest1ZHPR2 = do
  a <- Matrix.generateMutableDenseVector 10 (\_ -> 0.0:+0.0)
  x <- Matrix.generateMutableDenseVector 4 (\idx -> [1.0:+1.0, 1.0:+1.0, 1.0:+1.0, 1.0:+1.0] !! idx)
  y <- Matrix.generateMutableDenseVector 4 (\idx -> [1.0:+2.0, 1.0:+2.0, 1.0:+2.0, 1.0:+2.0] !! idx)
  BLAS.zhpr2 Matrix.SColumn Matrix.MatUpper 4 2.0 x y a
  resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector a
  resList `shouldBe` [12.0:+0.0, 12.0:+0.0, 12.0:+0.0, 12.0:+0.0, 12.0:+0.0, 12.0:+0.0, 12.0:+0.0, 12.0:+0.0, 12.0:+0.0, 12.0:+0.0]

-- [1 2 0 0]
-- [2 3 4 0]
-- [0 4 5 6]
-- [0 0 6 7]
--
-- [1 2]
-- [3 4]
-- [5 6]
-- [7 0]
--
-- [5]
-- [11]
-- [17]
-- [15]

sbmvSpec :: Spec
sbmvSpec =
  context "?SBMV" $ do
    describe "SSBMV" $ do
      it "4*4 a upper (row oriented)" $ do
        matvecTest1SSBMV
    describe "DSBMV" $ do
      it "4*4 a lower (column oriented)" $ do
        matvecTest1DSBMV

matvecTest1SSBMV :: IO ()
matvecTest1SSBMV = do
  a <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (4, 2) (\(x, y) -> [1, 2, 3, 4, 5, 6, 7, 0] !! (y * 4 + x))
  x <- Matrix.generateMutableDenseVector 4 (\_ -> 1)
  y <- Matrix.generateMutableDenseVector 4 (\_ -> 2)
  BLAS.ssbmv Matrix.MatUpper 1 1.0 a x 1.0 y
  resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector y
  resList `shouldBe` [5, 11, 17, 15]

-- [1 3 5 7]
-- [2 4 6 0]

matvecTest1DSBMV :: IO ()
matvecTest1DSBMV = do
  a <- Matrix.generateMutableDenseMatrix (Matrix.SColumn) (4, 2) (\(x, y) -> [1, 3, 5, 7, 2, 4, 6, 0] !! (y * 4 + x))
  x <- Matrix.generateMutableDenseVector 4 (\_ -> 1)
  y <- Matrix.generateMutableDenseVector 4 (\_ -> 2)
  BLAS.dsbmv Matrix.MatLower 1 1.0 a x 1.0 y
  resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector y
  resList `shouldBe` [5, 11, 17, 15]

spmvSpec :: Spec
spmvSpec =
  context "?SPMV" $ do
    describe "SSPMV" $ do
      it "3*3 a upper (row oriented)" $ do
        matvecTest1SSPMV
    describe "DSPMV" $ do
      it "3*3 a lower (column oriented)" $ do
        matvecTest1DSPMV

matvecTest1SSPMV :: IO ()
matvecTest1SSPMV = do
  a <- Matrix.generateMutableDenseVector 6 (\idx -> [1, 2, 3, 4, 5, 6] !! idx)
  x <- Matrix.generateMutableDenseVector 3 (\_ -> 1)
  y <- Matrix.generateMutableDenseVector 3 (\_ -> 2)
  BLAS.sspmv Matrix.SRow Matrix.MatUpper 3 1.0 a x 1.0 y
  resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector y
  resList `shouldBe` [8, 13, 16]

matvecTest1DSPMV :: IO ()
matvecTest1DSPMV = do
  a <- Matrix.generateMutableDenseVector 6 (\idx -> [1, 2, 3, 4, 5, 6] !! idx)
  x <- Matrix.generateMutableDenseVector 3 (\_ -> 1)
  y <- Matrix.generateMutableDenseVector 3 (\_ -> 2)
  BLAS.dspmv Matrix.SColumn Matrix.MatLower 3 1.0 a x 1.0 y
  resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector y
  resList `shouldBe` [8, 13, 16]

sprSpec :: Spec
sprSpec =
  context "?SPR" $ do
    describe "SSPR" $ do
      it "3*3 a upper (row oriented)" $ do
        matvecTest1SSPR
    describe "DSPR" $ do
      it "3*3 a upper (column oriented)" $ do
        matvecTest1DSPR

matvecTest1SSPR :: IO ()
matvecTest1SSPR = do
  a <- Matrix.generateMutableDenseVector 6 (\idx -> [1, 2, 3, 4, 5, 6] !! idx)
  x <- Matrix.generateMutableDenseVector 3 (\_ -> 1)
  BLAS.sspr Matrix.SRow Matrix.MatUpper 3 1.0 x a
  resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector a
  resList `shouldBe` [2, 3, 4, 5, 6, 7]

matvecTest1DSPR :: IO ()
matvecTest1DSPR = do
  a <- Matrix.generateMutableDenseVector 6 (\idx -> [1, 2, 3, 4, 5, 6] !! idx)
  x <- Matrix.generateMutableDenseVector 3 (\_ -> 1)
  BLAS.dspr Matrix.SRow Matrix.MatUpper 3 1.0 x a
  resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector a
  resList `shouldBe` [2, 3, 4, 5, 6, 7]

spr2Spec :: Spec
spr2Spec =
  context "?SPR2" $ do
    describe "SSPR2" $ do
      it "3*3 a upper (row oriented)" $ do
        matvecTest1SSPR2
    describe "DSPR2" $ do
      it "3*3 a upper (column oriented)" $ do
        matvecTest1DSPR2

matvecTest1SSPR2 :: IO ()
matvecTest1SSPR2 = do
  a <- Matrix.generateMutableDenseVector 6 (\idx -> [1, 2, 3, 4, 5, 6] !! idx)
  x <- Matrix.generateMutableDenseVector 3 (\_ -> 1)
  y <- Matrix.generateMutableDenseVector 3 (\idx -> [1, 2, 3] !! idx)
  BLAS.sspr2 Matrix.SRow Matrix.MatUpper 3 1.0 x y a
  resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector a
  resList `shouldBe` [3, 5, 7, 8, 10, 12]

matvecTest1DSPR2 :: IO ()
matvecTest1DSPR2 = do
  a <- Matrix.generateMutableDenseVector 6 (\idx -> [1, 2, 3, 4, 5, 6] !! idx)
  x <- Matrix.generateMutableDenseVector 3 (\_ -> 1)
  y <- Matrix.generateMutableDenseVector 3 (\idx -> [1, 2, 3] !! idx)
  BLAS.dspr2 Matrix.SRow Matrix.MatUpper 3 1.0 x y a
  resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector a
  resList `shouldBe` [3, 5, 7, 8, 10, 12]

symvSpec :: Spec
symvSpec =
  context "?SYMV" $ do
    describe "SSYMV" $ do
      it "3*3 a upper (row oriented)" $ do
        matvecTest1SSYMV
    describe "DSYMV" $ do
      it "3*3 a lower (column oriented)" $ do
        matvecTest1DSYMV

-- [1 2 3]
-- [0 5 6]
-- [0 0 7]
matvecTest1SSYMV :: IO ()
matvecTest1SSYMV = do
  a <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (3, 3) (\(x, y) -> [1, 2, 3, 0, 5, 6, 0, 0, 7] !! (y * 3 + x))
  x <- Matrix.generateMutableDenseVector 3 (\_ -> 1)
  y <- Matrix.generateMutableDenseVector 3 (\_ -> 2)
  BLAS.ssymv Matrix.MatUpper 1.0 a x 1.0 y
  resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector y
  resList `shouldBe` [8, 15, 18]

matvecTest1DSYMV :: IO ()
matvecTest1DSYMV = do
  a <- Matrix.generateMutableDenseMatrix (Matrix.SColumn) (3, 3) (\(x, y) -> [1, 0, 0, 2, 5, 0, 3, 6, 7] !! (y * 3 + x))
  x <- Matrix.generateMutableDenseVector 3 (\_ -> 1)
  y <- Matrix.generateMutableDenseVector 3 (\_ -> 2)
  BLAS.dsymv Matrix.MatLower 1.0 a x 1.0 y
  resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector y
  resList `shouldBe` [8, 15, 18]

syrSpec :: Spec
syrSpec =
  context "?SYR" $ do
    describe "SSYR" $ do
      it "3*3 a upper (row oriented)" $ do
        matvecTest1SSYR
    describe "DSYR" $ do
      it "3*3 a upper (column oriented)"$ do
        matvecTest1DSYR
-- [1 2 3]
-- [0 5 6]
-- [0 0 7]
matvecTest1SSYR :: IO ()
matvecTest1SSYR = do
  a <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (3, 3) (\(x, y) -> [1, 2, 3, 0, 5, 6, 0, 0, 7] !! (y * 3 + x))
  x <- Matrix.generateMutableDenseVector 3 (\_ -> 1)
  BLAS.ssyr Matrix.MatUpper 1.0 x a
  resList <- Matrix.mutableVectorToList $ _bufferDenMutMat a
  resList `shouldBe` [2, 3, 4, 0, 6, 7, 0, 0, 8]

matvecTest1DSYR :: IO ()
matvecTest1DSYR = do
  a <- Matrix.generateMutableDenseMatrix (Matrix.SColumn) (3, 3) (\(x, y) -> [1, 0, 0, 2, 5, 0, 3, 6, 7] !! (y * 3 + x))
  x <- Matrix.generateMutableDenseVector 3 (\_ -> 1)
  BLAS.dsyr Matrix.MatLower 1.0 x a
  resList <- Matrix.mutableVectorToList $ _bufferDenMutMat a
  resList `shouldBe` [2, 3, 4, 0, 6, 7, 0, 0, 8]

syr2Spec :: Spec
syr2Spec =
  context "?SYR2" $ do
    describe "SSYR2" $ do
      it "3*3 a upper (row oriented)" $ do
        matvecTest1SSYR2
    describe "DSYR2" $ do
      it "3*3 a upper (column oriented)" $ do
        matvecTest1DSYR2

matvecTest1SSYR2 :: IO ()
matvecTest1SSYR2 = do
  a <- Matrix.generateMutableDenseMatrix (Matrix.SRow) (3, 3) (\(x, y) -> [1, 2, 3, 0, 5, 6, 0, 0, 7] !! (y * 3 + x))
  x <- Matrix.generateMutableDenseVector 3 (\_ -> 1)
  y <- Matrix.generateMutableDenseVector 3 (\idx -> [1, 2, 3] !! idx)
  BLAS.ssyr2 Matrix.MatUpper 1.0 x y a
  resList <- Matrix.mutableVectorToList $ _bufferDenMutMat a
  resList `shouldBe` [3, 5, 7, 0, 9, 11, 0, 0, 13]

matvecTest1DSYR2 :: IO ()
matvecTest1DSYR2 = do
  a <- Matrix.generateMutableDenseMatrix (Matrix.SColumn) (3, 3) (\(x, y) -> [1, 0, 0, 2, 5, 0, 3, 6, 7] !! (y * 3 + x))
  x <- Matrix.generateMutableDenseVector 3 (\_ -> 1)
  y <- Matrix.generateMutableDenseVector 3 (\idx -> [1, 2, 3] !! idx)
  BLAS.dsyr2 Matrix.MatLower 1.0 x y a
  resList <- Matrix.mutableVectorToList $ _bufferDenMutMat a
  resList `shouldBe` [3, 5, 7, 0, 9, 11, 0, 0, 13]

tbmvSpec :: Spec
tbmvSpec =
  context "?TBMV" $ do
    describe "STBMV" $ do
      it "3x3 upper no trans (row oriented)" $ do
        matmatTest1STBMV
    describe "DTBMV" $ do
      it "3x3 lower trans (column oriented)" $ do
        matmatTest1DTBMV
    describe "CTBMV" $ do
      it "3x3 upper conj trans (row oriented)" $ do
        matmatTest1CTBMV
    describe "ZTBMV" $ do
      it "3x3 lower no trans (row oriented)" $ do
        matmatTest1ZTBMV

-- [ 1 2 0]   [-2]   [4]
-- [ 0 1 1] * [ 3] = [5]
-- [ 0 0 1]   [ 2]   [2]
matmatTest1STBMV:: IO ()
matmatTest1STBMV = do
    a  <- Matrix.generateMutableDenseMatrix (Matrix.SRow)  (3, 2) (\(x, y) -> [0, 2, 1, 1, 1, 1] !! (y * 3 + x))
    x  <- Matrix.generateMutableDenseVector 3 (\idx -> [-2, 3, 2] !! idx)
    BLAS.stbmv MatUpper NoTranspose MatUnit 1 a x
    resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector x
    resList `shouldBe` [4, 5, 2]

-- [1 1 0]   [-2]   [1]
-- [0 1 1] * [ 3] = [5]
-- [0 0 1]   [ 2]   [2]
matmatTest1DTBMV:: IO ()
matmatTest1DTBMV = do
    a  <- Matrix.generateMutableDenseMatrix (Matrix.SColumn)  (3, 2) (\(x, y) -> [0, 1, 1, 1, 1, 1] !! (y * 3 + x))
    x  <- Matrix.generateMutableDenseVector 3 (\idx -> [-2, 3, 2] !! idx)
    BLAS.dtbmv MatLower Transpose MatUnit 1 a x
    resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector x
    resList `shouldBe` [1, 5, 2]

-- [ 1:+0    0:+0    0:+0]   [-2]   [-2:+0   ]
-- [ 2:+(-2) 1:+0    0:+0] * [ 3] = [-1:+4   ]
-- [ 0:+0    1:+(-1) 1:+0]   [ 2]   [ 5:+(-3)]
matmatTest1CTBMV:: IO ()
matmatTest1CTBMV = do
    a  <- Matrix.generateMutableDenseMatrix (Matrix.SRow)  (3, 2) (\(x, y) -> [0:+0, 2:+2, 1:+1, 1:+1, 1:+1, 1:+1] !! (y * 3 + x))
    x  <- Matrix.generateMutableDenseVector 3 (\idx -> [-2, 3, 2] !! idx)
    BLAS.ctbmv MatUpper Matrix.ConjTranspose MatUnit 1 a x
    resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector x
    resList `shouldBe` [(-2):+0, (-1):+4, 5:+(-3)]

-- [1 1 0]   [2]   [5]
-- [0 1 1] * [3] = [5]
-- [0 0 1]   [2]   [2]
matmatTest1ZTBMV:: IO ()
matmatTest1ZTBMV = do
    a  <- Matrix.generateMutableDenseMatrix (Matrix.SColumn)  (3, 2) (\(x, y) -> [0:+0, 1:+0, 1:+0, 1:+0, 1:+0, 1:+0] !! (y * 3 + x))
    x  <- Matrix.generateMutableDenseVector 3 (\idx -> [2, 3, 2] !! idx)
    BLAS.ztbmv MatUpper NoTranspose MatNonUnit 1 a x -- TODO: NAN error when using Lower Transpose
    resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector x
    resList `shouldBe` [5, 5, 2]

tbsvSpec :: Spec
tbsvSpec =
  context "?TBSV" $ do
    describe "STBSV" $ do
      it "3x3 upper no trans (row oriented)" $ do
        matmatTest1STBSV
    describe "DTBSV" $ do
      it "3x3 lower trans (column oriented)" $ do
        matmatTest1DTBSV
    describe "CTBSV" $ do
      it "3x3 upper conj trans (row oriented)" $ do
        matmatTest1CTBSV
    describe "ZTBSV" $ do
      it "3x3 lower no trans (row oriented)" $ do
        matmatTest1ZTBSV

-- [ 1 2 0]   [-2]   [4]
-- [ 0 1 1] * [ 3] = [5]
-- [ 0 0 1]   [ 2]   [2]
matmatTest1STBSV:: IO ()
matmatTest1STBSV = do
    a  <- Matrix.generateMutableDenseMatrix (Matrix.SRow)  (3, 2) (\(x, y) -> [0, 2, 1, 1, 1, 1] !! (y * 3 + x))
    x  <- Matrix.generateMutableDenseVector 3 (\idx -> [4, 5, 2] !! idx)
    BLAS.stbsv MatUpper NoTranspose MatUnit 1 a x
    resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector x
    resList `shouldBe` [-2, 3, 2]

-- [1 1 0]   [-2]   [1]
-- [0 1 1] * [ 3] = [5]
-- [0 0 1]   [ 2]   [2]
matmatTest1DTBSV:: IO ()
matmatTest1DTBSV = do
    a  <- Matrix.generateMutableDenseMatrix (Matrix.SColumn)  (3, 2) (\(x, y) -> [0, 1, 1, 1, 1, 1] !! (y * 3 + x))
    x  <- Matrix.generateMutableDenseVector 3 (\idx -> [1, 5, 2] !! idx)
    BLAS.dtbsv MatLower Transpose MatUnit 1 a x
    resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector x
    resList `shouldBe` [-2, 3, 2]

-- [ 1:+0    0:+0    0:+0]   [-2]   [-2:+0   ]
-- [ 2:+(-2) 1:+0    0:+0] * [ 3] = [-1:+4   ]
-- [ 0:+0    1:+(-1) 1:+0]   [ 2]   [ 5:+(-3)]
matmatTest1CTBSV:: IO ()
matmatTest1CTBSV = do
    a  <- Matrix.generateMutableDenseMatrix (Matrix.SRow)  (3, 2) (\(x, y) -> [0:+0, 2:+2, 1:+1, 1:+1, 1:+1, 1:+1] !! (y * 3 + x))
    x  <- Matrix.generateMutableDenseVector 3 (\idx -> [(-2):+0, (-1):+4, 5:+(-3)] !! idx)
    BLAS.ctbsv MatUpper Matrix.ConjTranspose MatUnit 1 a x
    resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector x
    resList `shouldBe` [-2, 3, 2]

-- [1 1 0]   [2]   [5]
-- [0 1 1] * [3] = [5]
-- [0 0 1]   [2]   [2]
matmatTest1ZTBSV:: IO ()
matmatTest1ZTBSV = do
    a  <- Matrix.generateMutableDenseMatrix (Matrix.SColumn)  (3, 2) (\(x, y) -> [0:+0, 1:+0, 1:+0, 1:+0, 1:+0, 1:+0] !! (y * 3 + x))
    x  <- Matrix.generateMutableDenseVector 3 (\idx -> [5, 5, 2] !! idx)
    BLAS.ztbsv MatUpper NoTranspose MatNonUnit 1 a x -- TODO: NAN error when using Lower Transpose
    resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector x
    resList `shouldBe` [2, 3, 2]

tpmvSpec :: Spec
tpmvSpec =
  context "?TPMV" $ do
    describe "STPMV" $ do
      it "3x3 upper no trans (row oriented)" $ do
        matmatTest1STPMV
    describe "DTPMV" $ do
      it "3x3 lower trans (column oriented)" $ do
        matmatTest1DTPMV
    describe "CTPMV" $ do
      it "3x3 upper conj trans (row oriented)" $ do
        matmatTest1CTPMV
    describe "ZTPMV" $ do
      it "3x3 lower no trans (row oriented)" $ do
        matmatTest1ZTPMV

-- [ 1 2 0]   [-2]   [4]
-- [ 0 1 1] * [ 3] = [5]
-- [ 0 0 1]   [ 2]   [2]
matmatTest1STPMV:: IO ()
matmatTest1STPMV = do
    a  <- Matrix.generateMutableDenseVector 6 (\idx -> [1, 2, 0, 1, 1, 1] !! idx)
    x  <- Matrix.generateMutableDenseVector 3 (\idx -> [-2, 3, 2] !! idx)
    BLAS.stpmv Matrix.SRow MatUpper NoTranspose MatUnit 3 a x
    resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector x
    resList `shouldBe` [4, 5, 2]

-- [1 1 0]   [-2]   [1]
-- [0 1 1] * [ 3] = [5]
-- [0 0 1]   [ 2]   [2]
matmatTest1DTPMV:: IO ()
matmatTest1DTPMV = do
    a  <- Matrix.generateMutableDenseVector 6 (\idx -> [1, 1, 0, 1, 1, 1] !! idx)
    x  <- Matrix.generateMutableDenseVector 3 (\idx -> [-2, 3, 2] !! idx)
    BLAS.dtpmv Matrix.SColumn MatLower Transpose MatUnit 3 a x
    resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector x
    resList `shouldBe` [1, 5, 2]

-- [ 1:+0    0:+0    0:+0]   [-2]   [-2:+0   ]
-- [ 2:+(-2) 1:+0    0:+0] * [ 3] = [-1:+4   ]
-- [ 0:+0    1:+(-1) 1:+0]   [ 2]   [ 5:+(-3)]
matmatTest1CTPMV:: IO ()
matmatTest1CTPMV = do
    a  <- Matrix.generateMutableDenseVector 6 (\idx -> [1:+0, 2:+2, 0:+0, 1:+0, 1:+1, 1:+0] !! idx)
    x  <- Matrix.generateMutableDenseVector 3 (\idx -> [-2, 3, 2] !! idx)
    BLAS.ctpmv Matrix.SRow MatUpper Matrix.ConjTranspose MatUnit 3 a x
    resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector x
    resList `shouldBe` [(-2):+0, (-1):+4, 5:+(-3)]

-- [1 1 0]   [2]   [5]
-- [0 1 1] * [3] = [5]
-- [0 0 1]   [2]   [2]
matmatTest1ZTPMV:: IO ()
matmatTest1ZTPMV = do
    a  <- Matrix.generateMutableDenseVector 6 (\idx -> [1:+0, 1:+0, 1:+0, 0:+0, 1:+0, 1:+0] !! idx)
    x  <- Matrix.generateMutableDenseVector 3 (\idx -> [2, 3, 2] !! idx)
    BLAS.ztpmv Matrix.SColumn MatUpper NoTranspose MatNonUnit 3 a x -- TODO: NAN error when using Lower Transpose
    resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector x
    resList `shouldBe` [5, 5, 2]

tpsvSpec :: Spec
tpsvSpec =
  context "?TPSV" $ do
    describe "STPSV" $ do
      it "3x3 upper no trans (row oriented)" $ do
        matmatTest1STPSV
    describe "DTPSV" $ do
      it "3x3 lower trans (column oriented)" $ do
        matmatTest1DTPSV
    describe "CTPSV" $ do
      it "3x3 upper conj trans (row oriented)" $ do
        matmatTest1CTPSV
    describe "ZTPSV" $ do
      it "3x3 lower no trans (row oriented)" $ do
        matmatTest1ZTPSV

-- [ 1 2 0]   [-2]   [4]
-- [ 0 1 1] * [ 3] = [5]
-- [ 0 0 1]   [ 2]   [2]
matmatTest1STPSV:: IO ()
matmatTest1STPSV = do
    a  <- Matrix.generateMutableDenseVector 6 (\idx -> [1, 2, 0, 1, 1, 1] !! idx)
    x  <- Matrix.generateMutableDenseVector 3 (\idx -> [4, 5, 2] !! idx)
    BLAS.stpsv Matrix.SRow MatUpper NoTranspose MatUnit 3 a x
    resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector x
    resList `shouldBe` [-2, 3, 2]

-- [1 1 0]   [-2]   [1]
-- [0 1 1] * [ 3] = [5]
-- [0 0 1]   [ 2]   [2]
matmatTest1DTPSV:: IO ()
matmatTest1DTPSV = do
    a  <- Matrix.generateMutableDenseVector 6 (\idx -> [1, 1, 0, 1, 1, 1] !! idx)
    x  <- Matrix.generateMutableDenseVector 3 (\idx -> [1, 5, 2] !! idx)
    BLAS.dtpsv Matrix.SColumn MatLower Transpose MatUnit 3 a x
    resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector x
    resList `shouldBe` [-2, 3, 2]

-- [ 1:+0    0:+0    0:+0]   [-2]   [-2:+0   ]
-- [ 2:+(-2) 1:+0    0:+0] * [ 3] = [-1:+4   ]
-- [ 0:+0    1:+(-1) 1:+0]   [ 2]   [ 5:+(-3)]
matmatTest1CTPSV:: IO ()
matmatTest1CTPSV = do
    a  <- Matrix.generateMutableDenseVector 6 (\idx -> [1:+0, 2:+2, 0:+0, 1:+0, 1:+1, 1:+0] !! idx)
    x  <- Matrix.generateMutableDenseVector 3 (\idx -> [(-2):+0, (-1):+4, 5:+(-3)]!! idx)
    BLAS.ctpsv Matrix.SRow MatUpper Matrix.ConjTranspose MatUnit 3 a x
    resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector x
    resList `shouldBe` [-2, 3, 2]

-- [1 1 0]   [2]   [5]
-- [0 1 1] * [3] = [5]
-- [0 0 1]   [2]   [2]
matmatTest1ZTPSV:: IO ()
matmatTest1ZTPSV = do
    a  <- Matrix.generateMutableDenseVector 6 (\idx -> [1:+0, 1:+0, 1:+0, 0:+0, 1:+0, 1:+0] !! idx)
    x  <- Matrix.generateMutableDenseVector 3 (\idx -> [5, 5, 2] !! idx)
    BLAS.ztpsv Matrix.SColumn MatUpper NoTranspose MatNonUnit 3 a x -- TODO: NAN error when using Lower Transpose
    resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector x
    resList `shouldBe` [2, 3, 2]

trmvSpec :: Spec
trmvSpec =
  context "?TRMV" $ do
    describe "STRMV" $ do
      it "2x2 upper 1's" $ do
        matmatTest1STRMV
    describe "DTRMV" $ do
      it "2x2 upper 1's" $ do
        matmatTest1DTRMV
    describe "CTRMV" $ do
      it "2x2 upper 1's" $ do
        matmatTest1CTRMV
    describe "ZTRMV" $ do
      it "2x2 upper 1's" $ do
        matmatTest1ZTRMV

matmatTest1STRMV:: IO ()
matmatTest1STRMV = do
    left  <- Matrix.generateMutableDenseMatrix (Matrix.SRow)  (2,2)
            (\(i,j) -> if i >= j then (1.0::Float) else 0 )

    res  <- Matrix.generateMutableDenseVector  2 (\i -> if i == 0 then 2 else 1)
    BLAS.strmv MatUpper NoTranspose MatUnit left res
    resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector res
    resList `shouldBe` [3,1]

matmatTest1DTRMV:: IO ()
matmatTest1DTRMV = do
    left  <- Matrix.generateMutableDenseMatrix (Matrix.SRow)  (2,2)
                (\(i,j) -> if i >= j then (1.0::Double) else 0 )

    res  <- Matrix.generateMutableDenseVector  2 (\i -> if i == 0 then 2 else 1)
    BLAS.dtrmv MatUpper NoTranspose MatUnit left res
    resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector res
    resList `shouldBe` [3,1]

matmatTest1CTRMV:: IO ()
matmatTest1CTRMV = do
    left  <- Matrix.generateMutableDenseMatrix (Matrix.SRow)  (2,2)
                (\(i,j) -> if i >= j then (1.0::(Complex Float)) else 0 )

    res  <- Matrix.generateMutableDenseVector  2 (\i -> if i == 0 then 2 else 1)
    BLAS.ctrmv MatUpper NoTranspose MatUnit left res
    resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector res
    resList `shouldBe` [3,1]

matmatTest1ZTRMV:: IO ()
matmatTest1ZTRMV = do
    left  <- Matrix.generateMutableDenseMatrix (Matrix.SRow)  (2,2)
                (\(i,j) -> if i >= j then (1.0::(Complex Double )) else 0 )
    res  <- Matrix.generateMutableDenseVector  2 (\i -> if i == 0 then 2 else 1)
    BLAS.ztrmv MatUpper NoTranspose MatUnit left res
    resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector res
    resList `shouldBe` [3,1]

trsvSpec :: Spec
trsvSpec =
  context "?TRSV" $ do
    describe "STRSV" $ do
      it "2x2 upper 1's" $ do
        matmatTest1STRSV
    describe "DTRSV" $ do
      it "2x2 upper 1's" $ do
        matmatTest1DTRSV
    describe "CTRSV" $ do
      it "2x2 upper 1's" $ do
        matmatTest1CTRSV
    describe "ZTRSV" $ do
      it "2x2 upper 1's" $ do
        matmatTest1ZTRSV

matmatTest1STRSV:: IO ()
matmatTest1STRSV = do
    left  <- Matrix.generateMutableDenseMatrix (Matrix.SRow)  (2,2)
            (\(i,j) -> if i >= j then (1.0::Float) else 0 )

    res  <- Matrix.generateMutableDenseVector  2 (\i -> if i == 0 then 3 else 1)
    BLAS.strsv MatUpper NoTranspose MatUnit left res
    resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector res
    resList `shouldBe` [2,1]

matmatTest1DTRSV:: IO ()
matmatTest1DTRSV = do
    left  <- Matrix.generateMutableDenseMatrix (Matrix.SRow)  (2,2)
                (\(i,j) -> if i >= j then (1.0::Double) else 0 )

    res  <- Matrix.generateMutableDenseVector  2 (\i -> if i == 0 then 3 else 1)
    BLAS.dtrsv MatUpper NoTranspose MatUnit left res
    resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector res
    resList `shouldBe` [2,1]

matmatTest1CTRSV:: IO ()
matmatTest1CTRSV = do
    left  <- Matrix.generateMutableDenseMatrix (Matrix.SRow)  (2,2)
                (\(i,j) -> if i >= j then (1.0::(Complex Float)) else 0 )

    res  <- Matrix.generateMutableDenseVector  2 (\i -> if i == 0 then 3 else 1)
    BLAS.ctrsv MatUpper NoTranspose MatUnit left res
    resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector res
    resList `shouldBe` [2,1]

matmatTest1ZTRSV:: IO ()
matmatTest1ZTRSV = do
    left  <- Matrix.generateMutableDenseMatrix (Matrix.SRow)  (2,2)
                (\(i,j) -> if i >= j then (1.0::(Complex Double )) else 0 )
    res  <- Matrix.generateMutableDenseVector  2 (\i -> if i == 0 then 3 else 1)
    BLAS.ztrsv MatUpper NoTranspose MatUnit left res
    resList <- Matrix.mutableVectorToList $ _bufferMutDenseVector res
    resList `shouldBe` [2,1]

max_stars_repo_name: schnecki/hblas

max_issues_repo_name: schnecki/hblas

max_forks_repo_name: schnecki/hblas