kCalculateAmoebaCudaPmeFixedEField.h 18.4 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
/* -------------------------------------------------------------------------- *
 *                                   OpenMM                                   *
 * -------------------------------------------------------------------------- *
 * This is part of the OpenMM molecular simulation toolkit originating from   *
 * Simbios, the NIH National Center for Physics-Based Simulation of           *
 * Biological Structures at Stanford, funded under the NIH Roadmap for        *
 * Medical Research, grant U54 GM072970. See https://simtk.org.               *
 *                                                                            *
 * Portions copyright (c) 2009 Stanford University and the Authors.           *
 * Authors: Scott Le Grand, Peter Eastman                                     *
 * Contributors:                                                              *
 *                                                                            *
 * This program is free software: you can redistribute it and/or modify       *
 * it under the terms of the GNU Lesser General Public License as published   *
 * by the Free Software Foundation, either version 3 of the License, or       *
 * (at your option) any later version.                                        *
 *                                                                            *
 * This program is distributed in the hope that it will be useful,            *
 * but WITHOUT ANY WARRANTY; without even the implied warranty of             *
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              *
 * GNU Lesser General Public License for more details.                        *
 *                                                                            *
 * You should have received a copy of the GNU Lesser General Public License   *
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.      *
 * -------------------------------------------------------------------------- */

#include "amoebaScaleFactors.h"

__global__ 
#if (__CUDA_ARCH__ >= 200)
__launch_bounds__(GF1XX_NONBOND_THREADS_PER_BLOCK, 1)
#elif (__CUDA_ARCH__ >= 130)
__launch_bounds__(GT2XX_NONBOND_THREADS_PER_BLOCK, 1)
#else
__launch_bounds__(G8X_NONBOND_THREADS_PER_BLOCK, 1)
#endif
void METHOD_NAME(kCalculateAmoebaPmeDirectFixedE_Field, _kernel)(
                            unsigned int* workUnit,
                            float* outputEField,
                            float* outputEFieldPolar
#ifdef AMOEBA_DEBUG
                           , float4* debugArray, unsigned int targetAtom
#endif
){

#ifdef AMOEBA_DEBUG
Mark Friedrichs's avatar
Mark Friedrichs committed
47
    int maxPullIndex = 1;
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
    float4 pullBack[12];
#endif

    extern __shared__ FixedFieldParticle sA[];

    unsigned int totalWarps      = gridDim.x*blockDim.x/GRID;
    unsigned int warp            = (blockIdx.x*blockDim.x+threadIdx.x)/GRID;
    unsigned int numWorkUnits    = cSim.pInteractionCount[0];
    unsigned int pos             = warp*numWorkUnits/totalWarps;
    unsigned int end             = (warp+1)*numWorkUnits/totalWarps;
    unsigned int lasty           = 0xFFFFFFFF;

    while (pos < end)
    {

        unsigned int x;
        unsigned int y;
        bool bExclusionFlag;
Mark Friedrichs's avatar
Mark Friedrichs committed
66
67
68
69
        float dScaleValue;
        float pScaleValue;
        int  dScaleMask;
        int2 pScaleMask;
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101

        // extract cell coordinates

        decodeCell( workUnit[pos], &x, &y, &bExclusionFlag );

        unsigned int tgx           = threadIdx.x & (GRID - 1);
        unsigned int tbx           = threadIdx.x - tgx;
        unsigned int tj            = tgx;

        FixedFieldParticle* psA    = &sA[tbx];
        unsigned int atomI         = x + tgx;
        FixedFieldParticle localParticle;
        loadFixedFieldShared( &localParticle, atomI );

        float fieldSum[3];
        float fieldPolarSum[3];

        fieldSum[0]                = 0.0f;
        fieldSum[1]                = 0.0f;
        fieldSum[2]                = 0.0f;

        fieldPolarSum[0]           = 0.0f;
        fieldPolarSum[1]           = 0.0f;
        fieldPolarSum[2]           = 0.0f;

        if (x == y)
        {

            // load coordinates, charge, ...

            loadFixedFieldShared( &(sA[threadIdx.x]), atomI );

Mark Friedrichs's avatar
Mark Friedrichs committed
102
            if( bExclusionFlag ){
103
104
                unsigned int xi       = x >> GRIDBITS;
                unsigned int cell     = xi + xi*cAmoebaSim.paddedNumberOfAtoms/GRID-xi*(xi+1)/2;
Mark Friedrichs's avatar
Mark Friedrichs committed
105
106
107
108
                dScaleMask            = cAmoebaSim.pD_ScaleIndices[cAmoebaSim.pScaleIndicesIndex[cell]+tgx];
                pScaleMask            = cAmoebaSim.pP_ScaleIndices[cAmoebaSim.pScaleIndicesIndex[cell]+tgx];
            } else {
                dScaleValue = pScaleValue = 1.0f;
109

Mark Friedrichs's avatar
Mark Friedrichs committed
110
            }
111

Mark Friedrichs's avatar
Mark Friedrichs committed
112
113
            for (unsigned int j = 0; j < GRID; j++)
            {
114

Mark Friedrichs's avatar
Mark Friedrichs committed
115
                if( bExclusionFlag ){
116
117
                    getMaskedDScaleFactor( j, dScaleMask, &dScaleValue );
                    getMaskedPScaleFactor( j, pScaleMask, &pScaleValue );
Mark Friedrichs's avatar
Mark Friedrichs committed
118
                }
119

Mark Friedrichs's avatar
Mark Friedrichs committed
120
121
                float ijField[4][3];
                calculateFixedFieldRealSpacePairIxn_kernel( localParticle, psA[j], dScaleValue, pScaleValue, ijField
122
#ifdef AMOEBA_DEBUG
Mark Friedrichs's avatar
Mark Friedrichs committed
123
                                            , pullBack
124
#endif
Mark Friedrichs's avatar
Mark Friedrichs committed
125
                );
126

Mark Friedrichs's avatar
Mark Friedrichs committed
127
128
                // nan*0.0 = nan not 0.0, so explicitly exclude (atomI == atomJ) contribution
                // by setting match flag
129

Mark Friedrichs's avatar
Mark Friedrichs committed
130
                unsigned int match      = ( (atomI == (y + j)) || (atomI >= cAmoebaSim.numberOfAtoms) || ((y+j) >= cAmoebaSim.numberOfAtoms) ) ? 1 : 0;
131

Mark Friedrichs's avatar
Mark Friedrichs committed
132
                // add to field at atomI the field due atomJ's charge/dipole/quadrupole
133

Mark Friedrichs's avatar
Mark Friedrichs committed
134
135
136
                fieldSum[0]            += match ? 0.0f : ijField[0][0];
                fieldSum[1]            += match ? 0.0f : ijField[0][1];
                fieldSum[2]            += match ? 0.0f : ijField[0][2];
137

Mark Friedrichs's avatar
Mark Friedrichs committed
138
139
140
                fieldPolarSum[0]       += match ? 0.0f : ijField[2][0];
                fieldPolarSum[1]       += match ? 0.0f : ijField[2][1];
                fieldPolarSum[2]       += match ? 0.0f : ijField[2][2];
141
142

#ifdef AMOEBA_DEBUG
Mark Friedrichs's avatar
Mark Friedrichs committed
143
if( atomI == targetAtom || targetAtom == (y+j) ){
Mark Friedrichs's avatar
Mark Friedrichs committed
144
145
146
147
    unsigned int index                 = atomI == targetAtom ? (y + j) : atomI;
    unsigned int indexI                = 0;
    unsigned int indexJ                = indexI ? 0 : 2;
    unsigned int indices[4]            = { indexI, indexJ, indexI+1, indexJ+1 };
Mark Friedrichs's avatar
Mark Friedrichs committed
148
    float flag                         = 7.0f;
Mark Friedrichs's avatar
Mark Friedrichs committed
149
150
151
152
153
154

    debugArray[index].x                = (float) atomI;
    debugArray[index].y                = (float) (y + j);
    debugArray[index].z                = dScaleValue;
    debugArray[index].w                = pScaleValue;

Mark Friedrichs's avatar
Mark Friedrichs committed
155
156
157
158
159
160
161
162
163
164
165
166
    index                             += cAmoebaSim.paddedNumberOfAtoms;
    debugArray[index].x                = (float) bExclusionFlag;
    debugArray[index].y                = (float) (tgx);
    debugArray[index].z                = (float) j;
    debugArray[index].w                = flag;

    index                             += cAmoebaSim.paddedNumberOfAtoms;
    debugArray[index].x                = (float) dScaleMask;
    debugArray[index].y                = (float) pScaleMask.x;
    debugArray[index].z                = (float) pScaleMask.y;
    debugArray[index].w                = flag;

Mark Friedrichs's avatar
Mark Friedrichs committed
167
    for( int ii = 0; ii < 4; ii++ ){
Mark Friedrichs's avatar
Mark Friedrichs committed
168
        index                             += cAmoebaSim.paddedNumberOfAtoms;
Mark Friedrichs's avatar
Mark Friedrichs committed
169
170
171
172
173
        debugArray[index].x                = match ? 0.0f : ijField[indices[ii]][0];
        debugArray[index].y                = match ? 0.0f : ijField[indices[ii]][1];
        debugArray[index].z                = match ? 0.0f : ijField[indices[ii]][2];
        debugArray[index].w                = flag;
    }
Mark Friedrichs's avatar
Mark Friedrichs committed
174

Mark Friedrichs's avatar
Mark Friedrichs committed
175
176
177
178
179
180
181
    for( int pullIndex = 0; pullIndex < maxPullIndex; pullIndex++ ){
        index                             += cAmoebaSim.paddedNumberOfAtoms;
        debugArray[index].x                = pullBack[pullIndex].x;
        debugArray[index].y                = pullBack[pullIndex].y;
        debugArray[index].z                = pullBack[pullIndex].z;
        debugArray[index].w                = pullBack[pullIndex].w;
    }   
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199


}
#endif
            }

            // Write results

#ifdef USE_OUTPUT_BUFFER_PER_WARP
            unsigned int offset                 = 3*(x + tgx + warp*cAmoebaSim.paddedNumberOfAtoms);
            load3dArrayBufferPerWarp( offset, fieldSum,       outputEField );
            load3dArrayBufferPerWarp( offset, fieldPolarSum,  outputEFieldPolar );
#else
            unsigned int offset                 = 3*(x + tgx + (x >> GRIDBITS) * cAmoebaSim.paddedNumberOfAtoms);
            load3dArray( offset, fieldSum,       outputEField );
            load3dArray( offset, fieldPolarSum,  outputEFieldPolar );
#endif

Mark Friedrichs's avatar
Mark Friedrichs committed
200
        } else {
201

Mark Friedrichs's avatar
Mark Friedrichs committed
202
203
204
205
206
207
208
209
            if (lasty != y ) {
    
                // load coordinates, charge, ...
    
                loadFixedFieldShared( &(sA[threadIdx.x]), (y+tgx) );
    
            }

Mark Friedrichs's avatar
Mark Friedrichs committed
210
211
212
213
            unsigned int flags = cSim.pInteractionFlag[pos];
            if (flags == 0) {
                // No interactions in this block.
            } else {
214

Mark Friedrichs's avatar
Mark Friedrichs committed
215
                // zero shared fields
216

Mark Friedrichs's avatar
Mark Friedrichs committed
217
                zeroFixedFieldParticleSharedField( &(sA[threadIdx.x]) );
218

Mark Friedrichs's avatar
Mark Friedrichs committed
219
220
221
222
223
224
225
226
                if( bExclusionFlag ) {
                    unsigned int xi   = x >> GRIDBITS;
                    unsigned int yi   = y >> GRIDBITS;
                    unsigned int cell = xi+yi*cAmoebaSim.paddedNumberOfAtoms/GRID-yi*(yi+1)/2;
                    dScaleMask        = cAmoebaSim.pD_ScaleIndices[cAmoebaSim.pScaleIndicesIndex[cell]+tgx];
                    pScaleMask        = cAmoebaSim.pP_ScaleIndices[cAmoebaSim.pScaleIndicesIndex[cell]+tgx];
                } else {
                    dScaleValue = pScaleValue  = 1.0f;
227
228
                }

Mark Friedrichs's avatar
Mark Friedrichs committed
229
                for (unsigned int j = 0; j < GRID; j++){
230

Mark Friedrichs's avatar
Mark Friedrichs committed
231
232
233
234
235
                    unsigned int jIdx = (flags == 0xFFFFFFFF) ? tj : j;
                    if( bExclusionFlag ){
                        getMaskedDScaleFactor( jIdx, dScaleMask, &dScaleValue );
                        getMaskedPScaleFactor( jIdx, pScaleMask, &pScaleValue );
                    }
236
237

                    float ijField[4][3];
Mark Friedrichs's avatar
Mark Friedrichs committed
238
                    calculateFixedFieldRealSpacePairIxn_kernel( localParticle, psA[jIdx], dScaleValue, pScaleValue, ijField
239
240
241
242
243
#ifdef AMOEBA_DEBUG
                                                , pullBack
#endif
                    );

Mark Friedrichs's avatar
Mark Friedrichs committed
244
                    unsigned int outOfBounds     = ( (atomI >= cAmoebaSim.numberOfAtoms) || ((y+jIdx) >= cAmoebaSim.numberOfAtoms) ) ? 1 : 0;
Mark Friedrichs's avatar
Mark Friedrichs committed
245

246
                    // add to field at atomI the field due atomJ's charge/dipole/quadrupole
Mark Friedrichs's avatar
Mark Friedrichs committed
247
    
Mark Friedrichs's avatar
Mark Friedrichs committed
248
249
250
                    fieldSum[0]                 += outOfBounds ? 0.0f : ijField[0][0];
                    fieldSum[1]                 += outOfBounds ? 0.0f : ijField[0][1];
                    fieldSum[2]                 += outOfBounds ? 0.0f : ijField[0][2];
251

Mark Friedrichs's avatar
Mark Friedrichs committed
252
253
254
                    fieldPolarSum[0]            += outOfBounds ? 0.0f : ijField[2][0];
                    fieldPolarSum[1]            += outOfBounds ? 0.0f : ijField[2][1];
                    fieldPolarSum[2]            += outOfBounds ? 0.0f : ijField[2][2];
Mark Friedrichs's avatar
Mark Friedrichs committed
255
256
    
                    if( flags == 0xFFFFFFFF ){
257

Mark Friedrichs's avatar
Mark Friedrichs committed
258
259
260
261
262
263
264
265
266
267
268
                        // add to field at atomJ the field due atomI's charge/dipole/quadrupole
    
                        psA[jIdx].eField[0]        += outOfBounds ? 0.0f : ijField[1][0];
                        psA[jIdx].eField[1]        += outOfBounds ? 0.0f : ijField[1][1];
                        psA[jIdx].eField[2]        += outOfBounds ? 0.0f : ijField[1][2];
    
                        psA[jIdx].eFieldP[0]       += outOfBounds ? 0.0f : ijField[3][0];
                        psA[jIdx].eFieldP[1]       += outOfBounds ? 0.0f : ijField[3][1];
                        psA[jIdx].eFieldP[2]       += outOfBounds ? 0.0f : ijField[3][2];
 
                    } else {
269

Mark Friedrichs's avatar
Mark Friedrichs committed
270
271
272
                        sA[threadIdx.x].tempBuffer[0]  = outOfBounds ? 0.0f : ijField[1][0];
                        sA[threadIdx.x].tempBuffer[1]  = outOfBounds ? 0.0f : ijField[1][1];
                        sA[threadIdx.x].tempBuffer[2]  = outOfBounds ? 0.0f : ijField[1][2];
Mark Friedrichs's avatar
Mark Friedrichs committed
273
    
Mark Friedrichs's avatar
Mark Friedrichs committed
274
275
276
                        sA[threadIdx.x].tempBufferP[0] = outOfBounds ? 0.0f : ijField[3][0];
                        sA[threadIdx.x].tempBufferP[1] = outOfBounds ? 0.0f : ijField[3][1];
                        sA[threadIdx.x].tempBufferP[2] = outOfBounds ? 0.0f : ijField[3][2];
Mark Friedrichs's avatar
Mark Friedrichs committed
277
278

                        if( tgx % 2 == 0 ){
Mark Friedrichs's avatar
Mark Friedrichs committed
279
                            sumTempBuffer( sA[threadIdx.x], sA[threadIdx.x+1] ); 
Mark Friedrichs's avatar
Mark Friedrichs committed
280
281
                        } 
                        if( tgx % 4 == 0 ){
Mark Friedrichs's avatar
Mark Friedrichs committed
282
                            sumTempBuffer( sA[threadIdx.x], sA[threadIdx.x+2] ); 
Mark Friedrichs's avatar
Mark Friedrichs committed
283
284
                        } 
                        if( tgx % 8 == 0 ){
Mark Friedrichs's avatar
Mark Friedrichs committed
285
                            sumTempBuffer( sA[threadIdx.x], sA[threadIdx.x+4] ); 
Mark Friedrichs's avatar
Mark Friedrichs committed
286
287
                        } 
                        if( tgx % 16 == 0 ){
Mark Friedrichs's avatar
Mark Friedrichs committed
288
                            sumTempBuffer( sA[threadIdx.x], sA[threadIdx.x+8] ); 
Mark Friedrichs's avatar
Mark Friedrichs committed
289
290
291
292
                        } 

                        if (tgx == 0)
                        {
Mark Friedrichs's avatar
Mark Friedrichs committed
293
294
295
                            psA[jIdx].eField[0]  += sA[threadIdx.x].tempBuffer[0]  + sA[threadIdx.x+16].tempBuffer[0];
                            psA[jIdx].eField[1]  += sA[threadIdx.x].tempBuffer[1]  + sA[threadIdx.x+16].tempBuffer[1];
                            psA[jIdx].eField[2]  += sA[threadIdx.x].tempBuffer[2]  + sA[threadIdx.x+16].tempBuffer[2];
Mark Friedrichs's avatar
Mark Friedrichs committed
296

Mark Friedrichs's avatar
Mark Friedrichs committed
297
298
299
                            psA[jIdx].eFieldP[0] += sA[threadIdx.x].tempBufferP[0] + sA[threadIdx.x+16].tempBufferP[0];
                            psA[jIdx].eFieldP[1] += sA[threadIdx.x].tempBufferP[1] + sA[threadIdx.x+16].tempBufferP[1];
                            psA[jIdx].eFieldP[2] += sA[threadIdx.x].tempBufferP[2] + sA[threadIdx.x+16].tempBufferP[2];
Mark Friedrichs's avatar
Mark Friedrichs committed
300
301
                        }
                    }
302
303

#ifdef AMOEBA_DEBUG
Mark Friedrichs's avatar
Mark Friedrichs committed
304
if( (atomI == targetAtom || (y + jIdx) == targetAtom) ){
305

Mark Friedrichs's avatar
Mark Friedrichs committed
306
            unsigned int index                 = (atomI == targetAtom) ? (y + jIdx) : atomI;
307
308
309
310
            unsigned int indexI                = (atomI == targetAtom) ? 0 : 2;
            unsigned int indexJ                = (atomI == targetAtom) ? 2 : 0;

            debugArray[index].x                = (float) atomI;
Mark Friedrichs's avatar
Mark Friedrichs committed
311
312
313
            debugArray[index].y                = (float) (y + jIdx);
            debugArray[index].z                = dScaleValue;
            debugArray[index].w                = pScaleValue;
314
315

            float flag                         = 9.0f;
Mark Friedrichs's avatar
Mark Friedrichs committed
316
317
318
319
320
321
322
323
324
325
326
327
            index                             += cAmoebaSim.paddedNumberOfAtoms;
            debugArray[index].x                = (float) bExclusionFlag;
            debugArray[index].y                = (float) (tgx);
            debugArray[index].z                = (float) j;
            debugArray[index].w                = jIdx;
        
            index                             += cAmoebaSim.paddedNumberOfAtoms;
            debugArray[index].x                = (float) dScaleMask;
            debugArray[index].y                = (float) pScaleMask.x;
            debugArray[index].z                = (float) pScaleMask.y;
            debugArray[index].w                = (float) flags;
        
328
            index                             += cAmoebaSim.paddedNumberOfAtoms;
Mark Friedrichs's avatar
Mark Friedrichs committed
329
330
331
332
            debugArray[index].x                =  outOfBounds ? 0.0f : ijField[indexI][0];
            debugArray[index].y                =  outOfBounds ? 0.0f : ijField[indexI][1];
            debugArray[index].z                =  outOfBounds ? 0.0f : ijField[indexI][2];
            debugArray[index].w                =  flag;
333
334

            index                             += cAmoebaSim.paddedNumberOfAtoms;
Mark Friedrichs's avatar
Mark Friedrichs committed
335
336
337
            debugArray[index].x                =  outOfBounds ? 0.0f : ijField[indexJ][0];
            debugArray[index].y                =  outOfBounds ? 0.0f : ijField[indexJ][1];
            debugArray[index].z                =  outOfBounds ? 0.0f : ijField[indexJ][2];
338
339
340
            debugArray[index].w                = flag;

            index                             += cAmoebaSim.paddedNumberOfAtoms;
Mark Friedrichs's avatar
Mark Friedrichs committed
341
342
343
            debugArray[index].x                =  outOfBounds ? 0.0f : ijField[indexI+1][0];
            debugArray[index].y                =  outOfBounds ? 0.0f : ijField[indexI+1][1];
            debugArray[index].z                =  outOfBounds ? 0.0f : ijField[indexI+1][2];
344
345
346
            debugArray[index].w                = flag;

            index                             += cAmoebaSim.paddedNumberOfAtoms;
Mark Friedrichs's avatar
Mark Friedrichs committed
347
348
349
            debugArray[index].x                =  outOfBounds ? 0.0f : ijField[indexJ+1][0];
            debugArray[index].y                =  outOfBounds ? 0.0f : ijField[indexJ+1][1];
            debugArray[index].z                =  outOfBounds ? 0.0f : ijField[indexJ+1][2];
350
351
            debugArray[index].w                = flag;

Mark Friedrichs's avatar
Mark Friedrichs committed
352
353
354
355
356
357
358
359
            for( int pullIndex = 0; pullIndex < maxPullIndex; pullIndex++ ){
                index                             += cAmoebaSim.paddedNumberOfAtoms;
                debugArray[index].x                = pullBack[pullIndex].x;
                debugArray[index].y                = pullBack[pullIndex].y;
                debugArray[index].z                = pullBack[pullIndex].z;
                debugArray[index].w                = pullBack[pullIndex].w;
            }
}        
360
361
362
#endif
                    tj                  = (tj + 1) & (GRID - 1);

Mark Friedrichs's avatar
Mark Friedrichs committed
363
364
365
366
                } // j-loop block
    
                // Write results
    
367
#ifdef USE_OUTPUT_BUFFER_PER_WARP
Mark Friedrichs's avatar
Mark Friedrichs committed
368
369
370
371
372
373
374
375
                unsigned int offset                 = 3*(x + tgx + warp*cAmoebaSim.paddedNumberOfAtoms);
                load3dArrayBufferPerWarp( offset, fieldSum,       outputEField );
                load3dArrayBufferPerWarp( offset, fieldPolarSum,  outputEFieldPolar );
    
                offset                              = 3*(y + tgx + warp*cAmoebaSim.paddedNumberOfAtoms);
                load3dArrayBufferPerWarp( offset, sA[threadIdx.x].eField,  outputEField );
                load3dArrayBufferPerWarp( offset, sA[threadIdx.x].eFieldP, outputEFieldPolar );
    
376
#else
Mark Friedrichs's avatar
Mark Friedrichs committed
377
378
379
380
381
382
383
384
                unsigned int offset                 = 3*(x + tgx + (y >> GRIDBITS) * cAmoebaSim.paddedNumberOfAtoms);
                load3dArray( offset, fieldSum,       outputEField );
                load3dArray( offset, fieldPolarSum,  outputEFieldPolar );
    
                offset                              = 3*(y + tgx + (x >> GRIDBITS) * cAmoebaSim.paddedNumberOfAtoms);
                load3dArray( offset, sA[threadIdx.x].eField,  outputEField );
                load3dArray( offset, sA[threadIdx.x].eFieldP, outputEFieldPolar );
     
385
#endif
Mark Friedrichs's avatar
Mark Friedrichs committed
386
            } // end of pInteractionFlag block 
387
            lasty = y;
Mark Friedrichs's avatar
Mark Friedrichs committed
388
        } // x == y block
389
390
391
392

        pos++;
    }
}