PlasCom2  1.0
XPACC Multi-physics simluation application
Operators.f90
Go to the documentation of this file.
1 MODULE operators
2 
3  IMPLICIT NONE
4 
5 CONTAINS
6 
9 
34  SUBROUTINE applyoperator(numDim,dimSizes,numComponents,numPoints,opDir,opInterval,numStencils, &
35  stencilSizes,stencilStarts,numValues,stencilWeights,stencilOffsets,stencilID,U,dU)
36 
37  IMPLICIT NONE
38 
39  INTEGER(KIND=4), INTENT(IN) :: numDim, opDir,numStencils, numValues, numComponents
40  INTEGER(KIND=8), INTENT(IN) :: dimSizes(numdim), numPoints
41  INTEGER(KIND=8), INTENT(IN) :: opInterval(2*numdim)
42  INTEGER(KIND=4), INTENT(IN) :: stencilSizes(numstencils),stencilStarts(numstencils)
43  INTEGER(KIND=4), INTENT(IN) :: stencilOffsets(numvalues)
44  REAL(KIND=8), INTENT(IN) :: stencilWeights(numvalues)
45  INTEGER(KIND=4), INTENT(IN) :: stencilID(numpoints)
46 
47  REAL(KIND=8), INTENT(IN), TARGET :: U(numpoints*numcomponents)
48  REAL(KIND=8), INTENT(OUT), TARGET :: dU(numpoints*numcomponents)
49 
50  REAL(KIND=8) :: fac
51  INTEGER(KIND=4) :: iStencil, iWeight, iComp
52  INTEGER(KIND=8) :: plane, pointOffset, compOffset
53  INTEGER(KIND=8) :: I, J, K, jIndex, jkIndex, kIndex, iPoint
54 
55  REAL(KIND=8), DIMENSION(:), POINTER :: uCompPtr
56  REAL(KIND=8), DIMENSION(:), POINTER :: duCompPtr
57 
58  ! DO iComp = 1, numComponents
59 
60  ! compOffset = (iComp-1)*numPoints
61  ! uCompPtr => U(compOffset+1:compOffset+numPoints)
62  ! duCompPtr => dU(compOffset+1:compOffset+numPoints)
63 
64  IF(numdim == 1) THEN
65  DO i = opinterval(1), opinterval(2)
66  istencil = stencilid(i)
67  du(i) = 0.0_8
68  DO iweight = stencilstarts(istencil),stencilstarts(istencil)+stencilsizes(istencil) - 1
69  du(i) = du(i) + stencilweights(iweight)*u(i+stenciloffsets(iweight))
70  END DO
71  END DO
72  ELSE IF(numdim == 2) THEN
73  IF(opdir == 1) THEN
74  DO j = opinterval(3), opinterval(4)
75  jindex = (j-1)*dimsizes(1)
76  DO i = opinterval(1), opinterval(2)
77  ipoint = jindex + i
78  istencil = stencilid(ipoint)
79  du(ipoint) = 0.0_8
80  DO iweight = stencilstarts(istencil),stencilstarts(istencil)+stencilsizes(istencil) - 1
81  du(ipoint) = du(ipoint) + &
82  stencilweights(iweight)*u(ipoint+stenciloffsets(iweight))
83  END DO
84  END DO
85  END DO
86  ELSE IF(opdir == 2) THEN
87  plane = dimsizes(1)
88  DO j = opinterval(3), opinterval(4)
89  jindex = (j-1)*plane
90  DO i = opinterval(1), opinterval(2)
91  ipoint = jindex + i
92  istencil = stencilid(ipoint)
93  du(ipoint) = 0.0_8
94  DO iweight = stencilstarts(istencil),stencilstarts(istencil)+stencilsizes(istencil) - 1
95  du(ipoint) = du(ipoint) + &
96  stencilweights(iweight)*u(ipoint+stenciloffsets(iweight)*plane)
97  END DO
98  END DO
99  END DO
100  ENDIF
101  ELSE IF(numdim == 3) THEN
102  plane = dimsizes(1) * dimsizes(2)
103  IF(opdir == 1) THEN
104  DO k = opinterval(5), opinterval(6)
105  kindex = (k-1)*plane
106  DO j = opinterval(3), opinterval(4)
107  jkindex = kindex + (j-1)*dimsizes(1)
108  DO i = opinterval(1), opinterval(2)
109  ipoint = jkindex + i
110  istencil = stencilid(ipoint)
111  du(ipoint) = 0.0_8
112  DO iweight = stencilstarts(istencil),stencilstarts(istencil)+stencilsizes(istencil) - 1
113  du(ipoint) = du(ipoint) + &
114  stencilweights(iweight)*u(ipoint+stenciloffsets(iweight))
115  END DO
116  END DO
117  END DO
118  END DO
119  ELSE
120  pointoffset = 1
121  IF(opdir == 2) THEN
122  pointoffset = dimsizes(1)
123  ELSE
124  pointoffset = plane
125  ENDIF
126  DO k = opinterval(5), opinterval(6)
127  kindex = (k-1)*plane
128  DO j = opinterval(3), opinterval(4)
129  jkindex = kindex + (j-1)*dimsizes(1)
130  DO i = opinterval(1), opinterval(2)
131  ipoint = jkindex + i
132  istencil = stencilid(ipoint)
133  du(ipoint) = 0.0_8
134  DO iweight = stencilstarts(istencil),stencilstarts(istencil)+stencilsizes(istencil) - 1
135  du(ipoint) = du(ipoint) + &
136  stencilweights(iweight)*u(ipoint+stenciloffsets(iweight)*pointoffset)
137  END DO ! weight
138  END DO ! I
139  END DO ! J
140  END DO ! K
141  ENDIF ! (opDir)
142  ENDIF ! (numDim)
143  END SUBROUTINE applyoperator
144 
147 
172  SUBROUTINE applyoperatorv(numDim,dimSizes,numComponents,numPoints,opDir,opInterval,numStencils, &
173  stencilSizes,stencilStarts,numValues,stencilWeights,stencilOffsets,stencilID,U,dU)
174 
175  IMPLICIT NONE
176 
177  INTEGER(KIND=4), INTENT(IN) :: numDim, opDir,numStencils, numValues, numComponents
178  INTEGER(KIND=8), INTENT(IN) :: dimSizes(numdim), numPoints
179  INTEGER(KIND=8), INTENT(IN) :: opInterval(2*numdim)
180  INTEGER(KIND=4), INTENT(IN) :: stencilSizes(numstencils),stencilStarts(numstencils)
181  INTEGER(KIND=4), INTENT(IN) :: stencilOffsets(numvalues)
182  REAL(KIND=8), INTENT(IN) :: stencilWeights(numvalues)
183  INTEGER(KIND=4), INTENT(IN) :: stencilID(numpoints)
184 
185  REAL(KIND=8), INTENT(IN), TARGET :: U(numpoints*numcomponents)
186  REAL(KIND=8), INTENT(OUT), TARGET :: dU(numpoints*numcomponents)
187 
188  REAL(KIND=8) :: fac
189  INTEGER(KIND=4) :: iStencil, iWeight, iComp
190  INTEGER(KIND=8) :: plane, pointOffset, compOffset
191  INTEGER(KIND=8) :: I, J, K, jIndex, jkIndex, kIndex, iPoint
192 
193  REAL(KIND=8), DIMENSION(:), POINTER :: uCompPtr
194  REAL(KIND=8), DIMENSION(:), POINTER :: duCompPtr
195 
196  DO icomp = 1, numcomponents
197 
198  compoffset = (icomp-1)*numpoints
199  ucompptr => u(compoffset+1:compoffset+numpoints)
200  ducompptr => du(compoffset+1:compoffset+numpoints)
201 
202  IF(numdim == 1) THEN
203  DO i = opinterval(1), opinterval(2)
204  istencil = stencilid(i)
205  ducompptr(i) = 0.0_8
206  DO iweight = stencilstarts(istencil),stencilstarts(istencil)+stencilsizes(istencil) - 1
207  ducompptr(i) = ducompptr(i) + stencilweights(iweight)*ucompptr(i+stenciloffsets(iweight))
208  END DO
209  END DO
210  ELSE IF(numdim == 2) THEN
211  IF(opdir == 1) THEN
212  DO j = opinterval(3), opinterval(4)
213  jindex = (j-1)*dimsizes(1)
214  DO i = opinterval(1), opinterval(2)
215  ipoint = jindex + i
216  istencil = stencilid(ipoint)
217  ducompptr(ipoint) = 0.0_8
218  DO iweight = stencilstarts(istencil),stencilstarts(istencil)+stencilsizes(istencil) - 1
219  ducompptr(ipoint) = ducompptr(ipoint) + &
220  stencilweights(iweight)*ucompptr(ipoint+stenciloffsets(iweight))
221  END DO
222  END DO
223  END DO
224  ELSE IF(opdir == 2) THEN
225  plane = dimsizes(1)
226  DO j = opinterval(3), opinterval(4)
227  jindex = (j-1)*plane
228  DO i = opinterval(1), opinterval(2)
229  ipoint = jindex + i
230  istencil = stencilid(ipoint)
231  ducompptr(ipoint) = 0.0_8
232  DO iweight = stencilstarts(istencil),stencilstarts(istencil)+stencilsizes(istencil) - 1
233  ducompptr(ipoint) = ducompptr(ipoint) + &
234  stencilweights(iweight)*ucompptr(ipoint+stenciloffsets(iweight)*plane)
235  END DO
236  END DO
237  END DO
238  ENDIF
239  ELSE IF(numdim == 3) THEN
240  plane = dimsizes(1) * dimsizes(2)
241  IF(opdir == 1) THEN
242  DO k = opinterval(5), opinterval(6)
243  kindex = (k-1)*plane
244  DO j = opinterval(3), opinterval(4)
245  jkindex = kindex + (j-1)*dimsizes(1)
246  DO i = opinterval(1), opinterval(2)
247  ipoint = jkindex + i
248  istencil = stencilid(ipoint)
249  ducompptr(ipoint) = 0.0_8
250  DO iweight = stencilstarts(istencil),stencilstarts(istencil)+stencilsizes(istencil) - 1
251  ducompptr(ipoint) = ducompptr(ipoint) + &
252  stencilweights(iweight)*ucompptr(ipoint+stenciloffsets(iweight))
253  END DO
254  END DO
255  END DO
256  END DO
257  ELSE
258  pointoffset = 1
259  IF(opdir == 2) THEN
260  pointoffset = dimsizes(1)
261  ELSE
262  pointoffset = plane
263  ENDIF
264  DO k = opinterval(5), opinterval(6)
265  kindex = (k-1)*plane
266  DO j = opinterval(3), opinterval(4)
267  jkindex = kindex + (j-1)*dimsizes(1)
268  DO i = opinterval(1), opinterval(2)
269  ipoint = jkindex + i
270  istencil = stencilid(ipoint)
271  ducompptr(ipoint) = 0.0_8
272  DO iweight = stencilstarts(istencil),stencilstarts(istencil)+stencilsizes(istencil) - 1
273  ducompptr(ipoint) = ducompptr(ipoint) + &
274  stencilweights(iweight)*ucompptr(ipoint+stenciloffsets(iweight)*pointoffset)
275  END DO ! weight
276  END DO ! I
277  END DO ! J
278  END DO ! K
279  ENDIF ! (opDir)
280  ENDIF ! (numDim)
281  END DO ! (loop over components)
282  END SUBROUTINE applyoperatorv
283 
309  SUBROUTINE applyoperatorblobs(numDim,dimSizes,numComponents,numPointsBuffer,opDir,numStencils, &
310  stencilSizes,stencilStarts,numStencilValues,stencilWeights,stencilOffsets,numPointsStencil,&
311  numPointsApply,stencilPoints,U,dU)
312 
313  IMPLICIT NONE
314 
315  INTEGER(KIND=4), INTENT(IN) :: numDim, opDir,numStencils, numStencilValues, numComponents
316  INTEGER(KIND=8), INTENT(IN) :: numPointsApply, numPointsBuffer
317  INTEGER(KIND=8), INTENT(IN) :: dimSizes(numdim),numPointsStencil(numstencils)
318  INTEGER(KIND=8), INTENT(IN) :: stencilPoints(numpointsapply)
319  INTEGER(KIND=4), INTENT(IN) :: stencilSizes(numstencils),stencilStarts(numstencils)
320  INTEGER(KIND=4), INTENT(IN) :: stencilOffsets(numstencilvalues)
321  REAL(KIND=8), INTENT(IN) :: stencilWeights(numstencilvalues)
322  REAL(KIND=8), INTENT(IN) :: U(numpointsbuffer)
323  REAL(KIND=8), INTENT(OUT) :: dU(numpointsbuffer)
324 
325  INTEGER(KIND=8) :: stencilPointsOffset, iPoint, pointsStart,pointsEnd,numPointsThisStencil
326  INTEGER(KIND=4) :: iStencil, iWeight, stencilSize, stencilStart, stencilEnd
327 
328  ! dU(:) = 0.0_8
329  stencilpointsoffset = 1
330  ! WRITE(*,*) 'NumDim: ',numDim
331  ! WRITE(*,*) 'dimSizes:',dimSizes
332  ! WRITE(*,*) 'numComponents:',numComponents
333  ! WRITE(*,*) 'numPointsBuffer:',numPointsBuffer
334  ! WRITE(*,*) 'opDir:',opDir
335  ! WRITE(*,*) 'numStencils:',numStencils
336  ! WRITE(*,*) 'stencilSizes:',stencilSizes
337  ! WRITE(*,*) 'stencilStarts:',stencilStarts
338  ! WRITE(*,*) 'numStencilValues:',numStencilValues
339  ! WRITE(*,*) 'stencilWeights:',stencilWeights
340  ! WRITE(*,*) 'stencilOffsets:',stencilOffsets
341  ! WRITE(*,*) 'numPointsStencil:',numPointsStencil
342  ! WRITE(*,*) 'numPointsApply:',numPointsApply
343  DO istencil = 1, numstencils
344  numpointsthisstencil = numpointsstencil(istencil)
345  stencilsize = stencilsizes(istencil)
346  stencilstart = stencilstarts(istencil)
347  stencilend = stencilstart+stencilsize-1
348  pointsstart = stencilpointsoffset
349  pointsend = pointsstart + numpointsthisstencil-1
350  CALL applysinglestencil(numdim,dimsizes,numcomponents,numpointsbuffer,opdir,numpointsthisstencil, &
351  stencilpoints(pointsstart:pointsend),stencilsize,stencilweights(stencilstart:stencilend),&
352  stenciloffsets(stencilstart:stencilend),u,du)
353  stencilpointsoffset = stencilpointsoffset+numpointsthisstencil
354  END DO
355 
356  END SUBROUTINE applyoperatorblobs
357 
358 
359 
380  SUBROUTINE applysinglestencil(numDim,dimSizes,numComponents,numPointsBuffer,opDir,numPointsApply, &
381  applyPoints,stencilSize,stencilWeights,stencilOffsets,U,dU)
382 
383  IMPLICIT NONE
384 
385  INTEGER(KIND=4), INTENT(IN) :: numDim, opDir, stencilSize, numComponents
386  INTEGER(KIND=8), INTENT(IN) :: dimSizes(numdim), numPointsBuffer, numPointsApply
387  INTEGER(KIND=8), INTENT(IN) :: applyPoints(numpointsapply)
388  INTEGER(KIND=4), INTENT(IN) :: stencilOffsets(stencilsize)
389  REAL(KIND=8), INTENT(IN) :: stencilWeights(stencilsize)
390  REAL(KIND=8), INTENT(IN) :: U(numpointsbuffer)
391  REAL(KIND=8), INTENT(OUT) :: dU(numpointsbuffer)
392 
393  INTEGER(KIND=8) :: I, iPoint, plane
394  INTEGER(KIND=4) :: iStencil, iWeight
395 
396  IF(opdir == 1) THEN
397  DO i = 1,numpointsapply
398  ipoint = applypoints(i)
399  du(ipoint) = 0.0_8
400  DO iweight = 1, stencilsize
401  du(ipoint) = du(ipoint) + stencilweights(iweight)*u(ipoint+stenciloffsets(iweight))
402  END DO
403  END DO
404  ELSE
405  IF(opdir == 2) THEN
406  plane = dimsizes(1)
407  ELSE IF(opdir == 3) THEN
408  plane = dimsizes(1)*dimsizes(2)
409  ENDIF
410  DO i = 1, numpointsapply
411  ipoint = applypoints(i)
412  du(ipoint) = 0.0_8
413  DO iweight = 1,stencilsize
414  du(ipoint) = du(ipoint) + stencilweights(iweight)*u(ipoint+stenciloffsets(iweight)*plane)
415  END DO
416  END DO
417  END IF
418 
419  END SUBROUTINE applysinglestencil
420 
421 
423  SUBROUTINE yaxpy(numDim,numPoints,bufferSize,bufferInterval,a,X,Y)
424 
425  INTEGER(KIND=4), INTENT(IN) :: numDim
426  INTEGER(KIND=8), INTENT(IN) :: numPoints,bufferSize(numdim)
427  INTEGER(KIND=8), INTENT(IN) :: bufferInterval(2*numdim)
428  REAL(KIND=8), INTENT(IN) :: a,X(numpoints)
429  REAL(KIND=8), INTENT(INOUT) :: Y(numpoints)
430 
431  INTEGER(KIND=8) :: I, J, K
432  INTEGER(KIND=8) :: nPlane, zIndex, yIndex, yzIndex, bufferIndex, xSize
433  INTEGER(KIND=8) :: iStart,iEnd,jStart,jEnd,kStart,kEnd
434 
435  istart = bufferinterval(1)
436  iend = bufferinterval(2)
437  xsize = buffersize(1)
438 
439  IF(numdim == 1) THEN
440  DO i = istart, iend
441  y(i) = y(i) + a*x(i)
442  END DO
443  ELSE IF(numdim == 2) THEN
444  jstart = bufferinterval(3)
445  jend = bufferinterval(4)
446  DO j = jstart, jend
447  yindex = (j-1)*xsize
448  DO i = istart, iend
449  bufferindex = yindex + i
450  y(bufferindex) = y(bufferindex) + a*x(bufferindex)
451  END DO
452  END DO
453  ELSE IF(numdim == 3) THEN
454  jstart = bufferinterval(3)
455  jend = bufferinterval(4)
456  kstart = bufferinterval(5)
457  kend = bufferinterval(6)
458  nplane = xsize * buffersize(2)
459  DO k = kstart, kend
460  zindex = (k-1)*nplane
461  DO j = jstart, jend
462  yzindex = zindex + (j-1)*xsize
463  DO i = istart, iend
464  bufferindex = yzindex + i
465  y(bufferindex) = y(bufferindex) + a*x(bufferindex)
466  END DO
467  END DO
468  END DO
469  ENDIF
470 
471  END SUBROUTINE yaxpy
472 
474  SUBROUTINE zaxpy(numDim,numPoints,bufferSize,bufferInterval,a,X,Y,Z)
475 
476  INTEGER(KIND=4), INTENT(IN) :: numDim
477  INTEGER(KIND=8), INTENT(IN) :: numPoints,bufferSize(numdim)
478  INTEGER(KIND=8), INTENT(IN) :: bufferInterval(2*numdim)
479  REAL(KIND=8), INTENT(IN) :: a,X(numpoints),Y(numpoints)
480  REAL(KIND=8), INTENT(OUT) :: Z(numpoints)
481 
482  INTEGER(KIND=8) :: I, J, K
483  INTEGER(KIND=8) :: nPlane, zIndex, yIndex, yzIndex, bufferIndex, xSize
484  INTEGER(KIND=8) :: iStart,iEnd,jStart,jEnd,kStart,kEnd
485 
486  istart = bufferinterval(1)
487  iend = bufferinterval(2)
488  xsize = buffersize(1)
489 
490  IF(numdim == 1) THEN
491  DO i = istart, iend
492  z(i) = y(i) + a*x(i)
493  END DO
494  ELSE IF(numdim == 2) THEN
495  jstart = bufferinterval(3)
496  jend = bufferinterval(4)
497  DO j = jstart, jend
498  yindex = (j-1)*xsize
499  DO i = istart, iend
500  bufferindex = yindex + i
501  z(bufferindex) = y(bufferindex) + a*x(bufferindex)
502  END DO
503  END DO
504  ELSE IF(numdim == 3) THEN
505  jstart = bufferinterval(3)
506  jend = bufferinterval(4)
507  kstart = bufferinterval(5)
508  kend = bufferinterval(6)
509  nplane = xsize * buffersize(2)
510  DO k = kstart, kend
511  zindex = (k-1)*nplane
512  DO j = jstart, jend
513  yzindex = zindex + (j-1)*xsize
514  DO i = istart, iend
515  bufferindex = yzindex + i
516  z(bufferindex) = y(bufferindex) + a*x(bufferindex)
517  END DO
518  END DO
519  END DO
520  ENDIF
521 
522  END SUBROUTINE zaxpy
523 
525  SUBROUTINE zaxpby(numDim,numPoints,bufferSize,bufferInterval,a,b,X,Y,Z)
526 
527  INTEGER(KIND=4), INTENT(IN) :: numDim
528  INTEGER(KIND=8), INTENT(IN) :: numPoints,bufferSize(numdim)
529  INTEGER(KIND=8), INTENT(IN) :: bufferInterval(2*numdim)
530  REAL(KIND=8), INTENT(IN) :: a,b,X(numpoints),Y(numpoints)
531  REAL(KIND=8), INTENT(OUT) :: Z(numpoints)
532 
533  INTEGER(KIND=8) :: I, J, K
534  INTEGER(KIND=8) :: nPlane, zIndex, yIndex, yzIndex, bufferIndex, xSize
535  INTEGER(KIND=8) :: iStart,iEnd,jStart,jEnd,kStart,kEnd
536 
537  istart = bufferinterval(1)
538  iend = bufferinterval(2)
539  xsize = buffersize(1)
540 
541  IF(numdim == 1) THEN
542  DO i = istart, iend
543  z(i) = b*y(i) + a*x(i)
544  END DO
545  ELSE IF(numdim == 2) THEN
546  jstart = bufferinterval(3)
547  jend = bufferinterval(4)
548  DO j = jstart, jend
549  yindex = (j-1)*xsize
550  DO i = istart, iend
551  bufferindex = yindex + i
552  z(bufferindex) = b*y(bufferindex) + a*x(bufferindex)
553  END DO
554  END DO
555  ELSE IF(numdim == 3) THEN
556  jstart = bufferinterval(3)
557  jend = bufferinterval(4)
558  kstart = bufferinterval(5)
559  kend = bufferinterval(6)
560  nplane = xsize * buffersize(2)
561  DO k = kstart, kend
562  zindex = (k-1)*nplane
563  DO j = jstart, jend
564  yzindex = zindex + (j-1)*xsize
565  DO i = istart, iend
566  bufferindex = yzindex + i
567  z(bufferindex) = b*y(bufferindex) + a*x(bufferindex)
568  END DO
569  END DO
570  END DO
571  ENDIF
572 
573  END SUBROUTINE zaxpby
574 
576  SUBROUTINE yaxpby(numDim,numPoints,bufferSize,bufferInterval,a,b,X,Y)
577 
578  INTEGER(KIND=4), INTENT(IN) :: numDim
579  INTEGER(KIND=8), INTENT(IN) :: numPoints,bufferSize(numdim)
580  INTEGER(KIND=8), INTENT(IN) :: bufferInterval(2*numdim)
581  REAL(KIND=8), INTENT(IN) :: a,b,X(numpoints)
582  REAL(KIND=8), INTENT(INOUT) :: Y(numpoints)
583 
584  INTEGER(KIND=8) :: I, J, K
585  INTEGER(KIND=8) :: nPlane, zIndex, yIndex, yzIndex, bufferIndex, xSize
586  INTEGER(KIND=8) :: iStart,iEnd,jStart,jEnd,kStart,kEnd
587 
588  istart = bufferinterval(1)
589  iend = bufferinterval(2)
590  xsize = buffersize(1)
591 
592  IF(numdim == 1) THEN
593  DO i = istart, iend
594  y(i) = b*y(i) + a*x(i)
595  END DO
596  ELSE IF(numdim == 2) THEN
597  jstart = bufferinterval(3)
598  jend = bufferinterval(4)
599  DO j = jstart, jend
600  yindex = (j-1)*xsize
601  DO i = istart, iend
602  bufferindex = yindex + i
603  y(bufferindex) = b*y(bufferindex) + a*x(bufferindex)
604  END DO
605  END DO
606  ELSE IF(numdim == 3) THEN
607  jstart = bufferinterval(3)
608  jend = bufferinterval(4)
609  kstart = bufferinterval(5)
610  kend = bufferinterval(6)
611  nplane = xsize * buffersize(2)
612  DO k = kstart, kend
613  zindex = (k-1)*nplane
614  DO j = jstart, jend
615  yzindex = zindex + (j-1)*xsize
616  DO i = istart, iend
617  bufferindex = yzindex + i
618  y(bufferindex) = b*y(bufferindex) + a*x(bufferindex)
619  END DO
620  END DO
621  END DO
622  ENDIF
623 
624  END SUBROUTINE yaxpby
625 
627  SUBROUTINE yax(numDim,numPoints,bufferSize,bufferInterval,a,X,Y)
628 
629  INTEGER(KIND=4), INTENT(IN) :: numDim
630  INTEGER(KIND=8), INTENT(IN) :: numPoints,bufferSize(numdim)
631  INTEGER(KIND=8), INTENT(IN) :: bufferInterval(2*numdim)
632  REAL(KIND=8), INTENT(IN) :: a,X(numpoints)
633  REAL(KIND=8), INTENT(OUT) :: Y(numpoints)
634 
635  INTEGER(KIND=8) :: I, J, K
636  INTEGER(KIND=8) :: nPlane, zIndex, yIndex, yzIndex, bufferIndex, xSize
637  INTEGER(KIND=8) :: iStart,iEnd,jStart,jEnd,kStart,kEnd
638 
639  istart = bufferinterval(1)
640  iend = bufferinterval(2)
641  xsize = buffersize(1)
642 
643  IF(numdim == 1) THEN
644  DO i = istart, iend
645  y(i) = a*x(i)
646  END DO
647  ELSE IF(numdim == 2) THEN
648  jstart = bufferinterval(3)
649  jend = bufferinterval(4)
650  DO j = jstart, jend
651  yindex = (j-1)*xsize
652  DO i = istart, iend
653  bufferindex = yindex + i
654  y(bufferindex) = a*x(bufferindex)
655  END DO
656  END DO
657  ELSE IF(numdim == 3) THEN
658  jstart = bufferinterval(3)
659  jend = bufferinterval(4)
660  kstart = bufferinterval(5)
661  kend = bufferinterval(6)
662  nplane = xsize * buffersize(2)
663  DO k = kstart, kend
664  zindex = (k-1)*nplane
665  DO j = jstart, jend
666  yzindex = zindex + (j-1)*xsize
667  DO i = istart, iend
668  bufferindex = yzindex + i
669  y(bufferindex) = a*x(bufferindex)
670  END DO
671  END DO
672  END DO
673  ENDIF
674 
675  END SUBROUTINE yax
676 
678  SUBROUTINE zxy(numDim,numPoints,bufferSize,bufferInterval,X,Y,Z)
679 
680  INTEGER(KIND=4), INTENT(IN) :: numDim
681  INTEGER(KIND=8), INTENT(IN) :: numPoints,bufferSize(numdim)
682  INTEGER(KIND=8), INTENT(IN) :: bufferInterval(2*numdim)
683  REAL(KIND=8), INTENT(IN) :: X(numpoints)
684  REAL(KIND=8), INTENT(IN) :: Y(numpoints)
685  REAL(KIND=8), INTENT(OUT) :: Z(numpoints)
686 
687  INTEGER(KIND=8) :: I, J, K
688  INTEGER(KIND=8) :: nPlane, zIndex, yIndex, yzIndex, bufferIndex, xSize
689  INTEGER(KIND=8) :: iStart,iEnd,jStart,jEnd,kStart,kEnd
690 
691  istart = bufferinterval(1)
692  iend = bufferinterval(2)
693  xsize = buffersize(1)
694 
695  IF(numdim == 1) THEN
696  DO i = istart, iend
697  z(i) = x(i)*y(i)
698  END DO
699  ELSE IF(numdim == 2) THEN
700  jstart = bufferinterval(3)
701  jend = bufferinterval(4)
702  DO j = jstart, jend
703  yindex = (j-1)*xsize
704  DO i = istart, iend
705  bufferindex = yindex + i
706  z(bufferindex) = x(bufferindex)*y(bufferindex)
707  END DO
708  END DO
709  ELSE IF(numdim == 3) THEN
710  jstart = bufferinterval(3)
711  jend = bufferinterval(4)
712  kstart = bufferinterval(5)
713  kend = bufferinterval(6)
714  nplane = xsize * buffersize(2)
715  DO k = kstart, kend
716  zindex = (k-1)*nplane
717  DO j = jstart, jend
718  yzindex = zindex + (j-1)*xsize
719  DO i = istart, iend
720  bufferindex = yzindex + i
721  z(bufferindex) = x(bufferindex)*y(bufferindex)
722  END DO
723  END DO
724  END DO
725  ENDIF
726 
727  END SUBROUTINE zxy
728 
730  SUBROUTINE yxy(numDim,numPoints,bufferSize,bufferInterval,X,Y)
731 
732  INTEGER(KIND=4), INTENT(IN) :: numDim
733  INTEGER(KIND=8), INTENT(IN) :: numPoints,bufferSize(numdim)
734  INTEGER(KIND=8), INTENT(IN) :: bufferInterval(2*numdim)
735  REAL(KIND=8), INTENT(IN) :: X(numpoints)
736  REAL(KIND=8), INTENT(INOUT) :: Y(numpoints)
737 
738  INTEGER(KIND=8) :: I, J, K
739  INTEGER(KIND=8) :: nPlane, zIndex, yIndex, yzIndex, bufferIndex, xSize
740  INTEGER(KIND=8) :: iStart,iEnd,jStart,jEnd,kStart,kEnd
741 
742  istart = bufferinterval(1)
743  iend = bufferinterval(2)
744  xsize = buffersize(1)
745 
746  IF(numdim == 1) THEN
747  DO i = istart, iend
748  y(i) = x(i)*y(i)
749  END DO
750  ELSE IF(numdim == 2) THEN
751  jstart = bufferinterval(3)
752  jend = bufferinterval(4)
753  DO j = jstart, jend
754  yindex = (j-1)*xsize
755  DO i = istart, iend
756  bufferindex = yindex + i
757  y(bufferindex) = x(bufferindex)*y(bufferindex)
758  END DO
759  END DO
760  ELSE IF(numdim == 3) THEN
761  jstart = bufferinterval(3)
762  jend = bufferinterval(4)
763  kstart = bufferinterval(5)
764  kend = bufferinterval(6)
765  nplane = xsize * buffersize(2)
766  DO k = kstart, kend
767  zindex = (k-1)*nplane
768  DO j = jstart, jend
769  yzindex = zindex + (j-1)*xsize
770  DO i = istart, iend
771  bufferindex = yzindex + i
772  y(bufferindex) = x(bufferindex)*y(bufferindex)
773  END DO
774  END DO
775  END DO
776  ENDIF
777 
778  END SUBROUTINE yxy
779 
781  SUBROUTINE zwxpy(numDim,numPoints,bufferSize,bufferInterval,W,X,Y,Z)
782 
783  INTEGER(KIND=4), INTENT(IN) :: numDim
784  INTEGER(KIND=8), INTENT(IN) :: numPoints,bufferSize(numdim)
785  INTEGER(KIND=8), INTENT(IN) :: bufferInterval(2*numdim)
786  REAL(KIND=8), INTENT(IN) :: W(numpoints)
787  REAL(KIND=8), INTENT(IN) :: X(numpoints)
788  REAL(KIND=8), INTENT(IN) :: Y(numpoints)
789  REAL(KIND=8), INTENT(OUT) :: Z(numpoints)
790 
791  INTEGER(KIND=8) :: I, J, K
792  INTEGER(KIND=8) :: nPlane, zIndex, yIndex, yzIndex, bufferIndex, xSize
793  INTEGER(KIND=8) :: iStart,iEnd,jStart,jEnd,kStart,kEnd
794 
795  istart = bufferinterval(1)
796  iend = bufferinterval(2)
797  xsize = buffersize(1)
798 
799  IF(numdim == 1) THEN
800  DO i = istart, iend
801  z(i) = x(i)*w(i) + y(i)
802  END DO
803  ELSE IF(numdim == 2) THEN
804  jstart = bufferinterval(3)
805  jend = bufferinterval(4)
806  DO j = jstart, jend
807  yindex = (j-1)*xsize
808  DO i = istart, iend
809  bufferindex = yindex + i
810  z(bufferindex) = x(bufferindex)*w(bufferindex) + y(bufferindex)
811  END DO
812  END DO
813  ELSE IF(numdim == 3) THEN
814  jstart = bufferinterval(3)
815  jend = bufferinterval(4)
816  kstart = bufferinterval(5)
817  kend = bufferinterval(6)
818  nplane = xsize * buffersize(2)
819  DO k = kstart, kend
820  zindex = (k-1)*nplane
821  DO j = jstart, jend
822  yzindex = zindex + (j-1)*xsize
823  DO i = istart, iend
824  bufferindex = yzindex + i
825  z(bufferindex) = x(bufferindex)*w(bufferindex) + y(bufferindex)
826  END DO
827  END DO
828  END DO
829  ENDIF
830 
831  END SUBROUTINE zwxpy
832 
834  SUBROUTINE ywxpy(numDim,numPoints,bufferSize,bufferInterval,W,X,Y)
835 
836  INTEGER(KIND=4), INTENT(IN) :: numDim
837  INTEGER(KIND=8), INTENT(IN) :: numPoints,bufferSize(numdim)
838  INTEGER(KIND=8), INTENT(IN) :: bufferInterval(2*numdim)
839  REAL(KIND=8), INTENT(IN) :: W(numpoints)
840  REAL(KIND=8), INTENT(IN) :: X(numpoints)
841  REAL(KIND=8), INTENT(INOUT) :: Y(numpoints)
842 
843  INTEGER(KIND=8) :: I, J, K
844  INTEGER(KIND=8) :: nPlane, zIndex, yIndex, yzIndex, bufferIndex, xSize
845  INTEGER(KIND=8) :: iStart,iEnd,jStart,jEnd,kStart,kEnd
846 
847  istart = bufferinterval(1)
848  iend = bufferinterval(2)
849  xsize = buffersize(1)
850 
851  IF(numdim == 1) THEN
852  DO i = istart, iend
853  y(i) = x(i)*w(i) + y(i)
854  END DO
855  ELSE IF(numdim == 2) THEN
856  jstart = bufferinterval(3)
857  jend = bufferinterval(4)
858  DO j = jstart, jend
859  yindex = (j-1)*xsize
860  DO i = istart, iend
861  bufferindex = yindex + i
862  y(bufferindex) = x(bufferindex)*w(bufferindex) + y(bufferindex)
863  END DO
864  END DO
865  ELSE IF(numdim == 3) THEN
866  jstart = bufferinterval(3)
867  jend = bufferinterval(4)
868  kstart = bufferinterval(5)
869  kend = bufferinterval(6)
870  nplane = xsize * buffersize(2)
871  DO k = kstart, kend
872  zindex = (k-1)*nplane
873  DO j = jstart, jend
874  yzindex = zindex + (j-1)*xsize
875  DO i = istart, iend
876  bufferindex = yzindex + i
877  y(bufferindex) = x(bufferindex)*w(bufferindex) + y(bufferindex)
878  END DO
879  END DO
880  END DO
881  ENDIF
882 
883  END SUBROUTINE ywxpy
884 
886  SUBROUTINE zawpxy(numDim,numPoints,bufferSize,bufferInterval,a,W,X,Y,Z)
887 
888  INTEGER(KIND=4), INTENT(IN) :: numDim
889  INTEGER(KIND=8), INTENT(IN) :: numPoints,bufferSize(numdim)
890  INTEGER(KIND=8), INTENT(IN) :: bufferInterval(2*numdim)
891  REAL(KIND=8), INTENT(IN) :: a
892  REAL(KIND=8), INTENT(IN) :: W(numpoints)
893  REAL(KIND=8), INTENT(IN) :: X(numpoints)
894  REAL(KIND=8), INTENT(IN) :: Y(numpoints)
895  REAL(KIND=8), INTENT(OUT) :: Z(numpoints)
896 
897  INTEGER(KIND=8) :: I, J, K
898  INTEGER(KIND=8) :: nPlane, zIndex, yIndex, yzIndex, bufferIndex, xSize
899  INTEGER(KIND=8) :: iStart,iEnd,jStart,jEnd,kStart,kEnd
900 
901  istart = bufferinterval(1)
902  iend = bufferinterval(2)
903  xsize = buffersize(1)
904 
905  IF(numdim == 1) THEN
906  DO i = istart, iend
907  z(i) = a*w(i) + x(i)*y(i)
908  END DO
909  ELSE IF(numdim == 2) THEN
910  jstart = bufferinterval(3)
911  jend = bufferinterval(4)
912  DO j = jstart, jend
913  yindex = (j-1)*xsize
914  DO i = istart, iend
915  bufferindex = yindex + i
916  z(bufferindex) = a*w(bufferindex) + x(bufferindex)*y(bufferindex)
917  END DO
918  END DO
919  ELSE IF(numdim == 3) THEN
920  jstart = bufferinterval(3)
921  jend = bufferinterval(4)
922  kstart = bufferinterval(5)
923  kend = bufferinterval(6)
924  nplane = xsize * buffersize(2)
925  DO k = kstart, kend
926  zindex = (k-1)*nplane
927  DO j = jstart, jend
928  yzindex = zindex + (j-1)*xsize
929  DO i = istart, iend
930  bufferindex = yzindex + i
931  z(bufferindex) = a*w(bufferindex) + x(bufferindex)*y(bufferindex)
932  END DO
933  END DO
934  END DO
935  ENDIF
936 
937  END SUBROUTINE zawpxy
938 
940  SUBROUTINE zvwpxy(numDim,numPoints,bufferSize,bufferInterval,V,W,X,Y,Z)
941 
942  INTEGER(KIND=4), INTENT(IN) :: numDim
943  INTEGER(KIND=8), INTENT(IN) :: numPoints,bufferSize(numdim)
944  INTEGER(KIND=8), INTENT(IN) :: bufferInterval(2*numdim)
945  REAL(KIND=8), INTENT(IN) :: V(numpoints)
946  REAL(KIND=8), INTENT(IN) :: W(numpoints)
947  REAL(KIND=8), INTENT(IN) :: X(numpoints)
948  REAL(KIND=8), INTENT(IN) :: Y(numpoints)
949  REAL(KIND=8), INTENT(OUT) :: Z(numpoints)
950 
951  INTEGER(KIND=8) :: I, J, K
952  INTEGER(KIND=8) :: nPlane, zIndex, yIndex, yzIndex, bufferIndex, xSize
953  INTEGER(KIND=8) :: iStart,iEnd,jStart,jEnd,kStart,kEnd
954 
955  istart = bufferinterval(1)
956  iend = bufferinterval(2)
957  xsize = buffersize(1)
958 
959  IF(numdim == 1) THEN
960  DO i = istart, iend
961  z(i) = v(i)*w(i) + x(i)*y(i)
962  END DO
963  ELSE IF(numdim == 2) THEN
964  jstart = bufferinterval(3)
965  jend = bufferinterval(4)
966  DO j = jstart, jend
967  yindex = (j-1)*xsize
968  DO i = istart, iend
969  bufferindex = yindex + i
970  z(bufferindex) = v(bufferindex)*w(bufferindex) + &
971  x(bufferindex)*y(bufferindex)
972  END DO
973  END DO
974  ELSE IF(numdim == 3) THEN
975  jstart = bufferinterval(3)
976  jend = bufferinterval(4)
977  kstart = bufferinterval(5)
978  kend = bufferinterval(6)
979  nplane = xsize * buffersize(2)
980  DO k = kstart, kend
981  zindex = (k-1)*nplane
982  DO j = jstart, jend
983  yzindex = zindex + (j-1)*xsize
984  DO i = istart, iend
985  bufferindex = yzindex + i
986  z(bufferindex) = v(bufferindex)*w(bufferindex) + &
987  x(bufferindex)*y(bufferindex)
988  END DO
989  END DO
990  END DO
991  ENDIF
992 
993  END SUBROUTINE zvwpxy
994 
996  SUBROUTINE zwmxpy(numDim,numPoints,bufferSize,bufferInterval,W,X,Y,Z)
997 
998  INTEGER(KIND=4), INTENT(IN) :: numDim
999  INTEGER(KIND=8), INTENT(IN) :: numPoints,bufferSize(numdim)
1000  INTEGER(KIND=8), INTENT(IN) :: bufferInterval(2*numdim)
1001  REAL(KIND=8), INTENT(IN) :: W(numpoints)
1002  REAL(KIND=8), INTENT(IN) :: X(numpoints)
1003  REAL(KIND=8), INTENT(IN) :: Y(numpoints)
1004  REAL(KIND=8), INTENT(OUT) :: Z(numpoints)
1005 
1006  INTEGER(KIND=8) :: I, J, K
1007  INTEGER(KIND=8) :: nPlane, zIndex, yIndex, yzIndex, bufferIndex, xSize
1008  INTEGER(KIND=8) :: iStart,iEnd,jStart,jEnd,kStart,kEnd
1009 
1010  istart = bufferinterval(1)
1011  iend = bufferinterval(2)
1012  xsize = buffersize(1)
1013 
1014  IF(numdim == 1) THEN
1015  DO i = istart, iend
1016  z(i) = w(i)*(x(i)+y(i))
1017  END DO
1018  ELSE IF(numdim == 2) THEN
1019  jstart = bufferinterval(3)
1020  jend = bufferinterval(4)
1021  DO j = jstart, jend
1022  yindex = (j-1)*xsize
1023  DO i = istart, iend
1024  bufferindex = yindex + i
1025  z(bufferindex) = w(bufferindex)*(x(bufferindex)+y(bufferindex))
1026  END DO
1027  END DO
1028  ELSE IF(numdim == 3) THEN
1029  jstart = bufferinterval(3)
1030  jend = bufferinterval(4)
1031  kstart = bufferinterval(5)
1032  kend = bufferinterval(6)
1033  nplane = xsize * buffersize(2)
1034  DO k = kstart, kend
1035  zindex = (k-1)*nplane
1036  DO j = jstart, jend
1037  yzindex = zindex + (j-1)*xsize
1038  DO i = istart, iend
1039  bufferindex = yzindex + i
1040  z(bufferindex) = w(bufferindex)*(x(bufferindex) + y(bufferindex))
1041  END DO
1042  END DO
1043  END DO
1044  ENDIF
1045 
1046  END SUBROUTINE zwmxpy
1047 
1049  SUBROUTINE zxdoty(numDim,numPoints,bufferSize,bufferInterval,numComponents,X,Y,Z)
1050 
1051  INTEGER(KIND=4), INTENT(IN) :: numDim, numComponents
1052  INTEGER(KIND=8), INTENT(IN) :: numPoints,bufferSize(numdim)
1053  INTEGER(KIND=8), INTENT(IN) :: bufferInterval(2*numdim)
1054  REAL(KIND=8), INTENT(IN) :: X(numpoints*numcomponents)
1055  REAL(KIND=8), INTENT(IN) :: Y(numpoints*numcomponents)
1056  REAL(KIND=8), INTENT(OUT) :: Z(numpoints)
1057 
1058  INTEGER(KIND=8) :: I, J, K
1059  INTEGER(KIND=4) :: iComp
1060  INTEGER(KIND=8) :: nPlane, zIndex, yIndex, yzIndex, bufferIndex, xSize
1061  INTEGER(KIND=8) :: iStart,iEnd,jStart,jEnd,kStart,kEnd, compOffset
1062 
1063  REAL(KIND=8) :: zero = 0.0_8
1064 
1065  istart = bufferinterval(1)
1066  iend = bufferinterval(2)
1067  xsize = buffersize(1)
1068 
1069  CALL assignmentxa(numdim,numpoints,buffersize,bufferinterval, &
1070  zero,z)
1071 
1072  IF(numdim == 1) THEN
1073  DO icomp = 0, numcomponents-1
1074  compoffset = icomp*numpoints
1075  DO i = istart, iend
1076  z(i) = z(i) + x(i+compoffset)*y(i+compoffset)
1077  END DO
1078  END DO
1079  ELSE IF(numdim == 2) THEN
1080  jstart = bufferinterval(3)
1081  jend = bufferinterval(4)
1082  DO icomp = 0,numcomponents-1
1083  compoffset = icomp*numpoints
1084  DO j = jstart, jend
1085  yindex = (j-1)*xsize
1086  DO i = istart, iend
1087  bufferindex = yindex + i
1088  z(bufferindex) = z(bufferindex) + &
1089  x(bufferindex+compoffset)*y(bufferindex+compoffset)
1090  END DO
1091  END DO
1092  END DO
1093  ELSE IF(numdim == 3) THEN
1094  jstart = bufferinterval(3)
1095  jend = bufferinterval(4)
1096  kstart = bufferinterval(5)
1097  kend = bufferinterval(6)
1098  nplane = xsize * buffersize(2)
1099  DO icomp = 0,numcomponents-1
1100  compoffset = icomp*numpoints
1101  DO k = kstart, kend
1102  zindex = (k-1)*nplane
1103  DO j = jstart, jend
1104  yzindex = zindex + (j-1)*xsize
1105  DO i = istart, iend
1106  bufferindex = yzindex + i
1107  z(bufferindex) = z(bufferindex) + &
1108  x(bufferindex+compoffset)*y(bufferindex+compoffset)
1109  END DO
1110  END DO
1111  END DO
1112  END DO
1113  ENDIF
1114 
1115  END SUBROUTINE zxdoty
1116 
1118  SUBROUTINE xax(numDim,numPoints,bufferSize,bufferInterval,a,X)
1119 
1120  INTEGER(KIND=4), INTENT(IN) :: numDim
1121  INTEGER(KIND=8), INTENT(IN) :: numPoints,bufferSize(numdim)
1122  INTEGER(KIND=8), INTENT(IN) :: bufferInterval(2*numdim)
1123  REAL(KIND=8), INTENT(IN) :: a
1124  REAL(KIND=8), INTENT(INOUT) :: X(numpoints)
1125 
1126  INTEGER(KIND=8) :: I, J, K
1127  INTEGER(KIND=8) :: nPlane, zIndex, yIndex, yzIndex, bufferIndex, xSize
1128  INTEGER(KIND=8) :: iStart,iEnd,jStart,jEnd,kStart,kEnd
1129 
1130  istart = bufferinterval(1)
1131  iend = bufferinterval(2)
1132  xsize = buffersize(1)
1133 
1134  IF(numdim == 1) THEN
1135  DO i = istart, iend
1136  x(i) = a*x(i)
1137  END DO
1138  ELSE IF(numdim == 2) THEN
1139  jstart = bufferinterval(3)
1140  jend = bufferinterval(4)
1141  DO j = jstart, jend
1142  yindex = (j-1)*xsize
1143  DO i = istart, iend
1144  bufferindex = yindex + i
1145  x(bufferindex) = a*x(bufferindex)
1146  END DO
1147  END DO
1148  ELSE IF(numdim == 3) THEN
1149  jstart = bufferinterval(3)
1150  jend = bufferinterval(4)
1151  kstart = bufferinterval(5)
1152  kend = bufferinterval(6)
1153  nplane = xsize * buffersize(2)
1154  DO k = kstart, kend
1155  zindex = (k-1)*nplane
1156  DO j = jstart, jend
1157  yzindex = zindex + (j-1)*xsize
1158  DO i = istart, iend
1159  bufferindex = yzindex + i
1160  x(bufferindex) = a*x(bufferindex)
1161  END DO
1162  END DO
1163  END DO
1164  ENDIF
1165 
1166  END SUBROUTINE xax
1167 
1169  SUBROUTINE assignmentyx(numDim,numPoints,bufferSize,bufferInterval,X,Y)
1170 
1171  INTEGER(KIND=4), INTENT(IN) :: numDim
1172  INTEGER(KIND=8), INTENT(IN) :: numPoints
1173  INTEGER(KIND=8), INTENT(IN) :: bufferSize(numdim)
1174  INTEGER(KIND=8), INTENT(IN) :: bufferInterval(2*numdim)
1175  REAL(KIND=8), INTENT(IN) :: X(numpoints)
1176  REAL(KIND=8), INTENT(INOUT) :: Y(numpoints)
1177 
1178  INTEGER(KIND=8) :: I, J, K
1179  INTEGER(KIND=8) :: nPlane, zIndex, yIndex, yzIndex, bufferIndex, xSize
1180  INTEGER(KIND=8) :: iStart,iEnd,jStart,jEnd,kStart,kEnd
1181 
1182  istart = bufferinterval(1)
1183  iend = bufferinterval(2)
1184  xsize = buffersize(1)
1185 
1186  IF(numdim == 1) THEN
1187  DO i = istart, iend
1188  y(i) = x(i)
1189  END DO
1190  ELSE IF(numdim == 2) THEN
1191  jstart = bufferinterval(3)
1192  jend = bufferinterval(4)
1193  DO j = jstart, jend
1194  yindex = (j-1)*xsize
1195  DO i = istart, iend
1196  bufferindex = yindex + i
1197  y(bufferindex) = x(bufferindex)
1198  END DO
1199  END DO
1200  ELSE IF(numdim == 3) THEN
1201  jstart = bufferinterval(3)
1202  jend = bufferinterval(4)
1203  kstart = bufferinterval(5)
1204  kend = bufferinterval(6)
1205  nplane = xsize * buffersize(2)
1206  DO k = kstart, kend
1207  zindex = (k-1)*nplane
1208  DO j = jstart, jend
1209  yzindex = zindex + (j-1)*xsize
1210  DO i = istart, iend
1211  bufferindex = yzindex + i
1212  y(bufferindex) = x(bufferindex)
1213  END DO
1214  END DO
1215  END DO
1216  ENDIF
1217 
1218  END SUBROUTINE assignmentyx
1219 
1221  SUBROUTINE assignmentxa(numDim,numPoints,bufferSize,bufferInterval,a,X)
1222 
1223  INTEGER(KIND=4), INTENT(IN) :: numDim
1224  INTEGER(KIND=8), INTENT(IN) :: numPoints
1225  INTEGER(KIND=8), INTENT(IN) :: bufferSize(numdim)
1226  INTEGER(KIND=8), INTENT(IN) :: bufferInterval(2*numdim)
1227  REAL(KIND=8), INTENT(IN) :: a
1228  REAL(KIND=8), INTENT(INOUT) :: X(numpoints)
1229 
1230  INTEGER(KIND=8) :: I, J, K
1231  INTEGER(KIND=8) :: nPlane, zIndex, yIndex, yzIndex, bufferIndex, xSize
1232  INTEGER(KIND=8) :: iStart,iEnd,jStart,jEnd,kStart,kEnd
1233 
1234 
1235  istart = bufferinterval(1)
1236  iend = bufferinterval(2)
1237  xsize = buffersize(1)
1238 
1239  IF(numdim == 1) THEN
1240  DO i = istart, iend
1241  x(i) = a
1242  END DO
1243  ELSE IF(numdim == 2) THEN
1244  jstart = bufferinterval(3)
1245  jend = bufferinterval(4)
1246  DO j = jstart, jend
1247  yindex = (j-1)*xsize
1248  DO i = istart, iend
1249  bufferindex = yindex + i
1250  x(bufferindex) = a
1251  END DO
1252  END DO
1253  ELSE IF(numdim == 3) THEN
1254  jstart = bufferinterval(3)
1255  jend = bufferinterval(4)
1256  kstart = bufferinterval(5)
1257  kend = bufferinterval(6)
1258  nplane = xsize * buffersize(2)
1259  DO k = kstart, kend
1260  zindex = (k-1)*nplane
1261  DO j = jstart, jend
1262  yzindex = zindex + (j-1)*xsize
1263  DO i = istart, iend
1264  bufferindex = yzindex + i
1265  x(bufferindex) = a
1266  END DO
1267  END DO
1268  END DO
1269  ENDIF
1270 
1271  END SUBROUTINE assignmentxa
1272 
1274  SUBROUTINE assignmentyabsx(numDim,numPoints,bufferSize,bufferInterval,X,Y)
1275 
1276  INTEGER(KIND=4), INTENT(IN) :: numDim
1277  INTEGER(KIND=8), INTENT(IN) :: numPoints
1278  INTEGER(KIND=8), INTENT(IN) :: bufferSize(numdim)
1279  INTEGER(KIND=8), INTENT(IN) :: bufferInterval(2*numdim)
1280  REAL(KIND=8), INTENT(IN) :: X(numpoints)
1281  REAL(KIND=8), INTENT(OUT) :: Y(numpoints)
1282 
1283  INTEGER(KIND=8) :: I, J, K
1284  INTEGER(KIND=8) :: nPlane, zIndex, yIndex, yzIndex, bufferIndex, xSize
1285  INTEGER(KIND=8) :: iStart,iEnd,jStart,jEnd,kStart,kEnd
1286 
1287 
1288  istart = bufferinterval(1)
1289  iend = bufferinterval(2)
1290  xsize = buffersize(1)
1291 
1292  IF(numdim == 1) THEN
1293  DO i = istart, iend
1294  y(i) = abs(x(i))
1295  END DO
1296  ELSE IF(numdim == 2) THEN
1297  jstart = bufferinterval(3)
1298  jend = bufferinterval(4)
1299  DO j = jstart, jend
1300  yindex = (j-1)*xsize
1301  DO i = istart, iend
1302  bufferindex = yindex + i
1303  y(bufferindex) = abs(x(bufferindex))
1304  END DO
1305  END DO
1306  ELSE IF(numdim == 3) THEN
1307  jstart = bufferinterval(3)
1308  jend = bufferinterval(4)
1309  kstart = bufferinterval(5)
1310  kend = bufferinterval(6)
1311  nplane = xsize * buffersize(2)
1312  DO k = kstart, kend
1313  zindex = (k-1)*nplane
1314  DO j = jstart, jend
1315  yzindex = zindex + (j-1)*xsize
1316  DO i = istart, iend
1317  bufferindex = yzindex + i
1318  y(bufferindex) = abs(x(bufferindex))
1319  END DO
1320  END DO
1321  END DO
1322  ENDIF
1323 
1324  END SUBROUTINE assignmentyabsx
1325 
1327  SUBROUTINE veclen(numDim,numPoints,bufferSize,bufferInterval,numComp,V,lenV)
1328 
1329  INTEGER(KIND=4), INTENT(IN) :: numDim, numComp
1330  INTEGER(KIND=8), INTENT(IN) :: numPoints
1331  INTEGER(KIND=8), INTENT(IN) :: bufferSize(numdim)
1332  INTEGER(KIND=8), INTENT(IN) :: bufferInterval(2*numdim)
1333  REAL(KIND=8), INTENT(IN) :: V(numcomp*numpoints)
1334  REAL(KIND=8), INTENT(OUT) :: lenV(numpoints)
1335 
1336  INTEGER(KIND=8) :: I, J, K, L
1337  INTEGER(KIND=8) :: nPlane, zIndex, yIndex, yzIndex, bufferIndex, xSize
1338  INTEGER(KIND=8) :: iStart,iEnd,jStart,jEnd,kStart,kEnd
1339 
1340 
1341  istart = bufferinterval(1)
1342  iend = bufferinterval(2)
1343  xsize = buffersize(1)
1344 
1345  IF(numdim == 1) THEN
1346  DO i = istart, iend
1347  lenv(i) = 0.0_8
1348  DO l = 1, numcomp
1349  lenv(i) = lenv(i) + v((l-1)*numpoints+i)*v((l-1)*numpoints+i)
1350  END DO
1351  lenv(i) = sqrt(lenv(i))
1352  END DO
1353  ELSE IF(numdim == 2) THEN
1354  jstart = bufferinterval(3)
1355  jend = bufferinterval(4)
1356  DO j = jstart, jend
1357  yindex = (j-1)*xsize
1358  DO i = istart, iend
1359  bufferindex = yindex + i
1360  lenv(bufferindex) = 0.0_8
1361  DO l = 1,numcomp
1362  lenv(bufferindex) = lenv(bufferindex) + &
1363  v((l-1)*numpoints+bufferindex)*v((l-1)*numpoints+bufferindex)
1364  END DO
1365  lenv(bufferindex) = sqrt(lenv(bufferindex))
1366  END DO
1367  END DO
1368  ELSE IF(numdim == 3) THEN
1369  jstart = bufferinterval(3)
1370  jend = bufferinterval(4)
1371  kstart = bufferinterval(5)
1372  kend = bufferinterval(6)
1373  nplane = xsize * buffersize(2)
1374  DO k = kstart, kend
1375  zindex = (k-1)*nplane
1376  DO j = jstart, jend
1377  yzindex = zindex + (j-1)*xsize
1378  DO i = istart, iend
1379  bufferindex = yzindex + i
1380  lenv(bufferindex) = 0.0_8
1381  DO l = 1,numcomp
1382  lenv(bufferindex) = lenv(bufferindex) + &
1383  v((l-1)*numpoints+bufferindex)*v((l-1)*numpoints+bufferindex)
1384  END DO
1385  lenv(bufferindex) = sqrt(lenv(bufferindex))
1386  END DO
1387  END DO
1388  END DO
1389  ENDIF
1390 
1391  END SUBROUTINE veclen
1392 
1393 
1395  SUBROUTINE yaxm1(numDim,numPoints,bufferSize,bufferInterval,a,X,Y)
1396 
1397  INTEGER(KIND=4), INTENT(IN) :: numDim
1398  INTEGER(KIND=8), INTENT(IN) :: numPoints,bufferSize(numdim)
1399  INTEGER(KIND=8), INTENT(IN) :: bufferInterval(2*numdim)
1400  REAL(KIND=8), INTENT(IN) :: a
1401  REAL(KIND=8), INTENT(INOUT) :: X(numpoints)
1402  REAL(KIND=8), INTENT(INOUT) :: Y(numpoints)
1403 
1404  INTEGER(KIND=8) :: I, J, K
1405  INTEGER(KIND=8) :: nPlane, zIndex, yIndex, yzIndex, bufferIndex, xSize
1406  INTEGER(KIND=8) :: iStart,iEnd,jStart,jEnd,kStart,kEnd
1407 
1408  istart = bufferinterval(1)
1409  iend = bufferinterval(2)
1410  xsize = buffersize(1)
1411 
1412  IF(numdim == 1) THEN
1413  DO i = istart, iend
1414  y(i) = a/x(i)
1415  END DO
1416  ELSE IF(numdim == 2) THEN
1417  jstart = bufferinterval(3)
1418  jend = bufferinterval(4)
1419  DO j = jstart, jend
1420  yindex = (j-1)*xsize
1421  DO i = istart, iend
1422  bufferindex = yindex + i
1423  y(bufferindex) = a/x(bufferindex)
1424  END DO
1425  END DO
1426  ELSE IF(numdim == 3) THEN
1427  jstart = bufferinterval(3)
1428  jend = bufferinterval(4)
1429  kstart = bufferinterval(5)
1430  kend = bufferinterval(6)
1431  nplane = xsize * buffersize(2)
1432  DO k = kstart, kend
1433  zindex = (k-1)*nplane
1434  DO j = jstart, jend
1435  yzindex = zindex + (j-1)*xsize
1436  DO i = istart, iend
1437  bufferindex = yzindex + i
1438  y(bufferindex) = a/x(bufferindex)
1439  END DO
1440  END DO
1441  END DO
1442  ENDIF
1443 
1444  END SUBROUTINE yaxm1
1445 
1446 
1448  SUBROUTINE determinant3x3(numPoints,bufferSize,bufferInterval,inMatrix,matrixDeterminant)
1449 
1450  INTEGER(KIND=8), INTENT(IN) :: numPoints,bufferSize(3)
1451  INTEGER(KIND=8), INTENT(IN) :: bufferInterval(6)
1452  REAL(KIND=8), INTENT(IN) :: inMatrix(9*numpoints)
1453  REAL(KIND=8), INTENT(INOUT) :: matrixDeterminant(numpoints)
1454 
1455  INTEGER(KIND=8) :: I, J, K
1456  INTEGER(KIND=8) :: nPlane, zIndex, yIndex, yzIndex, bufferIndex, xSize
1457  INTEGER(KIND=8) :: iStart,iEnd,jStart,jEnd,kStart,kEnd
1458 
1459  istart = bufferinterval(1)
1460  iend = bufferinterval(2)
1461  xsize = buffersize(1)
1462  jstart = bufferinterval(3)
1463  jend = bufferinterval(4)
1464  kstart = bufferinterval(5)
1465  kend = bufferinterval(6)
1466  nplane = xsize * buffersize(2)
1467 
1468  DO k = kstart, kend
1469  zindex = (k-1)*nplane
1470  DO j = jstart, jend
1471  yzindex = zindex + (j-1)*xsize
1472  DO i = istart, iend
1473  bufferindex = yzindex + i
1474 
1475  matrixdeterminant(bufferindex) = &
1476  inmatrix(bufferindex)* &
1477  (inmatrix(bufferindex+8*numpoints) * &
1478  inmatrix(bufferindex+4*numpoints)- &
1479  inmatrix(bufferindex+7*numpoints)* &
1480  inmatrix(bufferindex+5*numpoints))- &
1481 
1482  inmatrix(bufferindex+3*numpoints)* &
1483  (inmatrix(bufferindex+8*numpoints)* &
1484  inmatrix(bufferindex+numpoints) - &
1485  inmatrix(bufferindex+7*numpoints) * &
1486  inmatrix(bufferindex+2*numpoints)) + &
1487 
1488  inmatrix(bufferindex+6*numpoints)* &
1489  (inmatrix(bufferindex+5*numpoints)* &
1490  inmatrix(bufferindex+numpoints)- &
1491  inmatrix(bufferindex+4*numpoints)* &
1492  inmatrix(bufferindex+2*numpoints))
1493 
1494  END DO
1495  END DO
1496  END DO
1497  END SUBROUTINE determinant3x3
1498 
1500  SUBROUTINE determinant2x2(numPoints,bufferSize,bufferInterval,inMatrix,matrixDeterminant)
1501 
1502  INTEGER(KIND=8), INTENT(IN) :: numPoints,bufferSize(2)
1503  INTEGER(KIND=8), INTENT(IN) :: bufferInterval(4)
1504  REAL(KIND=8), INTENT(IN) :: inMatrix(4*numpoints)
1505  REAL(KIND=8), INTENT(INOUT) :: matrixDeterminant(numpoints)
1506 
1507  INTEGER(KIND=8) :: I, J, K
1508  INTEGER(KIND=8) :: nPlane, zIndex, yIndex, yzIndex, bufferIndex, xSize
1509  INTEGER(KIND=8) :: iStart,iEnd,jStart,jEnd,kStart,kEnd
1510 
1511  istart = bufferinterval(1)
1512  iend = bufferinterval(2)
1513  xsize = buffersize(1)
1514  jstart = bufferinterval(3)
1515  jend = bufferinterval(4)
1516  nplane = xsize * buffersize(2)
1517 
1518  DO j = jstart, jend
1519  yindex = (j-1)*xsize
1520  DO i = istart, iend
1521  bufferindex = yindex + i
1522 
1523  matrixdeterminant(bufferindex) = &
1524  (inmatrix(bufferindex)* &
1525  inmatrix(bufferindex+3*numpoints)) - &
1526  (inmatrix(bufferindex+numpoints)* &
1527  inmatrix(bufferindex+2*numpoints))
1528  END DO
1529  END DO
1530 
1531  END SUBROUTINE determinant2x2
1532 
1534  SUBROUTINE metricsum4(numDim,numPoints,bufferSize,bufferInterval,&
1535  buf1,buf2,buf3,buf4,buf5,buf6,buf7,metricSum)
1536 
1537  INTEGER(KIND=4), INTENT(IN) :: numDim
1538  INTEGER(KIND=8), INTENT(IN) :: numPoints,bufferSize(numdim)
1539  INTEGER(KIND=8), INTENT(IN) :: bufferInterval(2*numdim)
1540  REAL(KIND=8), INTENT(IN) :: buf1(numpoints)
1541  REAL(KIND=8), INTENT(IN) :: buf2(numpoints)
1542  REAL(KIND=8), INTENT(IN) :: buf3(numpoints)
1543  REAL(KIND=8), INTENT(IN) :: buf4(numpoints)
1544  REAL(KIND=8), INTENT(IN) :: buf5(numpoints)
1545  REAL(KIND=8), INTENT(IN) :: buf6(numpoints)
1546  REAL(KIND=8), INTENT(IN) :: buf7(numpoints)
1547  REAL(KIND=8), INTENT(INOUT) :: metricSum(numpoints)
1548 
1549  INTEGER(KIND=8) :: I, J, K
1550  INTEGER(KIND=8) :: nPlane, zIndex, yIndex, yzIndex, bufferIndex, xSize
1551  INTEGER(KIND=8) :: iStart,iEnd,jStart,jEnd,kStart,kEnd
1552 
1553  istart = bufferinterval(1)
1554  iend = bufferinterval(2)
1555  xsize = buffersize(1)
1556 
1557  IF(numdim == 1) THEN
1558  DO i = istart, iend
1559  metricsum(i) = buf1(i)*buf4(i) - buf2(i)*buf3(i) + buf7(i)*(buf5(i)-buf6(i))
1560  END DO
1561  ELSE IF(numdim == 2) THEN
1562  jstart = bufferinterval(3)
1563  jend = bufferinterval(4)
1564  DO j = jstart, jend
1565  yindex = (j-1)*xsize
1566  DO i = istart, iend
1567  bufferindex = yindex + i
1568  metricsum(bufferindex) = buf1(bufferindex)*buf4(bufferindex) - &
1569  buf2(bufferindex)*buf3(bufferindex) + &
1570  buf7(bufferindex)*(buf5(bufferindex)-buf6(bufferindex))
1571  END DO
1572  END DO
1573  ELSE IF(numdim == 3) THEN
1574  jstart = bufferinterval(3)
1575  jend = bufferinterval(4)
1576  kstart = bufferinterval(5)
1577  kend = bufferinterval(6)
1578  nplane = xsize * buffersize(2)
1579  DO k = kstart, kend
1580  zindex = (k-1)*nplane
1581  DO j = jstart, jend
1582  yzindex = zindex + (j-1)*xsize
1583  DO i = istart, iend
1584  bufferindex = yzindex + i
1585  metricsum(bufferindex) = buf1(bufferindex)*buf4(bufferindex) - &
1586  buf2(bufferindex)*buf3(bufferindex) + &
1587  buf7(bufferindex)*(buf5(bufferindex)-buf6(bufferindex))
1588  END DO
1589  END DO
1590  END DO
1591  ENDIF
1592  END SUBROUTINE metricsum4
1593 
1594 END MODULE operators
operators::applyoperatorv
subroutine applyoperatorv(numDim, dimSizes, numComponents, numPoints, opDir, opInterval, numStencils, stencilSizes, stencilStarts, numValues, stencilWeights, stencilOffsets, stencilID, U, dU)
applyoperator applies an operator specified as a stencil set to the provided state data ...
Definition: Operators.f90:174
operators::assignmentyx
subroutine assignmentyx(numDim, numPoints, bufferSize, bufferInterval, X, Y)
ASSIGNMENTYX point-wise operator performing Y = X.
Definition: Operators.f90:1170
operators::ywxpy
subroutine ywxpy(numDim, numPoints, bufferSize, bufferInterval, W, X, Y)
YWXPY point-wise operator performing Y = WX + Y, where all are vectors.
Definition: Operators.f90:835
operators::assignmentyabsx
subroutine assignmentyabsx(numDim, numPoints, bufferSize, bufferInterval, X, Y)
ASSIGNMENTYABSX point-wise operator performing X = scalar a.
Definition: Operators.f90:1275
operators::applyoperator
subroutine applyoperator(numDim, dimSizes, numComponents, numPoints, opDir, opInterval, numStencils, stencilSizes, stencilStarts, numValues, stencilWeights, stencilOffsets, stencilID, U, dU)
applyoperator applies an operator specified as a stencil set to the provided state data ...
Definition: Operators.f90:36
operators::determinant3x3
subroutine determinant3x3(numPoints, bufferSize, bufferInterval, inMatrix, matrixDeterminant)
Computes determinant of 3x3 matrix.
Definition: Operators.f90:1449
operators::yaxpy
subroutine yaxpy(numDim, numPoints, bufferSize, bufferInterval, a, X, Y)
YAXPY point-wise operator performing Y = aX + Y (scalar a)
Definition: Operators.f90:424
operators::zwxpy
subroutine zwxpy(numDim, numPoints, bufferSize, bufferInterval, W, X, Y, Z)
ZWXPY point-wise operator performing Z = WX + Y, where all are vectors.
Definition: Operators.f90:782
operators::yxy
subroutine yxy(numDim, numPoints, bufferSize, bufferInterval, X, Y)
YXY point-wise operator performing Y = XY (all vectors)
Definition: Operators.f90:731
operators::veclen
subroutine veclen(numDim, numPoints, bufferSize, bufferInterval, numComp, V, lenV)
VECLEN point-wise operator returning the length of a numComp-dimensional vector.
Definition: Operators.f90:1328
operators::assignmentxa
subroutine assignmentxa(numDim, numPoints, bufferSize, bufferInterval, a, X)
ASSIGNMENTXA point-wise operator performing X = scalar a.
Definition: Operators.f90:1222
operators::applyoperatorblobs
subroutine applyoperatorblobs(numDim, dimSizes, numComponents, numPointsBuffer, opDir, numStencils, stencilSizes, stencilStarts, numStencilValues, stencilWeights, stencilOffsets, numPointsStencil, numPointsApply, stencilPoints, U, dU)
applyoperatorblobs applies an operator by applying each stencil in turn to all the points to which it...
Definition: Operators.f90:312
operators::yaxpby
subroutine yaxpby(numDim, numPoints, bufferSize, bufferInterval, a, b, X, Y)
YAXPBY point-wise operator performing Y = aX + bY (scalar a,b)
Definition: Operators.f90:577
operators::zxdoty
subroutine zxdoty(numDim, numPoints, bufferSize, bufferInterval, numComponents, X, Y, Z)
ZXDOTY numComponents-vector inner product Z = X * Y.
Definition: Operators.f90:1050
operators::zawpxy
subroutine zawpxy(numDim, numPoints, bufferSize, bufferInterval, a, W, X, Y, Z)
ZAWPXY point-wise operator performing Z = aW + XY.
Definition: Operators.f90:887
operators::applysinglestencil
subroutine applysinglestencil(numDim, dimSizes, numComponents, numPointsBuffer, opDir, numPointsApply, applyPoints, stencilSize, stencilWeights, stencilOffsets, U, dU)
applysinglestencil applies an operator by applying a given stencil to the specified points ...
Definition: Operators.f90:382
operators::determinant2x2
subroutine determinant2x2(numPoints, bufferSize, bufferInterval, inMatrix, matrixDeterminant)
Computes determinant of 2x2 matrix.
Definition: Operators.f90:1501
operators::zwmxpy
subroutine zwmxpy(numDim, numPoints, bufferSize, bufferInterval, W, X, Y, Z)
ZWMXPY point-wise operator performing Z = W(X+Y) where all are vectors.
Definition: Operators.f90:997
operators::yax
subroutine yax(numDim, numPoints, bufferSize, bufferInterval, a, X, Y)
YAX point-wise operator performing Y = aX (scalar a)
Definition: Operators.f90:628
operators::zxy
subroutine zxy(numDim, numPoints, bufferSize, bufferInterval, X, Y, Z)
ZXY point-wise operator performing Z = XY (all vectors)
Definition: Operators.f90:679
operators::zvwpxy
subroutine zvwpxy(numDim, numPoints, bufferSize, bufferInterval, V, W, X, Y, Z)
ZVWPXY point-wise operator performing Z = VW + XY.
Definition: Operators.f90:941
operators::zaxpy
subroutine zaxpy(numDim, numPoints, bufferSize, bufferInterval, a, X, Y, Z)
ZAXPY point-wise operator performing Z = aX + Y (scalar a)
Definition: Operators.f90:475
operators::yaxm1
subroutine yaxm1(numDim, numPoints, bufferSize, bufferInterval, a, X, Y)
YAXM1 point-wise operator performing Y = a/X (scalar a)
Definition: Operators.f90:1396
operators::metricsum4
subroutine metricsum4(numDim, numPoints, bufferSize, bufferInterval, buf1, buf2, buf3, buf4, buf5, buf6, buf7, metricSum)
Computes buf1*buf4 - buf2*buf3 + buf7*(buf5 - buf6)
Definition: Operators.f90:1536
operators::zaxpby
subroutine zaxpby(numDim, numPoints, bufferSize, bufferInterval, a, b, X, Y, Z)
ZAXPBY point-wise operator performing Z = aX + bY (scalar a,b)
Definition: Operators.f90:526
operators::xax
subroutine xax(numDim, numPoints, bufferSize, bufferInterval, a, X)
XAX point-wise operator performing X = aX (scalar a)
Definition: Operators.f90:1119