{-# LANGUAGE TypeFamilies #-}
module Test.Multiply (
   multiplySquare,
   squareSquare,
   power,
   ) where

import qualified Test.Utility as Util

import qualified Numeric.LAPACK.Matrix.Square as Square
import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix
import qualified Numeric.LAPACK.Matrix as Matrix
import qualified Numeric.LAPACK.Vector as Vector
import Numeric.LAPACK.Matrix (Matrix, ShapeInt, (#*##))
import Numeric.LAPACK.Scalar (RealOf)

import qualified Numeric.Netlib.Class as Class



multiplySquare ::
   (Matrix.Power typ, Matrix.MultiplySquare typ,
    Matrix.SquareShape typ, Matrix.HeightOf typ ~ ShapeInt,
    Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   Matrix typ a -> Bool
multiplySquare a =
   Util.approxArray -- (Scalar.selectReal 1e-1 1e-5)
      (Matrix.toSquare $ Matrix.square a)
      (a #*## Matrix.toSquare a)

squareSquare ::
   (Matrix.Power typ, Matrix.MultiplySquare typ,
    Matrix.SquareShape typ, Matrix.HeightOf typ ~ ShapeInt,
    Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   Matrix typ a -> Bool
squareSquare a =
   Util.approxArray -- (Scalar.selectReal 1e-1 1e-5)
      (Matrix.toSquare $ Matrix.square a)
      (Square.square $ Matrix.toSquare a)

power ::
   (Matrix.Power typ, Matrix.MultiplySquare typ,
    Matrix.SquareShape typ, Matrix.HeightOf typ ~ ShapeInt,
    Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   Int -> Matrix typ a -> Bool
power n a =
   let b = Matrix.toSquare (Matrix.power (n+1) a)
       c = a #*## Matrix.toSquare (Matrix.power n a)
       normInf1 = Vector.normInf1 . ArrMatrix.toVector
   in Util.approxArrayTol (1e-6 * (normInf1 b + normInf1 c)) b c