update to version from Seq no More paper
authorSimon Marlow <marlowsd@gmail.com>
Fri, 5 Nov 2010 13:26:32 +0000 (13:26 +0000)
committerSimon Marlow <marlowsd@gmail.com>
Fri, 5 Nov 2010 13:26:32 +0000 (13:26 +0000)
parallel/linsolv/CRA.hs
parallel/linsolv/LinSolv.hs [deleted file]
parallel/linsolv/Main.hs
parallel/linsolv/Makefile
parallel/linsolv/Matrix.hs
parallel/linsolv/ModArithm.hs
parallel/linsolv/ParForce.hs [deleted file]

index 55c8445..86c1612 100644 (file)
@@ -1,95 +1,64 @@
-{-
-$Id: CRA.hs,v 1.1 1996/01/08 20:07:34 partain Exp $
-
-This is revision: $Revision: 1.1 $
-
-Module realizing the Chinese Remainder Algorithm (on 2 inputs as well as
-on 2 input lists).
-
-Changelog:
-$Log: CRA.hs,v $
-Revision 1.1  1996/01/08 20:07:34  partain
-Initial revision
-
---# Revision 1.2  1994/12/13  19:53:03  hwloidl
---# tree_IMCRA0 is a tree-based version of CRA over lists.
---# Nomally it works over infinite lists (with -DNO_FILTER over finite lists;
---# NOTE: Depending on NO_FILTER different tree_IMCRA0 functions are compiled
---#       with *different typings*. So, any module using tree_IMCRA0 must be
---#       recompiled, too, when changing this flag).
---# s2IMCRA is the same as s1IMCRA (i.e. fold-based) but works on infinite lists.
---# tree_IMCRA0 handles the det-check now (therefore list of dets as input).
---# It tries to push the det-check as much backwards as possible by using the
---# handle fails fct (possible problem: only fails is immediately needed when calling
---# handle_fails)
---#
---# Revision 1.1  1994/12/08  19:24:21  hwloidl
---# Initial revision
---#
----------------------------------------------------------------------------- -}
-
-module CRA (binCRA,s1IMCRA,s2IMCRA,tree_IMCRA0,
-           isPrime) where 
-
-import ModArithm (modHom, modDif, modProd, modInv) 
+-- Time-stamp: <2010-11-03 10:31:32 simonmar>
+--
+-- Chinese Remainder Algorithm.
+-- Works over lists of Integers (2 input lists).
+----------------------------------------------------------------------------
 
-import Random  (randomInts)    -- Just for testing
-import ParForce 
+-- @node Chinese Remainder Algorithm, , ,
+-- @chapter Chinese Remainder Algorithm
 
-#ifdef SEQ
+module CRA (binCRA, seq_list_CRA0, seq_list_CRA,
+           -- par_binCRA, tree_IMCRA0,
+           isPrime) where 
 
-par :: a -> b -> b
-par x y = y
+-- @menu
+-- * Imports::                 
+-- * Auxiliary functions::     
+-- * Basic binary CRA operation::  
+-- * CRA over lists::          
+-- @end menu
 
-seq :: a -> b -> b
-seq x y = y
+-- @node Imports, Auxiliary functions, Chinese Remainder Algorithm, Chinese Remainder Algorithm
+-- @section Imports
 
-_parGlobal_ :: Int# -> a -> b -> b
-_parGlobal_ n x y = y
+import ModArithm (modHom, modDif, modProd, modInv) 
 
-_parLocal_ :: Int# -> a -> b -> b
-_parLocal_ n x y = y
+#if defined(STRATEGIES)
+import Control.Parallel
+import Control.Parallel.Strategies
 #endif
 
-#ifdef GUM
--- GUM has only par 
-
-import {-fool mkdependHS; ToDo: rm-}
-       Parallel
-
-_parGlobal_ :: Int# -> a -> b -> b
-_parGlobal_ n x y = par x y
-
-_parLocal_ :: Int# -> a -> b -> b
-_parLocal_ n x y = par x y
-#endif
+-- @node Auxiliary functions, Basic binary CRA operation, Imports, Chinese Remainder Algorithm
+-- @section Auxiliary functions
 
-isPrime :: Integer -> Bool
+--isPrime :: Integer -> Bool
 isPrime 2 = True
 isPrime n 
  | even n    = False
  | otherwise = isPrime' n 3
 
-isPrime' :: Integer -> Integer -> Bool
-
+--isPrime' :: Integer -> Integer -> Bool
 isPrime' n l1     | n < l1*l1        = True
                   | n `mod` l1 == 0  = False
                   | otherwise        = isPrime' n (l1+2)
 
+-- @node Basic binary CRA operation, CRA over lists, Auxiliary functions, Chinese Remainder Algorithm
+-- @section Basic binary CRA operation
+#if 0 && defined(STRATEGIES)
 
 par_binCRA :: Integer -> Integer -> Integer -> Integer -> Integer -> Integer 
 
 par_binCRA m1 m2 inv a1 a2 = 
                         (if d == 0 then a1
                                     else a)
+                         `sparking` rnf a
                          where ab = modHom m2 a1
-                               d  = _parGlobal_ 0# a
-                                       modDif m2 a2 a1
+                               d  = modDif m2 a2 a1
                                b  = modProd m2 d inv
                                b0 = if (b+b>m2) then b-m2
                                                     else b
-                               a  = m1*b0+a1             {- not HWL version -}
-
+                               a  = m1*b0+a1
+#endif
 
 binCRA :: Integer -> Integer -> Integer -> Integer -> Integer -> Integer 
 
@@ -102,25 +71,34 @@ binCRA m1 m2 inv a1 a2 = (if d == 0 then a1
                                                 else b
                                a  =  m1*b0+a1             {- not HWL version -}
 
-s1IMCRA :: [Integer] -> [Integer] -> (Integer,Integer)
+-- @node CRA over lists,  , Basic binary CRA operation, Chinese Remainder Algorithm
+-- @section CRA over lists
 
-s1IMCRA (m:ms) (r:rs) = foldl 
-                        (\ (m0,r0) (m1,r1) -> (m0*m1,binCRA m0 m1 (modInv m1 m0) r0 r1) )
-                        (m,r) (zip ms rs)
+-- @node Sequential,  , CRA over lists, CRA over lists
+-- @subsection Sequential
+
+seq_list_CRA0 :: [Integer] -> [Integer] -> (Integer,Integer)
+
+seq_list_CRA0 (m:ms) (r:rs) = 
+  foldl (\ (m0,r0) (m1,r1) -> (m0*m1, binCRA m0 m1 (modInv m1 m0) r0 r1) )
+        (m,r) 
+        (zip ms rs)
 
------------------------------------------------------------------------------
--- Same as s1IMCRA but applicable for infinite lists
------------------------------------------------------------------------------
 
-s2IMCRA :: Integer -> [Integer] -> [Integer] -> [Integer] -> (Integer,Integer)
-s2IMCRA mBound (m:ms) (r:rs) (d:ds)  = s2IMCRA' mBound ms rs ds m r 
+-- Same as above but applicable for infinite lists
 
-s2IMCRA' :: Integer -> [Integer] -> [Integer] -> [Integer] -> Integer -> Integer ->
+seq_list_CRA :: Integer -> [Integer] -> [Integer] -> [Integer] -> (Integer,Integer)
+seq_list_CRA mBound (m:ms) (r:rs) (d:ds)  = seq_list_CRA' mBound ms rs ds m r 
+
+seq_list_CRA' :: Integer -> [Integer] -> [Integer] -> [Integer] -> Integer -> Integer ->
            (Integer,Integer)
-s2IMCRA' mBound (m:ms) (r:rs) (d:ds) accM accR 
+seq_list_CRA' mBound (m:ms) (r:rs) (d:ds) accM accR 
   | accM > mBound = (accM, accR)
-  | d == 0     = s2IMCRA' mBound ms rs ds accM accR
-  | otherwise  = s2IMCRA' mBound ms rs ds (m*accM) (binCRA accM m (modInv m accM) accR r)
+  | d == 0               = seq_list_CRA' mBound ms rs ds accM accR
+  | otherwise            = -- trace ("seq_list_CRA': (m,accM) = " ++ (show (m,accM))) $
+                         seq_list_CRA' mBound ms rs ds (m*accM) 
+                                       (my_cra accM m (modInv m accM) accR r)
+                         where my_cra = binCRA {- or par_binCRA -}
 
 -----------------------------------------------------------------------------
 -- Note:
@@ -135,102 +113,55 @@ s2IMCRA' mBound (m:ms) (r:rs) (d:ds) accM accR
 -- Note: ms as ds are infinte lists 
 
 tree_IMCRA0 :: Integer -> [Integer] -> [Integer] -> [Integer] -> 
-               (Integer,Integer)
+              (Integer,Integer)
 
 
 tree_IMCRA0 n ms as ds = 
        let
-         res@(m, a, fails) = {-_parGlobal_ 11# (forcelist ms')
-                               (_parGlobal_ 12# (forcelist as')
-                                 (_parGlobal_ 13# (forcelist ds') -}
-                                    tree_IMCRA0' ms' as' ds' -- ))
-                             where ms' = take n ms 
-                                   as' = take n as 
-                                   ds' = take n ds
-         handle_fails n m a (m1:ms) (a1:as) (d1:ds) =
-            -- _parGlobal_ 46# m $
-            -- _parGlobal_ 47# a $
-                        if n == 0 then (m, a)
-                        else if d1 == 0 then (handle_fails n m a ms as ds)
-                        else handle_fails (n-1) m' a' ms as ds
-                             where
-                                       m'  = m * m1
-                                       a'  = par_binCRA m m1 inv a a1
-                                       inv = modInv m1 m
-       in
-         handle_fails fails m a ms as ds
+         res@(m, a, fails) = {-parGlobal 11 11 (forcelist ms')
+                                (parGlobal 12 12 (forcelist as')
+                                 (parGlobal 13 13 (forcelist ds') -}
+                             tree_IMCRA0' ms' as' ds' -- ))
+                             where ms' = take (fromInteger n) ms 
+                                   as' = take (fromInteger n) as 
+                                   ds' = take (fromInteger n) ds
+         handle_fails :: Integer -> Integer -> Integer -> 
+                         [Integer] -> [Integer] -> [Integer] -> (Integer, Integer)
+         handle_fails n m a (m1:ms) (a1:as) (d1:ds) =
+            -- parGlobal 46 46 1 0 m $
+            -- parGlobal 47 47 1 0 a $
+               (if n == 0 
+                  then (m, a)
+                  else if d1 == 0 
+                         then (handle_fails n m a ms as ds)
+                         else handle_fails (n-1) m' a' ms as ds
+               )
+#if defined(STRATEGIES)
+               `using` \ x -> do par m (pseq a (return x)) --{ m' <- rpar m; a' <- rseq a; return x }
+#endif
+               where
+                     m'  = m * m1
+                     a'  = {-par_binCRA-} binCRA m m1 inv a a1
+                     inv = modInv m1 m
+       in
+         handle_fails fails m a ms as ds
 
 tree_IMCRA0' [m] [a] [0] = (1, 1, 1)   -- FAIL due to unlucky prime
 tree_IMCRA0' [m] [a] [_] = (m, a, 0)   -- normal case
 tree_IMCRA0' ms as ds =
-       let 
-         n = length ms
-         (left_ms, right_ms) = splitAt (n `div` 2) ms
-         (left_as, right_as) = splitAt (n `div` 2) as
-         (left_ds, right_ds) = splitAt (n `div` 2) ds
-         left@(left_M, left_CRA, left_fails)  = 
-               tree_IMCRA0' left_ms left_as left_ds
-         right@(right_M, right_CRA, right_fails)  = 
-               tree_IMCRA0' right_ms right_as right_ds
-         inv = modInv right_M left_M
-         cra = par_binCRA left_M right_M inv left_CRA right_CRA
-       in
-#if defined(TU_CRA)
-         _parGlobal_ 5# left
-          (_parGlobal_ 6# right
-            (_parLocal_ 9#  inv
-            --seq cra
-               (left_M * right_M, 
-                cra,
-                left_fails + right_fails ) -- ))
-#elif defined(COMP_CRA)
-         _parGlobal_ 5# left_CRA
-          (_parGlobal_ 6# left_M
-           (_parGlobal_ 7#  right_CRA
-             (_parGlobal_ 8#  right_M
-               (_parGlobal_ 9#  inv
-                 (left_M * right_M, 
-                  cra,
-                  left_fails + right_fails)))))
-#elif defined(ASYM_CRA)
-         _parGlobal_ 5# left
-          (_parGlobal_ 7# right
-            (_parGlobal_ 9#  inv
-               (seq cra
-                 (left_M * right_M, 
-                  cra,
-                  left_fails + right_fails)) ) )
-#endif
-
--- +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
--- Experimental pat starts here
--- +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
--- This is a brain damaged version of doing a binary-tree CRA over lists.
-
-tree_IMCRA :: [Integer] -> [Integer] -> (Integer,Integer)
-
-tree_IMCRA ms as = _parGlobal_ 0# inv 
-                                 (m1*m2, par_binCRA m1 m2 inv a1 a2)
-                  where --( [(m1,a1),(m2,a2)] : _ ) =
-                        --(m1,a1) = (head foo) !! 0
-                        --(m2,a2) = (head foo) !! 1
-                        -- foo :: Integer -- (Integer,Integer) -> (Integer,Integer) -> Boolean
-                        [(m1, a1), (m2, a2)] =   
-                               head (dropWhile (\ l -> (length l) /= 2)
-                                         (iterate tree_IMCRA' (zip ms as)))
-
-                        inv = modInv m2 m1
-                        -- cra = par_binCRA m1 m2 inv a1 a2
-
-tree_IMCRA' :: [(Integer,Integer)] -> [(Integer,Integer)]
-
-tree_IMCRA' [] = []
-tree_IMCRA' [(m1,a1)] = [(m1,a1)]
-tree_IMCRA' ((m1,a1) : (m2,a2) : rest) = 
-   _parGlobal_ 11# cra ({-_parGlobal_ 12# (forcelist rest')-} (cra : rest') )
-   where cra = (m1*m2, binCRA m1 m2 (modInv m2 m1) a1 a2)
-        rest' = tree_IMCRA' rest
+       let 
+         n = length ms
+         (left_ms, right_ms) = splitAt (n `div` 2) ms
+         (left_as, right_as) = splitAt (n `div` 2) as
+         (left_ds, right_ds) = splitAt (n `div` 2) ds
+         left@(left_M, left_CRA, left_fails)  = 
+              tree_IMCRA0' left_ms left_as left_ds
+         right@(right_M, right_CRA, right_fails)  = 
+              tree_IMCRA0' right_ms right_as right_ds
+         inv = modInv right_M left_M
+         cra = {-par_binCRA-} binCRA left_M right_M inv left_CRA right_CRA
+       in
+       (left_M * right_M, cra, left_fails + right_fails )
 
 -- IMCRA with list of products of the modules as additional argument i.e.
 -- in s2IMCRA m M r the following holds M_i = m_1 * ... * m_{i-1}
@@ -251,72 +182,3 @@ s3IMCRA (m:ms) pp (r:rs) = foldl
                         (\ (m0,r0) (m1,r1,p) -> (p,binCRA m0 m1 (modInv m1 m0) r0 r1) )
                         (m,r) (zip3 ms rs pp)
 
------------------------------------------------------------------------------
--- Old cheating version
------------------------------------------------------------------------------
-
-#if 0   /* defined(NO_FILTER) */
-
--- Note: ms as are finte lists 
---      no ds needed as no test (== 0) is performed
-
-tree_IMCRA0 :: [Integer] -> [Integer] -> 
-               (Integer,Integer)
-
--- tree_IMCRA0 [] [] = (1, 1)
-tree_IMCRA0 [m] [a] = (m, a)
-tree_IMCRA0 ms as =
-       let 
-               n = length ms
-               (left_ms, right_ms) = splitAt (n `div` 2) ms
-               (left_as, right_as) = splitAt (n `div` 2) as
-               left@(left_M, left_CRA)  = tree_IMCRA0 left_ms left_as
-               right@(right_M, right_CRA)  = tree_IMCRA0 right_ms right_as
-               inv = modInv right_M left_M
-       in
-#if defined(TU_CRA)
-               _parGlobal_ 5# left
-                (_parGlobal_ 6# right
-                   (_parLocal_ 9#  inv
-                       (left_M * right_M, 
-                        par_binCRA left_M right_M inv left_CRA right_CRA)))
-#elif defined(COMP_CRA)
-               _parGlobal_ 5# left_CRA
-                (_parGlobal_ 6# left_M
-                 (_parGlobal_ 7#  right_CRA
-                  (_parGlobal_ 8#  right_M
-                   (_parGlobal_ 9#  inv
-                       (left_M * right_M, 
-                        par_binCRA left_M right_M inv left_CRA right_CRA)))))
-#elif defined(ASYM_CRA)
-               _parGlobal_ 5# left
-                --(_parLocal_ 6# right
-                   (_parGlobal_ 9#  inv
-                       (left_M * right_M, 
-                        par_binCRA left_M right_M inv left_CRA right_CRA)) --)
-#endif 
-        
-{-
-                               --(left_CRA + right_CRA, 
-                               --binCRA left_M right_M inv left_CRA right_CRA) )
--}
-
-{- Test check for equal length lists. 
-       | length ms /= length as = error "tree_IMCRA0: different lengths of input lists"
--}
-
-#endif
-
--- Old fail handler
-         {-
-         handle_fails 0 m a _       _       _      = (m, a)
-         handle_fails n m a (m1:ms) (a1:as) (0:ds) = 
-               handle_fails n m a ms as ds
-         handle_fails n m a (m1:ms) (a1:as) (d1:ds) =
-               handle_fails (n-1) m' a' ms as ds
-               where
-                 m'  = m * m1
-                 a'  = par_binCRA m m1 inv a a1
-                 inv = modInv m1 m
-
-         -}
diff --git a/parallel/linsolv/LinSolv.hs b/parallel/linsolv/LinSolv.hs
deleted file mode 100644 (file)
index 3ac1176..0000000
+++ /dev/null
@@ -1,346 +0,0 @@
-{-
-$Id: LinSolv.hs,v 1.1 1997/07/31 22:15:08 sof Exp $
-
-This is revision: $Revision: 1.1 $
-
-Module for computing the solution of a system of linear system of equations
-in several variables.
-The rhs of the system of equations is the first parameter (given as a 
-square matrix). The lhs of the system is the second parameter (given as a
-vector).
-
-ChangeLog:
-$Log: LinSolv.hs,v $
-Revision 1.1  1997/07/31 22:15:08  sof
-Renamed LinSolv-par to LinSolv
-
-Revision 1.1  1996/01/08 20:07:34  partain
-Initial revision
-
---# Revision 1.4  1994/12/13  19:49:38  hwloidl
---# Version based on infinite lists.
---# Has seq steps in all variants (even without det check)
---# -DNOFILTER uses finite part of inifinte xList and avoids any det check
---# -DSEQ_CRA uses a fold-based CRA instead of the tree-based one (both over
---#           infinte lists)
---#
---# Revision 1.3  1994/11/23  01:04:43  hwloidl
---# Test version of fwd mapping and hom sol (no CRA in this version)
---#  New: prime list computed in par and packet forced
---#       modDet computed in par to hom sol
---#       No modDet == 0 test in gen_xList to avoid seq in that step ->
---#       count and filter are used to detect number of necessary hom ims
---#       xList is now again an infinite list (it is not forced in this version)
---#  Note: The result of linsolv just triggers computation up to hom sols but
---#        not until CRAs. The value is not important, just the tiggering
---#  NO_FILTER flag: if set any modDet == 0 test is avoided (i.e. no filtering
---#                  either; this doesn't always yield a correct solution but
---#                  shows how this test seq the computation
---#
---# Revision 1.2  1994/11/19  21:49:38  hwloidl
---# Compute list of primes in parallel to fwd mapping and sol in hom ims
---#
---# Revision 1.1  1994/11/19  02:00:05  hwloidl
---# Initial revision
---#
-
----------------------------------------------------------------------------- -}
-
-module LinSolv(linSolv) where
-
-import ModArithm {- (modHom, Hom(hom))-}
-
-import Matrix (SqMatrix, Vector, vector, matBounds,
-               determinant, transp, replaceColumn, size,
-               maxElem, maxElemVec, scalarMult, vecScalarQuot,
-               matGcd, vecGcd, matHom, vecHom) 
-
-import CRA (s1IMCRA,s2IMCRA,tree_IMCRA0)
-  
-#ifndef SEQ
-
-import ParForce 
-import {-fool mkdependHS; ToDo: rm-}
-       Parallel
-
-#else
-
-parmap = map
-
-par :: a -> b -> b
-par x y = y
-
-seq :: a -> b -> b
-seq x y = y
-
-_parGlobal_ :: Int# -> a -> b -> b
-_parGlobal_ n x y = y
-
-_parLocal_ :: Int# -> a -> b -> b
-_parLocal_ n x y = y
-#endif
-
-#ifdef GUM
--- GUM has only par 
-_parGlobal_ :: Int# -> a -> b -> b
-_parGlobal_ n x y = par x y
-
-_parLocal_ :: Int# -> a -> b -> b
-_parLocal_ n x y = par x y
-#endif
-
--- ---------------------------------------------------------------------------
--- Auxiliary functions
--- ---------------------------------------------------------------------------
-
-primes :: [Integer]
-
-primes = 2 : oddPrimesFrom 3
-         where oddPrimesFrom n | isPrime n primes = n : oddPrimesFrom (n+2)
-                               | otherwise        =     oddPrimesFrom (n+2)
-
-isPrime :: Integer -> [Integer] -> Bool
-
-isPrime n (l1:ls) | n < l1*l1        = True
-                  | n `mod` l1 == 0  = False
-                  | otherwise        = isPrime n ls
-
-projection :: Integer -> [[Integer]] -> [Integer]
-
-projection i [] = [] 
-projection i (ls1:lss) = ls1!!(i-1) : projection i lss
-fact :: Integer -> Integer
-
-fact 0   = 1
-fact (n+1) = (n+1) * fact n
-
--- ---------------------------------------------------------------------------
--- Function for solving a system of linear equations
--- ---------------------------------------------------------------------------
-
--- linSolv :: (Integral a) => SqMatrix a -> Vector a -> (a,a,Vector a)   
-
-linSolv a b = 
-  let 
-   {- Step1: -}
-    n = toInteger (size a) 
-    s = max (maxElem a) (maxElemVec b)
-    pBound = 2 * s^n * fact n
-
-   {- Step2: -}
-    xList :: [[Integer]]
-    xList = _parGlobal_ 19# (forcelist1 100 primes) 
-               _parGlobal_ 18# noOfPrimes 
-                   (gen_xList primes)
-
-    gen_xList :: [Integer] -> [[Integer]]
-    gen_xList (p:ps)
-                   = let 
-                       b0 = vecHom p b
-                       a0 = matHom p a
-                       modDet = modHom p (determinant a0)
-
-                       {-  Parallelism inside a homomorphic solution -}
-
-                       homSol = _parGlobal_ 8# modDet 
-                                               (p : modDet : pmx)
-
-                       pmx = parmap ( \ j ->
-                                       let
-                                               a1 = replaceColumn j a0 b0
-                                       in
-                                               modHom p (determinant a1) )
-                                    [jLo..jHi] 
-                       ((iLo,jLo),(iHi,jHi)) = matBounds a
-
-                       restList = gen_xList ps 
-                     in
-                       _parGlobal_ 3#  homSol 
-                                       (homSol : restList)
-
-   {- Step3: -}
-
-   {- s1IMCRA variant -}
-
-    count (l1:ls) i accum = if accum>pBound
-                              then i
-                              else count ls (i+1) (l1*accum)
-
-    -- noOfPrimes = count (projection 1 xList) 0 1 
-
-    -- Just for TESTING:
-    --  noOfPrimes_0 = count primes 0 1 
-    --  xList_0 = take noOfPrimes xList
-
-#if defined(NO_FILTER)
-    {-
-    noOfPrimes = count (projection 1 xList) 0 1 
-    xList'' = take noOfPrimes xList
-
-    luckyPrimes :: [Integer]
-    luckyPrimes = projection 1 xList''
-    -}
-
-    -- 1. The lists passed to CRA are finite
-    -- 2. No test for det == 0 is ever performed
-
-    noOfPrimes = count (projection 1 xList) 0 1
-
-    primeList = take noOfPrimes (projection 1 xList)
-    detList = take noOfPrimes (projection 2 xList)
-
-    -- Force the lists?
-
-    det = snd (tree_IMCRA0 primeList detList) 
-
-    x_i i = snd (tree_IMCRA0 primeList x_i_List)
-           where
-                 x_i_List = take noOfPrimes (projection (i+2) xList)
-
-    x = vector (parmap x_i [1..n])
-
-#else
-    -- No more cheating ...
-
-    {- The naive, inefficient version
-    xList' = filter (\ (x1:x2:xs) -> x2 /= 0) xList
-    noOfPrimes = count (projection 1 xList') 0 1
-    xList'' = take noOfPrimes xList'
-    -}
-
-    noOfPrimes = count (projection 1 xList) 0 1
-    primeList = projection 1 xList
-    detList = projection 2 xList     -- list of dets in hom ims if det==0
-                                    -- we will throw the result away in 
-                                    -- tree_IMCRA0 => needed in there
-#if defined(SEQ_CRA)
-    det = a
-         where (_, a) = s2IMCRA pBound {-noOfPrimes-} primes detList detList
-
-    x_i i = a
-           where (_, a) = s2IMCRA pBound {-noOfPrimes-} primes x_i_List detList
-                 x_i_List = projection (i+2) xList
-#else
-    det = snd (tree_IMCRA0 noOfPrimes primeList detList detList)
-
-    x_i i = snd (tree_IMCRA0 noOfPrimes primeList x_i_List detList)
-           where
-                  x_i_List = projection (i+2) xList
-#endif
-
-    x = vector (parmap x_i [1..n])
-#endif
-
-   {- Step4: -}
-    gcdVec = vecGcd x
-    gcdVal = seq det 
-                  (gcd gcdVec det)
-
-    -- Solution i.e. a * (res_a/res_b) * res_x = b
-    res_a :: Integer
-    res_a = gcdVec `div` gcdVal
-    res_b :: Integer
-    res_b = det `div` gcdVal
-    res_x :: Vector Integer
-    res_x = vecScalarQuot gcdVec x
-  in 
-    -- (det, gcdVec, gcdVal, x)   -- Just for testing
-    -- (pBound , length xList, last xList) --, sum (map sum xList))
-    -- (noOfPrimes, pBound, sum (parmap sum xList)) --, head xList, last xList) --, sum (map sum xList))
-    --(noOfPrimes, pBound, sum (parmap sum xList''), det, x)
-    --(noOfPrimes, pBound, det)
-    --(take 20 xList)
-    --(res_x, res_a, res_b)
-    (vecScalarQuot gcdVec x, gcdVec `div` gcdVal, det `div` gcdVal)   
-    --(x, 99, 99)   
-    --(noOfPrimes, x)
-
--- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
-   {- s3IMCRA variant:  
-
-    --luckyPrimes :: [Integer]
-    --luckyPrimes = take noOfPrimes (projection 1 xList)
-
-    (noOfPrimes, luckyPrimes, primeProds) = count (projection 1 xList) 0 1 
-        where count :: [Integer] -> Integer -> Integer -> (Integer,[Integer],[Integer])
-              count (l1:ls) i accum =  
-                    let 
-                         new_accum = l1*accum
-                         res@(n,lP,pP) = count ls (i+1) new_accum
-                       in 
-                         if accum>pBound
-                           then (i, [], [])
-                           else (n, l1:lP, new_accum:pP)
-
-    det = snd( s3IMCRA luckyPrimes
-                       (take noOfPrimes (projection 2 xList))
-                       primeProds )
-
-    x_i i = snd( s3IMCRA luckyPrimes
-                         (take noOfPrimes (projection (i+2) xList))
-                         primeProds )
-
-    x = vector [ x_i i | i <- [1..n] ]              {- should be a vector -}
-   -}
-
-{-
-main =
- let 
-   a1 = sqMatrix (listArray ((1,1),(3,3)) [ 1, 1, 1, 
-                                            0, 1, 2,
-                                            2, 2, 1] )
-   l1 = [6, 8, 9]
-   b1 = vector l1
-   x = linSolv a1 b1
- in
-   appendChan stdout "Testing linSolv: \n" abort $
-   appendChan stdout "Solving a1*x=b1 with the following a1 and b1: \n a1 = " abort $
-   appendChan stdout (show a1) abort $
-   appendChan stdout "\nb1 = " abort $
-   appendChan stdout (show b1) abort $
-   appendChan stdout "\nx = " abort $
-   appendChan stdout (show x) abort $
-   done 
--}
-   {- Seq version: 
-    xList :: [[Integer]]
-    xList = [ let
-                 b0 = vecHom p b
-               in 
-                  p : modDet :
-                  [ let
-                      a1 = replaceColumn j a0 b0
-                    in
-                      modHom p (determinant a1)
-                  | j <- [jLo..jHi] ]
-             | p <- primes, a0 <- [matHom p a], 
-               modDet <- [modHom p (determinant a0)],
-               modDet /= 0 ]
-             where
-               ((iLo,jLo),(iHi,jHi)) = matBounds a
-   -}
-{- 
-  Old version of the body of gen_xList:
-                       if modDet == 0 
-                         then restList
-                         else   _parGlobal_ 3# homSol 
-                                       (_parGlobal_ 4# (forcelist restList)
-                                                       (homSol : restList) )
-
--}
-
-                       {-
-                       Seq. hom sol
-
-                       homSol = p : modDet :
-                                [ let
-                                       a1 = replaceColumn j a0 b0
-                                  in
-                                       modHom p (determinant a1)
-                                | j <- [jLo..jHi] ]
-                                where ((iLo,jLo),(iHi,jHi)) = matBounds a
-                       -}
index 54c8e78..4fe8994 100644 (file)
-{- Main module for testing linSolv.hs, the solution of a system of 
-   linear equations by using the modular method.                             -}
+{- 
+Time-stamp: <2010-11-03 09:26:39 simonmar>
+$Id: testLinSolv.hs,v 1.6.2.5 2002/06/21 02:39:07 hwloidl Exp $
+
+Main module for testing linSolv.hs, the solution of a system of 
+linear equations by using the modular method.                             -}
+
+-- @node Top, Modular Arithmetic, (dir), (dir)
+-- @top A Linear System Solver
+
+-- @node Test wrapper for LinSolv,  , Top, Top
+-- @chapter Test wrapper for LinSolv
 
 module Main where
 
+--Other  modules:
+-- @menu
+-- * Modular Arithmetic::      
+-- * ADT Matrix::              
+-- * LinSolv top level function::  
+-- @end menu
+
 -- ----------------------------------------------------------------------------
+-- @node Imports, Functions for I/O, Test wrapper for LinSolv, Test wrapper for LinSolv
+-- @section Imports
+--
 -- Using modular method for linear system solving
 -- ----------------------------------------------------------------------------
 
+import System.Environment(getArgs)
+
+-- @include ModArithm.hs
 import ModArithm
 
-import Matrix (SqMatrix, Vector, vector,
-              lolSqMatrix,listSqMatrix,maxElemVec,jaH)
+-- @include Matrix-list.hs
+import Matrix (SqMatrix, Vector, 
+               vector,
+              lolSqMatrix,listSqMatrix,maxElemVec)
 
-import ParForce 
+#if defined(STRATEGIES)
+import Control.Parallel
+import Control.Parallel.Strategies
+#endif
+import Control.DeepSeq
+
+-- @include LinSolv-par.hs
+import LinSolvS (linSolv)
 
-import LinSolv (linSolv)
+-- fixed input data
+import Inputs8
 
+-- get_input :: (Integral a, NFData a) =>
+--              Int -> [(Int, SqMatrix a, Vector a)] -> (SqMatrix a, Vector a)
+get_input n = case (lookup (fromIntegral n) all_inputs) of
+                Just x -> x
+                Nothing -> error $ "Unknown input " ++ (show n) ++ ". Possible inputs are " ++ (show (map fst all_inputs))
+              
 -- ----------------------------------------------------------------------------
+-- @node Functions for I/O, Main fct, Imports, Test wrapper for LinSolv
+-- @section Functions for I/O
+--
 -- Functions for I/O
 -- ----------------------------------------------------------------------------
 
-seq :: [String] -> IO ()
-
-seq []      = return ()
-seq (s1:ss) = putStr s1 >> seq ss
-
-{-
-getCpuTime fail succ resps = 
-  GetCpuTime : dblDispatch fail succ resps
-
-getTime fail succ resps = 
-  GetTime : dblDispatch fail succ resps
-
-dblDispatch fail succ ~(resp:resps) = 
-  case resp of Dbl t       -> succ t resps
-               Failure msg -> fail msg resps
--}
-
 -- Compute one number out of result; JUST FOR TESTING
 -- This demands the computation of the full result and is fast to print
 -- i.e. not much IO overhead in timing.
 
-compact :: (Integral a) => (Vector a, a, a) -> [a]
-compact (x', a, b) = if a == 0 
-                    then if b == 0 then bonzo x
-                                   else 1 : bonzo x
-                    else if b == 1 then bonzo x
-                                   else 1 : bonzo x
-                    where bonzo x = parmap (\ y -> if y==0 then 0 else 1) x
+--compact :: (Integral a) => (Vector a, a, a, a) -> a
+{-
+compact (x', a, b, _) = if a == 0 
+                      then if b == 0 then x
+                                     else 1 +  x
+                      else if b == 1 then  x
+                                     else 1 + x
+                    where -- bonzo x = map (\ y -> if y==0 then 0 else 1) x
                           x = jaH x'
+-}
 
 -- ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
-main = 
-       let
-         a1 = listSqMatrix ((1,1),(3,3)) [ 1, 1, 1, 
-                                                  0, 1, 2,
-                                                  2, 2, 1 ] 
-         b1 = vector [ 6, 8, 9 ]
-         x1 = linSolv a1 b1            {- x1 = 1/-1 * [ -1, -2, -3] -}
-
-         a2 = listSqMatrix ((1,1),(3,3)) [  1, -1,  4, 
-                                                  -3,  0,  1,
-                                                   2, -2,  1 ] 
-         b2 = vector [ 12, -7, 10 ]
-         x2 = linSolv a2 b2            {- x2 = 1/1 * [3, -1, 2] -}
-
-         a3 = listSqMatrix ((1,1),(2,2)) [  1, 3, 
-                                                  -2, 2 ] 
-         b3 = vector [ 11, 10 ]
-         x3 = linSolv a3 b3            {- x3 = 1/1 * [ -1, 4 ] -}
-
-         a4 = listSqMatrix ((1,1),(4,4)) [  1,  1,  0,  2,
-                                                   2, -1,  1, -1, 
-                                                   3,  2, -2,  4,
-                                                  -1,  2,  3,  0 ] 
-         b4 = vector [ 8, 1, 13, 3 ]
-         x4 = linSolv a4 b4            {- x4 = 1/-1 * [ -1, 1, -2, -4 ] -}
-
-         a5 = listSqMatrix ((1,1),(3,3)) [  0,  0, -2,
-                                                   7,  0,  5,
-                                                  -3,  4, -6 ] 
-         b5 = vector [ 7, 7, 7 ]
-         x5 = linSolv a5 b5            {- x5 = 7/-8 * [ -4, 1, 4 ] -}
-
-         a6 = listSqMatrix ((1,1),(2,2)) [  3,  0,
-                                                  -3,  3 ] 
-         b6 = vector [ 0, 2 ]
-         x6 = linSolv a6 b6            {- x6 = 2/3 * [ 0, 1 ] -}
-
-         a7 = listSqMatrix ((1,1),(4,4)) [ -1, -1,  0,  0,
-                                                  -1,  0,  0, -1,
-                                                   0, -1, -1,  1,
-                                                  -1,  1,  0,  0 ] 
-         b7 = vector [ -1, 0, 0, 0 ]
-         x7 = linSolv a7 b7            {- x7 = 1/2 * [ 1, 1, -2, -1 ] -}
-
-         a8 = listSqMatrix ((1,1),(8,8)) 
-                       [  1,  4,  0, 11, -4,  5, 23, -8, 
-                           0,  5, 17, -3, 29, 11,  7,  4,
-                         -2, -7, -9, -1,  0,  6,  9, 21,
-                         32,  6, 11, -6,  2,  5,  9, -9, 
-                         13, 17, 29, -9, 12, 35, 38,  0,
-                          7,  5, -7,  8, -1,  0,  1,  2,
-                          2,  3, -6, 24,  7, -1,  5, 19,
-                         -3,  2,  9, 13, 25, -7, -1,  5 ] 
-
-         b8 = vector [526, 244, 274, 320, 1051, 125, 276, 24]
-         x8 = linSolv a8 b8            {- x8 =  [ 4, 9, -3, 0, 2, 1, 23, 5 ] -}
-
-         a9 = listSqMatrix ((1,1),(5,5)) 
-                       [   1,  0,  5,  9,  2,
-                           3, -5, 13, -1,  0,
-                          -7,  6,  8,  1, 23,
-                          -1,  2,  5, 17,  1,
-                           0, -5, -2, -1,  7 ] 
-        b9 = vector [115, 156, 229, 140, 37]
-        x9 = linSolv a9 b9                     {- x9 = [ 1, -3, 11, 5, 7 ] -}
-         a10 = listSqMatrix ((1,1),(6,6)) 
-                       [  0,  0,  2,  1,  2, -1,
-                          3, -5, -1, -7,  0,  2,
-                          1,  0,  3, -2,  5,  0,
-                         -3,  1, -5,  0,  0, -4,
-                          1,  2,  0, -6,  2, -1,
-                          7, -2,  5, -3,  0,  5 ] 
-        b10 = vector [14, -18, 19, 18, 3, -24]
-        x10 = linSolv a10 b10             {- x10 = [ 1, 0, -1, 2, 5, -4 ] -}
-
-         a11 = listSqMatrix ((1,1),(7,7))
-               [ 13, -6,  5,  1,  9,  11, -21,
-                  0,  2, -1, 15, -1,   7,   6,
-                 27, -5,  3,  2,  0,  -1,  -8,
-                  3,  1, 11, -4,  8,   0, -13,
-                  5,  7,  9, 10, 13, -29,   6,
-                 17, -2,  7,  1,  6,  -5,   1,
-                  4, -1,  0, 13,  6, -19,   2 ]
-        b11 = vector [ 3, 1, -4, 9, 13, -11, 5 ]
-        x11 = linSolv a11 b11  {- x11 = 1/11138030 *
-                                     [ -1361273, 24524143, -5023241, 4646680, 
-                                       -4426356, -4257689, -14542705]  -}
-
-         a_huge = lolSqMatrix ((1,1),(16,16)) [
-                    [44924, 20447, 29276, 40951, 5939, 996, 828, 62553, 63695, 30588, 25619, 32394, 63129, 63694, 39185, 40378],
-                    [42921, 62054, 23399, 29333, 8566, 30909, 30055, 64891, 54909, 49172, 46320, 32652, 45621, 17632, 9264, 19098],
-                    [25815, 52465, 63565, 65486, 55552, 52004, 44895, 34511, 49784, 4985, 10311, 27617, 31058, 34162, 35039, 26088],
-                    [49670, 59116, 17022, 4689, 38334, 62487, 56107, 45233, 46942, 37146, 55954, 35711, 22658, 2164, 44732, 55607],
-                    [44168, 5116, 18693, 40788, 34627, 7578, 35898, 40475, 19712, 4726, 7064, 19388, 23497, 2134, 14499, 2139],
-                    [60549, 55183, 26167, 20043, 33321, 13385, 32805, 27000, 37236, 6897, 22401, 11099, 60617, 34354, 53034, 8541],
-                    [59027, 54672, 7370, 12417, 4486, 30214, 29519, 22169, 32584, 46623, 62480, 25354, 36661, 30636, 41936, 25868],
-                    [40590, 22341, 11608, 19798, 34183, 55625, 44401, 33781, 60677, 27455, 62873, 39102, 23854, 19340, 33328, 18615],
-                    [65221, 56835, 15302, 13304, 16, 24870, 62514, 24506, 20878, 45057, 2315, 12301, 29937, 43638, 36656, 10205],
-                    [27639, 33459, 56786, 29138, 51908, 15590, 44807, 23451, 42840, 33062, 61254, 14966, 18778, 62308, 8203, 45007],
-                    [47620, 3963, 27701, 20749, 40909, 18793, 25356, 61272, 31210, 52127, 59875, 23758, 61105, 35731, 44582, 5403],
-                    [4923, 59571, 27361, 32663, 54364, 50740, 47073, 46706, 11430, 39028, 62126, 15938, 13854, 21280, 41020, 44250],
-                    [36495, 39319, 22790, 8282, 21049, 53476, 58482, 63690, 24856, 36757, 61169, 28765, 61383, 15468, 27518, 48825],
-                    [6670, 22264, 30113, 51909, 7221, 31731, 63234, 7272, 22321, 44164, 7022, 19185, 28804, 53810, 10030, 8397],
-                    [23372, 39852, 37439, 8528, 6944, 9497, 40252, 4427, 52306, 14132, 5411, 12954, 23639, 54223, 59029, 57407],
-                    [5848, 700, 41324, 34164, 54152, 46356, 32955, 4917, 41869, 52246, 39441, 30706, 45053, 17938, 282, 4849]
-                    ]
-
-         b_huge = vector [36495, 39319, 22790, 8282, 21049, 53476, 58482, 63690, 24856, 36757, 61169, 28765, 61383, 15468, 27518, 48825]
-
-        x_huge = linSolv a_huge b_huge
-        -- This is the preferred test setting:
-
-#if (INPUT_SIZE==1)
-        a = a4
-        b = b4
-        x = x4
-#elif (INPUT_SIZE==2)
-        a = a9
-        b = b9
-        x = x9
-#elif (INPUT_SIZE==3)
-        a = a8
-        b = b8
-        x = x8
-#elif (INPUT_SIZE==4)
-        a = a11
-        b = b11
-        x = x11
-#elif (INPUT_SIZE==9)
-        a = a_huge
-        b = b_huge
-        x = x_huge
-#else
-        a = a9
-        b = b9
-        x = x9
+-- @node Main fct,  , Functions for I/O, Test wrapper for LinSolv
+-- @section Main fct
+
+main = do
+        args <- getArgs
+        let 
+         v = read (args!!0) :: Int -- version (i.e. strategy to use)
+         n = read (args!!1) :: Int -- input
+        -- Pick input based on command line argument
+         (a, b) = get_input n
+         x = linSolv v a b
+#ifdef TEST
+       putStrLn ( case determinant a of
+                       0 -> error "determinant 0, no solution"
+                       x -> "determinant: " ++ show x)
 #endif
-       in  
-         seq [ "Testing linSolv with various matrices:\n",
-               "\nSolving a*x=b for a = ",
-               show a,
-               "b = ",
-               show b,
-               "\nSolution: ",
-               show x,
-               "\n"
-             ]
-
+# if defined(OUTPUT)
+       putStr ("Testing linSolv with various matrices:\n" ++
+               "\nSolving a*x=b for a = " ++
+               show a ++
+               "b = " ++ 
+               show b ++ 
+               "\nSolution: " ++
+               show x ++
+               "\n")
+# else
+        putStrLn (x `deepseq` "done" )
+#endif /* OUTPUT */
+        -- putStr ( show (compact x) ++ "done")
index b829568..4c4a1fe 100644 (file)
@@ -1,7 +1,11 @@
 TOP = ../..
 include $(TOP)/mk/boilerplate.mk
 
-SRC_HC_OPTS += -cpp
+SRC_HC_OPTS += -cpp -DSTRATEGIES -package random -package parallel
+
+# 28 = version
+# 83 = input
+PROG_ARGS = 28 83
 
 include $(TOP)/mk/target.mk
 
index 8158aeb..d0763ed 100644 (file)
@@ -1,53 +1,73 @@
-{- 
-$Id: Matrix.hs,v 1.1 1997/07/31 22:14:00 sof Exp $
-
-This is revision: $Revision: 1.1 $
-
-Data Encapsulation of the ADT Matrix.
-Internal representation is a list of lists.
-
-Changelog:
-$Log: Matrix.hs,v $
-Revision 1.1  1997/07/31 22:14:00  sof
-Renamed Matrix-list to Matrix
-
-Revision 1.2  1996/07/25 21:20:27  partain
-Bulk of final changes for 2.01
-
-Revision 1.1  1996/01/08 20:07:34  partain
-Initial revision
-
---# Revision 1.3  1994/11/23  01:06:18  hwloidl
---# Version for testing fwd mapping and hom sol.
---# Unchanged from previous version
---#
---# Revision 1.2  1994/11/19  21:50:15  hwloidl
---# *** empty log message ***
---#
---# Revision 1.1  1994/11/19  02:00:05  hwloidl
---# Initial revision
---#
---# Revision 1.1  1994/11/19  02:00:05  hwloidl
---# Initial revision
---#
---------------------------------------------------------------------------  -}
+-- Time-stamp: <2010-11-03 09:27:34 simonmar>
+-- $Id: Matrix.hs,v 1.4.2.5 2002/06/15 01:34:29 hwloidl Exp $
+
+-- Data Encapsulation of the ADT Matrix.
+-- Internal representation is a list of lists.
+--
+-- LinSolv remark: default determinant is parallel (d&c over 1st line)
+-----------------------------------------------------------------------------
+
+-- @node ADT Matrix, , ,
+-- @chapter ADT Matrix
 
 module Matrix(SqMatrix, Vector, {- MatBounds, -}
-              listSqMatrix, lolSqMatrix, vector,
+              (!!-), (!-), sqMatrix, 
+              vecBounds, matBounds, vecCont, matCont,
+              listSqMatrix, lolSqMatrix, unSqMatrix, vector, unvector,
               determinant, transp, replaceColumn, size,
               maxElem, maxElemVec, scalarMult, vecScalarQuot, 
-              matGcd, vecGcd, matHom, vecHom, matBounds,
-             matMult, jaH) 
+              matGcd, vecGcd, matHom, vecHom, matBounds, matCont,
+              matMult, matCompact
+              ) 
               {- showsMatrix, matEqual, matSum, matMult)
                  matSum',matSum'',showIt,makeUnique) -}                   where
 
+-- @menu
+-- * Imports::                 
+-- * Data Types::              
+-- * Constants::               
+-- * Aux functions::           
+-- * Data type constructors::  
+-- * H.o. fcts::               
+-- * Misc operations::         
+-- * Arithmetic Operations::   
+-- * I/O Operations::          
+-- * Instances::               
+-- @end menu
+
+-- @node Imports, Data Types, ADT Matrix, ADT Matrix
+-- @section Imports
+
+import Data.List(transpose)
+import Data.Array
+
 import ModArithm ({- Hom(hom), -} modHom)
+-- only needed if we use array based LU Decomp later
 
-#ifndef SEQ
-import ParForce 
+#if defined(STRATEGIES)
+import Control.Parallel.Strategies
 #endif
+import Control.DeepSeq
+
+infixl 5 !!-
+infixl 5 !-
+
+m!!-(i,j) = (m'!!i')!!j'
+            where bds@((rl,cl),(rh,ch)) = matBounds m
+                  i' = i - rl
+                  j' = j - cl
+                  m' = matCont m
+
+v!-i = v'!!i'
+       where bds@(rl,rh) = vecBounds v
+             i' = i - rl
+             v' = vecCont v
+
 
 -- ----------------------------------------------------------------------------
+-- @node Data Types, Constants, Imports, ADT Matrix
+-- @section Data Types
+--
 -- Data Type definitions
 -- ----------------------------------------------------------------------------
 
@@ -57,107 +77,128 @@ data  (Integral a)  =>  Vector a   = VectorC VecBounds [a]
 type  MatBounds = ((Int,Int),(Int,Int))
 type  VecBounds = ((Int),(Int))
 
+instance (NFData a, Integral a) => NFData (Vector a) where
+  -- rnf x@(VectorC b l) = rnf b >> rnf l >>  return x
+  rnf (VectorC b l) = rnf b `seq` rnf l
+
+instance (NFData a, Integral a) => NFData (SqMatrix a) where
+  -- rnf x@(SqMatrixC b m) = rnf b >> rnf m >> return x
+  rnf (SqMatrixC b m) = rnf b `seq` rnf m
+
+-- ----------------------------------------------------------------------------
+-- @node Aux functions, Data type constructors, Constants, ADT Matrix
+-- @section Aux functions
 -- ----------------------------------------------------------------------------
 
 lol :: (Integral a) => Int -> [a] -> [[a]]
-
 lol _ [] = []
 lol n l = let 
            (line, rest) = splitAt n l
           in
            line : (lol n rest)
 
-#ifndef SEQ
-
-mat_map = parmap
-
-listCompwiseComp :: (a -> b -> c) -> [a] -> [b] -> [c]
-
-listCompwiseComp = par_zipWith 
-                                       {- parmap f' (parzip l l')
-                                       where f' (a,b) = a `f` b -}
-#else
 
 mat_map = map
 
 listCompwiseComp :: (a -> b -> c) -> [a] -> [b] -> [c]
-
 listCompwiseComp = zipWith 
                                        {- map f' (zip l l')
                                        where f' (a,b) = a `f` b -}
-#endif
 
 -- ----------------------------------------------------------------------------
--- Data type constructors
+-- @node Data type constructors, H.o. fcts, Aux functions, ADT Matrix
+-- @section Data type constructors
 -- ----------------------------------------------------------------------------
 
-{- ----------------------------------------------------------------------------
-sqMatrix, i.e. arrays removed for GUM
-   ----------------------------------------------------------------------------
-
 sqMatrix :: (Integral a)  =>  Array (Int,Int) a -> SqMatrix a
-
 sqMatrix arr = SqMatrixC b [ [ (arr!(i,j)) | j <- [jLo..jHi] ]
                            | i <- [iLo..iHi] ] 
               where b@((iLo,jLo),(iHi,jHi)) = (bounds arr)
--}
-listSqMatrix :: (Integral a)  =>  MatBounds -> [a] -> SqMatrix a
 
+unSqMatrix :: (Integral a)  =>  SqMatrix a -> Array (Int,Int) a 
+unSqMatrix (SqMatrixC b@((iLo,jLo),(iHi,jHi)) m)  
+ = array b (concat [ [ ((i,j), (m!!(i-1))!!(j-1)) | j <- [jLo..jHi] ]
+                   | i <- [iLo..iHi] ])
+
+listSqMatrix :: (Integral a)  =>  MatBounds -> [a] -> SqMatrix a
 listSqMatrix b@((iLo,jLo),(iHi,jHi)) l = SqMatrixC b (take m (lol n l))
                                         where m = iHi - iLo +1
                                               n = jHi - jLo + 1
 
 lolSqMatrix :: (Integral a)  =>  MatBounds -> [[a]] -> SqMatrix a
-
 lolSqMatrix b l = SqMatrixC b l
 
 matBounds (SqMatrixC b _) = b
+matCont   (SqMatrixC _ m) = m
 
 vector :: (Integral a)  =>  [a] -> Vector a
-
 vector l = VectorC ((1),(n)) l
            where n = length l
 
 vecBounds (VectorC b _) = b
+vecCont   (VectorC _ v) = v
 
--- only for testing
-
-jaH (VectorC b l) = l
+unvector :: (Integral a)  =>  Vector a -> Array (Int) a
+unvector (VectorC b@(x,y) l) = array b (zip [x..y] l)
 
 -- ----------------------------------------------------------------------------
+-- @node H.o. fcts, Misc operations, Data type constructors, ADT Matrix
+-- @section H.o. fcts
+--
 -- Mapping and other general operations
 -- ----------------------------------------------------------------------------
 
+#if defined(STRATEGIES)
+matMapUnary :: (Integral a, NFData a)  =>  
+#else
 matMapUnary :: (Integral a)  =>  
+#endif
                (a -> a) -> SqMatrix a -> SqMatrix a
 
 matMapUnary f (SqMatrixC b mat) =
        SqMatrixC b (mat_map (mat_map f) mat)
 
-matCompwiseComp :: (Integral a, Integral b, Integral c)  =>  
+matCompwiseComp :: (Integral a, Integral b, Integral c
+#if defined(STRATEGIES)
+                     ,NFData a, NFData b, NFData c
+#endif
+                   )  =>  
                    (a -> b -> c) -> SqMatrix a -> SqMatrix b -> SqMatrix c
-
 matCompwiseComp f (SqMatrixC bnds@((iLo,jLo),(iHi,jHi)) mat) (SqMatrixC bnds' mat') = 
        if (bnds==bnds')
          then SqMatrixC bnds [ listCompwiseComp f (mat!!(k-1)) (mat'!!(k-1))
                             | k <- [iLo..iHi] ]
          else error "matCompwiseComp: Matrices have different bounds\n"
 
+#if defined(STRATEGIES)
+matFold :: (Integral a, NFData a)  =>  (a -> a -> a) -> a -> SqMatrix a -> a
+#else
 matFold :: (Integral a)  =>  (a -> a -> a) -> a -> SqMatrix a -> a
-
+#endif
 matFold f init (SqMatrixC _ mat) = foldl f init (mat_map (foldl f init) mat)
 
 
+#if defined(STRATEGIES)
+vecFold :: (Integral a, NFData a)  =>  (a -> a -> a) -> a -> Vector a -> a
+#else
 vecFold :: (Integral a)  =>  (a -> a -> a) -> a -> Vector a -> a
-
+#endif
 vecFold f init (VectorC _ mat) = foldl f init mat 
 
 -- ----------------------------------------------------------------------------
+-- @node Misc operations, Arithmetic Operations, H.o. fcts, ADT Matrix
+-- @section Misc operations
+--
 -- Misc operations
 -- ----------------------------------------------------------------------------
  
-size :: (Integral a)  =>  SqMatrix a -> Int
+-- Just for testing; demands computation of all elems of the matrix
+
+matCompact x = matFold max 0 (matMapUnary signum x)
 
+-- ---------------------------------------------------------------------------
+
+size :: (Integral a)  =>  SqMatrix a -> Int
 size (SqMatrixC ((iLo,jLo),(iHi,jHi)) mat) = 
                        if (iLo==jLo) && (iHi==jHi)
                          then iHi-iLo+1
@@ -191,7 +232,6 @@ replaceColumn j (SqMatrixC b@((iLo,jLo),(iHi,jHi)) mat) (VectorC _ v) =
 -- transp :: (Ix a, Ix b) => Array (a,b) c -> Array (b,a) c
 
 transp :: (Integral a) => SqMatrix a -> SqMatrix a
-
 transp (SqMatrixC b@((iLo,jLo),(iHi,jHi)) mat) = SqMatrixC b (transpose mat)
        {- 
        SqMatrixC b [ [ line !! (j-1) | line <- mat ] | j <- [jLo..jHi] ]
@@ -199,40 +239,61 @@ transp (SqMatrixC b@((iLo,jLo),(iHi,jHi)) mat) = SqMatrixC b (transpose mat)
 
 -- maxElem :: (Ix a, Ix b, Ord c) => Array (a,b) c -> c
 
+#if defined(STRATEGIES)
+maxElem :: (Integral a, NFData a) => SqMatrix a -> a
+#else
 maxElem :: (Integral a) => SqMatrix a -> a
-
+#endif
 maxElem (SqMatrixC _ mat) = maximum ( mat_map maximum mat )
 
+#if defined(STRATEGIES)
+maxElemVec :: (Integral a, NFData a) => Vector a -> a
+#else
 maxElemVec :: (Integral a) => Vector a -> a
-
+#endif
 maxElemVec (VectorC _ vec) = maximum vec
 
 -- ----------------------------------------------------------------------------
--- Arithmetic Operations
+-- @node Arithmetic Operations, I/O Operations, Misc operations, ADT Matrix
+-- @section Arithmetic Operations
 -- ----------------------------------------------------------------------------
  
 -- scalarMult :: (Ix a, Ix b, Num c) => c -> Array (a,b) c -> Array (a,b) c
 
+#if defined(STRATEGIES)
+scalarMult :: (Integral a, NFData a) => a -> SqMatrix a -> SqMatrix a
+#else
 scalarMult :: (Integral a) => a -> SqMatrix a -> SqMatrix a
-
+#endif
 scalarMult x = matMapUnary (x*)
                {-
                SqMatrixC b [ mat_map (x*) line | line <- mat ]
                -}
 
+#if defined(STRATEGIES)
+vecScalarQuot :: (Integral a, NFData a) => a -> Vector a -> Vector a
+#else
 vecScalarQuot :: (Integral a) => a -> Vector a -> Vector a
-
+#endif
 vecScalarQuot x (VectorC b vec) = 
               VectorC b (mat_map (`div` x) vec)
 
+#if defined(STRATEGIES)
+crossProd :: (Integral a, NFData a) => Vector a -> Vector a -> a
+#else
 crossProd :: (Integral a) => Vector a -> Vector a -> a
+#endif
+crossProd (VectorC _ vec) (VectorC _ vec') = sum (zipWith (+) vec vec')
+       -- foldl (+) 0 (listCompwiseComp (*) vec vec')
 
-crossProd (VectorC _ vec) (VectorC _ vec') =
-       foldl (+) 0 (listCompwiseComp (*) vec vec')
+-- @cindex determinant
 
 -- determinant :: (Ix a, Ix b, Num c) => Array (a,b) c -> c
 
-determinant :: (Integral a) => SqMatrix a -> a
+determinant :: (
+                 Integral a 
+                , NFData a
+               ) => SqMatrix a -> a
 
 determinant (SqMatrixC ((iLo,jLo),(iHi,jHi)) mat) 
        | jHi-jLo+1 == 1 =  let 
@@ -244,52 +305,45 @@ determinant (SqMatrixC ((iLo,jLo),(iHi,jHi)) mat)
                               [mat_2_1,mat_2_2] ] = mat
                            in
                              mat_1_1 * mat_2_2 -  mat_1_2 * mat_2_1
-       | otherwise      =  
-#ifndef SEQ
-#ifdef PAR_DET_SZ
-             sum (if (size<PAR_DET_SZ) then l_seq else l_par)
-#else
-            sum l_par
-#endif
-#else
-            sum l_seq
-#endif
-             where
-              newLine _ [] = []
-              newLine j line = [ line !! (k-1) | k <- [jLo..j-1] ] ++
-                               [ line !! (k-1) | k <- [j+1..jHi] ]
-
-               size = jHi-jLo+1
-
-#ifndef SEQ
-               l_par =  parmap 
-                    ( \ j ->
-                             (let 
-                        sign = if (even (j-jLo)) then 1
-                                                 else -1
-                        mat' = SqMatrixC ((iLo,jLo),(iHi-1,jHi-1))
-                               (map (newLine j) (tail mat))
-                        x = sign * ( (head mat) !! (j-1) ) * (determinant mat')
-                       in
-                        x) )
-                     [jLo..jHi] 
+       | otherwise      =  sum l_par 
+            where
+             l_par =   map determine1 [jLo..jHi]
+             determine1 j = 
+                        (if pivot > 0 then
+                          sign*pivot*det' 
+                        else
+                          0) -- `sparking` rnf sign 
+                        where
+                           sign = if (even (j-jLo)) then 1 else -1
+                           pivot = (head mat) !! (j-1)
+                           mat_h' = (map (newLine j) (tail mat))
+                           mat' = SqMatrixC ((iLo,jLo),(iHi-1,jHi-1))
+                                            mat_h'
+                           det' = determinant mat'
+
+#if 0
+                           strategyD r = 
+                             parList (parList rnf) mat_h'  `par`
+                             rnf det'         `par` 
+                             r0 r
 #endif
-
-               l_seq =  [ let 
-                        sign = if (even (j-jLo)) then 1
-                                                 else -1
-                        mat' = SqMatrixC ((iLo,jLo),(iHi-1,jHi-1))
-                               (map (newLine j) (tail mat))
-                        x = sign * ( (head mat) !! (j-1) ) * (determinant mat')
-                      in
-                        x
-                    |
-                      j <- [jLo..jHi] ]
+             tree_sum [] = 0
+             tree_sum [x] = x
+             tree_sum xs = (left+right)
+                           where (l,r) = splitAt (length xs `div` 2) xs
+                                 left = tree_sum l
+                                 right = tree_sum r
+             newLine _ [] = []
+             newLine j line = (pre ++ post)
+                               where                             
+                                pre  = [ line !! (k-1) | k <- [jLo..j-1] ]
+                               post = [ line !! (k-1) | k <- [j+1..jHi] ]
+
+{- seq determinant! -}
 
 -- matEqual :: (Ix a, Ix b, Eq c) => Array (a,b) c -> Array (a,b) c -> Bool
 
-matEqual :: (Integral a) => SqMatrix a -> SqMatrix a -> Bool
-
+matEqual :: (Integral a, NFData a) => SqMatrix a -> SqMatrix a -> Bool
 matEqual (SqMatrixC bnds@((iLo,jLo),(iHi,jHi)) mat) (SqMatrixC bnds' mat') = 
        if (bnds==bnds')
          then foldl (&&) True 
@@ -299,8 +353,7 @@ matEqual (SqMatrixC bnds@((iLo,jLo),(iHi,jHi)) mat) (SqMatrixC bnds' mat') =
          else error "matEqual: Matrices have different bounds\n"
 
 
-vecEqual :: (Integral a) => Vector a -> Vector a -> Bool
-
+vecEqual :: (Integral a, NFData a) => Vector a -> Vector a -> Bool
 vecEqual (VectorC bnds vec) (VectorC bnds' vec') = 
        if (bnds==bnds')
          then foldl (&&) True (listCompwiseComp (==) vec vec')
@@ -308,94 +361,57 @@ vecEqual (VectorC bnds vec) (VectorC bnds' vec') =
 
 -- matSum :: (Ix a, Ix b, Num c) -> Array (a,b) c -> Array (a,b) c -> Array (a,b) c
 
-matSum :: (Integral a) => SqMatrix a -> SqMatrix a -> SqMatrix a
-
+matSum :: (Integral a, NFData a) => SqMatrix a -> SqMatrix a -> SqMatrix a
 matSum = matCompwiseComp (+)
 
-
-matDif :: (Integral a) => SqMatrix a -> SqMatrix a -> SqMatrix a
-
+matDif :: (Integral a, NFData a) => SqMatrix a -> SqMatrix a -> SqMatrix a
 matDif = matCompwiseComp (-)
 
-{-
-matSum :: (Num a) => SqMatrix a -> SqMatrix a -> SqMatrix a
-
-matSum (SqMatrixC mat) (SqMatrixC mat') = 
-       if (bnds==bnds')
-         then SqMatrixC (array bnds
-                         [ i := mat!i + mat'!i | i <- indices mat ] )
-         else error "matSum: Matrices have different bounds\n"
-       where bnds  = bounds mat
-             bnds' = bounds mat'
-             
--- matDif :: (Ix a, Ix b, Num c) -> Array (a,b) c -> Array (a,b) c -> Array (a,b) c
-
-matDif :: (Num a) => SqMatrix a -> SqMatrix a -> SqMatrix a
-
-matDif (SqMatrixC mat) (SqMatrixC mat') = 
-       if (bnds==bnds')
-         then SqMatrixC (array bnds
-                         [ i := mat!i - mat'!i | i <- indices mat ] )
-         else error "matDif: Matrices have different bounds\n"
-       where bnds  = bounds mat
-             bnds' = bounds mat'
--}
-
--- matMult :: (Ix a, Ix b, Ix c, Num d) => 
---            Array (a,b) d -> Array (b,c) d -> Array (a,c) d
-
--- matMult :: (Integral a) => SqMatrix a -> SqMatrix a -> SqMatrix a
+-- @cindex mat mult
 
+{- parallel matrix multiplication -}
 matMult (SqMatrixC bnds mat) (SqMatrixC bnds' mat') = 
         SqMatrixC resultBounds 
-       [ [ let 
+#if defined(__PARALLEL_HASKELL__) || defined(__GRANSIM__)
+        (parMap rwhnf
+          (\i -> 
+           parMap rnf
+             (\j ->
+#else
+       (map (\i -> map (\j ->
+#endif
+             let
                line =    (VectorC ((jLo),(jHi)) (getLine i mat))
                column =  (VectorC ((iLo'),(iHi')) (getColumn j mat'))
-           in 
+             in
                crossProd line column
-          | j <- [jLo..jHi] ] 
-        | i <- [iLo..iHi] ]
+              )
+              [iLo..iHi]
+           )
+           [jLo..jHi] 
+        )
        where getLine i mat = mat !! (i-1)
              getColumn j mat = [ line !! (j-1) | line <- mat ]
+              size = iHi - iLo + 1           
               ((iLo,jLo),(iHi,jHi)) = bnds
               ((iLo',jLo'),(iHi',jHi')) = bnds'
               resultBounds 
                | (jLo,jHi)==(iLo',iHi')  = ((iLo,jLo'),(iHi,jHi'))
                | otherwise               = error "matMult: incompatible bounds"
 
-
-matAbs :: (Integral a) => SqMatrix a -> SqMatrix a
-
+matAbs :: (Integral a, NFData a) => SqMatrix a -> SqMatrix a
 matAbs = matMapUnary abs
 
 
-matSignum :: (Integral a) => SqMatrix a -> SqMatrix a
-
+matSignum :: (Integral a, NFData a) => SqMatrix a -> SqMatrix a
 matSignum = matMapUnary signum
 
 
-{-
-matAbs :: (Num a) => SqMatrix a -> SqMatrix a
-
-matAbs (SqMatrixC mat) = 
-       SqMatrixC (array (bounds mat)
-                  [ i := abs (mat!i) | i <- indices mat ] )
-
-matSignum :: (Num a) => SqMatrix a -> SqMatrix a
-
-matSignum (SqMatrixC mat) = 
-          SqMatrixC (array (bounds mat)
-                      [ i := signum (mat!i) | i <- indices mat ] )
--}
-
-matGcd :: (Integral a)  =>  SqMatrix a -> a
-
+matGcd :: (Integral a, NFData a)  =>  SqMatrix a -> a
 matGcd m = matFold gcd (maxElem m) m
 
 
-vecGcd :: (Integral a)  =>  Vector a -> a
-
+vecGcd :: (Integral a, NFData a)  =>  Vector a -> a
 vecGcd m = vecFold gcd (maxElemVec m) m
 
 
@@ -416,28 +432,29 @@ matBounds (SqMatrixC mat) = bounds mat
 -}
 
 matFromInteger :: Integer -> SqMatrix Integer
-
 matFromInteger n = SqMatrixC ((1,1),(1,1)) [[n]]
 
 -- ----------------------------------------------------------------------------
--- I/O Operations
+-- @node I/O Operations, Instances, Arithmetic Operations, ADT Matrix
+-- @section I/O Operations
 -- ----------------------------------------------------------------------------
 
 -- showsMatrix :: (Ix a, Ix b, Text c) => Array (a,b) c -> ShowS
 
 showsMatrix :: (Integral a) => SqMatrix a -> ShowS
-
 showsMatrix (SqMatrixC _ mat) = ( (++) ("Matrix: \n" ++
                                   (foldl (++) "" [ show line ++ "\n" 
                                                  | line <- mat ] ) ) )
 
 
 showsVector :: (Integral a) => Vector a -> ShowS
-
 showsVector (VectorC _ vec) = 
        ( (++) ("Vector: " ++ show vec) ) 
 
 -- ----------------------------------------------------------------------------
+-- @node Instances,  , I/O Operations, ADT Matrix
+-- @section Instances
+--
 -- Instance definitions for the ADT of Square Matrices and Vectors
 -- ----------------------------------------------------------------------------
 
@@ -446,14 +463,15 @@ instance (Eq a) => Eq [a] where
  l == l' = foldl (&&) True (listCompwiseComp (==) l l')
 -}
 
-instance (Integral a) => Eq (SqMatrix a) where
+instance (Integral a, NFData a) => Eq (SqMatrix a) where
  (==) = matEqual
 
+instance (Integral a) => Read (SqMatrix a) where
+ readsPrec p  = error "readsPrec of Matrix: Not yet implemented!\n"
 instance (Integral a) => Show (SqMatrix a) where
- --readsPrec p  = error "readsPrec of Matrix: Not yet implemented!\n"
  showsPrec p  = showsMatrix
 
-instance (Integral a) => Num (SqMatrix a) where                
+instance (Integral a, NFData a) => Num (SqMatrix a) where                
  (+) = matSum
  (-) = matDif
  (*) = matMult
@@ -464,20 +482,12 @@ instance (Integral a) => Num (SqMatrix a) where
                {- matFromInteger -}
 
 
-instance (Integral a) => Eq (Vector a) where
+instance (Integral a, NFData a) => Eq (Vector a) where
  (==) = vecEqual
 
-instance (Integral a) => Show (Vector a) where
- --readsPrec p  = error "readsPrec of Vector: Not yet implemented!\n"
- showsPrec p  = showsVector
-
+instance (Integral a, NFData a) => Read (Vector a) where
+ readsPrec p  = error "readsPrec of Vector: Not yet implemented!\n"
 
-{-
-instance  (Integral a) => Hom (Vector a) where
- hom = vecHom
- -- hom p (VectorC v) = vector (map (modHom p) (elems v)) 
+instance (Integral a, NFData a) => Show (Vector a) where
+ showsPrec p  = showsVector
 
-instance  (Integral a) => Hom (SqMatrix a) where
- hom = matHom
- --hom p = matMapUnary (modHom p)
--}
index 31bc83b..44e2101 100644 (file)
@@ -1,50 +1,70 @@
-{-
-$Id: ModArithm.hs,v 1.1 1996/01/08 20:07:35 partain Exp $
+-- Time-stamp: <Sat Jun 05 2010 01:37:17 Stardate: Stardate: [-28]3175.12 hwloidl>
+--
+-- Modular Arithmetic over Z_p
+-----------------------------------------------------------------------------
 
-This is revision: $Revision: 1.1 $
+-- @node Modular Arithmetic, ADT Matrix, Top, Top
+-- @chapter Modular Arithmetic
 
-Modular Arithmetic over Integrals with definition of Hom class.
-Changelog:
-$Log: ModArithm.hs,v $
-Revision 1.1  1996/01/08 20:07:35  partain
-Initial revision
-
---# Revision 1.2  1994/11/19  21:50:17  hwloidl
---# *** empty log message ***
---#
---# Revision 1.1  1994/11/19  02:00:05  hwloidl
---# Initial revision
---#
---# Revision 1.1  1994/11/19  02:00:05  hwloidl
---# Initial revision
---#
-----------------------------------------------------------------------  -}
-
-module ModArithm(modHom, modSum, modDif, modProd, modInv  {-, Hom(hom) -} )  where
+module ModArithm (modHom, modSum, modDif, modProd, modQuot, modInv  
+                 {-, Hom(hom) -} )  where
 
+{-# SPECIALISE 
+    modHom  :: Int -> Int -> Int
+  #-}
 modHom :: (Integral a) => a -> a -> a
 modHom m x = x `mod` m
 
 mapMod :: (Integral a) => (a -> a -> a) -> a -> a -> a -> a
 mapMod f m x y = modHom m (f x y)
 
+{-# SPECIALISE 
+    modSum  :: Int -> Int -> Int -> Int
+  #-}
 modSum :: (Integral a) => a -> a -> a -> a
 modSum = mapMod (+)  
 
+{-# SPECIALISE 
+    modDif  :: Int -> Int -> Int -> Int
+  #-}
 modDif :: (Integral a) => a -> a -> a -> a
 modDif = mapMod (-)  
 
+{-# SPECIALISE 
+    modProd :: Int -> Int -> Int -> Int
+  #-}
 modProd :: (Integral a) => a -> a -> a -> a
 modProd = mapMod (*)
 
-modInv _ 0 = 0
+{-# SPECIALISE 
+    modQuot :: Int -> Int -> Int -> Int
+  #-}
+modQuot :: (Integral a) => a -> a -> a -> a
+modQuot m x y = modProd m x (modInv m y)
+
+{-# SPECIALISE 
+       modInv :: Int -> Int -> Int
+  #-}
+modInv :: (Integral a) => a -> a -> a
 modInv m x = let 
-               (g,_,inv) = gcdCF m x
+               (g,foo,inv) = gcdCF m x
              in
                if (g /= 1) 
-                 then error ("modInv: Input values are not relative prime:\n " ++ (show m) ++ "\n " ++ (show x) ++"\n") 
+                 then modHom m inv -- error $ "modInv: Input values " ++ (show (m,x)) ++ " are not relative prime!" ++ ("** Wrong GCD res: gcd= " ++ (show g) ++ "; but x*a+y*b= " ++ (show (m*foo+x*inv)))
                  else modHom m inv
 
+gcdCF_with_check x y =
+  let 
+    res@(g,a,b) = gcdCF x y
+  in
+  if check_gcdCF x y res
+    then res
+    else error ("** Wrong GCD res: gcd= " ++ (show g) ++ "; but x*a+y*b= " ++ (show (x*a+y*b)))
+
+{-# SPECIALISE 
+    gcdCF :: Int -> Int -> (Int,Int,Int)
+  #-}
+gcdCF :: (Integral a) => a -> a -> (a,a,a)
 gcdCF x y = gcdCF' x y 1 0 0 1 
             where gcdCF' x 0 x1 x2 _ _   = (x,x1,x2)
                   gcdCF' x y x1 x2 y1 y2 | x<y       = gcdCF' y x y1 y2 x1 x2 
@@ -57,30 +77,8 @@ gcdCF x y = gcdCF' x y 1 0 0 1
                                                          gcdCF' y z y1 y2 z1 z2
 
 
-{- Main for testing gcdCF!
-
-main = let  
-         l1 = [7, 12, 62, 54, 55]
-         l2 = [3,  9, 30, 48, 15]
-         l = map ( \ (x,y) -> gcdCF x y ) (zip l1 l2)
-         showTuple (z,x,y) = "(gcd: " ++ show z ++ " Cofactors: " ++ show x ++
-                             " , " ++ show y ++ " )\n"
-         lshow = foldl (++) "" (map showTuple l)
-       in 
-         appendChan stdout ("First list: " ++ (showList l1  
-                            "\nSecond List: " ++ (showList l2 
-                            "\nResult: \n" ++ lshow))) 
-         abort done
--}
-
--- ---------------------------------------------------------------------------
-{-
-
-class Hom b where
- hom :: Integer -> b -> b
-
-instance Hom Integer where
- hom = modHom
- --hom m x = x `mod` m
+check_gcdCF :: (Integral a) => a -> a -> (a,a,a) -> Bool
+check_gcdCF x y (g,a,b) = if (x*a+y*b)==g
+                            then True
+                            else False
 
--}
\ No newline at end of file
diff --git a/parallel/linsolv/ParForce.hs b/parallel/linsolv/ParForce.hs
deleted file mode 100644 (file)
index b003fc8..0000000
+++ /dev/null
@@ -1,99 +0,0 @@
-{- 
-$Id: ParForce.hs,v 1.1 1996/01/08 20:07:35 partain Exp $
-
-This is revision: $Revision: 1.1 $
-
-  Module for forcing parallelism using par and friends.
-
-Changelog:
-$Log: ParForce.hs,v $
-Revision 1.1  1996/01/08 20:07:35  partain
-Initial revision
-
---# Revision 1.3  1994/11/23  01:07:23  hwloidl
---# Version for testing fwd mapping and hom sol.
---#
---# Revision 1.2  1994/11/19  21:50:18  hwloidl
---# *** empty log message ***
---#
---# Revision 1.1  1994/11/19  02:00:05  hwloidl
---# Initial revision
---#
---# Revision 1.1  1994/11/19  02:00:05  hwloidl
---# Initial revision
---#
-----------------------------------------------------------------------  -}
-
-module ParForce (parmap, parmap1, forcelist, forcelist1, par_iterate,
-                par_zip, par_zipWith) where
-
-#ifdef SEQ
-
-par :: a -> b -> b
-par x y = y
-
-seq :: a -> b -> b
-seq x y = y
-
-_parGlobal_ :: Int# -> a -> b -> b
-_parGlobal_ n x y = y
-
-#else
-
-import {-fool mkdependHS; ToDo: rm-}
-       Parallel
-
-#endif
-
-#ifdef GUM
--- GUM has only par 
-_parGlobal_ :: Int# -> a -> b -> b
-_parGlobal_ n x y = par x y
-
-_parLocal_ :: Int# -> a -> b -> b
-_parLocal_ n x y = par x y
-#endif
-
-forcelist [] = ()
-forcelist (x:xs) = seq x (forcelist xs)
-
-forcelist1 0 (x:xs) = ()
-forcelist1 n (x:xs) = seq x (forcelist1 (n-1) xs)
-
-parmap f []     = []
-parmap f (x:xs) = _parGlobal_ 1# fx 
-                               ( _parGlobal_ 6# (forcelist pmxs)  
-                                                (fx:pmxs) )
-                  where fx = f x
-                       pmxs = parmap f xs
-                 
-
-parmap1 f l = parmap' l
-             where parmap' []     = []
-                   parmap' (x:xs) = _parGlobal_ 1# 
-                                               pmxs 
-                                               (_parGlobal_ 2# 
-                                                            fx
-                                                            (fx:pmxs) )
-                                   where fx = f x
-                                         pmxs = parmap' xs
-
-
-par_iterate :: (a -> a) -> a -> [a]
-par_iterate f x = _parGlobal_ 13# fx 
-                       {- _parGlobal_ 14# (forcelist rest) -}
-                                       (fx : rest)
-                 where fx = f x
-                       rest = par_iterate f (f x)
-
-par_zipWith :: (a -> b -> c) -> [a] -> [b] -> [c]
-par_zipWith z (a:as) (b:bs) = _parGlobal_ 15# zab  
-                               (_parGlobal_ 16# (forcelist rest) 
-                                                (zab : rest)
-                               )
-                             where zab = z a b
-                                   rest = par_zipWith z as bs
-par_zipWith _ _ _ = []
-
-par_zip :: [a] -> [b] -> [(a,b)]
-par_zip = par_zipWith (\ a b -> (a,b))