Minor optimization by exploiting demag kernel symmetry

cuda/conv_kernmul.go

@@ -19,20 +19,24 @@ func kernMulRSymm3D_async(fftM [3]*data.Slice, Kxx, Kyy, Kzz, Kyz, Kxz, Kxy *dat
 }
 
 // kernel multiplication for 2D demag convolution on X and Y, exploiting full kernel symmetry.
+// Launch grid covers only the non-redundant rows iy in [0, Ny/2]; each thread
+// also writes the mirrored row Ny-iy when distinct (see kernmulrsymm2dxy.cu).
 func kernMulRSymm2Dxy_async(fftMx, fftMy, Kxx, Kyy, Kxy *data.Slice, Nx, Ny int) {
 	util.Argument(fftMy.NComp() == 1 && Kxx.NComp() == 1)
 
-	cfg := make3DConf([3]int{Nx, Ny, 1})
+	cfg := make3DConf([3]int{Nx, Ny/2 + 1, 1})
 	k_kernmulRSymm2Dxy_async(fftMx.DevPtr(0), fftMy.DevPtr(0),
 		Kxx.DevPtr(0), Kyy.DevPtr(0), Kxy.DevPtr(0),
 		Nx, Ny, cfg)
 }
 
 // kernel multiplication for 2D demag convolution on Z, exploiting full kernel symmetry.
+// Launch grid covers only the non-redundant rows iy in [0, Ny/2]; each thread
+// also writes the mirrored row Ny-iy when distinct (see kernmulrsymm2dz.cu).
 func kernMulRSymm2Dz_async(fftMz, Kzz *data.Slice, Nx, Ny int) {
 	util.Argument(fftMz.NComp() == 1 && Kzz.NComp() == 1)
 
-	cfg := make3DConf([3]int{Nx, Ny, 1})
+	cfg := make3DConf([3]int{Nx, Ny/2 + 1, 1})
 	k_kernmulRSymm2Dz_async(fftMz.DevPtr(0), Kzz.DevPtr(0), Nx, Ny, cfg)
 }
 

cuda/kernmulrsymm2dxy.cu

@@ -1,7 +1,12 @@
 // 2D XY (in-plane) micromagnetic kernel multiplication:
 // |Mx| = |Kxx Kxy| * |Mx|
 // |My|   |Kyx Kyy|   |My|
-// Using the same symmetries as kernmulrsymm3d.cu
+//
+// fftKxx/fftKyy/fftKxy only store the non-redundant rows iy in [0, Ny/2].
+// Launch grid covers only that half range; each thread handles row iy
+// (using Kxy as-is) and, if it has a distinct mirror row Ny-iy, also
+// handles that row (using -Kxy), avoiding launching threads for rows
+// that would just re-read the same kernel value.
 extern "C" __global__ void
 kernmulRSymm2Dxy(float* __restrict__  fftMx,  float* __restrict__  fftMy,
                  float* __restrict__  fftKxx, float* __restrict__  fftKyy, float* __restrict__  fftKxy,
@@ -10,33 +15,43 @@ kernmulRSymm2Dxy(float* __restrict__  fftMx,  float* __restrict__  fftMy,
     int ix = blockIdx.x * blockDim.x + threadIdx.x;
     int iy = blockIdx.y * blockDim.y + threadIdx.y;
 
-    if(ix>= Nx || iy>=Ny) {
+    if(ix>=Nx || iy>Ny/2) {
         return;
     }
 
     int I = iy*Nx + ix;
-    int e = 2 * I;
-
-    float reMx = fftMx[e  ];
-    float imMx = fftMx[e+1];
-    float reMy = fftMy[e  ];
-    float imMy = fftMy[e+1];
-
-    // symmetry factor
-    float fxy = 1.0f;
-    if (iy > Ny/2) {
-        iy = Ny-iy;
-        fxy = -fxy;
-    }
-    I = iy*Nx + ix;
-
     float Kxx = fftKxx[I];
     float Kyy = fftKyy[I];
-    float Kxy = fxy * fftKxy[I];
+    float Kxy = fftKxy[I];
+
+    // row iy
+    {
+        int e = 2 * (iy*Nx + ix);
+        float reMx = fftMx[e  ];
+        float imMx = fftMx[e+1];
+        float reMy = fftMy[e  ];
+        float imMy = fftMy[e+1];
+
+        fftMx[e  ] = reMx * Kxx + reMy * Kxy;
+        fftMx[e+1] = imMx * Kxx + imMy * Kxy;
+        fftMy[e  ] = reMx * Kxy + reMy * Kyy;
+        fftMy[e+1] = imMx * Kxy + imMy * Kyy;
+    }
 
-    fftMx[e  ] = reMx * Kxx + reMy * Kxy;
-    fftMx[e+1] = imMx * Kxx + imMy * Kxy;
-    fftMy[e  ] = reMx * Kxy + reMy * Kyy;
-    fftMy[e+1] = imMx * Kxy + imMy * Kyy;
+    // mirror row Ny-iy, if distinct from iy (sign of Kxy flips)
+    if (iy != 0 && 2*iy != Ny) {
+        int my = Ny - iy;
+        int e = 2 * (my*Nx + ix);
+        float reMx = fftMx[e  ];
+        float imMx = fftMx[e+1];
+        float reMy = fftMy[e  ];
+        float imMy = fftMy[e+1];
+        float Kxym = -Kxy;
+
+        fftMx[e  ] = reMx * Kxx + reMy * Kxym;
+        fftMx[e+1] = imMx * Kxx + imMy * Kxym;
+        fftMy[e  ] = reMx * Kxym + reMy * Kyy;
+        fftMy[e+1] = imMx * Kxym + imMy * Kyy;
+    }
 }
 

cuda/kernmulrsymm2dz.cu

@@ -1,30 +1,41 @@
 // 2D Z (out-of-plane only) micromagnetic kernel multiplication:
 // Mz = Kzz * Mz
-// Using the same symmetries as kernmulrsymm3d.cu
+//
+// fftKzz only stores the non-redundant rows iy in [0, Ny/2].
+// Launch grid covers only that half range; each thread handles row iy
+// and, if it has a distinct mirror row Ny-iy (same Kzz value), handles
+// that row too, avoiding launching threads for rows that would just
+// re-read the same kernel value.
 extern "C" __global__ void
 kernmulRSymm2Dz(float* __restrict__  fftMz, float* __restrict__  fftKzz, int Nx, int Ny) {
 
     int ix = blockIdx.x * blockDim.x + threadIdx.x;
     int iy = blockIdx.y * blockDim.y + threadIdx.y;
 
-    if(ix>= Nx || iy>=Ny) {
+    if(ix>=Nx || iy>Ny/2) {
         return;
     }
 
     int I = iy*Nx + ix;
-    int e = 2 * I;
-
-    float reMz = fftMz[e  ];
-    float imMz = fftMz[e+1];
+    float Kzz = fftKzz[I];
 
-    if (iy > Ny/2) {
-        iy = Ny-iy;
+    // row iy
+    {
+        int e = 2 * (iy*Nx + ix);
+        float reMz = fftMz[e  ];
+        float imMz = fftMz[e+1];
+        fftMz[e  ] = reMz * Kzz;
+        fftMz[e+1] = imMz * Kzz;
     }
-    I = iy*Nx + ix;
 
-    float Kzz = fftKzz[I];
-
-    fftMz[e  ] = reMz * Kzz;
-    fftMz[e+1] = imMz * Kzz;
+    // mirror row Ny-iy, if distinct from iy
+    if (iy != 0 && 2*iy != Ny) {
+        int my = Ny - iy;
+        int e = 2 * (my*Nx + ix);
+        float reMz = fftMz[e  ];
+        float imMz = fftMz[e+1];
+        fftMz[e  ] = reMz * Kzz;
+        fftMz[e+1] = imMz * Kzz;
+    }
 }