V3FIT
parallel_vmec_module.f90
1 !------------------------------------------------
2 MODULE parallel_vmec_module
3  !------------------------------------------------
4 
5  USE stel_kinds
6  USE mpi_inc
7  IMPLICIT NONE
8  INTEGER :: TOFU
9  INTEGER :: grank=0, gnranks=1
10  INTEGER :: vrank=0, vnranks=1
11  INTEGER :: rank=0, nranks=1
12  INTEGER :: last_ns = -1
13  INTEGER :: par_ns, mynsnum, par_nzeta, par_ntheta3
14  INTEGER :: par_ntmax, par_ntor, par_mpol1, par_nznt, par_nuv, par_nuv3
15  INTEGER :: blocksize, ntmaxblocksize
16  INTEGER :: tlglob, trglob, t1lglob, t1rglob
17  INTEGER :: nuvmin, nuvmax, nuv3min, nuv3max
18  INTEGER :: MPI_ERR
19 #if defined(MPI_OPT)
20  INTEGER :: MPI_STAT(MPI_STATUS_SIZE)
21 #endif
22  INTEGER :: RUNVMEC_COMM_WORLD
23  INTEGER :: NS_COMM
24  INTEGER :: VAC_COMM
25  INTEGER :: TWODCOMM, px, py
26  INTEGER :: NS_RESLTN=0
27  INTEGER :: num_grids
28  INTEGER, ALLOCATABLE, DIMENSION(:) :: grid_procs
29  INTEGER, ALLOCATABLE, DIMENSION(:) :: grid_size
30  REAL(dp), ALLOCATABLE, DIMENSION(:) :: grid_time
31  REAL(dp), ALLOCATABLE, DIMENSION(:) :: f3d_time
32  INTEGER, ALLOCATABLE, DIMENSION(:) :: f3d_num
33 
34  REAL(dp), ALLOCATABLE, DIMENSION(:) :: vgrid_time
35 
36  INTEGER, ALLOCATABLE, DIMENSION(:) :: tlglob_arr, trglob_arr
37  INTEGER, ALLOCATABLE, DIMENSION(:) :: nuv3min_arr,nuv3max_arr
38  INTEGER, ALLOCATABLE, DIMENSION(:) :: trow, brow
39  INTEGER, ALLOCATABLE, DIMENSION(:) :: lcol, rcol
40 
41  INTEGER, ALLOCATABLE, DIMENSION (:) :: blkrcounts, blkdisp
42  INTEGER, ALLOCATABLE, DIMENSION (:) :: ntblkrcounts, ntblkdisp
43  INTEGER, ALLOCATABLE, DIMENSION (:) :: nsrcounts, nsdisp
44  INTEGER, PRIVATE :: mmax, nmax, tmax
45 
46  LOGICAL :: PARVMEC=.false.
47  LOGICAL :: LV3FITCALL=.false.
48  LOGICAL :: LIFFREEB=.false.
49  LOGICAL :: LPRECOND=.false.
50 
51  LOGICAL :: lactive=.false.
52  LOGICAL :: vlactive=.false.
53 
54  CHARACTER*100 :: envvar
55  CHARACTER*100 :: envval
56 
57  REAL(dp) :: bcyclic_comp_time
58  REAL(dp) :: bcyclic_comm_time
59  REAL(dp) :: waitall_time
60  REAL(dp) :: dgemm_time
61  REAL(dp) :: dgetrf_time
62  REAL(dp) :: dgetrs_time
63  REAL(dp) :: dgemv_time
64  REAL(dp) :: ForwardSolveLoop_time
65  REAL(dp), DIMENSION(12) :: t
66 
67  REAL(dp) :: allgather_time=0
68  REAL(dp) :: allreduce_time=0
69  REAL(dp) :: broadcast_time=0
70  REAL(dp) :: sendrecv_time=0
71  REAL(dp) :: scatter_time=0
72 
73 !
74 ! OVERLOADED FUNCTIONS
75 !
76  INTERFACE copylastntype
77  MODULE PROCEDURE copy4lastntype, copy1lastntype, &
78  copym4lastntype, copym1lastntype
79  END INTERFACE
80 
81 CONTAINS
82 
83  !--------------------------------------------------------------------------
84  ! Read in environment variables
85  !--------------------------------------------------------------------------
86  SUBROUTINE myenvvariables
87 
88  parvmec=.true.
89  envvar='PARVMEC'
90  CALL getenv(envvar,envval)
91  IF (envval.EQ.'FALSE'.OR.envval.EQ.'false' &
92  .OR.envval.EQ.'F'.OR.envval.EQ.'f') THEN
93  parvmec=.false.
94  END IF
95 
96  lv3fitcall=.false.
97  envvar='LV3FITCALL'
98  CALL getenv(envvar,envval)
99  IF (envval.EQ.'TRUE'.OR.envval.EQ.'true' &
100  .OR.envval.EQ.'T'.OR.envval.EQ.'t') THEN
101  lv3fitcall=.true.
102  END IF
103 
104  END SUBROUTINE myenvvariables
105  !--------------------------------------------------------------------------
106 
107 
108  !--------------------------------------------------------------------------
109  ! Declarations of all timers to be used in the parallel implementation
110  !--------------------------------------------------------------------------
111  SUBROUTINE initializeparallel
112 
113 #if defined(MPI_OPT)
114  CALL mpi_init(mpi_err)
115  CALL mpi_comm_rank(mpi_comm_world,grank,mpi_err)
116  CALL mpi_comm_size(mpi_comm_world,gnranks,mpi_err)
117 #endif
118 
119  END SUBROUTINE initializeparallel
120  !--------------------------------------------------------------------------
121 
122  !--------------------------------------------------------------------------
123  SUBROUTINE initrunvmec(INCOMM, LVALUE)
124  INTEGER, INTENT(IN) :: INCOMM
125  LOGICAL, INTENT(IN) :: LVALUE
126 
127 ! NS_RESLTN = 0 ! SAL 070619
128  liffreeb = lvalue
129 #if defined(MPI_OPT)
130  CALL mpi_comm_dup(incomm,runvmec_comm_world,mpi_err)
131  CALL mpi_comm_rank(runvmec_comm_world,grank,mpi_err)
132  CALL mpi_comm_size(runvmec_comm_world,gnranks,mpi_err)
133 #endif
134  END SUBROUTINE initrunvmec
135  !--------------------------------------------------------------------------
136 
137  !--------------------------------------------------------------------------
138  SUBROUTINE initsurfacecomm(ns, nzeta, ntheta3, ntmax, ntor, mpol1)
139 
140 INTEGER, INTENT(IN) :: ns, nzeta, ntheta3
141  INTEGER, INTENT(IN) :: ntmax, ntor, mpol1
142  INTEGER :: i
143  LOGICAL :: FIRSTPASS = .true.
144 
145  par_ns = ns
146  par_nzeta = nzeta
147  par_ntheta3 = ntheta3
148  par_ntmax = ntmax
149  par_ntor = ntor
150  par_mpol1 = mpol1
151  par_nznt = par_nzeta*par_ntheta3
152  blocksize = (par_mpol1 + 1)*(par_ntor+1)
153  ntmaxblocksize = 3*ntmax*blocksize
154  mmax = par_mpol1 + 1
155  nmax = par_ntor + 1
156  tmax = 3*par_ntmax
157  ns_resltn = ns_resltn + 1
158 
159  IF (lv3fitcall) THEN
160  IF (last_ns.NE.par_ns) THEN
161  CALL setsurfacecommunicator
162  CALL setsurfacepartitions
163  CALL setsurfacepartitionarrays
164  last_ns = par_ns
165  END IF
166  ELSE
167  CALL setsurfacecommunicator
168  CALL setsurfacepartitions
169  CALL setsurfacepartitionarrays
170  last_ns = par_ns
171 
172  IF (lactive) THEN
173  IF (grank.EQ.0) THEN
174  IF (firstpass) THEN
175  CALL setoutputfile(rank,nranks,'parvmecinfo')
176  WRITE(tofu,*)"============================================================"
177  WRITE(tofu,*) 'PARVMEC = ',parvmec
178  WRITE(tofu,*) "> available processor count:", gnranks
179  WRITE(tofu,*) '> global rank:', grank
180  WRITE(tofu,*) '> nzeta: ', par_nzeta
181  WRITE(tofu,*) '> ntheta3: ', par_ntheta3
182  WRITE(tofu,*) '> ntor: ', par_ntor
183  WRITE(tofu,*) '> mpol1: ', par_mpol1
184  WRITE(tofu,*) '> ntmax: ', par_ntmax
185  WRITE(tofu,*) '> blocksize: ',(par_mpol1+1)*(par_ntor+1)
186  WRITE(tofu,*)"============================================================"
187  WRITE(tofu,*)
188  CALL flush(tofu)
189  firstpass = .false.
190  END IF
191  WRITE(tofu,*) ">>> grid resolution: ",par_ns
192  WRITE(tofu,*) ">>> active processors: ",nranks
193  WRITE(tofu,*) ">>> xc/gc size: ", par_ns*(par_ntor+1)*(par_mpol1+1)*3*ntmax
194  WRITE(tofu,*)"------------------------------------------------------------"
195  WRITE(tofu,*)
196  CALL flush(tofu)
197  END IF
198  END IF
199  END IF
200 
201  END SUBROUTINE initsurfacecomm
202  !------------------------------------------------
203 
204  !------------------------------------------------
205  ! Setting up the communicator for parallel surface
206  ! computations
207  !------------------------------------------------
208  SUBROUTINE setsurfacecommunicator
209 
210  INTEGER :: num_active, color
211 
212  num_active = min(gnranks, par_ns/2)
213 
214  IF (grank.LT.num_active) THEN
215  color=1
216  ELSE
217  color=0
218  END IF
219 #if defined(MPI_OPT)
220  CALL mpi_comm_split(runvmec_comm_world, color, grank, ns_comm, &
221  mpi_err)
222 
223  IF (color .eq. 1) THEN
224  CALL mpi_comm_size(ns_comm,nranks,mpi_err)
225  IF (nranks.NE.num_active) THEN
226  stop 'num_active != nranks in InitSurfaceCommunicator!'
227  END IF
228  lactive = .true.
229  CALL mpi_comm_rank(ns_comm,rank,mpi_err)
230  ELSE
231  nranks = 1
232  rank = 0
233  lactive = .false.
234  END IF
235 #endif
236  END SUBROUTINE setsurfacecommunicator
237  !------------------------------------------------
238 
239  !------------------------------------------------
240  ! Setting up the partitions for parallel surface
241  ! computations
242  !------------------------------------------------
243  SUBROUTINE setsurfacepartitions
244 
245  INTEGER :: mypart, local_err
246 #if defined(MPI_OPT)
247  IF (par_ns.LT.nranks) THEN
248  IF(grank.EQ.0) print *,"NS is less than NRANKS. Aborting!"
249  CALL stopmpi(5645)
250  END IF
251 
252  mypart=par_ns/nranks
253  IF (rank.LT.mod(par_ns,nranks)) mypart = mypart + 1
254  IF (mod(par_ns,nranks).NE.0) THEN
255  IF (rank.LT.mod(par_ns,nranks)) THEN
256  tlglob = rank*mypart
257  ELSE IF (rank .GE. mod(par_ns,nranks)) THEN
258  tlglob = mod(par_ns, nranks)*(mypart + 1) &
259  + (rank - mod(par_ns, nranks))*mypart
260  END IF
261  ELSE
262  tlglob = rank*mypart
263  END IF
264 
265  tlglob = tlglob + 1
266  trglob = tlglob + mypart - 1
267 
268  t1lglob = tlglob - 1; IF (rank.EQ.0) t1lglob=1
269  t1rglob = trglob + 1; IF (rank.EQ.nranks-1) t1rglob=par_ns
270 
271  IF(mypart.LT.2) THEN
272  CALL mpi_barrier(ns_comm,mpi_err)
273  WRITE(tofu,*) '***********************************************************'
274  WRITE(tofu,*) '* This version is not yet tested for mynsnum <= 2. Aborting!'
275  WRITE(tofu,*) '***********************************************************'
276  IF (rank.EQ.0) THEN
277  WRITE(*,*)
278  WRITE(*,*) '***********************************************************'
279  WRITE(*,*) '* This version is not yet tested for mynsnum <= 2. Aborting!'
280  WRITE(*,*) '***********************************************************'
281  WRITE(*,*)
282  END IF
283  CALL mpi_abort(ns_comm,local_err,mpi_err)
284  END IF
285 #endif
286  END SUBROUTINE setsurfacepartitions
287  !------------------------------------------------
288 
289  !------------------------------------------------
290  ! Setting up the partition arrays for parallel surface
291  ! computations
292  !------------------------------------------------
293  SUBROUTINE setsurfacepartitionarrays
294 
295  INTEGER, ALLOCATABLE, DIMENSION(:) :: localpart
296  INTEGER :: i, smallpart, largepart
297  INTEGER, SAVE :: lastns
298 
299  smallpart = par_ns/nranks; largepart = smallpart
300  IF(mod(par_ns,nranks) .GT. 0) largepart = largepart + 1
301 
302  ALLOCATE (localpart(nranks))
303  DO i = 0, nranks - 1
304  localpart(i + 1) = smallpart
305 
306  IF (mod(par_ns,nranks) .GT. 0 .and. &
307  i .LT. mod(par_ns, nranks)) THEN
308  localpart(i + 1) = localpart(i + 1) + 1
309  END IF
310  END DO
311 
312  IF(ALLOCATED(tlglob_arr)) DEALLOCATE(tlglob_arr)
313  IF(ALLOCATED(trglob_arr)) DEALLOCATE(trglob_arr)
314  ALLOCATE (tlglob_arr(nranks), trglob_arr(nranks))
315 
316  tlglob_arr(1) = 1
317  DO i = 2, nranks
318  tlglob_arr(i) = tlglob_arr(i - 1) + localpart(i - 1)
319  END DO
320  DO i = 1, nranks
321  trglob_arr(i) = tlglob_arr(i) + localpart(i) - 1
322  END DO
323 
324  DEALLOCATE (localpart)
325 
326  CALL computensallgatherparameters(nranks)
327  CALL computeblockallgatherparameters(nranks)
328  CALL computentmaxblockallgatherparameters(nranks)
329 
330  END SUBROUTINE setsurfacepartitionarrays
331  !------------------------------------------------
332 
333  !------------------------------------------------
334  ! Compute AllGather vector variant parameters for
335  ! blocksized movements.
336  !------------------------------------------------
337  SUBROUTINE computentmaxblockallgatherparameters(activeranks)
338 
339  INTEGER :: activeranks
340  INTEGER :: i
341  IF(.NOT.ALLOCATED(ntblkrcounts)) ALLOCATE(ntblkrcounts(activeranks))
342  IF(.NOT.ALLOCATED(ntblkdisp)) ALLOCATE(ntblkdisp(activeranks))
343  DO i = 1, activeranks
344  ntblkrcounts(i) = (trglob_arr(i) - tlglob_arr(i) + 1)*ntmaxblocksize
345  END DO
346  ntblkdisp(1) = 0
347  DO i = 2, activeranks
348  ntblkdisp(i) = ntblkdisp(i - 1) + ntblkrcounts(i - 1)
349  END DO
350  END SUBROUTINE computentmaxblockallgatherparameters
351  !------------------------------------------------
352 
353  !------------------------------------------------
354  ! Compute AllGather vector variant parameters for
355  ! blocksized movements.
356  !------------------------------------------------
357  SUBROUTINE computeblockallgatherparameters(activeranks)
358 
359  INTEGER :: activeranks
360  INTEGER :: i
361  IF(.NOT.ALLOCATED(blkrcounts)) ALLOCATE(blkrcounts(activeranks))
362  IF(.NOT.ALLOCATED(blkdisp)) ALLOCATE(blkdisp(activeranks))
363  DO i = 1, activeranks
364  blkrcounts(i) = (trglob_arr(i) - tlglob_arr(i) + 1)*blocksize
365  END DO
366  blkdisp(1)=0
367  DO i = 2, activeranks
368  blkdisp(i) = blkdisp(i - 1) + blkrcounts(i - 1)
369  END DO
370  END SUBROUTINE computeblockallgatherparameters
371  !------------------------------------------------
372 
373  !------------------------------------------------
374  ! Compute AllGather vector variant parameters for
375  ! blocksized movements.
376  !------------------------------------------------
377  SUBROUTINE computensallgatherparameters(activeranks)
378 
379  INTEGER :: activeranks
380  INTEGER :: i
381  IF(.NOT.ALLOCATED(nsrcounts)) ALLOCATE(nsrcounts(activeranks))
382  IF(.NOT.ALLOCATED(nsdisp)) ALLOCATE(nsdisp(activeranks))
383  DO i = 1, activeranks
384  nsrcounts(i) = (trglob_arr(i) - tlglob_arr(i) + 1)
385  END DO
386  nsdisp(1) = 0
387  DO i = 2, activeranks
388  nsdisp(i) = nsdisp(i - 1) + nsrcounts(i - 1)
389  END DO
390  END SUBROUTINE computensallgatherparameters
391  !------------------------------------------------
392 
393  !--------------------------------------------------------------------------
394  SUBROUTINE finalizesurfacecomm(INCOMM)
395 
396  INTEGER, INTENT(INOUT) :: INCOMM
397 #if defined(MPI_OPT)
398  CALL mpi_comm_free(incomm, mpi_err)
399  lactive = .false.
400  IF (ALLOCATED(ntblkrcounts)) DEALLOCATE(ntblkrcounts)
401  IF (ALLOCATED(ntblkdisp)) DEALLOCATE(ntblkdisp)
402  IF (ALLOCATED(blkrcounts)) DEALLOCATE(blkrcounts)
403  IF (ALLOCATED(blkdisp)) DEALLOCATE(blkdisp)
404  IF (ALLOCATED(nsrcounts)) DEALLOCATE(nsrcounts)
405  IF (ALLOCATED(nsdisp)) DEALLOCATE(nsdisp)
406 #endif
407  END SUBROUTINE finalizesurfacecomm
408  !--------------------------------------------------------------------------
409 
410  !--------------------------------------------------------------------------
411  SUBROUTINE finalizerunvmec(INCOMM)
412 
413  INTEGER, INTENT(INOUT) :: INCOMM
414 #if defined(MPI_OPT)
415  CALL mpi_comm_free(incomm, mpi_err)
416  IF(liffreeb) CALL mpi_comm_free(vac_comm,mpi_err)
417  incomm=0
418  vac_comm = 0
419  rank = 0
420  IF (.not.lv3fitcall) par_ns = 0
421  nranks = 1
422  grank = 0
423  gnranks = 1
424  vrank = 0
425  vnranks = 1
426  last_ns = -1
427  ns_resltn=0
428  vlactive = .false.
429 #endif
430  END SUBROUTINE finalizerunvmec
431  !--------------------------------------------------------------------------
432 
433  !------------------------------------------------
434  ! Setting up the communicator for parallel vacuum
435  ! computations
436  ! nuv = nzeta*ntheta1
437  ! nuv3 = nzeta*ntheta3
438  !------------------------------------------------
439  SUBROUTINE setvacuumcommunicator(nuv, nuv3, mgs)
440 
441  INTEGER, INTENT(IN) :: nuv, nuv3, mgs
442  INTEGER :: num_active, color, mypart
443 
444  par_nuv3 = nuv3
445  par_nuv = nuv
446 #if defined(MPI_OPT)
447  num_active = min(gnranks, par_nuv3)
448 
449  IF(grank .LT. num_active) THEN
450  color = 1
451  ELSE
452  color = 0
453  END IF
454  CALL mpi_comm_split(runvmec_comm_world, color, grank, vac_comm, &
455  mpi_err)
456  IF (color .eq. 1) THEN
457  CALL mpi_comm_rank(vac_comm,vrank,mpi_err)
458  CALL mpi_comm_size(vac_comm,vnranks,mpi_err)
459  CALL setvacuumpartitions(nuv3,nuv3min, nuv3max)
460  CALL setnuv3partitionarrays
461  vlactive=.true.
462  ELSE
463  vnranks = 1
464  vrank = 0
465  vlactive=.false.
466  END IF
467 #else
468  nuv3min = 1
469  nuv3max = nuv3
470 #endif
471  END SUBROUTINE setvacuumcommunicator
472  !------------------------------------------------
473 
474 
475  !------------------------------------------------
476  ! Setting up the partitions for parallel vacuum
477  ! computations
478  !------------------------------------------------
479  SUBROUTINE setvacuumpartitions(num, left, right)
480 
481  INTEGER, INTENT(IN) :: num
482  INTEGER, INTENT(INOUT) :: left, right
483  INTEGER :: mypart
484 
485  IF (num.LT.vnranks) THEN
486  IF(grank.EQ.0) print *,"NUM is less than VNRANKS. Aborting!"
487  CALL stopmpi(456745)
488  END IF
489 
490  mypart = num/vnranks
491  IF (vrank .LT. mod(num,vnranks)) mypart = mypart + 1
492  IF (mod(num,vnranks).NE.0) THEN
493  IF (vrank.LT.mod(num,vnranks)) THEN
494  left = vrank*mypart
495  ELSE IF (vrank.GE.mod(num,vnranks)) THEN
496  left = mod(num,vnranks)*(mypart + 1) &
497  + (vrank - mod(num,vnranks))*mypart
498  END IF
499  ELSE
500  left = vrank*mypart
501  END IF
502 
503  left = left + 1
504  right = left + mypart - 1
505 
506  END SUBROUTINE setvacuumpartitions
507  !------------------------------------------------
508 
509  !------------------------------------------------
510  ! Setting up the partition arrays for parallel vacuum
511  ! computations
512  !------------------------------------------------
513  SUBROUTINE setnuv3partitionarrays
514 
515  INTEGER, ALLOCATABLE, DIMENSION(:) :: localpart
516  INTEGER :: i, smallpart, largepart
517 
518  IF (.NOT.ALLOCATED(nuv3min_arr)) THEN
519  ALLOCATE(nuv3min_arr(vnranks), nuv3max_arr(vnranks))
520  END IF
521 
522  smallpart = par_nuv3/vnranks
523  largepart = smallpart
524  IF (mod(par_nuv3,vnranks).GT.0) THEN
525  largepart = largepart+1
526  END IF
527 
528  ALLOCATE (localpart(vnranks))
529  DO i = 0, vnranks - 1
530 
531  localpart(i + 1) = smallpart
532 
533  IF (mod(par_nuv3, vnranks) .GT. 0 .AND. &
534  i .LT. mod(par_nuv3, vnranks)) THEN
535  localpart(i + 1) = localpart(i + 1) + 1
536  END IF
537  END DO
538 
539  nuv3min_arr(1) = 1
540  DO i = 2, vnranks
541  nuv3min_arr(i) = nuv3min_arr(i - 1) + localpart(i - 1)
542  END DO
543  DO i = 1, vnranks
544  nuv3max_arr(i) = nuv3min_arr(i) + localpart(i) - 1
545  END DO
546 
547  DEALLOCATE (localpart)
548 
549  END SUBROUTINE setnuv3partitionarrays
550  !------------------------------------------------
551 
552  !------------------------------------------------
553  SUBROUTINE finalizeparallel
554 
555  INTEGER :: istat
556 
557  envvar = 'LPRECOND'
558  CALL getenv(envvar, envval)
559  IF (envval .EQ. 'TRUE') THEN
560  lprecond = .true.
561  ELSE IF (envval .EQ. 'FALSE') THEN
562  lprecond=.false.
563  END IF
564 
565  IF (grank.EQ.0) THEN
566  WRITE(*,*)
567  WRITE(*,'(1x,a,i4)') 'NO. OF PROCS: ',gnranks
568  WRITE(*,100) 'PARVMEC : ',parvmec
569  WRITE(*,100) 'LPRECOND : ',lprecond
570  WRITE(*,100) 'LV3FITCALL : ',lv3fitcall
571  END IF
572  100 FORMAT(1x,a,l4)
573 
574 #if defined(MPI_OPT)
575  CALL mpi_finalize(istat)
576 #endif
577  END SUBROUTINE finalizeparallel
578  !------------------------------------------------
579 
580  !------------------------------------------------
581  ! Print debugging output to compare aray values
582  !------------------------------------------------
583  SUBROUTINE printnsarray(arr, left, right, fileno, ifstop, string)
584 
585  INTEGER, INTENT(IN) :: fileno, left, right
586  LOGICAL, INTENT(IN) :: ifstop
587  REAL(dp), DIMENSION (par_ns + 1), INTENT(IN) :: arr
588  CHARACTER*(*), INTENT(IN) :: string
589 
590  INTEGER :: i
591 
592 #if defined(_WIN32)
593  RETURN
594 #endif
595 
596  DO i=left, right
597  WRITE(fileno + rank,50) string, i, arr(i)
598  CALL flush(fileno+rank)
599  END DO
600  WRITE(fileno+rank,*)
601  CALL flush(fileno+rank)
602  IF(ifstop) stop 'STOPPED CODE'
603 
604  50 FORMAT(a6,1i6,1p,e20.6)
605 
606  END SUBROUTINE printnsarray
607  !------------------------------------------------
608 
609  !------------------------------------------------
610  SUBROUTINE stopmpi(code)
611 
612  INTEGER, INTENT(IN) :: code
613 
614 #if defined(MPI_OPT)
615  CALL mpi_barrier(mpi_comm_world, mpi_err)
616 #endif
617  WRITE(*,*) 'Stopping program with code:', code
618  stop
619  END SUBROUTINE stopmpi
620  !------------------------------------------------
621 
622  !------------------------------------------------
623  SUBROUTINE parallel2serial4x(inarr,outarr)
624 
625  REAL(dp), INTENT(IN) :: inarr(par_ns*(par_ntor+1)*(par_mpol1+1)*3*(par_ntmax))
626  REAL(dp), INTENT(OUT) :: outarr(par_ns*(par_ntor+1)*(par_mpol1+1)*3*(par_ntmax))
627  INTEGER :: i, j, k, l, lk, lks, lkp
628 
629  lks=0
630  DO l=1, 3*par_ntmax
631  DO k=0, par_mpol1
632  DO j=0, par_ntor
633  DO i=1, par_ns
634  lks = lks +1
635  lkp = j + (par_ntor + 1)*k &
636  + (par_ntor + 1)*(par_mpol1 + 1)*(i - 1) &
637  + (par_ntor + 1)*(par_mpol1 + 1)*par_ns*(l - 1) + 1
638  outarr(lks) = inarr(lkp)
639  END DO
640  END DO
641  END DO
642  END DO
643 
644  END SUBROUTINE parallel2serial4x
645  !------------------------------------------------
646 
647  !------------------------------------------------
648  SUBROUTINE serial2parallel4x (inarr,outarr)
649 
650  REAL(dp), INTENT(IN) :: inarr(par_ns*(par_ntor+1)*(par_mpol1+1)*3*(par_ntmax))
651  REAL(dp), INTENT(OUT) :: outarr(par_ns*(par_ntor+1)*(par_mpol1+1)*3*(par_ntmax))
652  INTEGER :: i, j, k, l, lk, lks, lkp
653 
654  lks=0
655  DO l=1, 3*par_ntmax
656  DO k=0, par_mpol1
657  DO j=0, par_ntor
658  DO i=1, par_ns
659  lks = lks+1
660  lkp = j + (par_ntor + 1)*k &
661  + (par_ntor + 1)*(par_mpol1 + 1)*(i - 1) &
662  + (par_ntor + 1)*(par_mpol1 + 1)*par_ns*(l - 1) + 1
663  outarr(lkp) = inarr(lks)
664  END DO
665  END DO
666  END DO
667  END DO
668 
669  END SUBROUTINE serial2parallel4x
670  !------------------------------------------------
671 
672  !------------------------------------------------
673  SUBROUTINE parallel2serial2x(inarr, outarr)
674 
675  REAL(dp), DIMENSION(par_nzeta,par_ntheta3,par_ns), INTENT(IN) :: inarr
676  REAL(dp), DIMENSION(par_ns*par_nzeta*par_ntheta3+1), INTENT(OUT) :: outarr
677  INTEGER :: i, j, k, l
678 
679  l=0
680  DO k = 1, par_ntheta3
681  DO j = 1, par_nzeta
682  DO i = 1, par_ns
683  l = l + 1
684  outarr(l) = inarr(j,k,i)
685  END DO
686  END DO
687  END DO
688 
689  END SUBROUTINE parallel2serial2x
690  !------------------------------------------------
691 
692  !------------------------------------------------
693  SUBROUTINE rprintoutlineararray(arr, left, right, flag, fileno)
694 
695  REAL(dp), DIMENSION(:) :: arr
696  INTEGER, INTENT(IN) :: fileno, left, right
697  LOGICAL, INTENT(IN) :: flag !flag: TRUE for parallel, FALSE for serial
698  INTEGER :: i, j, k, l, lk
699 
700  REAL(dp), ALLOCATABLE, DIMENSION(:) :: tmp
701  ALLOCATE(tmp(ntmaxblocksize*par_ns))
702 
703  CALL tolastntype(arr,tmp)
704  lk = 0
705  DO l = 1, 3*par_ntmax
706  DO k = 0, par_mpol1
707  DO j = 0, par_ntor
708  DO i = 1, par_ns
709  lk = lk + 1
710  IF (flag) THEN
711  lk = j + nmax*(k + mmax*((i - 1) + par_ns*(l - 1))) + 1
712  END IF
713  IF (left .LE. i .AND. i .LE. right) THEN
714  WRITE(fileno+rank,50) i, j, k, l, tmp(lk)
715  CALL flush(fileno+rank)
716  END IF
717  END DO
718  END DO
719  END DO
720  END DO
721  DEALLOCATE(tmp)
722 
723  50 FORMAT(4i6,1p,e24.14)
724 
725  END SUBROUTINE rprintoutlineararray
726  !------------------------------------------------
727 
728 
729  !------------------------------------------------
730  SUBROUTINE printoutlineararray(arr, left, right, flag, fileno)
731 
732  REAL(dp), DIMENSION(ntmaxblocksize*par_ns) :: arr
733  INTEGER, INTENT(IN) :: fileno, left, right
734  LOGICAL, INTENT(IN) :: flag !flag: TRUE for parallel, FALSE for serial
735  INTEGER :: i, j, k, l, lk
736 
737  lk=0
738  DO l = 1, 3*par_ntmax
739  DO k = 0, par_mpol1
740  DO j = 0, par_ntor
741  DO i = 1, par_ns
742  lk=lk + 1
743  IF(flag) THEN
744  lk = j + nmax*(k + mmax*((i - 1) + par_ns*(l-1)))+1
745  END IF
746  IF (left .LE. i .AND. i .LE. right) THEN
747  WRITE(fileno+rank,50) i, j, k, l, arr(lk)
748  CALL flush(fileno+rank)
749  END IF
750  END DO
751  END DO
752  END DO
753  END DO
754 
755  50 FORMAT(4i6,1p,e24.14)
756 
757  END SUBROUTINE printoutlineararray
758  !------------------------------------------------
759 
760  !------------------------------------------------
761  ! Prints out [nsmin, nsmax] part of a (nzeta, ntheta,ns) array
762  !------------------------------------------------
763  SUBROUTINE printparallelijsparray (arr, left, right, fileno, ifstop, string)
764 
765  INTEGER, INTENT(IN) :: fileno, left, right
766  LOGICAL, INTENT(IN) :: ifstop
767  REAL(dp), DIMENSION (par_nzeta,par_ntheta3,par_ns), INTENT(IN) :: arr
768  CHARACTER*(*), INTENT(IN) :: string
769 
770  INTEGER :: i, j, k
771 
772  DO i=left, right
773  DO j=1, par_nzeta
774  DO k=1, par_ntheta3
775  WRITE(fileno+rank,50) string, i, j, k, arr(j,k,i)
776  CALL flush(fileno+rank)
777  END DO
778  END DO
779  END DO
780  WRITE(fileno+rank,*)
781  CALL flush(fileno+rank)
782  IF(ifstop) stop 'STOPPED PARALLEL CODE'
783  50 FORMAT(a,3i6,1p,e24.14)
784 
785  END SUBROUTINE printparallelijsparray
786  !------------------------------------------------
787 
788  !------------------------------------------------
789  ! Prints out [nsmin, nsmax] part of a (ns, nzeta, ntheta) array
790  !------------------------------------------------
791  SUBROUTINE printserialijsparray (arr, left, right, fileno, ifstop, string)
792 
793  INTEGER, INTENT(IN) :: fileno, left, right
794  LOGICAL, INTENT(IN) :: ifstop
795  REAL(dp), DIMENSION (par_ns,par_nzeta,par_ntheta3), INTENT(IN) :: arr
796  CHARACTER*(*), INTENT(IN) :: string
797 
798  INTEGER :: i, j, k
799 
800  DO i = left, right
801  DO j = 1, par_nzeta
802  DO k = 1, par_ntheta3
803  WRITE(fileno+rank,50) string, i, j, k, arr(i, j, k)
804  CALL flush(fileno+rank)
805  END DO
806  END DO
807  END DO
808  WRITE(fileno+rank,*)
809  CALL flush(fileno+rank)
810  IF(ifstop) stop 'STOPPED SERIAL CODE'
811  50 FORMAT(a,3i6,1p,e24.14)
812 
813  END SUBROUTINE printserialijsparray
814  !------------------------------------------------
815 
816  !------------------------------------------------
817  ! Prints out [nsmin, nsmax] part of a (ntor,mpol1,ns) array
818  !------------------------------------------------
819  SUBROUTINE printparallelmnsparray (arr, left, right, fileno, ifstop, string)
820 
821  INTEGER, INTENT(IN) :: fileno, left, right
822  LOGICAL, INTENT(IN) :: ifstop
823  REAL(dp), DIMENSION (0:par_ntor,0:par_mpol1,1:par_ns), INTENT(IN) :: arr
824  CHARACTER*(*), INTENT(IN) :: string
825 
826  INTEGER :: i, j, k
827 
828  DO i=left, right
829  DO j=0, par_mpol1
830  DO k=0, par_ntor
831  WRITE(fileno+rank,50) string, i, j, k, arr(k,j,i)
832  CALL flush(fileno+rank)
833  END DO
834  END DO
835  END DO
836  WRITE(fileno+rank,*)
837  CALL flush(fileno+rank)
838  IF(ifstop) stop 'STOPPED PARALLEL CODE'
839  50 FORMAT(a,3i6,1p,e20.12)
840 
841  END SUBROUTINE printparallelmnsparray
842  !------------------------------------------------
843 
844  !------------------------------------------------
845  ! Prints out [nsmin, nsmax] part of a (ns, ntor, mpol1) array
846  !------------------------------------------------
847  SUBROUTINE printserialmnsparray (arr, left, right, fileno, ifstop, string)
848 
849  INTEGER, INTENT(IN) :: fileno, left, right
850  LOGICAL, INTENT(IN) :: ifstop
851  REAL(dp), DIMENSION (1:par_ns,0:par_ntor,0:par_mpol1), INTENT(IN) :: arr
852  CHARACTER*(*), INTENT(IN) :: string
853 
854  INTEGER :: i, j, k
855 
856  DO i=left, right
857  DO k=0, par_mpol1
858  DO j=0, par_ntor
859  WRITE(fileno+rank,50) string, i, j, k, arr(i, j, k)
860  CALL flush(fileno+rank)
861  END DO
862  END DO
863  END DO
864  WRITE(fileno+rank,*)
865  CALL flush(fileno+rank)
866  IF(ifstop) stop 'STOPPED SERIAL CODE'
867  50 FORMAT(a,3i6,1p,e20.12)
868 
869  END SUBROUTINE printserialmnsparray
870  !------------------------------------------------
871 
872  !------------------------------------------------
873  SUBROUTINE gather4xarray(arr)
874  !-----------------------------------------------
875 
876  !-----------------------------------------------
877  REAL(dp), DIMENSION(0:par_ntor,0:par_mpol1,par_ns,3*par_ntmax), INTENT(INOUT) :: arr
878  INTEGER :: i
879  INTEGER :: blksize, numjs
880  INTEGER, ALLOCATABLE, DIMENSION(:) :: counts, disps
881  REAL(dp), ALLOCATABLE, DIMENSION(:,:,:,:) :: sendbuf
882  REAL(dp), ALLOCATABLE, DIMENSION(:) :: recvbuf
883  REAL(dp) :: allgvton, allgvtoff
884  !-----------------------------------------------
885  IF (nranks.EQ.1 .OR. .NOT.lactive) THEN
886  RETURN
887  END IF
888 
889  CALL second0(allgvton)
890 
891  blksize = (par_ntor + 1)*(par_mpol1 + 1)*3*par_ntmax
892  numjs = trglob - tlglob + 1
893  ALLOCATE(sendbuf(0:par_ntor,0:par_mpol1,numjs,1:3*par_ntmax))
894  ALLOCATE(recvbuf(par_ns*blksize))
895  ALLOCATE(counts(nranks),disps(nranks))
896 
897  DO i = 1, nranks
898  counts(i) = (trglob_arr(i) - tlglob_arr(i) + 1)*blksize
899  END DO
900 
901  disps(1) = 0
902  DO i = 2, nranks
903  disps(i) = disps(i - 1) + counts(i - 1)
904  END DO
905 #if defined(MPI_OPT)
906  sendbuf(0:par_ntor,0:par_mpol1,1:numjs,1:3*par_ntmax) = &
907  arr(0:par_ntor,0:par_mpol1,tlglob:trglob,1:3*par_ntmax)
908  CALL mpi_allgatherv(sendbuf, numjs*blksize, mpi_real8, recvbuf, &
909  counts, disps, mpi_real8, ns_comm, mpi_err)
910  DO i = 1, nranks
911  numjs = trglob_arr(i) - tlglob_arr(i) + 1
912  arr(0:par_ntor,0:par_mpol1,tlglob_arr(i):trglob_arr(i),1:3*par_ntmax) = &
913  reshape(recvbuf(disps(i)+1:disps(i)+counts(i)), &
914  (/par_ntor+1,par_mpol1+1,numjs,3*par_ntmax/))
915  END DO
916 #endif
917  DEALLOCATE(sendbuf, recvbuf)
918  DEALLOCATE(counts, disps)
919  CALL second0(allgvtoff)
920  allgather_time = allgather_time + (allgvtoff-allgvton)
921  END SUBROUTINE gather4xarray
922  !------------------------------------------------
923 
924  !------------------------------------------------
925  !
926  !------------------------------------------------
927  SUBROUTINE gatherreordered4xarray(arr)
928  !-----------------------------------------------
929 
930  !-----------------------------------------------
931  REAL(dp), DIMENSION(0:par_ntor,0:par_mpol1,3*par_ntmax,par_ns), INTENT(INOUT) :: arr
932  INTEGER :: i
933  INTEGER :: blksize, numjs
934  INTEGER, ALLOCATABLE, DIMENSION(:) :: counts, disps
935  REAL(dp) :: allgvton, allgvtoff
936  !-----------------------------------------------
937  IF (nranks.EQ.1 .OR. .NOT.lactive) THEN
938  RETURN
939  END IF
940 
941  CALL second0(allgvton)
942 
943  blksize = (par_ntor + 1)*(par_mpol1 + 1)*3*par_ntmax
944  numjs = trglob - tlglob + 1
945  ALLOCATE(counts(nranks),disps(nranks))
946 
947  DO i = 1, nranks
948  counts(i) = (trglob_arr(i) - tlglob_arr(i) + 1)*blksize
949  END DO
950 
951  disps(1) = 0
952  DO i = 2, nranks
953  disps(i) = disps(i - 1) + counts(i - 1)
954  END DO
955 #if defined(MPI_OPT)
956  CALL mpi_allgatherv(mpi_in_place, numjs*blksize, mpi_real8, arr, &
957  counts, disps, mpi_real8, ns_comm, mpi_err)
958 #endif
959  DEALLOCATE(counts, disps)
960  CALL second0(allgvtoff)
961  allgather_time = allgather_time + (allgvtoff-allgvton)
962  END SUBROUTINE gatherreordered4xarray
963  !------------------------------------------------
964 
965 
966  !------------------------------------------------
967  SUBROUTINE gather2xarray(arr)
968 
969  REAL(dp), DIMENSION(par_nznt,par_ns), INTENT(INOUT) :: arr
970  INTEGER :: i
971  INTEGER :: par_nsmin, par_nsmax, blksize, numjs
972  INTEGER, ALLOCATABLE, DIMENSION(:) :: counts, disps
973  REAL(dp) :: allgvton, allgvtoff
974  !-----------------------------------------------
975  IF (nranks.EQ.1 .OR. .NOT.lactive) THEN
976  RETURN
977  END IF
978 
979  CALL second0(allgvton)
980 
981  blksize = par_nznt
982  numjs = trglob - tlglob+1
983  ALLOCATE(counts(nranks),disps(nranks))
984 
985  DO i = 1, nranks
986  counts(i) =(trglob_arr(i) - tlglob_arr(i) + 1)*blksize
987  END DO
988 
989  disps(1) = 0
990  DO i = 2, nranks
991  disps(i) = disps(i - 1) + counts(i - 1)
992  END DO
993 
994 #if defined(MPI_OPT)
995  CALL mpi_allgatherv(mpi_in_place, numjs*blksize, mpi_real8, arr, &
996  counts, disps, mpi_real8, ns_comm, mpi_err)
997 #endif
998  DEALLOCATE(counts, disps)
999  CALL second0(allgvtoff)
1000  allgather_time = allgather_time + (allgvtoff-allgvton)
1001  END SUBROUTINE gather2xarray
1002  !------------------------------------------------
1003 
1004  !------------------------------------------------
1005  SUBROUTINE gather1xarray(arr)
1006  !-----------------------------------------------
1007  REAL(dp), DIMENSION(par_ns), INTENT(INOUT) :: arr
1008 
1009  INTEGER :: i, numjs
1010  INTEGER, ALLOCATABLE, DIMENSION(:) :: counts, disps
1011  REAL(dp) :: allgvton, allgvtoff
1012  !-----------------------------------------------
1013  IF (nranks.EQ.1 .OR. .NOT.lactive) THEN
1014  RETURN
1015  END IF
1016 
1017  CALL second0(allgvton)
1018 
1019  numjs = trglob - tlglob + 1
1020  ALLOCATE(counts(nranks),disps(nranks))
1021 
1022  DO i=1,nranks
1023  counts(i) = (trglob_arr(i) - tlglob_arr(i) + 1)
1024  END DO
1025 
1026  disps(1)=0
1027  DO i = 2,nranks
1028  disps(i)= disps(i - 1) +counts(i - 1)
1029  END DO
1030 
1031 #if defined(MPI_OPT)
1032  CALL mpi_allgatherv(mpi_in_place, numjs, mpi_real8, arr, counts, &
1033  disps, mpi_real8, ns_comm, mpi_err)
1034 #endif
1035  DEALLOCATE(counts, disps)
1036  CALL second0(allgvtoff)
1037  allgather_time = allgather_time + (allgvtoff - allgvton)
1038  END SUBROUTINE gather1xarray
1039  !------------------------------------------------
1040 
1041  !------------------------------------------------
1042  !------------------------------------------------
1043 
1044  !------------------------------------------------
1045  ! Setting the output files
1046  !------------------------------------------------
1047  SUBROUTINE setoutputfile(iam, nprocs, prefix)
1048 
1049  INTEGER, INTENT(IN) :: iam, nprocs
1050  CHARACTER (*), INTENT(IN) :: prefix
1051  INTEGER :: istat
1052  CHARACTER(100) :: fname, cfname
1053  CHARACTER(50) :: ciam, cnprocs
1054  CHARACTER(25) :: cprefix
1055 
1056  WRITE(ciam,*) iam
1057  WRITE(cnprocs,*) nprocs
1058  WRITE(cprefix,*) prefix
1059  ciam = adjustl(ciam)
1060  cnprocs = adjustl(cnprocs)
1061  cprefix = adjustl(cprefix)
1062  tofu = 4*nprocs + iam + 1000
1063 
1064  fname=trim(cprefix)//'.txt'
1065  OPEN(unit=tofu, file=fname, status="REPLACE", action="WRITE", &
1066  form="FORMATTED",position="APPEND", iostat=istat)
1067 
1068  END SUBROUTINE setoutputfile
1069  !------------------------------------------------
1070 
1071  !------------------------------------------------
1072  SUBROUTINE tolastns(inarr, outarr)
1073  REAL(dp), INTENT(IN), DIMENSION(0:par_ntor,0:par_mpol1,par_ns,3*par_ntmax) :: &
1074  inarr
1075  INTEGER :: i, j, k, l, cnt
1076  REAL(dp), INTENT(INOUT), DIMENSION(ntmaxblocksize,par_ns) :: outarr
1077 
1078  DO i = t1lglob, t1rglob
1079  cnt = 0
1080  DO l = 1, 3*par_ntmax
1081  DO k = 0, par_mpol1
1082  DO j = 0, par_ntor
1083  cnt = cnt + 1
1084  outarr(cnt,i) = inarr(j,k,i,l)
1085  END DO
1086  END DO
1087  END DO
1088  END DO
1089 
1090  END SUBROUTINE tolastns
1091  !------------------------------------------------
1092 
1093  !------------------------------------------------
1094  SUBROUTINE tolastntype(inarr, outarr)
1095  REAL(dp), INTENT(INOUT), DIMENSION(0:par_ntor,0:par_mpol1,par_ns,3*par_ntmax) :: outarr
1096  REAL(dp), INTENT(IN), DIMENSION(ntmaxblocksize,par_ns) :: inarr
1097  INTEGER :: i, j, k, l, cnt
1098 
1099  DO i=t1lglob, t1rglob
1100  cnt=0
1101  DO l=1, 3*par_ntmax
1102  DO k=0, par_mpol1
1103  DO j=0, par_ntor
1104  cnt=cnt+1
1105  outarr(j,k,i,l)=inarr(cnt,i)
1106  END DO
1107  END DO
1108  END DO
1109  END DO
1110  END SUBROUTINE tolastntype
1111  !------------------------------------------------
1112 
1113  !------------------------------------------------
1114  !
1115  !------------------------------------------------
1116  SUBROUTINE copylastns(a1, a2)
1117  REAL(dp), INTENT(IN), DIMENSION(ntmaxblocksize,par_ns) :: a1
1118  REAL(dp), INTENT(INOUT), DIMENSION(ntmaxblocksize,par_ns) :: a2
1119  a2(:,t1lglob:t1rglob) = a1(:,t1lglob:t1rglob)
1120  END SUBROUTINE copylastns
1121  !------------------------------------------------
1122 
1123  !------------------------------------------------
1124  SUBROUTINE copy1lastntype(a1, a2)
1125  REAL(dp), INTENT(IN), DIMENSION(ntmaxblocksize*par_ns) :: a1
1126  REAL(dp), INTENT(INOUT), DIMENSION(ntmaxblocksize*par_ns) :: a2
1127  CALL copy4lastntype(a1, a2)
1128  END SUBROUTINE copy1lastntype
1129  !------------------------------------------------
1130 
1131  !------------------------------------------------
1132  SUBROUTINE copy4lastntype(a1, a2)
1133  REAL(dp), INTENT(IN), DIMENSION(0:par_ntor,0:par_mpol1,par_ns,3*par_ntmax) :: a1
1134  REAL(dp), INTENT(INOUT), DIMENSION(0:par_ntor,0:par_mpol1,par_ns,3*par_ntmax) :: a2
1135  a2(:,:,t1lglob:t1rglob,:) = a1(:,:,t1lglob:t1rglob,:)
1136  END SUBROUTINE copy4lastntype
1137  !------------------------------------------------
1138 
1139  !------------------------------------------------
1140  SUBROUTINE copym1lastntype(a1, a2, d1)
1141  REAL(dp), INTENT(IN), DIMENSION(ntmaxblocksize*par_ns) :: a1
1142  REAL(dp), INTENT(INOUT), DIMENSION(ntmaxblocksize*par_ns) :: a2
1143  REAL(dp), INTENT(IN) :: d1
1144  CALL copym4lastntype(a1, a2, d1)
1145  END SUBROUTINE copym1lastntype
1146  !------------------------------------------------
1147 
1148  !------------------------------------------------
1149  SUBROUTINE copym4lastntype(a1, a2, d1)
1150  REAL(dp), INTENT(IN), DIMENSION(0:par_ntor,0:par_mpol1,par_ns,3*par_ntmax) :: a1
1151  REAL(dp), INTENT(INOUT), DIMENSION(0:par_ntor,0:par_mpol1,par_ns,3*par_ntmax) :: a2
1152  REAL(dp), INTENT(IN) :: d1
1153  a2(:,:,t1lglob:t1rglob,:) = d1*a1(:,:,t1lglob:t1rglob,:)
1154  END SUBROUTINE copym4lastntype
1155  !------------------------------------------------
1156 
1157  !------------------------------------------------
1158  SUBROUTINE saxlastns(a, x, vec)
1159  REAL(dp), INTENT(IN), DIMENSION(ntmaxblocksize,par_ns) :: a, x
1160  REAL(dp), INTENT(INOUT), DIMENSION(ntmaxblocksize,par_ns) :: vec
1161  vec(:,tlglob:trglob) = a(:,tlglob:trglob)*x(:,tlglob:trglob)
1162  END SUBROUTINE saxlastns
1163  !------------------------------------------------
1164 
1165  !------------------------------------------------
1166  SUBROUTINE saxlastns1(a, x, vec)
1167  REAL(dp), INTENT(IN), DIMENSION(ntmaxblocksize,tlglob:trglob) :: a
1168  REAL(dp), INTENT(IN), DIMENSION(ntmaxblocksize,par_ns) :: x
1169  REAL(dp), INTENT(INOUT), DIMENSION(ntmaxblocksize,tlglob:trglob) :: vec
1170  vec(:,tlglob:trglob) = a(:,tlglob:trglob)*x(:,tlglob:trglob)
1171  END SUBROUTINE saxlastns1
1172 
1173  !------------------------------------------------
1174  SUBROUTINE saxlastntype(a, x, vec)
1175  REAL(dp), INTENT(IN), DIMENSION(blocksize,par_ns,3*par_ntmax) :: a, x
1176  REAL(dp), INTENT(INOUT), DIMENSION(blocksize,par_ns,3*par_ntmax) :: vec
1177  vec(:,t1lglob:t1rglob,:) = a(:,t1lglob:t1rglob,:)*x(:,t1lglob:t1rglob,:)
1178  END SUBROUTINE saxlastntype
1179  !------------------------------------------------
1180 
1181  !------------------------------------------------
1182  SUBROUTINE saxpbylastntype(a, x, b, y, vec)
1183  REAL(dp), INTENT(IN), DIMENSION(blocksize,par_ns,3*par_ntmax) :: x, y
1184  REAL(dp), INTENT(INOUT), DIMENSION(blocksize,par_ns,3*par_ntmax) :: vec
1185  REAL(dp), INTENT(IN) :: a, b
1186  vec(:,t1lglob:t1rglob,:) = a*x(:,t1lglob:t1rglob,:) + b*y(:,t1lglob:t1rglob,:)
1187  END SUBROUTINE saxpbylastntype
1188  !------------------------------------------------
1189 
1190  !------------------------------------------------
1191  SUBROUTINE saxpbylastns(a, x, b, y, vec)
1192  REAL(dp), INTENT(IN), DIMENSION(ntmaxblocksize,tlglob:trglob) :: x
1193  REAL(dp), INTENT(IN), DIMENSION(ntmaxblocksize,par_ns) :: y
1194  REAL(dp), INTENT(INOUT), DIMENSION(ntmaxblocksize,par_ns) :: vec
1195  REAL(dp), INTENT(IN) :: a, b
1196  vec(:,tlglob:trglob) = a*x(:,tlglob:trglob) + b*y(:,tlglob:trglob)
1197  END SUBROUTINE saxpbylastns
1198  !------------------------------------------------
1199 
1200  !------------------------------------------------
1201  SUBROUTINE saxpby1lastns(a, x, b, y, vec)
1202  REAL(dp), INTENT(IN), DIMENSION(ntmaxblocksize,par_ns) :: x, y
1203  REAL(dp), INTENT(INOUT), DIMENSION(ntmaxblocksize,par_ns) :: vec
1204  REAL(dp), INTENT(IN) :: a, b
1205  vec(:,tlglob:trglob) = a*x(:,tlglob:trglob) + b*y(:,tlglob:trglob)
1206  END SUBROUTINE saxpby1lastns
1207  !------------------------------------------------
1208 
1209  !------------------------------------------------
1210  SUBROUTINE getderivlastns(x1, x2, d1, vec)
1211  REAL(dp), INTENT(IN), DIMENSION(ntmaxblocksize,par_ns) :: x1, x2
1212  REAL(dp), INTENT(INOUT), DIMENSION(ntmaxblocksize,tlglob:trglob) :: vec
1213  REAL(dp), INTENT(IN) :: d1
1214  IF (d1 .eq. 0) RETURN
1215  vec = (x1(:,tlglob:trglob) - x2(:,tlglob:trglob))/d1
1216  END SUBROUTINE getderivlastns
1217  !------------------------------------------------
1218 
1219  !------------------------------------------------
1220  SUBROUTINE zerolastntype(a1)
1221  REAL(dp), INTENT(INOUT), DIMENSION(blocksize,par_ns,3*par_ntmax) :: a1
1222  a1(:,t1lglob:t1rglob,:) = 0
1223  END SUBROUTINE zerolastntype
1224  !------------------------------------------------
1225 
1226  !------------------------------------------------
1227  SUBROUTINE copyparallellinearsubarray(fromarr,toarr, first, last)
1228  INTEGER, INTENT(IN) :: first, last
1229  REAL(dp), INTENT(IN), DIMENSION(blocksize,par_ns,3*par_ntmax) :: fromarr
1230  REAL(dp), INTENT(INOUT), DIMENSION(blocksize,par_ns,3*par_ntmax) :: toarr
1231  toarr(:,first:last,:) = fromarr(:,first:last,:)
1232  END SUBROUTINE copyparallellinearsubarray
1233  !------------------------------------------------
1234 
1235  !------------------------------------------------
1236  SUBROUTINE padsides(arr)
1237  REAL(dp), INTENT(INOUT), DIMENSION(blocksize,par_ns,3*par_ntmax) :: arr
1238  INTEGER :: left, right, tag1=1
1239  REAL(dp) :: ton, toff
1240 
1241  CALL second0(ton)
1242 
1243 #if defined(MPI_OPT)
1244  left = rank - 1; IF(rank.EQ.0) left = mpi_proc_null
1245  right = rank + 1; IF(rank.EQ.nranks-1) right = mpi_proc_null
1246  CALL mpi_sendrecv(arr(:,tlglob,:), ntmaxblocksize, mpi_real8, &
1247  left, tag1, arr(:,t1rglob,:), ntmaxblocksize, mpi_real8, &
1248  right, tag1, ns_comm, mpi_stat,mpi_err)
1249  CALL mpi_sendrecv(arr(:,trglob,:),ntmaxblocksize,mpi_real8, &
1250  right, tag1, arr(:,t1lglob,:), ntmaxblocksize,mpi_real8, &
1251  left, tag1, ns_comm, mpi_stat,mpi_err)
1252 #endif
1253 
1254  CALL second0(toff)
1255  sendrecv_time = sendrecv_time + (toff - ton)
1256 
1257  END SUBROUTINE padsides
1258  !------------------------------------------------
1259 
1260  !------------------------------------------------
1261  SUBROUTINE padsides1x(arr)
1262  REAL(dp), INTENT(INOUT), DIMENSION(par_ns) :: arr
1263  INTEGER :: left, right, tag1
1264  REAL(dp) :: ton, toff
1265 
1266  tag1 = 1
1267 
1268  CALL second0(ton)
1269 
1270 #if defined(MPI_OPT)
1271  left=rank-1; IF(rank.EQ.0) left = mpi_proc_null
1272  right=rank+1; IF(rank.EQ.nranks - 1) right = mpi_proc_null
1273 
1274  IF (grank .LT. nranks) THEN
1275  CALL mpi_sendrecv(arr(tlglob), 1, mpi_real8, left, tag1, &
1276  arr(t1rglob), 1, mpi_real8, right, tag1, ns_comm, &
1277  mpi_stat, mpi_err)
1278  END IF
1279 #endif
1280 
1281  CALL second0(toff)
1282  sendrecv_time = sendrecv_time + (toff - ton)
1283 
1284  END SUBROUTINE padsides1x
1285  !------------------------------------------------
1286 
1287  !------------------------------------------------
1288  SUBROUTINE compareedgevalues(pxc, pxsave)
1289  REAL(dp), INTENT(IN), DIMENSION(blocksize,par_ns,3*par_ntmax) :: pxc, pxsave
1290 
1291  IF (rank .EQ. nranks - 1 .AND. &
1292  any(pxc(:,par_ns,:) .NE. pxsave(:,par_ns,:))) THEN
1293  print *,' xsave != xc at edge returning from GMRES'
1294  END IF
1295 
1296  END SUBROUTINE compareedgevalues
1297  !------------------------------------------------
1298 
1299 END MODULE parallel_vmec_module
1300 !------------------------------------------------
mpi_inc
Umbrella module avoid multiple inlcudes of the mpif.h header.
Definition: mpi_inc.f:11
parallel_vmec_module::copylastntype
Definition: parallel_vmec_module.f90:76
Generated on Thu Mar 5 2020 15:49:27 for V3FIT by   doxygen 1.8.17