### Merge Jacobi exercises from 20 and 22

parent fd87ac2f
 %% Cell type:markdown id: tags: # Exercises: linear algebra ## problem 1 Implement the following vector norms: - 1-norm - p-norm - infinitiy-norm For simplicity use Signatures like haskell norm1 :: Vec Double -> Double  %% Cell type:code id: tags:  haskell  %% Cell type:markdown id: tags: ## problem 2 - Generate a $15 \times 15$ Hilbert matrix $H = (h_{ij})_{i,j=1}^{15}$ with $h_{ij} = \frac{1}{i+j-1}$. - Compute the determinant of $H$. - Generate the vector $r = (1, 1, \ldots, 1)^T$ and the product $y = H r$ - Solve the equation using - $LU$ decomposition - singular value decomposition - <\> %% Cell type:code id: tags:  haskell  %% Cell type:markdown id: tags: ## problem 3 Which of the vectors $$v_1=\left(\begin{array}{c} 1\\2\\3\\4\\5 \end{array}\right),\quad v_2=\left(\begin{array}{c}-1\\27\\26\\1\\-27 \end{array}\right),\quad v_3=\left(\begin{array}{c} 160\\-48\\112\\-160\\48 \end{array}\right),\quad v_4=\left(\begin{array}{c} 120\\234\\-23\\-43\\29 \end{array}\right)$$ are orthogonal? %% Cell type:code id: tags:  haskell  %% Cell type:markdown id: tags: ## problem 4 The module Data.Array.Repa.FFTW provides functions for Fast Fourier Transforms (FFT) for Repa arrays of up to three dimensions. Use the 1-d version to implement a simple low-pass filter that works as follows: - First, an FFT is applied to the input signal. - The coefficients corresponding to frequencies over some pre-determined cutoff are set to 0, the other ones are left untouched. - Then an inverse FFT is performed. The map from index k of the resulting array to the corresponding frequency for an FFT of length n is haskell frq n k | 2*k < n = k | otherwise = k - n  Use the low-pass filter to denoise a signal created by adding uniformly distributed noise between $-0.03$ and $0.03$ to the function $x \mapsto exp(-x^2)$ for $x$ sampled on 200 points between $-2$ and $2$ and plot the results. Note: - The function, its noisy version and the filtered version should be real, but the FFT functions require and return Complex Double arrays. - The FFT functions work with arrays in representation F (foreign). To convert from and to D or U, copyS can be used. - A simple way to generate a list of n random numbers without using IO is haskell take n . randomRs (-0.03, 0.03) $mkStdGen 0  %% Cell type:code id: tags:  haskell  %% Cell type:markdown id: tags: ## problem 5 The *Jacobi* iterative algorithm for solving linear equations decomposes the matrix$A$into the diagonal part$D$and the remnant$R$: $$A = D + R,$$ and then iterates the following way: $$x^{( k + 1 )} = D^{ - 1} ( b - R x^{( k )} ),$$ where$x^{ ( k )}$is the$k$-th approximation of$x$. Implement a parallel Jacobi solver for linear systems using Repa arrays. Since Repa's mmultP (from Data.Array.Repa.Algorithms.Matrix) only works on DIM2 arrays, the right hand side should be a matrix as well, i.e. you should effectively solve $$A X = Y$$ for$X$, where$A$,$X$and$Y$are matrices. "Conventional" linear systems can then be solved by taking$X$and$Y$as matrices with only a single column. Iterations should stop either after a pre-determined maximum number of iterations, or when the difference between preceeding iterates is smaller in norm than some pre-defined tolerance. In order to actually use the parallelisation, this exercise should not be done in the notebook interface. You should also test your program on the 1d Poisson example from lecture 22 and examine its behaviour in ThreadScope. %% Cell type:code id: tags:  haskell  ... ...  %% Cell type:markdown id: tags: # Exercises: boundary value problems ## problem 1 implement a *Gauss-Seidel* algorithm for solving linear equations. The Gauss–Seidel method is an iterative technique for solving a square system of n linear equations with unknown x: $$A x = b$$ It is defined by the iteration $$L_∗ x^{( k + 1 )} = b − U x^{( k )}$$ where$x^{( k )}$is the$k$-th approximation or iteration of$x$,$x^{k + 1}$is the next value of$x$, and the matrix$A$is decomposed into a lower triangular component$L_∗$, and a strictly upper triangular component$U$:$A = L_∗ + U$We rewrite this a little bit to for the algorithm: $$x^{( k + 1 )} = T x^{( k )} + C$$ where $$T = - L_∗^{-1} U$$ and $$C = L_∗^{-1} b$$ Write a recursive function which makes this iteration. The iteration should stop both from a tolerance value of the difference of adjacent$x$and a maximum of iterations. For a test case use the matrix and right hand side from the lecture. %% Cell type:code id: tags:  haskell  %% Cell type:markdown id: tags: ## problem 2 Implement the *jacobi* iterative algorithm for solving linear equations. First decompose the matrix$A$into the diagonal part$D$and the remnant$R$: $$A = D + R$$ Then iterate the following $$x^{( k + 1 )} = D^{ - 1} ( b - R x^{( k )} )$$ where$x^{ ( k )}$is the$k$-th approximation or iteration of$x$and$x^{ ( k + 1 )}$is the next. Write a recursive function which makes this iteration. The iteration should stop both from a tolerance value of the difference of adjacent x and a maximum of iterations. For a testcase use the matrix and right hand side from the lecture. %% Cell type:code id: tags:  haskell  %% Cell type:markdown id: tags: ## problem 3 Write an implementation of the *conjugate gradient* (CG) method. This method solves $$A x = b$$ by an iteration of the form$x_0 = 0$,$r_0 = p_0 = b, \begin{align} \alpha_k &= \frac{{\lVert r_k \rVert}^2}{\langle p_k, A p_k \rangle} \\ x_{k+1} &= x_k + \alpha_k p_k \\ r_{k+1} &= r_k - \alpha_k A p_k \\ p_{k+1} &= r_{k+1} + \frac{{\lVert r_{k+1} \rVert}^2}{{\lVert r_k \rVert}^2} p_k \end{align}, stopping the iterations at the smallestK$for which$\lVert r_K \rVert < \text{tol}$. An interesting property of this method in contrast to the ones from the previous problems is that it does not need$A$in form of a matrix; it suffices to be able to apply$A$to any given vector$p$, so$A$can also be passed as a function Vec Double -> Vec Double. The CG implementation should therefore have a signature like this: haskell cg :: Double -> (Vec Double -> Vec Double) -> Vec Double -> Vec Double cg tol f rhs = [...]  Apply the method to the same test case as above, but this time reimplement it using the function haskell conv :: Vec t -> Vec t -> Vec t  from Numeric.LinearAlgebra.Repa, which implements a convolution without actually creating the matrix. Note that conv returns a vector that is larger than the input; use the boundary values to handle this. %% Cell type:code id: tags:  haskell  %% Cell type:markdown id: tags: ## problem 4 ## problem 3 Let A be a$20.000 \times 20.000$matrix, which entries are all zero besides the prime numbers$2, 3, 5, 7, \ldots, 224737$on the diagonal and the number 1 in all entries$a_{ij}$with $$|i − j| = 1, 2, 4, 8, \ldots , 16384.$$ What is the$(1, 1)$entry of$A^{−1}$? How can you check/visualize the structure of the matrix? Be sure to use a good method for saving the matrix. %% Cell type:code id: tags:  haskell  ... ...  %% Cell type:markdown id: tags: # Exercises: linear algebra ## problem 1 Implement the following vector norms: - 1-norm - p-norm - infinitiy-norm For simplicity use Signatures like haskell norm1 :: Vec Double -> Double  %% Cell type:code id: tags:  haskell import Data.Array.Repa hiding (map) import qualified Data.Array.Repa as R import Numeric.LinearAlgebra.Helpers (norm2, vec, Vec(..)) norm1 :: Vec Double -> Double norm1 v = foldAllS (\acc -> (acc +) . abs) 0 v normP :: Double -> Vec Double -> Double normP p v | p < 1 = error "normP: no such p-norm" | otherwise = (**(1/p))$ foldAllS (\acc -> (acc +) . (**p) . abs) 0 v normInf :: Vec Double -> Double normInf = maximum . map abs . toList -- test v = vec [1,2,-3] norm1 v norm2 v -- from the library normP 2 v normInf v  %% Cell type:markdown id: tags: ## problem 2 - Generate a $15 \times 15$ Hilbert matrix $H = (h_{ij})_{i,j=1}^{15}$ with $h_{ij} = \frac{1}{i+j-1}$. - Compute the determinant of $H$. - Generate the vector $r = (1, 1, \ldots, 1)^T$ and the product $y = H r$ - Solve the equation using - $LU$ decomposition - singular value decomposition - <\> %% Cell type:code id: tags:  haskell import Numeric.LinearAlgebra.Helpers import Numeric.LinearAlgebra.Repa import Data.Array.Repa hiding (map) import qualified Data.Array.Repa as Repa computeS $ones (ix1 2) :: Vec Double  %% Cell type:code id: tags:  haskell n = 15 hilb :: Mat Double hilb = computeS$ fromFunction (ix2 n n) $\(Z:. i :. j) -> 1 / fromIntegral (i+j+1) rhs :: Vec Double rhs = hilb app computeS (ones (ix1 n)) linearSolve hilb (asColumn rhs) linearSolveSVD hilb (asColumn rhs) hilb <\> rhs  %% Cell type:markdown id: tags: ## problem 3 Which of the vectors $$v_1=\left(\begin{array}{c} 1\\2\\3\\4\\5 \end{array}\right),\quad v_2=\left(\begin{array}{c}-1\\27\\26\\1\\-27 \end{array}\right),\quad v_3=\left(\begin{array}{c} 160\\-48\\112\\-160\\48 \end{array}\right),\quad v_4=\left(\begin{array}{c} 120\\234\\-23\\-43\\29 \end{array}\right)$$ are orthogonal? %% Cell type:code id: tags:  haskell let vs = (4><5) [ 1, 2, 3, 4, 5 , -1, 27, 26, 1, -27 , 160, -48, 112, -160, 48 , 120, 234, -23, -43, 29 ] in vs mul htrans vs  %% Cell type:markdown id: tags: ## problem 4 The module Data.Array.Repa.FFTW provides functions for Fast Fourier Transforms (FFT) for Repa arrays of up to three dimensions. Use the 1-d version to implement a simple low-pass filter that works as follows: - First, an FFT is applied to the input signal. - The coefficients corresponding to frequencies over some pre-determined cutoff are set to 0, the other ones are left untouched. - Then an inverse FFT is performed. The map from index k of the resulting array to the corresponding frequency for an FFT of length n is haskell frq n k | 2*k < n = k | otherwise = k - n  Use the low-pass filter to denoise a signal created by adding uniformly distributed noise between$-0.03$and$0.03$to the function$x \mapsto exp(-x^2)$for$x$sampled on 200 points between$-2$and$2$and plot the results. Note: - The function, its noisy version and the filtered version should be real, but the FFT functions require and return Complex Double arrays. - The FFT functions work with arrays in representation F (foreign). To convert from and to D or U, copyS can be used. - A simple way to generate a list of n random numbers without using IO is haskell take n . randomRs (-0.03, 0.03)$ mkStdGen 0  %% Cell type:code id: tags:  haskell {-# LANGUAGE FlexibleContexts #-} import Data.Complex import Data.Array.Repa hiding ((++), map) import qualified Data.Array.Repa as R import Data.Array.Repa.Eval import Data.Array.Repa.FFTW import Data.Array.Repa.Repr.ForeignPtr import System.Random  %% Cell type:code id: tags:  haskell frq n k | 2*k < n = k | otherwise = k - n  %% Cell type:code id: tags:  haskell cutoff :: Int -> Array F DIM1 (Complex Double) -> Array F DIM1 (Complex Double) cutoff cut x = ifft $computeS$ fromFunction (Z:.n) getElem where Z:.n = extent x fx = fft x getElem (Z:.k) | abs (frq n k) < cut = index fx (Z:.k) | otherwise = 0 cutoffReal :: Source r Double => Int -> Array r DIM1 Double -> Array D DIM1 Double cutoffReal cut = R.map realPart . cutoff cut . computeS . R.map (:+ 0)  %% Cell type:code id: tags:  haskell linspace :: Int -> Double -> Double -> Array D DIM1 Double linspace n a b = fromFunction (Z:.n) $\(Z:.i) -> a + fromIntegral i / (fromIntegral n - 1) * (b - a) n :: Int n = 200 xs :: Array U DIM1 Double xs = computeS$ linspace n (-2) 2 func :: Array U DIM1 Double func = computeS $R.map (\x -> exp (-x^2)) xs noise :: Array U DIM1 Double noise = fromList (Z:.n)$ take n $randomRs (-0.03, 0.03)$ mkStdGen 0 func' :: Array U DIM1 Double func' = computeS $func +^ noise  %% Cell type:code id: tags:  haskell import Graphics.Rendering.Chart.Easy import Graphics.Rendering.Chart.Backend.Cairo toRenderable$ plot (line "exact" [zip (toList xs) (toList func)]) toRenderable $plot (line "noisy" [zip (toList xs) (toList func')]) toRenderable$ plot (line "smoothed" [zip (toList xs) (toList $cutoffReal 30 func')])  %% Cell type:markdown id: tags: ## problem 5 The *Jacobi* iterative algorithm for solving linear equations decomposes the matrix$A$into the diagonal part$D$and the remnant$R$: $$A = D + R,$$ and then iterates the following way: $$x^{( k + 1 )} = D^{ - 1} ( b - R x^{( k )} ),$$ where$x^{ ( k )}$is the$k$-th approximation of$x$. Implement a parallel Jacobi solver for linear systems using Repa arrays. Since Repa's mmultP (from Data.Array.Repa.Algorithms.Matrix) only works on DIM2 arrays, the right hand side should be a matrix as well, i.e. you should effectively solve $$A X = Y$$ for$X$, where$A$,$X$and$Y$are matrices. "Conventional" linear systems can then be solved by taking$X$and$Y$as matrices with only a single column. Iterations should stop either after a pre-determined maximum number of iterations, or when the difference between preceeding iterates is smaller in norm than some pre-defined tolerance. In order to actually use the parallelisation, this exercise should not be done in the notebook interface. You should also test your program on the 1d Poisson example from lecture 22 and examine its behaviour in ThreadScope. %% Cell type:code id: tags:  haskell {-# LANGUAGE BangPatterns #-} import Data.Array.Repa as R import Data.Array.Repa.Algorithms.Matrix import Numeric.LinearAlgebra.Helpers jacobi :: Monad m => Array U DIM2 Double -> Array U DIM2 Double -> m (Array U DIM2 Double) jacobi :: Monad m => Array U DIM2 Double -> Array U DIM2 Double -> Double -> Int -> m (Array U DIM2 Double) -- simple stopping rule: always do 500 iterations jacobi a y = jacobi' 500$ computeUnboxedS $zeros (extent y) jacobi a y tol maxit = jacobi' maxit$ computeUnboxedS $zeros (extent y) where d = fromFunction (extent y)$ \(Z:.(!j):._) -> unsafeIndex a (Z:.j:.j) r = computeUnboxedS $R.traverse a id$ \get idx@(Z:.(!j):.(!k)) -> if j==k then 0 else get idx jacobi' 0 x = return x jacobi' n x = do rx <- r mmultP x jacobi' (n-1) =<< computeUnboxedP ((y -^ rx) /^ d) x' <- computeUnboxedP ((y -^ rx) /^ d) delta <- fmap sqrt . sumAllP . R.map (**2) $x -^ x' if delta < tol then return x' else jacobi' (n-1) x'  %% Cell type:code id: tags:  haskell a :: Array U DIM2 Double a = fromListUnboxed (Z:.4:.4) [ 2, -1, 0, 0 , -1, 2, -1, 0 , 0, -1, 2, -1 , 0, 0, -1, 2 ] x :: Array U DIM2 Double x = fromListUnboxed (Z:.4:.1) [1..4] y :: Array U DIM2 Double y = a mmultS x jacobi a y jacobi a y 1e-6 1000  ... ...  %% Cell type:markdown id: tags: # Exercises: boundary value problems ## problem 1 implement a *Gauss-Seidel* algorithm for solving linear equations. The Gauss–Seidel method is an iterative technique for solving a square system of n linear equations with unknown x: $$A x = b$$ It is defined by the iteration $$L_∗ x^{( k + 1 )} = b − U x^{( k )}$$ where$x^{( k )}$is the$k$-th approximation or iteration of$x$,$x^{k + 1}$is the next value of$x$, and the matrix$A$is decomposed into a lower triangular component$L_∗$, and a strictly upper triangular component$U$:$A = L_∗ + U$We rewrite this a little bit to for the algorithm: $$x^{( k + 1 )} = T x^{( k )} + C$$ where $$T = - L_∗^{-1} U$$ and $$C = L_∗^{-1} b$$ Write a recursive function which makes this iteration. The iteration should stop both from a tolerance value of the difference of adjacent$x$and a maximum of iterations. For a test case use the matrix and right hand side from the lecture. %% Cell type:code id: tags:  haskell {-# LANGUAGE FlexibleContexts #-} import Numeric.LinearAlgebra.Repa import Numeric.LinearAlgebra.Helpers import qualified Data.Array.Repa as Repa import Data.Array.Repa hiding (map, (++)) -- lets first create a matrix a and right hand side as in the lecture idxf (Z :. i :. j) | i == j = 2 | i == j-1 || i == j+1 = -1 | otherwise = 0 setupsystem :: Int -> (Mat Double, Vec Double, Vec Double) setupsystem n = (a,f,xi) where -- divide the interval x :: Vec Double x = computeS$ linspace n (0,1::Double) -- cut of boundary points (they are zero anyway) xi :: Vec Double xi = subVector 1 98 x -- build the matrix a :: Mat Double a = computeS $fromFunction (ix2 (n-2) (n-2)) (idxf) :: Mat Double f :: Vec Double f = computeS$ Repa.map (\xi -> (1/(fromIntegral n))^2*exp(xi)) xi  %% Cell type:code id: tags:  haskell gaussSeidel :: Mat Double -> Vec Double -> Int -> Double -> (Vec Double, Double) gaussSeidel a b tol = gsit (computeS $ones (ix1 size)) 1 where n = size(extent b) maxit = 20000 l = computeS$ fromFunction (ix2 (n-2) (n-2)) (\(Z :. i :. j) -> if (i-j) <= 0 then a ! (ix2 i j) else 0 ) u = computeS $fromFunction (ix2 (n-2) (n-2)) (\(Z :. i :. j) -> if (i-j) > 0 then a ! (ix2 i j) else 0 ) il = inv l t = computeS (Repa.map negate il) mul u c = il app b gsit x it | it > maxit || err < tol = (x,err) | otherwise = gsit xn (it+1) where xn :: Vec Double xn = computeS$ (t app x) +^ c err :: Double err = norm2 $computeS (xn-^x) tol = 1e-8 n = 100 (a,b,x) = setupsystem n (sol,_) = gaussSeidel a b n tol import Graphics.Rendering.Chart.Easy (toRenderable, layout_title, plot, line, (.=)) import Graphics.Rendering.Chart.Backend.Cairo () toRenderable$ do layout_title .= "2-point boundary value problem" plot (line "sol" [zip (toList x) (toList sol)])  %% Cell type:markdown id: tags: