(git:374b731)
Loading...
Searching...
No Matches
dbt_test.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief General methods for testing DBT tensors.
10!> \author Patrick Seewald
11! **************************************************************************************************
12MODULE dbt_test
13
14
15 USE dbt_tas_base, ONLY: dbt_tas_info
16 USE dbm_tests, ONLY: generate_larnv_seed
17 USE dbt_methods, ONLY: &
18 dbt_copy, dbt_get_block, dbt_iterator_type, dbt_iterator_blocks_left, &
19 dbt_iterator_next_block, dbt_iterator_start, dbt_iterator_stop, &
20 dbt_reserve_blocks, dbt_get_stored_coordinates, dbt_put_block, &
21 dbt_contract, dbt_inverse_order
22 USE dbt_block, ONLY: block_nd
23 USE dbt_types, ONLY: &
24 dbt_create, dbt_destroy, dbt_type, dbt_distribution_type, &
25 dbt_distribution_destroy, dims_tensor, ndims_tensor, dbt_distribution_new, &
26 mp_environ_pgrid, dbt_pgrid_type, dbt_pgrid_create, dbt_pgrid_destroy, dbt_get_info, &
27 dbt_default_distvec
28 USE dbt_io, ONLY: &
29 dbt_write_blocks, dbt_write_block_indices
30 USE kinds, ONLY: dp, default_string_length, int_8, dp
31 USE dbt_allocate_wrap, ONLY: allocate_any
32 USE dbt_index, ONLY: &
33 combine_tensor_index, get_2d_indices_tensor, dbt_get_mapping_info
34 USE dbt_tas_test, ONLY: dbt_tas_checksum
35 USE message_passing, ONLY: mp_comm_type
36
37#include "../base/base_uses.f90"
38
39 IMPLICIT NONE
40 PRIVATE
41 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'dbt_test'
42
43 PUBLIC :: &
44 dbt_setup_test_tensor, &
45 dbt_contract_test, &
46 dbt_test_formats, &
47 dbt_checksum, &
48 dbt_reset_randmat_seed
49
50 INTERFACE dist_sparse_tensor_to_repl_dense_array
51 MODULE PROCEDURE dist_sparse_tensor_to_repl_dense_2d_array
52 MODULE PROCEDURE dist_sparse_tensor_to_repl_dense_3d_array
53 MODULE PROCEDURE dist_sparse_tensor_to_repl_dense_4d_array
54 END INTERFACE
55
56 INTEGER, SAVE :: randmat_counter = 0
57 INTEGER, PARAMETER, PRIVATE :: rand_seed_init = 12341313
58
59CONTAINS
60
61! **************************************************************************************************
62!> \brief check if two (arbitrarily mapped and distributed) tensors are equal.
63!> \author Patrick Seewald
64! **************************************************************************************************
65 FUNCTION dbt_equal(tensor1, tensor2)
66 TYPE(dbt_type), INTENT(INOUT) :: tensor1, tensor2
67 LOGICAL :: dbt_equal
68
69 TYPE(dbt_type) :: tensor2_tmp
70 TYPE(dbt_iterator_type) :: iter
71 TYPE(block_nd) :: blk_data1, blk_data2
72 INTEGER, DIMENSION(ndims_tensor(tensor1)) :: blk_size, ind_nd
73 LOGICAL :: found
74
75 ! create a copy of tensor2 that has exact same data format as tensor1
76 CALL dbt_create(tensor1, tensor2_tmp)
77
78 CALL dbt_reserve_blocks(tensor1, tensor2_tmp)
79 CALL dbt_copy(tensor2, tensor2_tmp)
80
81 dbt_equal = .true.
82
83!$OMP PARALLEL DEFAULT(NONE) SHARED(tensor1,tensor2_tmp,dbt_equal) &
84!$OMP PRIVATE(iter,ind_nd, blk_size,blk_data1,blk_data2,found)
85 CALL dbt_iterator_start(iter, tensor1)
86
87 DO WHILE (dbt_iterator_blocks_left(iter))
88 CALL dbt_iterator_next_block(iter, ind_nd, blk_size=blk_size)
89 CALL dbt_get_block(tensor1, ind_nd, blk_data1, found)
90 IF (.NOT. found) cpabort("Tensor block 1 not found")
91 CALL dbt_get_block(tensor2_tmp, ind_nd, blk_data2, found)
92 IF (.NOT. found) cpabort("Tensor block 2 not found")
93
94 IF (.NOT. blocks_equal(blk_data1, blk_data2)) THEN
95!$OMP CRITICAL
96 dbt_equal = .false.
97!$OMP END CRITICAL
98 END IF
99 END DO
100
101 CALL dbt_iterator_stop(iter)
102!$OMP END PARALLEL
103
104 CALL dbt_destroy(tensor2_tmp)
105 END FUNCTION
106
107! **************************************************************************************************
108!> \brief check if two blocks are equal
109!> \author Patrick Seewald
110! **************************************************************************************************
111 PURE FUNCTION blocks_equal(block1, block2)
112 TYPE(block_nd), INTENT(IN) :: block1, block2
113 LOGICAL :: blocks_equal
114
115 blocks_equal = maxval(abs(block1%blk - block2%blk)) .LT. 1.0e-12_dp
116
117 END FUNCTION
118
119! **************************************************************************************************
120!> \brief Compute factorial
121!> \author Patrick Seewald
122! **************************************************************************************************
123 PURE FUNCTION factorial(n)
124 INTEGER, INTENT(IN) :: n
125 INTEGER :: k
126 INTEGER :: factorial
127 factorial = product((/(k, k=1, n)/))
128 END FUNCTION
129
130! **************************************************************************************************
131!> \brief Compute all permutations p of (1, 2, ..., n)
132!> \author Patrick Seewald
133! **************************************************************************************************
134 SUBROUTINE permute(n, p)
135 INTEGER, INTENT(IN) :: n
136 INTEGER :: i, c
137 INTEGER, DIMENSION(n) :: pp
138 INTEGER, DIMENSION(n, factorial(n)), INTENT(OUT) :: p
139
140 pp = [(i, i=1, n)]
141 c = 1
142 CALL perm(1)
143 CONTAINS
144 RECURSIVE SUBROUTINE perm(i)
145 INTEGER, INTENT(IN) :: i
146 INTEGER :: j, t
147 IF (i == n) THEN
148 p(:, c) = pp(:)
149 c = c + 1
150 ELSE
151 DO j = i, n
152 t = pp(i)
153 pp(i) = pp(j)
154 pp(j) = t
155 call perm(i + 1)
156 t = pp(i)
157 pp(i) = pp(j)
158 pp(j) = t
159 END DO
160 END IF
161 END SUBROUTINE
162 END SUBROUTINE
163
164! **************************************************************************************************
165!> \brief Test equivalence of all tensor formats, using a random distribution.
166!> \param blk_size_i block sizes along respective dimension
167!> \param blk_ind_i index along respective dimension of non-zero blocks
168!> \param ndims tensor rank
169!> \param unit_nr output unit, needs to be a valid unit number on all mpi ranks
170!> \param verbose if .TRUE., print all tensor blocks
171!> \author Patrick Seewald
172! **************************************************************************************************
173 SUBROUTINE dbt_test_formats(ndims, mp_comm, unit_nr, verbose, &
174 blk_size_1, blk_size_2, blk_size_3, blk_size_4, &
175 blk_ind_1, blk_ind_2, blk_ind_3, blk_ind_4)
176 INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: blk_size_1, blk_size_2, blk_size_3, blk_size_4
177 INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: blk_ind_1, blk_ind_2, blk_ind_3, blk_ind_4
178 INTEGER, INTENT(IN) :: ndims
179 INTEGER, INTENT(IN) :: unit_nr
180 LOGICAL, INTENT(IN) :: verbose
181 TYPE(mp_comm_type), INTENT(IN) :: mp_comm
182 TYPE(dbt_distribution_type) :: dist1, dist2
183 TYPE(dbt_type) :: tensor1, tensor2
184 INTEGER :: isep, iblk
185 INTEGER, DIMENSION(:), ALLOCATABLE :: dist1_1, dist1_2, dist1_3, dist1_4, &
186 dist2_1, dist2_2, dist2_3, dist2_4
187 INTEGER :: nblks, imap
188 INTEGER, DIMENSION(ndims) :: pdims, myploc
189 LOGICAL :: eql
190 INTEGER :: iperm, idist, icount
191 INTEGER, DIMENSION(:), ALLOCATABLE :: map1, map2, map1_ref, map2_ref
192 INTEGER, DIMENSION(ndims, factorial(ndims)) :: perm
193 INTEGER :: io_unit
194 INTEGER :: mynode
195 TYPE(dbt_pgrid_type) :: comm_nd
196 CHARACTER(LEN=default_string_length) :: tensor_name
197
198 ! Process grid
199 pdims(:) = 0
200 CALL dbt_pgrid_create(mp_comm, pdims, comm_nd)
201 mynode = mp_comm%mepos
202
203 io_unit = 0
204 IF (mynode .EQ. 0) io_unit = unit_nr
205
206 CALL permute(ndims, perm)
207 ALLOCATE (map1_ref, source=perm(1:ndims/2, 1))
208 ALLOCATE (map2_ref, source=perm(ndims/2 + 1:ndims, 1))
209
210 IF (io_unit > 0) THEN
211 WRITE (io_unit, *)
212 WRITE (io_unit, '(A)') repeat("-", 80)
213 WRITE (io_unit, '(A,1X,I1)') "Testing matrix representations of tensor rank", ndims
214 WRITE (io_unit, '(A)') repeat("-", 80)
215 WRITE (io_unit, '(A)') "Block sizes:"
216
217 IF (ndims >= 1) THEN
218 WRITE (io_unit, '(T4,A,1X,I1,A,1X)', advance='no') 'Dim', 1, ':'
219 DO iblk = 1, SIZE(blk_size_1)
220 WRITE (io_unit, '(I2,1X)', advance='no') blk_size_1(iblk)
221 END DO
222 WRITE (io_unit, *)
223 END IF
224 IF (ndims >= 2) THEN
225 WRITE (io_unit, '(T4,A,1X,I1,A,1X)', advance='no') 'Dim', 2, ':'
226 DO iblk = 1, SIZE(blk_size_2)
227 WRITE (io_unit, '(I2,1X)', advance='no') blk_size_2(iblk)
228 END DO
229 WRITE (io_unit, *)
230 END IF
231 IF (ndims >= 3) THEN
232 WRITE (io_unit, '(T4,A,1X,I1,A,1X)', advance='no') 'Dim', 3, ':'
233 DO iblk = 1, SIZE(blk_size_3)
234 WRITE (io_unit, '(I2,1X)', advance='no') blk_size_3(iblk)
235 END DO
236 WRITE (io_unit, *)
237 END IF
238 IF (ndims >= 4) THEN
239 WRITE (io_unit, '(T4,A,1X,I1,A,1X)', advance='no') 'Dim', 4, ':'
240 DO iblk = 1, SIZE(blk_size_4)
241 WRITE (io_unit, '(I2,1X)', advance='no') blk_size_4(iblk)
242 END DO
243 WRITE (io_unit, *)
244 END IF
245
246 WRITE (io_unit, '(A)') "Non-zero blocks:"
247 DO iblk = 1, SIZE(blk_ind_1)
248 IF (ndims == 2) THEN
249 WRITE (io_unit, '(T4,A, I3, A, 2I3, 1X, A)') &
250 'Block', iblk, ': (', blk_ind_1(iblk), blk_ind_2(iblk), ')'
251 END IF
252 IF (ndims == 3) THEN
253 WRITE (io_unit, '(T4,A, I3, A, 3I3, 1X, A)') &
254 'Block', iblk, ': (', blk_ind_1(iblk), blk_ind_2(iblk), blk_ind_3(iblk), ')'
255 END IF
256 IF (ndims == 4) THEN
257 WRITE (io_unit, '(T4,A, I3, A, 4I3, 1X, A)') &
258 'Block', iblk, ': (', blk_ind_1(iblk), blk_ind_2(iblk), blk_ind_3(iblk), blk_ind_4(iblk), ')'
259 END IF
260 END DO
261
262 WRITE (io_unit, *)
263 WRITE (io_unit, '(A,1X)', advance='no') "Reference map:"
264 WRITE (io_unit, '(A1,1X)', advance='no') "("
265 DO imap = 1, SIZE(map1_ref)
266 WRITE (io_unit, '(I1,1X)', advance='no') map1_ref(imap)
267 END DO
268 WRITE (io_unit, '(A1,1X)', advance='no') "|"
269 DO imap = 1, SIZE(map2_ref)
270 WRITE (io_unit, '(I1,1X)', advance='no') map2_ref(imap)
271 END DO
272 WRITE (io_unit, '(A1)') ")"
273
274 END IF
275
276 icount = 0
277 DO iperm = 1, factorial(ndims)
278 DO isep = 1, ndims - 1
279 icount = icount + 1
280
281 ALLOCATE (map1, source=perm(1:isep, iperm))
282 ALLOCATE (map2, source=perm(isep + 1:ndims, iperm))
283
284 mynode = mp_comm%mepos
285 CALL mp_environ_pgrid(comm_nd, pdims, myploc)
286
287 IF (1 <= ndims) THEN
288 nblks = SIZE(blk_size_1)
289 ALLOCATE (dist1_1(nblks))
290 ALLOCATE (dist2_1(nblks))
291 CALL dbt_default_distvec(nblks, pdims(1), blk_size_1, dist1_1)
292 CALL dbt_default_distvec(nblks, pdims(1), blk_size_1, dist2_1)
293 END IF
294 IF (2 <= ndims) THEN
295 nblks = SIZE(blk_size_2)
296 ALLOCATE (dist1_2(nblks))
297 ALLOCATE (dist2_2(nblks))
298 CALL dbt_default_distvec(nblks, pdims(2), blk_size_2, dist1_2)
299 CALL dbt_default_distvec(nblks, pdims(2), blk_size_2, dist2_2)
300 END IF
301 IF (3 <= ndims) THEN
302 nblks = SIZE(blk_size_3)
303 ALLOCATE (dist1_3(nblks))
304 ALLOCATE (dist2_3(nblks))
305 CALL dbt_default_distvec(nblks, pdims(3), blk_size_3, dist1_3)
306 CALL dbt_default_distvec(nblks, pdims(3), blk_size_3, dist2_3)
307 END IF
308 IF (4 <= ndims) THEN
309 nblks = SIZE(blk_size_4)
310 ALLOCATE (dist1_4(nblks))
311 ALLOCATE (dist2_4(nblks))
312 CALL dbt_default_distvec(nblks, pdims(4), blk_size_4, dist1_4)
313 CALL dbt_default_distvec(nblks, pdims(4), blk_size_4, dist2_4)
314 END IF
315
316 WRITE (tensor_name, '(A,1X,I3,1X)') "Test", icount
317
318 IF (io_unit > 0) THEN
319 WRITE (io_unit, *)
320 WRITE (io_unit, '(A,A,1X)', advance='no') trim(tensor_name), ':'
321 WRITE (io_unit, '(A1,1X)', advance='no') "("
322 DO imap = 1, SIZE(map1)
323 WRITE (io_unit, '(I1,1X)', advance='no') map1(imap)
324 END DO
325 WRITE (io_unit, '(A1,1X)', advance='no') "|"
326 DO imap = 1, SIZE(map2)
327 WRITE (io_unit, '(I1,1X)', advance='no') map2(imap)
328 END DO
329 WRITE (io_unit, '(A1)') ")"
330
331 WRITE (io_unit, '(T4,A)') "Reference distribution:"
332 IF (1 <= ndims) THEN
333 WRITE (io_unit, '(T7,A,1X)', advance='no') "Dist vec 1:"
334 DO idist = 1, SIZE(dist2_1)
335 WRITE (io_unit, '(I2,1X)', advance='no') dist2_1(idist)
336 END DO
337 WRITE (io_unit, *)
338 END IF
339 IF (2 <= ndims) THEN
340 WRITE (io_unit, '(T7,A,1X)', advance='no') "Dist vec 2:"
341 DO idist = 1, SIZE(dist2_2)
342 WRITE (io_unit, '(I2,1X)', advance='no') dist2_2(idist)
343 END DO
344 WRITE (io_unit, *)
345 END IF
346 IF (3 <= ndims) THEN
347 WRITE (io_unit, '(T7,A,1X)', advance='no') "Dist vec 3:"
348 DO idist = 1, SIZE(dist2_3)
349 WRITE (io_unit, '(I2,1X)', advance='no') dist2_3(idist)
350 END DO
351 WRITE (io_unit, *)
352 END IF
353 IF (4 <= ndims) THEN
354 WRITE (io_unit, '(T7,A,1X)', advance='no') "Dist vec 4:"
355 DO idist = 1, SIZE(dist2_4)
356 WRITE (io_unit, '(I2,1X)', advance='no') dist2_4(idist)
357 END DO
358 WRITE (io_unit, *)
359 END IF
360
361 WRITE (io_unit, '(T4,A)') "Test distribution:"
362 IF (1 <= ndims) THEN
363 WRITE (io_unit, '(T7,A,1X)', advance='no') "Dist vec 1:"
364 DO idist = 1, SIZE(dist2_1)
365 WRITE (io_unit, '(I2,1X)', advance='no') dist1_1(idist)
366 END DO
367 WRITE (io_unit, *)
368 END IF
369 IF (2 <= ndims) THEN
370 WRITE (io_unit, '(T7,A,1X)', advance='no') "Dist vec 2:"
371 DO idist = 1, SIZE(dist2_2)
372 WRITE (io_unit, '(I2,1X)', advance='no') dist1_2(idist)
373 END DO
374 WRITE (io_unit, *)
375 END IF
376 IF (3 <= ndims) THEN
377 WRITE (io_unit, '(T7,A,1X)', advance='no') "Dist vec 3:"
378 DO idist = 1, SIZE(dist2_3)
379 WRITE (io_unit, '(I2,1X)', advance='no') dist1_3(idist)
380 END DO
381 WRITE (io_unit, *)
382 END IF
383 IF (4 <= ndims) THEN
384 WRITE (io_unit, '(T7,A,1X)', advance='no') "Dist vec 4:"
385 DO idist = 1, SIZE(dist2_4)
386 WRITE (io_unit, '(I2,1X)', advance='no') dist1_4(idist)
387 END DO
388 WRITE (io_unit, *)
389 END IF
390 END IF
391
392 IF (ndims == 2) THEN
393 CALL dbt_distribution_new(dist2, comm_nd, dist2_1, dist2_2)
394 CALL dbt_create(tensor2, "Ref", dist2, map1_ref, map2_ref, &
395 blk_size_1, blk_size_2)
396 CALL dbt_setup_test_tensor(tensor2, comm_nd%mp_comm_2d, .true., blk_ind_1, blk_ind_2)
397 END IF
398 IF (ndims == 3) THEN
399 CALL dbt_distribution_new(dist2, comm_nd, dist2_1, dist2_2, dist2_3)
400 CALL dbt_create(tensor2, "Ref", dist2, map1_ref, map2_ref, &
401 blk_size_1, blk_size_2, blk_size_3)
402 CALL dbt_setup_test_tensor(tensor2, comm_nd%mp_comm_2d, .true., blk_ind_1, blk_ind_2, blk_ind_3)
403 END IF
404 IF (ndims == 4) THEN
405 CALL dbt_distribution_new(dist2, comm_nd, dist2_1, dist2_2, dist2_3, dist2_4)
406 CALL dbt_create(tensor2, "Ref", dist2, map1_ref, map2_ref, &
407 blk_size_1, blk_size_2, blk_size_3, blk_size_4)
408 CALL dbt_setup_test_tensor(tensor2, comm_nd%mp_comm_2d, .true., blk_ind_1, blk_ind_2, blk_ind_3, blk_ind_4)
409 END IF
410
411 IF (verbose) CALL dbt_write_blocks(tensor2, io_unit, unit_nr)
412
413 IF (ndims == 2) THEN
414 CALL dbt_distribution_new(dist1, comm_nd, dist1_1, dist1_2)
415 CALL dbt_create(tensor1, tensor_name, dist1, map1, map2, &
416 blk_size_1, blk_size_2)
417 CALL dbt_setup_test_tensor(tensor1, comm_nd%mp_comm_2d, .true., blk_ind_1, blk_ind_2)
418 END IF
419 IF (ndims == 3) THEN
420 CALL dbt_distribution_new(dist1, comm_nd, dist1_1, dist1_2, dist1_3)
421 CALL dbt_create(tensor1, tensor_name, dist1, map1, map2, &
422 blk_size_1, blk_size_2, blk_size_3)
423 CALL dbt_setup_test_tensor(tensor1, comm_nd%mp_comm_2d, .true., blk_ind_1, blk_ind_2, blk_ind_3)
424 END IF
425 IF (ndims == 4) THEN
426 CALL dbt_distribution_new(dist1, comm_nd, dist1_1, dist1_2, dist1_3, dist1_4)
427 CALL dbt_create(tensor1, tensor_name, dist1, map1, map2, &
428 blk_size_1, blk_size_2, blk_size_3, blk_size_4)
429 CALL dbt_setup_test_tensor(tensor1, comm_nd%mp_comm_2d, .true., blk_ind_1, blk_ind_2, blk_ind_3, blk_ind_4)
430 END IF
431
432 IF (verbose) CALL dbt_write_blocks(tensor1, io_unit, unit_nr)
433
434 eql = dbt_equal(tensor1, tensor2)
435
436 IF (.NOT. eql) THEN
437 IF (io_unit > 0) WRITE (io_unit, '(A,1X,A)') trim(tensor_name), 'Test failed!'
438 cpabort('')
439 ELSE
440 IF (io_unit > 0) WRITE (io_unit, '(A,1X,A)') trim(tensor_name), 'Test passed!'
441 END IF
442 DEALLOCATE (map1, map2)
443
444 CALL dbt_destroy(tensor1)
445 CALL dbt_distribution_destroy(dist1)
446
447 CALL dbt_destroy(tensor2)
448 CALL dbt_distribution_destroy(dist2)
449
450 IF (1 <= ndims) THEN
451 DEALLOCATE (dist1_1, dist2_1)
452 END IF
453 IF (2 <= ndims) THEN
454 DEALLOCATE (dist1_2, dist2_2)
455 END IF
456 IF (3 <= ndims) THEN
457 DEALLOCATE (dist1_3, dist2_3)
458 END IF
459 IF (4 <= ndims) THEN
460 DEALLOCATE (dist1_4, dist2_4)
461 END IF
462
463 END DO
464 END DO
465 CALL dbt_pgrid_destroy(comm_nd)
466 END SUBROUTINE
467
468! **************************************************************************************************
469!> \brief Allocate and fill test tensor - entries are enumerated by their index s.t. they only depend
470!> on global properties of the tensor but not on distribution, matrix representation, etc.
471!> \param mp_comm communicator
472!> \param blk_ind_i index along respective dimension of non-zero blocks
473!> \author Patrick Seewald
474! **************************************************************************************************
475 SUBROUTINE dbt_setup_test_tensor(tensor, mp_comm, enumerate, blk_ind_1, blk_ind_2, blk_ind_3, blk_ind_4)
476 TYPE(dbt_type), INTENT(INOUT) :: tensor
477 CLASS(mp_comm_type), INTENT(IN) :: mp_comm
478 LOGICAL, INTENT(IN) :: enumerate
479 INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: blk_ind_1, blk_ind_2, blk_ind_3, blk_ind_4
480 INTEGER :: mynode
481
482 INTEGER :: i, ib, my_nblks_alloc, nblks_alloc, proc, nze
483 INTEGER, ALLOCATABLE, DIMENSION(:) :: my_blk_ind_1, my_blk_ind_2, my_blk_ind_3, my_blk_ind_4
484 INTEGER, DIMENSION(ndims_tensor(tensor)) :: blk_index, blk_offset, blk_size, &
485 tensor_dims
486 INTEGER, DIMENSION(:, :), ALLOCATABLE :: ind_nd
487 REAL(kind=dp), ALLOCATABLE, &
488 DIMENSION(:,:) :: blk_values_2
489 REAL(kind=dp), ALLOCATABLE, &
490 DIMENSION(:,:,:) :: blk_values_3
491 REAL(kind=dp), ALLOCATABLE, &
492 DIMENSION(:,:,:,:) :: blk_values_4
493 TYPE(dbt_iterator_type) :: iterator
494 INTEGER, DIMENSION(4) :: iseed
495 INTEGER, DIMENSION(2) :: blk_index_2d, nblks_2d
496
497 nblks_alloc = SIZE(blk_ind_1)
498 mynode = mp_comm%mepos
499
500 IF (.NOT. enumerate) THEN
501 cpassert(randmat_counter .NE. 0)
502
503 randmat_counter = randmat_counter + 1
504 END IF
505
506 ALLOCATE (ind_nd(nblks_alloc, ndims_tensor(tensor)))
507
508!$OMP PARALLEL DEFAULT(NONE) SHARED(ind_nd,nblks_alloc,tensor,mynode,enumerate,randmat_counter) &
509!$OMP SHARED(blk_ind_1, blk_ind_2, blk_ind_3, blk_ind_4) &
510!$OMP PRIVATE(my_nblks_alloc,ib,proc,i,iterator,blk_offset,blk_size,blk_index) &
511!$OMP PRIVATE(blk_index_2d,nblks_2d,iseed,nze,tensor_dims) &
512!$OMP PRIVATE(blk_values_2, blk_values_3, blk_values_4) &
513!$OMP PRIVATE(my_blk_ind_1, my_blk_ind_2, my_blk_ind_3, my_blk_ind_4)
514
515 my_nblks_alloc = 0
516!$OMP DO
517 DO ib = 1, nblks_alloc
518 IF (ndims_tensor(tensor) == 2) THEN
519 ind_nd(ib, :) = [blk_ind_1(ib), blk_ind_2(ib)]
520 END IF
521 IF (ndims_tensor(tensor) == 3) THEN
522 ind_nd(ib, :) = [blk_ind_1(ib), blk_ind_2(ib), blk_ind_3(ib)]
523 END IF
524 IF (ndims_tensor(tensor) == 4) THEN
525 ind_nd(ib, :) = [blk_ind_1(ib), blk_ind_2(ib), blk_ind_3(ib), blk_ind_4(ib)]
526 END IF
527 CALL dbt_get_stored_coordinates(tensor, ind_nd(ib, :), proc)
528 IF (proc == mynode) THEN
529 my_nblks_alloc = my_nblks_alloc + 1
530 END IF
531 END DO
532!$OMP END DO
533
534 IF (ndims_tensor(tensor) >= 1) THEN
535 ALLOCATE (my_blk_ind_1(my_nblks_alloc))
536 END IF
537 IF (ndims_tensor(tensor) >= 2) THEN
538 ALLOCATE (my_blk_ind_2(my_nblks_alloc))
539 END IF
540 IF (ndims_tensor(tensor) >= 3) THEN
541 ALLOCATE (my_blk_ind_3(my_nblks_alloc))
542 END IF
543 IF (ndims_tensor(tensor) >= 4) THEN
544 ALLOCATE (my_blk_ind_4(my_nblks_alloc))
545 END IF
546
547 i = 0
548!$OMP DO
549 DO ib = 1, nblks_alloc
550 CALL dbt_get_stored_coordinates(tensor, ind_nd(ib, :), proc)
551 IF (proc == mynode) THEN
552 i = i + 1
553 IF (ndims_tensor(tensor) >= 1) THEN
554 my_blk_ind_1(i) = blk_ind_1(ib)
555 END IF
556 IF (ndims_tensor(tensor) >= 2) THEN
557 my_blk_ind_2(i) = blk_ind_2(ib)
558 END IF
559 IF (ndims_tensor(tensor) >= 3) THEN
560 my_blk_ind_3(i) = blk_ind_3(ib)
561 END IF
562 IF (ndims_tensor(tensor) >= 4) THEN
563 my_blk_ind_4(i) = blk_ind_4(ib)
564 END IF
565 END IF
566 END DO
567!$OMP END DO
568
569 IF (ndims_tensor(tensor) == 2) THEN
570 CALL dbt_reserve_blocks(tensor, my_blk_ind_1, my_blk_ind_2)
571 END IF
572 IF (ndims_tensor(tensor) == 3) THEN
573 CALL dbt_reserve_blocks(tensor, my_blk_ind_1, my_blk_ind_2, my_blk_ind_3)
574 END IF
575 IF (ndims_tensor(tensor) == 4) THEN
576 CALL dbt_reserve_blocks(tensor, my_blk_ind_1, my_blk_ind_2, my_blk_ind_3, my_blk_ind_4)
577 END IF
578
579 CALL dbt_iterator_start(iterator, tensor)
580 DO WHILE (dbt_iterator_blocks_left(iterator))
581 CALL dbt_iterator_next_block(iterator, blk_index, blk_size=blk_size, blk_offset=blk_offset)
582
583 IF (.NOT. enumerate) THEN
584 blk_index_2d = int(get_2d_indices_tensor(tensor%nd_index_blk, blk_index))
585 CALL dbt_get_mapping_info(tensor%nd_index_blk, dims_2d=nblks_2d)
586 iseed = generate_larnv_seed(blk_index_2d(1), nblks_2d(1), blk_index_2d(2), nblks_2d(2), randmat_counter)
587 nze = product(blk_size)
588 END IF
589
590 IF (ndims_tensor(tensor) == 2) THEN
591 CALL allocate_any(blk_values_2, shape_spec=blk_size)
592 CALL dims_tensor(tensor, tensor_dims)
593 IF (enumerate) THEN
594 CALL enumerate_block_elements(blk_size, blk_offset, tensor_dims, blk_2=blk_values_2)
595 ELSE
596 CALL dlarnv(1, iseed, nze, blk_values_2)
597 END IF
598 CALL dbt_put_block(tensor, blk_index, blk_size, blk_values_2)
599 DEALLOCATE (blk_values_2)
600 END IF
601 IF (ndims_tensor(tensor) == 3) THEN
602 CALL allocate_any(blk_values_3, shape_spec=blk_size)
603 CALL dims_tensor(tensor, tensor_dims)
604 IF (enumerate) THEN
605 CALL enumerate_block_elements(blk_size, blk_offset, tensor_dims, blk_3=blk_values_3)
606 ELSE
607 CALL dlarnv(1, iseed, nze, blk_values_3)
608 END IF
609 CALL dbt_put_block(tensor, blk_index, blk_size, blk_values_3)
610 DEALLOCATE (blk_values_3)
611 END IF
612 IF (ndims_tensor(tensor) == 4) THEN
613 CALL allocate_any(blk_values_4, shape_spec=blk_size)
614 CALL dims_tensor(tensor, tensor_dims)
615 IF (enumerate) THEN
616 CALL enumerate_block_elements(blk_size, blk_offset, tensor_dims, blk_4=blk_values_4)
617 ELSE
618 CALL dlarnv(1, iseed, nze, blk_values_4)
619 END IF
620 CALL dbt_put_block(tensor, blk_index, blk_size, blk_values_4)
621 DEALLOCATE (blk_values_4)
622 END IF
623 END DO
624 CALL dbt_iterator_stop(iterator)
625!$OMP END PARALLEL
626
627 END SUBROUTINE
628
629! **************************************************************************************************
630!> \brief Enumerate tensor entries in block
631!> \param blk_2 block values for 2 dimensions
632!> \param blk_3 block values for 3 dimensions
633!> \param blk_size size of block
634!> \param blk_offset block offset (indices of first element)
635!> \param tensor_size global tensor sizes
636!> \author Patrick Seewald
637! **************************************************************************************************
638 SUBROUTINE enumerate_block_elements(blk_size, blk_offset, tensor_size, blk_2, blk_3, blk_4)
639 INTEGER, DIMENSION(:), INTENT(IN) :: blk_size, blk_offset, tensor_size
640 REAL(kind=dp), DIMENSION(:,:), &
641 OPTIONAL, INTENT(OUT) :: blk_2
642 REAL(kind=dp), DIMENSION(:,:,:), &
643 OPTIONAL, INTENT(OUT) :: blk_3
644 REAL(kind=dp), DIMENSION(:,:,:,:), &
645 OPTIONAL, INTENT(OUT) :: blk_4
646 INTEGER :: ndim
647 INTEGER, DIMENSION(SIZE(blk_size)) :: arr_ind, tens_ind
648 INTEGER :: i_1, i_2, i_3, i_4
649
650 ndim = SIZE(tensor_size)
651
652 IF (ndim == 2) THEN
653 DO i_2 = 1, blk_size(2)
654 DO i_1 = 1, blk_size(1)
655 arr_ind(:) = [i_1, i_2]
656 tens_ind(:) = arr_ind(:) + blk_offset(:) - 1
657 blk_2(arr_ind(1), arr_ind(2)) = combine_tensor_index(tens_ind, tensor_size)
658 END DO
659 END DO
660 END IF
661 IF (ndim == 3) THEN
662 DO i_3 = 1, blk_size(3)
663 DO i_2 = 1, blk_size(2)
664 DO i_1 = 1, blk_size(1)
665 arr_ind(:) = [i_1, i_2, i_3]
666 tens_ind(:) = arr_ind(:) + blk_offset(:) - 1
667 blk_3(arr_ind(1), arr_ind(2), arr_ind(3)) = combine_tensor_index(tens_ind, tensor_size)
668 END DO
669 END DO
670 END DO
671 END IF
672 IF (ndim == 4) THEN
673 DO i_4 = 1, blk_size(4)
674 DO i_3 = 1, blk_size(3)
675 DO i_2 = 1, blk_size(2)
676 DO i_1 = 1, blk_size(1)
677 arr_ind(:) = [i_1, i_2, i_3, i_4]
678 tens_ind(:) = arr_ind(:) + blk_offset(:) - 1
679 blk_4(arr_ind(1), arr_ind(2), arr_ind(3), arr_ind(4)) = combine_tensor_index(tens_ind, tensor_size)
680 END DO
681 END DO
682 END DO
683 END DO
684 END IF
685
686 END SUBROUTINE
687
688! **************************************************************************************************
689!> \brief Transform a distributed sparse tensor to a replicated dense array. This is only useful for
690!> testing tensor contraction by matrix multiplication of dense arrays.
691!> \author Patrick Seewald
692! **************************************************************************************************
693 SUBROUTINE dist_sparse_tensor_to_repl_dense_2d_array(tensor, array)
694
695 TYPE(dbt_type), INTENT(INOUT) :: tensor
696 REAL(dp), ALLOCATABLE, DIMENSION(:,:), &
697 INTENT(OUT) :: array
698 REAL(dp), ALLOCATABLE, DIMENSION(:,:) :: block
699 INTEGER, DIMENSION(ndims_tensor(tensor)) :: dims_nd, ind_nd, blk_size, blk_offset
700 TYPE(dbt_iterator_type) :: iterator
701 INTEGER :: idim
702 INTEGER, DIMENSION(ndims_tensor(tensor)) :: blk_start, blk_end
703 LOGICAL :: found
704
705 cpassert(ndims_tensor(tensor) .EQ. 2)
706 CALL dbt_get_info(tensor, nfull_total=dims_nd)
707 CALL allocate_any(array, shape_spec=dims_nd)
708 array(:,:) = 0.0_dp
709
710!$OMP PARALLEL DEFAULT(NONE) SHARED(tensor,array) &
711!$OMP PRIVATE(iterator,ind_nd,blk_size,blk_offset,block,found,idim,blk_start,blk_end)
712 CALL dbt_iterator_start(iterator, tensor)
713 DO WHILE (dbt_iterator_blocks_left(iterator))
714 CALL dbt_iterator_next_block(iterator, ind_nd, blk_size=blk_size, blk_offset=blk_offset)
715 CALL dbt_get_block(tensor, ind_nd, block, found)
716 cpassert(found)
717
718 DO idim = 1, ndims_tensor(tensor)
719 blk_start(idim) = blk_offset(idim)
720 blk_end(idim) = blk_offset(idim) + blk_size(idim) - 1
721 END DO
722 array(blk_start(1):blk_end(1), blk_start(2):blk_end(2)) = &
723 block(:,:)
724
725 DEALLOCATE (block)
726 END DO
727 CALL dbt_iterator_stop(iterator)
728!$OMP END PARALLEL
729 CALL tensor%pgrid%mp_comm_2d%sum(array)
730
731 END SUBROUTINE
732! **************************************************************************************************
733!> \brief Transform a distributed sparse tensor to a replicated dense array. This is only useful for
734!> testing tensor contraction by matrix multiplication of dense arrays.
735!> \author Patrick Seewald
736! **************************************************************************************************
737 SUBROUTINE dist_sparse_tensor_to_repl_dense_3d_array(tensor, array)
738
739 TYPE(dbt_type), INTENT(INOUT) :: tensor
740 REAL(dp), ALLOCATABLE, DIMENSION(:,:,:), &
741 INTENT(OUT) :: array
742 REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) :: block
743 INTEGER, DIMENSION(ndims_tensor(tensor)) :: dims_nd, ind_nd, blk_size, blk_offset
744 TYPE(dbt_iterator_type) :: iterator
745 INTEGER :: idim
746 INTEGER, DIMENSION(ndims_tensor(tensor)) :: blk_start, blk_end
747 LOGICAL :: found
748
749 cpassert(ndims_tensor(tensor) .EQ. 3)
750 CALL dbt_get_info(tensor, nfull_total=dims_nd)
751 CALL allocate_any(array, shape_spec=dims_nd)
752 array(:,:,:) = 0.0_dp
753
754!$OMP PARALLEL DEFAULT(NONE) SHARED(tensor,array) &
755!$OMP PRIVATE(iterator,ind_nd,blk_size,blk_offset,block,found,idim,blk_start,blk_end)
756 CALL dbt_iterator_start(iterator, tensor)
757 DO WHILE (dbt_iterator_blocks_left(iterator))
758 CALL dbt_iterator_next_block(iterator, ind_nd, blk_size=blk_size, blk_offset=blk_offset)
759 CALL dbt_get_block(tensor, ind_nd, block, found)
760 cpassert(found)
761
762 DO idim = 1, ndims_tensor(tensor)
763 blk_start(idim) = blk_offset(idim)
764 blk_end(idim) = blk_offset(idim) + blk_size(idim) - 1
765 END DO
766 array(blk_start(1):blk_end(1), blk_start(2):blk_end(2), blk_start(3):blk_end(3)) = &
767 block(:,:,:)
768
769 DEALLOCATE (block)
770 END DO
771 CALL dbt_iterator_stop(iterator)
772!$OMP END PARALLEL
773 CALL tensor%pgrid%mp_comm_2d%sum(array)
774
775 END SUBROUTINE
776! **************************************************************************************************
777!> \brief Transform a distributed sparse tensor to a replicated dense array. This is only useful for
778!> testing tensor contraction by matrix multiplication of dense arrays.
779!> \author Patrick Seewald
780! **************************************************************************************************
781 SUBROUTINE dist_sparse_tensor_to_repl_dense_4d_array(tensor, array)
782
783 TYPE(dbt_type), INTENT(INOUT) :: tensor
784 REAL(dp), ALLOCATABLE, DIMENSION(:,:,:,:), &
785 INTENT(OUT) :: array
786 REAL(dp), ALLOCATABLE, DIMENSION(:,:,:,:) :: block
787 INTEGER, DIMENSION(ndims_tensor(tensor)) :: dims_nd, ind_nd, blk_size, blk_offset
788 TYPE(dbt_iterator_type) :: iterator
789 INTEGER :: idim
790 INTEGER, DIMENSION(ndims_tensor(tensor)) :: blk_start, blk_end
791 LOGICAL :: found
792
793 cpassert(ndims_tensor(tensor) .EQ. 4)
794 CALL dbt_get_info(tensor, nfull_total=dims_nd)
795 CALL allocate_any(array, shape_spec=dims_nd)
796 array(:,:,:,:) = 0.0_dp
797
798!$OMP PARALLEL DEFAULT(NONE) SHARED(tensor,array) &
799!$OMP PRIVATE(iterator,ind_nd,blk_size,blk_offset,block,found,idim,blk_start,blk_end)
800 CALL dbt_iterator_start(iterator, tensor)
801 DO WHILE (dbt_iterator_blocks_left(iterator))
802 CALL dbt_iterator_next_block(iterator, ind_nd, blk_size=blk_size, blk_offset=blk_offset)
803 CALL dbt_get_block(tensor, ind_nd, block, found)
804 cpassert(found)
805
806 DO idim = 1, ndims_tensor(tensor)
807 blk_start(idim) = blk_offset(idim)
808 blk_end(idim) = blk_offset(idim) + blk_size(idim) - 1
809 END DO
810 array(blk_start(1):blk_end(1), blk_start(2):blk_end(2), blk_start(3):blk_end(3), blk_start(4):blk_end(4)) = &
811 block(:,:,:,:)
812
813 DEALLOCATE (block)
814 END DO
815 CALL dbt_iterator_stop(iterator)
816!$OMP END PARALLEL
817 CALL tensor%pgrid%mp_comm_2d%sum(array)
818
819 END SUBROUTINE
820
821! **************************************************************************************************
822!> \brief test tensor contraction
823!> \note for testing/debugging, simply replace a call to dbt_contract with a call to this routine
824!> \author Patrick Seewald
825! **************************************************************************************************
826 SUBROUTINE dbt_contract_test(alpha, tensor_1, tensor_2, beta, tensor_3, &
827 contract_1, notcontract_1, &
828 contract_2, notcontract_2, &
829 map_1, map_2, &
830 unit_nr, &
831 bounds_1, bounds_2, bounds_3, &
832 log_verbose, write_int)
833
834 REAL(dp), INTENT(IN) :: alpha
835 TYPE(dbt_type), INTENT(INOUT) :: tensor_1, tensor_2, tensor_3
836 REAL(dp), INTENT(IN) :: beta
837 INTEGER, DIMENSION(:), INTENT(IN) :: contract_1, contract_2, &
838 notcontract_1, notcontract_2, &
839 map_1, map_2
840 INTEGER, INTENT(IN) :: unit_nr
841 INTEGER, DIMENSION(2, SIZE(contract_1)), &
842 OPTIONAL :: bounds_1
843 INTEGER, DIMENSION(2, SIZE(notcontract_1)), &
844 OPTIONAL :: bounds_2
845 INTEGER, DIMENSION(2, SIZE(notcontract_2)), &
846 OPTIONAL :: bounds_3
847 LOGICAL, INTENT(IN), OPTIONAL :: log_verbose
848 LOGICAL, INTENT(IN), OPTIONAL :: write_int
849 INTEGER :: io_unit, mynode
850 TYPE(mp_comm_type) :: mp_comm
851 INTEGER, DIMENSION(:), ALLOCATABLE :: size_1, size_2, size_3, &
852 order_t1, order_t2, order_t3
853 INTEGER, DIMENSION(2, ndims_tensor(tensor_1)) :: bounds_t1
854 INTEGER, DIMENSION(2, ndims_tensor(tensor_2)) :: bounds_t2
855
856 REAL(kind=dp), ALLOCATABLE, &
857 DIMENSION(:,:) :: array_1_2d, &
858 array_2_2d, &
859 array_3_2d, &
860 array_1_2d_full, &
861 array_2_2d_full, &
862 array_3_0_2d, &
863 array_1_rs2d, &
864 array_2_rs2d, &
865 array_3_rs2d, &
866 array_3_0_rs2d
867 REAL(kind=dp), ALLOCATABLE, &
868 DIMENSION(:,:,:) :: array_1_3d, &
869 array_2_3d, &
870 array_3_3d, &
871 array_1_3d_full, &
872 array_2_3d_full, &
873 array_3_0_3d, &
874 array_1_rs3d, &
875 array_2_rs3d, &
876 array_3_rs3d, &
877 array_3_0_rs3d
878 REAL(kind=dp), ALLOCATABLE, &
879 DIMENSION(:,:,:,:) :: array_1_4d, &
880 array_2_4d, &
881 array_3_4d, &
882 array_1_4d_full, &
883 array_2_4d_full, &
884 array_3_0_4d, &
885 array_1_rs4d, &
886 array_2_rs4d, &
887 array_3_rs4d, &
888 array_3_0_rs4d
889 REAL(kind=dp), ALLOCATABLE, &
890 DIMENSION(:, :) :: array_1_mm, &
891 array_2_mm, &
892 array_3_mm, &
893 array_3_test_mm
894 LOGICAL :: eql, notzero
895 LOGICAL, PARAMETER :: debug = .false.
896 REAL(kind=dp) :: cs_1, cs_2, cs_3, eql_diff
897 LOGICAL :: do_crop_1, do_crop_2
898
899 mp_comm = tensor_1%pgrid%mp_comm_2d
900 mynode = mp_comm%mepos
901 io_unit = -1
902 IF (mynode .EQ. 0) io_unit = unit_nr
903
904 cs_1 = dbt_checksum(tensor_1)
905 cs_2 = dbt_checksum(tensor_2)
906 cs_3 = dbt_checksum(tensor_3)
907
908 IF (io_unit > 0) THEN
909 WRITE (io_unit, *)
910 WRITE (io_unit, '(A)') repeat("-", 80)
911 WRITE (io_unit, '(A,1X,A,1X,A,1X,A,1X,A,1X,A)') "Testing tensor contraction", &
912 trim(tensor_1%name), "x", trim(tensor_2%name), "=", trim(tensor_3%name)
913 WRITE (io_unit, '(A)') repeat("-", 80)
914 END IF
915
916 IF (debug) THEN
917 IF (io_unit > 0) THEN
918 WRITE (io_unit, "(A, E9.2)") "checksum ", trim(tensor_1%name), cs_1
919 WRITE (io_unit, "(A, E9.2)") "checksum ", trim(tensor_2%name), cs_2
920 WRITE (io_unit, "(A, E9.2)") "checksum ", trim(tensor_3%name), cs_3
921 END IF
922 END IF
923
924 IF (debug) THEN
925 CALL dbt_write_block_indices(tensor_1, io_unit, unit_nr)
926 CALL dbt_write_blocks(tensor_1, io_unit, unit_nr, write_int)
927 END IF
928
929 SELECT CASE (ndims_tensor(tensor_3))
930 CASE (2)
931 CALL dist_sparse_tensor_to_repl_dense_array(tensor_3, array_3_0_2d)
932 CASE (3)
933 CALL dist_sparse_tensor_to_repl_dense_array(tensor_3, array_3_0_3d)
934 CASE (4)
935 CALL dist_sparse_tensor_to_repl_dense_array(tensor_3, array_3_0_4d)
936 END SELECT
937
938 CALL dbt_contract(alpha, tensor_1, tensor_2, beta, tensor_3, &
939 contract_1, notcontract_1, &
940 contract_2, notcontract_2, &
941 map_1, map_2, &
942 bounds_1=bounds_1, bounds_2=bounds_2, bounds_3=bounds_3, &
943 filter_eps=1.0e-12_dp, &
944 unit_nr=io_unit, log_verbose=log_verbose)
945
946 cs_3 = dbt_checksum(tensor_3)
947
948 IF (debug) THEN
949 IF (io_unit > 0) THEN
950 WRITE (io_unit, "(A, E9.2)") "checksum ", trim(tensor_3%name), cs_3
951 END IF
952 END IF
953
954 do_crop_1 = .false.; do_crop_2 = .false.!; do_crop_3 = .FALSE.
955
956 ! crop tensor as first step
957 bounds_t1(1, :) = 1
958 CALL dbt_get_info(tensor_1, nfull_total=bounds_t1(2, :))
959
960 bounds_t2(1, :) = 1
961 CALL dbt_get_info(tensor_2, nfull_total=bounds_t2(2, :))
962
963 IF (PRESENT(bounds_1)) THEN
964 bounds_t1(:, contract_1) = bounds_1
965 do_crop_1 = .true.
966 bounds_t2(:, contract_2) = bounds_1
967 do_crop_2 = .true.
968 END IF
969
970 IF (PRESENT(bounds_2)) THEN
971 bounds_t1(:, notcontract_1) = bounds_2
972 do_crop_1 = .true.
973 END IF
974
975 IF (PRESENT(bounds_3)) THEN
976 bounds_t2(:, notcontract_2) = bounds_3
977 do_crop_2 = .true.
978 END IF
979
980 ! Convert tensors to simple multidimensional arrays
981 SELECT CASE (ndims_tensor(tensor_1))
982 CASE (2)
983 CALL dist_sparse_tensor_to_repl_dense_array(tensor_1, array_1_2d_full)
984 CALL allocate_any(array_1_2d, shape_spec=shape(array_1_2d_full))
985 array_1_2d = 0.0_dp
986 array_1_2d(bounds_t1(1, 1):bounds_t1(2, 1), bounds_t1(1, 2):bounds_t1(2, 2)) = &
987 array_1_2d_full(bounds_t1(1, 1):bounds_t1(2, 1), bounds_t1(1, 2):bounds_t1(2, 2))
988
989 CASE (3)
990 CALL dist_sparse_tensor_to_repl_dense_array(tensor_1, array_1_3d_full)
991 CALL allocate_any(array_1_3d, shape_spec=shape(array_1_3d_full))
992 array_1_3d = 0.0_dp
993 array_1_3d(bounds_t1(1, 1):bounds_t1(2, 1), bounds_t1(1, 2):bounds_t1(2, 2), bounds_t1(1, 3):bounds_t1(2, 3)) = &
994 array_1_3d_full(bounds_t1(1, 1):bounds_t1(2, 1), bounds_t1(1, 2):bounds_t1(2, 2), bounds_t1(1, 3):bounds_t1(2, 3))
995
996 CASE (4)
997 CALL dist_sparse_tensor_to_repl_dense_array(tensor_1, array_1_4d_full)
998 CALL allocate_any(array_1_4d, shape_spec=shape(array_1_4d_full))
999 array_1_4d = 0.0_dp
1000 array_1_4d(bounds_t1(1, 1):bounds_t1(2, 1), bounds_t1(1, 2):bounds_t1(2, 2), bounds_t1(1, 3):bounds_t1(2, 3),&
1001 & bounds_t1(1, 4):bounds_t1(2, 4)) = &
1002 array_1_4d_full(bounds_t1(1, 1):bounds_t1(2, 1), bounds_t1(1, 2):bounds_t1(2, 2), bounds_t1(1, 3):bounds_t1(2, 3),&
1003 & bounds_t1(1, 4):bounds_t1(2, 4))
1004
1005 END SELECT
1006 SELECT CASE (ndims_tensor(tensor_2))
1007 CASE (2)
1008 CALL dist_sparse_tensor_to_repl_dense_array(tensor_2, array_2_2d_full)
1009 CALL allocate_any(array_2_2d, shape_spec=shape(array_2_2d_full))
1010 array_2_2d = 0.0_dp
1011 array_2_2d(bounds_t2(1, 1):bounds_t2(2, 1), bounds_t2(1, 2):bounds_t2(2, 2)) = &
1012 array_2_2d_full(bounds_t2(1, 1):bounds_t2(2, 1), bounds_t2(1, 2):bounds_t2(2, 2))
1013
1014 CASE (3)
1015 CALL dist_sparse_tensor_to_repl_dense_array(tensor_2, array_2_3d_full)
1016 CALL allocate_any(array_2_3d, shape_spec=shape(array_2_3d_full))
1017 array_2_3d = 0.0_dp
1018 array_2_3d(bounds_t2(1, 1):bounds_t2(2, 1), bounds_t2(1, 2):bounds_t2(2, 2), bounds_t2(1, 3):bounds_t2(2, 3)) = &
1019 array_2_3d_full(bounds_t2(1, 1):bounds_t2(2, 1), bounds_t2(1, 2):bounds_t2(2, 2), bounds_t2(1, 3):bounds_t2(2, 3))
1020
1021 CASE (4)
1022 CALL dist_sparse_tensor_to_repl_dense_array(tensor_2, array_2_4d_full)
1023 CALL allocate_any(array_2_4d, shape_spec=shape(array_2_4d_full))
1024 array_2_4d = 0.0_dp
1025 array_2_4d(bounds_t2(1, 1):bounds_t2(2, 1), bounds_t2(1, 2):bounds_t2(2, 2), bounds_t2(1, 3):bounds_t2(2, 3),&
1026 & bounds_t2(1, 4):bounds_t2(2, 4)) = &
1027 array_2_4d_full(bounds_t2(1, 1):bounds_t2(2, 1), bounds_t2(1, 2):bounds_t2(2, 2), bounds_t2(1, 3):bounds_t2(2, 3),&
1028 & bounds_t2(1, 4):bounds_t2(2, 4))
1029
1030 END SELECT
1031 SELECT CASE (ndims_tensor(tensor_3))
1032 CASE (2)
1033 CALL dist_sparse_tensor_to_repl_dense_array(tensor_3, array_3_2d)
1034
1035 CASE (3)
1036 CALL dist_sparse_tensor_to_repl_dense_array(tensor_3, array_3_3d)
1037
1038 CASE (4)
1039 CALL dist_sparse_tensor_to_repl_dense_array(tensor_3, array_3_4d)
1040
1041 END SELECT
1042
1043 ! Get array sizes
1044
1045 SELECT CASE (ndims_tensor(tensor_1))
1046 CASE (2)
1047 ALLOCATE (size_1, source=shape(array_1_2d))
1048
1049 CASE (3)
1050 ALLOCATE (size_1, source=shape(array_1_3d))
1051
1052 CASE (4)
1053 ALLOCATE (size_1, source=shape(array_1_4d))
1054
1055 END SELECT
1056 SELECT CASE (ndims_tensor(tensor_2))
1057 CASE (2)
1058 ALLOCATE (size_2, source=shape(array_2_2d))
1059
1060 CASE (3)
1061 ALLOCATE (size_2, source=shape(array_2_3d))
1062
1063 CASE (4)
1064 ALLOCATE (size_2, source=shape(array_2_4d))
1065
1066 END SELECT
1067 SELECT CASE (ndims_tensor(tensor_3))
1068 CASE (2)
1069 ALLOCATE (size_3, source=shape(array_3_2d))
1070
1071 CASE (3)
1072 ALLOCATE (size_3, source=shape(array_3_3d))
1073
1074 CASE (4)
1075 ALLOCATE (size_3, source=shape(array_3_4d))
1076
1077 END SELECT
1078
1079 ALLOCATE (order_t1(ndims_tensor(tensor_1)))
1080 ALLOCATE (order_t2(ndims_tensor(tensor_2)))
1081 ALLOCATE (order_t3(ndims_tensor(tensor_3)))
1082
1083 associate(map_t1_1 => notcontract_1, map_t1_2 => contract_1, &
1084 map_t2_1 => notcontract_2, map_t2_2 => contract_2, &
1085 map_t3_1 => map_1, map_t3_2 => map_2)
1086
1087 order_t1(:) = dbt_inverse_order([map_t1_1, map_t1_2])
1088
1089 SELECT CASE (ndims_tensor(tensor_1))
1090 CASE (2)
1091 CALL allocate_any(array_1_rs2d, source=array_1_2d, order=order_t1)
1092 CALL allocate_any(array_1_mm, sizes_2d(size_1, map_t1_1, map_t1_2))
1093 array_1_mm(:, :) = reshape(array_1_rs2d, shape(array_1_mm))
1094 CASE (3)
1095 CALL allocate_any(array_1_rs3d, source=array_1_3d, order=order_t1)
1096 CALL allocate_any(array_1_mm, sizes_2d(size_1, map_t1_1, map_t1_2))
1097 array_1_mm(:, :) = reshape(array_1_rs3d, shape(array_1_mm))
1098 CASE (4)
1099 CALL allocate_any(array_1_rs4d, source=array_1_4d, order=order_t1)
1100 CALL allocate_any(array_1_mm, sizes_2d(size_1, map_t1_1, map_t1_2))
1101 array_1_mm(:, :) = reshape(array_1_rs4d, shape(array_1_mm))
1102 END SELECT
1103 order_t2(:) = dbt_inverse_order([map_t2_1, map_t2_2])
1104
1105 SELECT CASE (ndims_tensor(tensor_2))
1106 CASE (2)
1107 CALL allocate_any(array_2_rs2d, source=array_2_2d, order=order_t2)
1108 CALL allocate_any(array_2_mm, sizes_2d(size_2, map_t2_1, map_t2_2))
1109 array_2_mm(:, :) = reshape(array_2_rs2d, shape(array_2_mm))
1110 CASE (3)
1111 CALL allocate_any(array_2_rs3d, source=array_2_3d, order=order_t2)
1112 CALL allocate_any(array_2_mm, sizes_2d(size_2, map_t2_1, map_t2_2))
1113 array_2_mm(:, :) = reshape(array_2_rs3d, shape(array_2_mm))
1114 CASE (4)
1115 CALL allocate_any(array_2_rs4d, source=array_2_4d, order=order_t2)
1116 CALL allocate_any(array_2_mm, sizes_2d(size_2, map_t2_1, map_t2_2))
1117 array_2_mm(:, :) = reshape(array_2_rs4d, shape(array_2_mm))
1118 END SELECT
1119 order_t3(:) = dbt_inverse_order([map_t3_1, map_t3_2])
1120
1121 SELECT CASE (ndims_tensor(tensor_3))
1122 CASE (2)
1123 CALL allocate_any(array_3_rs2d, source=array_3_2d, order=order_t3)
1124 CALL allocate_any(array_3_mm, sizes_2d(size_3, map_t3_1, map_t3_2))
1125 array_3_mm(:, :) = reshape(array_3_rs2d, shape(array_3_mm))
1126 CASE (3)
1127 CALL allocate_any(array_3_rs3d, source=array_3_3d, order=order_t3)
1128 CALL allocate_any(array_3_mm, sizes_2d(size_3, map_t3_1, map_t3_2))
1129 array_3_mm(:, :) = reshape(array_3_rs3d, shape(array_3_mm))
1130 CASE (4)
1131 CALL allocate_any(array_3_rs4d, source=array_3_4d, order=order_t3)
1132 CALL allocate_any(array_3_mm, sizes_2d(size_3, map_t3_1, map_t3_2))
1133 array_3_mm(:, :) = reshape(array_3_rs4d, shape(array_3_mm))
1134 END SELECT
1135
1136 SELECT CASE (ndims_tensor(tensor_3))
1137 CASE (2)
1138 CALL allocate_any(array_3_0_rs2d, source=array_3_0_2d, order=order_t3)
1139 CALL allocate_any(array_3_test_mm, sizes_2d(size_3, map_t3_1, map_t3_2))
1140 array_3_test_mm(:, :) = reshape(array_3_0_rs2d, shape(array_3_mm))
1141 CASE (3)
1142 CALL allocate_any(array_3_0_rs3d, source=array_3_0_3d, order=order_t3)
1143 CALL allocate_any(array_3_test_mm, sizes_2d(size_3, map_t3_1, map_t3_2))
1144 array_3_test_mm(:, :) = reshape(array_3_0_rs3d, shape(array_3_mm))
1145 CASE (4)
1146 CALL allocate_any(array_3_0_rs4d, source=array_3_0_4d, order=order_t3)
1147 CALL allocate_any(array_3_test_mm, sizes_2d(size_3, map_t3_1, map_t3_2))
1148 array_3_test_mm(:, :) = reshape(array_3_0_rs4d, shape(array_3_mm))
1149 END SELECT
1150
1151 array_3_test_mm(:, :) = beta*array_3_test_mm(:, :) + alpha*matmul(array_1_mm, transpose(array_2_mm))
1152
1153 END associate
1154
1155 eql_diff = maxval(abs(array_3_test_mm(:, :) - array_3_mm(:, :)))
1156 notzero = maxval(abs(array_3_test_mm(:, :))) .GT. 1.0e-12_dp
1157
1158 eql = eql_diff .LT. 1.0e-11_dp
1159
1160 IF (.NOT. eql .OR. .NOT. notzero) THEN
1161 IF (io_unit > 0) WRITE (io_unit, *) 'Test failed!', eql_diff
1162 cpabort('')
1163 ELSE
1164 IF (io_unit > 0) WRITE (io_unit, *) 'Test passed!', eql_diff
1165 END IF
1166
1167 END SUBROUTINE
1168
1169! **************************************************************************************************
1170!> \brief mapped sizes in 2d
1171!> \author Patrick Seewald
1172! **************************************************************************************************
1173 FUNCTION sizes_2d(nd_sizes, map1, map2)
1174 INTEGER, DIMENSION(:), INTENT(IN) :: nd_sizes, map1, map2
1175 INTEGER, DIMENSION(2) :: sizes_2d
1176 sizes_2d(1) = product(nd_sizes(map1))
1177 sizes_2d(2) = product(nd_sizes(map2))
1178 END FUNCTION
1179
1180! **************************************************************************************************
1181!> \brief checksum of a tensor consistent with block_checksum
1182!> \author Patrick Seewald
1183! **************************************************************************************************
1184 FUNCTION dbt_checksum(tensor)
1185 TYPE(dbt_type), INTENT(IN) :: tensor
1186 REAL(kind=dp) :: dbt_checksum
1187 dbt_checksum = dbt_tas_checksum(tensor%matrix_rep)
1188 END FUNCTION
1189
1190! **************************************************************************************************
1191!> \brief Reset the seed used for generating random matrices to default value
1192!> \author Patrick Seewald
1193! **************************************************************************************************
1194 SUBROUTINE dbt_reset_randmat_seed()
1195 randmat_counter = rand_seed_init
1196 END SUBROUTINE
1197
1198 END MODULE
tensor
struct tensor_ tensor
dbt_allocate_wrap::allocate_any
Definition dbt_allocate_wrap.F:26
dbt_block::dbt_get_block
Definition dbt_block.F:85
dbt_block::dbt_put_block
Definition dbt_block.F:78
dbt_block::dbt_reserve_blocks
Definition dbt_block.F:95
dbt_types::dbt_create
Definition dbt_types.F:166
dbm_tests
Definition dbm_tests.F:8
dbm_tests::generate_larnv_seed
integer function, dimension(4), public generate_larnv_seed(irow, nrow, icol, ncol, ival)
Generate a seed respecting the lapack constraints,.
Definition dbm_tests.F:435
dbt_allocate_wrap
Wrapper for allocating, copying and reshaping arrays.
Definition dbt_allocate_wrap.F:15
dbt_block
Methods to operate on n-dimensional tensor blocks.
Definition dbt_block.F:12
dbt_block::dbt_iterator_blocks_left
logical function, public dbt_iterator_blocks_left(iterator)
Generalization of block_iterator_blocks_left for tensors.
Definition dbt_block.F:197
dbt_block::dbt_iterator_stop
subroutine, public dbt_iterator_stop(iterator)
Generalization of block_iterator_stop for tensors.
Definition dbt_block.F:134
dbt_block::dbt_iterator_start
subroutine, public dbt_iterator_start(iterator, tensor)
Generalization of block_iterator_start for tensors.
Definition dbt_block.F:121
dbt_block::dbt_iterator_next_block
subroutine, public dbt_iterator_next_block(iterator, ind_nd, blk_size, blk_offset)
iterate over nd blocks of an nd rank tensor, index only (blocks must be retrieved by calling dbt_get_...
Definition dbt_block.F:161
dbt_index
tensor index and mapping to DBM index
Definition dbt_index.F:12
dbt_index::get_2d_indices_tensor
pure integer(kind=int_8) function, dimension(2), public get_2d_indices_tensor(map, ind_in)
transform nd index to 2d index, using info from index mapping.
Definition dbt_index.F:318
dbt_index::dbt_inverse_order
pure integer function, dimension(size(order)), public dbt_inverse_order(order)
Invert order.
Definition dbt_index.F:410
dbt_index::combine_tensor_index
pure integer(kind=int_8) function, public combine_tensor_index(ind_in, dims)
transform nd index to flat index
Definition dbt_index.F:235
dbt_index::dbt_get_mapping_info
pure subroutine, public dbt_get_mapping_info(map, ndim_nd, ndim1_2d, ndim2_2d, dims_2d_i8, dims_2d, dims_nd, dims1_2d, dims2_2d, map1_2d, map2_2d, map_nd, base, col_major)
get mapping info
Definition dbt_index.F:176
dbt_io
DBT tensor Input / Output.
Definition dbt_io.F:12
dbt_io::dbt_write_blocks
subroutine, public dbt_write_blocks(tensor, io_unit_master, io_unit_all, write_int)
Write all tensor blocks.
Definition dbt_io.F:212
dbt_io::dbt_write_block_indices
subroutine, public dbt_write_block_indices(tensor, io_unit_master, io_unit_all)
Definition dbt_io.F:365
dbt_methods
DBT tensor framework for block-sparse tensor contraction. Representation of n-rank tensors as DBT tal...
Definition dbt_methods.F:16
dbt_methods::dbt_copy
subroutine, public dbt_copy(tensor_in, tensor_out, order, summation, bounds, move_data, unit_nr)
Copy tensor data. Redistributes tensor data according to distributions of target and source tensor....
Definition dbt_methods.F:114
dbt_methods::dbt_contract
subroutine, public dbt_contract(alpha, tensor_1, tensor_2, beta, tensor_3, contract_1, notcontract_1, contract_2, notcontract_2, map_1, map_2, bounds_1, bounds_2, bounds_3, optimize_dist, pgrid_opt_1, pgrid_opt_2, pgrid_opt_3, filter_eps, flop, move_data, retain_sparsity, unit_nr, log_verbose)
Contract tensors by multiplying matrix representations. tensor_3(map_1, map_2) := alpha * tensor_1(no...
Definition dbt_methods.F:499
dbt_tas_base
Tall-and-skinny matrices: base routines similar to DBM API, mostly wrappers around existing DBM routi...
Definition dbt_tas_base.F:13
dbt_tas_base::dbt_tas_info
type(dbt_tas_split_info) function, public dbt_tas_info(matrix)
get info on mpi grid splitting
Definition dbt_tas_base.F:822
dbt_tas_test
testing infrastructure for tall-and-skinny matrices
Definition dbt_tas_test.F:12
dbt_tas_test::dbt_tas_checksum
real(kind=dp) function, public dbt_tas_checksum(matrix)
Calculate checksum of tall-and-skinny matrix consistent with dbm_checksum.
Definition dbt_tas_test.F:475
dbt_test
General methods for testing DBT tensors.
Definition dbt_test.F:12
dbt_test::dbt_checksum
real(kind=dp) function, public dbt_checksum(tensor)
checksum of a tensor consistent with block_checksum
Definition dbt_test.F:1185
dbt_test::dbt_reset_randmat_seed
subroutine, public dbt_reset_randmat_seed()
Reset the seed used for generating random matrices to default value.
Definition dbt_test.F:1195
dbt_test::dbt_test_formats
subroutine, public dbt_test_formats(ndims, mp_comm, unit_nr, verbose, blk_size_1, blk_size_2, blk_size_3, blk_size_4, blk_ind_1, blk_ind_2, blk_ind_3, blk_ind_4)
Test equivalence of all tensor formats, using a random distribution.
Definition dbt_test.F:176
dbt_test::dbt_contract_test
subroutine, public dbt_contract_test(alpha, tensor_1, tensor_2, beta, tensor_3, contract_1, notcontract_1, contract_2, notcontract_2, map_1, map_2, unit_nr, bounds_1, bounds_2, bounds_3, log_verbose, write_int)
test tensor contraction
Definition dbt_test.F:833
dbt_test::dbt_setup_test_tensor
subroutine, public dbt_setup_test_tensor(tensor, mp_comm, enumerate, blk_ind_1, blk_ind_2, blk_ind_3, blk_ind_4)
Allocate and fill test tensor - entries are enumerated by their index s.t. they only depend on global...
Definition dbt_test.F:476
dbt_types
DBT tensor framework for block-sparse tensor contraction: Types and create/destroy routines.
Definition dbt_types.F:12
dbt_types::dbt_pgrid_destroy
subroutine, public dbt_pgrid_destroy(pgrid, keep_comm)
destroy process grid
Definition dbt_types.F:905
dbt_types::dbt_distribution_new
subroutine, public dbt_distribution_new(dist, pgrid, nd_dist_1, nd_dist_2, nd_dist_3, nd_dist_4)
Create a tensor distribution.
Definition dbt_types.F:886
dbt_types::dims_tensor
subroutine, public dims_tensor(tensor, dims)
tensor dimensions
Definition dbt_types.F:1238
dbt_types::dbt_destroy
subroutine, public dbt_destroy(tensor)
Destroy a tensor.
Definition dbt_types.F:1410
dbt_types::dbt_get_info
subroutine, public dbt_get_info(tensor, nblks_total, nfull_total, nblks_local, nfull_local, pdims, my_ploc, blks_local_1, blks_local_2, blks_local_3, blks_local_4, proc_dist_1, proc_dist_2, proc_dist_3, proc_dist_4, blk_size_1, blk_size_2, blk_size_3, blk_size_4, blk_offset_1, blk_offset_2, blk_offset_3, blk_offset_4, distribution, name)
As block_get_info but for tensors.
Definition dbt_types.F:1656
dbt_types::ndims_tensor
pure integer function, public ndims_tensor(tensor)
tensor rank
Definition dbt_types.F:1227
dbt_types::dbt_default_distvec
subroutine, public dbt_default_distvec(nblk, nproc, blk_size, dist)
get a load-balanced and randomized distribution along one tensor dimension
Definition dbt_types.F:1876
dbt_types::mp_environ_pgrid
subroutine, public mp_environ_pgrid(pgrid, dims, task_coor)
as mp_environ but for special pgrid type
Definition dbt_types.F:768
dbt_types::dbt_get_stored_coordinates
subroutine, public dbt_get_stored_coordinates(tensor, ind_nd, processor)
Generalization of block_get_stored_coordinates for tensors.
Definition dbt_types.F:1510
dbt_types::dbt_pgrid_create
subroutine, public dbt_pgrid_create(mp_comm, dims, pgrid, tensor_dims)
Definition dbt_types.F:1525
dbt_types::dbt_distribution_destroy
subroutine, public dbt_distribution_destroy(dist)
Destroy tensor distribution.
Definition dbt_types.F:926
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::int_8
integer, parameter, public int_8
Definition kinds.F:54
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
dbt_block::block_nd
Definition dbt_block.F:68
dbt_block::dbt_iterator_type
Definition dbt_block.F:63
dbt_types::dbt_distribution_type
Definition dbt_types.F:129
dbt_types::dbt_pgrid_type
Definition dbt_types.F:96
dbt_types::dbt_type
Definition dbt_types.F:110
message_passing::mp_comm_type
Definition message_passing.F:189