Commit 80871ae2 authored by Jochen Schulz's avatar Jochen Schulz
Browse files

small fix (removed arr)

parent a8a5bd99
%% Cell type:markdown id: tags:
# poisson equation in 2d
- Poisson Problem also describes the stationary heat equation.
**Poisson equation** <br />
Find $u \in C^2(\Omega)\cap C(\overline{\Omega})$ with
$$\left \{ \begin{array}{rcll}- \triangle u & = & f & \mbox{in } \Omega\\
u & = & 0 & \mbox{auf } \partial \Omega\\ \end{array} \right.$$
for $\Omega=(0,1)^2$ and $f \in C(\Omega)$.
With $f = 0$ the equation is known as *laplace equation*
**Laplacian** <br />
$\triangle u := \sum_{i=1}^d \frac{\partial ^2 u}{\partial x_i^2}$
- equidistant gridsize: $h= \frac 1 N$,
$N \in \mathbb{N}$
- set of all gridpoints
$$Z_h := \left\{ (x,y) \in \overline{\Omega} \ \mid \ x=z_1h, \ y=z_2h \text{ mit } z_1,z_2 \in \mathbb{Z} \right\}.$$
- inner gridpoints: $\omega_h := Z_h \cap \Omega$
- approximation of $\frac{\partial ^2 u}{\partial x^2} (x,y)$
$$\frac{u(x -h,y) - 2 u(x,y) + u(x+h,y)}{h^2} = \frac{\partial ^2 u}{\partial
x^2} (x,y) + \mathcal{O}(h^2)$$
- approximation of $\frac{\partial ^2 u}{\partial y^2} (x,y)$
$$\frac{u(x ,y-h) - 2 u(x,y) + u(x,y+h)}{h^2} = \frac{\partial ^2 u}{\partial
y^2} (x,y) + \mathcal{O}(h^2)$$
- the sum gives the approximation for $\triangle u(x,y)$:
$$\frac{1}{h^2} \left( u(x,y-h) + u(x-h,y) - 4 u(x,y) + u(x,y+h) +
u(x+h,y) \right)$$
- definition $u_{i,j}:=u(ih,jh)$ ergibt an Gitterpunkten $(ih,jh)$
$$-u_{i,j-1} - u_{i-1,j} + 4 u_{i,j} - u_{i+1,j} - u_{i,j+1} = h^2 f_{ij}$$
with $i,j \in \{ 1, \dots , N-1 \}$ and $f_{ij}:=f(ih,jh)$.
- boundary conditions
$u_{0,i}=u_{N,i}=u_{i,0}=u_{i,N}=0$, $i=0, \dots ,N$.
- save 2D-values of $u$ in a vector: shape the inner variables
$$\begin{array}{cccc}
u(h,(N-1)h) & u(2h,(N-1)h) & \ldots & u((N-1)h,(N-1)h)\\
\vdots & \vdots & \vdots & \vdots \\
u(h,2h) & u(2h,2h) & \ldots & u((N-1)h,2h)\\
u(h,h), & u(2h,h) & \ldots & u((N-1)h,h)\\
\end{array}$$
this gives the vector $U_{i+(N-1)(j-1)}=u_{i,j}$.
Linear system for $U=(U_i)_{i=1}^{(N-1)^2}$
$$A U = F$$
with
- $F:=(f_i)_{i=1}^{(N-1)^2}$ mit $f_{i+(N-1)(j-1)}=f(ih,jh)$, $i,j \in \{1,
\dots ,N-1 \}$,
- $$A := \frac{1}{h^2} tridiag(-I_{N-1}, T, -I_{N-1}) \in \mathbb{R}^{(N-1)^2
\times (N-1)^2},$$
$$T := tridiag(-1,4,-1) \in \mathbb{R}^{(N-1)\times (N-1)}.$$
%% Cell type:code id: tags:
``` haskell
:extension NoMonomorphismRestriction
import Numeric.Container
import Numeric.LinearAlgebra.Algorithms
import Numeric.LinearAlgebra.Util
import Control.DeepSeq (force)
idx2grid :: Int -> Int -> (Int, Int)
idx2grid n i = (i `mod` n, i `div` n)
idxf :: Int -> (Int, Int) -> Double
idxf n (i, j)
| i == j = 4
| xi == xj && abs(yi - yj) == 1 = -1
| yi == yj && abs(xi - xj) == 1 = -1
| otherwise = 0
where (xi, yi) = idx2grid n i
(xj, yj) = idx2grid n j
poisson :: ([(Double, Double)], Matrix Double, Vector Double)
poisson = (mesh, reshape (n-1) $ head . toColumns $ loes, f)
where
n = 50
-- coordinates
x = linspace (n+1) (0,1::Double)
h = x @> 1
-- inner points
xInner = subVector 1 (n-1) x
mesh = [ (xj, yj) | xj <- toList xInner , yj <- toList xInner ]
-- laplace matrix
a = buildMatrix ((n-1)^2) ((n-1)^2) (idxf (n-1)) :: Matrix Double
-- right hand side
f = fromList $ map (\(x, y) -> (h^2)*x*(y^4)) mesh
-- solve linear system
loes = linearSolve a (asColumn f)
(me, loes, f) = poisson
```
%% Cell type:code id: tags:
``` haskell
import Graphics.Rendering.Plot
setPlots 1 1 >> withPlot (1, 1) (setDataset loes)
```
%%%% Output: display_data
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAASwAAAEsCAIAAAD2HxkiAAAABmJLR0QA/wD/AP+gvaeTAAAdOUlEQVR4nO2d2ZIbOa+ES1vbMxP/+7/rtJZzoTgOOpOdIIpVgnqc35WpBpdaIBdSIOrweDwWY0wdx+oFGPOnYyc0phg7oTHF2AmNKcZOaEwxdkJjirETGlOMndCYYuyExhRjJzSmGDuhMcXYCY0pxk5oTDF2QmOKsRMaU4yd0Jhi7ITGFGMnNKYYO6ExxdgJjSnGTmhMMXZCY4qxExpTzLl6AV9yv99fPMKfZj/DK+f6iuNx/X8h0JeHYgPdZWYxh12L/04OPr+27Ai72r/VYlbYb9t9ksPhsPqvbKCbC/lYaJDCj6PGFGMnNKYYO6ExxbyRMBPGGCsMNozZNp992+n2ODmTK9yVMIoLoz5tADHe4/FoDaD77Xb7+PhIrP53NnbCva9rtnvtjf5i+xf4bap7SDtCKKUAWWVF+xg41eFwYPm3NZiRYRg/jhpTjJ3QmGLshMYU81JhZj5m29U+7L73aDPjj4y26+Fk4e6psHDy5/XD4RBmwIgw8n6/n06ntjm+cmbWCSe9Ytvu88rE5I2+a/d3s59XswXzKS9aiTkej63nHA6/pY6Fsx+Pxw3FYT+OGlOMndCYYuyExhRTKcxMZv1vrkzsGjI9Ho/NByycrvvJ+F+zzCfEQMrL4/Foo8TuaGDfHhEcXbEwA7yzksH35fzyUvZ8qWaWNzL73qdLdB+xH2eF+AmftGLmswkZMOCT5zO6hriak07ox1FjirETGlOMndCYYoozZnbtvqvwsyIkY+NsjJdaAAQq2XhYd2f2Cwh/jdZGcff7HYQTobuw/f1+hwwYLeR017MVr3PCFSrIi2/K1Ph8ODz+pH12Pdseb+iEWdVnEi29cFOnoR2Px1aqgYyZ7oxAa/9e6igweeXmL/xb2c+rkTP22fVs3n0GdpKu26xeHqetvfLoHBMaU4yd0Jhi7ITryT4O7Tpd1j583Mqy+fG+bPDyBWwZE86HHPpQtZLBA+6tfExOt8JeGGTH7xpoe5hu+Z1sAMmk6vlmhRlIkTmdTq3B8Xi8XC4wo85ra4/oer2KlYcUF3p6ZyVj3kvhLp/0Uu2E8+Mv8vCzZy9sAqEHhsXRwAn1CDwaby9sPwnPxgx+HDWmGDuhMcXYCY0p5o32E44IA6koZTIkG4mpUkGUXg+MlrXfOyLNhpRhE1gREwol5nA4wN6lx+9paLpS03M0IeTAyf+DMmYm75vN7bNucLvd9HRZt9nWS/Xh7+2EwLaFmxZyQt4u2F6dp04DTtuOD4WevJ/QmO+NndCYYuyExhQzFRPyLqyWbQPC7ie6+2RIySHTTAy2SQgnBgynW2GQDSlFs/uJJvUeJf6xnqfTVWfC/YdwNlrj7KEBW1bg7t73KftdBb3NvQgMQnkThBntVCPLAyfc9nBG7EX3rsEyzAphBtRR/W6zhW5F8OHb7QZO2F4+0IEszBjzvbETGlOMndCYYt6o0FMYM2S7p0KUeSEnK2zokHIkoNXdUyHlyOHomDCr04ToIPAhCzFxDMlCi9ZdoBLUtkoM8FInXJJCTtYrJlNSsspH6EVamLndbin7SWFmhZfOOO2yszDD2ibsD0ylrcFuw8PhcL1e2wHP5zOM3x7+5H5CP44aU4yd0Jhi7ITGFPPqmFAwGexmdZ1JHWjEYGa0rM7Exikda0V8nh1Q24dvMgObroFe3rj9trpLyI5OmL1s/NdQf0vpNPPyYzZDBez1esLxtZCzQqcJ7fXRpa4Ow38Nq8hAU7/bjM82CC3tgmG7IAsz0ATVZ5JXv58wZa/9Nvzynhx/cj3hd8SIwbg9exEYhOpodvzsdyLTGuga2JyGpitwg3rJ689enSW692ZwTGhMMXZCY4qxExpTzEuFmTBS3zBG4k/mM04mhZNUxkxoHy5vUgdasdOqberL3f2TFmb0BkLOmOEqMpz41v67tWdd53a7tWIMHN17VeAG+FxnpQs9YMpLefZUtuSkffYuX/GlEHodTJe1T8mnS1LWSnngQolmfC3YvjU4nU7t+rt1olJCzgx+HDWmGDuhMcXYCY0pZrbQk9hzxaQeo7shxIxyo2OepRfk6KggJXXMx3gszLT22e5dA728VIg4EuELwqIy9983BPI76GFAUG6u1yvYaGFGnw049ixbFnrK2k86VWifVT7YaTf0oqxThfa3223GaUeEmVANbpvZq6PRrx9cKHGMhRb2YZBPwYArcLd/1d/Ik07ox1FjirETGlOMndCYYt6o0FNonIoYQ4PuX8dD1pEQMYyRtg1BJ0NEDvmywsx+MeGSLPT0eDwgaOSjg1edwY/1IPyEN48+HM2OW5n4MmgDtt82o2VboSVr/3g8ILmJ73It5FyvVz1dOP7mXto2Vzhh+0mYMQNOyMJM+wkYiw9/rYTVUZDB2uNlafqrkUfw46gxxdgJjSnGTmhMMW9U6AmYDHZXTBfqOpP24QLC8VMRNYfEM0JLGPJlM2a4qes46Yra2bMxYj/enKRSHQ1Pk+4+f1dldZpUBkyY0QJCC9h3pwN7vTzWaWD8bXWarDCj7+ORFBmRtgbvlF9I2uHibpfLpf0ry05w9uDkiGMJqSyDn23W2vNdlRp/snt3hNThhN2z60/9VzMynYDLOq042L0PZzWOCY0pxk5oTDF2QtNnW+3hzak92C1jwmxIk+rOI4RP7fNKTKhkaPtUBkwozIDQcr1es0JOankpYaZ7MtvmioAqVXKbhRkeDYQZ0Gba0/uswN3+VafIvJcws61UkIqVs7H4i+35NtX3cddeDDj5nRKOn03iG1EyhCtyBW62ge78G4Zo6rMXFnrSi8/ix1FjirETGlOMndCYYl73Y334DJ0VckZmzAo50EzFeCu2Ao0rJSy0hMLMCiEndTipOlHLXMYMb2WCtDWYju0/Pz9BuYHxocYMH913zZhZ4Xi6+WKhZVslg518Mo9MexHcKNm0uOzy+DtLC6TdT34ROuFDbuo9Ho/wHXQ+n9vpTqdTu36uwM1OLnSaSSf046gxxdgJjSnGTmhMMbMVuFMv785GgLp7OHgqwgwNQmGGg3URYoUGk/YjMV52q1QqJtxWmFl65Xp1oafURjM4umdMCNXl4eqPH0tIZQXuFc2U/KjdJrQPvSK0D4WT1h5KCY3cRlB6iOVQfdu19l0n3FaYYfWitc8Wejqfz0I4gd2GT/t2hPP5LPLUnjqNruwEzWUCP44aU4yd0Jhi7ITGFLPjVqbQYCZwnyer03D3MFgXf+0aiCZ/wjEYN7VQlNKZwoh675gQ+sKP6RA0sizEP76Ls/f8635KDFBZgXtFU992k8JMKD9mU0xCoSUlzIDQcr1etTDz+fkJ9np5KwpPtdNNOiGQTVs7n89amOEUGRgQhBm+WEKngTOfxY+jxhRjJzSmGDuhMcX8QU5YK/xkDSaXt6L7pE6W7Z5a4aS4sGK6vW+Plm9cBn9EnByXH0NdJ0zQCZMqdAZMKLSkhBkWWtiem2CvdZpsAk3Wq1lZaT85nU6tshK+s54NhDDDhZ748PXiU7yvEy7JXxG6f535vhy5b1L281rx+HfKiD14UWgwkyS45NXR9q/PZvsJ7ADsbgiE2eHdFXz4bXdenr73vjqQEf6gx1Fj3hM7oTHF2AmNKea9Xo02OZr+JMzb0kpMNyNE2+sgasXWIbAXQktXmIEEHRZy2uVBwo1OoOHD2VaY6cZ4kKfWvnT+drux7gJnA+o+ccZMu4BvXOiJSQW7k0LLtvZdnWNcKekqGVnlQ6uv3D37HaHtU9sd54UZ0Vx6yk1Y/U3sP4Qmf6Iv7iR+HDWmGDuhMcXYCY0p5n0rcGcHzI4fxpwrYrxxHSg0CKOUFTGe2IzTFX5mYsKujrX8jo4JQZhhY13oqdVpFgoRQdaCg32uVsty30OYWRG8vlKYWeGEOsVEy6FZYWaFE4bV07SaOu+EqcNfkk7IrxwUTng6ncKyV+CEXJBbH86GTujHUWOKsRMaU4yd0JhiKndRbPhz5/x02YAzO2AYgoYpJqE2oIUQLfx07UPhp7XXP9Z3Q+Lla7gCN6B/fF8oTmP78GRqgUCsLcvG1da2VWJS9iNCi4itR+xTwszIXSuUjzvlhcGALC2EFbVT9pzmxmlxcDgpJ1+SGTPsw7A1qR3tdDp9fn629py2BnlqM2lrkz7px1FjirETGlOMndCYYt63vMXesk0q/lxnsOFok/H2ivFDe+4+rmQ8mzMnMIznOWhcPdreVO4nnOyePY+Td3k4vk5D66o+QsnoJqBooSWV0QLKxP1+10pG9k1sLMzw6Vokqf2EMP7tdoM3Z8L2yPP53B4vVHZ6boxqx9fCzGQF7uL9hFnjlB++m/2u03X/59lv/LB71h7gZNFtp9P22UJPkzgmNKYYO6ExxdgJjSnmjdLWsg/Z8wGPDhJCoYVzOMSA995+QhHrcx4Wv0pNJ+joDBgoyM32LORMCjM6rQzgtDU+27y5qTWG/YT8orj2cLpnQ2fMiMVnKf6JYm8/TCkNIwPCv7U0MrO80H6kKb41QvtwQx17RSprr3v+f8H7CVPThfbZs6Ev1qRP+nHUmGLshMYUYyc0ppg3KvQU2mdHyHbnMCC1nlTMxp9oHSg7/jNG0uOzkNNON18gXEy3DAgz0IQMGFZuIEkNDke/6ozzkyBjRl+syTvzvQo98QjhJ6kRdPeUUvJICjOhE2qpILTvCiF6fG3P4q1Oi2N77eT69hipwN1++NwB2Brw8kSWH/y1W+gpdS+l8OOoMcXYCY0pxk5oTDHvu59wntRT+4r4c2b8kYA5JeR0Q9Bx4YeFk6xuFMaEqR/rF9rKBGEhbzVqm+HrssN0pZQSM/lj/X/KCbNekXIDvstDA6GL8Cehvb6tQ2HmLgtJsf1kBe4VwkzbZGGGLxakyIA6ynIoNGcqcPN30Ax+HDWmGDuhMcXYCY0p5j8VE86w4SP+JgsIhSK94EnVJ5xupLuIkEdmDBfcNrVMFa5ncm2TfGMn1OeRP1xxW+jpwuuavS3CFA0tFWihJcyA0UIOV+zWiWChkJMVZliBFHXy77SfUFfghlepPYUZ0LHE1ZlUR/04akwxdkJjirETGlOMndD0yUoRWSVjUuqYXN4k2472jYWZvdn1LlmhRoZ6INtrYUb8lWfUus4S7T9kYYaTTrQww2WdRKGnrm4EhZvC+uXjFbjh0LLYCRUb/m/QlV5Tfpgdf0SenRx/w+V1BxyfboUWHTazV2c1fhw1phg7oTHF2AmNKcYx4ZuyeQTFBlqY0Ylg3QI2wqA7vmguMiFmIaEFlBvOsEm9g/65Wn28Ym1Z7IRvxCuFmRH71HZHtk9t6oUmbCAMN+nC+CsqcLfN4/EYHj646DKBH0eNKcZOaEwxdkJjivnGMSGXFenafNVc0X1kSSn79yH723SoG2lhY+m96gyCLr46UKmpVWK6ykrKXq92v1/ql2/thEDoANqgW+N51/Vsjs75CIWZRd5qI0JLSpgJb/QWPpnPLX/taJwWB01R5+pOFbhHjndk5SP4cdSYYuyExhRjJzSmmB1jwhHlIxxhtf3m03UPB5Se8e4rFrAroTAzMsJ4c8SAlyeEllAomtSZ9OIneZ0w82Kl8XnTb3uyXjz+TiP/IivMaF1n84yZ1pgzYODkh69Gg7S1UJhxoSdj/iDshMYUYyc0ppjiH+v1a8ez9lo4GZluV3tu7r08TTj+N2JbGWnFgDNUOiFvV5mxH3nLuW5mxw/thXb6guUtGbI607waqfcT8ncWlHXiNDS910krKyy0dD/5arWT/HfS1phd1cs/kMlfHVL2/I2z4keOcbGXZ19xOKtxTGhMMXZCY4qxExpTzJYxYZinNp80M7MhkJUSoV4uVDuoayDsQRhgLRSUBl4ezw61jMLxIc7Jnj1AT7ctHIBpYYaFFtZphJAT1q3S52dSp3lfYabrgeJcTN4ToVw5P+CMfVa9nF9POSlhZiSrLpWUlx1/Bj+OGlOMndCYYuyExhTzRmlroVKy9LSQ1j7UUcCA7dsuusbzU3dpDbRyAAZdoSVcXvvJ8XgU9rw83X2JlKFQB9K61IuzJsLpdMSYDfkmj+6N0tb4r9o+1CpfbB+uP5xOuHRonz2c7pLGDz/bfXNW5M2NZ8A8/p+vuluYMeY/hZ3QmGLshMYU8+pCT2GkBPbCgO1TQk6oA41IEdpeCDNsoLuHQg5nhHDGzKQwE+pe42d7ie6E7tXckP22RKzgpYWeZpSM0D7lgcuYj006oe4ONaSPx6PuzvaiOzTZvjt+e3W4qceHAfV0S+96iRSZTRwyK8yI/Yfb4sdRY4qxExpTjJ3QmGLeKGOmG+uLZra7tt9DmEnFnDAg/JUTekbWs+Hytj0cXl6oxITTzZCN97aND98oY2bkombt2yYrDVo5APuF8s54wNb4QVvgdPN0Ogkxs9u9/eR0Op1OJ708MGgrTB8OB1hAa7zQDr2nfWsAFa9hNF4/7MGDywHo79OugUZ7Ef9V2/9n9xMasytCHT2sqiu1GseExhRjJzSmmI0fR7XyschH+VApyXbf1p67LxntoSskzEyXOrpwednxs91ZZ9ILzspsKwxSzV3ZN20NmqyspB6s9V3FQotOAeGbErQEFlqgO0gRWmhZSOo4n89Qa2j5HdBRzudzOyBs2Ft+j1Keb/YS9nA4/M736/Xa2kPho8Pv7yp7RJWR4FpzAsqIQdvUGTyhSvc+Hrhs+zg6ufTsidh7uqxBarrseibvqpHpZgYMu2ftJ0/+i71oEseExhRjJzSmGDuhMcVMCTMQDbPKMilXhnGC2FC3RFIEZ7RAoBJWatLTLaSUgDBzOp3E2bvdbufzuf3r6XSC411+p7W/3++hkAOyENhDnSs+HFjA7XaDE/LVXF8NqE8vAPceJAx1d0uK3ZsLabli6s2p3E/ITbDPjs+faDmUm8Keu6fqzGc34LHcl91PGO4P5D274/ZLL40O1hPuaRbdQX1l9Bdo92yPN1+MH0eNKcZOaEwxdkJjitkxbW3EYFKY4biibUJQoYUc3vvDOlN2K1MrFYwIM2wg7CFDZYmEGUiIeabgiOWtEGZEFMevpOeUnVSc1o2B275655fWaZ7/1orAhmz8fkKdhhYKM3DY7EXQzAbfQghhZYLz1FJOuNDuGN5xJ8TbJe+EoItA2lrrcgvtD+zaayeEw4EsPIDvDT69WSdsmyAdQ5PHD8tqvVKq8eOoMcXYCY0pxk5oTDH7/lgPMc+yszCzuT1rLbA8CDPAWCf0aOVgoa1MEMbcbjcIGiEDhseHAUN7GLztfqAaMzrEPdBOJa6po68+wCd/W2FGR4nbBo3vlTGTTXDRBlypqYUTo0J1VN9GcBfyYkC3YD0Q/spKBmzhg+lgeTACuBmLLpzWB5WgtCy0kDIE30H66rC4yjKPkMG47BUIUZDXxkIOnEx22uVrdIZdiB9HjSnGTmhMMXZCY4r5ToWethVmRqYLY3FtkFpPOB0biL0//Go0vSsifMd9V+qA5Y3bh/lJHLFzDKyFmZmzfYgyZr6NMNMVWsTeJZYueECwb5ucccKfiNGelmAvpIXu4Pou5zW09tfrNZwOzl47JtdlgrkOv6fsgDLxvMVFis/pdLper18tj9dzOp3a7pwW9/n5KfLgQIVaepWjwCtAvP34+GjtL5dLa3A+n1vlBoScp5Ic5rV91czyvhW4N/ymWTfdYToLLzW+Pt7wK6z7xa+zAnf9j5qzAlPTscFCaYB6NOibtdfdt8UxoTHF2AmNKcZOaEwxs4We2qaOoJ6I6BY+YX0MpuANgWCpQ/nuaiGHA4Ic3R1iMOg+KcywPSgfUGepu9MKhBzIYoP1sz0IM7qgONtfLpd2usvlIjJmWJjR+9pYHW2nW5bl4+OjNbhcLq1yc7lcWJiBylHQfFN1NKtkhKOlpssqE90uPMJXBodeZaTxQkxdA6HNLstyOp2EPXxnsRNyqifc9JCFdz6fWydnbXYmd5QNFvl7FR8v/xW+MXVe28h+Qv0fBgwu/hrix1FjirETGlOMndCYYvYt9ARNzjyC7iIP69eH3X937cFAZ2B0M2CEEvOgukywfi6stBDZDBuuCtNOp3UjttfThcIMB11gDxkzn5+frT3XoYL9jSlh5kBJc7CV6ePjo71eHx8frXJzuVza5rO7KOnNzWWCl2bMaGEmzCDRn/DgHMproQX8MEyLW0iuFGl33aw6Yc86DdizcgDfUDyjuG+6WXtiP2H3OkLamk5z0/sJH7/vZly+OP/taCDMwFfk5XJpDxnS1roup4Wcr1ayAj+OGlOMndCYYuyExhTzukJPSyTMdH+fbZtd8aDtzntzxLN7N8ATPx93a0iL9fPRhUoSTKelI5AiwhiP7fV0y++/70OUxTutDrS1Cuw592BGmAE4QgZhBmJCUGIgRHymIkAOzfjOpiyv20/In7BMEipg49M9jbXQwoPoDBiQCnhGuI3CNDcxftdFRdpaWNUfdKPb7dYKJ5wxE9r/+++/7fI+Pz+F/f1+b5sL7SecFGb4VgFhBtLWfvz4IdTR4/EIiWyg3EB6kDNmjPne2AmNKcZOaEwxszGhyHznZtY+FROO2EOOiBZywCC01yk73b+KI4K9QhyzaR3lQcV/U9P9GuFXE/o+X70m7KFLN6FHROAs5GhhhlfIO63aAbkWMJSo0SkycOkrM2Y4JUU3FzrvOtSG886bd0R3XgAnxLD+JtbP9py31dp3vai1h1eRwcvGHo8HKBnX61XLiSCcgP31egXhBKYDoYXttdDy77//wvJAN0odDru0TgPUWWzLslwul9aG1dG//vqr7f7z58/2erFOAy69TPBGj6PbpgLNE65H/8dbezjZxY902bD7CrY9vdlrt+vRvZETGvNnYic0ppjKCtxZYUZLC6G9Vm5WLG+RjzGhUKSfeUa6Q9JG2D21AD1+dj3Z7l2D5WtGRpt8BN0v3HhdjZnuolMJNHpvUXjXhkVZYEBWF3UiFaiXYUlpFmbAHpQMbT+i64RCjrAf0YG0PSwPDkcf3TK2P/MX7HIgpXAFbhBm/v7771Zu+fnz548fP341Pz4+2u5QVCpLZdra/IAz03XF2xaWQ7VXh3/Vgh5/xYxsaBTjgz2/go9/8BBfaoffq+J3lW2dhcdfkUJ8flDaWneD5fIFXSdsP4E0NHDRQ+8nCjgbsLavVjKCY0JjirETGlOMndCYYir3E4IBy486z4tjJBg8rOzEwgx0Dw10ExIptDDzoFfSQ8lqiJFYyGHlA6YLlZLWPpWgw4ejp+va67OhryYDNwOExJy21uouy7L8R4QZhrUQ3eQRQoOWMKsOciPBzVgJYJ1Gjx9q3Fp91coE36bwnQW3NewPDL0CytqH6muYhaenC8VkFsNTUnyYSgrV1mD/IaupIPN8tZIR/DhqTDF2QmOKsRMaU8yO+wmZ0GByOv1rNRtsK+SEQUso5EB33vqk7TmmgqBLx2B6/LAGTBjjsT1MB/b6ZOqYkNHVgZ9FZVr7v/76q73c30aYCRNiwqyR1PgrptPNMCeD95WGN4roHt5nutIRN7VXj2ThjTstLy/8UkjZL2Pl7UQT4FpvwicPh8OPHz9ag58/f7aeBk4IymoWP44aU4yd0Jhi7ITGFPO6/YTdT8b/Oj/diJCjDbbVgUJdJ6sDZYWfVIjYFX5S0838Fh+e2+4nLby/EbZNQEz48fHR2vOr1KDkjJg65NVbmcAGhBYeQdgzoVyWNZi0z943Kfvsfcn2KZ1psjvbh+vX43MzBVc3BxcFv+LtiK0Y87///W/1ShY/jhpTjp3QmGLshMYU89KtTAA/9GtCe22wovvMgFq3GDHQMRL/eh5OB7siwvEnu4sMJP4rG2dPYApdvWKh399BmPn4+Pjnn39+NSeFmamiL8aYefw4akwxdkJjirETGlOMndCYYuyExhRjJzSmGDuhMcXYCY0pxk5oTDF2QmOKsRMaU4yd0Jhi7ITGFGMnNKYYO6ExxdgJjSnGTmhMMXZCY4qxExpTzP8BkQQaOZyGbRYAAAAASUVORK5CYII=)
%% Cell type:code id: tags:
``` haskell
-- write out matrices to visualize with python
(xl,yl) = unzip me
xmat = reshape 49 (fromList xl)
ymat = reshape 49 (fromList yl)
-- create x and y matrix
saveMatrix "xmat.txt" "%f" xmat
saveMatrix "ymat.txt" "%f" ymat
saveMatrix "zmat.txt" "%f" loes
-- ../vis.py -x xmat.txt -y ymat.txt -z zmat.txt -t image
```
%% Cell type:markdown id: tags:
<img src="../images/poisson.png"></img>
%% Cell type:markdown id: tags:
## parallization
lets do the same in parallel using repa.
%% Cell type:code id: tags:
``` haskell
import Data.Array.Repa as R hiding ((++))
poissonOperator :: Int -> Array U DIM2 Double -> Array D DIM2 Double
poissonOperator n arr = fromFunction ext getElem
where ext = extent arr
get idx = if inShape ext idx then index arr idx else 0
getElem (Z:.i:.j) = ( - get (Z:.(i-1):.j) - get (Z:.(i+1):.j)
- get (Z:.i:.(j-1)) - get (Z:.i:.(j+1))
+ 4*get (Z:.i:.j))
* fromIntegral (n^2)
rhs :: Int -> Array U DIM2 Double
rhs n = computeS $ fromFunction (Z:.n:.n) $ \(Z:.i:.j) -> f (toCoord i) (toCoord j)
where h = 1 / fromIntegral n
toCoord k = h*fromIntegral k + 0.5
f x y = x*(y^4)
```
%% Cell type:markdown id: tags:
instead of HMatrix's `linearSolve`, we need a solver that works with Repa arrays.
%% Cell type:markdown id: tags:
## conjugate gradient solver
CG is a simple, efficient iterative solver for linear systems $A x = y$ with positive definite matrix A (i.e. $x^T A x \geq 0$ for all $x$).
%% Cell type:code id: tags:
``` haskell
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE ScopedTypeVariables #-}
import Control.Monad.Identity
import Data.Array.Repa as R hiding ((++))
data CGState sh = CGState { cgx :: Array U sh Double
, cgp :: Array U sh Double
, cgr :: Array U sh Double
, cgr2 :: Double
}
```
%% Cell type:markdown id: tags:
`Arr` is a type alias for convenience, `CGState` is the state of a CG iteration.
`CGState` is the state of a CG iteration.
CG takes a function implementing a linear operator (type `Arr sh U -> Arr sh D`), a right hand side and an initial guess and returns the (lazy, infinte) list of the `CGState`s of all iterations.
CG takes a function implementing a linear operator (type `Array U sh -> Array D sh `), a right hand side and an initial guess and returns the (lazy, infinte) list of the `CGState`s of all iterations.
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`.
%% Cell type:code id: tags:
``` haskell
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::Int))
scale a = R.map (* a)
-- the main iteration; maps a CGState to another CGState
cgStep (CGState x p r r2) =
runIdentity $ do
-- need inline type annotation here because computeP is polymorphic in its
-- output type (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'
```
%% Cell type:markdown id: tags:
Using lazyness, the iterations can be decoupled from stopping rule and output of information.
%% Cell type:code id: tags:
``` haskell
takeUntil :: (a -> Bool) -> [a] -> [a]
takeUntil _ [] = []
takeUntil predicate (x:xs)
| predicate x = [x]
| otherwise = x : takeUntil predicate xs
process :: Monad m => [a] -> (a -> m ()) -> m a
process xs f = foldM (\_ x -> f x >> return x) undefined xs
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
```
%% Cell type:markdown id: tags:
- `takeUntil` is like `takeWhile`, but also returns the final iterate (why waste it?).
- `process` runs over a list, performs a monadic action on each element and returns the last element. It does not retain the start of the list while iterating, so it gets garbage collected properly.
- `runCG` executes `cg` with initial guess 0, a stopping rule based on residuals relative to the first initial guess, and outputs the residual while iterating.
### Example
%% Cell type:code id: tags:
``` haskell
solution <- let n = 50 in runCG 1e-3 (poissonOperator n) (rhs n)
```
%%%% Output: display_data
%% Cell type:code id: tags:
``` haskell
import Data.Packed.Repa
setPlots 1 1 >> withPlot (1, 1) (setDataset $ repaToMatrix $ copyS solution)
```
%%%% Output: display_data
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAASwAAAEsCAIAAAD2HxkiAAAABmJLR0QA/wD/AP+gvaeTAAAb30lEQVR4nO2d23rbPK+EKVtO//u/3XqzDvp8XcxAGRACJTjpvEehRULUBo4wBqHl9Xo1IUQdl+oJCPGvIycUohg5oRDFyAmFKEZOKEQxckIhipETClGMnFCIYuSEQhQjJxSiGDmhEMXICYUoRk4oRDFyQiGKkRMKUYycUIhi5IRCFCMnFKIYOaEQxcgJhShGTihEMXJCIYqREwpRzFq14+fzmemQHG47HL273WPtVt7fLSSbPNJk/4lcLuxfiN0KnyzLAlv7DtDZbl3XT77z8fExNusNlkOL/3Lj7q7PHH7oZL7RVEc67O48F/AivpU325bXQf9+yLIsfYfn85lxQj2OClGMnFCIYuSE4hP/yLtJ5j6TJzlQmHEjGd4h2v/Mfb1er92T2XFcockfet5CW0c6hLCB3HjgB0Fd2woC+w62eb1ee2u9MHO/3zMx4WQnPPR2eR9rrhP+mCPd7DC+daTDIDs8kDjhsviSJPTvdeDL5dIPT0rEehwVohg5oRDFyAmFKGZmTBhVL+yTNMRd3HhoX9AhaQ1m7g5PHmmmw9zT6FrgnaOEfnB3g0CuxFyvV2j2Sszz+eyVmNfr1Vu73+/Dx7RB1gn55T9T/HBv9JDxkDV3+FxrUSec/nVGmpbdfujmuITUTnDCTWEG+sPW/iqM6Drj6HFUiGLkhEIUIycUopgDhRm7NRPjuUEdD5wKrW0eOO9/2lRH+pOm2yEUOLmrIjJKzLIsvdDSWoO1SKC1vF4v6A9774cnf6w/MGMmKQlEb76kY/yz1vhV2Oww3tz85CuiTsgXH9kVgKCsWPvEl7iuo4wZIb43ckIhipETClFMKiZ8Pp/8OT4UP0Qjn4mRkmvNbUajstDc+FSjw/ncMiFiJiBs8Z/juTBjY0IQWl6vF9d1YNfHraKYmTGz4wJDc6LbQH/3vk+6Wd9/rjU7/PF4kCNNfr+0qdk8IUJZaW1gQSBxQldrtcYhi41/+4TQ46gQxcgJhShGTihEMWU1ZkYEgIzAAE/t0aDODau48Yy1HcPHrW3275vuVSBbbYcoPCyM/jpPivk2E0Pe7/c+6ns8Hv1We8V7U2+UMcMvSV5KyYgfUWsgfryztdB5a7Pzb1wnHHfLUFZaaw3UTrsgkGfM2N2B/gnW+tMOnZUxI8T3Rk4oRDFyQiGKyWbMkDfjhGL0aLDhKgShgHNEtgkNJ0HdSMQYGh4NhjOBevIyWUJKDARyr88pLy+6Fgm2ti3pZVyJmcuBGTPu1olSyuv1Cokf/EZvxjEej8dua0fP7VCZ51BhJlO4qXlKzIgwA8OJMHO5XPoTBdaShZ70OCpEMXJCIYqREwpRTNnrsqO/6iZ1mlCokxyel3mSuk4mRIzqNLC1fWauMOMuZeL74rMFYSZ6pBnOqzHjXm97K4+rHSPiB6Q4cO3EFWbGrW0aDzlhdG4ZYcb9BikUZsAJQZhZ1xWEmV4dvV6v/Uz+6C4gvfRnxqo+vfqijBkhfhRyQiGKkRMKUcx5xX+T1qKRcaj/DpknOXy3ShS1Fo3xMsKM/cQ97aEMGAgao2fVdiZazqFKDJBNW4PlJD2hW9MKBvxmej6fRPwY0UKiyko/PGoNMipC1qKHBs2kMMN9uM0WZnjlJVtCGxYE9k67rmtfJNsqq3axEggzUFSGKDFJYebAV6Pxpr2cmeFRa/xmcv/bnGnNPdKMtbZ1lpLD22cm/g/JzBz6L+bdZplLnEQxoRDFyAmFKOZ9nXDiv/uTrUX3dai1qPG3Pe21Z/VQji301B/q3GSukXyafrhtjosfrwFlBeY21xpPoLnf7+PnzQ6PZsxY/Yw0Nz/5SyhFZlkWOG8gzKzrylNkYF+w1glOI9R94jcnnMMox+aOvq1c4Rp3pXnyFcC/IKLWeLOlVydya7VOSLa2rWW4JA3N1nGyJ1bCjBD/KHJCIYqREwpRzIFpa6Fo4fUffz/hSsyIMNPvC8LupJQSteYOP07mCZ23Fq8TlYmU3JiQl/d9fq4zZiNzOKsg89xuN0i46fd1v9/7hBt7UfrOxRkzwG4txH4Sbb6PlLLjRidex51wxFpGHeVO2Lwq+pzQy8msMAMvCbTW7KH1Y59mFa97D/R/S5gR4ucgJxSiGDmhEMXMrMDtBoQTY0I3dLFxORc/eAINV1bu9/u4MPOnOa7EuMKMmxtEIsZoiDhdmOFrl0BKAWHmZRYcwqHBxKBIzO1264fbgty9NbtADM7qyPF+xXlLmfJOmMkjiUopUfEjJKWE5FBXa01+vxQKM271NHBCGA5JMFaYsXXv+11HBS1oTnRCPY4KUYycUIhi5IRCFHPeO+vdziElxm7lSQw2ziYBwIh6AdYyUgrXacD4pkoEzVCoA8bhvLk6Tb/VXnH7CSn4a4NAHhO2z1GlrVhBYsI/Bw7CTN90FSk7md2ct5RprjDjZqVwL7XWQnJF1AkzyxGfz6cVV3trrteFKmKNyM59/5AfukXvIWMmKszAvuxphwOBI4XTDqeR/w9oCfQ4KkQxckIhipETClFM2avRgM0QkTejPyKTUGckEBpXO0a0EB5hhuLVUDWdaPR7dEzI61PY0tLQgTQfj0c/3FYVetKX1HMl5n2FGVdZcS8n6RBqtnNzXKwwE5JSeNaba+33798Tj/RMJ+QZM7ZIDFh2lz4RGfnPO+jBCftTYa8p3G8T/VCPo0IUIycUohg5oRDFvJEwwz+JCjPRADUUCLm6TiZe5TrNZjAcEmZCoXVoZVMbUCxgLRJs4u9CgyDQprwQYWbznM9KkUn+WD95PSGfa9RtJgozh+a4zBVmotZAmMmsTmxxYYbrNEA0YyYkzLhpa1aY4SeK3GBwWpLocVSIYuSEQhQjJxSimLKlTNGA3tVpoJnRdZIyjxvih4a71kKikXsgJ8eEpJgvZLTAmonN/kShGJGUiBITupmjnLqUKaMZgoSQVCCTGTNzpZSotd+/fxNrSWHmYQpDHeeEIJxYYaafm9VdrAXoAEosz5ixTWXMCPGvICcUohg5oRDFnJcxk3yGDg2fLvPMtZZRfaLDQ8LMlMm0YZ5mMRERZqxxd27uWJK+w5kYELYzX43mdubKCr+3uIrTBl4YFk2gIVLK4/Ho+1tlhcs8rk7TCzOu6uOuJySFWJpXedp1QqKUtK21SLxwUzPrCUkRGniz2stU4L7f75D1Nv6SephYZdqaJfkfJjP8TGvRmy/0FZ6cG2+6uw4N3zG3HnAzMG4LNyUnw60lmxkUEwpRjJxQiGK+jROGIsyoteSTRiYYnmv85N25xjPDQ0+2UeNRa4fyLuro5vM9CQCiqWFWiRlvvgZKM/HyJPzVaFaYOTT/JpoxwzNF3IwZG4n9hVdqGlnKRIQc+8p7mzED1xR0GrJ2iacNRTkvbW1k69uKH4XWRjqP959rzT3SzU8Gx47sK3PgjZ7n0NYk3+ZxVIifipxQiGLkhEIUU1boKfqEzeOBTZ0GOoBcQZpRYQaklB0FbELWojIPz5ixTXJiR1blta/JZMws//FVf9B17vf7uq79cHukkCXHbwkbUvZbvz5on1OdcKJcEbW2Y/isub1yr23jzRas9bZjHSY/luR6QtLfqqP2vMGu+YHzr+noNZUwI8TPQU4oRDFyQiGKOa/QU6jZPOllJHQJDc/EXVGZJ2rNTXkholG+zPGhMSE0QWgBYQb6X69XUKRIkZjNOlHQdGPIvvnVMe4gW4H7Z0gpI515/3EtpAW/ApLCDG+61kb23jdDTujiVluDXUePFIZHb4nNv3egx1EhipETClGMnFCIYsqWMo0IM7zJhRn+TM+j2aQ1N8SfK/PwtUjJ3CBXAOOBE2ALN5FSSy/zajRbnxe2jmfA2NmGlJhkEAi8b8YM7z/SOSSlkK3t+KyUQ4WZMytwQxPgTmjF0qS8Sc6btebeEuR+40ftosdRIYqREwpRjJxQiGLKXo3mdg6FwtG4OaQSbcYDZ8o80SBtXDTaMXyiMLP5SW/N/lgPFbuhSfSwzfMAw8eVmO8kzJC5une2HT4utLgdouJHyFpSmIlKKYXCzOaJal/jqqP98B3CDFmHCdY2h49fRPu9aY41gB5HhShGTihEMXJCIYp5l+K/bv+5Mk9ya0bXcYWZkE6TtObOLSTb2GPnPOm70GDrZgeurIyfh1qOfStTdGtGmIHOoUvClZVNa+O6jmvNaifcWjRPbW4Fbq6OhIQZqCIDCwihcBNUamqt9XWc7KFxoeVPxRrSwR443ADkMKPocVSIYuSEQhQjJxSimPfNmHlna4eKRhnVx+2wYzi3lhRmyI/1m3E+78+nSra6HKrivG/GDL9ZR+6GcSmFB+UtKMwkrfEcFyvMnJx/M1eYgbFwWkCYgVedue82Oy1jJokeR4UoRk4oRDFyQiGKOc8JDw1t3X1l9l5rLbM1OjwqtETPQ8Z4UgTKjHWHJ1dRlL0aLUpIxNv8ZKI1uCF4MheXYUasgWDArbmVmrj8EF3ZxOUxC2gtvZQCY22KzJOmrT0ej1654SWtmqkT5V41flwZKgs9RfvnPUfWzrfWY185GGpOn9txRxpCMaEQxcgJhShGTng2SZnnTH3L3d2hk4k+PSaNnzkceN931s89zuTTv1UgosO5tZCuE7UGK5tCMg9Ys29KCwkz0RozsPiICzNcc9o80tByRHJcSSrVUX6znrC78WZyuHukE60ljecn04YB6WXEpcnu+PCkNdt5InocFaIYOaEQxcgJfziHPuGLKfyQjBneOb/ruSGctZ9purEQORbbDOXf2P6hM8/XBz5NLSb+ajRXibHG+aGR45obLr61E44f+aafhE6NG8Tvnttmc65XhxSFjFzhDg+dcy7MwNbkZEbukDMdr0ePo0IUIycUohg5oRDFvHVM+I9waGRi5QduLarrkK0uyUJPIV3HFWbcQyPNJHLCQyCX0+2/o7lbmHGtNXMsE+/OQzNmRraG9j7X8Xr0OCpEMXJCIYqREwpRjGLCb8ahgYobgyUnY/uHhBnbf3xfLRjcujPvm/9KoSfxl6MdL9Qh86UAnblO4+4uJMz8+fs0CYqjx1EhipETClGMnFCIYt4oJoQCJPYTqBtLmhDBR/durY3MllgL7Z1PZmRux3Ho7tyrxoUc29kGgV81Nz8Z35rkjZwQIPdldGx+JpuFaCfaJ83Q2BNwtRPS3852fKu79x0ZMCFh5jj0OCpEMXJCIYqREwpRTGVMeLLAEBI/ksOjotFECQo6cOMjOShk79Z4rXITMjXFzhTKnDCjRpTvLi+9EkFirjAz7pCb1jYLvRwH0TPJh383RZN7MjIP2VEUPY4KUYycUIhi5IRCFHNsTMglAd6ZD29UMIiKH0lrUWEmNJno3JLCTOhIubWvDBLOzz14Bw50wqQWEpUQ+NYdgkT/tzucOyHp3Fq7XC7jOg3fl/0EXjB2uVz6wtUjDjwuQtjOGV0n6ZCHSilgMLmeUI+jQhQjJxSiGDmhEMW8dcbMacKMK6WcaS00PCnMNHMVThZmQrH35ifE8neReb5NxkxSrgDxw74P3XUDGJ60xud2ufz/E8rr9eqbbeDG5cIMGAdroDHwQ7OTgeFcmCHHtdnMONU7O6QeR4UoRk4oRDFyQiGKmRwTjkfhcwWGaOYHb1otJDQ8NNXoZLhOYztEJ1MozPDmn7/HT5S7r/fhPGEmdNeODLdXiFhz5YqodgJbibVNoYVYu16vkOphd9dvhf7wJnfYVzNlXUBKgTeKXa9Xkg5iLYdyR+xFIc1lWex5IPeAu+vxzkfzXR9HzzyJ7j+ud7b2zifqrTyhkO/qhEL8GOSEQhRTuZTpTPEjJFe41qADN+5KKaHJzBWNklfB1cMmCjM7+ketVfFGxX8XU3kFmlyJcS/AuDADzTYgzFyv179NW+wkJMzAVtBdlmWBfUGH6/Vq/eorlq0XIZHaKiDkQHPTCcnem/l+6U+FveJcyHFlHtjvm3hgO/RxNPktxa2F9pXk5O/j0LFEz9vc/xWHDk9a+0YoJhSiGDmhEMXICYUo5o3UUS5XPJ9PEsQ/n89oHomVN3rjfbMZdQSsPZ/PkLV1XaFCCYlnHo8HaEighYA1nrZi82PgE67EQH+Y+aYWSrQZe3uAMAOnkaciwXC7lWt1vHko3yZtLWktqk8cqihw6dW9G9ytcKeSBYFWBw7dqXbmoXpKmV271kL7cjnUIfU4KkQxckIhipETClFMKiZ0VwABXJhxQwJSecUKMzZccYWcr7a2oDDTTGVYkFLAmpU32meg87qu/aakMNNbs7sLCTPNxE7RjBlocmHmer2CEtP3t1ttk+g6oVAc7pYo51Xg5lvdwDeqrNjhRK5wv024lBKyFpVSRlbZja9OhGOBG9EeGv9GsDefXWEYUkfdIw315013MqQ5Fz2OClGMnFCIYuSEQhQzMyYMxbKbwgwJdjelF76V59/wArt8uRAIMzbyIRkzVkrhwgz8Ap4UZh6PB1mptDmB/sw8Hg+YqruUyV1L1f9NgsBlK42JCDMjustuYWYubyTMLDTfgkfhI8pKRvywwgw0rZOT/uAJPEUOrC3LAl8BYM0VZrjH2vWEtqj2pim71VprRnqFsdwJ+WnnFzGqvfHjmuuQehwVohg5oRDFyAmFKOZdljJtdh6P2pPCTPscmG0uTSJSis2Y4cIMtwZAyLeZf9P3D2XM8JjNDnGXMsFVcI3zOI2c1TYgzBCdZjNjJhQ0TgwLK5cyuaEwd0IedkfzSOB6hMQPyPxaTIkquHugyjV3wqg1cEI4USNOSGQh+xVgTxS3Rtihjvb9wSejaWtuwg1MtW8m09b0OCpEMXJCIYqREwpRzMylTC2Siu5mzIxIL2SrDRr7mVhhhlhrW8smYPIThRnXmpsxw2WevmnjSRvFwZkJ1ZixIehXE2tjP9bDNYUg0FVioDmeMeNKGBkmp63xfAjetDc6JLXwrTwrhWihzayaW8wLw+zryogCuZgX3IeyUvjXx+Z9TKQUe6fe7/d+oNUYIJsHvK534x3CzHhSlCvM3G63vv+6rrCe0F1eOJ5wc5w02vQ4KkQ5ckIhipETClHMse+sjwoz0Hn8x/rNfU2UeUL5N9za5vBxYSaazQMs5jVM7s/3xJr7Q/zEmLB57yHf0RwXZtp3yZhxhRke+NolM9CfqHavrQWBff91XcF+LyHAVKN5JPzQdkgp/dzcpY9cB3o8HmANfJifN6j/ncyY4Qk0m1pd3+TCzO126zvcbrfb7fa3ua5rrypfLheecOPm05CZR9HjqBDFyAmFKEZOKEQx5y1lgsDJCjO8f2jrZjMk8/D8m0NlHhtHEeNRa/bQbIjI82/6Jk+Isfk3XKRxY0IbDI8rMZu6y3iKjHvDZKisMWNva3IzWeDuiap84zdT29KcuJQCWSmu23AppbfmqqP3+x2sEeOv1wumer/fQZgB6ZXrNET9ap4TNu97FpaMretKhJl1XUGY6ZuXy2Vd1/66gHLjJtz0M3HvVc4/+jg68Wvs5F1HdblQf7s11N8dnhH9p1+yk3dH+EedUIj3QU4oRDHZmHD8ASMqzMy1tkPI+cq4teZG7cdZ22wmrcFso5oTHCm/TLuvePRIbXPE/ngzQ3Y9Yd9cjHpBDvJiXuMMAoOV6cZn4iap8JVQ7jI5OFIrpdh7sR/uSilc5uHrCX///k1SXqxx0Gl69aIZYcatOgVNkhjUBmoH2/PWfwJNWMoESsz1ev34+PjbvFwuHx8f/WWyGTYgzPCbuSU4r9BTFOvSsJWPtcPHrW1+43JrVf+77Fb7/cK/feCrMNT/Qt9Z/6dJvlijF8V+6Wf+E2b+s038N9gUEwpRjpxQiGLkhEIUc+B6QvsUzhOmbDoFeWq3NbBhbPsc9y/m3etRla+3D5VXrPhhk1r64ZCVYpukrIuVUiCPJDq8V33+DO+byYwZt7hOo5D7p3nCzPV6BWHm169ff5uXy+XXr1/9cBBmbrdbf2LtSqg3FWYW791mfKwrfoSC+LniRyg/0+q6Nsup7wAuDR7bzLeJlRzBGkz+q+OyM9/cHV85ydcT2lPhZq6RmXMn3FGBG46F3yFkbkn0OCpEMXJCIYqREwpRzIEVuDd/R/77t43RbaBFrLkr2ezEuMDgrtnrrdn1PrxCLuS4cJ3GtQbaye12IzKPTXkhxjf7h1JkuDATKsgNLP/x9xMI82xMCFJKL8wsywIZMx8fH31/0GlseW8QgcYPxHKeMOOOXWjBInsBYF82jCbiB2x19QPXGgni3Ry6Rt9KHy0btZi1jmTpY1Q02rGAcLcwA9jvdFeYIQW5baEnXtmJy6FaTyjE90ZOKEQxckIhijk2Y4ZnpfC4y1rjwgwEDKAo8Mm4BXa5MOPmkUAYZgvsci2EvAjpT3OWtc3+E4WZTEzYtiIxrtNAENgvZVqWBTJmuDAD6Tg2gBw/CsuxhZ54youb9QLNUKIQv94888PeGXbyIa+zKjHPmBlXbls8s+w4J2xbcii5CptfxF/Bxfa2pdOA20ChJ/5mNZsER7wumUCjx1EhipETClGMnFCIYmYWerJPxjzM4ykyjeo0r63KTgQb2/CiMjbugiYPnGApE19qkAnq3Mm4KS98bm5MCDFeKEUmFBM27/6xYkl/FZZl6WPCTWGGvNQJ8m94iBhlZqEnS0h6scuFiBayqZ1Ak8hCvHRK85bk2AybuXd2yFpSz3TnFnKz0Dvro8IMfBISSy+XS+9Ff5ywH/7r1y+o+wRiKan7pIwZIb43ckIhipETClHM5FUUvINbxwXUEYgfelNWmOF7X8y72mHvECmBNRsahYQZ/gIpCOp4GoqbpMIjxqg1N0TkEaObIjOeQBONCRfz4jQoxvPx8dH352/btj/998bBcpTzXo22mGoldqUSv8DWMUjTKjFcEsgsyYnKFaFb2XUz/o0w4ma757Y5nMytTXVCVywlOg2Ipa01vrwQmvBateR6Qj2OClGMnFCIYuSEQhQz+dVoPa524go5XKcJ/ULKYxU7VRjOA0i+0IGLQJvWeFQGEYg7nOe42EPjex9v2rm5J5bEhC39Xh1IoLE6DfmxPrSoIko2Y6Y/a9ap+FabthaK2uECh6y5+wLjfHj05kveyu4XyplzC523aJMT9UkQZuC7bF1XvrIJlhf2Lvq///1vfNoWPY4KUYycUIhi5IRCFDNTmLHSy3HGrTBjgxM+nEcjfFFVSOaxw6PW4Eh5UBeKGDfnBv13b93s0D5TFRM28ws7xIRu2dI3zZixq4dsh77pXhK+NdTZHX6oteS9+D7GRzrs7jzSgZBczupW+yUVa6wcFeLHPo4e9z/5za0ljb+ztYm7dj02tDXJj3VCIb4LckIhijmw+C9/ULZbT+4//hy/w1To1/bQ3KJjbZQVWgYRnQzg/pQfNUiwaSu8OihstQsjQHqBBBooj7J72q21/e9REkJMQY+jQhQjJxSiGDmhEMXICYUoRk4oRDFyQiGKkRMKUYycUIhi5IRCFCMnFKIYOaEQxcgJhShGTihEMXJCIYqREwpRjJxQiGLkhEIUIycUohg5oRDF/B9aBv6QzKY/3AAAAABJRU5ErkJggg==)
%% Cell type:markdown id: tags:
## Stencils
Operators like `poissonOperator` are a frequent pattern: each element of the resulting array is a linear combination of the surrounding elements of the input array.
Repa has special means for constructing such operations in 2-d: Stencils.
```haskell
makeStencil2 :: Num a => Int -> Int -> (DIM2 -> Maybe a) -> Stencil DIM2 a
```
%% Cell type:code id: tags:
``` haskell
import Data.Array.Repa.Stencil.Dim2
poissonStencil :: Integer -> Stencil DIM2 Double
poissonStencil n = makeStencil2 3 3 getElem
where val = fromIntegral n^2
getElem (Z:. -1 :. 0) = Just $ -val
getElem (Z:. 1 :. 0) = Just $ -val
getElem (Z:. 0 :. -1) = Just $ -val
getElem (Z:. 0 :. 1) = Just $ -val
getElem (Z:. 0 :. 0) = Just $ 4*val
getElem _ = Nothing
```
%% Cell type:markdown id: tags:
Stencils are applied to 2-d arrays using
```haskell
mapStencil2 :: Source r a => Boundary a -> Stencil DIM2 a -> Array r DIM2 a -> Array PC5 DIM2 a
```
- `Boundary a` is a type that handles the boundary conditions. It can be
- `BoundFixed x` returns a fixed value `x` at the boundary,
- `BoundConst x` assumes pixels outside the array have the value `x`,
- `BoundClamp` extends the boundary values to outside the boundary,
where `x :: a`
- `PC5` is a special representation for stencil results; it is defined as
```haskell
type PC5 = P C (P (S D) (P (S D) (P (S D) (P (S D) X))))
```
(No need to understand this.)
- `PC5` arrays can be converted to `D` arrays using the polymorphic function
```haskell
delay :: (Shape sh, Source r e) => Array r sh e -> Array D sh e
```
An alternative way to solve the Poisson problem with Repa is then
%% Cell type:code id: tags:
``` haskell
poissonOperator :: Integer -> Array U DIM2 Double -> Array D DIM2 Double
poissonOperator n = delay . mapStencil2 (BoundConst 0) (poissonStencil n)
solution <- let n = 50 in runCG 1e-3 (poissonOperator n) (rhs n)
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
### Stencil syntax
There is special syntax for stencils. The `poissonStencil` above can be produced more easily by
%% Cell type:code id: tags:
``` haskell
{-# LANGUAGE QuasiQuotes #-}
poissonStencil :: Stencil DIM2 Double
poissonStencil = [stencil2| 0 -1 0
-1 4 -1
0 -1 0 |]
```
%% Cell type:markdown id: tags:
The syntax is processed at compile time, so scaling is not possible in this version:
%% Cell type:code id: tags:
``` haskell
{-# LANGUAGE QuasiQuotes #-}
poissonStencil :: Stencil DIM2 Double
poissonStencil n = [stencil2| 0 (-val) 0
(-val) (4*val) (-val)
0 (-val) 0 |]
where val = fromIntegral n^2
```
%%%% Output: display_data
......
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