@@ -103,8 +103,8 @@ unsafeGetColPerm :: (WriteVector x m, Elem e) => Perm (n,p) e -> Int -> m (x n e) unsafeGetColPerm (Identity n) j = newBasisVector n j unsafeGetColPerm (Perm h sigma) j - | h == NoTrans = newBasisVector n (P.unsafeAt sigma j) - | otherwise = newBasisVector n (P.indexOf sigma j) + | h == NoTrans = newBasisVector n (P.indexOf sigma j) + | otherwise = newBasisVector n (P.unsafeAt sigma j) where n = P.size sigma {-# INLINE unsafeGetColPerm #-} @@ -153,11 +153,11 @@ unsafeDoSApplyAddPerm alpha p@(Perm h sigma) x beta y forM_ [0..(n-1)] $ \i -> let i' = P.unsafeAt sigma i in case h of - NoTrans -> do + ConjTrans -> do e <- unsafeReadElem x i f <- unsafeReadElem y i' unsafeWriteElem y i' (alpha*e + f) - ConjTrans -> do + NoTrans -> do e <- unsafeReadElem x i' f <- unsafeReadElem y i unsafeWriteElem y i (alpha*e + f) @@ -175,9 +175,9 @@ unsafeDoSApplyAddMatPerm alpha p@(Perm h sigma) b beta c = forM_ [0..(m-1)] $ \i -> let i' = P.unsafeAt sigma i in case h of - NoTrans -> unsafeAxpyVector alpha (unsafeRowView b i) + ConjTrans -> unsafeAxpyVector alpha (unsafeRowView b i) (unsafeRowView c i') - ConjTrans -> unsafeAxpyVector alpha (unsafeRowView b i') + NoTrans -> unsafeAxpyVector alpha (unsafeRowView b i') (unsafeRowView c i) {-# INLINE unsafeDoSApplyAddMatPerm #-} lib/Data/Matrix/Perm.hs @@ -9,7 +9,7 @@ import Orthogonal( tests_Orthogonal ) main :: IO () main = do args <- getArgs - let n = if null args then 1000 else read (head args) + let n = if null args then 100 else read (head args) printf "\nRunnings tests for field `%s'\n" field tests/Main.hs @@ -5,20 +5,126 @@ module Orthogonal import Driver import Monadic import Test.QuickCheck hiding ( vector ) -import Test.QuickCheck.BLAS( Pos(..), Nat2(..) ) +import qualified Test.QuickCheck as QC +import Test.QuickCheck.BLAS( Index(..), Pos(..), Nat(..), Nat2(..) ) import qualified Test.QuickCheck.BLAS as Test import Control.Monad +import Data.Permute( Permute ) +import qualified Data.Permute as P + import Data.Elem.BLAS import Data.Vector.Dense import Data.Vector.Dense.ST import Data.Matrix.Dense import Data.Matrix.Dense.ST import Data.Matrix.House +import Data.Matrix.Perm import Data.Matrix.QR -import Debug.Trace + +testPermute :: Int -> Gen Permute +testPermute n = do + es <- QC.vector n :: Gen [Int] + return $ P.order n es + +testPerm :: Int -> Gen (Perm (n,n) e) +testPerm n = do + p <- liftM permFromPermute $ testPermute n + elements [ identityPerm n, p, herm p ] + +prop_perm_herm (Nat n) = + forAll (testPerm n) $ \p -> + permuteFromPerm (herm p) == P.inverse (permuteFromPerm p) + +prop_perm_col (Index n i) = + forAll (testPerm n) $ \p -> + col p i + === + p <*> (basisVector n i :: V) + +prop_perm_apply_basis (Index n i) = + forAll (testPerm n) $ \p -> + p <*> (basisVector n i :: V) + === + basisVector n (P.indexOf (permuteFromPerm p) i) + +prop_perm_herm_apply (Nat n) = + forAll (testPerm n) $ \p -> + forAll (Test.vector n) $ \(x :: V) -> + p <*> herm p <*> x === x + +prop_herm_perm_apply (Nat n) = + forAll (testPerm n) $ \p -> + forAll (Test.vector n) $ \(x :: V) -> + herm p <*> p <*> x === x + +prop_perm_solve (Nat n) = + forAll (testPerm n) $ \p -> + forAll (Test.vector n) $ \(x :: V) -> + p <\> x === herm p <*> x + +prop_perm_applyMat_cols (Nat2 (m,n)) = + forAll (testPerm m) $ \p -> + forAll (Test.matrix (m,n)) $ \(a :: M) -> + cols (p <**> a) === map (p <*>) (cols a) + +prop_perm_herm_applyMat (Nat2 (m,n)) = + forAll (testPerm m) $ \p -> + forAll (Test.matrix (m,n)) $ \(a :: M) -> + p <**> herm p <**> a === a + +prop_herm_perm_applyMat (Nat2 (m,n)) = + forAll (testPerm m) $ \p -> + forAll (Test.matrix (m,n)) $ \(a :: M) -> + herm p <**> p <**> a === a + +prop_perm_solveMat_cols (Nat2 (m,n)) = + forAll (testPerm m) $ \p -> + forAll (Test.matrix (m,n)) $ \(a :: M) -> + cols (p <\\> a) === map (p <\>) (cols a) + +prop_perm_solveMat (Nat2 (m,n)) = + forAll (testPerm m) $ \p -> + forAll (Test.matrix (m,n)) $ \(a :: M) -> + p <\\> a === herm p <**> a + +prop_perm_doSApplyVector_ alpha (Nat n) = + forAll (testPerm n) $ \p -> + forAll (Test.vector n) $ \(x :: V) -> + monadicST $ do + x' <- run $ unsafeThawVector x + x'' <- run $ freezeVector x' + run $ doSApplyVector_ alpha p x' + assert $ x ~== p <*> (alpha *> x'') + +prop_perm_doSApplyMatrix_ alpha (Nat2 (m,n)) = + forAll (testPerm m) $ \p -> + forAll (Test.matrix (m,n)) $ \(b :: M) -> + monadicST $ do + b' <- run $ unsafeThawMatrix b + b'' <- run $ freezeMatrix b' + run $ doSApplyMatrix_ alpha p b' + assert $ b ~== p <**> (alpha *> b'') + +prop_perm_doSSolveVector alpha (Nat n) = + forAll (testPerm n) $ \p -> + forAll (Test.vector n) $ \(x :: V) -> + forAll (Test.vector n) $ \y -> + monadicST $ do + x' <- run $ unsafeThawVector x + run $ doSSolveVector alpha p y x' + assert $ x ~== p <\> (alpha *> y) + +prop_perm_doSSolveMatrix alpha (Nat2 (m,n)) = + forAll (testPerm m) $ \p -> + forAll (Test.matrix (m,n)) $ \(b :: M) -> + forAll (Test.matrix (m,n)) $ \c -> + monadicST $ do + b' <- run $ unsafeThawMatrix b + run $ doSSolveMatrix alpha p c b' + assert $ b ~== p <\\> (alpha *> c) prop_setReflector_snd (Pos n) = monadicST $ do @@ -96,7 +202,24 @@ prop_qrFactor_doSSolveMatrix alpha (Nat2 (n,p)) = assert $ b ~== qr <\\> (alpha *> c) tests_Orthogonal = - [ ("snd . setReflector", mytest prop_setReflector_snd) + [ + ("perm herm", mytest prop_perm_herm) + , ("perm col", mytest prop_perm_col) + , ("perm apply basis", mytest prop_perm_apply_basis) + , ("perm herm apply", mytest prop_perm_herm_apply) + , ("herm perm apply", mytest prop_herm_perm_apply) + , ("perm solve", mytest prop_perm_solve) + , ("perm applyMat cols", mytest prop_perm_applyMat_cols) + , ("perm herm applyMat", mytest prop_perm_herm_applyMat) + , ("herm perm applyMat", mytest prop_herm_perm_applyMat) + , ("perm solveMat cols", mytest prop_perm_solveMat_cols) + , ("perm solveMat", mytest prop_perm_solveMat) + , ("perm doApplyVector_", mytest prop_perm_doSApplyVector_) + , ("perm doApplyMatrix_", mytest prop_perm_doSApplyMatrix_) + , ("perm doSolveVector", mytest prop_perm_doSSolveVector) + , ("perm doSolveMatrix", mytest prop_perm_doSSolveMatrix) + + , ("snd . setReflector", mytest prop_setReflector_snd) , ("fst . reflector", mytest prop_reflector_fst) , ("reflector <*>", mytest prop_reflector_vector) , ("reflector <**>", mytest prop_reflector_matrix) tests/Orthogonal.hs 65e96e1cbfebad34e8fdb4aaf251640971c9be50 Patrick Perry patperry@stanford.edu http://github.com/patperry/lapack/commit/e47271400903cc1e770e19af6487518e43278204 e47271400903cc1e770e19af6487518e43278204 2009-01-24T20:30:39-08:00 2009-01-24T20:30:39-08:00 add tests for Perm and fix bugs 7fcb7b5fe20f81a9ecb2ca5e477b5726a5d80a03 Patrick Perry patperry@stanford.edu