Commit c3816570 authored by Christoph Ruegge's avatar Christoph Ruegge
Browse files

Move some stuff to lahelpers

parent 4d2de736
......@@ -3,11 +3,10 @@
{-# LANGUAGE BangPatterns #-}
import Data.Array.Repa
import Data.Array.Repa.Helpers (saveArrayAsPNG)
import Data.Maybe (fromJust)
import qualified Graphics.Rendering.Plot as P
import Numeric.LinearAlgebra.Helpers
import Numeric.LinearAlgebra.Repa (linearSolve)
import Numeric.LinearAlgebra.Repa.Conversion (repa2hm)
poisson :: Int -> (Double -> Double -> Double) -> Mat Double
poisson n f = computeS . reshape gridshape . fromJust . linearSolve laplace $ asColumn rhs
......@@ -50,8 +49,4 @@ poisson n f = computeS . reshape gridshape . fromJust . linearSolve laplace $ as
c = 1/h**2
main :: IO ()
main = do
let solution = poisson 50 (\x y -> x*y**4)
P.writeFigure P.PNG "01.png" (500, 500) $ do
P.setPlots 1 1
P.withPlot (1, 1) $ P.setDataset $ repa2hm solution
main = saveArrayAsPNG "01.png" (500, 500) $ poisson 50 (\x y -> x*y**4)
......@@ -2,13 +2,11 @@
{-# LANGUAGE BangPatterns #-}
import CG (runCG)
import Data.Array.Repa
import qualified Graphics.Rendering.Plot as P
import Numeric.LinearAlgebra.Repa.Conversion (repa2hm)
import Data.Array.Repa.Helpers (runCG, saveArrayAsPNG)
poisson :: Int -> (Double -> Double -> Double) -> IO (Array U DIM2 Double)
poisson n f = runCG 1e-6 laplace =<< computeP rhs
poisson n f = runCG laplace =<< computeP rhs
where
-- coordinates
h = 1 / fromIntegral n :: Double
......@@ -30,10 +28,5 @@ poisson n f = runCG 1e-6 laplace =<< computeP rhs
- get (Z:.jx:.(jy-1)) - get (Z:.jx:.(jy+1)) )
/ h**2
main :: IO ()
main = do
solution <- poisson 50 (\x y -> x*y**4)
P.writeFigure P.PNG "02.png" (500, 500) $ do
P.setPlots 1 1
P.withPlot (1, 1) $ P.setDataset $ repa2hm $ copyS solution
main = saveArrayAsPNG "02.png" (500, 500) =<< poisson 50 (\x y -> x*y**4)
......@@ -30,15 +30,13 @@
--
-- delay :: (Shape sh, Source r e) => Array r sh e -> Array D sh e
import CG (runCG)
import Data.Array.Repa as R
import Data.Array.Repa.Helpers (runCG, saveArrayAsPNG)
import Data.Array.Repa.Stencil
import Data.Array.Repa.Stencil.Dim2
import qualified Graphics.Rendering.Plot as P
import Numeric.LinearAlgebra.Repa.Conversion (repa2hm)
poisson :: Int -> (Double -> Double -> Double) -> IO (Array U DIM2 Double)
poisson n f = runCG 1e-6 laplace =<< computeP rhs
poisson n f = runCG laplace =<< computeP rhs
where
-- coordinates
h = 1 / fromIntegral n :: Double
......@@ -74,8 +72,4 @@ poisson n f = runCG 1e-6 laplace =<< computeP rhs
-- That's why we have to use R.map to do the division by h**2 here.
main :: IO ()
main = do
solution <- poisson 50 (\x y -> x*y**4)
P.writeFigure P.PNG "03.png" (500, 500) $ do
P.setPlots 1 1
P.withPlot (1, 1) $ P.setDataset $ repa2hm $ copyS solution
main = saveArrayAsPNG "03.png" (500, 500) =<< poisson 50 (\x y -> x*y**4)
{-# OPTIONS_GHC -Odph -rtsopts -threaded -fno-liberate-case -funfolding-use-threshold1000 -funfolding-keeness-factor1000 -fllvm -optlo-O3 #-}
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE ScopedTypeVariables #-}
-- CG is a simple, efficient iterative solver for linear systems
--
-- A x = y
--
-- with a positive definite matrix A, i.e. (x^T A x >= 0) for all x.
module CG where
import Control.Monad (foldM)
import Control.Monad.Identity (runIdentity)
import Data.Array.Repa as R hiding ((++))
-- CGState is the state of a CG iteration, i.e. everything that needs to be
-- passed to the next iteration.
data CGState sh = CGState { cgx :: !(Array U sh Double)
, cgp :: !(Array U sh Double)
, cgr :: !(Array U sh Double)
, cgr2 :: !Double
}
-- CG takes a function implementing a linear operator, a right hand side and an
-- initial guess and returns the (lazy, infinte) list of the CGStates of all
-- iterations. This way, the iterations can be decoupled from the stopping
-- criterion and the output of diagnostic information.
--
-- Computations use Repa's parallel computation in a straight-forward way and
-- needs to be performed in a monad. We use the Identity monad which does
-- essentially nothing except for providing the sequencing. All operations are
-- performed strictly using BangPatterns.
cg :: Shape sh => (Array U sh Double -> Array D sh Double) -> Array U sh Double -> Array U sh Double-> [CGState sh]
cg op rhs initial =
runIdentity $ do
!rInit <- computeP $ rhs -^ op initial
!r2Init <- normSquaredP rInit
return $ iterate cgStep (CGState initial rInit rInit r2Init)
where normSquaredP = sumAllP . R.map (**2)
scale a = R.map (* a)
cgStep (CGState !x !p !r !r2) =
-- We need a monad to avoid nested parallelism; the Identity
-- works fine.
runIdentity $ do
-- We need inline type annotation here because computeP is
-- polymorphic in its output type, and the type of q is not
-- fixed by any other expression below. Syntax like this
-- requires ScopedTypeVariables.
!(q :: Array U sh Double) <- computeP $ op p
!qp <- sumAllP $ q *^ p
let alpha = r2 / qp
!x' <- computeP $ x +^ scale alpha p
!r' <- computeP $ r -^ scale alpha q
!r2' <- normSquaredP r'
let beta = r2' / r2
!p' <- computeP $ r' +^ scale beta p
return $ CGState x' p' r' r2'
-- takeUntil is like takeWhile, but also returns the final iterate (no need to
-- waste it).
takeUntil :: (a -> Bool) -> [a] -> [a]
takeUntil _ [] = []
takeUntil predicate (x:xs)
| predicate x = [x]
| otherwise = x : takeUntil predicate xs
-- We want to be able to output some information about the iterates, i.e. run an
-- IO action over the list of iterates, before returning the final one. One
-- option would be something like
--
-- fmap last . forM iterates $ \it -> do [...]
--
-- However, this would perform the `fmap last` only *after* the IO actions are
-- done, therfore retaining the entire list of iterates in memory.
--
-- We want all but the final iterate to be garbage collected directly after
-- printing information. This can be done using a monadic fold, where we do not
-- actually accumulate a result but only return the last value.
process :: Monad m => [a] -> (a -> m ()) -> m a
process xs f = foldM (\_ x -> f x >> return x) undefined xs
-- runCG calls cg with initial guess 0, a stopping rule based on residuals
-- relative to the first initial guess, and outputs the residual while
-- iterating -- therefore, it returns a result in the IO monad.
runCG :: Shape sh => Double -> (Array U sh Double -> Array D sh Double) -> Array U sh Double -> IO (Array U sh Double)
runCG tol op rhs = do
let initial = computeS $ fromFunction (extent rhs) (const 0)
let steps' = cg op rhs initial
let r20 = cgr2 $ head steps'
let steps = takeUntil (\x -> sqrt (cgr2 x / r20) < tol) steps'
result <- process (zip [(1::Int)..] steps) $ \(n, cgs) ->
putStrLn $ show n ++ " " ++ show (sqrt $ cgr2 cgs / r20)
return $ cgx . snd $ result
......@@ -5,38 +5,24 @@ maintainer: c.ruegge@math.uni-goettingen.de
build-type: Simple
cabal-version: >=1.10
library
default-language: Haskell2010
hs-source-dirs: lib
exposed-modules: CG
build-depends: base >=4.8 && <4.9
, mtl
, repa
executable 01-hmatrix
default-language: Haskell2010
main-is: 01-hmatrix.hs
build-depends: base >=4.8 && <4.9
, lahelpers
, plot
, repa
, repa-linear-algebra
, poisson-equation
executable 02-repa
default-language: Haskell2010
main-is: 02-repa.hs
build-depends: base >=4.8 && <4.9
, plot
, lahelpers
, repa
, repa-linear-algebra
, poisson-equation
executable 03-stencils
default-language: Haskell2010
main-is: 03-stencils.hs
build-depends: base >=4.8 && <4.9
, plot
, lahelpers
, repa
, repa-linear-algebra
, poisson-equation
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment