~nickpapior/siesta/trunk-buds-format0.92
» Revision
616
: Src/buds/sources_mpi/src/mpi/bud_Dist1D_Array1D.f90
« back to all changes in this revision
Viewing changes to Src/buds/sources_mpi/src/mpi/bud_Dist1D_Array1D.f90
- browse files at revision 616
- compare with another revision
- download diff
- download tarball
- view history from revision 616
- Committer: Nick Papior
- Date: 2017-04-07 12:42:28 UTC
- Revision ID: nickpapior@gmail.com-20170407124228-u5t08yr2p4fhzfeo
Initial commit of buds merged into siesta
Currently I have only enabled buds compilation.
Currently I have only enabled buds compilation.
- files added:
- Src/buds
- Src/buds/.travis
- Src/buds/doc
- Src/buds/include
- Src/buds/setup.makes
- Src/buds/smeka
- Src/buds/smeka/vendor
- Src/buds/sources
- Src/buds/sources/src
- Src/buds/sources/src/mpi
- Src/buds/sources_mpi
- Src/buds/sources_mpi/src
- Src/buds/sources_mpi/src/mpi
- Src/buds/src
- Src/buds/src/mpi
- Src/buds/src/mpi/test
- Src/buds/src/test
- Src/buds/tools
- files modified:
- Src/Makefile
added
removed
Src/buds/sources_mpi/src/mpi/bud_Dist1D_Array1D.f90
1
module bud_iDist1D_b1D
2
use bud_iDist1D
3
use bud_bArray1D
4
use mpi
5
use bud_MP_Comm
6
implicit none
7
private
8
integer, parameter :: rr_ = selected_real_kind(p=6) ! single (real*4)
9
integer, parameter :: rd_ = selected_real_kind(p=15) ! double (real*8)
10
integer, parameter :: is_ = selected_int_kind(4) ! short (integer*2)
11
integer, parameter :: ii_ = selected_int_kind(9) ! int (integer*4)
12
integer, parameter :: il_ = selected_int_kind(18) ! long (integer*8)
13
integer, parameter, private :: BUD_ID_LEn = 36
14
character(len=*), parameter, private :: &
15
BUD_MOd = 'BUD_CC2(BUD_MOD,_)BUD_MOD_NAME'
16
character(len=*), parameter, private :: &
17
BUD_TYPe = "iD1D_b1D"
18
interface dist1d
19
module procedure get_elem1_
20
end interface
21
public :: dist1d
22
interface dist1d_p
23
module procedure get_elem1p_
24
end interface
25
public :: dist1d_p
26
interface array
27
module procedure get_elem2_
28
end interface
29
public :: array
30
interface array_p
31
module procedure get_elem2p_
32
end interface
33
public :: array_p
34
interface dimensions
35
module procedure dimensions_
36
end interface
37
public :: dimensions
38
interface distributed_index
39
module procedure distributed_index_
40
end interface
41
public :: distributed_index
42
interface distribute
43
module procedure distribute_
44
end interface
45
public :: distribute
46
interface new
47
module procedure new_dist_index_
48
end interface new
49
type iD1D_b1D
50
type(iD1D_b1D_), pointer :: D => null()
51
integer :: error_ = 0
52
end type iD1D_b1D
53
type iD1D_b1D_
54
type(iDist1D) :: e1
55
type(bArray1D) :: e2
56
integer :: refs_ = 0
57
character(len=BUD_ID_LEN) :: id_ = "null_id"
58
end type iD1D_b1D_
59
interface new
60
module procedure new_
61
module procedure new_data_
62
end interface
63
public :: new
64
interface assignment(=)
65
module procedure get_elem1_assign_
66
module procedure get_elem2_assign_
67
module procedure set_elem1_
68
module procedure set_elem2_
69
end interface
70
interface element
71
module procedure get_elem1_
72
module procedure get_elem2_
73
end interface
74
interface set_element
75
module procedure set_elem1_
76
module procedure set_elem2_
77
end interface
78
interface element1
79
module procedure get_elem1_
80
end interface
81
public :: element1
82
interface set_element1
83
module procedure set_elem1_
84
end interface
85
public :: set_element1
86
interface element1_p
87
module procedure get_elem1p_
88
end interface
89
public :: element1_p
90
interface element2
91
module procedure get_elem2_
92
end interface
93
public :: element2
94
interface set_element2
95
module procedure set_elem2_
96
end interface
97
public :: set_element2
98
interface element2_p
99
module procedure get_elem2p_
100
end interface
101
public :: element2_p
102
public :: iD1D_b1D
103
private :: iD1D_b1D_
104
interface assignment(=)
105
module procedure common_assign_
106
end interface
107
public :: assignment(=)
108
private :: common_assign_
109
interface initialize
110
module procedure common_initialize_
111
end interface
112
public :: initialize
113
private :: common_initialize_
114
interface is_initialized
115
module procedure common_is_initialized_
116
end interface
117
public :: is_initialized
118
private :: common_is_initialized_
119
interface initialized
120
module procedure common_is_initialized_
121
end interface
122
public :: initialized
123
interface is_initd
124
module procedure common_is_initialized_
125
end interface
126
public :: is_initd
127
interface is_same
128
module procedure common_is_same_
129
end interface
130
public :: is_same
131
private :: common_is_same_
132
interface same
133
module procedure common_is_same_
134
end interface
135
public :: same
136
interface delete
137
module procedure common_delete_
138
end interface
139
public :: delete
140
private :: common_delete_
141
interface nullify
142
module procedure common_nullify_
143
end interface
144
public :: nullify
145
private :: common_nullify_
146
interface copy
147
module procedure copy_
148
end interface
149
public :: copy
150
private :: common_copy_
151
interface print
152
module procedure print_
153
end interface
154
public :: print
155
interface references
156
module procedure common_references_
157
end interface
158
public :: references
159
private :: common_references_
160
interface refs
161
module procedure common_references_
162
end interface
163
public :: refs
164
interface set_error
165
module procedure common_set_error_is_
166
module procedure common_set_error_ii_
167
module procedure common_set_error_il_
168
end interface
169
public :: set_error
170
private :: common_set_error_is_
171
private :: common_set_error_ii_
172
private :: common_set_error_il_
173
interface error
174
module procedure common_error_
175
end interface
176
public :: error
177
private :: common_error_
178
contains
179
subroutine common_copy_(from, to)
180
type(iD1D_b1D), intent(in) :: from
181
type(iD1D_b1D), intent(inout) :: to
182
call set_error(to, error(from))
183
end subroutine common_copy_
184
subroutine common_initialize_(this)
185
type(iD1D_b1D), intent(inout) :: this
186
integer :: error
187
call delete(this)
188
allocate(this%D, stat=error)
189
call set_error(this, error)
190
if ( error /= 0 ) return
191
this%D%refs_ = 1
192
call common_tag_object_(this)
193
end subroutine common_initialize_
194
pure function common_is_initialized_(this) result(init)
195
type(iD1D_b1D), intent(in) :: this
196
logical :: init
197
init = associated(this%D)
198
end function common_is_initialized_
199
elemental function common_is_same_(lhs, rhs) result(same)
200
type(iD1D_b1D), intent(in) :: lhs, rhs
201
logical :: same
202
same = is_initd(lhs) .and. is_initd(rhs)
203
if ( .not. same ) return
204
same = associated(lhs%D, target=rhs%D)
205
end function common_is_same_
206
subroutine common_delete_(this)
207
type(iD1D_b1D), intent(inout) :: this
208
integer :: error
209
call set_error(this, 0)
210
if (.not. is_initd(this) ) return
211
this%D%refs_ = this%D%refs_ - 1
212
if ( 0 == this%D%refs_ ) then
213
call delete_(this)
214
deallocate(this%D, stat=error)
215
call set_error(this, error)
216
end if
217
nullify(this%D)
218
end subroutine common_delete_
219
elemental subroutine common_nullify_(this)
220
type(iD1D_b1D), intent(inout) :: this
221
if (.not. is_initd(this) ) return
222
nullify(this%D)
223
end subroutine common_nullify_
224
subroutine common_assign_(lhs, rhs)
225
type(iD1D_b1D), intent(inout) :: lhs
226
type(iD1D_b1D), intent(in) :: rhs
227
call delete(lhs)
228
if ( .not. is_initd(rhs) ) return
229
lhs%D => rhs%D
230
lhs%D%refs_ = rhs%D%refs_ + 1
231
end subroutine common_assign_
232
elemental function common_references_(this) result(refs)
233
type(iD1D_b1D), intent(in) :: this
234
integer :: refs
235
if ( is_initd(this) ) then
236
refs = this%D%refs_
237
else
238
refs = 0
239
end if
240
end function common_references_
241
elemental function common_error_(this) result(error)
242
type(iD1D_b1D), intent(in) :: this
243
integer :: error
244
if ( is_initd(this) ) then
245
error = this%error_
246
else
247
error = 0
248
end if
249
end function common_error_
250
elemental subroutine common_set_error_is_(this, error)
251
type(iD1D_b1D), intent(inout) :: this
252
integer(is_), intent(in) :: error
253
this%error_ = error
254
end subroutine common_set_error_is_
255
elemental subroutine common_set_error_ii_(this, error)
256
type(iD1D_b1D), intent(inout) :: this
257
integer(ii_), intent(in) :: error
258
this%error_ = error
259
end subroutine common_set_error_ii_
260
elemental subroutine common_set_error_il_(this, error)
261
type(iD1D_b1D), intent(inout) :: this
262
integer(il_), intent(in) :: error
263
this%error_ = error
264
end subroutine common_set_error_il_
265
elemental function common_id_(this) result(str)
266
type(iD1D_b1D), intent(in) :: this
267
character(len=BUD_ID_LEn) :: str
268
str = this%D%id_
269
end function common_id_
270
subroutine common_tag_object_(this)
271
type(iD1D_b1D), intent(inout) :: this
272
end subroutine common_tag_object_
273
subroutine delete_(this)
274
type(iD1D_b1D), intent(inout) :: this
275
call set_error(this, 0)
276
call delete(this%D%e1)
277
if ( 0 /= error(this%D%e1) ) &
278
call set_error(this, error(this%D%e1))
279
call delete(this%D%e2)
280
if ( 0 /= error(this%D%e2) ) &
281
call set_error(this, error(this%D%e2))
282
end subroutine delete_
283
subroutine copy_(from, to)
284
type(iD1D_b1D), intent(in) :: from
285
type(iD1D_b1D), intent(inout) :: to
286
call delete(to)
287
if ( .not. is_initd(from) ) return
288
call initialize(to)
289
call common_copy_(from, to)
290
call copy(from%D%e1, to%D%e1)
291
call copy(from%D%e2, to%D%e2)
292
end subroutine copy_
293
subroutine new_data_(this, a, b &
294
)
295
type(iD1D_b1D), intent(inout) :: this
296
type(iDist1D), intent(inout) :: a
297
type(bArray1D), intent(inout) :: b
298
call new(this)
299
this%D%e1 = a
300
this%D%e2 = b
301
end subroutine new_data_
302
subroutine new_(this)
303
type(iD1D_b1D), intent(inout) :: this
304
call initialize(this)
305
end subroutine new_
306
subroutine get_elem1_(this, item)
307
type(iD1D_b1D), intent(in) :: this
308
type(iDist1D), intent(inout) :: item
309
if ( .not. is_initd(this) ) then
310
call delete(item)
311
else
312
item = this%D% e1
313
end if
314
end subroutine
315
subroutine get_elem1_assign_(item, this)
316
type(iDist1D), intent(inout) :: item
317
type(iD1D_b1D), intent(in) :: this
318
if ( .not. is_initd(this) ) then
319
call delete(item)
320
else
321
item = this%D% e1
322
end if
323
end subroutine
324
subroutine set_elem1_(this, item)
325
type(iD1D_b1D), intent(inout) :: this
326
type(iDist1D), intent(in) :: item
327
if ( .not. is_initd(this) ) return
328
this%D% e1 = item
329
end subroutine
330
function get_elem1p_(this) result(p)
331
type(iD1D_b1D), intent(inout) :: this
332
type(iDist1D), pointer :: p
333
if ( .not. is_initd(this) ) then
334
nullify(p)
335
else
336
p => this%D% e1
337
end if
338
end function
339
subroutine get_elem2_(this, item)
340
type(iD1D_b1D), intent(in) :: this
341
type(bArray1D), intent(inout) :: item
342
if ( .not. is_initd(this) ) then
343
call delete(item)
344
else
345
item = this%D% e2
346
end if
347
end subroutine
348
subroutine get_elem2_assign_(item, this)
349
type(bArray1D), intent(inout) :: item
350
type(iD1D_b1D), intent(in) :: this
351
if ( .not. is_initd(this) ) then
352
call delete(item)
353
else
354
item = this%D% e2
355
end if
356
end subroutine
357
subroutine set_elem2_(this, item)
358
type(iD1D_b1D), intent(inout) :: this
359
type(bArray1D), intent(in) :: item
360
if ( .not. is_initd(this) ) return
361
this%D% e2 = item
362
end subroutine
363
function get_elem2p_(this) result(p)
364
type(iD1D_b1D), intent(inout) :: this
365
type(bArray1D), pointer :: p
366
if ( .not. is_initd(this) ) then
367
nullify(p)
368
else
369
p => this%D% e2
370
end if
371
end function
372
subroutine print_(this, info, indent)
373
type(iD1D_b1D), intent(in) :: this
374
character(len=*), intent(in), optional :: info
375
integer, intent(in), optional :: indent
376
integer :: lindent
377
character(len=32) :: fmt
378
character(len=256) :: name
379
name = "iD1D_b1D"
380
if ( present(info) ) name = info
381
lindent = 1
382
if ( present(indent) ) lindent = indent
383
write(fmt, '(a,i0,a)') '(t',lindent,',3a)'
384
if ( .not. is_initd(this) ) then
385
write(*,fmt) "<", trim(name), " not initialized>"
386
return
387
end if
388
write(fmt, '(a,i0,a)') '(t',lindent,',3a)'
389
lindent = lindent + 2 ! step indentation
390
write(*,fmt) "<<", trim(name), " coll>"
391
call print(this%D%e1, indent = lindent)
392
call print(this%D%e2, indent = lindent)
393
lindent = lindent - 2 ! go back to requested indentation
394
write(fmt, '(a,i0,a)') '(t',lindent,',a,i0,a)'
395
write(*,fmt) " <coll-refs: ", references(this), ">>"
396
end subroutine print_
397
subroutine new_dist_index_(this, dist, arr, dist_idx)
398
type(iD1D_b1D), intent(inout) :: this
399
type(iDist1D), intent(inout) :: dist
400
type(bArray1D), intent(inout) :: arr
401
integer, intent(in) :: dist_idx
402
call new(this, dist, arr)
403
end subroutine new_dist_index_
404
function distributed_index_(this) result(idx)
405
type(iD1D_b1D), intent(in) :: this
406
integer :: idx
407
if ( is_initd(this) ) then
408
idx = 1
409
else
410
idx = -1
411
end if
412
end function distributed_index_
413
pure function dimensions_(this) result(d)
414
type(iD1D_b1D), intent(in) :: this
415
integer :: d
416
if ( is_initd(this) ) then
417
d = 1
418
else
419
d = -1
420
end if
421
end function dimensions_
422
subroutine distribute_(this, parent, out_dist, out)
423
type(iD1D_b1D), intent(inout) :: this
424
type(MP_Comm), intent(inout) :: parent
425
type(iDist1D), intent(inout) :: out_dist
426
type(iD1D_b1D), intent(inout) :: out
427
type(MP_Comm) :: comm, out_comm
428
type(iDist1D) :: fake_dist, dist
429
type(bArray1D) :: arr
430
integer(ii_) :: ir, nranks, my_root, rank, in_rank
431
logical :: is_distr, my_distr
432
integer(ii_) :: dist_idx, dims
433
integer(ii_) :: nl, ng
434
integer(ii_) :: nout
435
integer(ii_), allocatable :: ranks(:)
436
integer(ii_), allocatable :: ashape(:), itmp1(:)
437
if ( .not. is_communicator(parent) ) return
438
rank = comm_rank(parent)
439
nranks = comm_size(parent)
440
dist = this
441
in_rank = comm_rank(dist)
442
if ( in_rank < 0 ) in_rank = in_rank - 1
443
comm = dist
444
arr = this
445
dims = dimensions(arr)
446
call AllReduce_Max(dims, ir, parent)
447
dims = ir
448
allocate(ashape(dims))
449
if ( is_initd(arr) ) then
450
dist_idx = distributed_index(this)
451
do ir = 1 , dims
452
ashape(ir) = size(arr, ir)
453
end do
454
else
455
dist_idx = 0
456
do ir = 1 , dims
457
ashape(ir) = 0
458
end do
459
end if
460
call AllReduce_Max(dist_idx, ir, parent)
461
dist_idx = ir
462
allocate(itmp1(dims))
463
call AllReduce_Max(ashape, itmp1, parent)
464
do ir = 1 , dims
465
if ( ir /= dist_idx ) then
466
ashape(ir) = itmp1(ir)
467
end if
468
end do
469
deallocate(itmp1)
470
out_comm = out_dist
471
if ( is_communicator(out_comm) ) then
472
if ( comm_rank(out_comm) == 0 ) then
473
ir = comm_rank(parent)
474
else
475
ir = -1
476
end if
477
call AllReduce_Max(ir, my_root, out_comm)
478
else
479
my_root = -1
480
end if
481
ng = size_global(dist)
482
call AllReduce_Max(ng, nl, parent)
483
ng = nl
484
nl = size_local(dist)
485
do ir = 0 , nranks - 1
486
if ( my_root == ir ) then
487
my_distr = .true.
488
else
489
my_distr = .false.
490
end if
491
call AllReduce_LOr(my_distr, is_distr, parent)
492
if ( .not. is_distr ) cycle
493
if ( my_distr ) then
494
call new_remote(parent, out_dist)
495
call create_ranks(out_dist)
496
call sub_dist(out_dist)
497
else
498
call new_remote(parent, fake_dist)
499
call create_ranks(fake_dist)
500
call sub_dist(fake_dist)
501
call delete(fake_dist)
502
end if
503
deallocate(ranks)
504
end do
505
call delete(dist)
506
call delete(comm)
507
call delete(arr)
508
call delete(out_comm)
509
contains
510
subroutine create_ranks(dist)
511
type(iDist1D), intent(inout) :: dist
512
type(MP_Comm) :: comm
513
comm = dist
514
call child_Bcast(parent, comm, size=nout)
515
allocate(ranks(-1:nout-1))
516
ranks(-1) = -1
517
call child_Bcast_ranks(parent, comm, nout, ranks(0:))
518
call delete(comm)
519
end subroutine create_ranks
520
subroutine sub_dist(out_dist)
521
use bud_Transfer, only: transfer_dim
522
type(iDist1D), intent(inout) :: out_dist
523
type(MP_Comm) :: out_comm
524
type(bArray1D) :: out_arr
525
logical :: run
526
logical, pointer :: dat (:)
527
logical, pointer :: odat (:)
528
integer(ii_) :: il, ig, oil
529
integer(ii_) :: out_nl
530
integer(ii_) :: i2
531
integer :: send_R, recv_R
532
integer(ii_), allocatable :: reqs(:), stats(:,:)
533
integer(ii_) :: stat(MPI_STATUS_SIZE)
534
run = .true.
535
if ( .not. run ) return
536
out_comm = out_dist
537
out_nl = size_local(out_dist)
538
if ( is_communicator(out_comm) ) then
539
ashape(dist_idx) = out_nl
540
call new(out_arr, ashape(:))
541
call new(out, out_dist, out_arr)
542
odat => array_p(out_arr)
543
call delete(out_arr)
544
end if
545
if ( is_initd(dist) ) then
546
allocate(reqs(nl), stats(MPI_STATUS_SIZE, nl))
547
do il = 1 , nl
548
reqs(il) = MPI_REQUEST_NULL
549
end do
550
end if
551
dat => array_p(arr)
552
if ( dist_idx == dims ) then
553
do ig = 1 , ng
554
send_R = g2rank(dist, ig)
555
recv_R = ranks(g2rank(out_dist, ig))
556
if ( recv_R == rank .and. &
557
send_R == in_rank ) then
558
il = g2l(dist, ig)
559
oil = g2l(out_dist, ig)
560
odat(oil) = dat(il)
561
else if ( recv_R == rank ) then
562
oil = g2l(out_dist, ig)
563
call Recv(odat(oil), MPI_ANY_SOURCE, ig, &
564
parent, stat)
565
else if ( send_R == in_rank ) then
566
il = g2l(dist, ig)
567
call ISSend(dat(il), recv_R, ig, &
568
parent, reqs(il))
569
end if
570
end do
571
if ( allocated(reqs) ) &
572
call WaitAll(nl, reqs, stats, parent)
573
else if ( dist_idx == 1 ) then
574
do i2 = 1 , ashape(2)
575
do ig = 1 , ng
576
send_R = g2rank(dist, ig)
577
recv_R = ranks(g2rank(out_dist, ig))
578
if ( recv_R == rank .and. &
579
send_R == in_rank ) then
580
il = g2l(dist, ig)
581
oil = g2l(out_dist, ig)
582
else if ( recv_R == rank ) then
583
oil = g2l(out_dist, ig)
584
else if ( send_R == in_rank ) then
585
il = g2l(dist, ig)
586
end if
587
end do ! ig
588
if ( allocated(reqs) ) &
589
call WaitAll(nl, reqs, stats, parent)
590
end do ! i2
591
end if
592
if ( allocated(reqs) ) deallocate(reqs)
593
call Barrier(parent)
594
call delete(out_comm)
595
end subroutine sub_dist
596
end subroutine distribute_
597
subroutine write_(f, this)
598
use bud_File
599
use bud_Transfer, only: transfer_dim
600
type( File ), intent(inout) :: f
601
type(iD1D_b1D), intent(in) :: this
602
type(iDist1D) :: dist
603
type(bArray1D) :: arr
604
type( MP_Comm ) :: comm
605
logical :: formatted, do_io
606
integer :: iu, io_rank, dat_rank, rank
607
character(len=64), parameter :: fmt_ = '(e20.16)'
608
logical, pointer :: dat (:)
609
logical, allocatable :: rdat(:)
610
integer :: dist_idx, l, n, ndat
611
integer :: d1, i1
612
integer :: status(MPI_STATUS_SIZE)
613
if ( .not. is_initd(this) ) return
614
dist = this
615
comm = dist
616
arr = this
617
rank = comm_rank(comm)
618
if ( is_open(f) ) then
619
io_rank = rank
620
else
621
io_rank = -1
622
end if
623
call AllReduce_Max(io_rank, iu, comm)
624
io_rank = iu
625
if ( io_rank < 0 ) then
626
call delete(dist)
627
call delete(comm)
628
call delete(arr)
629
return
630
end if
631
do_io = io_rank == rank
632
if ( do_io ) then
633
formatted = is_formatted(f)
634
iu = unit(f)
635
end if
636
d1 = size(arr, 1)
637
dat => array_p(arr)
638
call delete(arr)
639
dist_idx = distributed_index(this)
640
if ( dist_idx == 1 ) then
641
ndat = size_global(dist)
642
else
643
ndat = d1
644
end if
645
if ( do_io ) then
646
allocate(rdat(ndat))
647
if ( formatted ) then
648
write(iu, '(i16)') 1
649
write(iu, '(i16)') d1
650
else
651
write(iu) 1
652
write(iu) d1
653
end if
654
end if
655
if ( dist_idx == 1 ) then
656
i1 = 1
657
do while ( i1 <= d1 )
658
n = consecutive(dist, i1)
659
dat_rank = global2rank(dist, i1)
660
if ( do_io .and. rank == dat_rank ) then
661
l = global2local(dist, i1)
662
call transfer_dim(n, rdat(i1:), n, dat(l:))
663
else if ( do_io ) then
664
call Recv(rdat(i1), n, dat_rank, i1, comm, status)
665
else if ( rank == dat_rank ) then
666
l = global2local(dist, i1)
667
call SSend(dat(l), n, io_rank, i1, comm)
668
end if
669
i1 = i1 + n
670
end do
671
if ( do_io ) then
672
if ( formatted ) then
673
write(iu, fmt_) rdat(:)
674
else
675
write(iu) rdat(:)
676
end if
677
end if
678
end if
679
if ( do_io ) deallocate(rdat)
680
call delete(dist)
681
call delete(comm)
682
end subroutine write_
683
subroutine read_(f, dist, this, dist_idx)
684
use bud_File
685
use bud_Transfer, only: transfer_dim
686
type( File ), intent(inout) :: f
687
type(iDist1D), intent(inout) :: dist
688
type(iD1D_b1D), intent(inout) :: this
689
integer, intent(in), optional :: dist_idx
690
type(bArray1D) :: arr
691
type( MP_Comm ) :: comm
692
logical :: formatted, do_io
693
integer :: iu, io_rank, dat_rank, rank
694
character(len=64), parameter :: fmt_ = '(e20.16)'
695
integer(ii_) :: ashape(1)
696
logical, pointer :: dat (:)
697
logical, allocatable :: rdat(:)
698
integer :: ldist_idx, l, n
699
integer :: d1, i1
700
integer :: status(MPI_STATUS_SIZE)
701
if ( .not. is_initd(this) ) return
702
ldist_idx = 1
703
if ( present(dist_idx) ) ldist_idx = dist_idx
704
comm = dist
705
rank = comm_rank(comm)
706
if ( is_open(f) ) then
707
io_rank = rank
708
else
709
io_rank = -1
710
end if
711
call AllReduce_Max(io_rank, iu, comm)
712
io_rank = iu
713
if ( io_rank < 0 ) then
714
call delete(comm)
715
call delete(this)
716
call set_error(this, -1)
717
return
718
end if
719
do_io = io_rank == rank
720
if ( do_io ) then
721
formatted = is_formatted(f)
722
iu = unit(f)
723
if ( formatted ) then
724
read(iu, '(i16)') l
725
read(iu, '(i16)') ashape
726
d1 = ashape(1)
727
else
728
read(iu) l
729
read(iu) ashape
730
d1 = ashape(1)
731
end if
732
end if
733
select case ( ldist_idx )
734
case ( 1 )
735
ashape(1) = size_local(dist)
736
end select
737
call new(arr, ashape)
738
dat => array_p(arr)
739
call new(this, dist, arr, ldist_idx)
740
call delete(arr)
741
if ( do_io ) then
742
allocate(rdat(d1))
743
end if
744
if ( do_io ) then
745
if ( formatted ) then
746
read(iu, fmt_) rdat(:)
747
else
748
read(iu) rdat(:)
749
end if
750
end if
751
select case ( ldist_idx )
752
case ( 1 )
753
i1 = 1
754
do while ( i1 <= d1 )
755
n = consecutive(dist, i1)
756
dat_rank = global2rank(dist, i1)
757
if ( do_io .and. rank == dat_rank ) then
758
l = global2local(dist, i1)
759
call transfer_dim(n, dat(l:), n, rdat(i1:))
760
else if ( do_io ) then
761
call SSend(rdat(i1), n, dat_rank, i1, comm)
762
else if ( rank == dat_rank ) then
763
l = global2local(dist, i1)
764
call Recv(dat(l), n, io_rank, i1, comm, status)
765
end if
766
i1 = i1 + n
767
end do
768
end select
769
call delete(comm)
770
end subroutine read_
771
end module
772
module bud_iDist1D_r1D
773
use bud_iDist1D
774
use bud_rArray1D
775
use mpi
776
use bud_MP_Comm
777
implicit none
778
private
779
integer, parameter :: rr_ = selected_real_kind(p=6) ! single (real*4)
780
integer, parameter :: rd_ = selected_real_kind(p=15) ! double (real*8)
781
integer, parameter :: is_ = selected_int_kind(4) ! short (integer*2)
782
integer, parameter :: ii_ = selected_int_kind(9) ! int (integer*4)
783
integer, parameter :: il_ = selected_int_kind(18) ! long (integer*8)
784
integer, parameter, private :: BUD_ID_LEn = 36
785
character(len=*), parameter, private :: &
786
BUD_MOd = 'BUD_CC2(BUD_MOD,_)BUD_MOD_NAME'
787
character(len=*), parameter, private :: &
788
BUD_TYPe = "iD1D_r1D"
789
interface dist1d
790
module procedure get_elem1_
791
end interface
792
public :: dist1d
793
interface dist1d_p
794
module procedure get_elem1p_
795
end interface
796
public :: dist1d_p
797
interface array
798
module procedure get_elem2_
799
end interface
800
public :: array
801
interface array_p
802
module procedure get_elem2p_
803
end interface
804
public :: array_p
805
interface dimensions
806
module procedure dimensions_
807
end interface
808
public :: dimensions
809
interface distributed_index
810
module procedure distributed_index_
811
end interface
812
public :: distributed_index
813
interface distribute
814
module procedure distribute_
815
end interface
816
public :: distribute
817
interface new
818
module procedure new_dist_index_
819
end interface new
820
type iD1D_r1D
821
type(iD1D_r1D_), pointer :: D => null()
822
integer :: error_ = 0
823
end type iD1D_r1D
824
type iD1D_r1D_
825
type(iDist1D) :: e1
826
type(rArray1D) :: e2
827
integer :: refs_ = 0
828
character(len=BUD_ID_LEN) :: id_ = "null_id"
829
end type iD1D_r1D_
830
interface new
831
module procedure new_
832
module procedure new_data_
833
end interface
834
public :: new
835
interface assignment(=)
836
module procedure get_elem1_assign_
837
module procedure get_elem2_assign_
838
module procedure set_elem1_
839
module procedure set_elem2_
840
end interface
841
interface element
842
module procedure get_elem1_
843
module procedure get_elem2_
844
end interface
845
interface set_element
846
module procedure set_elem1_
847
module procedure set_elem2_
848
end interface
849
interface element1
850
module procedure get_elem1_
851
end interface
852
public :: element1
853
interface set_element1
854
module procedure set_elem1_
855
end interface
856
public :: set_element1
857
interface element1_p
858
module procedure get_elem1p_
859
end interface
860
public :: element1_p
861
interface element2
862
module procedure get_elem2_
863
end interface
864
public :: element2
865
interface set_element2
866
module procedure set_elem2_
867
end interface
868
public :: set_element2
869
interface element2_p
870
module procedure get_elem2p_
871
end interface
872
public :: element2_p
873
public :: iD1D_r1D
874
private :: iD1D_r1D_
875
interface assignment(=)
876
module procedure common_assign_
877
end interface
878
public :: assignment(=)
879
private :: common_assign_
880
interface initialize
881
module procedure common_initialize_
882
end interface
883
public :: initialize
884
private :: common_initialize_
885
interface is_initialized
886
module procedure common_is_initialized_
887
end interface
888
public :: is_initialized
889
private :: common_is_initialized_
890
interface initialized
891
module procedure common_is_initialized_
892
end interface
893
public :: initialized
894
interface is_initd
895
module procedure common_is_initialized_
896
end interface
897
public :: is_initd
898
interface is_same
899
module procedure common_is_same_
900
end interface
901
public :: is_same
902
private :: common_is_same_
903
interface same
904
module procedure common_is_same_
905
end interface
906
public :: same
907
interface delete
908
module procedure common_delete_
909
end interface
910
public :: delete
911
private :: common_delete_
912
interface nullify
913
module procedure common_nullify_
914
end interface
915
public :: nullify
916
private :: common_nullify_
917
interface copy
918
module procedure copy_
919
end interface
920
public :: copy
921
private :: common_copy_
922
interface print
923
module procedure print_
924
end interface
925
public :: print
926
interface references
927
module procedure common_references_
928
end interface
929
public :: references
930
private :: common_references_
931
interface refs
932
module procedure common_references_
933
end interface
934
public :: refs
935
interface set_error
936
module procedure common_set_error_is_
937
module procedure common_set_error_ii_
938
module procedure common_set_error_il_
939
end interface
940
public :: set_error
941
private :: common_set_error_is_
942
private :: common_set_error_ii_
943
private :: common_set_error_il_
944
interface error
945
module procedure common_error_
946
end interface
947
public :: error
948
private :: common_error_
949
contains
950
subroutine common_copy_(from, to)
951
type(iD1D_r1D), intent(in) :: from
952
type(iD1D_r1D), intent(inout) :: to
953
call set_error(to, error(from))
954
end subroutine common_copy_
955
subroutine common_initialize_(this)
956
type(iD1D_r1D), intent(inout) :: this
957
integer :: error
958
call delete(this)
959
allocate(this%D, stat=error)
960
call set_error(this, error)
961
if ( error /= 0 ) return
962
this%D%refs_ = 1
963
call common_tag_object_(this)
964
end subroutine common_initialize_
965
pure function common_is_initialized_(this) result(init)
966
type(iD1D_r1D), intent(in) :: this
967
logical :: init
968
init = associated(this%D)
969
end function common_is_initialized_
970
elemental function common_is_same_(lhs, rhs) result(same)
971
type(iD1D_r1D), intent(in) :: lhs, rhs
972
logical :: same
973
same = is_initd(lhs) .and. is_initd(rhs)
974
if ( .not. same ) return
975
same = associated(lhs%D, target=rhs%D)
976
end function common_is_same_
977
subroutine common_delete_(this)
978
type(iD1D_r1D), intent(inout) :: this
979
integer :: error
980
call set_error(this, 0)
981
if (.not. is_initd(this) ) return
982
this%D%refs_ = this%D%refs_ - 1
983
if ( 0 == this%D%refs_ ) then
984
call delete_(this)
985
deallocate(this%D, stat=error)
986
call set_error(this, error)
987
end if
988
nullify(this%D)
989
end subroutine common_delete_
990
elemental subroutine common_nullify_(this)
991
type(iD1D_r1D), intent(inout) :: this
992
if (.not. is_initd(this) ) return
993
nullify(this%D)
994
end subroutine common_nullify_
995
subroutine common_assign_(lhs, rhs)
996
type(iD1D_r1D), intent(inout) :: lhs
997
type(iD1D_r1D), intent(in) :: rhs
998
call delete(lhs)
999
if ( .not. is_initd(rhs) ) return
1000
lhs%D => rhs%D
1001
lhs%D%refs_ = rhs%D%refs_ + 1
1002
end subroutine common_assign_
1003
elemental function common_references_(this) result(refs)
1004
type(iD1D_r1D), intent(in) :: this
1005
integer :: refs
1006
if ( is_initd(this) ) then
1007
refs = this%D%refs_
1008
else
1009
refs = 0
1010
end if
1011
end function common_references_
1012
elemental function common_error_(this) result(error)
1013
type(iD1D_r1D), intent(in) :: this
1014
integer :: error
1015
if ( is_initd(this) ) then
1016
error = this%error_
1017
else
1018
error = 0
1019
end if
1020
end function common_error_
1021
elemental subroutine common_set_error_is_(this, error)
1022
type(iD1D_r1D), intent(inout) :: this
1023
integer(is_), intent(in) :: error
1024
this%error_ = error
1025
end subroutine common_set_error_is_
1026
elemental subroutine common_set_error_ii_(this, error)
1027
type(iD1D_r1D), intent(inout) :: this
1028
integer(ii_), intent(in) :: error
1029
this%error_ = error
1030
end subroutine common_set_error_ii_
1031
elemental subroutine common_set_error_il_(this, error)
1032
type(iD1D_r1D), intent(inout) :: this
1033
integer(il_), intent(in) :: error
1034
this%error_ = error
1035
end subroutine common_set_error_il_
1036
elemental function common_id_(this) result(str)
1037
type(iD1D_r1D), intent(in) :: this
1038
character(len=BUD_ID_LEn) :: str
1039
str = this%D%id_
1040
end function common_id_
1041
subroutine common_tag_object_(this)
1042
type(iD1D_r1D), intent(inout) :: this
1043
end subroutine common_tag_object_
1044
subroutine delete_(this)
1045
type(iD1D_r1D), intent(inout) :: this
1046
call set_error(this, 0)
1047
call delete(this%D%e1)
1048
if ( 0 /= error(this%D%e1) ) &
1049
call set_error(this, error(this%D%e1))
1050
call delete(this%D%e2)
1051
if ( 0 /= error(this%D%e2) ) &
1052
call set_error(this, error(this%D%e2))
1053
end subroutine delete_
1054
subroutine copy_(from, to)
1055
type(iD1D_r1D), intent(in) :: from
1056
type(iD1D_r1D), intent(inout) :: to
1057
call delete(to)
1058
if ( .not. is_initd(from) ) return
1059
call initialize(to)
1060
call common_copy_(from, to)
1061
call copy(from%D%e1, to%D%e1)
1062
call copy(from%D%e2, to%D%e2)
1063
end subroutine copy_
1064
subroutine new_data_(this, a, b &
1065
)
1066
type(iD1D_r1D), intent(inout) :: this
1067
type(iDist1D), intent(inout) :: a
1068
type(rArray1D), intent(inout) :: b
1069
call new(this)
1070
this%D%e1 = a
1071
this%D%e2 = b
1072
end subroutine new_data_
1073
subroutine new_(this)
1074
type(iD1D_r1D), intent(inout) :: this
1075
call initialize(this)
1076
end subroutine new_
1077
subroutine get_elem1_(this, item)
1078
type(iD1D_r1D), intent(in) :: this
1079
type(iDist1D), intent(inout) :: item
1080
if ( .not. is_initd(this) ) then
1081
call delete(item)
1082
else
1083
item = this%D% e1
1084
end if
1085
end subroutine
1086
subroutine get_elem1_assign_(item, this)
1087
type(iDist1D), intent(inout) :: item
1088
type(iD1D_r1D), intent(in) :: this
1089
if ( .not. is_initd(this) ) then
1090
call delete(item)
1091
else
1092
item = this%D% e1
1093
end if
1094
end subroutine
1095
subroutine set_elem1_(this, item)
1096
type(iD1D_r1D), intent(inout) :: this
1097
type(iDist1D), intent(in) :: item
1098
if ( .not. is_initd(this) ) return
1099
this%D% e1 = item
1100
end subroutine
1101
function get_elem1p_(this) result(p)
1102
type(iD1D_r1D), intent(inout) :: this
1103
type(iDist1D), pointer :: p
1104
if ( .not. is_initd(this) ) then
1105
nullify(p)
1106
else
1107
p => this%D% e1
1108
end if
1109
end function
1110
subroutine get_elem2_(this, item)
1111
type(iD1D_r1D), intent(in) :: this
1112
type(rArray1D), intent(inout) :: item
1113
if ( .not. is_initd(this) ) then
1114
call delete(item)
1115
else
1116
item = this%D% e2
1117
end if
1118
end subroutine
1119
subroutine get_elem2_assign_(item, this)
1120
type(rArray1D), intent(inout) :: item
1121
type(iD1D_r1D), intent(in) :: this
1122
if ( .not. is_initd(this) ) then
1123
call delete(item)
1124
else
1125
item = this%D% e2
1126
end if
1127
end subroutine
1128
subroutine set_elem2_(this, item)
1129
type(iD1D_r1D), intent(inout) :: this
1130
type(rArray1D), intent(in) :: item
1131
if ( .not. is_initd(this) ) return
1132
this%D% e2 = item
1133
end subroutine
1134
function get_elem2p_(this) result(p)
1135
type(iD1D_r1D), intent(inout) :: this
1136
type(rArray1D), pointer :: p
1137
if ( .not. is_initd(this) ) then
1138
nullify(p)
1139
else
1140
p => this%D% e2
1141
end if
1142
end function
1143
subroutine print_(this, info, indent)
1144
type(iD1D_r1D), intent(in) :: this
1145
character(len=*), intent(in), optional :: info
1146
integer, intent(in), optional :: indent
1147
integer :: lindent
1148
character(len=32) :: fmt
1149
character(len=256) :: name
1150
name = "iD1D_r1D"
1151
if ( present(info) ) name = info
1152
lindent = 1
1153
if ( present(indent) ) lindent = indent
1154
write(fmt, '(a,i0,a)') '(t',lindent,',3a)'
1155
if ( .not. is_initd(this) ) then
1156
write(*,fmt) "<", trim(name), " not initialized>"
1157
return
1158
end if
1159
write(fmt, '(a,i0,a)') '(t',lindent,',3a)'
1160
lindent = lindent + 2 ! step indentation
1161
write(*,fmt) "<<", trim(name), " coll>"
1162
call print(this%D%e1, indent = lindent)
1163
call print(this%D%e2, indent = lindent)
1164
lindent = lindent - 2 ! go back to requested indentation
1165
write(fmt, '(a,i0,a)') '(t',lindent,',a,i0,a)'
1166
write(*,fmt) " <coll-refs: ", references(this), ">>"
1167
end subroutine print_
1168
subroutine new_dist_index_(this, dist, arr, dist_idx)
1169
type(iD1D_r1D), intent(inout) :: this
1170
type(iDist1D), intent(inout) :: dist
1171
type(rArray1D), intent(inout) :: arr
1172
integer, intent(in) :: dist_idx
1173
call new(this, dist, arr)
1174
end subroutine new_dist_index_
1175
function distributed_index_(this) result(idx)
1176
type(iD1D_r1D), intent(in) :: this
1177
integer :: idx
1178
if ( is_initd(this) ) then
1179
idx = 1
1180
else
1181
idx = -1
1182
end if
1183
end function distributed_index_
1184
pure function dimensions_(this) result(d)
1185
type(iD1D_r1D), intent(in) :: this
1186
integer :: d
1187
if ( is_initd(this) ) then
1188
d = 1
1189
else
1190
d = -1
1191
end if
1192
end function dimensions_
1193
subroutine distribute_(this, parent, out_dist, out)
1194
type(iD1D_r1D), intent(inout) :: this
1195
type(MP_Comm), intent(inout) :: parent
1196
type(iDist1D), intent(inout) :: out_dist
1197
type(iD1D_r1D), intent(inout) :: out
1198
type(MP_Comm) :: comm, out_comm
1199
type(iDist1D) :: fake_dist, dist
1200
type(rArray1D) :: arr
1201
integer(ii_) :: ir, nranks, my_root, rank, in_rank
1202
logical :: is_distr, my_distr
1203
integer(ii_) :: dist_idx, dims
1204
integer(ii_) :: nl, ng
1205
integer(ii_) :: nout
1206
integer(ii_), allocatable :: ranks(:)
1207
integer(ii_), allocatable :: ashape(:), itmp1(:)
1208
if ( .not. is_communicator(parent) ) return
1209
rank = comm_rank(parent)
1210
nranks = comm_size(parent)
1211
dist = this
1212
in_rank = comm_rank(dist)
1213
if ( in_rank < 0 ) in_rank = in_rank - 1
1214
comm = dist
1215
arr = this
1216
dims = dimensions(arr)
1217
call AllReduce_Max(dims, ir, parent)
1218
dims = ir
1219
allocate(ashape(dims))
1220
if ( is_initd(arr) ) then
1221
dist_idx = distributed_index(this)
1222
do ir = 1 , dims
1223
ashape(ir) = size(arr, ir)
1224
end do
1225
else
1226
dist_idx = 0
1227
do ir = 1 , dims
1228
ashape(ir) = 0
1229
end do
1230
end if
1231
call AllReduce_Max(dist_idx, ir, parent)
1232
dist_idx = ir
1233
allocate(itmp1(dims))
1234
call AllReduce_Max(ashape, itmp1, parent)
1235
do ir = 1 , dims
1236
if ( ir /= dist_idx ) then
1237
ashape(ir) = itmp1(ir)
1238
end if
1239
end do
1240
deallocate(itmp1)
1241
out_comm = out_dist
1242
if ( is_communicator(out_comm) ) then
1243
if ( comm_rank(out_comm) == 0 ) then
1244
ir = comm_rank(parent)
1245
else
1246
ir = -1
1247
end if
1248
call AllReduce_Max(ir, my_root, out_comm)
1249
else
1250
my_root = -1
1251
end if
1252
ng = size_global(dist)
1253
call AllReduce_Max(ng, nl, parent)
1254
ng = nl
1255
nl = size_local(dist)
1256
do ir = 0 , nranks - 1
1257
if ( my_root == ir ) then
1258
my_distr = .true.
1259
else
1260
my_distr = .false.
1261
end if
1262
call AllReduce_LOr(my_distr, is_distr, parent)
1263
if ( .not. is_distr ) cycle
1264
if ( my_distr ) then
1265
call new_remote(parent, out_dist)
1266
call create_ranks(out_dist)
1267
call sub_dist(out_dist)
1268
else
1269
call new_remote(parent, fake_dist)
1270
call create_ranks(fake_dist)
1271
call sub_dist(fake_dist)
1272
call delete(fake_dist)
1273
end if
1274
deallocate(ranks)
1275
end do
1276
call delete(dist)
1277
call delete(comm)
1278
call delete(arr)
1279
call delete(out_comm)
1280
contains
1281
subroutine create_ranks(dist)
1282
type(iDist1D), intent(inout) :: dist
1283
type(MP_Comm) :: comm
1284
comm = dist
1285
call child_Bcast(parent, comm, size=nout)
1286
allocate(ranks(-1:nout-1))
1287
ranks(-1) = -1
1288
call child_Bcast_ranks(parent, comm, nout, ranks(0:))
1289
call delete(comm)
1290
end subroutine create_ranks
1291
subroutine sub_dist(out_dist)
1292
use bud_Transfer, only: transfer_dim
1293
type(iDist1D), intent(inout) :: out_dist
1294
type(MP_Comm) :: out_comm
1295
type(rArray1D) :: out_arr
1296
logical :: run
1297
real(rr_), pointer :: dat (:)
1298
real(rr_), pointer :: odat (:)
1299
integer(ii_) :: il, ig, oil
1300
integer(ii_) :: out_nl
1301
integer(ii_) :: i2
1302
integer :: send_R, recv_R
1303
integer(ii_), allocatable :: reqs(:), stats(:,:)
1304
integer(ii_) :: stat(MPI_STATUS_SIZE)
1305
run = .true.
1306
if ( .not. run ) return
1307
out_comm = out_dist
1308
out_nl = size_local(out_dist)
1309
if ( is_communicator(out_comm) ) then
1310
ashape(dist_idx) = out_nl
1311
call new(out_arr, ashape(:))
1312
call new(out, out_dist, out_arr)
1313
odat => array_p(out_arr)
1314
call delete(out_arr)
1315
end if
1316
if ( is_initd(dist) ) then
1317
allocate(reqs(nl), stats(MPI_STATUS_SIZE, nl))
1318
do il = 1 , nl
1319
reqs(il) = MPI_REQUEST_NULL
1320
end do
1321
end if
1322
dat => array_p(arr)
1323
if ( dist_idx == dims ) then
1324
do ig = 1 , ng
1325
send_R = g2rank(dist, ig)
1326
recv_R = ranks(g2rank(out_dist, ig))
1327
if ( recv_R == rank .and. &
1328
send_R == in_rank ) then
1329
il = g2l(dist, ig)
1330
oil = g2l(out_dist, ig)
1331
odat(oil) = dat(il)
1332
else if ( recv_R == rank ) then
1333
oil = g2l(out_dist, ig)
1334
call Recv(odat(oil), MPI_ANY_SOURCE, ig, &
1335
parent, stat)
1336
else if ( send_R == in_rank ) then
1337
il = g2l(dist, ig)
1338
call ISSend(dat(il), recv_R, ig, &
1339
parent, reqs(il))
1340
end if
1341
end do
1342
if ( allocated(reqs) ) &
1343
call WaitAll(nl, reqs, stats, parent)
1344
else if ( dist_idx == 1 ) then
1345
do i2 = 1 , ashape(2)
1346
do ig = 1 , ng
1347
send_R = g2rank(dist, ig)
1348
recv_R = ranks(g2rank(out_dist, ig))
1349
if ( recv_R == rank .and. &
1350
send_R == in_rank ) then
1351
il = g2l(dist, ig)
1352
oil = g2l(out_dist, ig)
1353
else if ( recv_R == rank ) then
1354
oil = g2l(out_dist, ig)
1355
else if ( send_R == in_rank ) then
1356
il = g2l(dist, ig)
1357
end if
1358
end do ! ig
1359
if ( allocated(reqs) ) &
1360
call WaitAll(nl, reqs, stats, parent)
1361
end do ! i2
1362
end if
1363
if ( allocated(reqs) ) deallocate(reqs)
1364
call Barrier(parent)
1365
call delete(out_comm)
1366
end subroutine sub_dist
1367
end subroutine distribute_
1368
subroutine write_(f, this)
1369
use bud_File
1370
use bud_Transfer, only: transfer_dim
1371
type( File ), intent(inout) :: f
1372
type(iD1D_r1D), intent(in) :: this
1373
type(iDist1D) :: dist
1374
type(rArray1D) :: arr
1375
type( MP_Comm ) :: comm
1376
logical :: formatted, do_io
1377
integer :: iu, io_rank, dat_rank, rank
1378
character(len=64), parameter :: fmt_ = '(e20.16)'
1379
real(rr_), pointer :: dat (:)
1380
real(rr_), allocatable :: rdat(:)
1381
integer :: dist_idx, l, n, ndat
1382
integer :: d1, i1
1383
integer :: status(MPI_STATUS_SIZE)
1384
if ( .not. is_initd(this) ) return
1385
dist = this
1386
comm = dist
1387
arr = this
1388
rank = comm_rank(comm)
1389
if ( is_open(f) ) then
1390
io_rank = rank
1391
else
1392
io_rank = -1
1393
end if
1394
call AllReduce_Max(io_rank, iu, comm)
1395
io_rank = iu
1396
if ( io_rank < 0 ) then
1397
call delete(dist)
1398
call delete(comm)
1399
call delete(arr)
1400
return
1401
end if
1402
do_io = io_rank == rank
1403
if ( do_io ) then
1404
formatted = is_formatted(f)
1405
iu = unit(f)
1406
end if
1407
d1 = size(arr, 1)
1408
dat => array_p(arr)
1409
call delete(arr)
1410
dist_idx = distributed_index(this)
1411
if ( dist_idx == 1 ) then
1412
ndat = size_global(dist)
1413
else
1414
ndat = d1
1415
end if
1416
if ( do_io ) then
1417
allocate(rdat(ndat))
1418
if ( formatted ) then
1419
write(iu, '(i16)') 1
1420
write(iu, '(i16)') d1
1421
else
1422
write(iu) 1
1423
write(iu) d1
1424
end if
1425
end if
1426
if ( dist_idx == 1 ) then
1427
i1 = 1
1428
do while ( i1 <= d1 )
1429
n = consecutive(dist, i1)
1430
dat_rank = global2rank(dist, i1)
1431
if ( do_io .and. rank == dat_rank ) then
1432
l = global2local(dist, i1)
1433
call transfer_dim(n, rdat(i1:), n, dat(l:))
1434
else if ( do_io ) then
1435
call Recv(rdat(i1), n, dat_rank, i1, comm, status)
1436
else if ( rank == dat_rank ) then
1437
l = global2local(dist, i1)
1438
call SSend(dat(l), n, io_rank, i1, comm)
1439
end if
1440
i1 = i1 + n
1441
end do
1442
if ( do_io ) then
1443
if ( formatted ) then
1444
write(iu, fmt_) rdat(:)
1445
else
1446
write(iu) rdat(:)
1447
end if
1448
end if
1449
end if
1450
if ( do_io ) deallocate(rdat)
1451
call delete(dist)
1452
call delete(comm)
1453
end subroutine write_
1454
subroutine read_(f, dist, this, dist_idx)
1455
use bud_File
1456
use bud_Transfer, only: transfer_dim
1457
type( File ), intent(inout) :: f
1458
type(iDist1D), intent(inout) :: dist
1459
type(iD1D_r1D), intent(inout) :: this
1460
integer, intent(in), optional :: dist_idx
1461
type(rArray1D) :: arr
1462
type( MP_Comm ) :: comm
1463
logical :: formatted, do_io
1464
integer :: iu, io_rank, dat_rank, rank
1465
character(len=64), parameter :: fmt_ = '(e20.16)'
1466
integer(ii_) :: ashape(1)
1467
real(rr_), pointer :: dat (:)
1468
real(rr_), allocatable :: rdat(:)
1469
integer :: ldist_idx, l, n
1470
integer :: d1, i1
1471
integer :: status(MPI_STATUS_SIZE)
1472
if ( .not. is_initd(this) ) return
1473
ldist_idx = 1
1474
if ( present(dist_idx) ) ldist_idx = dist_idx
1475
comm = dist
1476
rank = comm_rank(comm)
1477
if ( is_open(f) ) then
1478
io_rank = rank
1479
else
1480
io_rank = -1
1481
end if
1482
call AllReduce_Max(io_rank, iu, comm)
1483
io_rank = iu
1484
if ( io_rank < 0 ) then
1485
call delete(comm)
1486
call delete(this)
1487
call set_error(this, -1)
1488
return
1489
end if
1490
do_io = io_rank == rank
1491
if ( do_io ) then
1492
formatted = is_formatted(f)
1493
iu = unit(f)
1494
if ( formatted ) then
1495
read(iu, '(i16)') l
1496
read(iu, '(i16)') ashape
1497
d1 = ashape(1)
1498
else
1499
read(iu) l
1500
read(iu) ashape
1501
d1 = ashape(1)
1502
end if
1503
end if
1504
select case ( ldist_idx )
1505
case ( 1 )
1506
ashape(1) = size_local(dist)
1507
end select
1508
call new(arr, ashape)
1509
dat => array_p(arr)
1510
call new(this, dist, arr, ldist_idx)
1511
call delete(arr)
1512
if ( do_io ) then
1513
allocate(rdat(d1))
1514
end if
1515
if ( do_io ) then
1516
if ( formatted ) then
1517
read(iu, fmt_) rdat(:)
1518
else
1519
read(iu) rdat(:)
1520
end if
1521
end if
1522
select case ( ldist_idx )
1523
case ( 1 )
1524
i1 = 1
1525
do while ( i1 <= d1 )
1526
n = consecutive(dist, i1)
1527
dat_rank = global2rank(dist, i1)
1528
if ( do_io .and. rank == dat_rank ) then
1529
l = global2local(dist, i1)
1530
call transfer_dim(n, dat(l:), n, rdat(i1:))
1531
else if ( do_io ) then
1532
call SSend(rdat(i1), n, dat_rank, i1, comm)
1533
else if ( rank == dat_rank ) then
1534
l = global2local(dist, i1)
1535
call Recv(dat(l), n, io_rank, i1, comm, status)
1536
end if
1537
i1 = i1 + n
1538
end do
1539
end select
1540
call delete(comm)
1541
end subroutine read_
1542
end module
1543
module bud_iDist1D_d1D
1544
use bud_iDist1D
1545
use bud_dArray1D
1546
use mpi
1547
use bud_MP_Comm
1548
implicit none
1549
private
1550
integer, parameter :: rr_ = selected_real_kind(p=6) ! single (real*4)
1551
integer, parameter :: rd_ = selected_real_kind(p=15) ! double (real*8)
1552
integer, parameter :: is_ = selected_int_kind(4) ! short (integer*2)
1553
integer, parameter :: ii_ = selected_int_kind(9) ! int (integer*4)
1554
integer, parameter :: il_ = selected_int_kind(18) ! long (integer*8)
1555
integer, parameter, private :: BUD_ID_LEn = 36
1556
character(len=*), parameter, private :: &
1557
BUD_MOd = 'BUD_CC2(BUD_MOD,_)BUD_MOD_NAME'
1558
character(len=*), parameter, private :: &
1559
BUD_TYPe = "iD1D_d1D"
1560
interface dist1d
1561
module procedure get_elem1_
1562
end interface
1563
public :: dist1d
1564
interface dist1d_p
1565
module procedure get_elem1p_
1566
end interface
1567
public :: dist1d_p
1568
interface array
1569
module procedure get_elem2_
1570
end interface
1571
public :: array
1572
interface array_p
1573
module procedure get_elem2p_
1574
end interface
1575
public :: array_p
1576
interface dimensions
1577
module procedure dimensions_
1578
end interface
1579
public :: dimensions
1580
interface distributed_index
1581
module procedure distributed_index_
1582
end interface
1583
public :: distributed_index
1584
interface distribute
1585
module procedure distribute_
1586
end interface
1587
public :: distribute
1588
interface new
1589
module procedure new_dist_index_
1590
end interface new
1591
type iD1D_d1D
1592
type(iD1D_d1D_), pointer :: D => null()
1593
integer :: error_ = 0
1594
end type iD1D_d1D
1595
type iD1D_d1D_
1596
type(iDist1D) :: e1
1597
type(dArray1D) :: e2
1598
integer :: refs_ = 0
1599
character(len=BUD_ID_LEN) :: id_ = "null_id"
1600
end type iD1D_d1D_
1601
interface new
1602
module procedure new_
1603
module procedure new_data_
1604
end interface
1605
public :: new
1606
interface assignment(=)
1607
module procedure get_elem1_assign_
1608
module procedure get_elem2_assign_
1609
module procedure set_elem1_
1610
module procedure set_elem2_
1611
end interface
1612
interface element
1613
module procedure get_elem1_
1614
module procedure get_elem2_
1615
end interface
1616
interface set_element
1617
module procedure set_elem1_
1618
module procedure set_elem2_
1619
end interface
1620
interface element1
1621
module procedure get_elem1_
1622
end interface
1623
public :: element1
1624
interface set_element1
1625
module procedure set_elem1_
1626
end interface
1627
public :: set_element1
1628
interface element1_p
1629
module procedure get_elem1p_
1630
end interface
1631
public :: element1_p
1632
interface element2
1633
module procedure get_elem2_
1634
end interface
1635
public :: element2
1636
interface set_element2
1637
module procedure set_elem2_
1638
end interface
1639
public :: set_element2
1640
interface element2_p
1641
module procedure get_elem2p_
1642
end interface
1643
public :: element2_p
1644
public :: iD1D_d1D
1645
private :: iD1D_d1D_
1646
interface assignment(=)
1647
module procedure common_assign_
1648
end interface
1649
public :: assignment(=)
1650
private :: common_assign_
1651
interface initialize
1652
module procedure common_initialize_
1653
end interface
1654
public :: initialize
1655
private :: common_initialize_
1656
interface is_initialized
1657
module procedure common_is_initialized_
1658
end interface
1659
public :: is_initialized
1660
private :: common_is_initialized_
1661
interface initialized
1662
module procedure common_is_initialized_
1663
end interface
1664
public :: initialized
1665
interface is_initd
1666
module procedure common_is_initialized_
1667
end interface
1668
public :: is_initd
1669
interface is_same
1670
module procedure common_is_same_
1671
end interface
1672
public :: is_same
1673
private :: common_is_same_
1674
interface same
1675
module procedure common_is_same_
1676
end interface
1677
public :: same
1678
interface delete
1679
module procedure common_delete_
1680
end interface
1681
public :: delete
1682
private :: common_delete_
1683
interface nullify
1684
module procedure common_nullify_
1685
end interface
1686
public :: nullify
1687
private :: common_nullify_
1688
interface copy
1689
module procedure copy_
1690
end interface
1691
public :: copy
1692
private :: common_copy_
1693
interface print
1694
module procedure print_
1695
end interface
1696
public :: print
1697
interface references
1698
module procedure common_references_
1699
end interface
1700
public :: references
1701
private :: common_references_
1702
interface refs
1703
module procedure common_references_
1704
end interface
1705
public :: refs
1706
interface set_error
1707
module procedure common_set_error_is_
1708
module procedure common_set_error_ii_
1709
module procedure common_set_error_il_
1710
end interface
1711
public :: set_error
1712
private :: common_set_error_is_
1713
private :: common_set_error_ii_
1714
private :: common_set_error_il_
1715
interface error
1716
module procedure common_error_
1717
end interface
1718
public :: error
1719
private :: common_error_
1720
contains
1721
subroutine common_copy_(from, to)
1722
type(iD1D_d1D), intent(in) :: from
1723
type(iD1D_d1D), intent(inout) :: to
1724
call set_error(to, error(from))
1725
end subroutine common_copy_
1726
subroutine common_initialize_(this)
1727
type(iD1D_d1D), intent(inout) :: this
1728
integer :: error
1729
call delete(this)
1730
allocate(this%D, stat=error)
1731
call set_error(this, error)
1732
if ( error /= 0 ) return
1733
this%D%refs_ = 1
1734
call common_tag_object_(this)
1735
end subroutine common_initialize_
1736
pure function common_is_initialized_(this) result(init)
1737
type(iD1D_d1D), intent(in) :: this
1738
logical :: init
1739
init = associated(this%D)
1740
end function common_is_initialized_
1741
elemental function common_is_same_(lhs, rhs) result(same)
1742
type(iD1D_d1D), intent(in) :: lhs, rhs
1743
logical :: same
1744
same = is_initd(lhs) .and. is_initd(rhs)
1745
if ( .not. same ) return
1746
same = associated(lhs%D, target=rhs%D)
1747
end function common_is_same_
1748
subroutine common_delete_(this)
1749
type(iD1D_d1D), intent(inout) :: this
1750
integer :: error
1751
call set_error(this, 0)
1752
if (.not. is_initd(this) ) return
1753
this%D%refs_ = this%D%refs_ - 1
1754
if ( 0 == this%D%refs_ ) then
1755
call delete_(this)
1756
deallocate(this%D, stat=error)
1757
call set_error(this, error)
1758
end if
1759
nullify(this%D)
1760
end subroutine common_delete_
1761
elemental subroutine common_nullify_(this)
1762
type(iD1D_d1D), intent(inout) :: this
1763
if (.not. is_initd(this) ) return
1764
nullify(this%D)
1765
end subroutine common_nullify_
1766
subroutine common_assign_(lhs, rhs)
1767
type(iD1D_d1D), intent(inout) :: lhs
1768
type(iD1D_d1D), intent(in) :: rhs
1769
call delete(lhs)
1770
if ( .not. is_initd(rhs) ) return
1771
lhs%D => rhs%D
1772
lhs%D%refs_ = rhs%D%refs_ + 1
1773
end subroutine common_assign_
1774
elemental function common_references_(this) result(refs)
1775
type(iD1D_d1D), intent(in) :: this
1776
integer :: refs
1777
if ( is_initd(this) ) then
1778
refs = this%D%refs_
1779
else
1780
refs = 0
1781
end if
1782
end function common_references_
1783
elemental function common_error_(this) result(error)
1784
type(iD1D_d1D), intent(in) :: this
1785
integer :: error
1786
if ( is_initd(this) ) then
1787
error = this%error_
1788
else
1789
error = 0
1790
end if
1791
end function common_error_
1792
elemental subroutine common_set_error_is_(this, error)
1793
type(iD1D_d1D), intent(inout) :: this
1794
integer(is_), intent(in) :: error
1795
this%error_ = error
1796
end subroutine common_set_error_is_
1797
elemental subroutine common_set_error_ii_(this, error)
1798
type(iD1D_d1D), intent(inout) :: this
1799
integer(ii_), intent(in) :: error
1800
this%error_ = error
1801
end subroutine common_set_error_ii_
1802
elemental subroutine common_set_error_il_(this, error)
1803
type(iD1D_d1D), intent(inout) :: this
1804
integer(il_), intent(in) :: error
1805
this%error_ = error
1806
end subroutine common_set_error_il_
1807
elemental function common_id_(this) result(str)
1808
type(iD1D_d1D), intent(in) :: this
1809
character(len=BUD_ID_LEn) :: str
1810
str = this%D%id_
1811
end function common_id_
1812
subroutine common_tag_object_(this)
1813
type(iD1D_d1D), intent(inout) :: this
1814
end subroutine common_tag_object_
1815
subroutine delete_(this)
1816
type(iD1D_d1D), intent(inout) :: this
1817
call set_error(this, 0)
1818
call delete(this%D%e1)
1819
if ( 0 /= error(this%D%e1) ) &
1820
call set_error(this, error(this%D%e1))
1821
call delete(this%D%e2)
1822
if ( 0 /= error(this%D%e2) ) &
1823
call set_error(this, error(this%D%e2))
1824
end subroutine delete_
1825
subroutine copy_(from, to)
1826
type(iD1D_d1D), intent(in) :: from
1827
type(iD1D_d1D), intent(inout) :: to
1828
call delete(to)
1829
if ( .not. is_initd(from) ) return
1830
call initialize(to)
1831
call common_copy_(from, to)
1832
call copy(from%D%e1, to%D%e1)
1833
call copy(from%D%e2, to%D%e2)
1834
end subroutine copy_
1835
subroutine new_data_(this, a, b &
1836
)
1837
type(iD1D_d1D), intent(inout) :: this
1838
type(iDist1D), intent(inout) :: a
1839
type(dArray1D), intent(inout) :: b
1840
call new(this)
1841
this%D%e1 = a
1842
this%D%e2 = b
1843
end subroutine new_data_
1844
subroutine new_(this)
1845
type(iD1D_d1D), intent(inout) :: this
1846
call initialize(this)
1847
end subroutine new_
1848
subroutine get_elem1_(this, item)
1849
type(iD1D_d1D), intent(in) :: this
1850
type(iDist1D), intent(inout) :: item
1851
if ( .not. is_initd(this) ) then
1852
call delete(item)
1853
else
1854
item = this%D% e1
1855
end if
1856
end subroutine
1857
subroutine get_elem1_assign_(item, this)
1858
type(iDist1D), intent(inout) :: item
1859
type(iD1D_d1D), intent(in) :: this
1860
if ( .not. is_initd(this) ) then
1861
call delete(item)
1862
else
1863
item = this%D% e1
1864
end if
1865
end subroutine
1866
subroutine set_elem1_(this, item)
1867
type(iD1D_d1D), intent(inout) :: this
1868
type(iDist1D), intent(in) :: item
1869
if ( .not. is_initd(this) ) return
1870
this%D% e1 = item
1871
end subroutine
1872
function get_elem1p_(this) result(p)
1873
type(iD1D_d1D), intent(inout) :: this
1874
type(iDist1D), pointer :: p
1875
if ( .not. is_initd(this) ) then
1876
nullify(p)
1877
else
1878
p => this%D% e1
1879
end if
1880
end function
1881
subroutine get_elem2_(this, item)
1882
type(iD1D_d1D), intent(in) :: this
1883
type(dArray1D), intent(inout) :: item
1884
if ( .not. is_initd(this) ) then
1885
call delete(item)
1886
else
1887
item = this%D% e2
1888
end if
1889
end subroutine
1890
subroutine get_elem2_assign_(item, this)
1891
type(dArray1D), intent(inout) :: item
1892
type(iD1D_d1D), intent(in) :: this
1893
if ( .not. is_initd(this) ) then
1894
call delete(item)
1895
else
1896
item = this%D% e2
1897
end if
1898
end subroutine
1899
subroutine set_elem2_(this, item)
1900
type(iD1D_d1D), intent(inout) :: this
1901
type(dArray1D), intent(in) :: item
1902
if ( .not. is_initd(this) ) return
1903
this%D% e2 = item
1904
end subroutine
1905
function get_elem2p_(this) result(p)
1906
type(iD1D_d1D), intent(inout) :: this
1907
type(dArray1D), pointer :: p
1908
if ( .not. is_initd(this) ) then
1909
nullify(p)
1910
else
1911
p => this%D% e2
1912
end if
1913
end function
1914
subroutine print_(this, info, indent)
1915
type(iD1D_d1D), intent(in) :: this
1916
character(len=*), intent(in), optional :: info
1917
integer, intent(in), optional :: indent
1918
integer :: lindent
1919
character(len=32) :: fmt
1920
character(len=256) :: name
1921
name = "iD1D_d1D"
1922
if ( present(info) ) name = info
1923
lindent = 1
1924
if ( present(indent) ) lindent = indent
1925
write(fmt, '(a,i0,a)') '(t',lindent,',3a)'
1926
if ( .not. is_initd(this) ) then
1927
write(*,fmt) "<", trim(name), " not initialized>"
1928
return
1929
end if
1930
write(fmt, '(a,i0,a)') '(t',lindent,',3a)'
1931
lindent = lindent + 2 ! step indentation
1932
write(*,fmt) "<<", trim(name), " coll>"
1933
call print(this%D%e1, indent = lindent)
1934
call print(this%D%e2, indent = lindent)
1935
lindent = lindent - 2 ! go back to requested indentation
1936
write(fmt, '(a,i0,a)') '(t',lindent,',a,i0,a)'
1937
write(*,fmt) " <coll-refs: ", references(this), ">>"
1938
end subroutine print_
1939
subroutine new_dist_index_(this, dist, arr, dist_idx)
1940
type(iD1D_d1D), intent(inout) :: this
1941
type(iDist1D), intent(inout) :: dist
1942
type(dArray1D), intent(inout) :: arr
1943
integer, intent(in) :: dist_idx
1944
call new(this, dist, arr)
1945
end subroutine new_dist_index_
1946
function distributed_index_(this) result(idx)
1947
type(iD1D_d1D), intent(in) :: this
1948
integer :: idx
1949
if ( is_initd(this) ) then
1950
idx = 1
1951
else
1952
idx = -1
1953
end if
1954
end function distributed_index_
1955
pure function dimensions_(this) result(d)
1956
type(iD1D_d1D), intent(in) :: this
1957
integer :: d
1958
if ( is_initd(this) ) then
1959
d = 1
1960
else
1961
d = -1
1962
end if
1963
end function dimensions_
1964
subroutine distribute_(this, parent, out_dist, out)
1965
type(iD1D_d1D), intent(inout) :: this
1966
type(MP_Comm), intent(inout) :: parent
1967
type(iDist1D), intent(inout) :: out_dist
1968
type(iD1D_d1D), intent(inout) :: out
1969
type(MP_Comm) :: comm, out_comm
1970
type(iDist1D) :: fake_dist, dist
1971
type(dArray1D) :: arr
1972
integer(ii_) :: ir, nranks, my_root, rank, in_rank
1973
logical :: is_distr, my_distr
1974
integer(ii_) :: dist_idx, dims
1975
integer(ii_) :: nl, ng
1976
integer(ii_) :: nout
1977
integer(ii_), allocatable :: ranks(:)
1978
integer(ii_), allocatable :: ashape(:), itmp1(:)
1979
if ( .not. is_communicator(parent) ) return
1980
rank = comm_rank(parent)
1981
nranks = comm_size(parent)
1982
dist = this
1983
in_rank = comm_rank(dist)
1984
if ( in_rank < 0 ) in_rank = in_rank - 1
1985
comm = dist
1986
arr = this
1987
dims = dimensions(arr)
1988
call AllReduce_Max(dims, ir, parent)
1989
dims = ir
1990
allocate(ashape(dims))
1991
if ( is_initd(arr) ) then
1992
dist_idx = distributed_index(this)
1993
do ir = 1 , dims
1994
ashape(ir) = size(arr, ir)
1995
end do
1996
else
1997
dist_idx = 0
1998
do ir = 1 , dims
1999
ashape(ir) = 0
2000
end do
2001
end if
2002
call AllReduce_Max(dist_idx, ir, parent)
2003
dist_idx = ir
2004
allocate(itmp1(dims))
2005
call AllReduce_Max(ashape, itmp1, parent)
2006
do ir = 1 , dims
2007
if ( ir /= dist_idx ) then
2008
ashape(ir) = itmp1(ir)
2009
end if
2010
end do
2011
deallocate(itmp1)
2012
out_comm = out_dist
2013
if ( is_communicator(out_comm) ) then
2014
if ( comm_rank(out_comm) == 0 ) then
2015
ir = comm_rank(parent)
2016
else
2017
ir = -1
2018
end if
2019
call AllReduce_Max(ir, my_root, out_comm)
2020
else
2021
my_root = -1
2022
end if
2023
ng = size_global(dist)
2024
call AllReduce_Max(ng, nl, parent)
2025
ng = nl
2026
nl = size_local(dist)
2027
do ir = 0 , nranks - 1
2028
if ( my_root == ir ) then
2029
my_distr = .true.
2030
else
2031
my_distr = .false.
2032
end if
2033
call AllReduce_LOr(my_distr, is_distr, parent)
2034
if ( .not. is_distr ) cycle
2035
if ( my_distr ) then
2036
call new_remote(parent, out_dist)
2037
call create_ranks(out_dist)
2038
call sub_dist(out_dist)
2039
else
2040
call new_remote(parent, fake_dist)
2041
call create_ranks(fake_dist)
2042
call sub_dist(fake_dist)
2043
call delete(fake_dist)
2044
end if
2045
deallocate(ranks)
2046
end do
2047
call delete(dist)
2048
call delete(comm)
2049
call delete(arr)
2050
call delete(out_comm)
2051
contains
2052
subroutine create_ranks(dist)
2053
type(iDist1D), intent(inout) :: dist
2054
type(MP_Comm) :: comm
2055
comm = dist
2056
call child_Bcast(parent, comm, size=nout)
2057
allocate(ranks(-1:nout-1))
2058
ranks(-1) = -1
2059
call child_Bcast_ranks(parent, comm, nout, ranks(0:))
2060
call delete(comm)
2061
end subroutine create_ranks
2062
subroutine sub_dist(out_dist)
2063
use bud_Transfer, only: transfer_dim
2064
type(iDist1D), intent(inout) :: out_dist
2065
type(MP_Comm) :: out_comm
2066
type(dArray1D) :: out_arr
2067
logical :: run
2068
real(rd_), pointer :: dat (:)
2069
real(rd_), pointer :: odat (:)
2070
integer(ii_) :: il, ig, oil
2071
integer(ii_) :: out_nl
2072
integer(ii_) :: i2
2073
integer :: send_R, recv_R
2074
integer(ii_), allocatable :: reqs(:), stats(:,:)
2075
integer(ii_) :: stat(MPI_STATUS_SIZE)
2076
run = .true.
2077
if ( .not. run ) return
2078
out_comm = out_dist
2079
out_nl = size_local(out_dist)
2080
if ( is_communicator(out_comm) ) then
2081
ashape(dist_idx) = out_nl
2082
call new(out_arr, ashape(:))
2083
call new(out, out_dist, out_arr)
2084
odat => array_p(out_arr)
2085
call delete(out_arr)
2086
end if
2087
if ( is_initd(dist) ) then
2088
allocate(reqs(nl), stats(MPI_STATUS_SIZE, nl))
2089
do il = 1 , nl
2090
reqs(il) = MPI_REQUEST_NULL
2091
end do
2092
end if
2093
dat => array_p(arr)
2094
if ( dist_idx == dims ) then
2095
do ig = 1 , ng
2096
send_R = g2rank(dist, ig)
2097
recv_R = ranks(g2rank(out_dist, ig))
2098
if ( recv_R == rank .and. &
2099
send_R == in_rank ) then
2100
il = g2l(dist, ig)
2101
oil = g2l(out_dist, ig)
2102
odat(oil) = dat(il)
2103
else if ( recv_R == rank ) then
2104
oil = g2l(out_dist, ig)
2105
call Recv(odat(oil), MPI_ANY_SOURCE, ig, &
2106
parent, stat)
2107
else if ( send_R == in_rank ) then
2108
il = g2l(dist, ig)
2109
call ISSend(dat(il), recv_R, ig, &
2110
parent, reqs(il))
2111
end if
2112
end do
2113
if ( allocated(reqs) ) &
2114
call WaitAll(nl, reqs, stats, parent)
2115
else if ( dist_idx == 1 ) then
2116
do i2 = 1 , ashape(2)
2117
do ig = 1 , ng
2118
send_R = g2rank(dist, ig)
2119
recv_R = ranks(g2rank(out_dist, ig))
2120
if ( recv_R == rank .and. &
2121
send_R == in_rank ) then
2122
il = g2l(dist, ig)
2123
oil = g2l(out_dist, ig)
2124
else if ( recv_R == rank ) then
2125
oil = g2l(out_dist, ig)
2126
else if ( send_R == in_rank ) then
2127
il = g2l(dist, ig)
2128
end if
2129
end do ! ig
2130
if ( allocated(reqs) ) &
2131
call WaitAll(nl, reqs, stats, parent)
2132
end do ! i2
2133
end if
2134
if ( allocated(reqs) ) deallocate(reqs)
2135
call Barrier(parent)
2136
call delete(out_comm)
2137
end subroutine sub_dist
2138
end subroutine distribute_
2139
subroutine write_(f, this)
2140
use bud_File
2141
use bud_Transfer, only: transfer_dim
2142
type( File ), intent(inout) :: f
2143
type(iD1D_d1D), intent(in) :: this
2144
type(iDist1D) :: dist
2145
type(dArray1D) :: arr
2146
type( MP_Comm ) :: comm
2147
logical :: formatted, do_io
2148
integer :: iu, io_rank, dat_rank, rank
2149
character(len=64), parameter :: fmt_ = '(e20.16)'
2150
real(rd_), pointer :: dat (:)
2151
real(rd_), allocatable :: rdat(:)
2152
integer :: dist_idx, l, n, ndat
2153
integer :: d1, i1
2154
integer :: status(MPI_STATUS_SIZE)
2155
if ( .not. is_initd(this) ) return
2156
dist = this
2157
comm = dist
2158
arr = this
2159
rank = comm_rank(comm)
2160
if ( is_open(f) ) then
2161
io_rank = rank
2162
else
2163
io_rank = -1
2164
end if
2165
call AllReduce_Max(io_rank, iu, comm)
2166
io_rank = iu
2167
if ( io_rank < 0 ) then
2168
call delete(dist)
2169
call delete(comm)
2170
call delete(arr)
2171
return
2172
end if
2173
do_io = io_rank == rank
2174
if ( do_io ) then
2175
formatted = is_formatted(f)
2176
iu = unit(f)
2177
end if
2178
d1 = size(arr, 1)
2179
dat => array_p(arr)
2180
call delete(arr)
2181
dist_idx = distributed_index(this)
2182
if ( dist_idx == 1 ) then
2183
ndat = size_global(dist)
2184
else
2185
ndat = d1
2186
end if
2187
if ( do_io ) then
2188
allocate(rdat(ndat))
2189
if ( formatted ) then
2190
write(iu, '(i16)') 1
2191
write(iu, '(i16)') d1
2192
else
2193
write(iu) 1
2194
write(iu) d1
2195
end if
2196
end if
2197
if ( dist_idx == 1 ) then
2198
i1 = 1
2199
do while ( i1 <= d1 )
2200
n = consecutive(dist, i1)
2201
dat_rank = global2rank(dist, i1)
2202
if ( do_io .and. rank == dat_rank ) then
2203
l = global2local(dist, i1)
2204
call transfer_dim(n, rdat(i1:), n, dat(l:))
2205
else if ( do_io ) then
2206
call Recv(rdat(i1), n, dat_rank, i1, comm, status)
2207
else if ( rank == dat_rank ) then
2208
l = global2local(dist, i1)
2209
call SSend(dat(l), n, io_rank, i1, comm)
2210
end if
2211
i1 = i1 + n
2212
end do
2213
if ( do_io ) then
2214
if ( formatted ) then
2215
write(iu, fmt_) rdat(:)
2216
else
2217
write(iu) rdat(:)
2218
end if
2219
end if
2220
end if
2221
if ( do_io ) deallocate(rdat)
2222
call delete(dist)
2223
call delete(comm)
2224
end subroutine write_
2225
subroutine read_(f, dist, this, dist_idx)
2226
use bud_File
2227
use bud_Transfer, only: transfer_dim
2228
type( File ), intent(inout) :: f
2229
type(iDist1D), intent(inout) :: dist
2230
type(iD1D_d1D), intent(inout) :: this
2231
integer, intent(in), optional :: dist_idx
2232
type(dArray1D) :: arr
2233
type( MP_Comm ) :: comm
2234
logical :: formatted, do_io
2235
integer :: iu, io_rank, dat_rank, rank
2236
character(len=64), parameter :: fmt_ = '(e20.16)'
2237
integer(ii_) :: ashape(1)
2238
real(rd_), pointer :: dat (:)
2239
real(rd_), allocatable :: rdat(:)
2240
integer :: ldist_idx, l, n
2241
integer :: d1, i1
2242
integer :: status(MPI_STATUS_SIZE)
2243
if ( .not. is_initd(this) ) return
2244
ldist_idx = 1
2245
if ( present(dist_idx) ) ldist_idx = dist_idx
2246
comm = dist
2247
rank = comm_rank(comm)
2248
if ( is_open(f) ) then
2249
io_rank = rank
2250
else
2251
io_rank = -1
2252
end if
2253
call AllReduce_Max(io_rank, iu, comm)
2254
io_rank = iu
2255
if ( io_rank < 0 ) then
2256
call delete(comm)
2257
call delete(this)
2258
call set_error(this, -1)
2259
return
2260
end if
2261
do_io = io_rank == rank
2262
if ( do_io ) then
2263
formatted = is_formatted(f)
2264
iu = unit(f)
2265
if ( formatted ) then
2266
read(iu, '(i16)') l
2267
read(iu, '(i16)') ashape
2268
d1 = ashape(1)
2269
else
2270
read(iu) l
2271
read(iu) ashape
2272
d1 = ashape(1)
2273
end if
2274
end if
2275
select case ( ldist_idx )
2276
case ( 1 )
2277
ashape(1) = size_local(dist)
2278
end select
2279
call new(arr, ashape)
2280
dat => array_p(arr)
2281
call new(this, dist, arr, ldist_idx)
2282
call delete(arr)
2283
if ( do_io ) then
2284
allocate(rdat(d1))
2285
end if
2286
if ( do_io ) then
2287
if ( formatted ) then
2288
read(iu, fmt_) rdat(:)
2289
else
2290
read(iu) rdat(:)
2291
end if
2292
end if
2293
select case ( ldist_idx )
2294
case ( 1 )
2295
i1 = 1
2296
do while ( i1 <= d1 )
2297
n = consecutive(dist, i1)
2298
dat_rank = global2rank(dist, i1)
2299
if ( do_io .and. rank == dat_rank ) then
2300
l = global2local(dist, i1)
2301
call transfer_dim(n, dat(l:), n, rdat(i1:))
2302
else if ( do_io ) then
2303
call SSend(rdat(i1), n, dat_rank, i1, comm)
2304
else if ( rank == dat_rank ) then
2305
l = global2local(dist, i1)
2306
call Recv(dat(l), n, io_rank, i1, comm, status)
2307
end if
2308
i1 = i1 + n
2309
end do
2310
end select
2311
call delete(comm)
2312
end subroutine read_
2313
end module
2314
module bud_iDist1D_c1D
2315
use bud_iDist1D
2316
use bud_cArray1D
2317
use mpi
2318
use bud_MP_Comm
2319
implicit none
2320
private
2321
integer, parameter :: rr_ = selected_real_kind(p=6) ! single (real*4)
2322
integer, parameter :: rd_ = selected_real_kind(p=15) ! double (real*8)
2323
integer, parameter :: is_ = selected_int_kind(4) ! short (integer*2)
2324
integer, parameter :: ii_ = selected_int_kind(9) ! int (integer*4)
2325
integer, parameter :: il_ = selected_int_kind(18) ! long (integer*8)
2326
integer, parameter, private :: BUD_ID_LEn = 36
2327
character(len=*), parameter, private :: &
2328
BUD_MOd = 'BUD_CC2(BUD_MOD,_)BUD_MOD_NAME'
2329
character(len=*), parameter, private :: &
2330
BUD_TYPe = "iD1D_c1D"
2331
interface dist1d
2332
module procedure get_elem1_
2333
end interface
2334
public :: dist1d
2335
interface dist1d_p
2336
module procedure get_elem1p_
2337
end interface
2338
public :: dist1d_p
2339
interface array
2340
module procedure get_elem2_
2341
end interface
2342
public :: array
2343
interface array_p
2344
module procedure get_elem2p_
2345
end interface
2346
public :: array_p
2347
interface dimensions
2348
module procedure dimensions_
2349
end interface
2350
public :: dimensions
2351
interface distributed_index
2352
module procedure distributed_index_
2353
end interface
2354
public :: distributed_index
2355
interface distribute
2356
module procedure distribute_
2357
end interface
2358
public :: distribute
2359
interface new
2360
module procedure new_dist_index_
2361
end interface new
2362
type iD1D_c1D
2363
type(iD1D_c1D_), pointer :: D => null()
2364
integer :: error_ = 0
2365
end type iD1D_c1D
2366
type iD1D_c1D_
2367
type(iDist1D) :: e1
2368
type(cArray1D) :: e2
2369
integer :: refs_ = 0
2370
character(len=BUD_ID_LEN) :: id_ = "null_id"
2371
end type iD1D_c1D_
2372
interface new
2373
module procedure new_
2374
module procedure new_data_
2375
end interface
2376
public :: new
2377
interface assignment(=)
2378
module procedure get_elem1_assign_
2379
module procedure get_elem2_assign_
2380
module procedure set_elem1_
2381
module procedure set_elem2_
2382
end interface
2383
interface element
2384
module procedure get_elem1_
2385
module procedure get_elem2_
2386
end interface
2387
interface set_element
2388
module procedure set_elem1_
2389
module procedure set_elem2_
2390
end interface
2391
interface element1
2392
module procedure get_elem1_
2393
end interface
2394
public :: element1
2395
interface set_element1
2396
module procedure set_elem1_
2397
end interface
2398
public :: set_element1
2399
interface element1_p
2400
module procedure get_elem1p_
2401
end interface
2402
public :: element1_p
2403
interface element2
2404
module procedure get_elem2_
2405
end interface
2406
public :: element2
2407
interface set_element2
2408
module procedure set_elem2_
2409
end interface
2410
public :: set_element2
2411
interface element2_p
2412
module procedure get_elem2p_
2413
end interface
2414
public :: element2_p
2415
public :: iD1D_c1D
2416
private :: iD1D_c1D_
2417
interface assignment(=)
2418
module procedure common_assign_
2419
end interface
2420
public :: assignment(=)
2421
private :: common_assign_
2422
interface initialize
2423
module procedure common_initialize_
2424
end interface
2425
public :: initialize
2426
private :: common_initialize_
2427
interface is_initialized
2428
module procedure common_is_initialized_
2429
end interface
2430
public :: is_initialized
2431
private :: common_is_initialized_
2432
interface initialized
2433
module procedure common_is_initialized_
2434
end interface
2435
public :: initialized
2436
interface is_initd
2437
module procedure common_is_initialized_
2438
end interface
2439
public :: is_initd
2440
interface is_same
2441
module procedure common_is_same_
2442
end interface
2443
public :: is_same
2444
private :: common_is_same_
2445
interface same
2446
module procedure common_is_same_
2447
end interface
2448
public :: same
2449
interface delete
2450
module procedure common_delete_
2451
end interface
2452
public :: delete
2453
private :: common_delete_
2454
interface nullify
2455
module procedure common_nullify_
2456
end interface
2457
public :: nullify
2458
private :: common_nullify_
2459
interface copy
2460
module procedure copy_
2461
end interface
2462
public :: copy
2463
private :: common_copy_
2464
interface print
2465
module procedure print_
2466
end interface
2467
public :: print
2468
interface references
2469
module procedure common_references_
2470
end interface
2471
public :: references
2472
private :: common_references_
2473
interface refs
2474
module procedure common_references_
2475
end interface
2476
public :: refs
2477
interface set_error
2478
module procedure common_set_error_is_
2479
module procedure common_set_error_ii_
2480
module procedure common_set_error_il_
2481
end interface
2482
public :: set_error
2483
private :: common_set_error_is_
2484
private :: common_set_error_ii_
2485
private :: common_set_error_il_
2486
interface error
2487
module procedure common_error_
2488
end interface
2489
public :: error
2490
private :: common_error_
2491
contains
2492
subroutine common_copy_(from, to)
2493
type(iD1D_c1D), intent(in) :: from
2494
type(iD1D_c1D), intent(inout) :: to
2495
call set_error(to, error(from))
2496
end subroutine common_copy_
2497
subroutine common_initialize_(this)
2498
type(iD1D_c1D), intent(inout) :: this
2499
integer :: error
2500
call delete(this)
2501
allocate(this%D, stat=error)
2502
call set_error(this, error)
2503
if ( error /= 0 ) return
2504
this%D%refs_ = 1
2505
call common_tag_object_(this)
2506
end subroutine common_initialize_
2507
pure function common_is_initialized_(this) result(init)
2508
type(iD1D_c1D), intent(in) :: this
2509
logical :: init
2510
init = associated(this%D)
2511
end function common_is_initialized_
2512
elemental function common_is_same_(lhs, rhs) result(same)
2513
type(iD1D_c1D), intent(in) :: lhs, rhs
2514
logical :: same
2515
same = is_initd(lhs) .and. is_initd(rhs)
2516
if ( .not. same ) return
2517
same = associated(lhs%D, target=rhs%D)
2518
end function common_is_same_
2519
subroutine common_delete_(this)
2520
type(iD1D_c1D), intent(inout) :: this
2521
integer :: error
2522
call set_error(this, 0)
2523
if (.not. is_initd(this) ) return
2524
this%D%refs_ = this%D%refs_ - 1
2525
if ( 0 == this%D%refs_ ) then
2526
call delete_(this)
2527
deallocate(this%D, stat=error)
2528
call set_error(this, error)
2529
end if
2530
nullify(this%D)
2531
end subroutine common_delete_
2532
elemental subroutine common_nullify_(this)
2533
type(iD1D_c1D), intent(inout) :: this
2534
if (.not. is_initd(this) ) return
2535
nullify(this%D)
2536
end subroutine common_nullify_
2537
subroutine common_assign_(lhs, rhs)
2538
type(iD1D_c1D), intent(inout) :: lhs
2539
type(iD1D_c1D), intent(in) :: rhs
2540
call delete(lhs)
2541
if ( .not. is_initd(rhs) ) return
2542
lhs%D => rhs%D
2543
lhs%D%refs_ = rhs%D%refs_ + 1
2544
end subroutine common_assign_
2545
elemental function common_references_(this) result(refs)
2546
type(iD1D_c1D), intent(in) :: this
2547
integer :: refs
2548
if ( is_initd(this) ) then
2549
refs = this%D%refs_
2550
else
2551
refs = 0
2552
end if
2553
end function common_references_
2554
elemental function common_error_(this) result(error)
2555
type(iD1D_c1D), intent(in) :: this
2556
integer :: error
2557
if ( is_initd(this) ) then
2558
error = this%error_
2559
else
2560
error = 0
2561
end if
2562
end function common_error_
2563
elemental subroutine common_set_error_is_(this, error)
2564
type(iD1D_c1D), intent(inout) :: this
2565
integer(is_), intent(in) :: error
2566
this%error_ = error
2567
end subroutine common_set_error_is_
2568
elemental subroutine common_set_error_ii_(this, error)
2569
type(iD1D_c1D), intent(inout) :: this
2570
integer(ii_), intent(in) :: error
2571
this%error_ = error
2572
end subroutine common_set_error_ii_
2573
elemental subroutine common_set_error_il_(this, error)
2574
type(iD1D_c1D), intent(inout) :: this
2575
integer(il_), intent(in) :: error
2576
this%error_ = error
2577
end subroutine common_set_error_il_
2578
elemental function common_id_(this) result(str)
2579
type(iD1D_c1D), intent(in) :: this
2580
character(len=BUD_ID_LEn) :: str
2581
str = this%D%id_
2582
end function common_id_
2583
subroutine common_tag_object_(this)
2584
type(iD1D_c1D), intent(inout) :: this
2585
end subroutine common_tag_object_
2586
subroutine delete_(this)
2587
type(iD1D_c1D), intent(inout) :: this
2588
call set_error(this, 0)
2589
call delete(this%D%e1)
2590
if ( 0 /= error(this%D%e1) ) &
2591
call set_error(this, error(this%D%e1))
2592
call delete(this%D%e2)
2593
if ( 0 /= error(this%D%e2) ) &
2594
call set_error(this, error(this%D%e2))
2595
end subroutine delete_
2596
subroutine copy_(from, to)
2597
type(iD1D_c1D), intent(in) :: from
2598
type(iD1D_c1D), intent(inout) :: to
2599
call delete(to)
2600
if ( .not. is_initd(from) ) return
2601
call initialize(to)
2602
call common_copy_(from, to)
2603
call copy(from%D%e1, to%D%e1)
2604
call copy(from%D%e2, to%D%e2)
2605
end subroutine copy_
2606
subroutine new_data_(this, a, b &
2607
)
2608
type(iD1D_c1D), intent(inout) :: this
2609
type(iDist1D), intent(inout) :: a
2610
type(cArray1D), intent(inout) :: b
2611
call new(this)
2612
this%D%e1 = a
2613
this%D%e2 = b
2614
end subroutine new_data_
2615
subroutine new_(this)
2616
type(iD1D_c1D), intent(inout) :: this
2617
call initialize(this)
2618
end subroutine new_
2619
subroutine get_elem1_(this, item)
2620
type(iD1D_c1D), intent(in) :: this
2621
type(iDist1D), intent(inout) :: item
2622
if ( .not. is_initd(this) ) then
2623
call delete(item)
2624
else
2625
item = this%D% e1
2626
end if
2627
end subroutine
2628
subroutine get_elem1_assign_(item, this)
2629
type(iDist1D), intent(inout) :: item
2630
type(iD1D_c1D), intent(in) :: this
2631
if ( .not. is_initd(this) ) then
2632
call delete(item)
2633
else
2634
item = this%D% e1
2635
end if
2636
end subroutine
2637
subroutine set_elem1_(this, item)
2638
type(iD1D_c1D), intent(inout) :: this
2639
type(iDist1D), intent(in) :: item
2640
if ( .not. is_initd(this) ) return
2641
this%D% e1 = item
2642
end subroutine
2643
function get_elem1p_(this) result(p)
2644
type(iD1D_c1D), intent(inout) :: this
2645
type(iDist1D), pointer :: p
2646
if ( .not. is_initd(this) ) then
2647
nullify(p)
2648
else
2649
p => this%D% e1
2650
end if
2651
end function
2652
subroutine get_elem2_(this, item)
2653
type(iD1D_c1D), intent(in) :: this
2654
type(cArray1D), intent(inout) :: item
2655
if ( .not. is_initd(this) ) then
2656
call delete(item)
2657
else
2658
item = this%D% e2
2659
end if
2660
end subroutine
2661
subroutine get_elem2_assign_(item, this)
2662
type(cArray1D), intent(inout) :: item
2663
type(iD1D_c1D), intent(in) :: this
2664
if ( .not. is_initd(this) ) then
2665
call delete(item)
2666
else
2667
item = this%D% e2
2668
end if
2669
end subroutine
2670
subroutine set_elem2_(this, item)
2671
type(iD1D_c1D), intent(inout) :: this
2672
type(cArray1D), intent(in) :: item
2673
if ( .not. is_initd(this) ) return
2674
this%D% e2 = item
2675
end subroutine
2676
function get_elem2p_(this) result(p)
2677
type(iD1D_c1D), intent(inout) :: this
2678
type(cArray1D), pointer :: p
2679
if ( .not. is_initd(this) ) then
2680
nullify(p)
2681
else
2682
p => this%D% e2
2683
end if
2684
end function
2685
subroutine print_(this, info, indent)
2686
type(iD1D_c1D), intent(in) :: this
2687
character(len=*), intent(in), optional :: info
2688
integer, intent(in), optional :: indent
2689
integer :: lindent
2690
character(len=32) :: fmt
2691
character(len=256) :: name
2692
name = "iD1D_c1D"
2693
if ( present(info) ) name = info
2694
lindent = 1
2695
if ( present(indent) ) lindent = indent
2696
write(fmt, '(a,i0,a)') '(t',lindent,',3a)'
2697
if ( .not. is_initd(this) ) then
2698
write(*,fmt) "<", trim(name), " not initialized>"
2699
return
2700
end if
2701
write(fmt, '(a,i0,a)') '(t',lindent,',3a)'
2702
lindent = lindent + 2 ! step indentation
2703
write(*,fmt) "<<", trim(name), " coll>"
2704
call print(this%D%e1, indent = lindent)
2705
call print(this%D%e2, indent = lindent)
2706
lindent = lindent - 2 ! go back to requested indentation
2707
write(fmt, '(a,i0,a)') '(t',lindent,',a,i0,a)'
2708
write(*,fmt) " <coll-refs: ", references(this), ">>"
2709
end subroutine print_
2710
subroutine new_dist_index_(this, dist, arr, dist_idx)
2711
type(iD1D_c1D), intent(inout) :: this
2712
type(iDist1D), intent(inout) :: dist
2713
type(cArray1D), intent(inout) :: arr
2714
integer, intent(in) :: dist_idx
2715
call new(this, dist, arr)
2716
end subroutine new_dist_index_
2717
function distributed_index_(this) result(idx)
2718
type(iD1D_c1D), intent(in) :: this
2719
integer :: idx
2720
if ( is_initd(this) ) then
2721
idx = 1
2722
else
2723
idx = -1
2724
end if
2725
end function distributed_index_
2726
pure function dimensions_(this) result(d)
2727
type(iD1D_c1D), intent(in) :: this
2728
integer :: d
2729
if ( is_initd(this) ) then
2730
d = 1
2731
else
2732
d = -1
2733
end if
2734
end function dimensions_
2735
subroutine distribute_(this, parent, out_dist, out)
2736
type(iD1D_c1D), intent(inout) :: this
2737
type(MP_Comm), intent(inout) :: parent
2738
type(iDist1D), intent(inout) :: out_dist
2739
type(iD1D_c1D), intent(inout) :: out
2740
type(MP_Comm) :: comm, out_comm
2741
type(iDist1D) :: fake_dist, dist
2742
type(cArray1D) :: arr
2743
integer(ii_) :: ir, nranks, my_root, rank, in_rank
2744
logical :: is_distr, my_distr
2745
integer(ii_) :: dist_idx, dims
2746
integer(ii_) :: nl, ng
2747
integer(ii_) :: nout
2748
integer(ii_), allocatable :: ranks(:)
2749
integer(ii_), allocatable :: ashape(:), itmp1(:)
2750
if ( .not. is_communicator(parent) ) return
2751
rank = comm_rank(parent)
2752
nranks = comm_size(parent)
2753
dist = this
2754
in_rank = comm_rank(dist)
2755
if ( in_rank < 0 ) in_rank = in_rank - 1
2756
comm = dist
2757
arr = this
2758
dims = dimensions(arr)
2759
call AllReduce_Max(dims, ir, parent)
2760
dims = ir
2761
allocate(ashape(dims))
2762
if ( is_initd(arr) ) then
2763
dist_idx = distributed_index(this)
2764
do ir = 1 , dims
2765
ashape(ir) = size(arr, ir)
2766
end do
2767
else
2768
dist_idx = 0
2769
do ir = 1 , dims
2770
ashape(ir) = 0
2771
end do
2772
end if
2773
call AllReduce_Max(dist_idx, ir, parent)
2774
dist_idx = ir
2775
allocate(itmp1(dims))
2776
call AllReduce_Max(ashape, itmp1, parent)
2777
do ir = 1 , dims
2778
if ( ir /= dist_idx ) then
2779
ashape(ir) = itmp1(ir)
2780
end if
2781
end do
2782
deallocate(itmp1)
2783
out_comm = out_dist
2784
if ( is_communicator(out_comm) ) then
2785
if ( comm_rank(out_comm) == 0 ) then
2786
ir = comm_rank(parent)
2787
else
2788
ir = -1
2789
end if
2790
call AllReduce_Max(ir, my_root, out_comm)
2791
else
2792
my_root = -1
2793
end if
2794
ng = size_global(dist)
2795
call AllReduce_Max(ng, nl, parent)
2796
ng = nl
2797
nl = size_local(dist)
2798
do ir = 0 , nranks - 1
2799
if ( my_root == ir ) then
2800
my_distr = .true.
2801
else
2802
my_distr = .false.
2803
end if
2804
call AllReduce_LOr(my_distr, is_distr, parent)
2805
if ( .not. is_distr ) cycle
2806
if ( my_distr ) then
2807
call new_remote(parent, out_dist)
2808
call create_ranks(out_dist)
2809
call sub_dist(out_dist)
2810
else
2811
call new_remote(parent, fake_dist)
2812
call create_ranks(fake_dist)
2813
call sub_dist(fake_dist)
2814
call delete(fake_dist)
2815
end if
2816
deallocate(ranks)
2817
end do
2818
call delete(dist)
2819
call delete(comm)
2820
call delete(arr)
2821
call delete(out_comm)
2822
contains
2823
subroutine create_ranks(dist)
2824
type(iDist1D), intent(inout) :: dist
2825
type(MP_Comm) :: comm
2826
comm = dist
2827
call child_Bcast(parent, comm, size=nout)
2828
allocate(ranks(-1:nout-1))
2829
ranks(-1) = -1
2830
call child_Bcast_ranks(parent, comm, nout, ranks(0:))
2831
call delete(comm)
2832
end subroutine create_ranks
2833
subroutine sub_dist(out_dist)
2834
use bud_Transfer, only: transfer_dim
2835
type(iDist1D), intent(inout) :: out_dist
2836
type(MP_Comm) :: out_comm
2837
type(cArray1D) :: out_arr
2838
logical :: run
2839
complex(rr_), pointer :: dat (:)
2840
complex(rr_), pointer :: odat (:)
2841
integer(ii_) :: il, ig, oil
2842
integer(ii_) :: out_nl
2843
integer(ii_) :: i2
2844
integer :: send_R, recv_R
2845
integer(ii_), allocatable :: reqs(:), stats(:,:)
2846
integer(ii_) :: stat(MPI_STATUS_SIZE)
2847
run = .true.
2848
if ( .not. run ) return
2849
out_comm = out_dist
2850
out_nl = size_local(out_dist)
2851
if ( is_communicator(out_comm) ) then
2852
ashape(dist_idx) = out_nl
2853
call new(out_arr, ashape(:))
2854
call new(out, out_dist, out_arr)
2855
odat => array_p(out_arr)
2856
call delete(out_arr)
2857
end if
2858
if ( is_initd(dist) ) then
2859
allocate(reqs(nl), stats(MPI_STATUS_SIZE, nl))
2860
do il = 1 , nl
2861
reqs(il) = MPI_REQUEST_NULL
2862
end do
2863
end if
2864
dat => array_p(arr)
2865
if ( dist_idx == dims ) then
2866
do ig = 1 , ng
2867
send_R = g2rank(dist, ig)
2868
recv_R = ranks(g2rank(out_dist, ig))
2869
if ( recv_R == rank .and. &
2870
send_R == in_rank ) then
2871
il = g2l(dist, ig)
2872
oil = g2l(out_dist, ig)
2873
odat(oil) = dat(il)
2874
else if ( recv_R == rank ) then
2875
oil = g2l(out_dist, ig)
2876
call Recv(odat(oil), MPI_ANY_SOURCE, ig, &
2877
parent, stat)
2878
else if ( send_R == in_rank ) then
2879
il = g2l(dist, ig)
2880
call ISSend(dat(il), recv_R, ig, &
2881
parent, reqs(il))
2882
end if
2883
end do
2884
if ( allocated(reqs) ) &
2885
call WaitAll(nl, reqs, stats, parent)
2886
else if ( dist_idx == 1 ) then
2887
do i2 = 1 , ashape(2)
2888
do ig = 1 , ng
2889
send_R = g2rank(dist, ig)
2890
recv_R = ranks(g2rank(out_dist, ig))
2891
if ( recv_R == rank .and. &
2892
send_R == in_rank ) then
2893
il = g2l(dist, ig)
2894
oil = g2l(out_dist, ig)
2895
else if ( recv_R == rank ) then
2896
oil = g2l(out_dist, ig)
2897
else if ( send_R == in_rank ) then
2898
il = g2l(dist, ig)
2899
end if
2900
end do ! ig
2901
if ( allocated(reqs) ) &
2902
call WaitAll(nl, reqs, stats, parent)
2903
end do ! i2
2904
end if
2905
if ( allocated(reqs) ) deallocate(reqs)
2906
call Barrier(parent)
2907
call delete(out_comm)
2908
end subroutine sub_dist
2909
end subroutine distribute_
2910
subroutine write_(f, this)
2911
use bud_File
2912
use bud_Transfer, only: transfer_dim
2913
type( File ), intent(inout) :: f
2914
type(iD1D_c1D), intent(in) :: this
2915
type(iDist1D) :: dist
2916
type(cArray1D) :: arr
2917
type( MP_Comm ) :: comm
2918
logical :: formatted, do_io
2919
integer :: iu, io_rank, dat_rank, rank
2920
character(len=64), parameter :: fmt_ = '(e20.16)'
2921
complex(rr_), pointer :: dat (:)
2922
complex(rr_), allocatable :: rdat(:)
2923
integer :: dist_idx, l, n, ndat
2924
integer :: d1, i1
2925
integer :: status(MPI_STATUS_SIZE)
2926
if ( .not. is_initd(this) ) return
2927
dist = this
2928
comm = dist
2929
arr = this
2930
rank = comm_rank(comm)
2931
if ( is_open(f) ) then
2932
io_rank = rank
2933
else
2934
io_rank = -1
2935
end if
2936
call AllReduce_Max(io_rank, iu, comm)
2937
io_rank = iu
2938
if ( io_rank < 0 ) then
2939
call delete(dist)
2940
call delete(comm)
2941
call delete(arr)
2942
return
2943
end if
2944
do_io = io_rank == rank
2945
if ( do_io ) then
2946
formatted = is_formatted(f)
2947
iu = unit(f)
2948
end if
2949
d1 = size(arr, 1)
2950
dat => array_p(arr)
2951
call delete(arr)
2952
dist_idx = distributed_index(this)
2953
if ( dist_idx == 1 ) then
2954
ndat = size_global(dist)
2955
else
2956
ndat = d1
2957
end if
2958
if ( do_io ) then
2959
allocate(rdat(ndat))
2960
if ( formatted ) then
2961
write(iu, '(i16)') 1
2962
write(iu, '(i16)') d1
2963
else
2964
write(iu) 1
2965
write(iu) d1
2966
end if
2967
end if
2968
if ( dist_idx == 1 ) then
2969
i1 = 1
2970
do while ( i1 <= d1 )
2971
n = consecutive(dist, i1)
2972
dat_rank = global2rank(dist, i1)
2973
if ( do_io .and. rank == dat_rank ) then
2974
l = global2local(dist, i1)
2975
call transfer_dim(n, rdat(i1:), n, dat(l:))
2976
else if ( do_io ) then
2977
call Recv(rdat(i1), n, dat_rank, i1, comm, status)
2978
else if ( rank == dat_rank ) then
2979
l = global2local(dist, i1)
2980
call SSend(dat(l), n, io_rank, i1, comm)
2981
end if
2982
i1 = i1 + n
2983
end do
2984
if ( do_io ) then
2985
if ( formatted ) then
2986
write(iu, fmt_) rdat(:)
2987
else
2988
write(iu) rdat(:)
2989
end if
2990
end if
2991
end if
2992
if ( do_io ) deallocate(rdat)
2993
call delete(dist)
2994
call delete(comm)
2995
end subroutine write_
2996
subroutine read_(f, dist, this, dist_idx)
2997
use bud_File
2998
use bud_Transfer, only: transfer_dim
2999
type( File ), intent(inout) :: f
3000
type(iDist1D), intent(inout) :: dist
3001
type(iD1D_c1D), intent(inout) :: this
3002
integer, intent(in), optional :: dist_idx
3003
type(cArray1D) :: arr
3004
type( MP_Comm ) :: comm
3005
logical :: formatted, do_io
3006
integer :: iu, io_rank, dat_rank, rank
3007
character(len=64), parameter :: fmt_ = '(e20.16)'
3008
integer(ii_) :: ashape(1)
3009
complex(rr_), pointer :: dat (:)
3010
complex(rr_), allocatable :: rdat(:)
3011
integer :: ldist_idx, l, n
3012
integer :: d1, i1
3013
integer :: status(MPI_STATUS_SIZE)
3014
if ( .not. is_initd(this) ) return
3015
ldist_idx = 1
3016
if ( present(dist_idx) ) ldist_idx = dist_idx
3017
comm = dist
3018
rank = comm_rank(comm)
3019
if ( is_open(f) ) then
3020
io_rank = rank
3021
else
3022
io_rank = -1
3023
end if
3024
call AllReduce_Max(io_rank, iu, comm)
3025
io_rank = iu
3026
if ( io_rank < 0 ) then
3027
call delete(comm)
3028
call delete(this)
3029
call set_error(this, -1)
3030
return
3031
end if
3032
do_io = io_rank == rank
3033
if ( do_io ) then
3034
formatted = is_formatted(f)
3035
iu = unit(f)
3036
if ( formatted ) then
3037
read(iu, '(i16)') l
3038
read(iu, '(i16)') ashape
3039
d1 = ashape(1)
3040
else
3041
read(iu) l
3042
read(iu) ashape
3043
d1 = ashape(1)
3044
end if
3045
end if
3046
select case ( ldist_idx )
3047
case ( 1 )
3048
ashape(1) = size_local(dist)
3049
end select
3050
call new(arr, ashape)
3051
dat => array_p(arr)
3052
call new(this, dist, arr, ldist_idx)
3053
call delete(arr)
3054
if ( do_io ) then
3055
allocate(rdat(d1))
3056
end if
3057
if ( do_io ) then
3058
if ( formatted ) then
3059
read(iu, fmt_) rdat(:)
3060
else
3061
read(iu) rdat(:)
3062
end if
3063
end if
3064
select case ( ldist_idx )
3065
case ( 1 )
3066
i1 = 1
3067
do while ( i1 <= d1 )
3068
n = consecutive(dist, i1)
3069
dat_rank = global2rank(dist, i1)
3070
if ( do_io .and. rank == dat_rank ) then
3071
l = global2local(dist, i1)
3072
call transfer_dim(n, dat(l:), n, rdat(i1:))
3073
else if ( do_io ) then
3074
call SSend(rdat(i1), n, dat_rank, i1, comm)
3075
else if ( rank == dat_rank ) then
3076
l = global2local(dist, i1)
3077
call Recv(dat(l), n, io_rank, i1, comm, status)
3078
end if
3079
i1 = i1 + n
3080
end do
3081
end select
3082
call delete(comm)
3083
end subroutine read_
3084
end module
3085
module bud_iDist1D_z1D
3086
use bud_iDist1D
3087
use bud_zArray1D
3088
use mpi
3089
use bud_MP_Comm
3090
implicit none
3091
private
3092
integer, parameter :: rr_ = selected_real_kind(p=6) ! single (real*4)
3093
integer, parameter :: rd_ = selected_real_kind(p=15) ! double (real*8)
3094
integer, parameter :: is_ = selected_int_kind(4) ! short (integer*2)
3095
integer, parameter :: ii_ = selected_int_kind(9) ! int (integer*4)
3096
integer, parameter :: il_ = selected_int_kind(18) ! long (integer*8)
3097
integer, parameter, private :: BUD_ID_LEn = 36
3098
character(len=*), parameter, private :: &
3099
BUD_MOd = 'BUD_CC2(BUD_MOD,_)BUD_MOD_NAME'
3100
character(len=*), parameter, private :: &
3101
BUD_TYPe = "iD1D_z1D"
3102
interface dist1d
3103
module procedure get_elem1_
3104
end interface
3105
public :: dist1d
3106
interface dist1d_p
3107
module procedure get_elem1p_
3108
end interface
3109
public :: dist1d_p
3110
interface array
3111
module procedure get_elem2_
3112
end interface
3113
public :: array
3114
interface array_p
3115
module procedure get_elem2p_
3116
end interface
3117
public :: array_p
3118
interface dimensions
3119
module procedure dimensions_
3120
end interface
3121
public :: dimensions
3122
interface distributed_index
3123
module procedure distributed_index_
3124
end interface
3125
public :: distributed_index
3126
interface distribute
3127
module procedure distribute_
3128
end interface
3129
public :: distribute
3130
interface new
3131
module procedure new_dist_index_
3132
end interface new
3133
type iD1D_z1D
3134
type(iD1D_z1D_), pointer :: D => null()
3135
integer :: error_ = 0
3136
end type iD1D_z1D
3137
type iD1D_z1D_
3138
type(iDist1D) :: e1
3139
type(zArray1D) :: e2
3140
integer :: refs_ = 0
3141
character(len=BUD_ID_LEN) :: id_ = "null_id"
3142
end type iD1D_z1D_
3143
interface new
3144
module procedure new_
3145
module procedure new_data_
3146
end interface
3147
public :: new
3148
interface assignment(=)
3149
module procedure get_elem1_assign_
3150
module procedure get_elem2_assign_
3151
module procedure set_elem1_
3152
module procedure set_elem2_
3153
end interface
3154
interface element
3155
module procedure get_elem1_
3156
module procedure get_elem2_
3157
end interface
3158
interface set_element
3159
module procedure set_elem1_
3160
module procedure set_elem2_
3161
end interface
3162
interface element1
3163
module procedure get_elem1_
3164
end interface
3165
public :: element1
3166
interface set_element1
3167
module procedure set_elem1_
3168
end interface
3169
public :: set_element1
3170
interface element1_p
3171
module procedure get_elem1p_
3172
end interface
3173
public :: element1_p
3174
interface element2
3175
module procedure get_elem2_
3176
end interface
3177
public :: element2
3178
interface set_element2
3179
module procedure set_elem2_
3180
end interface
3181
public :: set_element2
3182
interface element2_p
3183
module procedure get_elem2p_
3184
end interface
3185
public :: element2_p
3186
public :: iD1D_z1D
3187
private :: iD1D_z1D_
3188
interface assignment(=)
3189
module procedure common_assign_
3190
end interface
3191
public :: assignment(=)
3192
private :: common_assign_
3193
interface initialize
3194
module procedure common_initialize_
3195
end interface
3196
public :: initialize
3197
private :: common_initialize_
3198
interface is_initialized
3199
module procedure common_is_initialized_
3200
end interface
3201
public :: is_initialized
3202
private :: common_is_initialized_
3203
interface initialized
3204
module procedure common_is_initialized_
3205
end interface
3206
public :: initialized
3207
interface is_initd
3208
module procedure common_is_initialized_
3209
end interface
3210
public :: is_initd
3211
interface is_same
3212
module procedure common_is_same_
3213
end interface
3214
public :: is_same
3215
private :: common_is_same_
3216
interface same
3217
module procedure common_is_same_
3218
end interface
3219
public :: same
3220
interface delete
3221
module procedure common_delete_
3222
end interface
3223
public :: delete
3224
private :: common_delete_
3225
interface nullify
3226
module procedure common_nullify_
3227
end interface
3228
public :: nullify
3229
private :: common_nullify_
3230
interface copy
3231
module procedure copy_
3232
end interface
3233
public :: copy
3234
private :: common_copy_
3235
interface print
3236
module procedure print_
3237
end interface
3238
public :: print
3239
interface references
3240
module procedure common_references_
3241
end interface
3242
public :: references
3243
private :: common_references_
3244
interface refs
3245
module procedure common_references_
3246
end interface
3247
public :: refs
3248
interface set_error
3249
module procedure common_set_error_is_
3250
module procedure common_set_error_ii_
3251
module procedure common_set_error_il_
3252
end interface
3253
public :: set_error
3254
private :: common_set_error_is_
3255
private :: common_set_error_ii_
3256
private :: common_set_error_il_
3257
interface error
3258
module procedure common_error_
3259
end interface
3260
public :: error
3261
private :: common_error_
3262
contains
3263
subroutine common_copy_(from, to)
3264
type(iD1D_z1D), intent(in) :: from
3265
type(iD1D_z1D), intent(inout) :: to
3266
call set_error(to, error(from))
3267
end subroutine common_copy_
3268
subroutine common_initialize_(this)
3269
type(iD1D_z1D), intent(inout) :: this
3270
integer :: error
3271
call delete(this)
3272
allocate(this%D, stat=error)
3273
call set_error(this, error)
3274
if ( error /= 0 ) return
3275
this%D%refs_ = 1
3276
call common_tag_object_(this)
3277
end subroutine common_initialize_
3278
pure function common_is_initialized_(this) result(init)
3279
type(iD1D_z1D), intent(in) :: this
3280
logical :: init
3281
init = associated(this%D)
3282
end function common_is_initialized_
3283
elemental function common_is_same_(lhs, rhs) result(same)
3284
type(iD1D_z1D), intent(in) :: lhs, rhs
3285
logical :: same
3286
same = is_initd(lhs) .and. is_initd(rhs)
3287
if ( .not. same ) return
3288
same = associated(lhs%D, target=rhs%D)
3289
end function common_is_same_
3290
subroutine common_delete_(this)
3291
type(iD1D_z1D), intent(inout) :: this
3292
integer :: error
3293
call set_error(this, 0)
3294
if (.not. is_initd(this) ) return
3295
this%D%refs_ = this%D%refs_ - 1
3296
if ( 0 == this%D%refs_ ) then
3297
call delete_(this)
3298
deallocate(this%D, stat=error)
3299
call set_error(this, error)
3300
end if
3301
nullify(this%D)
3302
end subroutine common_delete_
3303
elemental subroutine common_nullify_(this)
3304
type(iD1D_z1D), intent(inout) :: this
3305
if (.not. is_initd(this) ) return
3306
nullify(this%D)
3307
end subroutine common_nullify_
3308
subroutine common_assign_(lhs, rhs)
3309
type(iD1D_z1D), intent(inout) :: lhs
3310
type(iD1D_z1D), intent(in) :: rhs
3311
call delete(lhs)
3312
if ( .not. is_initd(rhs) ) return
3313
lhs%D => rhs%D
3314
lhs%D%refs_ = rhs%D%refs_ + 1
3315
end subroutine common_assign_
3316
elemental function common_references_(this) result(refs)
3317
type(iD1D_z1D), intent(in) :: this
3318
integer :: refs
3319
if ( is_initd(this) ) then
3320
refs = this%D%refs_
3321
else
3322
refs = 0
3323
end if
3324
end function common_references_
3325
elemental function common_error_(this) result(error)
3326
type(iD1D_z1D), intent(in) :: this
3327
integer :: error
3328
if ( is_initd(this) ) then
3329
error = this%error_
3330
else
3331
error = 0
3332
end if
3333
end function common_error_
3334
elemental subroutine common_set_error_is_(this, error)
3335
type(iD1D_z1D), intent(inout) :: this
3336
integer(is_), intent(in) :: error
3337
this%error_ = error
3338
end subroutine common_set_error_is_
3339
elemental subroutine common_set_error_ii_(this, error)
3340
type(iD1D_z1D), intent(inout) :: this
3341
integer(ii_), intent(in) :: error
3342
this%error_ = error
3343
end subroutine common_set_error_ii_
3344
elemental subroutine common_set_error_il_(this, error)
3345
type(iD1D_z1D), intent(inout) :: this
3346
integer(il_), intent(in) :: error
3347
this%error_ = error
3348
end subroutine common_set_error_il_
3349
elemental function common_id_(this) result(str)
3350
type(iD1D_z1D), intent(in) :: this
3351
character(len=BUD_ID_LEn) :: str
3352
str = this%D%id_
3353
end function common_id_
3354
subroutine common_tag_object_(this)
3355
type(iD1D_z1D), intent(inout) :: this
3356
end subroutine common_tag_object_
3357
subroutine delete_(this)
3358
type(iD1D_z1D), intent(inout) :: this
3359
call set_error(this, 0)
3360
call delete(this%D%e1)
3361
if ( 0 /= error(this%D%e1) ) &
3362
call set_error(this, error(this%D%e1))
3363
call delete(this%D%e2)
3364
if ( 0 /= error(this%D%e2) ) &
3365
call set_error(this, error(this%D%e2))
3366
end subroutine delete_
3367
subroutine copy_(from, to)
3368
type(iD1D_z1D), intent(in) :: from
3369
type(iD1D_z1D), intent(inout) :: to
3370
call delete(to)
3371
if ( .not. is_initd(from) ) return
3372
call initialize(to)
3373
call common_copy_(from, to)
3374
call copy(from%D%e1, to%D%e1)
3375
call copy(from%D%e2, to%D%e2)
3376
end subroutine copy_
3377
subroutine new_data_(this, a, b &
3378
)
3379
type(iD1D_z1D), intent(inout) :: this
3380
type(iDist1D), intent(inout) :: a
3381
type(zArray1D), intent(inout) :: b
3382
call new(this)
3383
this%D%e1 = a
3384
this%D%e2 = b
3385
end subroutine new_data_
3386
subroutine new_(this)
3387
type(iD1D_z1D), intent(inout) :: this
3388
call initialize(this)
3389
end subroutine new_
3390
subroutine get_elem1_(this, item)
3391
type(iD1D_z1D), intent(in) :: this
3392
type(iDist1D), intent(inout) :: item
3393
if ( .not. is_initd(this) ) then
3394
call delete(item)
3395
else
3396
item = this%D% e1
3397
end if
3398
end subroutine
3399
subroutine get_elem1_assign_(item, this)
3400
type(iDist1D), intent(inout) :: item
3401
type(iD1D_z1D), intent(in) :: this
3402
if ( .not. is_initd(this) ) then
3403
call delete(item)
3404
else
3405
item = this%D% e1
3406
end if
3407
end subroutine
3408
subroutine set_elem1_(this, item)
3409
type(iD1D_z1D), intent(inout) :: this
3410
type(iDist1D), intent(in) :: item
3411
if ( .not. is_initd(this) ) return
3412
this%D% e1 = item
3413
end subroutine
3414
function get_elem1p_(this) result(p)
3415
type(iD1D_z1D), intent(inout) :: this
3416
type(iDist1D), pointer :: p
3417
if ( .not. is_initd(this) ) then
3418
nullify(p)
3419
else
3420
p => this%D% e1
3421
end if
3422
end function
3423
subroutine get_elem2_(this, item)
3424
type(iD1D_z1D), intent(in) :: this
3425
type(zArray1D), intent(inout) :: item
3426
if ( .not. is_initd(this) ) then
3427
call delete(item)
3428
else
3429
item = this%D% e2
3430
end if
3431
end subroutine
3432
subroutine get_elem2_assign_(item, this)
3433
type(zArray1D), intent(inout) :: item
3434
type(iD1D_z1D), intent(in) :: this
3435
if ( .not. is_initd(this) ) then
3436
call delete(item)
3437
else
3438
item = this%D% e2
3439
end if
3440
end subroutine
3441
subroutine set_elem2_(this, item)
3442
type(iD1D_z1D), intent(inout) :: this
3443
type(zArray1D), intent(in) :: item
3444
if ( .not. is_initd(this) ) return
3445
this%D% e2 = item
3446
end subroutine
3447
function get_elem2p_(this) result(p)
3448
type(iD1D_z1D), intent(inout) :: this
3449
type(zArray1D), pointer :: p
3450
if ( .not. is_initd(this) ) then
3451
nullify(p)
3452
else
3453
p => this%D% e2
3454
end if
3455
end function
3456
subroutine print_(this, info, indent)
3457
type(iD1D_z1D), intent(in) :: this
3458
character(len=*), intent(in), optional :: info
3459
integer, intent(in), optional :: indent
3460
integer :: lindent
3461
character(len=32) :: fmt
3462
character(len=256) :: name
3463
name = "iD1D_z1D"
3464
if ( present(info) ) name = info
3465
lindent = 1
3466
if ( present(indent) ) lindent = indent
3467
write(fmt, '(a,i0,a)') '(t',lindent,',3a)'
3468
if ( .not. is_initd(this) ) then
3469
write(*,fmt) "<", trim(name), " not initialized>"
3470
return
3471
end if
3472
write(fmt, '(a,i0,a)') '(t',lindent,',3a)'
3473
lindent = lindent + 2 ! step indentation
3474
write(*,fmt) "<<", trim(name), " coll>"
3475
call print(this%D%e1, indent = lindent)
3476
call print(this%D%e2, indent = lindent)
3477
lindent = lindent - 2 ! go back to requested indentation
3478
write(fmt, '(a,i0,a)') '(t',lindent,',a,i0,a)'
3479
write(*,fmt) " <coll-refs: ", references(this), ">>"
3480
end subroutine print_
3481
subroutine new_dist_index_(this, dist, arr, dist_idx)
3482
type(iD1D_z1D), intent(inout) :: this
3483
type(iDist1D), intent(inout) :: dist
3484
type(zArray1D), intent(inout) :: arr
3485
integer, intent(in) :: dist_idx
3486
call new(this, dist, arr)
3487
end subroutine new_dist_index_
3488
function distributed_index_(this) result(idx)
3489
type(iD1D_z1D), intent(in) :: this
3490
integer :: idx
3491
if ( is_initd(this) ) then
3492
idx = 1
3493
else
3494
idx = -1
3495
end if
3496
end function distributed_index_
3497
pure function dimensions_(this) result(d)
3498
type(iD1D_z1D), intent(in) :: this
3499
integer :: d
3500
if ( is_initd(this) ) then
3501
d = 1
3502
else
3503
d = -1
3504
end if
3505
end function dimensions_
3506
subroutine distribute_(this, parent, out_dist, out)
3507
type(iD1D_z1D), intent(inout) :: this
3508
type(MP_Comm), intent(inout) :: parent
3509
type(iDist1D), intent(inout) :: out_dist
3510
type(iD1D_z1D), intent(inout) :: out
3511
type(MP_Comm) :: comm, out_comm
3512
type(iDist1D) :: fake_dist, dist
3513
type(zArray1D) :: arr
3514
integer(ii_) :: ir, nranks, my_root, rank, in_rank
3515
logical :: is_distr, my_distr
3516
integer(ii_) :: dist_idx, dims
3517
integer(ii_) :: nl, ng
3518
integer(ii_) :: nout
3519
integer(ii_), allocatable :: ranks(:)
3520
integer(ii_), allocatable :: ashape(:), itmp1(:)
3521
if ( .not. is_communicator(parent) ) return
3522
rank = comm_rank(parent)
3523
nranks = comm_size(parent)
3524
dist = this
3525
in_rank = comm_rank(dist)
3526
if ( in_rank < 0 ) in_rank = in_rank - 1
3527
comm = dist
3528
arr = this
3529
dims = dimensions(arr)
3530
call AllReduce_Max(dims, ir, parent)
3531
dims = ir
3532
allocate(ashape(dims))
3533
if ( is_initd(arr) ) then
3534
dist_idx = distributed_index(this)
3535
do ir = 1 , dims
3536
ashape(ir) = size(arr, ir)
3537
end do
3538
else
3539
dist_idx = 0
3540
do ir = 1 , dims
3541
ashape(ir) = 0
3542
end do
3543
end if
3544
call AllReduce_Max(dist_idx, ir, parent)
3545
dist_idx = ir
3546
allocate(itmp1(dims))
3547
call AllReduce_Max(ashape, itmp1, parent)
3548
do ir = 1 , dims
3549
if ( ir /= dist_idx ) then
3550
ashape(ir) = itmp1(ir)
3551
end if
3552
end do
3553
deallocate(itmp1)
3554
out_comm = out_dist
3555
if ( is_communicator(out_comm) ) then
3556
if ( comm_rank(out_comm) == 0 ) then
3557
ir = comm_rank(parent)
3558
else
3559
ir = -1
3560
end if
3561
call AllReduce_Max(ir, my_root, out_comm)
3562
else
3563
my_root = -1
3564
end if
3565
ng = size_global(dist)
3566
call AllReduce_Max(ng, nl, parent)
3567
ng = nl
3568
nl = size_local(dist)
3569
do ir = 0 , nranks - 1
3570
if ( my_root == ir ) then
3571
my_distr = .true.
3572
else
3573
my_distr = .false.
3574
end if
3575
call AllReduce_LOr(my_distr, is_distr, parent)
3576
if ( .not. is_distr ) cycle
3577
if ( my_distr ) then
3578
call new_remote(parent, out_dist)
3579
call create_ranks(out_dist)
3580
call sub_dist(out_dist)
3581
else
3582
call new_remote(parent, fake_dist)
3583
call create_ranks(fake_dist)
3584
call sub_dist(fake_dist)
3585
call delete(fake_dist)
3586
end if
3587
deallocate(ranks)
3588
end do
3589
call delete(dist)
3590
call delete(comm)
3591
call delete(arr)
3592
call delete(out_comm)
3593
contains
3594
subroutine create_ranks(dist)
3595
type(iDist1D), intent(inout) :: dist
3596
type(MP_Comm) :: comm
3597
comm = dist
3598
call child_Bcast(parent, comm, size=nout)
3599
allocate(ranks(-1:nout-1))
3600
ranks(-1) = -1
3601
call child_Bcast_ranks(parent, comm, nout, ranks(0:))
3602
call delete(comm)
3603
end subroutine create_ranks
3604
subroutine sub_dist(out_dist)
3605
use bud_Transfer, only: transfer_dim
3606
type(iDist1D), intent(inout) :: out_dist
3607
type(MP_Comm) :: out_comm
3608
type(zArray1D) :: out_arr
3609
logical :: run
3610
complex(rd_), pointer :: dat (:)
3611
complex(rd_), pointer :: odat (:)
3612
integer(ii_) :: il, ig, oil
3613
integer(ii_) :: out_nl
3614
integer(ii_) :: i2
3615
integer :: send_R, recv_R
3616
integer(ii_), allocatable :: reqs(:), stats(:,:)
3617
integer(ii_) :: stat(MPI_STATUS_SIZE)
3618
run = .true.
3619
if ( .not. run ) return
3620
out_comm = out_dist
3621
out_nl = size_local(out_dist)
3622
if ( is_communicator(out_comm) ) then
3623
ashape(dist_idx) = out_nl
3624
call new(out_arr, ashape(:))
3625
call new(out, out_dist, out_arr)
3626
odat => array_p(out_arr)
3627
call delete(out_arr)
3628
end if
3629
if ( is_initd(dist) ) then
3630
allocate(reqs(nl), stats(MPI_STATUS_SIZE, nl))
3631
do il = 1 , nl
3632
reqs(il) = MPI_REQUEST_NULL
3633
end do
3634
end if
3635
dat => array_p(arr)
3636
if ( dist_idx == dims ) then
3637
do ig = 1 , ng
3638
send_R = g2rank(dist, ig)
3639
recv_R = ranks(g2rank(out_dist, ig))
3640
if ( recv_R == rank .and. &
3641
send_R == in_rank ) then
3642
il = g2l(dist, ig)
3643
oil = g2l(out_dist, ig)
3644
odat(oil) = dat(il)
3645
else if ( recv_R == rank ) then
3646
oil = g2l(out_dist, ig)
3647
call Recv(odat(oil), MPI_ANY_SOURCE, ig, &
3648
parent, stat)
3649
else if ( send_R == in_rank ) then
3650
il = g2l(dist, ig)
3651
call ISSend(dat(il), recv_R, ig, &
3652
parent, reqs(il))
3653
end if
3654
end do
3655
if ( allocated(reqs) ) &
3656
call WaitAll(nl, reqs, stats, parent)
3657
else if ( dist_idx == 1 ) then
3658
do i2 = 1 , ashape(2)
3659
do ig = 1 , ng
3660
send_R = g2rank(dist, ig)
3661
recv_R = ranks(g2rank(out_dist, ig))
3662
if ( recv_R == rank .and. &
3663
send_R == in_rank ) then
3664
il = g2l(dist, ig)
3665
oil = g2l(out_dist, ig)
3666
else if ( recv_R == rank ) then
3667
oil = g2l(out_dist, ig)
3668
else if ( send_R == in_rank ) then
3669
il = g2l(dist, ig)
3670
end if
3671
end do ! ig
3672
if ( allocated(reqs) ) &
3673
call WaitAll(nl, reqs, stats, parent)
3674
end do ! i2
3675
end if
3676
if ( allocated(reqs) ) deallocate(reqs)
3677
call Barrier(parent)
3678
call delete(out_comm)
3679
end subroutine sub_dist
3680
end subroutine distribute_
3681
subroutine write_(f, this)
3682
use bud_File
3683
use bud_Transfer, only: transfer_dim
3684
type( File ), intent(inout) :: f
3685
type(iD1D_z1D), intent(in) :: this
3686
type(iDist1D) :: dist
3687
type(zArray1D) :: arr
3688
type( MP_Comm ) :: comm
3689
logical :: formatted, do_io
3690
integer :: iu, io_rank, dat_rank, rank
3691
character(len=64), parameter :: fmt_ = '(e20.16)'
3692
complex(rd_), pointer :: dat (:)
3693
complex(rd_), allocatable :: rdat(:)
3694
integer :: dist_idx, l, n, ndat
3695
integer :: d1, i1
3696
integer :: status(MPI_STATUS_SIZE)
3697
if ( .not. is_initd(this) ) return
3698
dist = this
3699
comm = dist
3700
arr = this
3701
rank = comm_rank(comm)
3702
if ( is_open(f) ) then
3703
io_rank = rank
3704
else
3705
io_rank = -1
3706
end if
3707
call AllReduce_Max(io_rank, iu, comm)
3708
io_rank = iu
3709
if ( io_rank < 0 ) then
3710
call delete(dist)
3711
call delete(comm)
3712
call delete(arr)
3713
return
3714
end if
3715
do_io = io_rank == rank
3716
if ( do_io ) then
3717
formatted = is_formatted(f)
3718
iu = unit(f)
3719
end if
3720
d1 = size(arr, 1)
3721
dat => array_p(arr)
3722
call delete(arr)
3723
dist_idx = distributed_index(this)
3724
if ( dist_idx == 1 ) then
3725
ndat = size_global(dist)
3726
else
3727
ndat = d1
3728
end if
3729
if ( do_io ) then
3730
allocate(rdat(ndat))
3731
if ( formatted ) then
3732
write(iu, '(i16)') 1
3733
write(iu, '(i16)') d1
3734
else
3735
write(iu) 1
3736
write(iu) d1
3737
end if
3738
end if
3739
if ( dist_idx == 1 ) then
3740
i1 = 1
3741
do while ( i1 <= d1 )
3742
n = consecutive(dist, i1)
3743
dat_rank = global2rank(dist, i1)
3744
if ( do_io .and. rank == dat_rank ) then
3745
l = global2local(dist, i1)
3746
call transfer_dim(n, rdat(i1:), n, dat(l:))
3747
else if ( do_io ) then
3748
call Recv(rdat(i1), n, dat_rank, i1, comm, status)
3749
else if ( rank == dat_rank ) then
3750
l = global2local(dist, i1)
3751
call SSend(dat(l), n, io_rank, i1, comm)
3752
end if
3753
i1 = i1 + n
3754
end do
3755
if ( do_io ) then
3756
if ( formatted ) then
3757
write(iu, fmt_) rdat(:)
3758
else
3759
write(iu) rdat(:)
3760
end if
3761
end if
3762
end if
3763
if ( do_io ) deallocate(rdat)
3764
call delete(dist)
3765
call delete(comm)
3766
end subroutine write_
3767
subroutine read_(f, dist, this, dist_idx)
3768
use bud_File
3769
use bud_Transfer, only: transfer_dim
3770
type( File ), intent(inout) :: f
3771
type(iDist1D), intent(inout) :: dist
3772
type(iD1D_z1D), intent(inout) :: this
3773
integer, intent(in), optional :: dist_idx
3774
type(zArray1D) :: arr
3775
type( MP_Comm ) :: comm
3776
logical :: formatted, do_io
3777
integer :: iu, io_rank, dat_rank, rank
3778
character(len=64), parameter :: fmt_ = '(e20.16)'
3779
integer(ii_) :: ashape(1)
3780
complex(rd_), pointer :: dat (:)
3781
complex(rd_), allocatable :: rdat(:)
3782
integer :: ldist_idx, l, n
3783
integer :: d1, i1
3784
integer :: status(MPI_STATUS_SIZE)
3785
if ( .not. is_initd(this) ) return
3786
ldist_idx = 1
3787
if ( present(dist_idx) ) ldist_idx = dist_idx
3788
comm = dist
3789
rank = comm_rank(comm)
3790
if ( is_open(f) ) then
3791
io_rank = rank
3792
else
3793
io_rank = -1
3794
end if
3795
call AllReduce_Max(io_rank, iu, comm)
3796
io_rank = iu
3797
if ( io_rank < 0 ) then
3798
call delete(comm)
3799
call delete(this)
3800
call set_error(this, -1)
3801
return
3802
end if
3803
do_io = io_rank == rank
3804
if ( do_io ) then
3805
formatted = is_formatted(f)
3806
iu = unit(f)
3807
if ( formatted ) then
3808
read(iu, '(i16)') l
3809
read(iu, '(i16)') ashape
3810
d1 = ashape(1)
3811
else
3812
read(iu) l
3813
read(iu) ashape
3814
d1 = ashape(1)
3815
end if
3816
end if
3817
select case ( ldist_idx )
3818
case ( 1 )
3819
ashape(1) = size_local(dist)
3820
end select
3821
call new(arr, ashape)
3822
dat => array_p(arr)
3823
call new(this, dist, arr, ldist_idx)
3824
call delete(arr)
3825
if ( do_io ) then
3826
allocate(rdat(d1))
3827
end if
3828
if ( do_io ) then
3829
if ( formatted ) then
3830
read(iu, fmt_) rdat(:)
3831
else
3832
read(iu) rdat(:)
3833
end if
3834
end if
3835
select case ( ldist_idx )
3836
case ( 1 )
3837
i1 = 1
3838
do while ( i1 <= d1 )
3839
n = consecutive(dist, i1)
3840
dat_rank = global2rank(dist, i1)
3841
if ( do_io .and. rank == dat_rank ) then
3842
l = global2local(dist, i1)
3843
call transfer_dim(n, dat(l:), n, rdat(i1:))
3844
else if ( do_io ) then
3845
call SSend(rdat(i1), n, dat_rank, i1, comm)
3846
else if ( rank == dat_rank ) then
3847
l = global2local(dist, i1)
3848
call Recv(dat(l), n, io_rank, i1, comm, status)
3849
end if
3850
i1 = i1 + n
3851
end do
3852
end select
3853
call delete(comm)
3854
end subroutine read_
3855
end module
Loggerhead is a web-based interface for Breezy
Version: 2.0.1