57
58
59
65
66
67
68#include "implicit_f.inc"
69
70
71
72#include "param_c.inc"
73#include "sphcom.inc"
74#include "task_c.inc"
75
76
77
78 INTEGER NSN, NSNR,NSP2SORTF,NSP2SORTL
79
80 INTEGER KXSP(NISP,*), IXSP(KVOISPH,*), NOD2SP(*),
81 . MA(*), JVOIS(*), JSTOR(*), JPERM(*), IREDUCE,,
82 . KREDUCE(*),LGAUGE(3,*)
83
85 . x(3,*),spbuf(nspbuf,*),dvois(*), bminma(12), gauge(llgauge,*)
86
87
88
89 INTEGER :: I, N, INOD,JNOD,NVOIS,NVOIS1,NVOIS2,K,M,NN,NS,MS,JK,IERROR,JTASK,MY_ADRV,NVOISS
90 my_real :: marge, aaa, aaa2, xi,yi,zi,xj,yj,zj, d2,d1x,d1y,d1z,dk
91 INTEGER NBX,NBY,NBZ
92 INTEGER (KIND=8) :: NBX8,NBY8,NBZ8,RES8,LVOXEL8
93
94
95
96 aaa = sqrt(nsn /
97 . ((bminma(7)-bminma(10))*(bminma(8)-bminma(11))
98 . +(bminma(8)-bminma(11))*(bminma(9)-bminma(12))
99 . +(bminma(9)-bminma(12))*(bminma(7)-bminma(10))))
100
101 aaa = 0.75*aaa
102
103 nbx = nint(aaa*(bminma(7)-bminma(10)))
104 nby = nint(aaa*(bminma(8)-bminma(11)))
105 nbz = nint(aaa*(bminma(9)-bminma(12)))
109
110 nbx8=nbx
111 nby8=nby
112 nbz8=nbz
113 res8=(nbx8+2)*(nby8+2)*(nbz8+2)
115
116 IF(res8 > lvoxel8) THEN
118 aaa = aaa/((nbx8+2)*(nby8+2)*(nbz8+2))
119 aaa = aaa**(third)
120 nbx = int((nbx+2)*aaa)-2
121 nby = int((nby+2)*aaa)-2
122 nbz = int((nbz+2)*aaa)-2
126 ENDIF
127
128 nbx8=nbx
129 nby8=nby
130 nbz8=nbz
131 res8=(nbx8+2)*(nby8+2)*(nbz8+2)
132
133 IF(res8 > lvoxel8) THEN
134 nbx =
min(100,
max(nbx8,1))
135 nby =
min(100,
max(nby8,1))
136 nbz =
min(100,
max(nbz8,1))
137 ENDIF
138
139
140
141 DO i=
inivoxel,(nbx+2)*(nby+2)*(nbz+2)
143 ENDDO
145
146 IF(itask==0)THEN
147 ALLOCATE(
kxspr(nsnr),
ixspr(kvoisph,nsnr),stat=ierror)
148 IF(ierror/=0)THEN
149 CALL ancmsg(msgid=248,anmode=aninfo_blind)
151 END IF
153 END IF
154
156
157
158
160 1 nsn ,nsnr ,x ,bminma ,nod2sp ,
161 2 nbx ,nby ,nbz ,marge ,itask ,
162 3 ma ,spbuf ,jvois ,jstor ,jperm ,
163 4 dvois ,ireduce,nsp2sortf,nsp2sortl,
voxel1 ,
164 5 kxsp ,ixsp ,kreduce ,lgauge ,gauge ,
166
167
168
169
170
171 IF(itask==0)THEN
172
173 ALLOCATE(
nv(numsph*nthread),
iv(numsph*nthread),
ivs(numsph),
174 .
kv(numsph),
iaux(kvoisph*(numsph+nsnr)),
175 . stat=ierror)
176
177 IF(ierror/=0)THEN
178 CALL ancmsg(msgid=248,anmode=aninfo_blind)
180 END IF
181 END IF
182
184
185 nv(itask*numsph+1:(itask+1)*numsph)=0
186
187 DO ns=nsp2sortf,nsp2sortl
188 n =ma(ns)
190 END DO
191
192
194
195
196 my_adrv=itask*numsph
197 DO ns=nsp2sortf,nsp2sortl
198 n =ma(ns)
199 inod=kxsp(3,n)
200 nvois=kxsp(5,n)
201 DO k=1,nvois
202 jnod=ixsp(k,n)
203 IF(jnod > 0) THEN
204 m=nod2sp(jnod)
205 IF(spbuf(1,m) < spbuf(1,n) .OR.
206 . (spbuf(1,m) == spbuf(1,n) .AND. kxsp(8,m) < kxsp(8,n)) )THEN
207
209 nv(my_adrv+ms)=
nv(my_adrv+ms)+1
210
211 END IF
212 ELSE
213
214 ENDIF
215 END DO
216 END DO
217
218 DO n=itask+1,nsnr,nthread
220 DO k=1,nvois
222 IF(jnod > 0) THEN
223 m=nod2sp(jnod)
224 ELSE
225 print *,'internal error'
226 END IF
227
229 nv(my_adrv+ms)=
nv(my_adrv+ms)+1
230
231 END DO
232 END DO
233
234
236
237
238 IF(itask==0)THEN
239
241 DO ns=1,nsp2sort-1
243 DO jtask=1,nthread
244 iv(ns+1)=
iv(ns+1)+
nv((jtask-1)*numsph+ns)
245 END DO
246 END DO
247
248 END IF
249
250
252
253
254 DO ns=nsp2sortf,nsp2sortl
256 END DO
257
258 DO jtask=1,nthread-1
259 DO ns=nsp2sortf,nsp2sortl
260 iv(jtask*numsph+ns)=
iv((jtask-1)*numsph+ns)+
nv((jtask-1)*numsph+ns)
261 END DO
262 END DO
263
264
266
267
268 DO ns=nsp2sortf,nsp2sortl
269 n =ma(ns)
270 inod=kxsp(3,n)
271 nvois=kxsp(5,n)
272 DO k=1,nvois
273 jnod=ixsp(k,n)
274 IF(jnod > 0) THEN
275 m=nod2sp(jnod)
276 IF(spbuf(1,m) < spbuf(1,n) .OR.
277 . (spbuf(1,m) == spbuf(1,n) .AND. kxsp(8,m) < kxsp(8,n)) )THEN
278
280 iaux(
iv(itask*numsph+ms))=inod
281 iv(itask*numsph+ms)=
iv(itask*numsph+ms)+1
282
283 END IF
284 ELSE
285
286 ENDIF
287 END DO
288 END DO
289
290 DO n=itask+1,nsnr,nthread
292 DO k=1,nvois
294 IF(jnod > 0) THEN
295 m=nod2sp(jnod)
296 ELSE
297 print *,'internal error'
298 END IF
299
301 iaux(
iv(itask*numsph+ms))=-n
302 iv(itask*numsph+ms)=
iv(itask*numsph+ms)+1
303
304 END DO
305 END DO
306
307
309
310
311 DO ns=nsp2sortf,nsp2sortl
312 n =ma(ns)
313
314 inod=kxsp(3,n)
315 xi = x(1,inod)
316 yi = x(2,inod)
317 zi = x(3,inod)
318
319 nvois=kxsp(5,n)
320 DO k=1,nvois
321 jnod=ixsp(k,n)
322 IF(jnod > 0) THEN
323 m=nod2sp(jnod)
324 xj = x(1,jnod)
325 yj = x(2,jnod)
326 zj = x(3,jnod)
327 aaa = spbuf(1,n)+spbuf(1,m)
328 ELSE
329 nn=-jnod
330 xj = xsphr(3,nn)
331 yj = xsphr(4,nn)
332 zj = xsphr(5,nn)
333 aaa = spbuf(1,n)+xsphr(2,nn)
334 ENDIF
335
336 aaa2 = aaa*aaa
337
338
339 d1x = xi - xj
340 d1y = yi - yj
341 d1z = zi - zj
342 d2 = d1x*d1x+d1y*d1y+d1z*d1z
343
344 jvois(k)=jnod
345 dvois(k)=d2/aaa2
346
347 END DO
348
349 nvoiss=0
350 DO jtask=1,nthread
351 nvoiss=nvoiss+
nv((jtask-1)*numsph+ns)
352 END DO
353 nvois=nvois+nvoiss
354 IF(nvois>kvoisph)THEN
355 ireduce=1
356 kreduce(n)=1
357 END IF
358
359 DO k=1,nvoiss
361 IF(jnod > 0) THEN
362 m=nod2sp(jnod)
363 xj = x(1,jnod)
364 yj = x(2,jnod)
365 zj = x(3,jnod)
366 aaa = spbuf(1,n)+spbuf(1,m)
367 ELSE
368 nn=-jnod
369 xj = xsphr(3,nn)
370 yj = xsphr(4,nn)
371 zj = xsphr(5,nn)
372 aaa = spbuf(1,n)+xsphr(2,nn)
373 ENDIF
374
375 aaa2 = aaa*aaa
376
377
378 d1x = xi - xj
379 d1y = yi - yj
380 d1z = zi - zj
381 d2 = d1x*d1x+d1y*d1y+d1z*d1z
382
383 jvois(kxsp(5,n)+k)=jnod
384 dvois(kxsp(5,n)+k)=d2/aaa2
385
386 END DO
387
389 IF(kreduce(n)/=0 .AND. nvois > kvoisph) THEN
391
392 CALL myqsort(nvois,dvois,jperm,ierror)
393 DO k=1,nvois
394 jstor(k)=jvois(k)
395 ENDDO
396 DO k=1,kvoisph
397 jvois(k)=jstor(jperm(k))
398 ENDDO
399
400 dk=dvois(kvoisph)
401
402
403 nvois=0
404 DO k=1,kvoisph
405 IF(dvois(k)<dk)THEN
406 nvois=nvois+1
407 END IF
408 END DO
409 END IF
410
411 nvois=
min(nvois,kvoisph)
412 kxsp(5,n)=nvois
413
414 nvois1=0
415 nvois2=nvois+1
416 DO k=1,nvois
417 jk =jvois(k)
418 dk =dvois(k)
419
420 IF(dk < one)THEN
421 nvois1=nvois1+1
422 ixsp(nvois1,n)=jk
423 ELSE
424 nvois2=nvois2-1
425 ixsp(nvois2,n)=jk
426 END IF
427 ENDDO
428 kxsp(4,n)=nvois1
429 IF(nvois1>lvoisph)ireduce=1
430
431 END DO
432
433
435
437
438 RETURN
subroutine myqsort(n, a, perm, error)
integer *8, dimension(:), allocatable iv
integer, dimension(:), allocatable iaux
integer, dimension(:), allocatable kv
integer, dimension(:,:), allocatable ixspr
integer, dimension(:), allocatable kxspr
integer, dimension(:), allocatable ivs
integer, dimension(:), allocatable nv
logical, dimension(:), allocatable bool_sph_sort
integer, dimension(lvoxel) voxel1
subroutine sptrivox(nsn, x, bminma, nod2sp, nbx, nby, nbz, nlist, spbuf, jvois, jstor, jperm, dvois, ireduce, nsphactf, nsphactl, voxel, kxsp, ixsp, kreduce, ipartsp, sz_intp_dist, max_intp_dist_part, pre_search)
subroutine ancmsg(msgid, msgtype, anmode, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13, i14, i15, i16, i17, i18, i19, i20, r1, r2, r3, r4, r5, r6, r7, r8, r9, c1, c2, c3, c4, c5, c6, c7, c8, c9, prmode)