-
Notifications
You must be signed in to change notification settings - Fork 1
/
smeagol_matrix_utils.F
1933 lines (1701 loc) · 103 KB
/
smeagol_matrix_utils.F
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
! !
! SPDX-License-Identifier: GPL-2.0-or-later !
!--------------------------------------------------------------------------------------------------!
! **************************************************************************************************
!> \brief Routines to convert sparse matrices between DBCSR (distributed-blocks compressed sparse rows)
!> and SIESTA (distributed compressed sparse columns) formats.
!> \author Sergey Chulkov
!> \author Christian Ahart
!> \author Clotilde Cucinotta
! **************************************************************************************************
MODULE smeagol_matrix_utils
USE cell_types, ONLY: cell_type, &
real_to_scaled, &
scaled_to_real
USE cp_dbcsr_api, ONLY: dbcsr_get_block_p, &
dbcsr_get_info, &
dbcsr_p_type, &
dbcsr_set
USE kinds, ONLY: dp, &
dp_size, &
int_8
USE message_passing, ONLY: mp_para_env_type, &
mp_request_type, &
mp_waitall
USE negf_matrix_utils, ONLY: get_index_by_cell
#if defined(__SMEAGOL)
USE parallel, ONLY: GetNodeOrbs, &
GlobalToLocalOrb, &
LocalToGlobalOrb, &
WhichNodeOrb
#endif
USE particle_types, ONLY: particle_type
USE qs_neighbor_list_types, ONLY: get_iterator_info, &
get_neighbor_list_set_p, &
neighbor_list_iterate, &
neighbor_list_iterator_create, &
neighbor_list_iterator_p_type, &
neighbor_list_iterator_release, &
neighbor_list_set_p_type
USE qs_subsys_types, ONLY: qs_subsys_get, &
qs_subsys_type
USE util, ONLY: sort
#include "./base/base_uses.f90"
IMPLICIT NONE
PRIVATE
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'smeagol_matrix_utils'
LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE.
INTEGER, PARAMETER, PRIVATE :: neighbor_list_iatom_index = 1
INTEGER, PARAMETER, PRIVATE :: neighbor_list_jatom_index = 2
INTEGER, PARAMETER, PRIVATE :: neighbor_list_dbcsr_image_index = 3
INTEGER, PARAMETER, PRIVATE :: neighbor_list_siesta_image_index = 4
INTEGER, PARAMETER, PRIVATE :: neighbor_list_siesta_transp_image_index = 5
INTEGER, PARAMETER, PRIVATE :: neighbor_list_dim1 = neighbor_list_siesta_transp_image_index
PUBLIC :: siesta_distrib_csc_struct_type
PUBLIC :: siesta_struct_create, siesta_struct_release
PUBLIC :: convert_dbcsr_to_distributed_siesta, convert_distributed_siesta_to_dbcsr
PRIVATE :: get_negf_cell_ijk, index_in_canonical_enumeration, number_from_canonical_enumeration, pbc_0_1
PRIVATE :: get_number_of_mpi_sendrecv_requests, assign_nonzero_elements_to_requests
!> number of DBCSR matrix elements to receive from a given rank
INTEGER, PARAMETER, PRIVATE :: nelements_dbcsr_recv = 1
!> number of DBCSR matrix elements to send to a given rank
INTEGER, PARAMETER, PRIVATE :: nelements_dbcsr_send = 2
INTEGER, PARAMETER, PRIVATE :: nelements_dbcsr_dim2 = nelements_dbcsr_send
INTEGER(kind=int_8), PARAMETER, PRIVATE :: max_mpi_packet_size_bytes = 134217728 ! 128 MiB (to limit memory usage for matrix redistribution)
INTEGER(kind=int_8), PARAMETER, PRIVATE :: max_mpi_packet_size_dp = max_mpi_packet_size_bytes/INT(dp_size, kind=int_8)
! a portable way to determine the upper bound for tag value is to call
! MPI_COMM_GET_ATTR(comm, MPI_TAG_UB, max_mpi_rank, flag, ierror).
! The MPI specification guarantees a value of 32767
INTEGER, PARAMETER, PRIVATE :: max_mpi_rank = 32767
! **************************************************************************************************
!> \brief Sparsity pattern of replicated SIESTA compressed sparse column (CSC) matrices
! **************************************************************************************************
TYPE siesta_distrib_csc_struct_type
!> gather all non-zero matrix elements on the given MPI process.
!> Distribute the elements across MPI processes if gather_root < 0.
!> The matrix elements should be located on I/O process in case of bulk transport,
!> and should be distributed in case of SMEAGOL calculation.
INTEGER :: gather_root = 0
!> Based of full (.FALSE.) or upper-triangular (.TRUE.) DBCSR matrix.
!> It is used in CPASSERTs to allow access to lower-triangular matrix elements of non-symmetric DBCSR matrices
!> In case there is no bugs in this module, these CPASSERTs should newer trigger.
!> Therefore the 'symmetric' variable alongside with relevant CPASSERT calls are excessive and can be removed.
LOGICAL :: symmetric = .TRUE.
!> number of neighbour list nodes for each MPI process (0:num_pe-1).
!> If do_merge == .TRUE., nodes for different cell images along k cell vector are merged into one node
INTEGER, ALLOCATABLE, DIMENSION(:) :: nnodes_per_proc
!> replicated neighbour list (1:neighbor_list_dim1, 1:SUM(nnodes_per_proc)).
!> Neighbour list nodes are ordered according to their MPI ranks.
!> Thus, the first nnodes_per_proc(0) nodes are stored on MPI rank 0,
!> the next nnodes_per_proc(1) nodes reside on MPI rank 1, etc
!> Nodes for cell images along transport direction are merged into one node.
!> The number of non-zero DBCSR matrix blocks and their DBCSR cell image indices
!> are stored into 'n_dbcsr_cell_images_to_merge' and 'dbcsr_cell_image_to_merge' arrays
INTEGER, ALLOCATABLE, DIMENSION(:, :) :: nl_repl
!> number of DBCSR images for each local merged neighbour list node (1:nnodes_per_proc(para_env%mepos))
INTEGER, ALLOCATABLE, DIMENSION(:) :: n_dbcsr_cell_images_to_merge
!> list of DBCSR image indices to merge; (1:SUM(n_dbcsr_cell_images_to_merge))
INTEGER, ALLOCATABLE, DIMENSION(:) :: dbcsr_cell_image_to_merge
!> number of DBCSR non-zero matrix elements that should be received/sent from each MPI rank
!> (0:num_pe-1, 1:nelements_dbcsr_dim2)
INTEGER(kind=int_8), ALLOCATABLE, DIMENSION(:, :) :: nelements_per_proc
!> number of non-zero matrix elements local to this MPI rank.
INTEGER(kind=int_8) :: n_nonzero_elements = 0_int_8
INTEGER :: nrows = 0, ncols = 0
!> Number of non-zero matrix elements (columns) on each row
!> n_nonzero_cols(1:nrows); same as 'numh' in SMEAGOL code.
INTEGER, ALLOCATABLE, DIMENSION(:) :: n_nonzero_cols
!> offset of the first non-zero matrix elements on each row.
!> column_offset(1:nrows); same as 'listhptr' in SMEAGOL code.
!> It should be declared as INTEGER(kind=int_8), but SMEAGOL expects it to be INTEGER
INTEGER, ALLOCATABLE, DIMENSION(:) :: row_offset
!> column index of each non-zero matrix element.
!> col_index(1:n_nonzero_elements); same as 'listh' in SMEAGOL code
!> col index of the first non-zero matrix elements of irow row is col_index(row_offset(irow)+1)
INTEGER, ALLOCATABLE, DIMENSION(:) :: col_index
!> index of the non-zero matrix element in a communication buffer for each rank
INTEGER, ALLOCATABLE, DIMENSION(:) :: packed_index
!> R_atom_row - R_atom_col
REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: xij
!> equivalent atomic orbitals
INTEGER, ALLOCATABLE, DIMENSION(:) :: indxuo
!> atomic index on which the orbital is centred
INTEGER, ALLOCATABLE, DIMENSION(:) :: iaorb
!> coordinates of all atoms in the supercell
REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: xa
END TYPE siesta_distrib_csc_struct_type
CONTAINS
! **************************************************************************************************
!> \brief Map non-zero matrix blocks between sparse matrices in DBCSR and SIESTA formats.
!> \param siesta_struct structure that stores metadata (sparsity pattern) of sparse SIESTA matrices
!> \param matrix_dbcsr_kp DBCSR matrices for each cell image
!> \param subsys QuickStep molecular system
!> \param cell_to_index array to convert 3-D cell indices to 1-D DBCSR image indices
!> \param sab_nl pair-wise neighbour list
!> \param para_env MPI parallel environment
!> \param max_ij_cell_image largest index of cell images along i and j cell vectors (e.g. (2,0) in
!> case of 5 cell images (0,0), (1,0), (-1,0), (2,0), and (-2,0))
!> \param do_merge merge DBCSR images along transport direction (k cell vector)
!> \param gather_root distribute non-zero matrix elements of SIESTA matrices across all
!> parallel processes (-1), or gather them on the given MPI rank (>= 0).
! **************************************************************************************************
SUBROUTINE siesta_struct_create(siesta_struct, matrix_dbcsr_kp, subsys, cell_to_index, &
sab_nl, para_env, max_ij_cell_image, do_merge, gather_root)
TYPE(siesta_distrib_csc_struct_type), &
INTENT(inout) :: siesta_struct
TYPE(dbcsr_p_type), DIMENSION(:), INTENT(in) :: matrix_dbcsr_kp
TYPE(qs_subsys_type), POINTER :: subsys
INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
POINTER :: sab_nl
TYPE(mp_para_env_type), POINTER :: para_env
INTEGER, DIMENSION(2), INTENT(inout) :: max_ij_cell_image
LOGICAL, INTENT(in) :: do_merge
INTEGER, INTENT(in) :: gather_root
CHARACTER(len=*), PARAMETER :: routineN = 'siesta_struct_create'
CHARACTER(len=20) :: str_nelem, str_nelem_max
INTEGER :: handle, iatom, icol, icol_blk, icol_local, image, image_j, image_k, irow, &
irow_local, natoms, ncells_siesta_total, ncols_blk, ncols_total, nrows_local, &
nrows_total, offset
INTEGER(kind=int_8) :: n_nonzero_elements_local
INTEGER, DIMENSION(3) :: max_ijk_cell_image, ncells_siesta
INTEGER, DIMENSION(:), POINTER :: col_blk_offset, col_blk_size
LOGICAL :: do_distribute, is_root_rank
REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: particle_coords
REAL(kind=dp), DIMENSION(3) :: real_cell_shift, scaled_cell_shift
TYPE(cell_type), POINTER :: cell
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
CALL timeset(routineN, handle)
do_distribute = gather_root < 0
is_root_rank = gather_root == para_env%mepos
! here row_blk_offset / col_blk_offset are global indices of the first row / column of a given non-zero block.
! They are not offsets (index-1) but the actual indices.
CALL dbcsr_get_info(matrix=matrix_dbcsr_kp(1)%matrix, &
nfullrows_total=nrows_total, nfullcols_total=ncols_total, &
nblkcols_total=ncols_blk, col_blk_size=col_blk_size, col_blk_offset=col_blk_offset)
IF (debug_this_module) THEN
CPASSERT(nrows_total == ncols_total)
CPASSERT(gather_root < para_env%num_pe)
END IF
siesta_struct%gather_root = gather_root
CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=siesta_struct%symmetric)
! apply periodic boundary conditions to atomic coordinates
CALL qs_subsys_get(subsys, cell=cell, particle_set=particle_set, nparticle=natoms)
ALLOCATE (particle_coords(3, natoms))
DO iatom = 1, natoms
CALL pbc_0_1(particle_coords(1:3, iatom), particle_set(iatom)%r(1:3), cell)
END DO
! Note: in case we would like to limit the number of cell images along transport direction (k cell vector)
! by enabling 'BulkTransvCellSizeZ' keyword, we need to pass max_ijk_cell_image(1:3) vector
! to the subroutine instead of the reduced vector max_ij_cell_image(1:2)
max_ijk_cell_image(1:2) = max_ij_cell_image(1:2)
! determine the actual number of cell images along k cell vector. Unless the third element is also passed
! via subroutine arguments, an extra MPI_Allreduce operation is needed each time we call
! replicate_neighbour_list() / get_nnodes_local().
max_ijk_cell_image(3) = -1
IF (.NOT. do_merge) max_ijk_cell_image(3) = 1 ! bulk-transport calculation expects exactly 3 cell images along transport direction
! replicate pair-wise neighbour list. Identical non-zero matrix blocks from cell image along transport direction
! are grouped together if do_merge == .TRUE.
ALLOCATE (siesta_struct%nnodes_per_proc(0:para_env%num_pe - 1))
CALL replicate_neighbour_list(siesta_struct%nl_repl, &
siesta_struct%n_dbcsr_cell_images_to_merge, &
siesta_struct%dbcsr_cell_image_to_merge, &
siesta_struct%nnodes_per_proc, &
max_ijk_cell_image, sab_nl, para_env, particle_coords, cell, cell_to_index, do_merge)
max_ij_cell_image(1:2) = max_ijk_cell_image(1:2)
! count number of non-zero matrix elements that need to be send to and received from other parallel processes
ALLOCATE (siesta_struct%nelements_per_proc(0:para_env%num_pe - 1, nelements_dbcsr_dim2))
CALL count_remote_dbcsr_elements(siesta_struct%nelements_per_proc, siesta_struct%nnodes_per_proc, &
siesta_struct%nl_repl, matrix_dbcsr_kp, siesta_struct%symmetric, para_env, gather_root)
! number of SIESTA non-zero matrix elements that are going to be stored on this parallel process
n_nonzero_elements_local = SUM(siesta_struct%nelements_per_proc(:, nelements_dbcsr_recv))
siesta_struct%n_nonzero_elements = n_nonzero_elements_local
! as SMEAGOL uses 32-bits integers, the number of non-zero matrix elements is limited by 2^31 per MPI rank.
! Abort CP2K if we are about to exceed this limit.
IF (n_nonzero_elements_local > INT(HUGE(0), kind=int_8)) THEN
WRITE (str_nelem, '(I0)') n_nonzero_elements_local
WRITE (str_nelem_max, '(I0)') HUGE(0)
CALL cp_abort(__LOCATION__, &
"The number of non-zero matrix elements per MPI process "//TRIM(str_nelem)// &
" cannot exceed "//TRIM(str_nelem_max)// &
". Please increase the number of MPI processes to satisfy this SMEAGOL limitation.")
END IF
! in case there is no nonzero matrix element stored on this process, allocate arrays with one element to avoid SEGFAULT
IF (n_nonzero_elements_local == 0) n_nonzero_elements_local = 1
! number of SIESTA-matrix rows local to the given parallel process
IF (do_distribute) THEN
#if defined(__SMEAGOL)
CALL GetNodeOrbs(nrows_total, para_env%mepos, para_env%num_pe, nrows_local)
#else
CALL cp_abort(__LOCATION__, &
"CP2K was compiled with no SMEAGOL support.")
#endif
ELSE
IF (is_root_rank) THEN
nrows_local = nrows_total
ELSE
nrows_local = 0
END IF
END IF
! number of cell images along each cell vector. It is 2*m+1, as SIESTA images are ordered as 0, 1, -1, ..., m, -m
ncells_siesta(1:3) = 2*max_ijk_cell_image(1:3) + 1
! in case of merged cell images along the transport direction, there will be just 1 'merged' cell image along it
IF (do_merge) ncells_siesta(3) = 1
ncells_siesta_total = ncells_siesta(1)*ncells_siesta(2)*ncells_siesta(3)
! number of rows local to the given parallel process. Rows are distributed in a block-cyclic manner
siesta_struct%nrows = nrows_local
! number of columns of the matrix in its dense form. SIESTA uses 1-D (rows) block-cyclic distribution.
! All non-zero matrix elements on a given row are stored on the same parallel process
siesta_struct%ncols = nrows_total*ncells_siesta_total
! allocate at least one array element to avoid SIGFAULT when passing unallocated arrays to subroutines
IF (nrows_local == 0) nrows_local = 1
ALLOCATE (siesta_struct%n_nonzero_cols(nrows_local))
ALLOCATE (siesta_struct%row_offset(nrows_local))
ALLOCATE (siesta_struct%col_index(n_nonzero_elements_local))
ALLOCATE (siesta_struct%packed_index(n_nonzero_elements_local))
! restore the actual number of local rows
nrows_local = siesta_struct%nrows
! get number of non-zero matrix element on each local row (n_nonzero_cols),
! offset of the first non-zero matrix element for each local row (row_offset),
! global column indices of all local non-zero matrix elements (col_index), and
! the indices of all local non-zero matrix elements in the communication buffer (packed_index)
CALL get_nonzero_element_indices(siesta_struct%n_nonzero_cols, siesta_struct%row_offset, &
siesta_struct%col_index, siesta_struct%packed_index, &
siesta_struct%nl_repl, matrix_dbcsr_kp, &
siesta_struct%symmetric, para_env, gather_root)
! indices of equivalent atomic orbitals
ALLOCATE (siesta_struct%indxuo(siesta_struct%ncols))
DO icol = 1, ncols_total
siesta_struct%indxuo(icol) = icol
END DO
DO image = 2, ncells_siesta_total
siesta_struct%indxuo((image - 1)*ncols_total + 1:image*ncols_total) = siesta_struct%indxuo(1:ncols_total)
END DO
! particle index on which the orbital is centred
ALLOCATE (siesta_struct%iaorb(siesta_struct%ncols))
DO icol_blk = 1, ncols_blk
! col_blk_offset() is not an offset but the index of the first atomic orbital in the column block
siesta_struct%iaorb(col_blk_offset(icol_blk):col_blk_offset(icol_blk) + col_blk_size(icol_blk) - 1) = icol_blk
END DO
DO image = 2, ncells_siesta_total
siesta_struct%iaorb((image - 1)*ncols_total + 1:image*ncols_total) = siesta_struct%iaorb(1:ncols_total) + (image - 1)*natoms
END DO
! coordinates of all particles in each cell images
ALLOCATE (siesta_struct%xa(3, natoms*ncells_siesta_total))
DO image_k = 1, ncells_siesta(3)
!icell_siesta(3) = image_k
scaled_cell_shift(3) = REAL(number_from_canonical_enumeration(image_k), kind=dp) ! SIESTA -> actual cell index
DO image_j = 1, ncells_siesta(2)
scaled_cell_shift(2) = REAL(number_from_canonical_enumeration(image_j), kind=dp)
DO image = 1, ncells_siesta(1)
scaled_cell_shift(1) = REAL(number_from_canonical_enumeration(image), kind=dp)
CALL scaled_to_real(real_cell_shift, scaled_cell_shift, cell)
offset = (((image_k - 1)*ncells_siesta(2) + image_j - 1)*ncells_siesta(1) + image - 1)*natoms
DO iatom = 1, natoms
siesta_struct%xa(1:3, offset + iatom) = particle_set(iatom)%r(1:3) + real_cell_shift(1:3)
END DO
END DO
END DO
END DO
! inter-atomic distance
ALLOCATE (siesta_struct%xij(3, n_nonzero_elements_local))
DO irow_local = 1, nrows_local
IF (do_distribute) THEN
#if defined(__SMEAGOL)
CALL LocalToGlobalOrb(irow_local, para_env%mepos, para_env%num_pe, irow)
#else
CALL cp_abort(__LOCATION__, &
"CP2K was compiled with no SMEAGOL support.")
#endif
ELSE
irow = irow_local
IF (debug_this_module) THEN
CPASSERT(is_root_rank)
END IF
END IF
offset = siesta_struct%row_offset(irow_local)
DO icol_local = offset + 1, offset + siesta_struct%n_nonzero_cols(irow_local)
icol = siesta_struct%col_index(icol_local)
siesta_struct%xij(1:3, icol_local) = siesta_struct%xa(1:3, siesta_struct%iaorb(icol)) - &
siesta_struct%xa(1:3, siesta_struct%iaorb(irow))
END DO
END DO
DEALLOCATE (particle_coords)
CALL timestop(handle)
END SUBROUTINE siesta_struct_create
! **************************************************************************************************
!> \brief Release a SIESTA matrix structure
!> \param siesta_struct structure to release
! **************************************************************************************************
SUBROUTINE siesta_struct_release(siesta_struct)
TYPE(siesta_distrib_csc_struct_type), &
INTENT(inout) :: siesta_struct
CHARACTER(len=*), PARAMETER :: routineN = 'siesta_struct_release'
INTEGER :: handle
CALL timeset(routineN, handle)
siesta_struct%gather_root = -1
IF (ALLOCATED(siesta_struct%nnodes_per_proc)) DEALLOCATE (siesta_struct%nnodes_per_proc)
IF (ALLOCATED(siesta_struct%nl_repl)) DEALLOCATE (siesta_struct%nl_repl)
IF (ALLOCATED(siesta_struct%n_dbcsr_cell_images_to_merge)) DEALLOCATE (siesta_struct%n_dbcsr_cell_images_to_merge)
IF (ALLOCATED(siesta_struct%dbcsr_cell_image_to_merge)) DEALLOCATE (siesta_struct%dbcsr_cell_image_to_merge)
IF (ALLOCATED(siesta_struct%nelements_per_proc)) DEALLOCATE (siesta_struct%nelements_per_proc)
siesta_struct%n_nonzero_elements = 0
siesta_struct%nrows = 0
siesta_struct%ncols = 0
IF (ALLOCATED(siesta_struct%n_nonzero_cols)) DEALLOCATE (siesta_struct%n_nonzero_cols)
IF (ALLOCATED(siesta_struct%row_offset)) DEALLOCATE (siesta_struct%row_offset)
IF (ALLOCATED(siesta_struct%col_index)) DEALLOCATE (siesta_struct%col_index)
IF (ALLOCATED(siesta_struct%packed_index)) DEALLOCATE (siesta_struct%packed_index)
IF (ALLOCATED(siesta_struct%xij)) DEALLOCATE (siesta_struct%xij)
IF (ALLOCATED(siesta_struct%indxuo)) DEALLOCATE (siesta_struct%indxuo)
IF (ALLOCATED(siesta_struct%iaorb)) DEALLOCATE (siesta_struct%iaorb)
IF (ALLOCATED(siesta_struct%xa)) DEALLOCATE (siesta_struct%xa)
CALL timestop(handle)
END SUBROUTINE siesta_struct_release
! **************************************************************************************************
!> \brief Convert matrix from DBCSR to sparse SIESTA format.
!> \param matrix_siesta matrix in SIESTA format [out]
!> \param matrix_dbcsr_kp DBCSR matrix [in]
!> \param siesta_struct structure to map matrix blocks between formats
!> \param para_env MPI parallel environment
! **************************************************************************************************
SUBROUTINE convert_dbcsr_to_distributed_siesta(matrix_siesta, matrix_dbcsr_kp, siesta_struct, para_env)
REAL(kind=dp), DIMENSION(:), INTENT(out) :: matrix_siesta
TYPE(dbcsr_p_type), DIMENSION(:), INTENT(in) :: matrix_dbcsr_kp
TYPE(siesta_distrib_csc_struct_type), INTENT(in) :: siesta_struct
TYPE(mp_para_env_type), INTENT(in), POINTER :: para_env
CHARACTER(len=*), PARAMETER :: routineN = 'convert_dbcsr_to_distributed_siesta'
INTEGER :: first_col_minus_one, first_row_minus_one, handle, icol_blk, icol_local, &
image_dbcsr, image_ind, image_ind_offset, image_siesta, image_siesta_transp, inode, &
inode_proc, iproc, irequest, irow_blk, irow_local, irow_proc, mepos, n_image_ind, &
ncols_blk, ncols_local, nnodes_proc, node_offset, nprocs, nrequests_recv, &
nrequests_total, nrows_blk, nrows_local
INTEGER(kind=int_8) :: n_nonzero_elements_dbcsr, &
n_nonzero_elements_siesta, &
offset_recv_mepos, offset_send_mepos
INTEGER(kind=int_8), ALLOCATABLE, DIMENSION(:) :: n_packed_elements_per_proc, &
nelements_per_request, &
offset_per_proc, offset_per_request
INTEGER, ALLOCATABLE, DIMENSION(:) :: next_nonzero_element_offset, peer_rank, &
request_tag
INTEGER, DIMENSION(:), POINTER :: col_blk_offset, col_blk_size, &
row_blk_offset, row_blk_size
LOGICAL :: do_distribute, found, is_root_rank, &
symmetric
REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: recv_buffer, reorder_recv_buffer, &
send_buffer
REAL(kind=dp), DIMENSION(:, :), POINTER :: sm_block, sm_block_merged
TYPE(mp_request_type), ALLOCATABLE, DIMENSION(:) :: requests
CALL timeset(routineN, handle)
matrix_siesta(:) = 0.0_dp
mepos = para_env%mepos
nprocs = para_env%num_pe
do_distribute = siesta_struct%gather_root < 0
is_root_rank = siesta_struct%gather_root == mepos
CALL dbcsr_get_info(matrix=matrix_dbcsr_kp(1)%matrix, &
nblkrows_total=nrows_blk, nblkcols_total=ncols_blk, &
row_blk_size=row_blk_size, col_blk_size=col_blk_size, &
row_blk_offset=row_blk_offset, col_blk_offset=col_blk_offset)
symmetric = siesta_struct%symmetric
! number of locally stored SIESTA non-zero matrix elements
n_nonzero_elements_siesta = SUM(siesta_struct%nelements_per_proc(:, nelements_dbcsr_recv))
! number of locally stored DBCSR non-zero matrix elements
n_nonzero_elements_dbcsr = SUM(siesta_struct%nelements_per_proc(:, nelements_dbcsr_send))
! number of concurrent MPI isend / irecv operations
nrequests_recv = get_number_of_mpi_sendrecv_requests(mepos, siesta_struct%nelements_per_proc(:, nelements_dbcsr_recv), &
max_mpi_packet_size_dp)
nrequests_total = get_number_of_mpi_sendrecv_requests(mepos, siesta_struct%nelements_per_proc(:, nelements_dbcsr_send), &
max_mpi_packet_size_dp) + nrequests_recv
IF (nrequests_total > 0) THEN
! allocate MPI-related arrays. request_tag is not actually needed, as MPI standard guarantees the order of
! peer-to-peer messages with the same tag between same processes
ALLOCATE (requests(nrequests_total))
ALLOCATE (peer_rank(nrequests_total), request_tag(nrequests_total))
ALLOCATE (offset_per_request(nrequests_total), nelements_per_request(nrequests_total))
!requests(:) = mp_request_null
! split large messages into a number of smaller messages. It is not really needed
! unless we are going to send > 2^31 matrix elements per MPI request
IF (nrequests_recv > 0) THEN
CALL assign_nonzero_elements_to_requests(offset_per_request(1:nrequests_recv), &
nelements_per_request(1:nrequests_recv), &
peer_rank(1:nrequests_recv), &
request_tag(1:nrequests_recv), &
mepos, &
siesta_struct%nelements_per_proc(:, nelements_dbcsr_recv), &
max_mpi_packet_size_dp)
END IF
IF (nrequests_total > nrequests_recv) THEN
CALL assign_nonzero_elements_to_requests(offset_per_request(nrequests_recv + 1:nrequests_total), &
nelements_per_request(nrequests_recv + 1:nrequests_total), &
peer_rank(nrequests_recv + 1:nrequests_total), &
request_tag(nrequests_recv + 1:nrequests_total), &
mepos, &
siesta_struct%nelements_per_proc(:, nelements_dbcsr_send), &
max_mpi_packet_size_dp)
END IF
END IF
! point-to-point recv/send can be replaced with alltoallv, if data to distribute per rank is < 2^31 elements (should be OK due to SMEAGOL limitation)
! in principle, it is possible to overcome this limit by using derived datatypes (which will require additional wrapper functions, indeed).
!
! pre-post non-blocking receive operations
IF (n_nonzero_elements_siesta > 0) THEN
ALLOCATE (recv_buffer(n_nonzero_elements_siesta))
END IF
DO irequest = 1, nrequests_recv
CALL para_env%irecv(recv_buffer(offset_per_request(irequest) + 1: &
offset_per_request(irequest) + nelements_per_request(irequest)), &
peer_rank(irequest), requests(irequest), request_tag(irequest))
END DO
! pack local DBCSR non-zero matrix elements ordering by their target parallel process
ALLOCATE (offset_per_proc(0:nprocs - 1), n_packed_elements_per_proc(0:nprocs - 1))
offset_per_proc(0) = 0
DO iproc = 1, nprocs - 1
offset_per_proc(iproc) = offset_per_proc(iproc - 1) + siesta_struct%nelements_per_proc(iproc - 1, nelements_dbcsr_send)
END DO
n_packed_elements_per_proc(:) = 0
! number of local neighbour-list nodes and offset of the first local neighbour-list node
nnodes_proc = siesta_struct%nnodes_per_proc(mepos)
!node_offset = SUM(siesta_struct%nnodes_per_proc(0:mepos)) - siesta_struct%nnodes_per_proc(mepos)
node_offset = SUM(siesta_struct%nnodes_per_proc(0:mepos)) - nnodes_proc
! if do_distribute == .FALSE., send all matrix elements to MPI process with rank gather_root
! in case of do_distribute == .TRUE., iproc is determined by calling WhichNodeOrb()
iproc = siesta_struct%gather_root
IF (n_nonzero_elements_dbcsr > 0) THEN
ALLOCATE (send_buffer(n_nonzero_elements_dbcsr))
send_buffer(:) = 0.0_dp
! iterate over locally-stored DBCSR matrix blocks.
! inode_proc is the target parallel process (where data are going to be sent)
image_ind_offset = 0
DO inode_proc = 1, nnodes_proc
n_image_ind = siesta_struct%n_dbcsr_cell_images_to_merge(inode_proc)
IF (n_image_ind > 0) THEN
inode = node_offset + inode_proc
irow_blk = siesta_struct%nl_repl(neighbor_list_iatom_index, inode)
icol_blk = siesta_struct%nl_repl(neighbor_list_jatom_index, inode)
CPASSERT(irow_blk <= icol_blk .OR. (.NOT. symmetric))
image_siesta = siesta_struct%nl_repl(neighbor_list_siesta_image_index, inode)
image_siesta_transp = siesta_struct%nl_repl(neighbor_list_siesta_transp_image_index, inode)
nrows_local = row_blk_size(irow_blk)
ncols_local = col_blk_size(icol_blk)
first_row_minus_one = row_blk_offset(irow_blk) - 1
first_col_minus_one = col_blk_offset(icol_blk) - 1
! merging cell images along transport direction
IF (n_image_ind == 1) THEN
! the most common case. Nothing to merge, so there is no need to allocate memory for a merged block
image_dbcsr = siesta_struct%dbcsr_cell_image_to_merge(image_ind_offset + 1)
CALL dbcsr_get_block_p(matrix=matrix_dbcsr_kp(image_dbcsr)%matrix, &
row=irow_blk, col=icol_blk, block=sm_block_merged, found=found)
CPASSERT(found)
ELSE ! n_image_ind > 1
ALLOCATE (sm_block_merged(nrows_local, ncols_local))
DO image_ind = 1, n_image_ind
image_dbcsr = siesta_struct%dbcsr_cell_image_to_merge(image_ind + image_ind_offset)
CALL dbcsr_get_block_p(matrix=matrix_dbcsr_kp(image_dbcsr)%matrix, &
row=irow_blk, col=icol_blk, block=sm_block, found=found)
CPASSERT(found)
sm_block_merged(1:nrows_local, 1:ncols_local) = sm_block_merged(1:nrows_local, 1:ncols_local) + &
sm_block(1:nrows_local, 1:ncols_local)
END DO
END IF
! pack matrix elements for the 'normal' SIESTA matrix block
IF (image_siesta > 0) THEN
DO irow_local = 1, nrows_local
IF (do_distribute) THEN
#if defined(__SMEAGOL)
CALL WhichNodeOrb(irow_local + first_row_minus_one, nprocs, iproc)
#else
CALL cp_abort(__LOCATION__, &
"CP2K was compiled with no SMEAGOL support.")
#endif
END IF
! CPASSERT
IF (debug_this_module) THEN
CPASSERT(iproc >= 0 .AND. iproc < nprocs)
IF (n_packed_elements_per_proc(iproc) + ncols_local > &
siesta_struct%nelements_per_proc(iproc, nelements_dbcsr_send)) THEN
CALL cp__a(__SHORT_FILE__, __LINE__)
END IF
END IF
offset_send_mepos = offset_per_proc(iproc) + n_packed_elements_per_proc(iproc)
send_buffer(offset_send_mepos + 1:offset_send_mepos + ncols_local) = sm_block_merged(irow_local, 1:ncols_local)
n_packed_elements_per_proc(iproc) = n_packed_elements_per_proc(iproc) + ncols_local
END DO
END IF
! pack matrix elements of the transposed SIESTA matrix block
IF (image_siesta_transp > 0) THEN
DO icol_local = 1, ncols_local
IF (do_distribute) THEN
#if defined(__SMEAGOL)
CALL WhichNodeOrb(icol_local + first_col_minus_one, nprocs, iproc) ! iproc_orb
#else
CALL cp_abort(__LOCATION__, &
"CP2K was compiled with no SMEAGOL support.")
#endif
END IF
! CPASSERT
IF (debug_this_module) THEN
CPASSERT(iproc >= 0 .AND. iproc < nprocs)
IF (n_packed_elements_per_proc(iproc) + nrows_local > &
siesta_struct%nelements_per_proc(iproc, nelements_dbcsr_send)) THEN
CALL cp__a(__SHORT_FILE__, __LINE__)
END IF
END IF
offset_send_mepos = offset_per_proc(iproc) + n_packed_elements_per_proc(iproc)
send_buffer(offset_send_mepos + 1:offset_send_mepos + nrows_local) = sm_block_merged(1:nrows_local, icol_local)
n_packed_elements_per_proc(iproc) = n_packed_elements_per_proc(iproc) + nrows_local
END DO
END IF
IF (n_image_ind > 1) THEN
DEALLOCATE (sm_block_merged)
END IF
END IF
image_ind_offset = image_ind_offset + siesta_struct%n_dbcsr_cell_images_to_merge(inode_proc)
END DO
IF (debug_this_module) THEN
DO iproc = 0, nprocs - 1
IF (n_packed_elements_per_proc(iproc) /= siesta_struct%nelements_per_proc(iproc, nelements_dbcsr_send)) THEN
CALL cp__a(__SHORT_FILE__, __LINE__)
END IF
END DO
END IF
! send packed data to other parallel processes
DO irequest = nrequests_recv + 1, nrequests_total
CALL para_env%isend(send_buffer(offset_per_request(irequest) + 1: &
offset_per_request(irequest) + nelements_per_request(irequest)), &
peer_rank(irequest), requests(irequest), request_tag(irequest))
END DO
! copy data locally that stay on the same process.
IF (mepos > 0) THEN
offset_recv_mepos = SUM(siesta_struct%nelements_per_proc(0:mepos - 1, nelements_dbcsr_recv))
ELSE
offset_recv_mepos = 0
END IF
offset_send_mepos = offset_per_proc(mepos)
IF (debug_this_module) THEN
IF (n_packed_elements_per_proc(mepos) /= siesta_struct%nelements_per_proc(mepos, nelements_dbcsr_recv)) THEN
CALL cp__a(__SHORT_FILE__, __LINE__)
END IF
END IF
IF (n_packed_elements_per_proc(mepos) > 0) THEN
recv_buffer(offset_recv_mepos + 1:offset_recv_mepos + n_packed_elements_per_proc(mepos)) = &
send_buffer(offset_send_mepos + 1:offset_send_mepos + n_packed_elements_per_proc(mepos))
END IF
END IF
IF (nrequests_total > 0) THEN
! wait for pending isend/irecv requests
CALL mp_waitall(requests)
DEALLOCATE (nelements_per_request, offset_per_request, peer_rank, requests, request_tag)
END IF
! release send buffers
IF (ALLOCATED(send_buffer)) DEALLOCATE (send_buffer)
DEALLOCATE (offset_per_proc, n_packed_elements_per_proc)
! non-zero matrix elements in 'recv_buffer' array are grouped by their source MPI rank, local row index, and column index (in this order).
! Reorder the matrix elements ('reorder_recv_buffer') so they are grouped by their local row index, source MPI rank, and column index.
! (column indices are in the ascending order within each (row index, source MPI rank) block).
! The array 'packed_index' allows mapping matrix element between these intermediate order and SIESTA order : local row index, column index
IF (n_nonzero_elements_siesta > 0) THEN
ALLOCATE (reorder_recv_buffer(n_nonzero_elements_siesta))
ALLOCATE (next_nonzero_element_offset(siesta_struct%nrows))
next_nonzero_element_offset(:) = 0
offset_recv_mepos = 0
DO inode = 1, SIZE(siesta_struct%nl_repl, 2)
irow_blk = siesta_struct%nl_repl(neighbor_list_iatom_index, inode)
icol_blk = siesta_struct%nl_repl(neighbor_list_jatom_index, inode)
CPASSERT(irow_blk <= icol_blk .OR. (.NOT. symmetric))
image_siesta = siesta_struct%nl_repl(neighbor_list_siesta_image_index, inode)
image_siesta_transp = siesta_struct%nl_repl(neighbor_list_siesta_transp_image_index, inode)
nrows_local = row_blk_size(irow_blk)
ncols_local = col_blk_size(icol_blk)
first_row_minus_one = row_blk_offset(irow_blk) - 1
first_col_minus_one = col_blk_offset(icol_blk) - 1
! normal block
IF (image_siesta > 0) THEN
DO irow_local = 1, nrows_local
IF (do_distribute) THEN
#if defined(__SMEAGOL)
CALL GlobalToLocalOrb(irow_local + first_row_minus_one, mepos, nprocs, irow_proc)
#else
CALL cp_abort(__LOCATION__, &
"CP2K was compiled with no SMEAGOL support.")
#endif
ELSE
IF (is_root_rank) THEN
irow_proc = irow_local + first_row_minus_one
ELSE
irow_proc = 0
END IF
END IF
IF (irow_proc > 0) THEN
offset_send_mepos = siesta_struct%row_offset(irow_proc) + next_nonzero_element_offset(irow_proc)
reorder_recv_buffer(offset_send_mepos + 1:offset_send_mepos + ncols_local) = &
recv_buffer(offset_recv_mepos + 1:offset_recv_mepos + ncols_local)
offset_recv_mepos = offset_recv_mepos + ncols_local
next_nonzero_element_offset(irow_proc) = next_nonzero_element_offset(irow_proc) + ncols_local
END IF
END DO
END IF
! transposed block
IF (image_siesta_transp > 0) THEN
DO icol_local = 1, ncols_local
IF (do_distribute) THEN
#if defined(__SMEAGOL)
CALL GlobalToLocalOrb(icol_local + first_col_minus_one, mepos, nprocs, irow_proc)
#else
CALL cp_abort(__LOCATION__, &
"CP2K was compiled with no SMEAGOL support.")
#endif
ELSE
IF (is_root_rank) THEN
irow_proc = icol_local + first_col_minus_one
ELSE
irow_proc = 0
END IF
END IF
IF (irow_proc > 0) THEN
offset_send_mepos = siesta_struct%row_offset(irow_proc) + next_nonzero_element_offset(irow_proc)
reorder_recv_buffer(offset_send_mepos + 1:offset_send_mepos + nrows_local) = &
recv_buffer(offset_recv_mepos + 1:offset_recv_mepos + nrows_local)
offset_recv_mepos = offset_recv_mepos + nrows_local
next_nonzero_element_offset(irow_proc) = next_nonzero_element_offset(irow_proc) + nrows_local
END IF
END DO
END IF
END DO
IF (debug_this_module) THEN
DO irow_local = 1, siesta_struct%nrows
IF (siesta_struct%n_nonzero_cols(irow_local) /= next_nonzero_element_offset(irow_local)) THEN
CALL cp__a(__SHORT_FILE__, __LINE__)
END IF
END DO
END IF
DEALLOCATE (next_nonzero_element_offset)
DEALLOCATE (recv_buffer)
! Map non-zero matrix element between the intermediate order and SIESTA order
DO irow_local = 1, siesta_struct%nrows
offset_recv_mepos = siesta_struct%row_offset(irow_local)
DO icol_local = 1, siesta_struct%n_nonzero_cols(irow_local)
matrix_siesta(offset_recv_mepos + icol_local) = &
reorder_recv_buffer(offset_recv_mepos + siesta_struct%packed_index(offset_recv_mepos + icol_local))
END DO
END DO
DEALLOCATE (reorder_recv_buffer)
END IF
CALL timestop(handle)
END SUBROUTINE convert_dbcsr_to_distributed_siesta
! **************************************************************************************************
!> \brief Convert matrix from DBCSR to sparse SIESTA format.
!> \param matrix_dbcsr_kp DBCSR matrix [out]. The matrix is declared as INTENT(in) as pointers to
!> dbcsr matrices remain intact. However we have intention to update matrix elements
!> \param matrix_siesta matrix in SIESTA format [in]
!> \param siesta_struct structure to map matrix blocks between formats
!> \param para_env MPI parallel environment
! **************************************************************************************************
SUBROUTINE convert_distributed_siesta_to_dbcsr(matrix_dbcsr_kp, matrix_siesta, siesta_struct, para_env)
TYPE(dbcsr_p_type), DIMENSION(:), INTENT(in) :: matrix_dbcsr_kp
REAL(kind=dp), DIMENSION(:), INTENT(in) :: matrix_siesta
TYPE(siesta_distrib_csc_struct_type), INTENT(in) :: siesta_struct
TYPE(mp_para_env_type), INTENT(in), POINTER :: para_env
CHARACTER(len=*), PARAMETER :: routineN = 'convert_distributed_siesta_to_dbcsr'
INTEGER :: first_col_minus_one, first_row_minus_one, handle, icol_blk, icol_local, &
image_dbcsr, image_siesta, image_siesta_transp, inode, inode_proc, iproc, irequest, &
irow_blk, irow_local, irow_proc, mepos, n_image_ind, ncols_blk, ncols_local, nnodes_proc, &
node_offset, nprocs, nrequests_recv, nrequests_total, nrows_blk, nrows_local
INTEGER(kind=int_8) :: n_nonzero_elements_dbcsr, &
n_nonzero_elements_siesta, &
offset_recv_mepos, offset_send_mepos
INTEGER(kind=int_8), ALLOCATABLE, DIMENSION(:) :: n_packed_elements_per_proc, &
nelements_per_request, &
offset_per_proc, offset_per_request
INTEGER, ALLOCATABLE, DIMENSION(:) :: next_nonzero_element_offset, peer_rank, &
request_tag
INTEGER, DIMENSION(:), POINTER :: col_blk_offset, col_blk_size, &
row_blk_offset, row_blk_size
LOGICAL :: do_distribute, found, is_root_rank, &
symmetric
REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: recv_buffer, reorder_send_buffer, &
send_buffer
REAL(kind=dp), DIMENSION(:, :), POINTER :: sm_block
TYPE(mp_request_type), ALLOCATABLE, DIMENSION(:) :: requests
CALL timeset(routineN, handle)
DO image_dbcsr = 1, SIZE(matrix_dbcsr_kp)
CALL dbcsr_set(matrix_dbcsr_kp(image_dbcsr)%matrix, 0.0_dp)
END DO
mepos = para_env%mepos
nprocs = para_env%num_pe
do_distribute = siesta_struct%gather_root < 0
is_root_rank = siesta_struct%gather_root == mepos
CALL dbcsr_get_info(matrix=matrix_dbcsr_kp(1)%matrix, &
nblkrows_total=nrows_blk, nblkcols_total=ncols_blk, &
row_blk_size=row_blk_size, col_blk_size=col_blk_size, &
row_blk_offset=row_blk_offset, col_blk_offset=col_blk_offset)
symmetric = siesta_struct%symmetric
n_nonzero_elements_siesta = SUM(siesta_struct%nelements_per_proc(:, nelements_dbcsr_recv))
n_nonzero_elements_dbcsr = SUM(siesta_struct%nelements_per_proc(:, nelements_dbcsr_send))
nrequests_recv = get_number_of_mpi_sendrecv_requests(mepos, siesta_struct%nelements_per_proc(:, nelements_dbcsr_send), &
max_mpi_packet_size_dp)
nrequests_total = get_number_of_mpi_sendrecv_requests(mepos, siesta_struct%nelements_per_proc(:, nelements_dbcsr_recv), &
max_mpi_packet_size_dp) + nrequests_recv
IF (nrequests_total > 0) THEN
ALLOCATE (requests(nrequests_total))
ALLOCATE (peer_rank(nrequests_total), request_tag(nrequests_total))
ALLOCATE (offset_per_request(nrequests_total), nelements_per_request(nrequests_total))
!requests(:) = mp_request_null
IF (nrequests_recv > 0) THEN
CALL assign_nonzero_elements_to_requests(offset_per_request(1:nrequests_recv), &
nelements_per_request(1:nrequests_recv), &
peer_rank(1:nrequests_recv), &
request_tag(1:nrequests_recv), &
mepos, &
siesta_struct%nelements_per_proc(:, nelements_dbcsr_send), &
max_mpi_packet_size_dp)
END IF
IF (nrequests_total > nrequests_recv) THEN
CALL assign_nonzero_elements_to_requests(offset_per_request(nrequests_recv + 1:nrequests_total), &
nelements_per_request(nrequests_recv + 1:nrequests_total), &
peer_rank(nrequests_recv + 1:nrequests_total), &
request_tag(nrequests_recv + 1:nrequests_total), &
mepos, &
siesta_struct%nelements_per_proc(:, nelements_dbcsr_recv), &
max_mpi_packet_size_dp)
END IF
END IF
IF (n_nonzero_elements_dbcsr > 0) THEN
ALLOCATE (recv_buffer(n_nonzero_elements_dbcsr))
END IF
DO irequest = 1, nrequests_recv
CALL para_env%irecv(recv_buffer(offset_per_request(irequest) + 1: &
offset_per_request(irequest) + nelements_per_request(irequest)), &
peer_rank(irequest), requests(irequest), request_tag(irequest))
END DO
ALLOCATE (offset_per_proc(0:nprocs - 1), n_packed_elements_per_proc(0:nprocs - 1))
offset_per_proc(0) = 0
DO iproc = 1, nprocs - 1
offset_per_proc(iproc) = offset_per_proc(iproc - 1) + siesta_struct%nelements_per_proc(iproc - 1, nelements_dbcsr_send)
END DO
n_packed_elements_per_proc(:) = 0
IF (mepos > 0) THEN
node_offset = SUM(siesta_struct%nnodes_per_proc(0:mepos - 1))
ELSE
node_offset = 0
END IF
nnodes_proc = siesta_struct%nnodes_per_proc(mepos)
IF (n_nonzero_elements_siesta > 0) THEN
ALLOCATE (send_buffer(n_nonzero_elements_siesta))
ALLOCATE (reorder_send_buffer(n_nonzero_elements_siesta))
DO irow_local = 1, siesta_struct%nrows
offset_send_mepos = siesta_struct%row_offset(irow_local)
DO icol_local = 1, siesta_struct%n_nonzero_cols(irow_local)
reorder_send_buffer(offset_send_mepos + siesta_struct%packed_index(offset_send_mepos + icol_local)) = &
matrix_siesta(offset_send_mepos + icol_local)
END DO
END DO
ALLOCATE (next_nonzero_element_offset(siesta_struct%nrows))
next_nonzero_element_offset(:) = 0
offset_send_mepos = 0
DO inode = 1, SIZE(siesta_struct%nl_repl, 2)
irow_blk = siesta_struct%nl_repl(neighbor_list_iatom_index, inode)
icol_blk = siesta_struct%nl_repl(neighbor_list_jatom_index, inode)
CPASSERT(irow_blk <= icol_blk .OR. (.NOT. symmetric))
image_siesta = siesta_struct%nl_repl(neighbor_list_siesta_image_index, inode)
image_siesta_transp = siesta_struct%nl_repl(neighbor_list_siesta_transp_image_index, inode)
nrows_local = row_blk_size(irow_blk)
ncols_local = col_blk_size(icol_blk)
first_row_minus_one = row_blk_offset(irow_blk) - 1
first_col_minus_one = col_blk_offset(icol_blk) - 1
IF (image_siesta > 0) THEN
DO irow_local = 1, nrows_local
IF (do_distribute) THEN
#if defined(__SMEAGOL)
CALL GlobalToLocalOrb(irow_local + first_row_minus_one, mepos, nprocs, irow_proc)
#else
CALL cp_abort(__LOCATION__, &
"CP2K was compiled with no SMEAGOL support.")
#endif
ELSE
IF (is_root_rank) THEN
irow_proc = irow_local + first_row_minus_one
ELSE
irow_proc = 0
END IF
END IF
IF (irow_proc > 0) THEN
offset_recv_mepos = siesta_struct%row_offset(irow_proc) + next_nonzero_element_offset(irow_proc)
send_buffer(offset_send_mepos + 1:offset_send_mepos + ncols_local) = &
reorder_send_buffer(offset_recv_mepos + 1:offset_recv_mepos + ncols_local)
offset_send_mepos = offset_send_mepos + ncols_local
next_nonzero_element_offset(irow_proc) = next_nonzero_element_offset(irow_proc) + ncols_local
END IF
END DO
END IF
! transposed block
IF (image_siesta_transp > 0) THEN
DO icol_local = 1, ncols_local
IF (do_distribute) THEN
#if defined(__SMEAGOL)
CALL GlobalToLocalOrb(icol_local + first_col_minus_one, mepos, nprocs, irow_proc)
#else
CALL cp_abort(__LOCATION__, &
"CP2K was compiled with no SMEAGOL support.")
#endif
ELSE
IF (is_root_rank) THEN
irow_proc = icol_local + first_col_minus_one
ELSE
irow_proc = 0
END IF
END IF
IF (irow_proc > 0) THEN
offset_recv_mepos = siesta_struct%row_offset(irow_proc) + next_nonzero_element_offset(irow_proc)
send_buffer(offset_send_mepos + 1:offset_send_mepos + nrows_local) = &
reorder_send_buffer(offset_recv_mepos + 1:offset_recv_mepos + nrows_local)
offset_send_mepos = offset_send_mepos + nrows_local
next_nonzero_element_offset(irow_proc) = next_nonzero_element_offset(irow_proc) + nrows_local
END IF
END DO
END IF
END DO
IF (debug_this_module) THEN
DO irow_local = 1, siesta_struct%nrows
IF (siesta_struct%n_nonzero_cols(irow_local) /= next_nonzero_element_offset(irow_local)) THEN
CALL cp__a(__SHORT_FILE__, __LINE__)
END IF
END DO
END IF
DEALLOCATE (next_nonzero_element_offset)
DEALLOCATE (reorder_send_buffer)
DO irequest = nrequests_recv + 1, nrequests_total
CALL para_env%isend(send_buffer(offset_per_request(irequest) + 1: &
offset_per_request(irequest) + nelements_per_request(irequest)), &
peer_rank(irequest), requests(irequest), request_tag(irequest))
END DO
! copy data locally that stay on the same process.
IF (mepos > 0) THEN
offset_send_mepos = SUM(siesta_struct%nelements_per_proc(0:mepos - 1, nelements_dbcsr_recv))
ELSE