Red-Black Gauss-seidel with mir

Christoph alt.c4 at web.de
Mon Sep 14 09:50:16 UTC 2020


Hi Ilya,

On Sunday, 13 September 2020 at 19:29:31 UTC, 9il wrote:
> More details are required. What compiler and command line has 
> been used?

I have tested it with dmd and ldc and called them just with
$ dub build --compiler=ldc(dmd)
with no more configurations in the dub.json file.

I have compared them with a 100 x 100 example problem and the 
version with the for-loops from above takes around 1.8s compiled 
with ldc and 0.8s compiled with dmd.
The slow version from below takes with the same problem around 
18s(dmd) and 12s (ldc) on my maschine.

The driver function for my Gaussseidel method looks like this:
```
/++
This is a Gauss Seidel Red Black implementation
it solves AU = F, with A being a poisson matrix like this
         1 1 1 1 .. 1
         1 4 1 0 .. 1
         1 1 4 1 .. 1
         .          .
         . 0..1 4 1 .
         1 .. 1 1 1 1
so the borders of U remain unchanged
Params:
     U = slice of dimension Dim
     F = slice of dimension Dim
     h = the distance between the grid points
Returns: U
+/

Slice!(T*, Dim) GS_RB(T, size_t Dim, size_t max_iter = 10_000_000,
         size_t norm_iter = 1_000, double eps = 1e-8)
        (in Slice!(T*, Dim) F, Slice!(T*, Dim) U, in T h)
         if (1 <= Dim && Dim <= 3 && isFloatingPoint!T)
{

     const T h2 = h * h;

     foreach (it; 1 .. max_iter + 1)
     {
         if (it % norm_iter == 0)
         {
             const auto norm = compute_residual!(T, Dim)(F, U, 
h).nrmL2;
             if (norm <= eps)
             {
                 break;
             }

         }
         // rote Halbiteration
         sweep!(T, Dim, Color.red)(F, U, h2);
         // schwarze Halbiteration
         sweep!(T, Dim, Color.black)(F, U, h2);
     }
     return U;
}

```
And the full slow version looks like this:

```
void sweep(T, size_t Dim : 2, Color color)
      (in Slice!(T*, 2) F, Slice!(T*, 2) U, T h2)
{
     const auto m = F.shape[0];
     const auto n = F.shape[1];
     static if (color == Color.red)
     {
         U[1 .. m - 1, 1 .. n - 1].strided!(0, 1)(2, 2)[] = U[0 .. 
m - 2, 1 .. n - 1].strided!(0, 1)(2, 2);
         U[1 .. m - 1, 1 .. n - 1].strided!(0, 1)(2, 2)[] += U[2 
.. m, 1 .. n - 1].strided!(0, 1)(2, 2);
         U[1 .. m - 1, 1 .. n - 1].strided!(0, 1)(2, 2)[] += U[1 
.. m - 1, 0 .. n - 2].strided!(0, 1)(2, 2);
         U[1 .. m - 1, 1 .. n - 1].strided!(0, 1)(2, 2)[] += U[1 
.. m - 1, 2 .. n].strided!(0, 1)(2, 2);
         U[1 .. m - 1, 1 .. n - 1].strided!(0, 1)(2, 2)[] -= F[1 
.. m - 1, 1 .. n - 1].strided!(0, 1)(2, 2) * h2;
         U[1 .. m - 1, 1 .. n - 1].strided!(0, 1)(2, 2)[] /= 
cast(T) 4.0;

         U[2 .. m - 1, 2 .. n - 1].strided!(0, 1)(2, 2)[] = U[1 .. 
m - 2, 2 .. n - 1].strided!(0, 1)(2, 2);
         U[2 .. m - 1, 2 .. n - 1].strided!(0, 1)(2, 2)[] += U[3 
.. m, 2 .. n - 1].strided!(0, 1)(2, 2);
         U[2 .. m - 1, 2 .. n - 1].strided!(0, 1)(2, 2)[] += U[2 
.. m - 1, 1 .. n - 2].strided!(0, 1)(2, 2);
         U[2 .. m - 1, 2 .. n - 1].strided!(0, 1)(2, 2)[] += U[2 
.. m - 1, 3 .. n].strided!(0, 1)(2, 2);
         U[2 .. m - 1, 2 .. n - 1].strided!(0, 1)(2, 2)[] -= F[2 
.. m - 1, 2 .. n - 1].strided!(0, 1)(2, 2) * h2;
         U[2 .. m - 1, 2 .. n - 1].strided!(0, 1)(2, 2)[] /= 
cast(T) 4.0;
     }
     else static if (color == Color.black)
     {
         U[1 .. m - 1, 2 .. n - 1].strided!(0, 1)(2, 2)[] = U[0 .. 
m - 2, 2 .. n - 1].strided!(0, 1)(2, 2) / 4.0;
         U[1 .. m - 1, 2 .. n - 1].strided!(0, 1)(2, 2)[] += U[2 
.. m, 2 .. n - 1].strided!(0, 1)(2, 2) / 4.0;
         U[1 .. m - 1, 2 .. n - 1].strided!(0, 1)(2, 2)[] += U[1 
.. m - 1, 1 .. n - 2].strided!(0, 1)(2, 2) / 4.0;
         U[1 .. m - 1, 2 .. n - 1].strided!(0, 1)(2, 2)[] += U[1 
.. m - 1, 3 .. n].strided!(0, 1)(2, 2) / 4.0;
         U[1 .. m - 1, 2 .. n - 1].strided!(0, 1)(2, 2)[] -= F[1 
.. m - 1, 2 .. n - 1].strided!(0, 1)(2, 2) * h2 / 4.0;

         U[2 .. m - 1, 1 .. n - 1].strided!(0, 1)(2, 2)[] = U[1 .. 
m - 2, 1 .. n - 1].strided!(0, 1)(2, 2) / 4.0;
         U[2 .. m - 1, 1 .. n - 1].strided!(0, 1)(2, 2)[] += U[3 
.. m, 1 .. n - 1].strided!(0, 1)(2, 2) / 4.0;
         U[2 .. m - 1, 1 .. n - 1].strided!(0, 1)(2, 2)[] += U[2 
.. m - 1, 0 .. n - 2].strided!(0, 1)(2, 2) / 4.0;
         U[2 .. m - 1, 1 .. n - 1].strided!(0, 1)(2, 2)[] += U[2 
.. m - 1, 2 .. n].strided!(0, 1)(2, 2) / 4.0;
         U[2 .. m - 1, 1 .. n - 1].strided!(0, 1)(2, 2)[] -= F[2 
.. m - 1, 1 .. n - 1].strided!(0, 1)(2, 2) * h2 / 4.0;
     }
     else
     {
         static assert(false, color.stringof ~ "invalid color");
     }
}
```
Thank you very much for your time.

Christoph


More information about the Digitalmars-d-learn mailing list