(git:0e87843)
Loading...
Searching...
No Matches
cp_dbcsr_operations.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief DBCSR operations in CP2K
10!> \author Urban Borstnik
11!> \date 2009-05-12
12!> \version 0.8
13!>
14!> <b>Modification history:</b>
15!> - Created 2009-05-12
16!> - Generalized sm_fm_mulitply for matrices w/ different row/col block size (A. Bussy, 11.2018)
17! **************************************************************************************************
20 USE cp_dbcsr_api, ONLY: &
21 dbcsr_add, dbcsr_complete_redistribute, dbcsr_convert_sizes_to_offsets, dbcsr_copy, &
27 dbcsr_scale, dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, &
28 dbcsr_type_symmetric, dbcsr_valid_index, dbcsr_verify_matrix
35 USE cp_fm_types, ONLY: cp_fm_create,&
42 USE kinds, ONLY: default_string_length,&
43 dp,&
44 int_8
45 USE mathlib, ONLY: gcd,&
46 lcm
48
49!$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads
50#include "base/base_uses.f90"
51
52 IMPLICIT NONE
53 PRIVATE
54
55 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_dbcsr_operations'
56 LOGICAL, PARAMETER :: debug_mod = .false.
57
58 INTEGER, SAVE, PUBLIC :: max_elements_per_block = 32
59
60 ! A plan caches S/A routing for a source block layout and FM distribution.
61 ! Use check_plan for sources whose layout can change. The owner must release
62 ! the plan before releasing the FM structure. All owned storage is allocatable;
63 ! no matrix data pointers or MPI requests are retained. Calls must be serialized.
64 TYPE dbcsr_fm_thread_plan_type
65 INTEGER, ALLOCATABLE :: blocks(:, :), slots(:)
66 END TYPE dbcsr_fm_thread_plan_type
67
69 PRIVATE
70 TYPE(cp_fm_struct_type), POINTER :: matrix_struct => null()
71 CHARACTER :: matrix_type = " "
72 TYPE(dbcsr_fm_thread_plan_type), ALLOCATABLE :: thread(:)
73 INTEGER, ALLOCATABLE :: send_count(:), send_disp(:), recv_count(:), recv_disp(:), recv_index(:, :)
74 REAL(kind=dp), ALLOCATABLE :: send_values(:), recv_values(:)
76
77 PUBLIC :: dbcsr_multiply_local
78
79 ! CP2K API emulation
85
86 ! distribution_2d_type compatibility
88
90
91 ! matrix set
94
96 MODULE PROCEDURE allocate_dbcsr_matrix_set_1d
97 MODULE PROCEDURE allocate_dbcsr_matrix_set_2d
98 MODULE PROCEDURE allocate_dbcsr_matrix_set_3d
99 MODULE PROCEDURE allocate_dbcsr_matrix_set_4d
100 MODULE PROCEDURE allocate_dbcsr_matrix_set_5d
101 END INTERFACE
102
104 MODULE PROCEDURE deallocate_dbcsr_matrix_set_1d
105 MODULE PROCEDURE deallocate_dbcsr_matrix_set_2d
106 MODULE PROCEDURE deallocate_dbcsr_matrix_set_3d
107 MODULE PROCEDURE deallocate_dbcsr_matrix_set_4d
108 MODULE PROCEDURE deallocate_dbcsr_matrix_set_5d
109 END INTERFACE
110
111CONTAINS
112
113! **************************************************************************************************
114!> \brief Copy a BLACS matrix to a dbcsr matrix.
115!>
116!> real_matrix=beta*real_matrix+alpha*fm
117!> beta defaults to 0, alpha to 1
118!> \param[in] fm full matrix
119!> \param[out] matrix DBCSR matrix
120!> \param[in] keep_sparsity (optional) retains the sparsity of the input
121!> matrix
122!> \date 2009-10-13
123!> \par History
124!> 2009-10-13 rewritten based on copy_dbcsr_to_fm
125!> \author Urban Borstnik
126!> \version 2.0
127! **************************************************************************************************
128 SUBROUTINE copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
129 TYPE(cp_fm_type), INTENT(IN) :: fm
130 TYPE(dbcsr_type), INTENT(INOUT) :: matrix
131 LOGICAL, INTENT(IN), OPTIONAL :: keep_sparsity
132
133 CHARACTER(LEN=*), PARAMETER :: routinen = 'copy_fm_to_dbcsr'
134
135 INTEGER :: handle
136 LOGICAL :: my_keep_sparsity
137 TYPE(dbcsr_type) :: bc_mat, redist_mat
138
139 CALL timeset(routinen, handle)
140
141 my_keep_sparsity = .false.
142 IF (PRESENT(keep_sparsity)) my_keep_sparsity = keep_sparsity
143
144 CALL copy_fm_to_dbcsr_bc(fm, bc_mat)
145
146 IF (my_keep_sparsity) THEN
147 CALL dbcsr_create(redist_mat, template=matrix)
148 CALL dbcsr_complete_redistribute(bc_mat, redist_mat)
149 CALL dbcsr_copy(matrix, redist_mat, keep_sparsity=.true.)
150 CALL dbcsr_release(redist_mat)
151 ELSE
152 CALL dbcsr_complete_redistribute(bc_mat, matrix)
153 END IF
154
155 CALL dbcsr_release(bc_mat)
156
157 CALL timestop(handle)
158 END SUBROUTINE copy_fm_to_dbcsr
159
160! **************************************************************************************************
161!> \brief Copy a BLACS matrix to a dbcsr matrix with a special block-cyclic distribution,
162!> which requires no complete redistribution.
163!> \param fm ...
164!> \param bc_mat ...
165! **************************************************************************************************
166 SUBROUTINE copy_fm_to_dbcsr_bc(fm, bc_mat)
167 TYPE(cp_fm_type), INTENT(IN) :: fm
168 TYPE(dbcsr_type), INTENT(INOUT) :: bc_mat
169
170 CHARACTER(LEN=*), PARAMETER :: routinen = 'copy_fm_to_dbcsr_bc'
171
172 INTEGER :: col, handle, ncol_block, ncol_global, &
173 nrow_block, nrow_global, row
174 INTEGER, ALLOCATABLE, DIMENSION(:) :: first_col, first_row, last_col, last_row
175 INTEGER, DIMENSION(:), POINTER :: col_blk_size, row_blk_size
176 INTEGER, DIMENSION(:, :), POINTER :: pgrid
177 REAL(kind=dp), DIMENSION(:, :), POINTER :: dbcsr_block, fm_block
178 TYPE(dbcsr_distribution_type) :: bc_dist
179 TYPE(dbcsr_iterator_type) :: iter
180
181 CALL timeset(routinen, handle)
182
183 ! Create processor grid
184 pgrid => fm%matrix_struct%context%blacs2mpi
185
186 ! Create a block-cyclic distribution compatible with the FM matrix.
187 nrow_block = fm%matrix_struct%nrow_block
188 ncol_block = fm%matrix_struct%ncol_block
189 nrow_global = fm%matrix_struct%nrow_global
190 ncol_global = fm%matrix_struct%ncol_global
191 NULLIFY (col_blk_size, row_blk_size)
192 CALL dbcsr_create_dist_block_cyclic(bc_dist, &
193 nrows=nrow_global, ncolumns=ncol_global, & ! Actual full matrix size
194 nrow_block=nrow_block, ncol_block=ncol_block, & ! BLACS parameters
195 group_handle=fm%matrix_struct%para_env%get_handle(), pgrid=pgrid, &
196 row_blk_sizes=row_blk_size, col_blk_sizes=col_blk_size) ! block-cyclic row/col sizes
197
198 ! Create the block-cyclic DBCSR matrix
199 CALL dbcsr_create(bc_mat, "Block-cyclic ", bc_dist, &
200 dbcsr_type_no_symmetry, row_blk_size, col_blk_size, reuse_arrays=.true.)
201 CALL dbcsr_distribution_release(bc_dist)
202
203 ! allocate all blocks
204 CALL dbcsr_reserve_all_blocks(bc_mat)
205
206 CALL calculate_fm_block_ranges(bc_mat, first_row, last_row, first_col, last_col)
207
208 ! Copy the FM data to the block-cyclic DBCSR matrix. This step
209 ! could be skipped with appropriate DBCSR index manipulation.
210 fm_block => fm%local_data
211!$OMP PARALLEL DEFAULT(NONE) PRIVATE(iter, row, col, dbcsr_block) &
212!$OMP SHARED(bc_mat, last_row, first_row, last_col, first_col, fm_block)
213 CALL dbcsr_iterator_start(iter, bc_mat)
214 DO WHILE (dbcsr_iterator_blocks_left(iter))
215 CALL dbcsr_iterator_next_block(iter, row, col, dbcsr_block)
216 dbcsr_block(:, :) = fm_block(first_row(row):last_row(row), first_col(col):last_col(col))
217 END DO
218 CALL dbcsr_iterator_stop(iter)
219!$OMP END PARALLEL
220
221 CALL timestop(handle)
222 END SUBROUTINE copy_fm_to_dbcsr_bc
223
224! **************************************************************************************************
225!> \brief Copy a DBCSR matrix to a BLACS matrix
226!> \param[in] matrix DBCSR matrix
227!> \param[out] fm full matrix
228!> \param plan optional S/A routing plan; discard before releasing the FM structure
229!> \param check_plan check for layout changes and collectively rebuild the plan;
230!> defaults to false for callers that guarantee fixed layouts and threads
231! **************************************************************************************************
232 SUBROUTINE copy_dbcsr_to_fm(matrix, fm, plan, check_plan)
233 TYPE(dbcsr_type), INTENT(IN) :: matrix
234 TYPE(cp_fm_type), INTENT(INOUT) :: fm
235 TYPE(dbcsr_fm_copy_plan_type), INTENT(INOUT), &
236 OPTIONAL :: plan
237 LOGICAL, INTENT(IN), OPTIONAL :: check_plan
238
239 CHARACTER(LEN=*), PARAMETER :: routinen = 'copy_dbcsr_to_fm'
240
241 CHARACTER(len=default_string_length) :: name
242 INTEGER :: group_handle, handle, ncol_block, &
243 nfullcols_total, nfullrows_total, &
244 nrow_block
245 INTEGER, DIMENSION(:), POINTER :: col_blk_size, row_blk_size
246 INTEGER, DIMENSION(:, :), POINTER :: pgrid
247 LOGICAL :: validate_plan
248 TYPE(dbcsr_distribution_type) :: bc_dist, dist
249 TYPE(dbcsr_fm_copy_plan_type) :: local_plan
250 TYPE(dbcsr_type) :: bc_mat, matrix_nosym
251
252 CALL timeset(routinen, handle)
253 validate_plan = .false.
254 IF (PRESENT(check_plan)) validate_plan = check_plan
255
256 ! check compatibility
257 CALL dbcsr_get_info(matrix, &
258 name=name, &
259 distribution=dist, &
260 nfullrows_total=nfullrows_total, &
261 nfullcols_total=nfullcols_total)
262
263 cpassert(fm%matrix_struct%nrow_global == nfullrows_total)
264 cpassert(fm%matrix_struct%ncol_global == nfullcols_total)
265
266 SELECT CASE (dbcsr_get_matrix_type(matrix))
267 CASE (dbcsr_type_symmetric, dbcsr_type_antisymmetric)
268 IF (PRESENT(plan)) THEN
269 CALL copy_dbcsr_to_fm_sym(matrix, fm, plan, validate_plan)
270 ELSE
271 CALL copy_dbcsr_to_fm_sym(matrix, fm, local_plan, .false.)
272 END IF
273 CALL timestop(handle)
274 RETURN
275 END SELECT
276
277 ! Other matrix types use the uncached redistribution path.
278 IF (PRESENT(plan)) plan = dbcsr_fm_copy_plan_type()
279
280 ! info about the full matrix
281 nrow_block = fm%matrix_struct%nrow_block
282 ncol_block = fm%matrix_struct%ncol_block
283
284 ! Convert DBCSR to a block-cyclic
285 NULLIFY (col_blk_size, row_blk_size)
286 CALL dbcsr_distribution_get(dist, group=group_handle, pgrid=pgrid)
287 CALL dbcsr_create_dist_block_cyclic(bc_dist, &
288 nrows=nfullrows_total, ncolumns=nfullcols_total, &
289 nrow_block=nrow_block, ncol_block=ncol_block, &
290 group_handle=group_handle, pgrid=pgrid, &
291 row_blk_sizes=row_blk_size, col_blk_sizes=col_blk_size)
292
293 CALL dbcsr_create(bc_mat, "Block-cyclic"//name, bc_dist, &
294 dbcsr_type_no_symmetry, row_blk_size, col_blk_size, reuse_arrays=.true.)
295 CALL dbcsr_distribution_release(bc_dist)
296
297 IF (dbcsr_get_matrix_type(matrix) == dbcsr_type_no_symmetry) THEN
298 CALL dbcsr_complete_redistribute(matrix, bc_mat)
299 ELSE
300 CALL dbcsr_create(matrix_nosym, template=matrix, matrix_type=dbcsr_type_no_symmetry)
301 CALL dbcsr_desymmetrize(matrix, matrix_nosym)
302 CALL dbcsr_complete_redistribute(matrix_nosym, bc_mat)
303 CALL dbcsr_release(matrix_nosym)
304 END IF
305
306 CALL copy_dbcsr_to_fm_bc(bc_mat, fm)
307
308 CALL dbcsr_release(bc_mat)
309
310 CALL timestop(handle)
311 END SUBROUTINE copy_dbcsr_to_fm
312
313! **************************************************************************************************
314!> \brief Send S/A values using a routing plan, checking the block sequence while packing.
315!> Mutable sources require a collective check before using the cached routing.
316!> \param matrix symmetric or antisymmetric DBCSR workspace
317!> \param fm full matrix on the source communicator
318!> \param plan reusable routing and buffers; no outstanding communication on return
319!> \param check_plan rebuild collectively if the source layout, FM structure or threads change
320! **************************************************************************************************
321 SUBROUTINE copy_dbcsr_to_fm_sym(matrix, fm, plan, check_plan)
322 TYPE(dbcsr_type), INTENT(IN) :: matrix
323 TYPE(cp_fm_type), INTENT(INOUT) :: fm
324 TYPE(dbcsr_fm_copy_plan_type), INTENT(INOUT) :: plan
325 LOGICAL, INTENT(IN) :: check_plan
326
327 INTEGER :: attempt, col, col_offset, i, iblock, j, &
328 k, mirror, ncol, nrecv, nrow, &
329 nthreads, pos, rebuild, row, &
330 row_offset, thread
331 LOGICAL :: block_matches, layout_matches, transposed
332 REAL(kind=dp) :: symmetry_sign, value
333 REAL(kind=dp), DIMENSION(:, :), POINTER :: block
334 TYPE(dbcsr_iterator_type) :: iter
335 TYPE(mp_para_env_type), POINTER :: para_env
336
337 IF (.NOT. ALLOCATED(plan%thread)) CALL dbcsr_to_fm_plan_create(matrix, fm, plan)
338 para_env => fm%matrix_struct%context%para_env
339 symmetry_sign = 1.0_dp
340 IF (dbcsr_get_matrix_type(matrix) == dbcsr_type_antisymmetric) symmetry_sign = -1.0_dp
341
342 DO attempt = 1, 2
343 layout_matches = ASSOCIATED(plan%matrix_struct, fm%matrix_struct) .AND. &
344 plan%matrix_type == dbcsr_get_matrix_type(matrix)
345 rebuild = 0
346!$OMP PARALLEL DEFAULT(NONE) SHARED(matrix, plan, symmetry_sign, layout_matches) &
347!$OMP PRIVATE(iter, row, col, block, transposed, row_offset, col_offset, nrow, ncol, &
348!$OMP thread, nthreads, iblock, pos, mirror, i, j, k, value, block_matches) REDUCTION(MAX:rebuild)
349 thread = 0
350 nthreads = 1
351!$ thread = omp_get_thread_num()
352!$ nthreads = omp_get_num_threads()
353 IF (.NOT. layout_matches .OR. nthreads /= SIZE(plan%thread)) THEN
354 rebuild = 1
355 ELSE
356 iblock = 0
357 pos = 0
358 CALL dbcsr_iterator_readonly_start(iter, matrix, dynamic=.false.)
359 DO WHILE (dbcsr_iterator_blocks_left(iter))
360 CALL dbcsr_iterator_next_block(iter, row, col, block, transposed=transposed, &
361 row_size=nrow, col_size=ncol, row_offset=row_offset, col_offset=col_offset)
362 iblock = iblock + 1
363 block_matches = iblock <= SIZE(plan%thread(thread)%blocks, 2)
364 IF (block_matches) block_matches = all(plan%thread(thread)%blocks(:, iblock) == &
365 [row, col, nrow, ncol, row_offset, col_offset])
366 IF (.NOT. block_matches) THEN
367 rebuild = 1
368 EXIT
369 END IF
370 DO j = 1, ncol
371 DO i = 1, nrow
372 IF (transposed) THEN
373 value = symmetry_sign*block(j, i)
374 ELSE
375 value = block(i, j)
376 END IF
377 DO mirror = 1, 2
378 IF (mirror == 2 .AND. row == col) cycle
379 pos = pos + 1
380 k = plan%thread(thread)%slots(pos)
381 plan%send_values(k) = value
382 value = symmetry_sign*value
383 END DO
384 END DO
385 END DO
386 END DO
387 CALL dbcsr_iterator_stop(iter)
388 IF (iblock /= SIZE(plan%thread(thread)%blocks, 2) .OR. &
389 pos /= SIZE(plan%thread(thread)%slots)) rebuild = 1
390 END IF
391!$OMP END PARALLEL
392 IF (check_plan) CALL para_env%max(rebuild)
393 IF (rebuild == 0) EXIT
394 cpassert(check_plan .AND. attempt == 1)
395 CALL dbcsr_to_fm_plan_create(matrix, fm, plan)
396 END DO
397
398 ! Counts and destination indices were exchanged once when the plan was built.
399 CALL para_env%alltoall(plan%send_values, plan%send_count, plan%send_disp, &
400 plan%recv_values, plan%recv_count, plan%recv_disp)
401 nrecv = sum(plan%recv_count)
402 fm%local_data = 0.0_dp
403!$OMP PARALLEL DO DEFAULT(NONE) SHARED(fm, plan, nrecv) PRIVATE(k)
404 DO k = 1, nrecv
405 fm%local_data(plan%recv_index(1, k), plan%recv_index(2, k)) = plan%recv_values(k)
406 END DO
407!$OMP END PARALLEL DO
408
409 END SUBROUTINE copy_dbcsr_to_fm_sym
410
411! **************************************************************************************************
412!> \brief Cache source-to-packet slots and destination FM indices, not source pointers.
413!> Counting and indexing use the same static DBCSR thread assignment. The
414!> plan is valid only while that assignment and both matrix layouts are fixed.
415!> \param matrix symmetric or antisymmetric DBCSR workspace
416!> \param fm destination full matrix
417!> \param plan new routing plan
418! **************************************************************************************************
419 SUBROUTINE dbcsr_to_fm_plan_create(matrix, fm, plan)
420 TYPE(dbcsr_type), INTENT(IN) :: matrix
421 TYPE(cp_fm_type), INTENT(IN) :: fm
422 TYPE(dbcsr_fm_copy_plan_type), INTENT(OUT) :: plan
423
424 CHARACTER(LEN=*), PARAMETER :: routinen = 'dbcsr_to_fm_plan_create'
425
426 INTEGER :: col, col_offset, dest, handle, i, &
427 iblock, j, k, mirror, ncol, nproc, &
428 nrow, nthreads, pass, pos, row, &
429 row_offset, thread, x, y
430 INTEGER(KIND=int_8) :: count, offset
431 INTEGER(KIND=int_8), ALLOCATABLE, DIMENSION(:) :: local_pos
432 INTEGER(KIND=int_8), ALLOCATABLE, DIMENSION(:, :) :: thread_pos
433 INTEGER, ALLOCATABLE, DIMENSION(:) :: nblocks
434 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: send_index
435 LOGICAL :: transposed
436 TYPE(cp_fm_struct_type), POINTER :: fms
437 TYPE(dbcsr_iterator_type) :: iter
438 TYPE(mp_para_env_type), POINTER :: para_env
439
440 CALL timeset(routinen, handle)
441 fms => fm%matrix_struct
442 plan%matrix_struct => fms
443 plan%matrix_type = dbcsr_get_matrix_type(matrix)
444 para_env => fms%context%para_env
445 nproc = para_env%num_pe
446 nthreads = 1
447!$ nthreads = omp_get_max_threads()
448 ALLOCATE (thread_pos(0:nproc - 1, 0:nthreads - 1), nblocks(0:nthreads - 1))
449 ALLOCATE (plan%send_count(0:nproc - 1), plan%send_disp(0:nproc - 1), &
450 plan%recv_count(0:nproc - 1), plan%recv_disp(0:nproc - 1))
451 thread_pos(:, :) = 0
452 nblocks(:) = 0
453
454!$OMP PARALLEL DEFAULT(NONE) &
455!$OMP SHARED(matrix, fms, nproc, nthreads, &
456!$OMP thread_pos, nblocks, send_index, plan) &
457!$OMP PRIVATE(iter, row, col, row_offset, col_offset, nrow, ncol, &
458!$OMP thread, iblock, pos, pass, mirror, i, j, x, y, dest, k, offset, count, local_pos, transposed)
459 thread = 0
460!$ thread = omp_get_thread_num()
461 ALLOCATE (local_pos(0:nproc - 1))
462 local_pos(:) = 0
463 DO pass = 1, 2
464 iblock = 0
465 pos = 0
466 IF (pass == 2) THEN
467!$OMP SINGLE
468 cpassert(sum(thread_pos) <= int(huge(0), int_8)/2)
469 plan%send_count(:) = int(sum(thread_pos, dim=2))
470!$ nthreads = omp_get_num_threads()
471 ALLOCATE (plan%thread(0:nthreads - 1))
472 DO k = 0, nthreads - 1
473 ALLOCATE (plan%thread(k)%blocks(6, nblocks(k)), &
474 plan%thread(k)%slots(int(sum(thread_pos(:, k)))))
475 END DO
476 offset = 0
477 DO dest = 0, nproc - 1
478 plan%send_disp(dest) = int(offset)
479 DO k = 0, nthreads - 1
480 count = thread_pos(dest, k)
481 thread_pos(dest, k) = offset
482 offset = offset + count
483 END DO
484 END DO
485 ALLOCATE (send_index(2, max(1, sum(plan%send_count))), &
486 plan%send_values(max(1, sum(plan%send_count))))
487 send_index(:, :) = 0
488 plan%send_values(:) = 0.0_dp
489!$OMP END SINGLE
490 local_pos(:) = thread_pos(:, thread)
491 END IF
492 CALL dbcsr_iterator_readonly_start(iter, matrix, dynamic=.false.)
493 DO WHILE (dbcsr_iterator_blocks_left(iter))
494 CALL dbcsr_iterator_next_block(iter, row, col, row_size=nrow, col_size=ncol, &
495 row_offset=row_offset, col_offset=col_offset, transposed=transposed)
496 iblock = iblock + 1
497 IF (pass == 2) THEN
498 plan%thread(thread)%blocks(:, iblock) = [row, col, nrow, ncol, row_offset, col_offset]
499 END IF
500 DO j = 1, ncol
501 DO i = 1, nrow
502 DO mirror = 1, 2
503 IF (mirror == 2 .AND. row == col) cycle
504 x = row_offset + i - 1
505 y = col_offset + j - 1
506 IF (mirror == 2) THEN
507 x = col_offset + j - 1
508 y = row_offset + i - 1
509 END IF
510 dest = fms%context%blacs2mpi(fms%g2p_row(x), fms%g2p_col(y))
511 local_pos(dest) = local_pos(dest) + 1
512 IF (pass == 1) cycle
513 k = int(local_pos(dest))
514 send_index(:, k) = [fms%g2l_row(x), fms%g2l_col(y)]
515 pos = pos + 1
516 plan%thread(thread)%slots(pos) = k
517 END DO
518 END DO
519 END DO
520 END DO
521 CALL dbcsr_iterator_stop(iter)
522 IF (pass == 1) THEN
523 thread_pos(:, thread) = local_pos(:)
524 nblocks(thread) = iblock
525 END IF
526!$OMP BARRIER
527 END DO
528 DEALLOCATE (local_pos)
529!$OMP END PARALLEL
530 DEALLOCATE (thread_pos, nblocks)
531
532 CALL para_env%alltoall(plan%send_count, plan%recv_count, 1)
533 cpassert(sum(int(plan%recv_count, int_8)) <= int(huge(0), int_8)/2)
534 plan%recv_disp(0) = 0
535 DO dest = 1, nproc - 1
536 plan%recv_disp(dest) = plan%recv_disp(dest - 1) + plan%recv_count(dest - 1)
537 END DO
538 ALLOCATE (plan%recv_index(2, max(1, sum(plan%recv_count))), &
539 plan%recv_values(max(1, sum(plan%recv_count))))
540 plan%recv_index(:, :) = 0
541 plan%recv_values(:) = 0.0_dp
542 CALL para_env%alltoall(send_index, 2*plan%send_count, 2*plan%send_disp, &
543 plan%recv_index, 2*plan%recv_count, 2*plan%recv_disp)
544 CALL timestop(handle)
545
546 END SUBROUTINE dbcsr_to_fm_plan_create
547
548! **************************************************************************************************
549!> \brief Copy a DBCSR_BLACS matrix to a BLACS matrix
550!> \param bc_mat DBCSR matrix
551!> \param[out] fm full matrix
552! **************************************************************************************************
553 SUBROUTINE copy_dbcsr_to_fm_bc(bc_mat, fm)
554 TYPE(dbcsr_type), INTENT(IN) :: bc_mat
555 TYPE(cp_fm_type), INTENT(INOUT) :: fm
556
557 CHARACTER(LEN=*), PARAMETER :: routinen = 'copy_dbcsr_to_fm_bc'
558
559 INTEGER :: col, handle, row
560 INTEGER, ALLOCATABLE, DIMENSION(:) :: first_col, first_row, last_col, last_row
561 REAL(kind=dp), DIMENSION(:, :), POINTER :: dbcsr_block, fm_block
562 TYPE(dbcsr_iterator_type) :: iter
563
564 CALL timeset(routinen, handle)
565
566 CALL calculate_fm_block_ranges(bc_mat, first_row, last_row, first_col, last_col)
567
568 ! Now copy data to the FM matrix
569 fm_block => fm%local_data
570 fm_block = real(0.0, kind=dp)
571!$OMP PARALLEL DEFAULT(NONE) PRIVATE(iter, row, col, dbcsr_block) &
572!$OMP SHARED(bc_mat, last_row, first_row, last_col, first_col, fm_block)
573 CALL dbcsr_iterator_readonly_start(iter, bc_mat)
574 DO WHILE (dbcsr_iterator_blocks_left(iter))
575 CALL dbcsr_iterator_next_block(iter, row, col, dbcsr_block)
576 fm_block(first_row(row):last_row(row), first_col(col):last_col(col)) = dbcsr_block(:, :)
577 END DO
578 CALL dbcsr_iterator_stop(iter)
579!$OMP END PARALLEL
580
581 CALL timestop(handle)
582 END SUBROUTINE copy_dbcsr_to_fm_bc
583
584! **************************************************************************************************
585!> \brief Helper routine used to copy blocks from DBCSR into FM matrices and vice versa
586!> \param bc_mat ...
587!> \param first_row ...
588!> \param last_row ...
589!> \param first_col ...
590!> \param last_col ...
591!> \author Ole Schuett
592! **************************************************************************************************
593 SUBROUTINE calculate_fm_block_ranges(bc_mat, first_row, last_row, first_col, last_col)
594 TYPE(dbcsr_type), INTENT(IN) :: bc_mat
595 INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: first_row, last_row, first_col, last_col
596
597 INTEGER :: col, nblkcols_local, nblkcols_total, &
598 nblkrows_local, nblkrows_total, row
599 INTEGER, ALLOCATABLE, DIMENSION(:) :: local_col_sizes, local_row_sizes
600 INTEGER, DIMENSION(:), POINTER :: col_blk_size, local_cols, local_rows, &
601 row_blk_size
602
603 CALL dbcsr_get_info(bc_mat, &
604 nblkrows_total=nblkrows_total, &
605 nblkcols_total=nblkcols_total, &
606 nblkrows_local=nblkrows_local, &
607 nblkcols_local=nblkcols_local, &
608 local_rows=local_rows, &
609 local_cols=local_cols, &
610 row_blk_size=row_blk_size, &
611 col_blk_size=col_blk_size)
612
613 ! calculate first_row and last_row
614 ALLOCATE (local_row_sizes(nblkrows_total))
615 local_row_sizes(:) = 0
616 IF (nblkrows_local >= 1) THEN
617 DO row = 1, nblkrows_local
618 local_row_sizes(local_rows(row)) = row_blk_size(local_rows(row))
619 END DO
620 END IF
621 ALLOCATE (first_row(nblkrows_total), last_row(nblkrows_total))
622 CALL dbcsr_convert_sizes_to_offsets(local_row_sizes, first_row, last_row)
623 DEALLOCATE (local_row_sizes)
624
625 ! calculate first_col and last_col
626 ALLOCATE (local_col_sizes(nblkcols_total))
627 local_col_sizes(:) = 0
628 IF (nblkcols_local >= 1) THEN
629 DO col = 1, nblkcols_local
630 local_col_sizes(local_cols(col)) = col_blk_size(local_cols(col))
631 END DO
632 END IF
633 ALLOCATE (first_col(nblkcols_total), last_col(nblkcols_total))
634 CALL dbcsr_convert_sizes_to_offsets(local_col_sizes, first_col, last_col)
635 DEALLOCATE (local_col_sizes)
636
637 END SUBROUTINE calculate_fm_block_ranges
638
639! **************************************************************************************************
640!> \brief hack for dbcsr_copy_columns
641!> \param matrix_b ...
642!> \param matrix_a ...
643!> \param ncol ...
644!> \param source_start ...
645!> \param target_start ...
646!> \param para_env ...
647!> \param blacs_env ...
648!> \author vw
649! **************************************************************************************************
650 SUBROUTINE dbcsr_copy_columns_hack(matrix_b, matrix_a, &
651 ncol, source_start, target_start, para_env, blacs_env)
652
653 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_b
654 TYPE(dbcsr_type), INTENT(IN) :: matrix_a
655 INTEGER, INTENT(IN) :: ncol, source_start, target_start
656 TYPE(mp_para_env_type), POINTER :: para_env
657 TYPE(cp_blacs_env_type), POINTER :: blacs_env
658
659 INTEGER :: nfullcols_total, nfullrows_total
660 TYPE(cp_fm_struct_type), POINTER :: fm_struct
661 TYPE(cp_fm_type) :: fm_matrix_a, fm_matrix_b
662
663 NULLIFY (fm_struct)
664 CALL dbcsr_get_info(matrix_a, nfullrows_total=nfullrows_total, nfullcols_total=nfullcols_total)
665 CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=nfullrows_total, &
666 ncol_global=nfullcols_total, para_env=para_env)
667 CALL cp_fm_create(fm_matrix_a, fm_struct, name="fm_matrix_a")
668 CALL cp_fm_struct_release(fm_struct)
669
670 CALL dbcsr_get_info(matrix_b, nfullrows_total=nfullrows_total, nfullcols_total=nfullcols_total)
671 CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=nfullrows_total, &
672 ncol_global=nfullcols_total, para_env=para_env)
673 CALL cp_fm_create(fm_matrix_b, fm_struct, name="fm_matrix_b")
674 CALL cp_fm_struct_release(fm_struct)
675
676 CALL copy_dbcsr_to_fm(matrix_a, fm_matrix_a)
677 CALL copy_dbcsr_to_fm(matrix_b, fm_matrix_b)
678
679 CALL cp_fm_to_fm(fm_matrix_a, fm_matrix_b, ncol, source_start, target_start)
680
681 CALL copy_fm_to_dbcsr(fm_matrix_b, matrix_b)
682
683 CALL cp_fm_release(fm_matrix_a)
684 CALL cp_fm_release(fm_matrix_b)
685
686 END SUBROUTINE dbcsr_copy_columns_hack
687
688! **************************************************************************************************
689!> \brief Creates a DBCSR distribution from a distribution_2d
690!> \param[in] dist2d distribution_2d
691!> \param[out] dist DBCSR distribution
692!> \par History
693!> move form dbcsr_operation 01.2010
694! **************************************************************************************************
695 SUBROUTINE cp_dbcsr_dist2d_to_dist(dist2d, dist)
696 TYPE(distribution_2d_type), INTENT(IN), TARGET :: dist2d
697 TYPE(dbcsr_distribution_type), INTENT(OUT) :: dist
698
699 INTEGER, DIMENSION(:), POINTER :: col_dist, row_dist
700 INTEGER, DIMENSION(:, :), POINTER :: col_dist_data, pgrid, row_dist_data
701 TYPE(cp_blacs_env_type), POINTER :: blacs_env
702 TYPE(distribution_2d_type), POINTER :: dist2d_p
703 TYPE(mp_para_env_type), POINTER :: para_env
704
705 dist2d_p => dist2d
706 CALL distribution_2d_get(dist2d_p, &
707 row_distribution=row_dist_data, &
708 col_distribution=col_dist_data, &
709 blacs_env=blacs_env)
710 CALL blacs_env%get(para_env=para_env, blacs2mpi=pgrid)
711
712 ! map to 1D arrays
713 row_dist => row_dist_data(:, 1)
714 col_dist => col_dist_data(:, 1)
715 !row_cluster => row_dist_data(:, 2)
716 !col_cluster => col_dist_data(:, 2)
717
718 CALL dbcsr_distribution_new(dist, &
719 group=para_env%get_handle(), pgrid=pgrid, &
720 row_dist=row_dist, &
721 col_dist=col_dist)
722
723 END SUBROUTINE cp_dbcsr_dist2d_to_dist
724
725! **************************************************************************************************
726!> \brief multiply a dbcsr with a replicated array
727!> c = alpha_scalar * A (dbscr) * b + c
728!> \param[in] matrix_a DBSCR matrxx
729!> \param[in] vec_b vectors b
730!> \param[inout] vec_c vectors c
731!> \param[in] ncol nbr of columns
732!> \param[in] alpha alpha
733!>
734! **************************************************************************************************
735 SUBROUTINE dbcsr_multiply_local(matrix_a, vec_b, vec_c, ncol, alpha)
736 TYPE(dbcsr_type), INTENT(IN) :: matrix_a
737 REAL(dp), DIMENSION(:, :), INTENT(IN) :: vec_b
738 REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: vec_c
739 INTEGER, INTENT(in), OPTIONAL :: ncol
740 REAL(dp), INTENT(IN), OPTIONAL :: alpha
741
742 CHARACTER(LEN=*), PARAMETER :: routinen = 'dbcsr_multiply_local'
743
744 INTEGER :: col, coloff, my_ncol, row, rowoff, &
745 timing_handle
746 LOGICAL :: has_symm
747 REAL(dp) :: my_alpha, my_alpha2
748 REAL(dp), DIMENSION(:, :), POINTER :: data_d
749 TYPE(dbcsr_iterator_type) :: iter
750
751 CALL timeset(routinen, timing_handle)
752
753 my_alpha = 1.0_dp
754 IF (PRESENT(alpha)) my_alpha = alpha
755
756 my_ncol = SIZE(vec_b, 2)
757 IF (PRESENT(ncol)) my_ncol = ncol
758
759 my_alpha2 = 0.0_dp
760 IF (dbcsr_get_matrix_type(matrix_a) == dbcsr_type_symmetric) my_alpha2 = my_alpha
761 IF (dbcsr_get_matrix_type(matrix_a) == dbcsr_type_antisymmetric) my_alpha2 = -my_alpha
762
763 has_symm = (dbcsr_get_matrix_type(matrix_a) == dbcsr_type_symmetric .OR. &
764 dbcsr_get_matrix_type(matrix_a) == dbcsr_type_antisymmetric)
765
766!$OMP PARALLEL DEFAULT(NONE) SHARED(matrix_a,vec_b,vec_c,ncol,my_alpha2,my_alpha,my_ncol,has_symm) &
767!$OMP PRIVATE(iter,row,col,data_d,rowoff,coloff)
768 CALL dbcsr_iterator_readonly_start(iter, matrix_a, dynamic=.true., dynamic_byrows=.true.)
769 DO WHILE (dbcsr_iterator_blocks_left(iter))
770 CALL dbcsr_iterator_next_block(iter, row, col, data_d, row_offset=rowoff, col_offset=coloff)
771 IF (my_ncol /= 1) THEN
772 CALL dgemm('N', 'N', &
773 SIZE(data_d, 1), my_ncol, SIZE(data_d, 2), &
774 my_alpha, data_d(1, 1), SIZE(data_d, 1), &
775 vec_b(coloff, 1), SIZE(vec_b, 1), &
776 1.0_dp, vec_c(rowoff, 1), SIZE(vec_c, 1))
777 ELSE
778 CALL dgemv('N', SIZE(data_d, 1), SIZE(data_d, 2), &
779 my_alpha, data_d(1, 1), SIZE(data_d, 1), &
780 vec_b(coloff, 1), 1, &
781 1.0_dp, vec_c(rowoff, 1), 1)
782 END IF
783 END DO
784 CALL dbcsr_iterator_stop(iter)
785!$OMP END PARALLEL
786
787 ! FIXME ... in the symmetric case, the writes to vec_c depend on the column, not the row. This makes OMP-ing more difficult
788 ! needs e.g. a buffer for vec_c and a reduction of that buffer.
789 IF (has_symm) THEN
790 CALL dbcsr_iterator_readonly_start(iter, matrix_a)
791 DO WHILE (dbcsr_iterator_blocks_left(iter))
792 CALL dbcsr_iterator_next_block(iter, row, col, data_d, row_offset=rowoff, col_offset=coloff)
793 IF (row /= col) THEN
794 IF (my_ncol /= 1) THEN
795 CALL dgemm('T', 'N', &
796 SIZE(data_d, 2), my_ncol, SIZE(data_d, 1), &
797 my_alpha2, data_d(1, 1), SIZE(data_d, 1), &
798 vec_b(rowoff, 1), SIZE(vec_b, 1), &
799 1.0_dp, vec_c(coloff, 1), SIZE(vec_c, 1))
800 ELSE
801 CALL dgemv('T', SIZE(data_d, 1), SIZE(data_d, 2), &
802 my_alpha2, data_d(1, 1), SIZE(data_d, 1), &
803 vec_b(rowoff, 1), 1, &
804 1.0_dp, vec_c(coloff, 1), 1)
805 END IF
806 END IF
807 END DO
808 CALL dbcsr_iterator_stop(iter)
809 END IF
810
811 CALL timestop(timing_handle)
812 END SUBROUTINE dbcsr_multiply_local
813
814! **************************************************************************************************
815!> \brief multiply a dbcsr with a fm matrix
816!>
817!> For backwards compatibility with BLAS XGEMM, this routine supports
818!> the multiplication of matrices with incompatible dimensions.
819!>
820!> \param[in] matrix DBCSR matrix
821!> \param fm_in full matrix
822!> \param fm_out full matrix
823!> \param[in] ncol nbr of columns
824!> \param[in] alpha alpha
825!> \param[in] beta beta
826!>
827! **************************************************************************************************
828 SUBROUTINE cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
829 TYPE(dbcsr_type), INTENT(IN) :: matrix
830 TYPE(cp_fm_type), INTENT(IN) :: fm_in
831 TYPE(cp_fm_type), INTENT(INOUT) :: fm_out
832 INTEGER, INTENT(IN) :: ncol
833 REAL(dp), INTENT(IN), OPTIONAL :: alpha, beta
834
835 CHARACTER(LEN=*), PARAMETER :: routinen = 'cp_dbcsr_sm_fm_multiply'
836
837 INTEGER :: a_ncol, a_nrow, b_ncol, b_nrow, c_ncol, &
838 c_nrow, k_in, k_out, timing_handle, &
839 timing_handle_mult
840 INTEGER, DIMENSION(:), POINTER :: col_blk_size, col_blk_size_right_in, &
841 col_blk_size_right_out, col_dist, &
842 row_blk_size, row_dist
843 TYPE(dbcsr_type) :: in, out
844 TYPE(dbcsr_distribution_type) :: dist, dist_right_in, product_dist
845 REAL(dp) :: my_alpha, my_beta
846
847 CALL timeset(routinen, timing_handle)
848
849 my_alpha = 1.0_dp
850 my_beta = 0.0_dp
851 IF (PRESENT(alpha)) my_alpha = alpha
852 IF (PRESENT(beta)) my_beta = beta
853
854 ! TODO
855 CALL cp_fm_get_info(fm_in, ncol_global=b_ncol, nrow_global=b_nrow)
856 CALL cp_fm_get_info(fm_out, ncol_global=c_ncol, nrow_global=c_nrow)
857 CALL dbcsr_get_info(matrix, nfullrows_total=a_nrow, nfullcols_total=a_ncol)
858 !WRITE(*,*) "cp_dbcsr_sm_fm_multiply: A ", a_nrow, "x", a_ncol
859 !WRITE(*,*) "cp_dbcsr_sm_fm_multiply: B ", b_nrow, "x", b_ncol
860 !WRITE(*,*) "cp_dbcsr_sm_fm_multiply: C ", c_nrow, "x", c_ncol
861
862 CALL cp_fm_get_info(fm_out, ncol_global=k_out)
863
864 CALL cp_fm_get_info(fm_in, ncol_global=k_in)
865 !write(*,*)routineN//" -----------------------------------"
866 !IF (k_in /= k_out) &
867 ! WRITE(*,'(3(A,I5,1X),2(A,F5.2,1X))')&
868 ! routineN//" ncol", ncol,'k_in',k_in,'k_out',k_out,&
869 ! 'alpha',my_alpha,'beta',my_beta
870
871 IF (ncol > 0 .AND. k_out > 0 .AND. k_in > 0) THEN
872 CALL dbcsr_get_info(matrix, row_blk_size=row_blk_size, col_blk_size=col_blk_size, distribution=dist)
873 CALL dbcsr_create_dist_r_unrot(dist_right_in, dist, k_in, col_blk_size_right_in)
874
875 CALL dbcsr_create(in, "D", dist_right_in, dbcsr_type_no_symmetry, &
876 col_blk_size, col_blk_size_right_in)
877
878 CALL dbcsr_distribution_get(dist, row_dist=row_dist)
879 CALL dbcsr_distribution_get(dist_right_in, col_dist=col_dist)
880 CALL dbcsr_distribution_new(product_dist, template=dist, &
881 row_dist=row_dist, col_dist=col_dist)
882 ALLOCATE (col_blk_size_right_out(SIZE(col_blk_size_right_in)))
883 col_blk_size_right_out = col_blk_size_right_in
884 CALL match_col_sizes(col_blk_size_right_out, col_blk_size_right_in, k_out)
885
886 !if (k_in .ne. k_out) then
887 ! write(*,*)routineN//" in cs", col_blk_size_right_in
888 ! write(*,*)routineN//" out cs", col_blk_size_right_out
889 !endif
890
891 CALL dbcsr_create(out, "D", product_dist, dbcsr_type_no_symmetry, &
892 row_blk_size, col_blk_size_right_out)
893
894 CALL copy_fm_to_dbcsr(fm_in, in)
895 IF (ncol /= k_out .OR. my_beta /= 0.0_dp) THEN
896 CALL copy_fm_to_dbcsr(fm_out, out)
897 END IF
898
899 CALL timeset(routinen//'_core', timing_handle_mult)
900 CALL dbcsr_multiply("N", "N", my_alpha, matrix, in, my_beta, out, &
901 last_column=ncol)
902 CALL timestop(timing_handle_mult)
903
904 CALL copy_dbcsr_to_fm(out, fm_out)
905
906 CALL dbcsr_release(in)
907 CALL dbcsr_release(out)
908 DEALLOCATE (col_blk_size_right_in, col_blk_size_right_out)
909 CALL dbcsr_distribution_release(dist_right_in)
910 CALL dbcsr_distribution_release(product_dist)
911
912 END IF
913
914 CALL timestop(timing_handle)
915
916 END SUBROUTINE cp_dbcsr_sm_fm_multiply
917
918! **************************************************************************************************
919!> \brief ...
920!> \param sizes1 ...
921!> \param sizes2 ...
922!> \param full_num ...
923! **************************************************************************************************
924 SUBROUTINE match_col_sizes(sizes1, sizes2, full_num)
925 INTEGER, DIMENSION(:), INTENT(INOUT) :: sizes1
926 INTEGER, DIMENSION(:), INTENT(IN) :: sizes2
927 INTEGER, INTENT(IN) :: full_num
928
929 INTEGER :: left, n1, n2, p, rm, used
930
931 n1 = SIZE(sizes1)
932 n2 = SIZE(sizes2)
933 IF (n1 /= n2) THEN
934 cpabort("distributions must be equal!")
935 END IF
936 sizes1(1:n1) = sizes2(1:n1)
937 used = sum(sizes1(1:n1))
938 ! If sizes1 does not cover everything, then we increase the
939 ! size of the last block; otherwise we reduce the blocks
940 ! (from the end) until it is small enough.
941 IF (used < full_num) THEN
942 sizes1(n1) = sizes1(n1) + full_num - used
943 ELSE
944 left = used - full_num
945 p = n1
946 DO WHILE (left > 0 .AND. p > 0)
947 rm = min(left, sizes1(p))
948 sizes1(p) = sizes1(p) - rm
949 left = left - rm
950 p = p - 1
951 END DO
952 END IF
953 END SUBROUTINE match_col_sizes
954
955! **************************************************************************************************
956!> \brief performs the multiplication sparse_matrix+dense_mat*dens_mat^T
957!> if matrix_g is not explicitly given, matrix_v^T will be used
958!> this can be important to save the necessary redistribute for a
959!> different matrix_g and increase performance.
960!> \param sparse_matrix ...
961!> \param matrix_v ...
962!> \param matrix_g ...
963!> \param ncol ...
964!> \param alpha ...
965!> \param keep_sparsity Determines if the sparsity of sparse_matrix is retained
966!> by default it is TRUE
967!> \param symmetry_mode There are the following modes
968!> 1: sparse_matrix += 0.5*alpha*(v*g^T+g^T*v) (symmetric update)
969!> -1: sparse_matrix += 0.5*alpha*(v*g^T-g^T*v) (skewsymmetric update)
970!> else: sparse_matrix += alpha*v*g^T (no symmetry, default)
971!> saves some redistribution steps
972! **************************************************************************************************
973 SUBROUTINE cp_dbcsr_plus_fm_fm_t(sparse_matrix, matrix_v, matrix_g, ncol, alpha, keep_sparsity, symmetry_mode)
974 TYPE(dbcsr_type), INTENT(INOUT) :: sparse_matrix
975 TYPE(cp_fm_type), INTENT(IN) :: matrix_v
976 TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: matrix_g
977 INTEGER, INTENT(IN) :: ncol
978 REAL(kind=dp), INTENT(IN), OPTIONAL :: alpha
979 LOGICAL, INTENT(IN), OPTIONAL :: keep_sparsity
980 INTEGER, INTENT(IN), OPTIONAL :: symmetry_mode
981
982 CHARACTER(LEN=*), PARAMETER :: routinen = 'cp_dbcsr_plus_fm_fm_t'
983
984 INTEGER :: k, my_symmetry_mode, nao, npcols, &
985 timing_handle
986 INTEGER, DIMENSION(:), POINTER :: col_blk_size_left, col_dist_left, &
987 row_blk_size, row_dist
988 LOGICAL :: check_product, my_keep_sparsity
989 REAL(kind=dp) :: my_alpha, norm
990 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
991 TYPE(cp_fm_type) :: fm_matrix
992 TYPE(dbcsr_distribution_type) :: dist_left, sparse_dist
993 TYPE(dbcsr_type) :: mat_g, mat_v, sparse_matrix2, &
994 sparse_matrix3
995
996 check_product = .false.
997
998 CALL timeset(routinen, timing_handle)
999
1000 my_keep_sparsity = .true.
1001 IF (PRESENT(keep_sparsity)) my_keep_sparsity = keep_sparsity
1002
1003 my_symmetry_mode = 0
1004 IF (PRESENT(symmetry_mode)) my_symmetry_mode = symmetry_mode
1005
1006 NULLIFY (col_dist_left)
1007
1008 IF (ncol > 0) THEN
1009 IF (.NOT. dbcsr_valid_index(sparse_matrix)) THEN
1010 cpabort("sparse_matrix must pre-exist")
1011 END IF
1012 !
1013 ! Setup matrix_v
1014 CALL cp_fm_get_info(matrix_v, ncol_global=k)
1015 !WRITE(*,*)routineN//'truncated mult k, ncol',k,ncol,' PRESENT (matrix_g)',PRESENT (matrix_g)
1016 CALL dbcsr_get_info(sparse_matrix, distribution=sparse_dist)
1017 CALL dbcsr_distribution_get(sparse_dist, npcols=npcols, row_dist=row_dist)
1018 CALL create_bl_distribution(col_dist_left, col_blk_size_left, k, npcols)
1019 CALL dbcsr_distribution_new(dist_left, template=sparse_dist, &
1020 row_dist=row_dist, col_dist=col_dist_left)
1021 DEALLOCATE (col_dist_left)
1022 CALL dbcsr_get_info(sparse_matrix, row_blk_size=row_blk_size)
1023 CALL dbcsr_create(mat_v, "DBCSR matrix_v", dist_left, dbcsr_type_no_symmetry, &
1024 row_blk_size, col_blk_size_left)
1025 CALL copy_fm_to_dbcsr(matrix_v, mat_v)
1026 CALL dbcsr_verify_matrix(mat_v)
1027 !
1028 ! Setup matrix_g
1029 IF (PRESENT(matrix_g)) THEN
1030 CALL dbcsr_create(mat_g, "DBCSR matrix_g", dist_left, dbcsr_type_no_symmetry, &
1031 row_blk_size, col_blk_size_left)
1032 CALL copy_fm_to_dbcsr(matrix_g, mat_g)
1033 END IF
1034 !
1035 DEALLOCATE (col_blk_size_left)
1036 CALL dbcsr_distribution_release(dist_left)
1037 !
1038 !
1039 IF (check_product) THEN
1040 CALL cp_fm_get_info(matrix_v, nrow_global=nao)
1041 CALL cp_fm_struct_create(fm_struct_tmp, context=matrix_v%matrix_struct%context, nrow_global=nao, &
1042 ncol_global=nao, para_env=matrix_v%matrix_struct%para_env)
1043 CALL cp_fm_create(fm_matrix, fm_struct_tmp, name="fm matrix")
1044 CALL cp_fm_struct_release(fm_struct_tmp)
1045 CALL copy_dbcsr_to_fm(sparse_matrix, fm_matrix)
1046 CALL dbcsr_copy(sparse_matrix3, sparse_matrix)
1047 END IF
1048 !
1049 my_alpha = 1.0_dp
1050 IF (PRESENT(alpha)) my_alpha = alpha
1051 IF (PRESENT(matrix_g)) THEN
1052 IF (my_symmetry_mode == 1) THEN
1053 ! Symmetric mode
1054 CALL dbcsr_multiply("N", "T", 0.5_dp*my_alpha, mat_v, mat_g, &
1055 1.0_dp, sparse_matrix, &
1056 retain_sparsity=my_keep_sparsity, &
1057 last_k=ncol)
1058 CALL dbcsr_multiply("N", "T", 0.5_dp*my_alpha, mat_g, mat_v, &
1059 1.0_dp, sparse_matrix, &
1060 retain_sparsity=my_keep_sparsity, &
1061 last_k=ncol)
1062 ELSE IF (my_symmetry_mode == -1) THEN
1063 ! Skewsymmetric mode
1064 CALL dbcsr_multiply("N", "T", 0.5_dp*my_alpha, mat_v, mat_g, &
1065 1.0_dp, sparse_matrix, &
1066 retain_sparsity=my_keep_sparsity, &
1067 last_k=ncol)
1068 CALL dbcsr_multiply("N", "T", -0.5_dp*my_alpha, mat_g, mat_v, &
1069 1.0_dp, sparse_matrix, &
1070 retain_sparsity=my_keep_sparsity, &
1071 last_k=ncol)
1072 ELSE
1073 ! Normal mode
1074 CALL dbcsr_multiply("N", "T", my_alpha, mat_v, mat_g, &
1075 1.0_dp, sparse_matrix, &
1076 retain_sparsity=my_keep_sparsity, &
1077 last_k=ncol)
1078 END IF
1079 ELSE
1080 CALL dbcsr_multiply("N", "T", my_alpha, mat_v, mat_v, &
1081 1.0_dp, sparse_matrix, &
1082 retain_sparsity=my_keep_sparsity, &
1083 last_k=ncol)
1084 END IF
1085
1086 IF (check_product) THEN
1087 IF (PRESENT(matrix_g)) THEN
1088 IF (my_symmetry_mode == 1) THEN
1089 CALL cp_fm_gemm("N", "T", nao, nao, ncol, 0.5_dp*my_alpha, matrix_v, matrix_g, &
1090 1.0_dp, fm_matrix)
1091 CALL cp_fm_gemm("N", "T", nao, nao, ncol, 0.5_dp*my_alpha, matrix_g, matrix_v, &
1092 1.0_dp, fm_matrix)
1093 ELSE IF (my_symmetry_mode == -1) THEN
1094 CALL cp_fm_gemm("N", "T", nao, nao, ncol, 0.5_dp*my_alpha, matrix_v, matrix_g, &
1095 1.0_dp, fm_matrix)
1096 CALL cp_fm_gemm("N", "T", nao, nao, ncol, -0.5_dp*my_alpha, matrix_g, matrix_v, &
1097 1.0_dp, fm_matrix)
1098 ELSE
1099 CALL cp_fm_gemm("N", "T", nao, nao, ncol, my_alpha, matrix_v, matrix_g, &
1100 1.0_dp, fm_matrix)
1101 END IF
1102 ELSE
1103 CALL cp_fm_gemm("N", "T", nao, nao, ncol, my_alpha, matrix_v, matrix_v, &
1104 1.0_dp, fm_matrix)
1105 END IF
1106
1107 CALL dbcsr_copy(sparse_matrix2, sparse_matrix)
1108 CALL dbcsr_scale(sparse_matrix2, alpha_scalar=0.0_dp)
1109 CALL copy_fm_to_dbcsr(fm_matrix, sparse_matrix2, keep_sparsity=my_keep_sparsity)
1110 CALL dbcsr_add(sparse_matrix2, sparse_matrix, alpha_scalar=1.0_dp, &
1111 beta_scalar=-1.0_dp)
1112 norm = dbcsr_frobenius_norm(sparse_matrix2)
1113 WRITE (*, *) 'nao=', nao, ' k=', k, ' ncol=', ncol, ' my_alpha=', my_alpha
1114 WRITE (*, *) 'PRESENT (matrix_g)', PRESENT(matrix_g)
1115 WRITE (*, *) 'matrix_type=', dbcsr_get_matrix_type(sparse_matrix)
1116 WRITE (*, *) 'norm(sm+alpha*v*g^t - fm+alpha*v*g^t)/n=', norm/real(nao, dp)
1117 CALL dbcsr_release(sparse_matrix2)
1118 CALL dbcsr_release(sparse_matrix3)
1119 CALL cp_fm_release(fm_matrix)
1120 END IF
1121 CALL dbcsr_release(mat_v)
1122 IF (PRESENT(matrix_g)) CALL dbcsr_release(mat_g)
1123 END IF
1124 CALL timestop(timing_handle)
1125
1126 END SUBROUTINE cp_dbcsr_plus_fm_fm_t
1127
1128! **************************************************************************************************
1129!> \brief Utility function to copy a specially shaped fm to dbcsr_matrix
1130!> The result matrix will be the matrix in dbcsr format
1131!> with the row blocks sizes according to the block_sizes of the template
1132!> and the col blocks sizes evenly blocked with the internal dbcsr conversion
1133!> size (32 is the current default)
1134!> \param matrix ...
1135!> \param fm_in ...
1136!> \param template ...
1137! **************************************************************************************************
1138 SUBROUTINE cp_fm_to_dbcsr_row_template(matrix, fm_in, template)
1139 TYPE(dbcsr_type), INTENT(INOUT) :: matrix
1140 TYPE(cp_fm_type), INTENT(IN) :: fm_in
1141 TYPE(dbcsr_type), INTENT(IN) :: template
1142
1143 INTEGER :: k_in
1144 INTEGER, DIMENSION(:), POINTER :: col_blk_size_right_in, row_blk_size
1145 TYPE(dbcsr_distribution_type) :: dist_right_in, tmpl_dist
1146
1147 CALL cp_fm_get_info(fm_in, ncol_global=k_in)
1148
1149 CALL dbcsr_get_info(template, distribution=tmpl_dist)
1150 CALL dbcsr_create_dist_r_unrot(dist_right_in, tmpl_dist, k_in, col_blk_size_right_in)
1151 CALL dbcsr_get_info(template, row_blk_size=row_blk_size)
1152 CALL dbcsr_create(matrix, "D", dist_right_in, dbcsr_type_no_symmetry, &
1153 row_blk_size, col_blk_size_right_in)
1154
1155 CALL copy_fm_to_dbcsr(fm_in, matrix)
1156 DEALLOCATE (col_blk_size_right_in)
1157 CALL dbcsr_distribution_release(dist_right_in)
1158
1159 END SUBROUTINE cp_fm_to_dbcsr_row_template
1160
1161! **************************************************************************************************
1162!> \brief Utility function to create an arbitrary shaped dbcsr matrix
1163!> with the same processor grid as the template matrix
1164!> both row sizes and col sizes are evenly blocked with the internal
1165!> dbcsr_conversion size (32 is the current default)
1166!> \param matrix dbcsr matrix to be created
1167!> \param template template dbcsr matrix giving its mp_env
1168!> \param m global row size of output matrix
1169!> \param n global col size of output matrix
1170!> \param sym ...
1171! **************************************************************************************************
1172 SUBROUTINE cp_dbcsr_m_by_n_from_template(matrix, template, m, n, sym)
1173 TYPE(dbcsr_type), INTENT(INOUT) :: matrix, template
1174 INTEGER, INTENT(IN) :: m, n
1175 CHARACTER, INTENT(IN), OPTIONAL :: sym
1176
1177 CHARACTER :: mysym
1178 INTEGER :: npcols, nprows
1179 INTEGER, DIMENSION(:), POINTER :: col_blk_size, col_dist, row_blk_size, &
1180 row_dist
1181 TYPE(dbcsr_distribution_type) :: dist_m_n, tmpl_dist
1182
1183 CALL dbcsr_get_info(template, matrix_type=mysym, distribution=tmpl_dist)
1184
1185 IF (PRESENT(sym)) mysym = sym
1186
1187 NULLIFY (row_dist, col_dist)
1188 NULLIFY (row_blk_size, col_blk_size)
1189 !NULLIFY (row_cluster, col_cluster)
1190
1191 CALL dbcsr_distribution_get(tmpl_dist, nprows=nprows, npcols=npcols)
1192 CALL create_bl_distribution(row_dist, row_blk_size, m, nprows)
1193 CALL create_bl_distribution(col_dist, col_blk_size, n, npcols)
1194 CALL dbcsr_distribution_new(dist_m_n, template=tmpl_dist, &
1195 row_dist=row_dist, col_dist=col_dist, &
1196 !row_cluster=row_cluster, col_cluster=col_cluster, &
1197 reuse_arrays=.true.)
1198
1199 CALL dbcsr_create(matrix, "m_n_template", dist_m_n, mysym, &
1200 row_blk_size, col_blk_size, reuse_arrays=.true.)
1201 CALL dbcsr_distribution_release(dist_m_n)
1202
1203 END SUBROUTINE cp_dbcsr_m_by_n_from_template
1204
1205! **************************************************************************************************
1206!> \brief Utility function to create dbcsr matrix, m x n matrix (n arbitrary)
1207!> with the same processor grid and row distribution as the template matrix
1208!> col sizes are evenly blocked with the internal
1209!> dbcsr_conversion size (32 is the current default)
1210!> \param matrix dbcsr matrix to be created
1211!> \param template template dbcsr matrix giving its mp_env
1212!> \param n global col size of output matrix
1213!> \param sym ...
1214! **************************************************************************************************
1215 SUBROUTINE cp_dbcsr_m_by_n_from_row_template(matrix, template, n, sym)
1216 TYPE(dbcsr_type), INTENT(INOUT) :: matrix, template
1217 INTEGER :: n
1218 CHARACTER, OPTIONAL :: sym
1219
1220 CHARACTER :: mysym
1221 INTEGER :: npcols
1222 INTEGER, DIMENSION(:), POINTER :: col_blk_size, col_dist, row_blk_size, &
1223 row_dist
1224 TYPE(dbcsr_distribution_type) :: dist_m_n, tmpl_dist
1225
1226 mysym = dbcsr_get_matrix_type(template)
1227 IF (PRESENT(sym)) mysym = sym
1228
1229 CALL dbcsr_get_info(template, distribution=tmpl_dist)
1230 CALL dbcsr_distribution_get(tmpl_dist, &
1231 npcols=npcols, &
1232 row_dist=row_dist)
1233
1234 NULLIFY (col_dist, col_blk_size)
1235 CALL create_bl_distribution(col_dist, col_blk_size, n, npcols)
1236 CALL dbcsr_distribution_new(dist_m_n, template=tmpl_dist, &
1237 row_dist=row_dist, col_dist=col_dist)
1238
1239 CALL dbcsr_get_info(template, row_blk_size=row_blk_size)
1240 CALL dbcsr_create(matrix, "m_n_template", dist_m_n, mysym, row_blk_size, col_blk_size)
1241
1242 DEALLOCATE (col_dist, col_blk_size)
1243 CALL dbcsr_distribution_release(dist_m_n)
1244
1246
1247! **************************************************************************************************
1248!> \brief Distributes elements into blocks and into bins
1249!>
1250!> \param[out] block_distribution block distribution to bins
1251!> \param[out] block_size sizes of blocks
1252!> \param[in] nelements number of elements to bin
1253!> \param[in] nbins number of bins
1254!> \par Term clarification
1255!> An example: blocks are atom blocks and bins are process rows/columns.
1256! **************************************************************************************************
1257 SUBROUTINE create_bl_distribution(block_distribution, &
1258 block_size, nelements, nbins)
1259 INTEGER, DIMENSION(:), INTENT(OUT), POINTER :: block_distribution, block_size
1260 INTEGER, INTENT(IN) :: nelements, nbins
1261
1262 CHARACTER(len=*), PARAMETER :: routinen = 'create_bl_distribution', &
1263 routinep = modulen//':'//routinen
1264
1265 INTEGER :: bin, blk_layer, element_stack, els, &
1266 estimated_blocks, max_blocks_per_bin, &
1267 nblks, nblocks, stat
1268 INTEGER, DIMENSION(:), POINTER :: blk_dist, blk_sizes
1269
1270! ---------------------------------------------------------------------------
1271
1272 NULLIFY (block_distribution)
1273 NULLIFY (block_size)
1274 ! Define the sizes on which we build the distribution.
1275 IF (nelements > 0) THEN
1276
1277 nblocks = ceiling(real(nelements, kind=dp)/real(max_elements_per_block, kind=dp))
1278 max_blocks_per_bin = ceiling(real(nblocks, kind=dp)/real(nbins, kind=dp))
1279
1280 IF (debug_mod) THEN
1281 WRITE (*, '(1X,A,1X,A,I7,A,I7,A)') routinep, "For", nelements, &
1282 " elements and", nbins, " bins"
1283 WRITE (*, '(1X,A,1X,A,I7,A)') routinep, "There are", &
1284 max_elements_per_block, " max elements per block"
1285 WRITE (*, '(1X,A,1X,A,I7,A)') routinep, "There are", &
1286 nblocks, " blocks"
1287 WRITE (*, '(1X,A,1X,A,I7,A)') routinep, "There are", &
1288 max_blocks_per_bin, " max blocks/bin"
1289 END IF
1290
1291 estimated_blocks = max_blocks_per_bin*nbins
1292 ALLOCATE (blk_dist(estimated_blocks), stat=stat)
1293 IF (stat /= 0) THEN
1294 cpabort("blk_dist")
1295 END IF
1296 ALLOCATE (blk_sizes(estimated_blocks), stat=stat)
1297 IF (stat /= 0) THEN
1298 cpabort("blk_sizes")
1299 END IF
1300 element_stack = 0
1301 nblks = 0
1302 DO blk_layer = 1, max_blocks_per_bin
1303 DO bin = 0, nbins - 1
1304 els = min(max_elements_per_block, nelements - element_stack)
1305 IF (els > 0) THEN
1306 element_stack = element_stack + els
1307 nblks = nblks + 1
1308 blk_dist(nblks) = bin
1309 blk_sizes(nblks) = els
1310 IF (debug_mod) WRITE (*, '(1X,A,I5,A,I5,A,I5)') routinep//" Assigning", &
1311 els, " elements as block", nblks, " to bin", bin
1312 END IF
1313 END DO
1314 END DO
1315 ! Create the output arrays.
1316 IF (nblks == estimated_blocks) THEN
1317 block_distribution => blk_dist
1318 block_size => blk_sizes
1319 ELSE
1320 ALLOCATE (block_distribution(nblks), stat=stat)
1321 IF (stat /= 0) THEN
1322 cpabort("blk_dist")
1323 END IF
1324 block_distribution(:) = blk_dist(1:nblks)
1325 DEALLOCATE (blk_dist)
1326 ALLOCATE (block_size(nblks), stat=stat)
1327 IF (stat /= 0) THEN
1328 cpabort("blk_sizes")
1329 END IF
1330 block_size(:) = blk_sizes(1:nblks)
1331 DEALLOCATE (blk_sizes)
1332 END IF
1333 ELSE
1334 ALLOCATE (block_distribution(0), stat=stat)
1335 IF (stat /= 0) THEN
1336 cpabort("blk_dist")
1337 END IF
1338 ALLOCATE (block_size(0), stat=stat)
1339 IF (stat /= 0) THEN
1340 cpabort("blk_sizes")
1341 END IF
1342 END IF
13431579 FORMAT(i5, 1x, i5, 1x, i5, 1x, i5, 1x, i5, 1x, i5, 1x, i5, 1x, i5, 1x, i5, 1x, i5)
1344 IF (debug_mod) THEN
1345 WRITE (*, '(1X,A,A)') routinep//" Distribution"
1346 WRITE (*, 1579) block_distribution(:)
1347 WRITE (*, '(1X,A,A)') routinep//" Sizes"
1348 WRITE (*, 1579) block_size(:)
1349 END IF
1350 END SUBROUTINE create_bl_distribution
1351
1352! **************************************************************************************************
1353!> \brief Creates a new distribution for the right matrix in a matrix
1354!> multiplication with unrotated grid.
1355!> \param[out] dist_right new distribution for the right matrix
1356!> \param[in] dist_left the distribution of the left matrix
1357!> \param[in] ncolumns number of columns in right matrix
1358!> \param[out] right_col_blk_sizes sizes of blocks in the created column
1359!> \par The new row distribution for the right matrix is the same as the row
1360!> distribution of the left matrix, while the column distribution is
1361!> created so that it is appropriate to the parallel environment.
1362! **************************************************************************************************
1363 SUBROUTINE dbcsr_create_dist_r_unrot(dist_right, dist_left, ncolumns, &
1364 right_col_blk_sizes)
1365 TYPE(dbcsr_distribution_type), INTENT(OUT) :: dist_right
1366 TYPE(dbcsr_distribution_type), INTENT(IN) :: dist_left
1367 INTEGER, INTENT(IN) :: ncolumns
1368 INTEGER, DIMENSION(:), INTENT(OUT), POINTER :: right_col_blk_sizes
1369
1370 INTEGER :: multiplicity, ncols, nimages, npcols, &
1371 nprows
1372 INTEGER, ALLOCATABLE, DIMENSION(:) :: tmp_images
1373 INTEGER, DIMENSION(:), POINTER :: old_col_dist, right_col_dist, &
1374 right_row_dist
1375
1376 CALL dbcsr_distribution_get(dist_left, &
1377 ncols=ncols, &
1378 col_dist=old_col_dist, &
1379 nprows=nprows, &
1380 npcols=npcols)
1381
1382 ! Create the column distribution
1383 CALL create_bl_distribution(right_col_dist, right_col_blk_sizes, ncolumns, npcols)
1384 ! Create an even row distribution.
1385 ALLOCATE (right_row_dist(ncols), tmp_images(ncols))
1386 nimages = lcm(nprows, npcols)/nprows
1387 multiplicity = nprows/gcd(nprows, npcols)
1388 CALL rebin_distribution(right_row_dist, tmp_images, old_col_dist, nprows, multiplicity, nimages)
1389
1390 CALL dbcsr_distribution_new(dist_right, &
1391 template=dist_left, &
1392 row_dist=right_row_dist, &
1393 col_dist=right_col_dist, &
1394 !row_cluster=dummy,&
1395 !col_cluster=dummy,&
1396 reuse_arrays=.true.)
1397 DEALLOCATE (tmp_images)
1398 END SUBROUTINE dbcsr_create_dist_r_unrot
1399
1400! **************************************************************************************************
1401!> \brief Makes new distribution with decimation and multiplicity
1402!> \param[out] new_bins new real distribution
1403!> \param[out] images new image distribution
1404!> \param[in] source_bins Basis for the new distribution and images
1405!> \param[in] nbins number of bins in the new real distribution
1406!> \param[in] multiplicity multiplicity
1407!> \param[in] nimages number of images in the new distribution
1408!> \par Definition of multiplicity and nimages
1409!> Multiplicity and decimation (number of images) are used to
1410!> match process grid coordinates on non-square process
1411!> grids. Given source_nbins and target_nbins, their relation is
1412!> source_nbins * target_multiplicity
1413!> = target_nbins * target_nimages.
1414!> It is best when both multiplicity and nimages are small. To
1415!> get these two factors, then, one can use the following formulas:
1416!> nimages = lcm(source_nbins, target_nbins) / target_nbins
1417!> multiplicity = target_nbins / gcd(source_nbins, target_nbins)
1418!> from the target's point of view (nimages = target_nimages).
1419!> \par Mapping
1420!> The new distribution comprises of real bins and images within
1421!> bins. These can be view as target_nbins*nimages virtual
1422!> columns. These same virtual columns are also
1423!> source_nbins*multiplicity in number. Therefore these virtual
1424!> columns are mapped from source_nbins*multiplicity onto
1425!> target_bins*nimages (each target bin has nimages images):
1426!> Source 4: |1 2 3|4 5 6|7 8 9|A B C| (4*3)
1427!> Target 6: |1 2|3 4|5 6|7 8|9 A|B C| (6*2)
1428!> multiplicity=3, nimages=2, 12 virtual columns (1-C).
1429!> Source bin elements are evenly mapped into one of multiplicity
1430!> virtual columns. Other (non-even, block-size aware) mappings
1431!> could be better.
1432! **************************************************************************************************
1433 SUBROUTINE rebin_distribution(new_bins, images, source_bins, &
1434 nbins, multiplicity, nimages)
1435 INTEGER, DIMENSION(:), INTENT(OUT) :: new_bins, images
1436 INTEGER, DIMENSION(:), INTENT(IN) :: source_bins
1437 INTEGER, INTENT(IN) :: nbins, multiplicity, nimages
1438
1439 INTEGER :: bin, i, old_nbins, virtual_bin
1440 INTEGER, ALLOCATABLE, DIMENSION(:) :: bin_multiplier
1441
1442! ---------------------------------------------------------------------------
1443
1444 IF (mod(nbins*nimages, multiplicity) /= 0) THEN
1445 cpwarn("mulitplicity is not divisor of new process grid coordinate")
1446 END IF
1447 old_nbins = (nbins*nimages)/multiplicity
1448 ALLOCATE (bin_multiplier(0:old_nbins - 1))
1449 bin_multiplier(:) = 0
1450 DO i = 1, SIZE(new_bins)
1451 IF (i <= SIZE(source_bins)) THEN
1452 bin = source_bins(i)
1453 ELSE
1454 ! Fill remainder with a cyclic distribution
1455 bin = mod(i, old_nbins)
1456 END IF
1457 virtual_bin = bin*multiplicity + bin_multiplier(bin)
1458 new_bins(i) = virtual_bin/nimages
1459 images(i) = 1 + mod(virtual_bin, nimages)
1460 bin_multiplier(bin) = bin_multiplier(bin) + 1
1461 IF (bin_multiplier(bin) >= multiplicity) THEN
1462 bin_multiplier(bin) = 0
1463 END IF
1464 END DO
1465 END SUBROUTINE rebin_distribution
1466
1467! **************************************************************************************************
1468!> \brief Creates a block-cyclic compatible distribution
1469!>
1470!> All blocks in a dimension, except for possibly the last
1471!> block, have the same size.
1472!> \param[out] dist the elemental distribution
1473!> \param[in] nrows number of full rows
1474!> \param[in] ncolumns number of full columns
1475!> \param[in] nrow_block size of row blocks
1476!> \param[in] ncol_block size of column blocks
1477!> \param group_handle ...
1478!> \param pgrid ...
1479!> \param[out] row_blk_sizes row block sizes
1480!> \param[out] col_blk_sizes column block sizes
1481! **************************************************************************************************
1482 SUBROUTINE dbcsr_create_dist_block_cyclic(dist, nrows, ncolumns, &
1483 nrow_block, ncol_block, group_handle, pgrid, row_blk_sizes, col_blk_sizes)
1484 TYPE(dbcsr_distribution_type), INTENT(OUT) :: dist
1485 INTEGER, INTENT(IN) :: nrows, ncolumns, nrow_block, ncol_block, &
1486 group_handle
1487 INTEGER, DIMENSION(:, :), POINTER :: pgrid
1488 INTEGER, DIMENSION(:), INTENT(OUT), POINTER :: row_blk_sizes, col_blk_sizes
1489
1490 CHARACTER(len=*), PARAMETER :: routinen = 'dbcsr_create_dist_block_cyclic'
1491
1492 INTEGER :: nblkcols, nblkrows, npcols, nprows, &
1493 pdim, sz
1494 INTEGER, DIMENSION(:), POINTER :: cd_data, rd_data
1495
1496 ! Row sizes
1497 IF (nrow_block == 0) THEN
1498 nblkrows = 0
1499 sz = 0
1500 ELSE
1501 nblkrows = nrows/nrow_block
1502 sz = mod(nrows, nrow_block)
1503 END IF
1504 IF (sz > 0) nblkrows = nblkrows + 1
1505 ALLOCATE (row_blk_sizes(nblkrows), rd_data(nblkrows))
1506 row_blk_sizes = nrow_block
1507 IF (sz /= 0) row_blk_sizes(nblkrows) = sz
1508
1509 ! Column sizes
1510 IF (ncol_block == 0) THEN
1511 nblkcols = 0
1512 sz = 0
1513 ELSE
1514 nblkcols = ncolumns/ncol_block
1515 sz = mod(ncolumns, ncol_block)
1516 END IF
1517 IF (sz > 0) nblkcols = nblkcols + 1
1518 ALLOCATE (col_blk_sizes(nblkcols), cd_data(nblkcols))
1519 col_blk_sizes = ncol_block
1520 IF (sz /= 0) col_blk_sizes(nblkcols) = sz
1521 !
1522 IF (debug_mod) THEN
1523 WRITE (*, *) routinen//" nrows,nrow_block,nblkrows=", &
1524 nrows, nrow_block, nblkrows
1525 WRITE (*, *) routinen//" ncols,ncol_block,nblkcols=", &
1526 ncolumns, ncol_block, nblkcols
1527 END IF
1528 ! Calculate process row distribution
1529 nprows = SIZE(pgrid, 1)
1530 DO pdim = 0, min(nprows - 1, nblkrows - 1)
1531 rd_data(1 + pdim:nblkrows:nprows) = pdim
1532 END DO
1533 ! Calculate process column distribution
1534 npcols = SIZE(pgrid, 2)
1535 DO pdim = 0, min(npcols - 1, nblkcols - 1)
1536 cd_data(1 + pdim:nblkcols:npcols) = pdim
1537 END DO
1538 !
1539 IF (debug_mod) THEN
1540 WRITE (*, *) routinen//" row_dist", &
1541 rd_data
1542 WRITE (*, *) routinen//" col_dist", &
1543 cd_data
1544 END IF
1545 !
1546 CALL dbcsr_distribution_new(dist, &
1547 group=group_handle, pgrid=pgrid, &
1548 row_dist=rd_data, &
1549 col_dist=cd_data, &
1550 reuse_arrays=.true.)
1551
1552 END SUBROUTINE dbcsr_create_dist_block_cyclic
1553
1554! **************************************************************************************************
1555!> \brief Allocate and initialize a real matrix 1-dimensional set.
1556!> \param[in,out] matrix_set Set containing the DBCSR matrices
1557!> \param[in] nmatrix Size of set
1558!> \par History
1559!> 2009-08-17 Adapted from sparse_matrix_type for DBCSR
1560! **************************************************************************************************
1561 SUBROUTINE allocate_dbcsr_matrix_set_1d(matrix_set, nmatrix)
1562 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_set
1563 INTEGER, INTENT(IN) :: nmatrix
1564
1565 INTEGER :: imatrix
1566
1567 IF (ASSOCIATED(matrix_set)) CALL dbcsr_deallocate_matrix_set(matrix_set)
1568 ALLOCATE (matrix_set(nmatrix))
1569 DO imatrix = 1, nmatrix
1570 NULLIFY (matrix_set(imatrix)%matrix)
1571 END DO
1572 END SUBROUTINE allocate_dbcsr_matrix_set_1d
1573
1574! **************************************************************************************************
1575!> \brief Allocate and initialize a real matrix 2-dimensional set.
1576!> \param[in,out] matrix_set Set containing the DBCSR matrix pointer type
1577!> \param[in] nmatrix Size of set
1578!> \param mmatrix ...
1579!> \par History
1580!> 2009-08-17 Adapted from sparse_matrix_type for DBCSR
1581! **************************************************************************************************
1582 SUBROUTINE allocate_dbcsr_matrix_set_2d(matrix_set, nmatrix, mmatrix)
1583 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_set
1584 INTEGER, INTENT(IN) :: nmatrix, mmatrix
1585
1586 INTEGER :: imatrix, jmatrix
1587
1588 IF (ASSOCIATED(matrix_set)) CALL dbcsr_deallocate_matrix_set(matrix_set)
1589 ALLOCATE (matrix_set(nmatrix, mmatrix))
1590 DO jmatrix = 1, mmatrix
1591 DO imatrix = 1, nmatrix
1592 NULLIFY (matrix_set(imatrix, jmatrix)%matrix)
1593 END DO
1594 END DO
1595 END SUBROUTINE allocate_dbcsr_matrix_set_2d
1596
1597! **************************************************************************************************
1598!> \brief Allocate and initialize a real matrix 3-dimensional set.
1599!> \param[in,out] matrix_set Set containing the DBCSR matrix pointer type
1600!> \param[in] nmatrix Size of set
1601!> \param mmatrix ...
1602!> \param pmatrix ...
1603!> \par History
1604!> 2009-08-17 Adapted from sparse_matrix_type for DBCSR
1605! **************************************************************************************************
1606 SUBROUTINE allocate_dbcsr_matrix_set_3d(matrix_set, nmatrix, mmatrix, pmatrix)
1607 TYPE(dbcsr_p_type), DIMENSION(:, :, :), POINTER :: matrix_set
1608 INTEGER, INTENT(IN) :: nmatrix, mmatrix, pmatrix
1609
1610 INTEGER :: imatrix, jmatrix, kmatrix
1611
1612 IF (ASSOCIATED(matrix_set)) CALL dbcsr_deallocate_matrix_set(matrix_set)
1613 ALLOCATE (matrix_set(nmatrix, mmatrix, pmatrix))
1614 DO kmatrix = 1, pmatrix
1615 DO jmatrix = 1, mmatrix
1616 DO imatrix = 1, nmatrix
1617 NULLIFY (matrix_set(imatrix, jmatrix, kmatrix)%matrix)
1618 END DO
1619 END DO
1620 END DO
1621 END SUBROUTINE allocate_dbcsr_matrix_set_3d
1622
1623! **************************************************************************************************
1624!> \brief Allocate and initialize a real matrix 4-dimensional set.
1625!> \param[in,out] matrix_set Set containing the DBCSR matrix pointer type
1626!> \param[in] nmatrix Size of set
1627!> \param mmatrix ...
1628!> \param pmatrix ...
1629!> \param qmatrix ...
1630!> \par History
1631!> 2009-08-17 Adapted from sparse_matrix_type for DBCSR
1632! **************************************************************************************************
1633 SUBROUTINE allocate_dbcsr_matrix_set_4d(matrix_set, nmatrix, mmatrix, pmatrix, qmatrix)
1634 TYPE(dbcsr_p_type), DIMENSION(:, :, :, :), POINTER :: matrix_set
1635 INTEGER, INTENT(IN) :: nmatrix, mmatrix, pmatrix, qmatrix
1636
1637 INTEGER :: imatrix, jmatrix, kmatrix, lmatrix
1638
1639 IF (ASSOCIATED(matrix_set)) CALL dbcsr_deallocate_matrix_set(matrix_set)
1640 ALLOCATE (matrix_set(nmatrix, mmatrix, pmatrix, qmatrix))
1641 DO lmatrix = 1, qmatrix
1642 DO kmatrix = 1, pmatrix
1643 DO jmatrix = 1, mmatrix
1644 DO imatrix = 1, nmatrix
1645 NULLIFY (matrix_set(imatrix, jmatrix, kmatrix, lmatrix)%matrix)
1646 END DO
1647 END DO
1648 END DO
1649 END DO
1650 END SUBROUTINE allocate_dbcsr_matrix_set_4d
1651
1652! **************************************************************************************************
1653!> \brief Allocate and initialize a real matrix 5-dimensional set.
1654!> \param[in,out] matrix_set Set containing the DBCSR matrix pointer type
1655!> \param[in] nmatrix Size of set
1656!> \param mmatrix ...
1657!> \param pmatrix ...
1658!> \param qmatrix ...
1659!> \param smatrix ...
1660!> \par History
1661!> 2009-08-17 Adapted from sparse_matrix_type for DBCSR
1662! **************************************************************************************************
1663 SUBROUTINE allocate_dbcsr_matrix_set_5d(matrix_set, nmatrix, mmatrix, pmatrix, qmatrix, smatrix)
1664 TYPE(dbcsr_p_type), DIMENSION(:, :, :, :, :), &
1665 POINTER :: matrix_set
1666 INTEGER, INTENT(IN) :: nmatrix, mmatrix, pmatrix, qmatrix, &
1667 smatrix
1668
1669 INTEGER :: hmatrix, imatrix, jmatrix, kmatrix, &
1670 lmatrix
1671
1672 IF (ASSOCIATED(matrix_set)) CALL dbcsr_deallocate_matrix_set(matrix_set)
1673 ALLOCATE (matrix_set(nmatrix, mmatrix, pmatrix, qmatrix, smatrix))
1674 DO hmatrix = 1, smatrix
1675 DO lmatrix = 1, qmatrix
1676 DO kmatrix = 1, pmatrix
1677 DO jmatrix = 1, mmatrix
1678 DO imatrix = 1, nmatrix
1679 NULLIFY (matrix_set(imatrix, jmatrix, kmatrix, lmatrix, hmatrix)%matrix)
1680 END DO
1681 END DO
1682 END DO
1683 END DO
1684 END DO
1685 END SUBROUTINE allocate_dbcsr_matrix_set_5d
1686
1687 ! **************************************************************************************************
1688!> \brief Deallocate a real matrix set and release all of the member matrices.
1689!> \param[in,out] matrix_set Set containing the DBCSR matrix pointer type
1690!> \par History
1691!> 2009-08-17 Adapted from sparse_matrix_type for DBCSR
1692! **************************************************************************************************
1693 SUBROUTINE deallocate_dbcsr_matrix_set_1d(matrix_set)
1694
1695 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_set
1696
1697 INTEGER :: imatrix
1698
1699 IF (ASSOCIATED(matrix_set)) THEN
1700 DO imatrix = lbound(matrix_set, 1), ubound(matrix_set, 1)
1701 CALL dbcsr_deallocate_matrix(matrix_set(imatrix)%matrix)
1702 END DO
1703 DEALLOCATE (matrix_set)
1704 END IF
1705
1706 END SUBROUTINE deallocate_dbcsr_matrix_set_1d
1707
1708! **************************************************************************************************
1709!> \brief Deallocate a real matrix set and release all of the member matrices.
1710!> \param[in,out] matrix_set Set containing the DBCSR matrix pointer type
1711!> \par History
1712!> 2009-08-17 Adapted from sparse_matrix_type for DBCSR
1713! **************************************************************************************************
1714 SUBROUTINE deallocate_dbcsr_matrix_set_2d(matrix_set)
1715
1716 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_set
1717
1718 INTEGER :: imatrix, jmatrix
1719
1720 IF (ASSOCIATED(matrix_set)) THEN
1721 DO jmatrix = lbound(matrix_set, 2), ubound(matrix_set, 2)
1722 DO imatrix = lbound(matrix_set, 1), ubound(matrix_set, 1)
1723 CALL dbcsr_deallocate_matrix(matrix_set(imatrix, jmatrix)%matrix)
1724 END DO
1725 END DO
1726 DEALLOCATE (matrix_set)
1727 END IF
1728 END SUBROUTINE deallocate_dbcsr_matrix_set_2d
1729
1730! **************************************************************************************************
1731!> \brief Deallocate a real matrix set and release all of the member matrices.
1732!> \param[in,out] matrix_set Set containing the DBCSR matrix pointer type
1733!> \par History
1734!> 2009-08-17 Adapted from sparse_matrix_type for DBCSR
1735! **************************************************************************************************
1736 SUBROUTINE deallocate_dbcsr_matrix_set_3d(matrix_set)
1737
1738 TYPE(dbcsr_p_type), DIMENSION(:, :, :), POINTER :: matrix_set
1739
1740 INTEGER :: imatrix, jmatrix, kmatrix
1741
1742 IF (ASSOCIATED(matrix_set)) THEN
1743 DO kmatrix = lbound(matrix_set, 3), ubound(matrix_set, 3)
1744 DO jmatrix = lbound(matrix_set, 2), ubound(matrix_set, 2)
1745 DO imatrix = lbound(matrix_set, 1), ubound(matrix_set, 1)
1746 CALL dbcsr_deallocate_matrix(matrix_set(imatrix, jmatrix, kmatrix)%matrix)
1747 END DO
1748 END DO
1749 END DO
1750 DEALLOCATE (matrix_set)
1751 END IF
1752 END SUBROUTINE deallocate_dbcsr_matrix_set_3d
1753
1754! **************************************************************************************************
1755!> \brief Deallocate a real matrix set and release all of the member matrices.
1756!> \param[in,out] matrix_set Set containing the DBCSR matrix pointer type
1757!> \par History
1758!> 2009-08-17 Adapted from sparse_matrix_type for DBCSR
1759! **************************************************************************************************
1760 SUBROUTINE deallocate_dbcsr_matrix_set_4d(matrix_set)
1761
1762 TYPE(dbcsr_p_type), DIMENSION(:, :, :, :), POINTER :: matrix_set
1763
1764 INTEGER :: imatrix, jmatrix, kmatrix, lmatrix
1765
1766 IF (ASSOCIATED(matrix_set)) THEN
1767 DO lmatrix = lbound(matrix_set, 4), ubound(matrix_set, 4)
1768 DO kmatrix = lbound(matrix_set, 3), ubound(matrix_set, 3)
1769 DO jmatrix = lbound(matrix_set, 2), ubound(matrix_set, 2)
1770 DO imatrix = lbound(matrix_set, 1), ubound(matrix_set, 1)
1771 CALL dbcsr_deallocate_matrix(matrix_set(imatrix, jmatrix, kmatrix, lmatrix)%matrix)
1772 END DO
1773 END DO
1774 END DO
1775 END DO
1776 DEALLOCATE (matrix_set)
1777 END IF
1778 END SUBROUTINE deallocate_dbcsr_matrix_set_4d
1779
1780! **************************************************************************************************
1781!> \brief Deallocate a real matrix set and release all of the member matrices.
1782!> \param[in,out] matrix_set Set containing the DBCSR matrix pointer type
1783!> \par History
1784!> 2009-08-17 Adapted from sparse_matrix_type for DBCSR
1785! **************************************************************************************************
1786 SUBROUTINE deallocate_dbcsr_matrix_set_5d(matrix_set)
1787
1788 TYPE(dbcsr_p_type), DIMENSION(:, :, :, :, :), &
1789 POINTER :: matrix_set
1790
1791 INTEGER :: hmatrix, imatrix, jmatrix, kmatrix, &
1792 lmatrix
1793
1794 IF (ASSOCIATED(matrix_set)) THEN
1795 DO hmatrix = lbound(matrix_set, 5), ubound(matrix_set, 5)
1796 DO lmatrix = lbound(matrix_set, 4), ubound(matrix_set, 4)
1797 DO kmatrix = lbound(matrix_set, 3), ubound(matrix_set, 3)
1798 DO jmatrix = lbound(matrix_set, 2), ubound(matrix_set, 2)
1799 DO imatrix = lbound(matrix_set, 1), ubound(matrix_set, 1)
1800 CALL dbcsr_deallocate_matrix(matrix_set(imatrix, jmatrix, kmatrix, lmatrix, hmatrix)%matrix)
1801 END DO
1802 END DO
1803 END DO
1804 END DO
1805 END DO
1806 DEALLOCATE (matrix_set)
1807 END IF
1808 END SUBROUTINE deallocate_dbcsr_matrix_set_5d
1809
1810END MODULE cp_dbcsr_operations
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
methods related to the blacs parallel environment
subroutine, public dbcsr_verify_matrix(matrix, verbosity, local)
...
subroutine, public dbcsr_distribution_release(dist)
...
subroutine, public dbcsr_scale(matrix, alpha_scalar)
...
subroutine, public dbcsr_distribution_new(dist, template, group, pgrid, row_dist, col_dist, reuse_arrays)
...
subroutine, public dbcsr_deallocate_matrix(matrix)
...
character function, public dbcsr_get_matrix_type(matrix)
...
logical function, public dbcsr_iterator_blocks_left(iterator)
...
subroutine, public dbcsr_iterator_stop(iterator)
...
logical function, public dbcsr_valid_index(matrix)
...
subroutine, public dbcsr_desymmetrize(matrix_a, matrix_b)
...
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset, transposed)
...
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
subroutine, public dbcsr_release(matrix)
...
subroutine, public dbcsr_complete_redistribute(matrix, redist)
...
subroutine, public dbcsr_iterator_readonly_start(iterator, matrix, shared, dynamic, dynamic_byrows)
Like dbcsr_iterator_start() but with matrix being INTENT(IN). When invoking this routine,...
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
subroutine, public dbcsr_distribution_get(dist, row_dist, col_dist, nrows, ncols, has_threads, group, mynode, numnodes, nprows, npcols, myprow, mypcol, pgrid, subgroups_defined, prow_group, pcol_group)
...
subroutine, public dbcsr_reserve_all_blocks(matrix)
Reserves all blocks.
real(dp) function, public dbcsr_frobenius_norm(matrix)
Compute the frobenius norm of a dbcsr matrix.
DBCSR operations in CP2K.
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
integer, save, public max_elements_per_block
subroutine, public copy_dbcsr_to_fm(matrix, fm, plan, check_plan)
Copy a DBCSR matrix to a BLACS matrix.
subroutine, public dbcsr_create_dist_r_unrot(dist_right, dist_left, ncolumns, right_col_blk_sizes)
Creates a new distribution for the right matrix in a matrix multiplication with unrotated grid.
subroutine, public dbcsr_multiply_local(matrix_a, vec_b, vec_c, ncol, alpha)
multiply a dbcsr with a replicated array c = alpha_scalar * A (dbscr) * b + c
subroutine, public cp_dbcsr_dist2d_to_dist(dist2d, dist)
Creates a DBCSR distribution from a distribution_2d.
subroutine, public copy_fm_to_dbcsr_bc(fm, bc_mat)
Copy a BLACS matrix to a dbcsr matrix with a special block-cyclic distribution, which requires no com...
subroutine, public copy_dbcsr_to_fm_bc(bc_mat, fm)
Copy a DBCSR_BLACS matrix to a BLACS matrix.
subroutine, public cp_fm_to_dbcsr_row_template(matrix, fm_in, template)
Utility function to copy a specially shaped fm to dbcsr_matrix The result matrix will be the matrix i...
subroutine, public cp_dbcsr_m_by_n_from_row_template(matrix, template, n, sym)
Utility function to create dbcsr matrix, m x n matrix (n arbitrary) with the same processor grid and ...
subroutine, public cp_dbcsr_plus_fm_fm_t(sparse_matrix, matrix_v, matrix_g, ncol, alpha, keep_sparsity, symmetry_mode)
performs the multiplication sparse_matrix+dense_mat*dens_mat^T if matrix_g is not explicitly given,...
subroutine, public cp_dbcsr_m_by_n_from_template(matrix, template, m, n, sym)
Utility function to create an arbitrary shaped dbcsr matrix with the same processor grid as the templ...
subroutine, public dbcsr_copy_columns_hack(matrix_b, matrix_a, ncol, source_start, target_start, para_env, blacs_env)
hack for dbcsr_copy_columns
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Basic linear algebra operations for full matrices.
subroutine, public cp_fm_gemm(transa, transb, m, n, k, alpha, matrix_a, matrix_b, beta, matrix_c, a_first_col, a_first_row, b_first_col, b_first_row, c_first_col, c_first_row)
computes matrix_c = beta * matrix_c + alpha * ( matrix_a ** transa ) * ( matrix_b ** transb )
represent the structure of a full matrix
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
stores a mapping of 2D info (e.g. matrix) on a 2D processor distribution (i.e. blacs grid) where cpus...
subroutine, public distribution_2d_get(distribution_2d, row_distribution, col_distribution, n_row_distribution, n_col_distribution, n_local_rows, n_local_cols, local_rows, local_cols, flat_local_rows, flat_local_cols, n_flat_local_rows, n_flat_local_cols, blacs_env)
returns various attributes about the distribution_2d
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public int_8
Definition kinds.F:54
integer, parameter, public dp
Definition kinds.F:34
integer, parameter, public default_string_length
Definition kinds.F:57
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
elemental integer function, public lcm(a, b)
computes the least common multiplier of two numbers
Definition mathlib.F:1324
elemental integer function, public gcd(a, b)
computes the greatest common divisor of two number
Definition mathlib.F:1289
Interface to the message passing library MPI.
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
keeps the information about the structure of a full matrix
represent a full matrix
distributes pairs on a 2d grid of processors
stores all the informations relevant to an mpi environment