src/gromacs/ewald/pme_spread.clh

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35
  36 /*! \internal \file
  37  *  \brief Implements PME OpenCL spline parameter computation and charge spread kernels.
  38  * When including this and other PME OpenCL kernel files, plenty of common
  39  * constants/macros are expected to be defined (such as "order" which is PME interpolation order).
  40  * For details, please see how pme_program.cl is compiled in pme_gpu_program_impl_ocl.cpp.
  41  *
  42  * This file's kernels specifically expect the following definitions:
  43  *
  44  * - atomsPerBlock which expresses how many atoms are processed by a single work group
  45  * - order which is a PME interpolation order
  46  * - computeSplines and spreadCharges must evaluate to either true or false to specify which
  47  * kernel falvor is being compiled
  48  * - wrapX and wrapY must evaluate to either true or false to specify whether the grid overlap
  49  * in dimension X/Y is to be used
  50  *
  51  *  \author Aleksei Iupinov <a.yupinov@gmail.com>
  52  */
  53
  54 #include "gromacs/gpu_utils/vectype_ops.clh"
  55
  56 #include "pme_gpu_types.h"
  57 #include "pme_gpu_calculate_splines.clh"
  58
  59 /*
  60  * This define affects the spline calculation behaviour in the kernel.
  61  * 0: a single GPU thread handles a single dimension of a single particle (calculating and storing
  62  * (order) spline values and derivatives). 1: (order) threads do redundant work on this same task,
  63  * each one stores only a single theta and single dtheta into global arrays. The only efficiency
  64  * difference is less global store operations, countered by more redundant spline computation.
  65  *
  66  * TODO: estimate if this should be a boolean parameter (and add it to the unit test if so).
  67  */
  68 #define PME_GPU_PARALLEL_SPLINE 0
  69
  70 #ifndef COMPILE_SPREAD_HELPERS_ONCE
  71 #    define COMPILE_SPREAD_HELPERS_ONCE
  72
  73 /*! \brief
  74  * General purpose function for loading atom-related data from global to shared memory.
  75  *
  76  * \param[out] sm_destination    Local memory array for output.
  77  * \param[in]  gm_source         Global memory array for input.
  78  * \param[in] dataCountPerAtom   Number of data elements per single atom (e.g. DIM for an rvec coordinates array).
  79  *
  80  */
  81 inline void pme_gpu_stage_atom_data(__local float* __restrict__ sm_destination,
  82                                     __global const float* __restrict__ gm_source,
  83                                     const int dataCountPerAtom)
  84 {
  85     assert((get_local_id(2) * get_local_size(1) + get_local_id(1)) * get_local_size(0) + get_local_id(0)
  86            < MAX_INT);
  87     const int threadLocalIndex =
  88             (int)((get_local_id(2) * get_local_size(1) + get_local_id(1)) * get_local_size(0)
  89                   + get_local_id(0));
  90     const int localIndex      = threadLocalIndex;
  91     const int globalIndexBase = (int)get_group_id(XX) * atomsPerBlock * dataCountPerAtom;
  92     const int globalIndex     = globalIndexBase + localIndex;
  93     if (localIndex < atomsPerBlock * dataCountPerAtom)
  94     {
  95         assert(isfinite(float(gm_source[globalIndex])));
  96         sm_destination[localIndex] = gm_source[globalIndex];
  97     }
  98 }
  99
 100 /*! \brief
 101  * PME GPU spline parameter and gridline indices calculation.
 102  * This corresponds to the CPU functions calc_interpolation_idx() and make_bsplines().
 103  * First stage of the whole kernel.
 104  *
 105  * \param[in]  kernelParams             Input PME GPU data in constant memory.
 106  * \param[in]  atomIndexOffset          Starting atom index for the execution block w.r.t. global memory.
 107  * \param[in]  sm_coordinates           Atom coordinates in the shared memory (rvec)
 108  * \param[in]  sm_coefficients          Atom charges/coefficients in the shared memory.
 109  * \param[out] sm_theta                 Atom spline values in the shared memory.
 110  * \param[out] sm_gridlineIndices       Atom gridline indices in the shared memory.
 111  * \param[out] sm_fractCoords           Atom fractional coordinates (rvec)
 112  * \param[out] gm_theta                 Atom spline parameter values
 113  * \param[out] gm_dtheta                Atom spline parameter derivatives
 114  * \param[out] gm_gridlineIndices       Same as \p sm_gridlineIndices but in global memory.
 115  * \param[in]  gm_fractShiftsTable      Atom fractional coordinates correction table
 116  * \param[in]  gm_gridlineIndicesTable  Atom fractional coordinates correction table
 117  */
 118 gmx_opencl_inline void calculate_splines(const struct PmeOpenCLKernelParams kernelParams,
 119                                          const int                          atomIndexOffset,
 120                                          __local const float* __restrict__ sm_coordinates,
 121                                          __local const float* __restrict__ sm_coefficients,
 122                                          __local float* __restrict__ sm_theta,
 123                                          __local int* __restrict__ sm_gridlineIndices,
 124                                          __local float* __restrict__ sm_fractCoords,
 125                                          __global float* __restrict__ gm_theta,
 126                                          __global float* __restrict__ gm_dtheta,
 127                                          __global int* __restrict__ gm_gridlineIndices,
 128                                          __global const float* __restrict__ gm_fractShiftsTable,
 129                                          __global const int* __restrict__ gm_gridlineIndicesTable)
 130 {
 131     /* Thread index w.r.t. block */
 132     assert((get_local_id(2) * get_local_size(1) + get_local_id(1)) * get_local_size(0) + get_local_id(0)
 133            < MAX_INT);
 134     const int threadLocalIndex =
 135             (int)((get_local_id(2) * get_local_size(1) + get_local_id(1)) * get_local_size(0)
 136                   + get_local_id(0));
 137     /* Warp index w.r.t. block - could probably be obtained easier? */
 138     const int warpIndex = threadLocalIndex / warp_size;
 139     /* Thread index w.r.t. warp */
 140     const int threadWarpIndex = threadLocalIndex % warp_size;
 141     /* Atom index w.r.t. warp - alternating 0 1 0 1 .. */
 142     const int atomWarpIndex = threadWarpIndex % atomsPerWarp;
 143     /* Atom index w.r.t. block/shared memory */
 144     const int atomIndexLocal = warpIndex * atomsPerWarp + atomWarpIndex;
 145
 146     /* Spline contribution index in one dimension */
 147     const int orderIndex = threadWarpIndex / (atomsPerWarp * DIM);
 148     /* Dimension index */
 149     const int dimIndex = (threadWarpIndex - orderIndex * (atomsPerWarp * DIM)) / atomsPerWarp;
 150
 151     /* Multi-purpose index of rvec/ivec atom data */
 152     const int sharedMemoryIndex = atomIndexLocal * DIM + dimIndex;
 153
 154     /* Spline parameters need a working buffer.
 155      * With PME_GPU_PARALLEL_SPLINE == 0 it is just a local array of (order) values for some of the threads, which is fine;
 156      * With PME_GPU_PARALLEL_SPLINE == 1 (order) times more threads are involved, so the shared memory is used to avoid overhead.
 157      * The buffer's size, striding and indexing are adjusted accordingly.
 158      * The buffer is accessed with SPLINE_DATA_PTR and SPLINE_DATA macros.
 159      */
 160 #    if PME_GPU_PARALLEL_SPLINE
 161 #        define splineDataStride (atomsPerBlock * DIM)
 162     const int      splineDataIndex = sharedMemoryIndex;
 163     __local float  sm_splineData[splineDataStride * order];
 164     __local float* splineDataPtr = sm_splineData;
 165 #    else
 166 #        define splineDataStride 1
 167     const int splineDataIndex = 0;
 168     float     splineData[splineDataStride * order];
 169     float*    splineDataPtr = splineData;
 170 #    endif
 171
 172 #    define SPLINE_DATA_PTR(i) (splineDataPtr + ((i)*splineDataStride + splineDataIndex))
 173 #    define SPLINE_DATA(i) (*SPLINE_DATA_PTR(i))
 174
 175     const int localCheck = (dimIndex < DIM) && (orderIndex < (PME_GPU_PARALLEL_SPLINE ? order : 1));
 176
 177     if (localCheck)
 178     {
 179         /* Indices interpolation */
 180         if (orderIndex == 0)
 181         {
 182             int          tableIndex;
 183             int          tInt;
 184             float        n;
 185             float        t;
 186             const float3 x = vload3(atomIndexLocal, sm_coordinates);
 187
 188             /* Accessing fields in fshOffset/nXYZ/recipbox/... with dimIndex offset
 189              * puts them into local memory(!) instead of accessing the constant memory directly.
 190              * That's the reason for the switch, to unroll explicitly.
 191              * The commented parts correspond to the 0 components of the recipbox.
 192              */
 193             switch (dimIndex)
 194             {
 195                 case XX:
 196                     tableIndex = kernelParams.grid.tablesOffsets[XX];
 197                     n          = kernelParams.grid.realGridSizeFP[XX];
 198                     t          = x.x * kernelParams.current.recipBox[dimIndex][XX]
 199                         + x.y * kernelParams.current.recipBox[dimIndex][YY]
 200                         + x.z * kernelParams.current.recipBox[dimIndex][ZZ];
 201                     break;
 202
 203                 case YY:
 204                     tableIndex = kernelParams.grid.tablesOffsets[YY];
 205                     n          = kernelParams.grid.realGridSizeFP[YY];
 206                     t          = /*x.x * kernelParams.current.recipbox[dimIndex][XX] + */ x.y
 207                                 * kernelParams.current.recipBox[dimIndex][YY]
 208                         + x.z * kernelParams.current.recipBox[dimIndex][ZZ];
 209                     break;
 210
 211                 case ZZ:
 212                     tableIndex = kernelParams.grid.tablesOffsets[ZZ];
 213                     n          = kernelParams.grid.realGridSizeFP[ZZ];
 214                     t          = /*x.x * kernelParams.current.recipbox[dimIndex][XX] + x.y * kernelParams.current.recipbox[dimIndex][YY] + */ x
 215                                 .z
 216                         * kernelParams.current.recipBox[dimIndex][ZZ];
 217                     break;
 218                 default:
 219                     assert(false);
 220                     return;
 221                     break;
 222             }
 223             const float shift = c_pmeMaxUnitcellShift;
 224
 225             /* Fractional coordinates along box vectors, adding a positive shift to ensure t is positive for triclinic boxes */
 226             t                                 = (t + shift) * n;
 227             tInt                              = (int)t;
 228             sm_fractCoords[sharedMemoryIndex] = t - (float)tInt;
 229             tableIndex += tInt;
 230             assert(tInt >= 0);
 231             assert(tInt < c_pmeNeighborUnitcellCount * n);
 232             sm_fractCoords[sharedMemoryIndex] += gm_fractShiftsTable[tableIndex];
 233             sm_gridlineIndices[sharedMemoryIndex] = gm_gridlineIndicesTable[tableIndex];
 234             gm_gridlineIndices[atomIndexOffset * DIM + sharedMemoryIndex] =
 235                     sm_gridlineIndices[sharedMemoryIndex];
 236         }
 237
 238         /* B-spline calculation */
 239
 240         const int chargeCheck = pme_gpu_check_atom_charge(sm_coefficients[atomIndexLocal]);
 241         if (chargeCheck)
 242         {
 243             float div;
 244             int o = orderIndex; // This is an index that is set once for PME_GPU_PARALLEL_SPLINE == 1
 245
 246             const float dr = sm_fractCoords[sharedMemoryIndex];
 247             assert(isfinite(dr));
 248
 249             /* dr is relative offset from lower cell limit */
 250             *SPLINE_DATA_PTR(order - 1) = 0.0F;
 251             *SPLINE_DATA_PTR(1)         = dr;
 252             *SPLINE_DATA_PTR(0)         = 1.0F - dr;
 253
 254 #    pragma unroll order
 255             for (int k = 3; k < order; k++)
 256             {
 257                 div                     = 1.0F / ((float)k - 1.0F);
 258                 *SPLINE_DATA_PTR(k - 1) = div * dr * SPLINE_DATA(k - 2);
 259 #    pragma unroll
 260                 for (int l = 1; l < (k - 1); l++)
 261                 {
 262                     *SPLINE_DATA_PTR(k - l - 1) =
 263                             div
 264                             * ((dr + (float)l) * SPLINE_DATA(k - l - 2)
 265                                + ((float)k - (float)l - dr) * SPLINE_DATA(k - l - 1));
 266                 }
 267                 *SPLINE_DATA_PTR(0) = div * (1.0F - dr) * SPLINE_DATA(0);
 268             }
 269
 270             const int thetaIndexBase        = getSplineParamIndexBase(warpIndex, atomWarpIndex);
 271             const int thetaGlobalOffsetBase = atomIndexOffset * DIM * order;
 272
 273             /* Differentiation and storing the spline derivatives (dtheta) */
 274 #    if !PME_GPU_PARALLEL_SPLINE
 275             // With PME_GPU_PARALLEL_SPLINE == 1, o is already set to orderIndex;
 276             // With PME_GPU_PARALLEL_SPLINE == 0, we loop o over range(order).
 277 #        pragma unroll order
 278             for (o = 0; o < order; o++)
 279 #    endif
 280             {
 281                 const int thetaIndex       = getSplineParamIndex(thetaIndexBase, dimIndex, o);
 282                 const int thetaGlobalIndex = thetaGlobalOffsetBase + thetaIndex;
 283
 284                 const float dtheta = ((o > 0) ? SPLINE_DATA(o - 1) : 0.0F) - SPLINE_DATA(o);
 285                 assert(isfinite(dtheta));
 286                 gm_dtheta[thetaGlobalIndex] = dtheta;
 287             }
 288
 289             div                         = 1.0F / (order - 1.0F);
 290             *SPLINE_DATA_PTR(order - 1) = div * dr * SPLINE_DATA(order - 2);
 291 #    pragma unroll
 292             for (int k = 1; k < (order - 1); k++)
 293             {
 294                 *SPLINE_DATA_PTR(order - k - 1) = div
 295                                                   * ((dr + (float)k) * SPLINE_DATA(order - k - 2)
 296                                                      + (order - k - dr) * SPLINE_DATA(order - k - 1));
 297             }
 298             *SPLINE_DATA_PTR(0) = div * (1.0F - dr) * SPLINE_DATA(0);
 299
 300             /* Storing the spline values (theta) */
 301 #    if !PME_GPU_PARALLEL_SPLINE
 302             // See comment for the loop above
 303 #        pragma unroll order
 304             for (o = 0; o < order; o++)
 305 #    endif
 306             {
 307                 const int thetaIndex       = getSplineParamIndex(thetaIndexBase, dimIndex, o);
 308                 const int thetaGlobalIndex = thetaGlobalOffsetBase + thetaIndex;
 309
 310                 sm_theta[thetaIndex] = SPLINE_DATA(o);
 311                 assert(isfinite(sm_theta[thetaIndex]));
 312                 gm_theta[thetaGlobalIndex] = SPLINE_DATA(o);
 313             }
 314         }
 315     }
 316 }
 317
 318 /*! \brief
 319  * Charge spreading onto the grid.
 320  * This corresponds to the CPU function spread_coefficients_bsplines_thread().
 321  * Second stage of the whole kernel.
 322  *
 323  * \param[in]  kernelParams         Input PME GPU data in constant memory.
 324  * \param[in]  sm_coefficients      Atom coefficents/charges in the shared memory.
 325  * \param[in]  sm_gridlineIndices   Atom gridline indices in the shared memory.
 326  * \param[in]  sm_theta             Atom spline values in the shared memory.
 327  * \param[out] gm_grid              Global 3D grid for spreading.
 328  */
 329 gmx_opencl_inline void spread_charges(const struct PmeOpenCLKernelParams kernelParams,
 330                                       __local const float* __restrict__ sm_coefficients,
 331                                       __local const int* __restrict__ sm_gridlineIndices,
 332                                       __local const float* __restrict__ sm_theta,
 333                                       __global float* __restrict__ gm_grid)
 334 {
 335     const int nx  = kernelParams.grid.realGridSize[XX];
 336     const int ny  = kernelParams.grid.realGridSize[YY];
 337     const int nz  = kernelParams.grid.realGridSize[ZZ];
 338     const int pny = kernelParams.grid.realGridSizePadded[YY];
 339     const int pnz = kernelParams.grid.realGridSizePadded[ZZ];
 340
 341     const int offx = 0;
 342     const int offy = 0;
 343     const int offz = 0;
 344
 345     const int atomIndexLocal = get_local_id(ZZ);
 346
 347     const int chargeCheck = pme_gpu_check_atom_charge(sm_coefficients[atomIndexLocal]);
 348     if (chargeCheck)
 349     {
 350         // Spline Y/Z coordinates
 351         const int ithy   = get_local_id(YY);
 352         const int ithz   = get_local_id(XX);
 353         const int ixBase = sm_gridlineIndices[atomIndexLocal * DIM + XX] - offx;
 354         int       iy     = sm_gridlineIndices[atomIndexLocal * DIM + YY] - offy + ithy;
 355         if (wrapY & (iy >= ny))
 356         {
 357             iy -= ny;
 358         }
 359         int iz = sm_gridlineIndices[atomIndexLocal * DIM + ZZ] - offz + ithz;
 360         if (iz >= nz)
 361         {
 362             iz -= nz;
 363         }
 364
 365         /* Atom index w.r.t. warp - alternating 0 1 0 1 .. */
 366         const int atomWarpIndex = atomIndexLocal % atomsPerWarp;
 367         /* Warp index w.r.t. block - could probably be obtained easier? */
 368         const int warpIndex = atomIndexLocal / atomsPerWarp;
 369
 370         const int   splineIndexBase = getSplineParamIndexBase(warpIndex, atomWarpIndex);
 371         const int   splineIndexZ    = getSplineParamIndex(splineIndexBase, ZZ, ithz);
 372         const float thetaZ          = sm_theta[splineIndexZ];
 373         const int   splineIndexY    = getSplineParamIndex(splineIndexBase, YY, ithy);
 374         const float thetaY          = sm_theta[splineIndexY];
 375         const float constVal        = thetaZ * thetaY * sm_coefficients[atomIndexLocal];
 376         assert(isfinite(constVal));
 377         const int constOffset = iy * pnz + iz;
 378
 379 #    pragma unroll order
 380         for (int ithx = 0; (ithx < order); ithx++)
 381         {
 382             int ix = ixBase + ithx;
 383             if (wrapX & (ix >= nx))
 384             {
 385                 ix -= nx;
 386             }
 387             const int   gridIndexGlobal = ix * pny * pnz + constOffset;
 388             const int   splineIndexX    = getSplineParamIndex(splineIndexBase, XX, ithx);
 389             const float thetaX          = sm_theta[splineIndexX];
 390             assert(isfinite(thetaX));
 391             assert(isfinite(gm_grid[gridIndexGlobal]));
 392             atomicAdd_g_f(gm_grid + gridIndexGlobal, thetaX * constVal);
 393         }
 394     }
 395 }
 396
 397 #endif // COMPILE_SPREAD_HELPERS_ONCE
 398
 399 /*! \brief
 400  * A spline computation and/or charge spreading kernel function.
 401  * Please see the file description for additional defines which this kernel expects.
 402  *
 403  * \param[in]     kernelParams             Input PME GPU data in constant memory.
 404  * \param[in,out] gm_theta                 Atom spline parameter values
 405  * \param[out]    gm_dtheta                Atom spline parameter derivatives
 406  * \param[in,out] gm_gridlineIndices       Atom gridline indices (ivec)
 407  * \param[out]    gm_grid                  Global 3D grid for charge spreading.
 408  * \param[in]     gm_fractShiftsTable      Atom fractional coordinates correction table
 409  * \param[in]     gm_gridlineIndicesTable  Atom fractional coordinates correction table
 410  * \param[in]     gm_coefficients          Atom charges/coefficients.
 411  * \param[in]     gm_coordinates           Atom coordinates (rvec)
 412  */
 413 __attribute__((reqd_work_group_size(order, order, atomsPerBlock))) __kernel void CUSTOMIZED_KERNEL_NAME(
 414         pme_spline_and_spread_kernel)(const struct PmeOpenCLKernelParams kernelParams,
 415                                       __global float* __restrict__ gm_theta,
 416                                       __global float* __restrict__ gm_dtheta,
 417                                       __global int* __restrict__ gm_gridlineIndices,
 418                                       __global float* __restrict__ gm_grid,
 419                                       __global const float* __restrict__ gm_fractShiftsTable,
 420                                       __global const int* __restrict__ gm_gridlineIndicesTable,
 421                                       __global const float* __restrict__ gm_coefficients,
 422                                       __global const float* __restrict__ gm_coordinates)
 423 {
 424     // Gridline indices, ivec
 425     __local int sm_gridlineIndices[atomsPerBlock * DIM];
 426     // Charges
 427     __local float sm_coefficients[atomsPerBlock];
 428     // Spline values
 429     __local float sm_theta[atomsPerBlock * DIM * order];
 430     // Fractional coordinates - only for spline computation
 431     __local float sm_fractCoords[atomsPerBlock * DIM];
 432     // Staging coordinates - only for spline computation
 433     __local float sm_coordinates[DIM * atomsPerBlock];
 434
 435     const int atomIndexOffset = (int)get_group_id(XX) * atomsPerBlock;
 436
 437     /* Staging coefficients/charges for both spline and spread */
 438     pme_gpu_stage_atom_data(sm_coefficients, gm_coefficients, 1);
 439
 440     if (computeSplines)
 441     {
 442         /* Staging coordinates */
 443         pme_gpu_stage_atom_data(sm_coordinates, gm_coordinates, DIM);
 444
 445         barrier(CLK_LOCAL_MEM_FENCE);
 446         calculate_splines(kernelParams, atomIndexOffset, sm_coordinates, sm_coefficients, sm_theta,
 447                           sm_gridlineIndices, sm_fractCoords, gm_theta, gm_dtheta,
 448                           gm_gridlineIndices, gm_fractShiftsTable, gm_gridlineIndicesTable);
 449 #if !defined(_AMD_SOURCE_) && !defined(_NVIDIA_SOURCE_)
 450         /* This is only here for execution of e.g. 32-sized warps on 16-wide hardware; this was
 451          * __syncwarp() in CUDA. #2519
 452          */
 453         barrier(CLK_LOCAL_MEM_FENCE);
 454 #endif
 455     }
 456     else
 457     {
 458         /* Staging the data for spread
 459          * (the data is assumed to be in GPU global memory with proper layout already,
 460          * as in after running the spline kernel)
 461          */
 462         /* Spline data - only thetas (dthetas will only be needed in gather) */
 463         pme_gpu_stage_atom_data(sm_theta, gm_theta, DIM * order);
 464         /* Gridline indices - they're actually int and not float, but C99 is angry about overloads */
 465         pme_gpu_stage_atom_data((__local float*)sm_gridlineIndices,
 466                                 (__global const float*)gm_gridlineIndices, DIM);
 467
 468         barrier(CLK_LOCAL_MEM_FENCE);
 469     }
 470     /* Spreading */
 471     if (spreadCharges)
 472     {
 473         spread_charges(kernelParams, sm_coefficients, sm_gridlineIndices, sm_theta, gm_grid);
 474     }
 475 }