39
40
41
44
45
46
47#include "implicit_f.inc"
48
49
50
51#include "com01_c.inc"
52#include "com04_c.inc"
53#include "sphcom.inc"
54#include "param_c.inc"
55#include "parit_c.inc"
56#include "task_c.inc"
57
58
59
60 INTEGER KXSP(NISP,*),IXSP(KVOISPH,*),NOD2SP(*),
61 . ISPSYM(NSPCOND,*),WSP2SORT(*), ITASK, IREDUCE, KREDUCE(*),
62 . LGAUGE(3,*),ISORTSP
63
65 . x(3,*),spbuf(nspbuf,*),xspsym(3,*), myspatrue, gauge(llgauge,*)
66
67
68
69 integer
70 . n,inod,jnod,j,nvois,m,ncand,k1,k2,nvois1,nvois2,
71 . nvoiss,nvoiss1,nvoiss2, iaux, ierror,
72 . k, jk, nc, js, ns, nn, nb,jj1,jj2, jj, jjj,
73 . mwa(2*kvoisph),jstor(kvoisph), jperm(kvoisph),
74 . lvoired, ig
76 . dms,dms2,dk,
77 . xi,yi,zi,di,xj,yj,zj,dj,dd,dvois(kvoisph),
78 . dwa(numsph)
79 SAVE lvoired
80 LOGICAL :: SORTING_CONDITION
81
82 lvoired = 0
83 IF(ireduce==0)GO TO 100
84
85
86
87
89
90
91 DO ns=itask+1,nsp2sort,nthread
92 n=wsp2sort(ns)
93 dwa(n)=one
94 nvois1 =kxsp(4,n)
95 nvoiss1=kxsp(6,n)
96 IF(kreduce(n)/=0.OR.nvois1+nvoiss1>lvoisph)THEN
97
98 IF(nvois1+nvoiss1>lvoisph)THEN
99 kreduce(n)=kreduce(n)+10
100 lvoired = 1
101 END IF
102
103 inod=kxsp(3,n)
104 xi=x(1,inod)
105 yi=x(2,inod)
106 zi=x(3,inod)
107 di=spbuf(1,n)
108 nvois=kxsp(5,n)
109 ncand=kxsp(5,n)+kxsp(7,n)
110 DO k=1,nvois
111 jnod = ixsp(k,n)
112 IF(jnod>0)THEN
113 m =nod2sp(jnod)
114 xj=x(1,jnod)
115 yj=x(2,jnod)
116 zj=x(3,jnod)
117 dj=spbuf(1,m)
118 ELSE
119 nn = -jnod
120 xj=xsphr(3,nn)
121 yj=xsphr(4,nn)
122 zj=xsphr(5,nn)
123 dj=xsphr(2,nn)
124 END IF
125 dms =di+dj
126 dms2=dms*dms
127 dvois(k)=(xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj)
128 dvois(k)=dvois(k)/dms2
129 END DO
130 DO k=nvois+1,ncand
131 jk = ixsp(k,n)
132 IF(jk>0)THEN
133 nc=mod(jk,nspcond+1)
134 m=jk/(nspcond+1)
135 js=ispsym(nc,m)
136 dj=spbuf(1,m)
137 ELSE
138 nc=mod(-jk,nspcond+1)
139 m =-jk/(nspcond+1)
141 dj =xsphr(2,m)
142 END IF
143 xj =xspsym(1,js)
144 yj =xspsym(2,js)
145 zj =xspsym(3,js)
146 dms =di+dj
147 dms2=dms*dms
148 dvois(k)=(xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj)
149 dvois(k)=dvois(k)/dms2
150 END DO
151
152 sorting_condition = (.NOT.(
bool_sph_sort(n)).OR.isortsp==0.OR.nvois/=ncand)
153 IF(sorting_condition) THEN
154 CALL myqsort(ncand,dvois,jperm,ierror)
155 ELSE
156 DO k=1,kxsp(4,n)
157 jperm(k) = k
158 ENDDO
159 DO k=1,kxsp(5,n)-kxsp(4,n)+1
160 jperm(kxsp(4,n)+k) = kxsp(5,n)-k+1
161 ENDDO
162 ENDIF
163
164 DO k=1,ncand
165 jstor(k) = ixsp(k,n)
166 END DO
167
168 IF(kreduce(n) >= 10)dwa(n)=sqrt(dvois(lvoisph))
169
170 k1=0
171 k2=0
172 DO k=1,ncand
173 jk=jstor(jperm(k))
174 IF(jperm(k) <= nvois) THEN
175 k1=k1+1
176 ixsp(k1,n) = jk
177 ELSE
178 k2=k2+1
179 ixsp(nvois+k2,n) = jk
180 END IF
181 END DO
182
183 END IF
184 END DO
185
186
187
188
190
191
192 IF(lvoired /= 0)THEN
193
194 DO ns=itask+1,nsp2sort,nthread
195 n=wsp2sort(ns)
196 spbuf(1,n)=
min(spbuf(1,n),dwa(n)*spbuf(1,n))
197 spbuf(8,n)=spbuf(1,n)
198 END DO
199 END IF
200
201 IF(nspmd > 1)THEN
202
203
205
206 IF(itask==0) THEN
207
209
210
212 END IF
213 END IF
214
215
217
218
219 DO ns=itask+1,nsp2sort,nthread
220 n=wsp2sort(ns)
221
222 IF(mod(kreduce(n),10)/=0)THEN
223
224 nvois1 =kxsp(4,n)
225 nvois =kxsp(5,n)
226 nvoiss1=kxsp(6,n)
227 nvoiss =kxsp(7,n)
228 inod=kxsp(3,n)
229 xi=x(1,inod)
230 yi=x(2,inod)
231 zi=x(3,inod)
232 di=spbuf(1,n)
233
234
235 jnod = ixsp(nvois,n)
236 IF(jnod>0)THEN
237 m =nod2sp(jnod)
238 xj=x(1,jnod)
239 yj=x(2,jnod)
240 zj=x(3,jnod)
241 dj=spbuf(1,m)
242 ELSE
243 nn = -jnod
244 xj=xsphr(3,nn)
245 yj=xsphr(4,nn)
246 zj=xsphr(5,nn)
247 dj=xsphr(2,nn)
248 END IF
249 dd =(xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj)
250 dms =di+dj
251 dms2=dms*dms
252 dk=dd/dms2
253 myspatrue=
max(zero,
min(myspatrue,dk-one))
254 END IF
255
256 END DO
257
258 100 CONTINUE
259 IF(nspcond==0) THEN
260 DO n = itask+1,
nsphr,nthread
261
263 END DO
264 ELSE
265
266
267 DO n = itask+1,
nsphr,nthread
269 END DO
270 END IF
271
273
274
275 IF(iparit/=0)THEN
276 DO ns=itask+1,nsp2sort,nthread
277 n=wsp2sort(ns)
278 inod=kxsp(3,n)
279 xi=x(1,inod)
280 yi=x(2,inod)
281 zi=x(3,inod)
282 ncand=kxsp(5,n)
283 nvois1=0
284 nvois2=0
285 DO j=1,ncand
286 jnod=ixsp(j,n)
287 IF(jnod>0)THEN
288 m=nod2sp(jnod)
289 xj=x(1,jnod)
290 yj=x(2,jnod)
291 zj=x(3,jnod)
292 dd =(xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj)
293 dms =spbuf(1,n)+spbuf(1,m)
294 dms2=dms*dms
295 IF (kxsp(2,m)/=0.AND.dd<dms2) THEN
296 nvois1=nvois1+1
297 mwa(nvois1)=jnod
298 ELSE
299 nvois2=nvois2+1
300 mwa(kvoisph+nvois2)=jnod
301 END IF
302 ELSE
303 nn = -jnod
304 xj=xsphr(3,nn)
305 yj=xsphr(4,nn)
306 zj=xsphr(5,nn)
307 dd =(xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj)
308 dms =spbuf(1,n)+xsphr(2,nn)
309 dms2=dms*dms
310 IF (nint(xsphr(13,nn))/=0.AND.dd<dms2) THEN
311 nvois1=nvois1+1
312 mwa(nvois1)=jnod
314 ELSE
315 nvois2=nvois2+1
316 mwa(kvoisph+nvois2)=jnod
317 ENDIF
318 END IF
319 ENDDO
320
321 kxsp(4,n)=nvois1
322 DO j=1,nvois1
323 ixsp(j,n)=mwa(j)
324 ENDDO
325 DO j=1,nvois2
326 ixsp(nvois1+j,n)=mwa(kvoisph+j)
327 ENDDO
328
329
330
331 DO k = 1, nvois1
332 jk = ixsp(k,n)
333 IF(jk>0)THEN
334 dvois(k) = kxsp(8,nod2sp(jk))
335 ELSE
336 dvois(k) = nint(xsphr(6,-jk))
337 END IF
338 END DO
339 CALL myqsort(nvois1,dvois,jperm,ierror)
340 DO k=1,nvois1
341 jstor(k) = ixsp(k,n)
342 END DO
343 DO k=1,nvois1
344 ixsp(k,n) = jstor(jperm(k))
345 END DO
346 ENDDO
347
348
349 DO ns=itask+1,nsp2sort,nthread
350 n=wsp2sort(ns)
351 inod=kxsp(3,n)
352 xi =x(1,inod)
353 yi =x(2,inod)
354 zi =x(3,inod)
355 di =spbuf(1,n)
356 nvois2 =kxsp(5,n)
357 nvoiss =kxsp(7,n)
358 nvoiss1=0
359 nvoiss2=0
360 DO k=nvois2+1,nvois2+nvoiss
361 jk=ixsp(k,n)
362 IF(jk>0)THEN
363 nc=mod(jk,nspcond+1)
364 m =jk/(nspcond+1)
365 js=ispsym(nc,m)
366 dj =spbuf(1,m)
367 xj =xspsym(1,js)
368 yj =xspsym(2,js)
369 zj =xspsym(3,js)
370 dms =di+dj
371 dms2=dms*dms
372 dd=(xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj)
373 IF (kxsp(2,m)/=0.AND.dd<dms2) THEN
374 nvoiss1=nvoiss1+1
375 mwa(nvoiss1)=jk
376 ELSE
377 nvoiss2=nvoiss2+1
378 mwa(kvoisph+nvoiss2)=jk
379 ENDIF
380 ELSE
381 nc=mod(-jk,nspcond+1)
382 m =-jk/(nspcond+1)
384 dj =xsphr(2,m)
385 xj =xspsym(1,js)
386 yj =xspsym(2,js)
387 zj =xspsym(3,js)
388 dms =di+dj
389 dms2=dms*dms
390 dd=(xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj)
391 IF (nint(xsphr(13,m))/=0.AND.dd<dms2) THEN
392 nvoiss1=nvoiss1+1
393 mwa(nvoiss1)=jk
394 ELSE
395 nvoiss2=nvoiss2+1
396 mwa(kvoisph+nvoiss2)=jk
397 ENDIF
398 END IF
399 ENDDO
400 kxsp(6,n)=nvoiss1
401 DO j=1,nvoiss1
402 ixsp(nvois2+j,n)=mwa(j)
403 ENDDO
404 DO j=1,nvoiss2
405 ixsp(nvois2+nvoiss1+j,n)=mwa(kvoisph+j)
406 ENDDO
407
408
409 DO k = 1, nvoiss1
410 jk = ixsp(nvois2+k,n)
411 IF(jk>0)THEN
412 m=jk/(nspcond+1)
413 nc=mod(jk,nspcond+1)
414 dvois(k) = kxsp(8,m)
415
416 mwa(k) = nc
417 ELSE
418 m=-jk/(nspcond+1)
419 nc=mod(-jk,nspcond+1)
420 dvois(k) = xsphr(6,m)
421 mwa(k) = nc
422 END IF
423 END DO
424 CALL myqsort(nvoiss1,dvois,jperm,ierror)
425 DO k=1,nvoiss1
426 jstor(k) = ixsp(nvois2+k,n)
427 END DO
428 DO k=1,nvoiss1
429 ixsp(nvois2+k,n) = jstor(jperm(k))
430 END DO
431 DO k=1,nvoiss1
432 jstor(k) = mwa(k)
433 END DO
434 DO k=1,nvoiss1
435 mwa(k) = jstor(jperm(k))
436 END DO
437 IF(nspcond>1) THEN
438
439 m = nint(dvois(1))
440 nb = 1
441 DO k = 2, nvoiss1
442 IF(nint(dvois(k))/=m) THEN
443 IF(nb>1)THEN
444 jj1 = k-nb
445 jj2 = k-1
446
447 DO jj = jj1, jj2-1
448 DO jjj = jj+1, jj2
449 IF(mwa(jj)>mwa(jjj))THEN
450 iaux = mwa(jj)
451 mwa(jj) = mwa(jjj)
452 mwa(jjj) = iaux
453 iaux = ixsp(nvois2+jj,n)
454 ixsp(nvois2+jj,n) = ixsp(nvois2+jjj,n)
455 ixsp(nvois2+jjj,n) = iaux
456 END IF
457 END DO
458 END DO
459 END IF
460 m = nint(dvois(k))
461 nb = 1
462 ELSE
463 nb = nb + 1
464 END IF
465 END DO
466
467 IF(nb>1)THEN
468 jj1 = nvoiss1-nb+1
469 jj2 = nvoiss1
470
471 DO jj = jj1, jj2-1
472 DO jjj = jj+1, jj2
473 IF(mwa(jj)>mwa(jjj))THEN
474 iaux = mwa(jj)
475 mwa(jj) = mwa(jjj)
476 mwa(jjj) = iaux
477 iaux = ixsp(nvois2+jj,n)
478 ixsp(nvois2+jj,n) = ixsp(nvois2+jjj,n)
479 ixsp(nvois2+jjj,n) = iaux
480 END IF
481 END DO
482 END DO
483 END IF
484 END IF
485
486 ENDDO
487
488 ELSE
489 DO ns=itask+1,nsp2sort,nthread
490 n=wsp2sort(ns)
491 inod=kxsp(3,n)
492 xi=x(1,inod)
493 yi=x(2,inod)
494 zi=x(3,inod)
495 ncand=kxsp(5,n)
496 nvois1=0
497 nvois2=0
498 DO j=1,ncand
499 jnod=ixsp(j,n)
500 IF(jnod>0)THEN
501 m=nod2sp(jnod)
502 xj=x(1,jnod)
503 yj=x(2,jnod)
504 zj=x(3,jnod)
505 dd =(xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj)
506 dms =spbuf(1,n)+spbuf(1,m)
507 dms2=dms*dms
508 IF (kxsp(2,m)/=0.AND.dd<dms2) THEN
509 nvois1=nvois1+1
510 mwa(nvois1)=jnod
511 ELSE
512 nvois2=nvois2+1
513 mwa(kvoisph+nvois2)=jnod
514 END IF
515 ELSE
516 nn = -jnod
517 xj=xsphr(3,nn)
518 yj=xsphr(4,nn)
519 zj=xsphr(5,nn)
520 dd =(xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj)
521 dms =spbuf(1,n)+xsphr(2,nn)
522 dms2=dms*dms
523 IF (nint(xsphr(13,nn))/=0.AND.dd<dms2) THEN
524 nvois1=nvois1+1
525 mwa(nvois1)=jnod
527 ELSE
528 nvois2=nvois2+1
529 mwa(kvoisph+nvois2)=jnod
530 ENDIF
531 END IF
532 ENDDO
533
534 kxsp(4,n)=nvois1
535 DO j=1,nvois1
536 ixsp(j,n)=mwa(j)
537 ENDDO
538 DO j=1,nvois2
539 ixsp(nvois1+j,n)=mwa(kvoisph+j)
540 ENDDO
541 ENDDO
542
543
544 DO ns=itask+1,nsp2sort,nthread
545 n=wsp2sort(ns)
546 inod=kxsp(3,n)
547 xi =x(1,inod)
548 yi =x(2,inod)
549 zi =x(3,inod)
550 di =spbuf(1,n)
551 nvois2 =kxsp(5,n)
552 nvoiss =kxsp(7,n)
553 nvoiss1=0
554 nvoiss2=0
555 DO k=nvois2+1,nvois2+nvoiss
556 jk=ixsp(k,n)
557 IF(jk>0)THEN
558 nc=mod(jk,nspcond+1)
559 m =jk/(nspcond+1)
560 js=ispsym(nc,m)
561 dj =spbuf(1,m)
562 xj =xspsym(1,js)
563 yj =xspsym(2,js)
564 zj =xspsym(3,js)
565 dms =di+dj
566 dms2=dms*dms
567 dd=(xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj)
568 IF (kxsp(2,m)/=0.AND.dd<dms2) THEN
569 nvoiss1=nvoiss1+1
570 mwa(nvoiss1)=jk
571 ELSE
572 nvoiss2=nvoiss2+1
573 mwa(kvoisph+nvoiss2)=jk
574 ENDIF
575 ELSE
576 nc=mod(-jk,nspcond+1)
577 m =-jk/(nspcond+1)
579 dj =xsphr(2,m)
580 xj =xspsym(1,js)
581 yj =xspsym(2,js)
582 zj =xspsym(3,js)
583 dms =di+dj
584 dms2=dms*dms
585 dd=(xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj)
586 IF (nint(xsphr(13,m))/=0.AND.dd<dms2) THEN
587 nvoiss1=nvoiss1+1
588 mwa(nvoiss1)=jk
589 ELSE
590 nvoiss2=nvoiss2+1
591 mwa(kvoisph+nvoiss2)=jk
592 ENDIF
593 END IF
594 ENDDO
595 kxsp(6,n)=nvoiss1
596 DO j=1,nvoiss1
597 ixsp(nvois2+j,n)=mwa(j)
598 ENDDO
599 DO j=1,nvoiss2
600 ixsp(nvois2+nvoiss1+j,n)=mwa(kvoisph+j)
601 ENDDO
602 ENDDO
603 END IF
604
605
606 DO ig=1,nbgauge
607 IF(lgauge(1,ig) > -(numels+1))cycle
608 n=numsph+ig
609 xi =gauge(2,ig)
610 yi =gauge(3,ig)
611 zi =gauge(4,ig)
612 ncand=kxsp(5,n)
613 nvois1=0
614 nvois2=0
615 DO j=1,ncand
616 jnod=ixsp(j,n)
617 IF(jnod>0)THEN
618 m=nod2sp(jnod)
619 xj=x(1,jnod)
620 yj=x(2,jnod)
621 zj=x(3,jnod)
622 dd =(xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj)
623 dms =two*spbuf(1,m)
624 dms2=dms*dms
625 IF (kxsp(2,m)/=0.AND.dd<dms2) THEN
626 nvois1=nvois1+1
627 mwa(nvois1)=jnod
628 ELSE
629 nvois2=nvois2+1
630 mwa(kvoisph+nvois2)=jnod
631 END IF
632 ELSE
633 nn = -jnod
634 xj=xsphr(3,nn)
635 yj=xsphr(4,nn)
636 zj=xsphr(5,nn)
637 dd =(xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj)
638 dms =two*xsphr(2,nn)
639 dms2=dms*dms
640 IF (nint(xsphr(13,nn))/=0.AND.dd<dms2) THEN
641 nvois1=nvois1+1
642 mwa(nvois1)=jnod
644 ELSE
645 nvois2=nvois2+1
646 mwa(kvoisph+nvois2)=jnod
647 ENDIF
648 END IF
649 ENDDO
650
651 kxsp(4,n)=nvois1
652 DO j=1,nvois1
653 ixsp(j,n)=mwa(j)
654 ENDDO
655 DO j=1,nvois2
656 ixsp(nvois1+j,n)=mwa(kvoisph+j)
657 ENDDO
658
659
660 nvois2 =kxsp(5,n)
661 nvoiss =kxsp(7,n)
662 nvoiss1=0
663 nvoiss2=0
664 DO k=nvois2+1,nvois2+nvoiss
665 jk=ixsp(k,n)
666 IF(jk>0)THEN
667 nc=mod(jk,nspcond+1)
668 m =jk/(nspcond+1)
669 js=ispsym(nc,m)
670 dj =spbuf(1,m)
671 xj =xspsym(1,js)
672 yj =xspsym(2,js)
673 zj =xspsym(3,js)
674 dms =two*dj
675 dms2=dms*dms
676 dd=(xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj)
677 IF (kxsp(2,m)/=0.AND.dd<dms2) THEN
678 nvoiss1=nvoiss1+1
679 mwa(nvoiss1)=jk
680 ELSE
681 nvoiss2=nvoiss2+1
682 mwa(kvoisph+nvoiss2)=jk
683 ENDIF
684 ELSE
685 nc=mod(-jk,nspcond+1)
686 m =-jk/(nspcond+1)
688 dj =xsphr(2,m)
689 xj =xspsym(1,js)
690 yj =xspsym(2,js)
691 zj =xspsym(3,js)
692 dms =two*dj
693 dms2=dms*dms
694 dd=(xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj)
695 IF (nint(xsphr(13,m))/=0.AND.dd<dms2) THEN
696 nvoiss1=nvoiss1+1
697 mwa(nvoiss1)=jk
698 ELSE
699 nvoiss2=nvoiss2+1
700 mwa(kvoisph+nvoiss2)=jk
701 ENDIF
702 END IF
703 ENDDO
704 kxsp(6,n)=nvoiss1
705 DO j=1,nvoiss1
706 ixsp(nvois2+j,n)=mwa(j)
707 ENDDO
708 DO j=1,nvoiss2
709 ixsp(nvois2+nvoiss1+j,n)=mwa(kvoisph+j)
710 ENDDO
711 ENDDO
712
713
714 RETURN
subroutine myqsort(n, a, perm, error)
logical, dimension(:), allocatable bool_sph_sort
integer, dimension(:), allocatable isphr
integer, dimension(:,:), allocatable ispsymr
subroutine spmd_allglob_isum9(v, len)
subroutine spmd_sphgeth(kxsp, spbuf)