packages feed

repa-examples (empty) → 1.0.0.0

raw patch · 16 files changed

+1270/−0 lines, 16 filesdep +basedep +dph-basedep +dph-prim-parsetup-changed

Dependencies added: base, dph-base, dph-prim-par, random, repa

Files

+ FFT/src/DFT.hs view
@@ -0,0 +1,61 @@+{-# LANGUAGE TypeOperators #-}++module DFT +	( dft+	, idft+	, dftWithRoots)+where+import Data.Array.Repa		as A+import Data.Ratio+import StrictComplex+import Roots++-- | Compute the (non-fast) Discrete Fourier Transform of a vector.+dft 	:: Shape sh+	=> Array (sh :. Int) Complex+	-> Array (sh :. Int) Complex++dft v+ = let	rofu	= calcRofu (extent v)+   in	force $ dftWithRoots rofu v+++-- | Compute the inverse (non-fast) Discrete Fourier Transform of a vector.+idft 	:: Shape sh+	=> Array (sh :. Int) Complex+	-> Array (sh :. Int) Complex++idft v+ = let	_ :. len	= extent v+	scale		= fromIntegral len :*: 0+	rofu		= calcInverseRofu (extent v)+   in	force $ A.map (/ scale) $ dftWithRoots rofu v+++-- | Generic function for computation of forward or inverse Discrete Fourier Transforms.+dftWithRoots+	:: Shape sh+	=> Array (sh :. Int) Complex		-- ^ Roots of unity for this vector length.+	-> Array (sh :. Int) Complex		-- ^ Input vector.+	-> Array (sh :. Int) Complex++dftWithRoots rofu arr+ = traverse arr id (\_ k -> dftK rofu arr k)+ ++-- | Compute one value of the DFT.+dftK	:: Shape sh+	=> Array (sh :. Int) Complex 		-- ^ Roots of unity for this vector length.+	-> Array (sh :. Int) Complex		-- ^ Input vector.+	-> (sh :. Int)				-- ^ Index of the value we want.+	-> Complex++dftK rofu arrX (_ :. k)+ = A.sumAll $ A.zipWith (*) arrX wroots+ where	sh@(_ :. len)	= extent arrX++	-- All the roots we need to multiply with.+	wroots		= fromFunction sh elemFn+	elemFn (sh :. n) +		= rofu !: (sh :. (k * n) `mod` len)+
+ FFT/src/FFT.hs view
@@ -0,0 +1,136 @@+{-# LANGUAGE TypeOperators, PatternGuards #-}++-- | Computation of Fast Fourier Transforms using the Cooley-Tuckey algorithm.+module FFT +	( fft+	, ifft+	, fft2d+	, fftWithRoots )+where+import Data.Array.Repa		as A+import Data.Ratio+import StrictComplex+import Roots++-- Vector Transform -------------------------------------------------------------------------------+-- | Compute the (fast) Discrete Fourier Transform of a vector.+fft	:: Shape sh+	=> Array (sh :. Int) Complex+	-> Array (sh :. Int) Complex++fft v+ = let	rofu	= calcRofu (extent v)+   in	force $ fftWithRoots rofu v+++-- | Compute the (fast) Inverse Discrete Fourier Transform of a vector.+ifft	:: Shape sh+	=> Array (sh :. Int) Complex+	-> Array (sh :. Int) Complex++ifft v+ = let	_ :. len	= extent v+	scale		= fromIntegral len :*: 0+	rofu		= calcInverseRofu (extent v)+   in	force $ A.map (/ scale) $ fftWithRoots rofu v+++-- Matrix Transform -------------------------------------------------------------------------------+-- | Compute the (fast) Discrete Fourier Transform of a square matrix.+fft2d 	:: Array DIM2 Complex+	-> Array DIM2 Complex++fft2d arr+ = let	rofu		= calcRofu (extent arr)+  	fftTrans	= transpose . fftWithRoots rofu+   in	fftTrans $ fftTrans arr+++-- Cube Transform ---------------------------------------------------------------------------------+-- | Compute the (fast) Discrete Fourier Transform of a 3d array.+fft3d 	:: Array DIM3 Complex+	-> Array DIM3 Complex++fft3d arr+ = let	rofu		= calcRofu (extent arr)++	transpose3 arr+	 = traverse arr +        	(\(Z :. k :. l :. m)   -> (Z :. l :. m :. k)) +            	(\f (Z :. l :. m :. k) -> f (Z :. k :. l :. m)) ++	fftTrans	= transpose3 . fftWithRoots rofu+	+  in	fftTrans $ fftTrans $ fftTrans arr++	+-- Worker -----------------------------------------------------------------------------------------+-- | Generic function for computation of forward or inverse Discrete Fourier Transforms.+--	The length of the roots vector must be the same as the values vector.+--	The length of these vectors must be a power of two.+fftWithRoots	+	:: Shape sh+	=> Array (sh :. Int) Complex		-- ^ Roots of unity.+	-> Array (sh :. Int) Complex		-- ^ Input values.+        -> Array (sh :. Int) Complex++fftWithRoots rofu v+	| not $ (denominator $ toRational (logBase 2 $ fromIntegral vLen)) == 1+	= error $ "fft: vector length of " ++ show vLen ++ " is not a power of 2"+	+	| rLen /= vLen+	= error $ "fft: length of vector is not the length of the roots"+	+	| otherwise+	= fftWithRoots' rofu v++	where	_ :. rLen	= extent rofu+		_ :. vLen	= extent v+++fftWithRoots'+	:: Shape sh+	=> Array (sh :. Int) Complex+	-> Array (sh :. Int) Complex+        -> Array (sh :. Int) Complex++{-# INLINE fftWithRoots' #-}+fftWithRoots' rofu v+ = case extent v of+	_ :. 2	-> fft_two   rofu v+	dim	-> fft_split rofu v dim++{-# INLINE fft_two #-}+fft_two rofu v+ = let	vFn' vFn (sh :. 0)  = vFn (sh :. 0) + vFn (sh :. 1)+	vFn' vFn (sh :. 1)  = vFn (sh :. 0) - vFn (sh :. 1)+   in	traverse v id vFn'+	+{-# INLINE fft_split #-}+fft_split rofu v vLen+ = let 	fft_lr = force $ fftWithRoots' (splitRofu rofu) (splitVector v)++	fft_l  = traverse2 fft_lr rofu + 		   (\(sh :. 2 :. n) _ -> sh :. n)+		   (\f r (sh :. i)    -> f (sh :. 0 :. i) + r (sh :. i) * f (sh :. 1 :. i))++	fft_r  = traverse2 fft_lr rofu +		   (\(sh :. 2 :. n) _ -> sh :. n)+		   (\f r (sh :. i)    -> f (sh :. 0 :. i) - r (sh :. i) * f (sh :. 1 :. i))++   in	fft_l +:+ fft_r++{-# INLINE splitRofu #-}+splitRofu rofu+ = traverse rofu+	(\(rSh :. rLen) 	-> rSh :. (2::Int) :. (rLen `div` 2))+	(\rFn (sh :. _ :. i) 	-> rFn (sh :. 2*i))++{-# INLINE splitVector #-}+splitVector v + = let	vFn' vFn (sh :. 0 :. i) = vFn (sh :. 2*i)+	vFn' vFn (sh :. 1 :. i) = vFn (sh :. 2*i+1)+   in	traverse v+		(\(vSh :. vLen)    -> vSh :. 2 :. (vLen `div` 2)) +		vFn'+
+ FFT/src/Main.hs view
@@ -0,0 +1,194 @@+{-# LANGUAGE ParallelListComp, PatternGuards, ScopedTypeVariables #-}++module Main where++import FFT+import DFT+import Roots+import StrictComplex+import Vector+import Matrix+import PPM+import ColorRamp++import Data.Array.Repa	as A+import Data.List	as L+import Data.Maybe	+import Prelude		as P+import Control.Monad+import System.Environment+++-- Arg Parsing ------------------------------------------------------------------------------------+data Arg+	-- | Use DFT instead of FFT+	= ArgDFT			++	-- | Compute the inverse transform.+	| ArgInverse+	+	-- | Use a vector step function for input.+	| ArgVectorRealStep	+		Int		-- length of 'on' part+		Int		-- total length of vector++	-- | Read the a real valued vector from this file.+	| ArgVectorReal+		FilePath++	-- | Read a PPM file as input.+	| ArgPPM +		FilePath++	-- | Write the result to this file.+	| ArgOutMagnitude+		FilePath++	-- | Write the magnitude of the transformed matrix as a PPM file.+	| ArgOutPPMMagnitude+		FilePath+		+	-- | Clip tranformed values to this level.+	| ArgOutPPMClip+		Double+		+	deriving (Show, Eq)+++parseArgs []		= []+parseArgs (flag:xx)+	| "-dft"		<- flag+	= ArgDFT : parseArgs xx++	| "-inverse"		<- flag+	= ArgInverse : parseArgs xx++	| "-vector-real-step"	<- flag+	, onlen:len:rest	<- xx+	= ArgVectorRealStep (read onlen) (read len) : parseArgs rest++	| "-vector-real"   	<- flag+	, fileName:rest		<- xx+	= ArgVectorReal fileName : parseArgs rest+		+	| "-ppm"		<- flag+	, fileName:rest		<- xx+	= ArgPPM fileName : parseArgs rest++	| "-ppm-clip"		<- flag+	, num:rest		<- xx+	= ArgOutPPMClip (read num) : parseArgs rest+	+	| "-out-magnitude"	<- flag+	, fileName:rest		<- xx+	= ArgOutMagnitude fileName : parseArgs rest+	+	| "-out-ppm-magnitude"	<- flag+	, fileName:rest		<- xx+	= ArgOutPPMMagnitude fileName : parseArgs rest+	+	| "-out-ppm-clip"	<- flag+	, level:rest		<- xx+	= ArgOutPPMClip (read level) : parseArgs rest+	+	| otherwise	+	= error $ "bad arg " ++ flag ++ "\n"+++help	= unlines+	[ "Usage: fft [args..]"+	, ""+	, "  -dft                                      Use (slow) DFT instead of (fast) FFT."+	, "  -inverse                                  Compute the inverse transform."+	, ""+	, "INPUT:"+	, "  -vector-real      <filename>              Read a real valued input vector from this file."+	, "  -vector-real-step <onlen::Int> <len::Int> Use a real valued step function for input."+	, "  -ppm              <filename>              Use a PPM file for input."+	, ""+	, "OUTPUT:"+	, "  -out-vector-magnitude <file-name>         Write the magnitute of the transformed vector to file."+	, "  -out-ppm-magnitude    <file-name>         Write transformed matrix to a ppm file."+	, "  -out-ppm-clip         <val::Int>           ... while clipping transformed values to this level."                  +	, ""+	, "NOTE:" +	, "  - For the fast algorithm, the length of the input vector/dimensions of the array"+	, "    must be powers of two."+	, ""+	, "  - When using PPM input files, pixels are converted to grey-scale and then used as"+	, "    the real values of the input matrix." ]+	+-- Main -------------------------------------------------------------------------------------------+main :: IO ()+main + = do	args	<- liftM parseArgs $ getArgs+	+	-- Decide which algorithm to use+	let alg	| elem ArgDFT args	+		, elem ArgInverse args	= idft++		| elem ArgDFT args	= dft+		+		| elem ArgInverse args	= ifft+			+		| otherwise		= fft++	main' args alg+++main' args alg++	-- | Transform a real-valued step function +	| [(onLength, vecLength)]	<- [(ol, vl) | ArgVectorRealStep ol vl <- args]+	= let	offLength	= vecLength - onLength+		step_real	= P.replicate onLength 1 ++ P.replicate offLength 0+		step		= P.map (:*: 0) step_real+		+		arr	= A.fromList (Z :. vecLength) step+		arrT	= alg arr+	  in	outVector args arrT+	+	-- | Transform some vector from a file+	| [fileName]	<- [f	| ArgVectorReal f <- args]+	= do	arr_real :: Array DIM1 Double	<- readVectorFromTextFile fileName +		let arr		= A.map (\r -> r :*: 0) arr_real+		let arrT	= alg arr+		outVector args arrT+	+	-- | Transform a PPM file.+	| [fileName]	<- [f	| ArgPPM f <- args]+	= do	let loadPixel r g b = sqrt (fromIntegral r^2 + fromIntegral g^2 + fromIntegral b^2)+		arr_double	<- readPPMAsMatrix loadPixel fileName+		let arr_real	= (A.map (\r -> r :*: 0) arr_double) :: Array DIM2 Complex+		let arrT 	= fft2d arr_real++		outPPM args arrT+	+	-- Not sure what you mean...+	| otherwise+	= putStr help+	+	+outVector args vec+	| [fileName]	<- [f | ArgOutMagnitude f <- args ]+	= writeVectorAsTextFile (A.map mag vec) fileName+			+	| otherwise+	= return ()++outPPM :: [Arg] -> Array DIM2 Complex -> IO ()+outPPM args arr+	| [fileName]	<- [f | ArgOutPPMMagnitude f <- args ]+	, mClipLevel	<- listToMaybe [l | ArgOutPPMClip l <- args]+	= do	let arr_mag	= A.map mag arr+		let arr_clipped	= maybe arr_mag+					(\level -> A.map (\x -> if x > level then level else x) arr_mag)+					mClipLevel+					+		writeMatrixAsNormalisedPPM fileName +		 	pixelGrey arr_clipped+		+		+pixelColor x = (x, x, x)+pixelGrey  x = (x, x, x)+
+ FFT/src/Roots.hs view
@@ -0,0 +1,43 @@+{-# LANGUAGE TypeOperators #-}++module Roots+	( calcRofu+	, calcInverseRofu)+where+import Data.Array.Repa+import StrictComplex+++-- Roots of Unity ---------------------------------------------------------------------------------++-- | Fill a vector with roots of unity (Rofu)+calcRofu +	:: Shape sh+	=> (sh :. Int) 			-- ^ Length of resulting vector.+	-> Array (sh :. Int) Complex++calcRofu sh@(_ :. n) + = force $ fromFunction sh f+ where+    f :: Shape sh => (sh :. Int) -> Complex+    f (_ :. i) =      (cos  (2 * pi * (fromIntegral i) / len))+		:*: (- sin  (2 * pi * (fromIntegral i) / len))++    len	= fromIntegral n+++-- | Fill a vector with roots of unity (Rofu)+--	for the inverse transform.+calcInverseRofu+	:: Shape sh+	=> (sh :. Int) 			-- ^ Length of resulting vector.+	-> Array (sh :. Int) Complex++calcInverseRofu sh@(_ :. n) + = force $ fromFunction sh f+ where+    f :: Shape sh => (sh :. Int) -> Complex+    f (_ :. i) =      (cos  (2 * pi * (fromIntegral i) / len))+		:*:   (sin  (2 * pi * (fromIntegral i) / len))++    len	= fromIntegral n
+ FFT/src/StrictComplex.hs view
@@ -0,0 +1,32 @@+{-# LANGUAGE TypeOperators, TypeSynonymInstances #-}++module StrictComplex+	( Complex(..)+	, mag+	, (:*:)(..))+where+import 	Data.Array.Parallel.Base ((:*:)(..))++-- | Strict complex doubles.+type Complex +	= Double :*: Double++instance Num Complex where+  (r :*: i) + (r' :*: i')	= r+r' :*: i+i'+  (r :*: i) - (r' :*: i')	= r-r' :*: i-i'+  (r :*: i) * (r' :*: i')	= r*r' - i*i' :*: r*i' + r'*i+  fromInteger n			= fromInteger n :*: 0.0++instance Fractional Complex where+  (a :*: b) / (c :*: d)		+ 	= let	den	= c^2 + d^2+		re	= (a * c + b * d) / den+		im	= (b * c - a * d) / den+	  in	re :*: im+	+	+-- | Take the magnitude of a complex number.+mag :: Complex -> Double+mag (r :*: i)	= sqrt (r * r + i * i)++
+ LICENSE view
@@ -0,0 +1,13 @@+Copyright (c) 2010 The DPH Team++ Permission is hereby granted, free of charge, to any person+ obtaining a copy of this software and associated documentation+ files (the "Software"), to deal in the Software without+ restriction, including without limitation the rights to use,+ copy, modify, merge, publish, distribute, sublicense, and/or sell+ copies of the Software, and to permit persons to whom the+ Software is furnished to do so, subject to the following+ condition:++ The above copyright notice and this permission notice shall be+ included in all copies or substantial portions of the Software.
+ Laplace/src/Main.hs view
@@ -0,0 +1,91 @@+{-# LANGUAGE BangPatterns #-}++-- | Solver for the Laplace equation+--	Writes a PPM file of the completed solution.+--+--	You can use the ImageMagick convert program to make a png+--	with	"convert out.ppm out.png"+--+import Solver+import Data.Array.Repa		as A+import ColorRamp+import PPM+import System.Environment++main :: IO ()+main + = do	args	<- getArgs+	case args of+	  [steps, fileInput, fileOutput]	+	    -> laplace (read steps) fileInput fileOutput++	  _ -> do+		putStr usage+		return ()+++-- | Command line usage information.+usage	:: String+usage	= unlines+	[ "Usage: laplace <iterations> <input.ppm> <output.ppm>"+	, ""+	, "  iterations  :: Int       Number of iterations to use in the solver."+	, "  input.ppm   :: FileName  ASCII 8 bit RGB PPM file for initial and boundary values."+	, "  output.ppm  :: FileName  PPM file to write output to."+	, "" +	, "  Format of input file:"+	, "      Boundary values are indicated in greyscale,"+	, "        ie from the list [(x, x, x) | x <- [0 .. 255]]"+	, "      Non-boundary values are indicated in blue,"+	, "        ie (0, 0, 255)"+	, "      Any other pixel value is an error." +	, ""+	]+			++-- | Solve it.+laplace :: Int			-- ^ Number of iterations to use.+	-> FilePath 		-- ^ Input file.+	-> FilePath		-- ^ Output file+	-> IO ()++laplace steps fileInput fileOutput+ = do+	-- Load up the file containing boundary conditions.+	(matBoundMask, matBoundValue)	+		<- readPPMAsMatrix2 loadPixel fileInput++	-- Use the boundary condition values as the initial matrix.+	let matInitial	= matBoundValue++	-- Run the solver.+	let matFinal	= solveLaplace+				steps+				matBoundMask+				matBoundValue+				matInitial++	matFinal `deepSeqArray` return ()++	-- Write out the matrix as a colorised PPM image	+	writeMatrixAsNormalisedPPM+		fileOutput+		(rampColorHotToCold 0.0 1.0)+		matFinal+++-- | Extract boundary mask and value from a pixel value.+loadPixel :: Int -> Int -> Int -> (Double, Double)+loadPixel r g b+	-- A non-boundary value.+ 	| r == 0 && g == 0 && b == 255	+	= (1, 0)++	-- A boundary value.+	| (r == g) && (r == b) +	= (0, (fromIntegral r) / 255)+	+	| otherwise+	= error $ "Unhandled pixel value in input " ++ show (r, g, b)+	+
+ Laplace/src/Solver.hs view
@@ -0,0 +1,82 @@+{-# LANGUAGE BangPatterns #-}+module Solver +	(solveLaplace)+where	+import Data.Array.Repa		as A+import qualified Data.Array.Repa.Shape	as S+import Data.Array.Parallel.Unlifted				(Elt)+import qualified Data.Array.Parallel.Unlifted			as U+++-- | Solver for the Laplace equation.+solveLaplace+	:: Int			-- ^ Number of iterations to use.+	-> Array DIM2 Double	-- ^ Boundary value mask.+	-> Array DIM2 Double	-- ^ Boundary values.+	-> Array DIM2 Double	-- ^ Initial state.+	-> Array DIM2 Double++{-# INLINE solveLaplace #-}+solveLaplace steps arrBoundMask@Manifest{} arrBoundValue@Manifest{} arrInit+ = go steps arrInit+ where	go i arr@Manifest{}+	   | i == 0	= arr+	   | otherwise	= go (i - 1) +			$ force+			$ applyBoundary arrBoundMask arrBoundValue+			$ relaxLaplace arr+++-- | Perform matrix relaxation for the Laplace equation,+--	using a stencil function.+--+--   Computation fn is+--	u'(i,j) = (u(i-1,j) + u(i+1,j) + u(i,j-1) + u(i,j+1)) / 4+--+relaxLaplace+	:: Array DIM2 Double+	-> Array DIM2 Double++{-# INLINE relaxLaplace #-}+relaxLaplace !arr+ = traverse arr id elemFn+ where+	_ :. height :. width	+		= extent arr++	{-# INLINE elemFn #-}+	elemFn get d@(sh :. i :. j)+	 = if isBorder i j+		 then  get d+		 else (get (sh :. (i-1) :. j)+		   +   get (sh :. i     :. (j-1))+		   +   get (sh :. (i+1) :. j)+	 	   +   get (sh :. i     :. (j+1))) / 4++	-- Check if this element is on the border of the matrix.+	-- If so we can't apply the stencil because we don't have all the neighbours.+	{-# INLINE isBorder #-}+	isBorder i j+	 	=  (i == 0) || (i >= width  - 1) +	 	|| (j == 0) || (j >= height - 1) +++-- | Apply the boundary conditions to this matrix.+--	The mask  matrix has 0 in places where boundary conditions hold+--	and 1 otherwise.+--+--	The value matrix has the boundary condition value in places where it holds,+--	and 0 otherwise.+-- +applyBoundary+	:: Array DIM2 Double	-- ^ Boundary condition mask.+	-> Array DIM2 Double	-- ^ Boundary condition values.+	-> Array DIM2 Double	-- ^ Initial matrix.+	-> Array DIM2 Double	-- ^ Matrix with boundary conditions applied.++{-# INLINE applyBoundary #-}+applyBoundary arrBoundMask arrBoundValue arr+	= A.zipWith (+) arrBoundValue+	$ A.zipWith (*) arrBoundMask  arr++
+ MMult/src/Main.hs view
@@ -0,0 +1,102 @@+{-# LANGUAGE PatternGuards #-}++import Solver+import Data.Array.Repa	+import Data.Maybe+import Matrix+import System.Environment+import Control.Monad++-- Arg Parsing ------------------------------------------------------------------------------------+data Arg+	= ArgSolver       String+	| ArgMatrixRandom Int Int+	| ArgMatrixFile   FilePath+	| ArgOutFile	  FilePath+	deriving Show++isArgMatrix arg+ = case arg of+	ArgMatrixRandom{}	-> True+	ArgMatrixFile{}		-> True+	_			-> False++parseArgs []		= []+parseArgs (flag:xx)+	| "-file"	<- flag+	, f:rest	<- xx+	= ArgMatrixFile f : parseArgs rest++	| "-out"	<- flag+	, f:rest	<- xx+	= ArgOutFile f	: parseArgs rest+	+	| "-random"	<- flag+	, x:y:rest	<- xx+	= ArgMatrixRandom (read x) (read y) : parseArgs rest+	+	| otherwise	+	= error $ "bad arg " ++ flag ++ "\n"++printHelp+	= putStr 	+	$ unlines+	[ "Usage: mmult [args..]"+	, ""+	, "  -random <height> <width>   Use a random matrix of this size."+	, "  -file   <filename>         Read a matrix from this file."+	, "  -out    <filename>         Write resulting matrix to this file."+	, ""+	, "  Format of matrix file:"+	, "    MATRIX"+	, "    <width> <height>"+	, "    <whitespace separated values..>"+	, "" ]+++-- | Get a matrix from a file, or generate a random one.+getMatrix :: Arg -> IO (Array DIM2 Double)+getMatrix arg+ = case arg of+	ArgMatrixFile   fileName	+	 -> readMatrixFromTextFile fileName++	ArgMatrixRandom height width	+	 -> genRandomMatrix (Z :. height :. width)	++			+-- Main -------------------------------------------------------------------------------------------+main :: IO ()+main + = do	args	<- liftM parseArgs $ getArgs+	main' args++main' args+	| [argMat1, argMat2]	<- filter isArgMatrix args+	, mArgOut		<- listToMaybe [s | ArgOutFile s <- args]+	= do	+		-- Get matrices from files, +		-- or generate random ones we were asked to.+		mat1		<- getMatrix argMat1+		mat2		<- getMatrix argMat2++        	mat1+          	 `deepSeqArray` mat2+          	 `deepSeqArray` return ()+		+		-- Run the solver.+		let matResult	= mmMult mat1 mat2++		matResult `deepSeqArray` return ()++		-- Write the output to file if requested.+		case mArgOut of +		 Nothing	-> return ()+		 Just fileOut	-> writeMatrixAsTextFile matResult fileOut+					+	| otherwise+	= printHelp++++
+ MMult/src/Solver.hs view
@@ -0,0 +1,14 @@+{-# LANGUAGE TypeOperators #-}++module Solver where+import Data.Array.Repa++mmMult	:: Array DIM2 Double+	-> Array DIM2 Double+	-> Array DIM2 Double++mmMult arr1@Manifest{} arr2@Manifest{}+	= fold (+) 0 +  	$ traverse2 arr1 (force $ transpose arr2)+      		(\(sh :. m1 :. n1) -> \(_ :. n2 :. m2) -> (sh :. m1 :. n2 :. n1))+		(\f1 -> \f2 -> \(sh :. i :. j :. k) -> f1 (sh :. i :. k) * f2 (sh :. j :. k))
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ lib/ColorRamp.hs view
@@ -0,0 +1,47 @@++module ColorRamp+	(rampColorHotToCold)+where+++-- Color Ramps  -----------------------------------------------------------------------------------+-- | Standard Hot -> Cold hypsometric color ramp.+--	Sequence is red, yellow, green, cyan, blue.+rampColorHotToCold +	:: (Ord a, Floating a) +	=> a +	-> a +	-> a +	-> (a, a, a)+	+{-# INLINE rampColorHotToCold #-}+rampColorHotToCold vmin vmax vNotNorm+ = let	+	v	| vNotNorm < vmin	= vmin+	 	| vNotNorm > vmax	= vmax+		| otherwise		= vNotNorm+	+	dv	= vmax - vmin	++	result	| v < vmin + 0.25 * dv+		= ( 0+		  , 4 * (v - vmin) / dv+		  , 1.0)+		+		| v < vmin + 0.5 * dv+		= ( 0+		  , 1.0+		  , 1 + 4 * (vmin + 0.25 * dv - v) / dv)+		+		| v < vmin + 0.75 * dv+		= ( 4 * (v - vmin - 0.5 * dv) / dv+		  , 1.0+		  , 0.0)+		+		| otherwise+		= ( 1.0+		  , 1 + 4 * (vmin + 0.75 * dv - v) / dv+		  , 0)+		+  in	result+
+ lib/Matrix.hs view
@@ -0,0 +1,86 @@++-- | Reading and writing matricies as ASCII files.+--	We use ASCII so we can generate and check simple test data by hand,+--	and we don't want to fool around with byte order issues.+--+--   Matrix file format is like:+--+--	MATRIX			-- header+--	100 100			-- width and height+--	1.23 1.56 1.23 ...	-- data, separated by whitespace+--	....+--+-- TODO: Merge this with PPM.hs+-- TODO: Merge this with Vector.hs. +--	 We should really have fns that read and write arrays of arbitrary dimension.+--+module Matrix+	( readMatrixFromTextFile+	, writeMatrixAsTextFile+	, genRandomMatrix)+where+import Data.List				as L+import Data.Array.Repa				as A+import qualified Data.Array.Parallel.Unlifted	as U+import Prelude					as P+import System.IO+import Control.Monad+import Data.Char+import System.Random+import Vector+++-- Reading ----------------------------------------------------------------------------------------+-- | Read a matrix from a text file.+readMatrixFromTextFile+	:: (U.Elt a, Num a, Read a)+	=> FilePath			-- ^ File name of matrix file.+	-> IO (Array DIM2 a)	++readMatrixFromTextFile fileName+ = do	handle		<- openFile fileName ReadMode+	+	"MATRIX"	<- hGetLine handle+	[width, height]	<- liftM (P.map read . words) $ hGetLine handle+	str		<- hGetContents handle+	let vals	= readValues str++	let dims	= Z :. width :. height+	let mat		= fromList dims vals++	return mat++++-- | Write a matrix as a text file.+writeMatrixAsTextFile +	:: (U.Elt a, Show a)+	=> Array DIM2 a			-- ^ Matrix to write.+	-> FilePath			-- ^ File name of output file.+	-> IO ()++writeMatrixAsTextFile arr fileName+ = do	file	<- openFile fileName WriteMode	++	hPutStrLn file "MATRIX"++	let Z :. width :. height	+		= extent arr++	hPutStrLn file $ show width ++ " " ++ show height+		+	hWriteValues file $ toList arr+	hClose file+++-- | Generate a random(ish) matrix.+genRandomMatrix +	:: DIM2 +	-> IO (Array DIM2 Double)++genRandomMatrix sh+ = do	uarr	<- genRandomUArray (A.size sh)+	return	$ fromUArray sh uarr+++
+ lib/PPM.hs view
@@ -0,0 +1,171 @@+{-# LANGUAGE ScopedTypeVariables, BangPatterns #-}++-- | Writing out matricies as PPM image files.+module PPM+	( writeMatrixAsNormalisedPPM+	, writeMatrixAsPPM+	, readPPMAsMatrix+	, readPPMAsMatrix2)+where+import qualified Data.Array.Parallel.Unlifted 	as U+import Data.List				as L+import Data.Array.Repa				as A+import Prelude					as P+import System.IO+import Control.Monad+import Data.Char+++-- Write ------------------------------------------------------------------------------------------+-- | Convert a matrix to a PPM image,+--	while normalising it to the maximum value present in the matrix.+writeMatrixAsNormalisedPPM +	:: FilePath				-- ^ Filename of output file.+	-> (Double -> (Double, Double, Double))	-- ^ Function for producing colors from data values,+	-> Array DIM2 Double			-- ^ Matrix of values (need not be normalised).+	-> IO ()++writeMatrixAsNormalisedPPM fileName colorFn arr+ = let	-- Use the maximum elem in the array as the white value.+	vals	= U.toList $ toUArray arr+	maxVal	= maximum vals++	-- Normalise the array to the range [0..1] for display.+	arrNorm	= A.map (/ maxVal) arr++  in	writeMatrixAsPPM fileName colorFn arrNorm+++-- | Convert a matrix to a PPM image.+--	Matrix elements should be normalised to [0..1]+writeMatrixAsPPM +	:: FilePath+	-> (Double -> (Double, Double, Double))	-- ^ Function for producing colors from data values.+	-> Array DIM2 Double 			-- ^ Matrix of values, normalised to [0..1]+	-> IO ()++writeMatrixAsPPM fileName colorFn arr+ = let		+	-- Break flat array data into individual rows+	Z :. width :. height	 +		= extent arr++   in do+	file	<- openFile fileName WriteMode+	hPutStrLn file $ "P3"+	hPutStrLn file $ show width ++ " " ++ show height+	hPutStrLn file $ "255"+		+	hWritePixels file colorFn $ toList arr+	hClose file+++-- | Write out pixel values to a file.+hWritePixels +	:: Handle +	-> (Double -> (Double, Double, Double))	-- ^ Function for producing colors from data values.+	-> [Double] 				-- ^ Data values.+	-> IO ()++hWritePixels h colorFn !xx+ = go xx+ where+	go []		= return ()+	go (x:xs)+	 = do	let (r, g, b)	= colorFn x+		hPutStr h $ showInt $ truncate (r * 255)+		hPutStr h $ " "+		hPutStr h $ showInt $ truncate (g * 255)+		hPutStr h $ " "+		hPutStr h $ showInt $ truncate (b * 255)+		hPutStr h $ "\n"+		go xs+	+showInt :: Int -> String+showInt i	= show i+++-- Read -------------------------------------------------------------------------------------------+readPPMAsMatrix +	:: (Int -> Int -> Int -> Double)	-- ^ Function for producing array values from RGB pixel values+	-> FilePath				-- ^ File name of ppm file.+	-> IO (Array DIM2 Double)		-- ^ Loaded matrix.+	+readPPMAsMatrix pointFn fileName+ = do	file	<- openFile fileName ReadMode+	+	"P3"		<- hGetLine file+	[width, height]	<- liftM (P.map read . words) $ hGetLine file+	_maxVal		<- hGetLine file+	vals		<- loadPixels pointFn file++	let dims	= Z :. width :. height+	let mat		= fromList dims vals++	return mat+++-- | Read the values of two separate matricies that are encoded in a single PPM file.+readPPMAsMatrix2 +	:: (Int -> Int -> Int -> (Double, Double))	+					-- ^ Function for producing array values from RGB pixel values+	-> FilePath			-- ^ File name of ppm file.+	-> IO   ( Array DIM2 Double	-- ^ Loaded matrix.+		, Array DIM2 Double)+		+readPPMAsMatrix2 pointFn fileName+ = do	file	<- openFile fileName ReadMode+	+	"P3"		<- hGetLine file+	[width, height]	<- liftM (P.map read . words) $ hGetLine file+	_maxVal		<- hGetLine file+	vals		<- loadPixels pointFn file++	let dims	= Z :. width :. height+	let mat1	= fromList dims $ P.map fst vals+	let mat2	= fromList dims $ P.map snd vals++	return (mat1, mat2)+	++-- | Load list of array data from the data part of a PPM file.+loadPixels +	:: (Int -> Int -> Int -> a)	-- ^ Function for producing array values from RGB pixel values.+	-> Handle 			-- ^ Handle of file to load from.+	-> IO [a]+	+loadPixels pointFn handle+ = do	str		<- hGetContents handle+	let ints	= readInts str+	let vals	= convertLine pointFn ints+	return vals+++-- | Convert groups of RGB components to values.+convertLine +	:: (Int -> Int -> Int -> a)	-- ^ Function for producing array values from RGB pixel values.+	-> [Int] 			-- ^ Ints read from the PPM file.+	-> [a]				-- ^ Output values.++convertLine fn vs+ = case vs of+	[]			-> []+	r : g : b : rest	-> fn r g b : convertLine fn rest+	_			-> error "PPM.convertLine: bad pixel"+	+-- | Read a string containing ints separated by whitespace.	+readInts :: String -> [Int]+readInts cs	= readInts' [] cs+ where	readInts' _ []	= []+	readInts' acc (c : rest)+		| isSpace c+		= if null acc +			then readInts' [] rest+			else read (reverse acc) : readInts' [] rest++		| isDigit c+		= readInts' (c : acc) rest++		| otherwise+		= error $ "unexpected char in PPM file " ++ show (ord c)+
+ lib/Vector.hs view
@@ -0,0 +1,137 @@++-- | Reading and writing vectors as ASCII files.+--	We use ASCII so we can generate and check simple test data by hand,+--	and we don't want to fool around with byte order issues.+--+--   Vector file format is like:+--+--	VECTOR			-- header+--	100			-- length of vector+--	1.23 1.56 1.23 ...	-- data, separated by whitespace+--	....+--+--+module Vector+	( readVectorFromTextFile+	, writeVectorAsTextFile+	, genRandomVector++	, readValues+	, hWriteValues+	, genRandomUArray)+where+import Data.List				as L+import Data.Array.Repa				as A+import qualified Data.Array.Parallel.Unlifted	as U+import Prelude					as P+import System.IO+import Control.Monad+import Data.Char+import System.Random+++-- | Read a vector from a text file.+readVectorFromTextFile+	:: (U.Elt a, Num a, Read a)+	=> FilePath			-- ^ File name of vector file.+	-> IO (Array DIM1 a)	++readVectorFromTextFile fileName+ = do	handle		<- openFile fileName ReadMode+	+	"VECTOR"	<- hGetLine handle+	[len]		<- liftM (P.map readInt . words) $ hGetLine handle+	str		<- hGetContents handle+	let vals	= readValues str++	let dims	= Z :. len+	let vec		= fromList dims vals++	return vec++readInt :: String -> Int+readInt str+	| and $ P.map isDigit str+	= read str+	+	+-- | Write a vector as a text file.+writeVectorAsTextFile +	:: (U.Elt a, Show a)+	=> Array DIM1 a			-- ^ Vector to write.+	-> FilePath			-- ^ File name of output file.+	-> IO ()++writeVectorAsTextFile arr fileName+ = do	file	<- openFile fileName WriteMode	++	hPutStrLn file "VECTOR"++	let Z :. len+		= extent arr++	hPutStrLn file 	$ show len+	hWriteValues file $ toList arr+	hClose file+	++-- | Generate a random(ish) vector.+genRandomVector +	:: DIM2 +	-> IO (Array DIM2 Double)++genRandomVector sh+ = do	uarr	<- genRandomUArray (A.size sh)+	return	$ fromUArray sh uarr+			+++-- Stuff shared with Matrix module -------------------------------------------------------------+-- | Write out values to a file.+hWriteValues+	:: Show a+	=> Handle +	-> [a] 				-- ^ Data values.+	-> IO ()++hWriteValues handle xx+ = go xx+ where+	go []		= return ()+	go (x:xs)+	 = do	hPutStr handle $ show x+		hPutStr handle $ "\n"+		go xs+++-- | Read a string containing ints separated by whitespace.	+readValues :: (Num a, Read a) => String -> [a]+readValues cs	= readValues' [] cs+ where	readValues' _ []	= []+	readValues' acc (c : rest)+		| isSpace c+		= if null acc +			then readValues' [] rest+			else read (reverse acc) : readValues' [] rest++		| isDigit c || c == '.' || c == 'e' || c == '-'+		= readValues' (c : acc) rest++		| otherwise+		= error $ "unexpected char in Matrix file " ++ show (ord c)+++-- | Generate a random(ish) UArray of doubles.+-- The std random function is too slow to generate really big vectors+-- with.  Instead, we generate a short random vector and repeat that.+genRandomUArray :: Int -> IO (U.Array Double)+genRandomUArray n + = do	let k		= 1000+    	rg		<- newStdGen+    	let randvec	= U.randomRs k (-100, 100) rg+	let vec		= U.map (\i -> randvec U.!: (i `mod` k)) (U.enumFromTo 0 (n-1))+	return vec++++
+ repa-examples.cabal view
@@ -0,0 +1,59 @@+Name:                repa-examples+Version:             1.0.0.0+License:             MIT+License-file:        LICENSE+Author:              The DPH Team+Maintainer:          Ben Lippmeier <benl@ouroborus.net>+Build-Type:          Simple+Cabal-Version:       >=1.6+Stability:           experimental+Category:            Data Structures+Homepage:            http://trac.haskell.org/repa+Bug-reports:         http://trac.haskell.org/repa/newticket+Description:+        Examples using the Repa array library.++Synopsis:+        Examples using the Repa array library.+++Executable repa-laplace+  Build-depends: +        base            >= 3   && < 5, +        random          >= 1.0 && < 1.1,+        dph-prim-par    >= 0.4 && < 0.5,+        repa            >= 1.0 && < 1.1++  Main-is: Laplace/src/Main.hs+  other-modules: Solver Matrix Vector ColorRamp PPM+  hs-source-dirs: Laplace/src . lib+  ghc-options: -Odph -threaded +++Executable repa-mmult+  Build-depends: +        base            >= 3   && < 5, +        random          >= 1.0 && < 1.1,+        dph-prim-par    >= 0.4 && < 0.5,+        repa            >= 1.0 && < 1.1++  Main-is: MMult/src/Main.hs+  other-modules: Solver Matrix Vector ColorRamp PPM+  hs-source-dirs: MMult/src . lib+  ghc-options: -Odph -threaded +++Executable repa-fft+  Build-depends: +        base            >= 3   && < 5, +        random          >= 1.0 && < 1.1,+        dph-prim-par    >= 0.4 && < 0.5,+        dph-base        >= 0.4 && < 0.5,+        repa            >= 1.0 && < 1.1++  Main-is: FFT/src/Main.hs+  other-modules: FFT DFT Roots StrictComplex Matrix Vector ColorRamp PPM+  hs-source-dirs: FFT/src . lib+  ghc-options: -Odph -threaded ++