packages feed

lapack-0.4: test/Numeric/LAPACK/Matrix/BandedHermitian/Naive.hs

{-# LANGUAGE TypeFamilies #-}
module Numeric.LAPACK.Matrix.BandedHermitian.Naive where

import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix
import qualified Numeric.LAPACK.Matrix.Layout as Layout
import qualified Numeric.LAPACK.Matrix.Extent as Extent
import qualified Numeric.LAPACK.Matrix.Banded as Banded
import qualified Numeric.LAPACK.Matrix as Matrix
import Numeric.LAPACK.Matrix.Layout (Order)
import Numeric.LAPACK.Matrix ((#!))

import qualified Numeric.Netlib.Class as Class

import qualified Type.Data.Num.Unary as Unary

import qualified Data.Array.Comfort.Storable as Array
import qualified Data.Array.Comfort.Shape as Shape
import Data.Array.Comfort.Storable ((!))


toHermitian ::
   (Unary.Natural offDiag, Shape.Indexed size, Class.Floating a) =>
   Banded.Hermitian offDiag size a -> Matrix.Hermitian size a
toHermitian = ArrMatrix.lift1 $ \a ->
   case Array.shape a of
      shape@(Layout.BandedHermitian _offDiag order size) ->
         Array.sample (Layout.hermitian order size) $ \ix ->
            let bix = uncurry Layout.InsideBox ix
            in if Shape.inBounds shape bix then a!bix else 0

toBanded ::
   (Unary.Natural offDiag, Shape.Indexed size, Class.Floating a) =>
   Banded.Hermitian offDiag size a -> Banded.Square offDiag offDiag size a
toBanded a = flip ArrMatrix.lift1 a $ \av ->
   case Array.shape av of
      shape@(Layout.BandedHermitian offDiag order size) ->
         Array.sample
            (Layout.Banded (offDiag,offDiag) order (Extent.square size)) $
         \ix ->
            case ix of
               Layout.InsideBox r c -> a#!(r,c)
               _ -> if Shape.inBounds shape ix then av!ix else 0

forceOrder ::
   (Unary.Natural offDiag, Shape.Indexed size, Class.Floating a) =>
   Order ->
   Banded.Hermitian offDiag size a ->
   Banded.Hermitian offDiag size a
forceOrder newOrder = ArrMatrix.lift1 $ \a ->
   case Array.shape a of
      shape ->
         Array.sample (shape{Layout.bandedHermitianOrder = newOrder}) $
         \ix ->
            case ix of
               Layout.InsideBox _ _ -> a!ix
               _ ->
                  if newOrder == Layout.bandedHermitianOrder shape
                     then a!ix
                     else 0