LCOV - code coverage report
Current view: top level - src - cp_dbcsr_operations.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 79.6 % 452 360
Test Date: 2025-07-25 12:55:17 Functions: 75.0 % 28 21

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 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              : ! **************************************************************************************************
      18              : MODULE cp_dbcsr_operations
      19              :    USE cp_blacs_env,                    ONLY: cp_blacs_env_type
      20              :    USE cp_dbcsr_api,                    ONLY: &
      21              :         dbcsr_add, dbcsr_complete_redistribute, dbcsr_convert_sizes_to_offsets, dbcsr_copy, &
      22              :         dbcsr_create, dbcsr_deallocate_matrix, dbcsr_desymmetrize, dbcsr_distribution_get, &
      23              :         dbcsr_distribution_new, dbcsr_distribution_release, dbcsr_distribution_type, &
      24              :         dbcsr_get_info, dbcsr_get_matrix_type, dbcsr_iterator_blocks_left, &
      25              :         dbcsr_iterator_next_block, dbcsr_iterator_readonly_start, dbcsr_iterator_start, &
      26              :         dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, dbcsr_p_type, dbcsr_release, &
      27              :         dbcsr_scale, dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, &
      28              :         dbcsr_type_symmetric, dbcsr_valid_index, dbcsr_verify_matrix
      29              :    USE cp_dbcsr_contrib,                ONLY: dbcsr_frobenius_norm,&
      30              :                                               dbcsr_reserve_all_blocks
      31              :    USE cp_fm_basic_linalg,              ONLY: cp_fm_gemm
      32              :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      33              :                                               cp_fm_struct_release,&
      34              :                                               cp_fm_struct_type
      35              :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      36              :                                               cp_fm_get_info,&
      37              :                                               cp_fm_release,&
      38              :                                               cp_fm_to_fm,&
      39              :                                               cp_fm_type
      40              :    USE distribution_2d_types,           ONLY: distribution_2d_get,&
      41              :                                               distribution_2d_type
      42              :    USE kinds,                           ONLY: default_string_length,&
      43              :                                               dp
      44              :    USE mathlib,                         ONLY: gcd,&
      45              :                                               lcm
      46              :    USE message_passing,                 ONLY: mp_para_env_type
      47              : 
      48              : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads
      49              : #include "base/base_uses.f90"
      50              : 
      51              :    IMPLICIT NONE
      52              :    PRIVATE
      53              : 
      54              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_dbcsr_operations'
      55              :    LOGICAL, PARAMETER :: debug_mod = .FALSE.
      56              : 
      57              :    INTEGER, SAVE, PUBLIC :: max_elements_per_block = 32
      58              : 
      59              :    PUBLIC :: dbcsr_multiply_local
      60              : 
      61              :    ! CP2K API emulation
      62              :    PUBLIC :: copy_fm_to_dbcsr, copy_dbcsr_to_fm, &
      63              :              cp_dbcsr_sm_fm_multiply, cp_dbcsr_plus_fm_fm_t, &
      64              :              copy_dbcsr_to_fm_bc, copy_fm_to_dbcsr_bc, cp_fm_to_dbcsr_row_template, &
      65              :              cp_dbcsr_m_by_n_from_template, cp_dbcsr_m_by_n_from_row_template, &
      66              :              dbcsr_create_dist_r_unrot
      67              : 
      68              :    ! distribution_2d_type compatibility
      69              :    PUBLIC :: cp_dbcsr_dist2d_to_dist
      70              : 
      71              :    PUBLIC :: dbcsr_copy_columns_hack
      72              : 
      73              :    ! matrix set
      74              :    PUBLIC :: dbcsr_allocate_matrix_set
      75              :    PUBLIC :: dbcsr_deallocate_matrix_set
      76              : 
      77              :    INTERFACE dbcsr_allocate_matrix_set
      78              :       MODULE PROCEDURE allocate_dbcsr_matrix_set_1d
      79              :       MODULE PROCEDURE allocate_dbcsr_matrix_set_2d
      80              :       MODULE PROCEDURE allocate_dbcsr_matrix_set_3d
      81              :       MODULE PROCEDURE allocate_dbcsr_matrix_set_4d
      82              :       MODULE PROCEDURE allocate_dbcsr_matrix_set_5d
      83              :    END INTERFACE
      84              : 
      85              :    INTERFACE dbcsr_deallocate_matrix_set
      86              :       MODULE PROCEDURE deallocate_dbcsr_matrix_set_1d
      87              :       MODULE PROCEDURE deallocate_dbcsr_matrix_set_2d
      88              :       MODULE PROCEDURE deallocate_dbcsr_matrix_set_3d
      89              :       MODULE PROCEDURE deallocate_dbcsr_matrix_set_4d
      90              :       MODULE PROCEDURE deallocate_dbcsr_matrix_set_5d
      91              :    END INTERFACE
      92              : 
      93              : CONTAINS
      94              : 
      95              : ! **************************************************************************************************
      96              : !> \brief   Copy a BLACS matrix to a dbcsr matrix.
      97              : !>
      98              : !>          real_matrix=beta*real_matrix+alpha*fm
      99              : !>          beta defaults to 0, alpha to 1
     100              : !> \param[in] fm              full matrix
     101              : !> \param[out] matrix         DBCSR matrix
     102              : !> \param[in] keep_sparsity   (optional) retains the sparsity of the input
     103              : !>                            matrix
     104              : !> \date    2009-10-13
     105              : !> \par History
     106              : !>          2009-10-13 rewritten based on copy_dbcsr_to_fm
     107              : !> \author  Urban Borstnik
     108              : !> \version 2.0
     109              : ! **************************************************************************************************
     110      1093891 :    SUBROUTINE copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
     111              :       TYPE(cp_fm_type), INTENT(IN)                       :: fm
     112              :       TYPE(dbcsr_type), INTENT(INOUT)                    :: matrix
     113              :       LOGICAL, INTENT(IN), OPTIONAL                      :: keep_sparsity
     114              : 
     115              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'copy_fm_to_dbcsr'
     116              : 
     117              :       INTEGER                                            :: handle
     118              :       LOGICAL                                            :: my_keep_sparsity
     119              :       TYPE(dbcsr_type)                                   :: bc_mat, redist_mat
     120              : 
     121      1093891 :       CALL timeset(routineN, handle)
     122              : 
     123      1093891 :       my_keep_sparsity = .FALSE.
     124      1093891 :       IF (PRESENT(keep_sparsity)) my_keep_sparsity = keep_sparsity
     125              : 
     126      1093891 :       CALL copy_fm_to_dbcsr_bc(fm, bc_mat)
     127              : 
     128      1093891 :       IF (my_keep_sparsity) THEN
     129        90488 :          CALL dbcsr_create(redist_mat, template=matrix)
     130        90488 :          CALL dbcsr_complete_redistribute(bc_mat, redist_mat)
     131        90488 :          CALL dbcsr_copy(matrix, redist_mat, keep_sparsity=.TRUE.)
     132        90488 :          CALL dbcsr_release(redist_mat)
     133              :       ELSE
     134      1003403 :          CALL dbcsr_complete_redistribute(bc_mat, matrix)
     135              :       END IF
     136              : 
     137      1093891 :       CALL dbcsr_release(bc_mat)
     138              : 
     139      1093891 :       CALL timestop(handle)
     140      1093891 :    END SUBROUTINE copy_fm_to_dbcsr
     141              : 
     142              : ! **************************************************************************************************
     143              : !> \brief   Copy a BLACS matrix to a dbcsr matrix with a special block-cyclic distribution,
     144              : !>           which requires no complete redistribution.
     145              : !> \param fm ...
     146              : !> \param bc_mat ...
     147              : ! **************************************************************************************************
     148      1094261 :    SUBROUTINE copy_fm_to_dbcsr_bc(fm, bc_mat)
     149              :       TYPE(cp_fm_type), INTENT(IN)                       :: fm
     150              :       TYPE(dbcsr_type), INTENT(INOUT)                    :: bc_mat
     151              : 
     152              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'copy_fm_to_dbcsr_bc'
     153              : 
     154              :       INTEGER                                            :: col, handle, ncol_block, ncol_global, &
     155              :                                                             nrow_block, nrow_global, row
     156      1094261 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: first_col, first_row, last_col, last_row
     157      1094261 :       INTEGER, DIMENSION(:), POINTER                     :: col_blk_size, row_blk_size
     158      1094261 :       INTEGER, DIMENSION(:, :), POINTER                  :: pgrid
     159      1094261 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: dbcsr_block, fm_block
     160              :       TYPE(dbcsr_distribution_type)                      :: bc_dist
     161              :       TYPE(dbcsr_iterator_type)                          :: iter
     162              : 
     163      1094261 :       CALL timeset(routineN, handle)
     164              : 
     165      1094261 :       IF (fm%use_sp) CPABORT("copy_fm_to_dbcsr_bc: single precision not supported")
     166              : 
     167              :       ! Create processor grid
     168      1094261 :       pgrid => fm%matrix_struct%context%blacs2mpi
     169              : 
     170              :       ! Create a block-cyclic distribution compatible with the FM matrix.
     171      1094261 :       nrow_block = fm%matrix_struct%nrow_block
     172      1094261 :       ncol_block = fm%matrix_struct%ncol_block
     173      1094261 :       nrow_global = fm%matrix_struct%nrow_global
     174      1094261 :       ncol_global = fm%matrix_struct%ncol_global
     175      1094261 :       NULLIFY (col_blk_size, row_blk_size)
     176              :       CALL dbcsr_create_dist_block_cyclic(bc_dist, &
     177              :                                           nrows=nrow_global, ncolumns=ncol_global, & ! Actual full matrix size
     178              :                                           nrow_block=nrow_block, ncol_block=ncol_block, & ! BLACS parameters
     179              :                                           group_handle=fm%matrix_struct%para_env%get_handle(), pgrid=pgrid, &
     180      1094261 :                                           row_blk_sizes=row_blk_size, col_blk_sizes=col_blk_size) ! block-cyclic row/col sizes
     181              : 
     182              :       ! Create the block-cyclic DBCSR matrix
     183              :       CALL dbcsr_create(bc_mat, "Block-cyclic ", bc_dist, &
     184      1094261 :                         dbcsr_type_no_symmetry, row_blk_size, col_blk_size, reuse_arrays=.TRUE.)
     185      1094261 :       CALL dbcsr_distribution_release(bc_dist)
     186              : 
     187              :       ! allocate all blocks
     188      1094261 :       CALL dbcsr_reserve_all_blocks(bc_mat)
     189              : 
     190      1094261 :       CALL calculate_fm_block_ranges(bc_mat, first_row, last_row, first_col, last_col)
     191              : 
     192              :       ! Copy the FM data to the block-cyclic DBCSR matrix.  This step
     193              :       ! could be skipped with appropriate DBCSR index manipulation.
     194      1094261 :       fm_block => fm%local_data
     195              : !$OMP PARALLEL DEFAULT(NONE) PRIVATE(iter, row, col, dbcsr_block) &
     196      1094261 : !$OMP SHARED(bc_mat, last_row, first_row, last_col, first_col, fm_block)
     197              :       CALL dbcsr_iterator_start(iter, bc_mat)
     198              :       DO WHILE (dbcsr_iterator_blocks_left(iter))
     199              :          CALL dbcsr_iterator_next_block(iter, row, col, dbcsr_block)
     200              :          dbcsr_block(:, :) = fm_block(first_row(row):last_row(row), first_col(col):last_col(col))
     201              :       END DO
     202              :       CALL dbcsr_iterator_stop(iter)
     203              : !$OMP END PARALLEL
     204              : 
     205      1094261 :       CALL timestop(handle)
     206      2188522 :    END SUBROUTINE copy_fm_to_dbcsr_bc
     207              : 
     208              : ! **************************************************************************************************
     209              : !> \brief Copy a DBCSR matrix to a BLACS matrix
     210              : !> \param[in] matrix          DBCSR matrix
     211              : !> \param[out] fm             full matrix
     212              : ! **************************************************************************************************
     213      4140728 :    SUBROUTINE copy_dbcsr_to_fm(matrix, fm)
     214              :       TYPE(dbcsr_type), INTENT(IN)                       :: matrix
     215              :       TYPE(cp_fm_type), INTENT(INOUT)                    :: fm
     216              : 
     217              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'copy_dbcsr_to_fm'
     218              : 
     219              :       CHARACTER(len=default_string_length)               :: name
     220              :       INTEGER                                            :: group_handle, handle, ncol_block, &
     221              :                                                             nfullcols_total, nfullrows_total, &
     222              :                                                             nrow_block
     223      1035182 :       INTEGER, DIMENSION(:), POINTER                     :: col_blk_size, row_blk_size
     224      1035182 :       INTEGER, DIMENSION(:, :), POINTER                  :: pgrid
     225              :       TYPE(dbcsr_distribution_type)                      :: bc_dist, dist
     226              :       TYPE(dbcsr_type)                                   :: bc_mat, matrix_nosym
     227              : 
     228      1035182 :       CALL timeset(routineN, handle)
     229              : 
     230              :       ! check compatibility
     231              :       CALL dbcsr_get_info(matrix, &
     232              :                           name=name, &
     233              :                           distribution=dist, &
     234              :                           nfullrows_total=nfullrows_total, &
     235      1035182 :                           nfullcols_total=nfullcols_total)
     236              : 
     237      1035182 :       CPASSERT(fm%matrix_struct%nrow_global == nfullrows_total)
     238      1035182 :       CPASSERT(fm%matrix_struct%ncol_global == nfullcols_total)
     239              : 
     240              :       ! info about the full matrix
     241      1035182 :       nrow_block = fm%matrix_struct%nrow_block
     242      1035182 :       ncol_block = fm%matrix_struct%ncol_block
     243              : 
     244              :       ! Convert DBCSR to a block-cyclic
     245      1035182 :       NULLIFY (col_blk_size, row_blk_size)
     246      1035182 :       CALL dbcsr_distribution_get(dist, group=group_handle, pgrid=pgrid)
     247              :       CALL dbcsr_create_dist_block_cyclic(bc_dist, &
     248              :                                           nrows=nfullrows_total, ncolumns=nfullcols_total, &
     249              :                                           nrow_block=nrow_block, ncol_block=ncol_block, &
     250              :                                           group_handle=group_handle, pgrid=pgrid, &
     251      1035182 :                                           row_blk_sizes=row_blk_size, col_blk_sizes=col_blk_size)
     252              : 
     253              :       CALL dbcsr_create(bc_mat, "Block-cyclic"//name, bc_dist, &
     254      1035182 :                         dbcsr_type_no_symmetry, row_blk_size, col_blk_size, reuse_arrays=.TRUE.)
     255      1035182 :       CALL dbcsr_distribution_release(bc_dist)
     256              : 
     257      1035182 :       CALL dbcsr_create(matrix_nosym, template=matrix, matrix_type="N")
     258      1035182 :       CALL dbcsr_desymmetrize(matrix, matrix_nosym)
     259      1035182 :       CALL dbcsr_complete_redistribute(matrix_nosym, bc_mat)
     260      1035182 :       CALL dbcsr_release(matrix_nosym)
     261              : 
     262      1035182 :       CALL copy_dbcsr_to_fm_bc(bc_mat, fm)
     263              : 
     264      1035182 :       CALL dbcsr_release(bc_mat)
     265              : 
     266      1035182 :       CALL timestop(handle)
     267      1035182 :    END SUBROUTINE copy_dbcsr_to_fm
     268              : 
     269              : ! **************************************************************************************************
     270              : !> \brief Copy a DBCSR_BLACS matrix to a BLACS matrix
     271              : !> \param bc_mat DBCSR matrix
     272              : !> \param[out] fm             full matrix
     273              : ! **************************************************************************************************
     274      1035182 :    SUBROUTINE copy_dbcsr_to_fm_bc(bc_mat, fm)
     275              :       TYPE(dbcsr_type), INTENT(IN)                       :: bc_mat
     276              :       TYPE(cp_fm_type), INTENT(INOUT)                    :: fm
     277              : 
     278              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'copy_dbcsr_to_fm_bc'
     279              : 
     280              :       INTEGER                                            :: col, handle, row
     281      1035182 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: first_col, first_row, last_col, last_row
     282      1035182 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: dbcsr_block, fm_block
     283              :       TYPE(dbcsr_iterator_type)                          :: iter
     284              : 
     285      1035182 :       CALL timeset(routineN, handle)
     286              : 
     287      1035182 :       IF (fm%use_sp) CPABORT("copy_dbcsr_to_fm_bc: single precision not supported")
     288              : 
     289      1035182 :       CALL calculate_fm_block_ranges(bc_mat, first_row, last_row, first_col, last_col)
     290              : 
     291              :       ! Now copy data to the FM matrix
     292      1035182 :       fm_block => fm%local_data
     293    432002471 :       fm_block = REAL(0.0, KIND=dp)
     294              : !$OMP PARALLEL DEFAULT(NONE) PRIVATE(iter, row, col, dbcsr_block) &
     295      1035182 : !$OMP SHARED(bc_mat, last_row, first_row, last_col, first_col, fm_block)
     296              :       CALL dbcsr_iterator_readonly_start(iter, bc_mat)
     297              :       DO WHILE (dbcsr_iterator_blocks_left(iter))
     298              :          CALL dbcsr_iterator_next_block(iter, row, col, dbcsr_block)
     299              :          fm_block(first_row(row):last_row(row), first_col(col):last_col(col)) = dbcsr_block(:, :)
     300              :       END DO
     301              :       CALL dbcsr_iterator_stop(iter)
     302              : !$OMP END PARALLEL
     303              : 
     304      1035182 :       CALL timestop(handle)
     305      2070364 :    END SUBROUTINE copy_dbcsr_to_fm_bc
     306              : 
     307              : ! **************************************************************************************************
     308              : !> \brief Helper routine used to copy blocks from DBCSR into FM matrices and vice versa
     309              : !> \param bc_mat ...
     310              : !> \param first_row ...
     311              : !> \param last_row ...
     312              : !> \param first_col ...
     313              : !> \param last_col ...
     314              : !> \author Ole Schuett
     315              : ! **************************************************************************************************
     316      2129443 :    SUBROUTINE calculate_fm_block_ranges(bc_mat, first_row, last_row, first_col, last_col)
     317              :       TYPE(dbcsr_type), INTENT(IN)                       :: bc_mat
     318              :       INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT)    :: first_row, last_row, first_col, last_col
     319              : 
     320              :       INTEGER                                            :: col, nblkcols_local, nblkcols_total, &
     321              :                                                             nblkrows_local, nblkrows_total, row
     322              :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: local_col_sizes, local_row_sizes
     323      2129443 :       INTEGER, DIMENSION(:), POINTER                     :: col_blk_size, local_cols, local_rows, &
     324      2129443 :                                                             row_blk_size
     325              : 
     326              :       CALL dbcsr_get_info(bc_mat, &
     327              :                           nblkrows_total=nblkrows_total, &
     328              :                           nblkcols_total=nblkcols_total, &
     329              :                           nblkrows_local=nblkrows_local, &
     330              :                           nblkcols_local=nblkcols_local, &
     331              :                           local_rows=local_rows, &
     332              :                           local_cols=local_cols, &
     333              :                           row_blk_size=row_blk_size, &
     334      2129443 :                           col_blk_size=col_blk_size)
     335              : 
     336              :       ! calculate first_row and last_row
     337      6387785 :       ALLOCATE (local_row_sizes(nblkrows_total))
     338      8435421 :       local_row_sizes(:) = 0
     339      2129443 :       IF (nblkrows_local .GE. 1) THEN
     340      5318095 :          DO row = 1, nblkrows_local
     341      5318095 :             local_row_sizes(local_rows(row)) = row_blk_size(local_rows(row))
     342              :          END DO
     343              :       END IF
     344     10645039 :       ALLOCATE (first_row(nblkrows_total), last_row(nblkrows_total))
     345      2129443 :       CALL dbcsr_convert_sizes_to_offsets(local_row_sizes, first_row, last_row)
     346      2129443 :       DEALLOCATE (local_row_sizes)
     347              : 
     348              :       ! calculate first_col and last_col
     349      6386629 :       ALLOCATE (local_col_sizes(nblkcols_total))
     350      6376771 :       local_col_sizes(:) = 0
     351      2129443 :       IF (nblkcols_local .GE. 1) THEN
     352      6375071 :          DO col = 1, nblkcols_local
     353      6375071 :             local_col_sizes(local_cols(col)) = col_blk_size(local_cols(col))
     354              :          END DO
     355              :       END IF
     356     10640415 :       ALLOCATE (first_col(nblkcols_total), last_col(nblkcols_total))
     357      2129443 :       CALL dbcsr_convert_sizes_to_offsets(local_col_sizes, first_col, last_col)
     358      2129443 :       DEALLOCATE (local_col_sizes)
     359              : 
     360      2129443 :    END SUBROUTINE calculate_fm_block_ranges
     361              : 
     362              : ! **************************************************************************************************
     363              : !> \brief hack for dbcsr_copy_columns
     364              : !> \param matrix_b ...
     365              : !> \param matrix_a ...
     366              : !> \param ncol ...
     367              : !> \param source_start ...
     368              : !> \param target_start ...
     369              : !> \param para_env ...
     370              : !> \param blacs_env ...
     371              : !> \author vw
     372              : ! **************************************************************************************************
     373         9480 :    SUBROUTINE dbcsr_copy_columns_hack(matrix_b, matrix_a, &
     374              :                                       ncol, source_start, target_start, para_env, blacs_env)
     375              : 
     376              :       TYPE(dbcsr_type), INTENT(INOUT)                    :: matrix_b
     377              :       TYPE(dbcsr_type), INTENT(IN)                       :: matrix_a
     378              :       INTEGER, INTENT(IN)                                :: ncol, source_start, target_start
     379              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     380              :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
     381              : 
     382              :       INTEGER                                            :: nfullcols_total, nfullrows_total
     383              :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct
     384              :       TYPE(cp_fm_type)                                   :: fm_matrix_a, fm_matrix_b
     385              : 
     386         2370 :       NULLIFY (fm_struct)
     387         2370 :       CALL dbcsr_get_info(matrix_a, nfullrows_total=nfullrows_total, nfullcols_total=nfullcols_total)
     388              :       CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=nfullrows_total, &
     389         2370 :                                ncol_global=nfullcols_total, para_env=para_env)
     390         2370 :       CALL cp_fm_create(fm_matrix_a, fm_struct, name="fm_matrix_a")
     391         2370 :       CALL cp_fm_struct_release(fm_struct)
     392              : 
     393         2370 :       CALL dbcsr_get_info(matrix_b, nfullrows_total=nfullrows_total, nfullcols_total=nfullcols_total)
     394              :       CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=nfullrows_total, &
     395         2370 :                                ncol_global=nfullcols_total, para_env=para_env)
     396         2370 :       CALL cp_fm_create(fm_matrix_b, fm_struct, name="fm_matrix_b")
     397         2370 :       CALL cp_fm_struct_release(fm_struct)
     398              : 
     399         2370 :       CALL copy_dbcsr_to_fm(matrix_a, fm_matrix_a)
     400         2370 :       CALL copy_dbcsr_to_fm(matrix_b, fm_matrix_b)
     401              : 
     402         2370 :       CALL cp_fm_to_fm(fm_matrix_a, fm_matrix_b, ncol, source_start, target_start)
     403              : 
     404         2370 :       CALL copy_fm_to_dbcsr(fm_matrix_b, matrix_b)
     405              : 
     406         2370 :       CALL cp_fm_release(fm_matrix_a)
     407         2370 :       CALL cp_fm_release(fm_matrix_b)
     408              : 
     409         2370 :    END SUBROUTINE dbcsr_copy_columns_hack
     410              : 
     411              : ! **************************************************************************************************
     412              : !> \brief Creates a DBCSR distribution from a distribution_2d
     413              : !> \param[in] dist2d          distribution_2d
     414              : !> \param[out] dist           DBCSR distribution
     415              : !> \par History
     416              : !>    move form dbcsr_operation 01.2010
     417              : ! **************************************************************************************************
     418         9550 :    SUBROUTINE cp_dbcsr_dist2d_to_dist(dist2d, dist)
     419              :       TYPE(distribution_2d_type), INTENT(IN), TARGET     :: dist2d
     420              :       TYPE(dbcsr_distribution_type), INTENT(OUT)         :: dist
     421              : 
     422         9550 :       INTEGER, DIMENSION(:), POINTER                     :: col_dist, row_dist
     423         9550 :       INTEGER, DIMENSION(:, :), POINTER                  :: col_dist_data, pgrid, row_dist_data
     424              :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
     425              :       TYPE(distribution_2d_type), POINTER                :: dist2d_p
     426              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     427              : 
     428         9550 :       dist2d_p => dist2d
     429              :       CALL distribution_2d_get(dist2d_p, &
     430              :                                row_distribution=row_dist_data, &
     431              :                                col_distribution=col_dist_data, &
     432         9550 :                                blacs_env=blacs_env)
     433         9550 :       CALL blacs_env%get(para_env=para_env, blacs2mpi=pgrid)
     434              : 
     435              :       ! map to 1D arrays
     436         9550 :       row_dist => row_dist_data(:, 1)
     437         9550 :       col_dist => col_dist_data(:, 1)
     438              :       !row_cluster => row_dist_data(:, 2)
     439              :       !col_cluster => col_dist_data(:, 2)
     440              : 
     441              :       CALL dbcsr_distribution_new(dist, &
     442              :                                   group=para_env%get_handle(), pgrid=pgrid, &
     443              :                                   row_dist=row_dist, &
     444         9550 :                                   col_dist=col_dist)
     445              : 
     446         9550 :    END SUBROUTINE cp_dbcsr_dist2d_to_dist
     447              : 
     448              : ! **************************************************************************************************
     449              : !> \brief multiply a dbcsr with a replicated array
     450              : !>        c = alpha_scalar * A (dbscr) * b + c
     451              : !> \param[in] matrix_a DBSCR matrxx
     452              : !> \param[in]  vec_b        vectors b
     453              : !> \param[inout] vec_c      vectors c
     454              : !> \param[in]  ncol         nbr of columns
     455              : !> \param[in]  alpha        alpha
     456              : !>
     457              : ! **************************************************************************************************
     458            0 :    SUBROUTINE dbcsr_multiply_local(matrix_a, vec_b, vec_c, ncol, alpha)
     459              :       TYPE(dbcsr_type), INTENT(IN)                       :: matrix_a
     460              :       REAL(dp), DIMENSION(:, :), INTENT(IN)              :: vec_b
     461              :       REAL(dp), DIMENSION(:, :), INTENT(INOUT)           :: vec_c
     462              :       INTEGER, INTENT(in), OPTIONAL                      :: ncol
     463              :       REAL(dp), INTENT(IN), OPTIONAL                     :: alpha
     464              : 
     465              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'dbcsr_multiply_local'
     466              : 
     467              :       INTEGER                                            :: col, coloff, my_ncol, row, rowoff, &
     468              :                                                             timing_handle
     469              :       LOGICAL                                            :: has_symm
     470              :       REAL(dp)                                           :: my_alpha, my_alpha2
     471            0 :       REAL(dp), DIMENSION(:, :), POINTER                 :: data_d
     472              :       TYPE(dbcsr_iterator_type)                          :: iter
     473              : 
     474            0 :       CALL timeset(routineN, timing_handle)
     475              : 
     476            0 :       my_alpha = 1.0_dp
     477            0 :       IF (PRESENT(alpha)) my_alpha = alpha
     478              : 
     479            0 :       my_ncol = SIZE(vec_b, 2)
     480            0 :       IF (PRESENT(ncol)) my_ncol = ncol
     481              : 
     482            0 :       my_alpha2 = 0.0_dp
     483            0 :       IF (dbcsr_get_matrix_type(matrix_a) .EQ. dbcsr_type_symmetric) my_alpha2 = my_alpha
     484            0 :       IF (dbcsr_get_matrix_type(matrix_a) .EQ. dbcsr_type_antisymmetric) my_alpha2 = -my_alpha
     485              : 
     486              :       has_symm = (dbcsr_get_matrix_type(matrix_a) .EQ. dbcsr_type_symmetric .OR. &
     487            0 :                   dbcsr_get_matrix_type(matrix_a) .EQ. dbcsr_type_antisymmetric)
     488              : 
     489              : !$OMP     PARALLEL DEFAULT(NONE) SHARED(matrix_a,vec_b,vec_c,ncol,my_alpha2,my_alpha,my_ncol,has_symm) &
     490            0 : !$OMP              PRIVATE(iter,row,col,data_d,rowoff,coloff)
     491              :       CALL dbcsr_iterator_readonly_start(iter, matrix_a, dynamic=.TRUE., dynamic_byrows=.TRUE.)
     492              :       DO WHILE (dbcsr_iterator_blocks_left(iter))
     493              :          CALL dbcsr_iterator_next_block(iter, row, col, data_d, row_offset=rowoff, col_offset=coloff)
     494              :          IF (my_ncol .NE. 1) THEN
     495              :             CALL dgemm('N', 'N', &
     496              :                        SIZE(data_d, 1), my_ncol, SIZE(data_d, 2), &
     497              :                        my_alpha, data_d(1, 1), SIZE(data_d, 1), &
     498              :                        vec_b(coloff, 1), SIZE(vec_b, 1), &
     499              :                        1.0_dp, vec_c(rowoff, 1), SIZE(vec_c, 1))
     500              :          ELSE
     501              :             CALL dgemv('N', SIZE(data_d, 1), SIZE(data_d, 2), &
     502              :                        my_alpha, data_d(1, 1), SIZE(data_d, 1), &
     503              :                        vec_b(coloff, 1), 1, &
     504              :                        1.0_dp, vec_c(rowoff, 1), 1)
     505              :          END IF
     506              :       END DO
     507              :       CALL dbcsr_iterator_stop(iter)
     508              : !$OMP     END PARALLEL
     509              : 
     510              :       ! FIXME ... in the symmetric case, the writes to vec_c depend on the column, not the row. This makes OMP-ing more difficult
     511              :       ! needs e.g. a buffer for vec_c and a reduction of that buffer.
     512            0 :       IF (has_symm) THEN
     513            0 :          CALL dbcsr_iterator_readonly_start(iter, matrix_a)
     514            0 :          DO WHILE (dbcsr_iterator_blocks_left(iter))
     515            0 :             CALL dbcsr_iterator_next_block(iter, row, col, data_d, row_offset=rowoff, col_offset=coloff)
     516            0 :             IF (row .NE. col) THEN
     517            0 :                IF (my_ncol .NE. 1) THEN
     518              :                   CALL dgemm('T', 'N', &
     519              :                              SIZE(data_d, 2), my_ncol, SIZE(data_d, 1), &
     520              :                              my_alpha2, data_d(1, 1), SIZE(data_d, 1), &
     521              :                              vec_b(rowoff, 1), SIZE(vec_b, 1), &
     522            0 :                              1.0_dp, vec_c(coloff, 1), SIZE(vec_c, 1))
     523              :                ELSE
     524              :                   CALL dgemv('T', SIZE(data_d, 1), SIZE(data_d, 2), &
     525              :                              my_alpha2, data_d(1, 1), SIZE(data_d, 1), &
     526              :                              vec_b(rowoff, 1), 1, &
     527            0 :                              1.0_dp, vec_c(coloff, 1), 1)
     528              :                END IF
     529              :             END IF
     530              :          END DO
     531            0 :          CALL dbcsr_iterator_stop(iter)
     532              :       END IF
     533              : 
     534            0 :       CALL timestop(timing_handle)
     535            0 :    END SUBROUTINE dbcsr_multiply_local
     536              : 
     537              : ! **************************************************************************************************
     538              : !> \brief multiply a dbcsr with a fm matrix
     539              : !>
     540              : !> For backwards compatibility with BLAS XGEMM, this routine supports
     541              : !> the multiplication of matrices with incompatible dimensions.
     542              : !>
     543              : !> \param[in]  matrix         DBCSR matrix
     544              : !> \param fm_in full matrix
     545              : !> \param fm_out full matrix
     546              : !> \param[in]  ncol           nbr of columns
     547              : !> \param[in]  alpha          alpha
     548              : !> \param[in]  beta           beta
     549              : !>
     550              : ! **************************************************************************************************
     551      2709024 :    SUBROUTINE cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
     552              :       TYPE(dbcsr_type), INTENT(IN)                       :: matrix
     553              :       TYPE(cp_fm_type), INTENT(IN)                       :: fm_in
     554              :       TYPE(cp_fm_type), INTENT(INOUT)                    :: fm_out
     555              :       INTEGER, INTENT(IN)                                :: ncol
     556              :       REAL(dp), INTENT(IN), OPTIONAL                     :: alpha, beta
     557              : 
     558              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_dbcsr_sm_fm_multiply'
     559              : 
     560              :       INTEGER                                            :: a_ncol, a_nrow, b_ncol, b_nrow, c_ncol, &
     561              :                                                             c_nrow, k_in, k_out, timing_handle, &
     562              :                                                             timing_handle_mult
     563       451504 :       INTEGER, DIMENSION(:), POINTER                     :: col_blk_size, col_blk_size_right_in, &
     564       451504 :                                                             col_blk_size_right_out, col_dist, &
     565       451504 :                                                             row_blk_size, row_dist
     566              :       TYPE(dbcsr_type)                                   :: in, out
     567              :       TYPE(dbcsr_distribution_type)                      :: dist, dist_right_in, product_dist
     568              :       REAL(dp)                                           :: my_alpha, my_beta
     569              : 
     570       451504 :       CALL timeset(routineN, timing_handle)
     571              : 
     572       451504 :       my_alpha = 1.0_dp
     573       451504 :       my_beta = 0.0_dp
     574       451504 :       IF (PRESENT(alpha)) my_alpha = alpha
     575       451504 :       IF (PRESENT(beta)) my_beta = beta
     576              : 
     577              :       ! TODO
     578       451504 :       CALL cp_fm_get_info(fm_in, ncol_global=b_ncol, nrow_global=b_nrow)
     579       451504 :       CALL cp_fm_get_info(fm_out, ncol_global=c_ncol, nrow_global=c_nrow)
     580       451504 :       CALL dbcsr_get_info(matrix, nfullrows_total=a_nrow, nfullcols_total=a_ncol)
     581              :       !WRITE(*,*) "cp_dbcsr_sm_fm_multiply: A ", a_nrow, "x", a_ncol
     582              :       !WRITE(*,*) "cp_dbcsr_sm_fm_multiply: B ", b_nrow, "x", b_ncol
     583              :       !WRITE(*,*) "cp_dbcsr_sm_fm_multiply: C ", c_nrow, "x", c_ncol
     584              : 
     585       451504 :       CALL cp_fm_get_info(fm_out, ncol_global=k_out)
     586              : 
     587       451504 :       CALL cp_fm_get_info(fm_in, ncol_global=k_in)
     588              :       !write(*,*)routineN//" -----------------------------------"
     589              :       !IF (k_in .NE. k_out) &
     590              :       !   WRITE(*,'(3(A,I5,1X),2(A,F5.2,1X))')&
     591              :       !   routineN//" ncol", ncol,'k_in',k_in,'k_out',k_out,&
     592              :       !   'alpha',my_alpha,'beta',my_beta
     593              : 
     594       451504 :       IF (ncol .GT. 0 .AND. k_out .GT. 0 .AND. k_in .GT. 0) THEN
     595       450466 :          CALL dbcsr_get_info(matrix, row_blk_size=row_blk_size, col_blk_size=col_blk_size, distribution=dist)
     596       450466 :          CALL dbcsr_create_dist_r_unrot(dist_right_in, dist, k_in, col_blk_size_right_in)
     597              : 
     598              :          CALL dbcsr_create(in, "D", dist_right_in, dbcsr_type_no_symmetry, &
     599       450466 :                            col_blk_size, col_blk_size_right_in)
     600              : 
     601       450466 :          CALL dbcsr_distribution_get(dist, row_dist=row_dist)
     602       450466 :          CALL dbcsr_distribution_get(dist_right_in, col_dist=col_dist)
     603              :          CALL dbcsr_distribution_new(product_dist, template=dist, &
     604       450466 :                                      row_dist=row_dist, col_dist=col_dist)
     605      1351398 :          ALLOCATE (col_blk_size_right_out(SIZE(col_blk_size_right_in)))
     606      1830048 :          col_blk_size_right_out = col_blk_size_right_in
     607       450466 :          CALL match_col_sizes(col_blk_size_right_out, col_blk_size_right_in, k_out)
     608              : 
     609              :          !if (k_in .ne. k_out) then
     610              :          !   write(*,*)routineN//" in cs", col_blk_size_right_in
     611              :          !   write(*,*)routineN//" out cs", col_blk_size_right_out
     612              :          !endif
     613              : 
     614              :          CALL dbcsr_create(out, "D", product_dist, dbcsr_type_no_symmetry, &
     615       450466 :                            row_blk_size, col_blk_size_right_out)
     616              : 
     617       450466 :          CALL copy_fm_to_dbcsr(fm_in, in)
     618       450466 :          IF (ncol .NE. k_out .OR. my_beta .NE. 0.0_dp) &
     619        99604 :             CALL copy_fm_to_dbcsr(fm_out, out)
     620              : 
     621       450466 :          CALL timeset(routineN//'_core', timing_handle_mult)
     622              :          CALL dbcsr_multiply("N", "N", my_alpha, matrix, in, my_beta, out, &
     623       450466 :                              last_column=ncol)
     624       450466 :          CALL timestop(timing_handle_mult)
     625              : 
     626       450466 :          CALL copy_dbcsr_to_fm(out, fm_out)
     627              : 
     628       450466 :          CALL dbcsr_release(in)
     629       450466 :          CALL dbcsr_release(out)
     630       450466 :          DEALLOCATE (col_blk_size_right_in, col_blk_size_right_out)
     631       450466 :          CALL dbcsr_distribution_release(dist_right_in)
     632      1801864 :          CALL dbcsr_distribution_release(product_dist)
     633              : 
     634              :       END IF
     635              : 
     636       451504 :       CALL timestop(timing_handle)
     637              : 
     638       451504 :    END SUBROUTINE cp_dbcsr_sm_fm_multiply
     639              : 
     640              : ! **************************************************************************************************
     641              : !> \brief ...
     642              : !> \param sizes1 ...
     643              : !> \param sizes2 ...
     644              : !> \param full_num ...
     645              : ! **************************************************************************************************
     646       450466 :    SUBROUTINE match_col_sizes(sizes1, sizes2, full_num)
     647              :       INTEGER, DIMENSION(:), INTENT(INOUT)               :: sizes1
     648              :       INTEGER, DIMENSION(:), INTENT(IN)                  :: sizes2
     649              :       INTEGER, INTENT(IN)                                :: full_num
     650              : 
     651              :       INTEGER                                            :: left, n1, n2, p, rm, used
     652              : 
     653       450466 :       n1 = SIZE(sizes1)
     654       450466 :       n2 = SIZE(sizes2)
     655       450466 :       IF (n1 .NE. n2) &
     656            0 :          CPABORT("distributions must be equal!")
     657       915024 :       sizes1(1:n1) = sizes2(1:n1)
     658       915024 :       used = SUM(sizes1(1:n1))
     659              :       ! If sizes1 does not cover everything, then we increase the
     660              :       ! size of the last block; otherwise we reduce the blocks
     661              :       ! (from the end) until it is small enough.
     662       450466 :       IF (used .LT. full_num) THEN
     663            0 :          sizes1(n1) = sizes1(n1) + full_num - used
     664              :       ELSE
     665       450466 :          left = used - full_num
     666       450466 :          p = n1
     667       450466 :          DO WHILE (left .GT. 0 .AND. p .GT. 0)
     668            0 :             rm = MIN(left, sizes1(p))
     669            0 :             sizes1(p) = sizes1(p) - rm
     670            0 :             left = left - rm
     671            0 :             p = p - 1
     672              :          END DO
     673              :       END IF
     674       450466 :    END SUBROUTINE match_col_sizes
     675              : 
     676              : ! **************************************************************************************************
     677              : !> \brief performs the multiplication sparse_matrix+dense_mat*dens_mat^T
     678              : !>        if matrix_g is not explicitly given, matrix_v^T will be used
     679              : !>        this can be important to save the necessary redistribute for a
     680              : !>        different matrix_g and increase performance.
     681              : !> \param sparse_matrix ...
     682              : !> \param matrix_v ...
     683              : !> \param matrix_g ...
     684              : !> \param ncol ...
     685              : !> \param alpha ...
     686              : !> \param keep_sparsity Determines if the sparsity of sparse_matrix is retained
     687              : !>        by default it is TRUE
     688              : !> \param symmetry_mode There are the following modes
     689              : !>        1:     sparse_matrix += 0.5*alpha*(v*g^T+g^T*v)    (symmetric update)
     690              : !>        -1:    sparse_matrix += 0.5*alpha*(v*g^T-g^T*v)    (skewsymmetric update)
     691              : !>        else:  sparse_matrix += alpha*v*g^T                (no symmetry, default)
     692              : !>        saves some redistribution steps
     693              : ! **************************************************************************************************
     694       203770 :    SUBROUTINE cp_dbcsr_plus_fm_fm_t(sparse_matrix, matrix_v, matrix_g, ncol, alpha, keep_sparsity, symmetry_mode)
     695              :       TYPE(dbcsr_type), INTENT(INOUT)                    :: sparse_matrix
     696              :       TYPE(cp_fm_type), INTENT(IN)                       :: matrix_v
     697              :       TYPE(cp_fm_type), INTENT(IN), OPTIONAL             :: matrix_g
     698              :       INTEGER, INTENT(IN)                                :: ncol
     699              :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: alpha
     700              :       LOGICAL, INTENT(IN), OPTIONAL                      :: keep_sparsity
     701              :       INTEGER, INTENT(IN), OPTIONAL                      :: symmetry_mode
     702              : 
     703              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_dbcsr_plus_fm_fm_t'
     704              : 
     705              :       INTEGER                                            :: k, my_symmetry_mode, nao, npcols, &
     706              :                                                             timing_handle
     707       203770 :       INTEGER, DIMENSION(:), POINTER                     :: col_blk_size_left, col_dist_left, &
     708       203770 :                                                             row_blk_size, row_dist
     709              :       LOGICAL                                            :: check_product, my_keep_sparsity
     710              :       REAL(KIND=dp)                                      :: my_alpha, norm
     711              :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct_tmp
     712              :       TYPE(cp_fm_type)                                   :: fm_matrix
     713              :       TYPE(dbcsr_distribution_type)                      :: dist_left, sparse_dist
     714              :       TYPE(dbcsr_type)                                   :: mat_g, mat_v, sparse_matrix2, &
     715              :                                                             sparse_matrix3
     716              : 
     717       203770 :       check_product = .FALSE.
     718              : 
     719       203770 :       CALL timeset(routineN, timing_handle)
     720              : 
     721       203770 :       my_keep_sparsity = .TRUE.
     722       203770 :       IF (PRESENT(keep_sparsity)) my_keep_sparsity = keep_sparsity
     723              : 
     724       203770 :       my_symmetry_mode = 0
     725       203770 :       IF (PRESENT(symmetry_mode)) my_symmetry_mode = symmetry_mode
     726              : 
     727       203770 :       NULLIFY (col_dist_left)
     728              : 
     729       203770 :       IF (ncol .GT. 0) THEN
     730       202326 :          IF (.NOT. dbcsr_valid_index(sparse_matrix)) &
     731            0 :             CPABORT("sparse_matrix must pre-exist")
     732              :          !
     733              :          ! Setup matrix_v
     734       202326 :          CALL cp_fm_get_info(matrix_v, ncol_global=k)
     735              :          !WRITE(*,*)routineN//'truncated mult k, ncol',k,ncol,' PRESENT (matrix_g)',PRESENT (matrix_g)
     736       202326 :          CALL dbcsr_get_info(sparse_matrix, distribution=sparse_dist)
     737       202326 :          CALL dbcsr_distribution_get(sparse_dist, npcols=npcols, row_dist=row_dist)
     738       202326 :          CALL create_bl_distribution(col_dist_left, col_blk_size_left, k, npcols)
     739              :          CALL dbcsr_distribution_new(dist_left, template=sparse_dist, &
     740       202326 :                                      row_dist=row_dist, col_dist=col_dist_left)
     741       202326 :          DEALLOCATE (col_dist_left)
     742       202326 :          CALL dbcsr_get_info(sparse_matrix, row_blk_size=row_blk_size)
     743              :          CALL dbcsr_create(mat_v, "DBCSR matrix_v", dist_left, dbcsr_type_no_symmetry, &
     744       202326 :                            row_blk_size, col_blk_size_left)
     745       202326 :          CALL copy_fm_to_dbcsr(matrix_v, mat_v)
     746       202326 :          CALL dbcsr_verify_matrix(mat_v)
     747              :          !
     748              :          ! Setup matrix_g
     749       202326 :          IF (PRESENT(matrix_g)) THEN
     750              :             CALL dbcsr_create(mat_g, "DBCSR matrix_g", dist_left, dbcsr_type_no_symmetry, &
     751        89863 :                               row_blk_size, col_blk_size_left)
     752        89863 :             CALL copy_fm_to_dbcsr(matrix_g, mat_g)
     753              :          END IF
     754              :          !
     755       202326 :          DEALLOCATE (col_blk_size_left)
     756       202326 :          CALL dbcsr_distribution_release(dist_left)
     757              :          !
     758              :          !
     759              :          IF (check_product) THEN
     760              :             CALL cp_fm_get_info(matrix_v, nrow_global=nao)
     761              :             CALL cp_fm_struct_create(fm_struct_tmp, context=matrix_v%matrix_struct%context, nrow_global=nao, &
     762              :                                      ncol_global=nao, para_env=matrix_v%matrix_struct%para_env)
     763              :             CALL cp_fm_create(fm_matrix, fm_struct_tmp, name="fm matrix")
     764              :             CALL cp_fm_struct_release(fm_struct_tmp)
     765              :             CALL copy_dbcsr_to_fm(sparse_matrix, fm_matrix)
     766              :             CALL dbcsr_copy(sparse_matrix3, sparse_matrix)
     767              :          END IF
     768              :          !
     769       202326 :          my_alpha = 1.0_dp
     770       202326 :          IF (PRESENT(alpha)) my_alpha = alpha
     771       202326 :          IF (PRESENT(matrix_g)) THEN
     772        89863 :             IF (my_symmetry_mode == 1) THEN
     773              :                ! Symmetric mode
     774              :                CALL dbcsr_multiply("N", "T", 0.5_dp*my_alpha, mat_v, mat_g, &
     775              :                                    1.0_dp, sparse_matrix, &
     776              :                                    retain_sparsity=my_keep_sparsity, &
     777        36062 :                                    last_k=ncol)
     778              :                CALL dbcsr_multiply("N", "T", 0.5_dp*my_alpha, mat_g, mat_v, &
     779              :                                    1.0_dp, sparse_matrix, &
     780              :                                    retain_sparsity=my_keep_sparsity, &
     781        36062 :                                    last_k=ncol)
     782        53801 :             ELSE IF (my_symmetry_mode == -1) THEN
     783              :                ! Skewsymmetric mode
     784              :                CALL dbcsr_multiply("N", "T", 0.5_dp*my_alpha, mat_v, mat_g, &
     785              :                                    1.0_dp, sparse_matrix, &
     786              :                                    retain_sparsity=my_keep_sparsity, &
     787         2188 :                                    last_k=ncol)
     788              :                CALL dbcsr_multiply("N", "T", -0.5_dp*my_alpha, mat_g, mat_v, &
     789              :                                    1.0_dp, sparse_matrix, &
     790              :                                    retain_sparsity=my_keep_sparsity, &
     791         2188 :                                    last_k=ncol)
     792              :             ELSE
     793              :                ! Normal mode
     794              :                CALL dbcsr_multiply("N", "T", my_alpha, mat_v, mat_g, &
     795              :                                    1.0_dp, sparse_matrix, &
     796              :                                    retain_sparsity=my_keep_sparsity, &
     797        51613 :                                    last_k=ncol)
     798              :             END IF
     799              :          ELSE
     800              :             CALL dbcsr_multiply("N", "T", my_alpha, mat_v, mat_v, &
     801              :                                 1.0_dp, sparse_matrix, &
     802              :                                 retain_sparsity=my_keep_sparsity, &
     803       112463 :                                 last_k=ncol)
     804              :          END IF
     805              : 
     806              :          IF (check_product) THEN
     807              :             IF (PRESENT(matrix_g)) THEN
     808              :                IF (my_symmetry_mode == 1) THEN
     809              :                   CALL cp_fm_gemm("N", "T", nao, nao, ncol, 0.5_dp*my_alpha, matrix_v, matrix_g, &
     810              :                                   1.0_dp, fm_matrix)
     811              :                   CALL cp_fm_gemm("N", "T", nao, nao, ncol, 0.5_dp*my_alpha, matrix_g, matrix_v, &
     812              :                                   1.0_dp, fm_matrix)
     813              :                ELSE IF (my_symmetry_mode == -1) THEN
     814              :                   CALL cp_fm_gemm("N", "T", nao, nao, ncol, 0.5_dp*my_alpha, matrix_v, matrix_g, &
     815              :                                   1.0_dp, fm_matrix)
     816              :                   CALL cp_fm_gemm("N", "T", nao, nao, ncol, -0.5_dp*my_alpha, matrix_g, matrix_v, &
     817              :                                   1.0_dp, fm_matrix)
     818              :                ELSE
     819              :                   CALL cp_fm_gemm("N", "T", nao, nao, ncol, my_alpha, matrix_v, matrix_g, &
     820              :                                   1.0_dp, fm_matrix)
     821              :                END IF
     822              :             ELSE
     823              :                CALL cp_fm_gemm("N", "T", nao, nao, ncol, my_alpha, matrix_v, matrix_v, &
     824              :                                1.0_dp, fm_matrix)
     825              :             END IF
     826              : 
     827              :             CALL dbcsr_copy(sparse_matrix2, sparse_matrix)
     828              :             CALL dbcsr_scale(sparse_matrix2, alpha_scalar=0.0_dp)
     829              :             CALL copy_fm_to_dbcsr(fm_matrix, sparse_matrix2, keep_sparsity=my_keep_sparsity)
     830              :             CALL dbcsr_add(sparse_matrix2, sparse_matrix, alpha_scalar=1.0_dp, &
     831              :                            beta_scalar=-1.0_dp)
     832              :             norm = dbcsr_frobenius_norm(sparse_matrix2)
     833              :             WRITE (*, *) 'nao=', nao, ' k=', k, ' ncol=', ncol, ' my_alpha=', my_alpha
     834              :             WRITE (*, *) 'PRESENT (matrix_g)', PRESENT(matrix_g)
     835              :             WRITE (*, *) 'matrix_type=', dbcsr_get_matrix_type(sparse_matrix)
     836              :             WRITE (*, *) 'norm(sm+alpha*v*g^t - fm+alpha*v*g^t)/n=', norm/REAL(nao, dp)
     837              :             IF (norm/REAL(nao, dp) .GT. 1e-12_dp) THEN
     838              :                !WRITE(*,*) 'fm_matrix'
     839              :                !DO j=1,SIZE(fm_matrix%local_data,2)
     840              :                !   DO i=1,SIZE(fm_matrix%local_data,1)
     841              :                !      WRITE(*,'(A,I3,A,I3,A,E26.16,A)') 'a(',i,',',j,')=',fm_matrix%local_data(i,j),';'
     842              :                !   ENDDO
     843              :                !ENDDO
     844              :                !WRITE(*,*) 'mat_v'
     845              :                !CALL dbcsr_print(mat_v)
     846              :                !WRITE(*,*) 'mat_g'
     847              :                !CALL dbcsr_print(mat_g)
     848              :                !WRITE(*,*) 'sparse_matrix'
     849              :                !CALL dbcsr_print(sparse_matrix)
     850              :                !WRITE(*,*) 'sparse_matrix2 (-sm + sparse(fm))'
     851              :                !CALL dbcsr_print(sparse_matrix2)
     852              :                !WRITE(*,*) 'sparse_matrix3 (copy of sm input)'
     853              :                !CALL dbcsr_print(sparse_matrix3)
     854              :                !stop
     855              :             END IF
     856              :             CALL dbcsr_release(sparse_matrix2)
     857              :             CALL dbcsr_release(sparse_matrix3)
     858              :             CALL cp_fm_release(fm_matrix)
     859              :          END IF
     860       202326 :          CALL dbcsr_release(mat_v)
     861       202326 :          IF (PRESENT(matrix_g)) CALL dbcsr_release(mat_g)
     862              :       END IF
     863       203770 :       CALL timestop(timing_handle)
     864              : 
     865       203770 :    END SUBROUTINE cp_dbcsr_plus_fm_fm_t
     866              : 
     867              : ! **************************************************************************************************
     868              : !> \brief Utility function to copy a specially shaped fm to dbcsr_matrix
     869              : !>        The result matrix will be the matrix in dbcsr format
     870              : !>        with the row blocks sizes according to the block_sizes of the template
     871              : !>        and the col blocks sizes evenly blocked with the internal dbcsr conversion
     872              : !>        size (32 is the current default)
     873              : !> \param matrix ...
     874              : !> \param fm_in ...
     875              : !> \param template ...
     876              : ! **************************************************************************************************
     877        14439 :    SUBROUTINE cp_fm_to_dbcsr_row_template(matrix, fm_in, template)
     878              :       TYPE(dbcsr_type), INTENT(INOUT)                    :: matrix
     879              :       TYPE(cp_fm_type), INTENT(IN)                       :: fm_in
     880              :       TYPE(dbcsr_type), INTENT(IN)                       :: template
     881              : 
     882              :       INTEGER                                            :: k_in
     883         4813 :       INTEGER, DIMENSION(:), POINTER                     :: col_blk_size_right_in, row_blk_size
     884              :       TYPE(dbcsr_distribution_type)                      :: dist_right_in, tmpl_dist
     885              : 
     886         4813 :       CALL cp_fm_get_info(fm_in, ncol_global=k_in)
     887              : 
     888         4813 :       CALL dbcsr_get_info(template, distribution=tmpl_dist)
     889         4813 :       CALL dbcsr_create_dist_r_unrot(dist_right_in, tmpl_dist, k_in, col_blk_size_right_in)
     890         4813 :       CALL dbcsr_get_info(template, row_blk_size=row_blk_size)
     891              :       CALL dbcsr_create(matrix, "D", dist_right_in, dbcsr_type_no_symmetry, &
     892         4813 :                         row_blk_size, col_blk_size_right_in)
     893              : 
     894         4813 :       CALL copy_fm_to_dbcsr(fm_in, matrix)
     895         4813 :       DEALLOCATE (col_blk_size_right_in)
     896         4813 :       CALL dbcsr_distribution_release(dist_right_in)
     897              : 
     898         4813 :    END SUBROUTINE cp_fm_to_dbcsr_row_template
     899              : 
     900              : ! **************************************************************************************************
     901              : !> \brief Utility function to create an arbitrary shaped dbcsr matrix
     902              : !>        with the same processor grid as the template matrix
     903              : !>        both row sizes and col sizes are evenly blocked with the internal
     904              : !>        dbcsr_conversion size (32 is the current default)
     905              : !> \param matrix dbcsr matrix to be created
     906              : !> \param template template dbcsr matrix giving its mp_env
     907              : !> \param m global row size of output matrix
     908              : !> \param n global col size of output matrix
     909              : !> \param sym ...
     910              : ! **************************************************************************************************
     911       283184 :    SUBROUTINE cp_dbcsr_m_by_n_from_template(matrix, template, m, n, sym)
     912              :       TYPE(dbcsr_type), INTENT(INOUT)                    :: matrix, template
     913              :       INTEGER, INTENT(IN)                                :: m, n
     914              :       CHARACTER, INTENT(IN), OPTIONAL                    :: sym
     915              : 
     916              :       CHARACTER                                          :: mysym
     917              :       INTEGER                                            :: npcols, nprows
     918       141592 :       INTEGER, DIMENSION(:), POINTER                     :: col_blk_size, col_dist, row_blk_size, &
     919       141592 :                                                             row_dist
     920              :       TYPE(dbcsr_distribution_type)                      :: dist_m_n, tmpl_dist
     921              : 
     922       141592 :       CALL dbcsr_get_info(template, matrix_type=mysym, distribution=tmpl_dist)
     923              : 
     924       141592 :       IF (PRESENT(sym)) mysym = sym
     925              : 
     926       141592 :       NULLIFY (row_dist, col_dist)
     927       141592 :       NULLIFY (row_blk_size, col_blk_size)
     928              :       !NULLIFY (row_cluster, col_cluster)
     929              : 
     930       141592 :       CALL dbcsr_distribution_get(tmpl_dist, nprows=nprows, npcols=npcols)
     931       141592 :       CALL create_bl_distribution(row_dist, row_blk_size, m, nprows)
     932       141592 :       CALL create_bl_distribution(col_dist, col_blk_size, n, npcols)
     933              :       CALL dbcsr_distribution_new(dist_m_n, template=tmpl_dist, &
     934              :                                   row_dist=row_dist, col_dist=col_dist, &
     935              :                                   !row_cluster=row_cluster, col_cluster=col_cluster, &
     936       141592 :                                   reuse_arrays=.TRUE.)
     937              : 
     938              :       CALL dbcsr_create(matrix, "m_n_template", dist_m_n, mysym, &
     939       141592 :                         row_blk_size, col_blk_size, reuse_arrays=.TRUE.)
     940       141592 :       CALL dbcsr_distribution_release(dist_m_n)
     941              : 
     942       141592 :    END SUBROUTINE cp_dbcsr_m_by_n_from_template
     943              : 
     944              : ! **************************************************************************************************
     945              : !> \brief Utility function to create dbcsr matrix, m x n matrix (n arbitrary)
     946              : !>        with the same processor grid and row distribution  as the template matrix
     947              : !>        col sizes are evenly blocked with the internal
     948              : !>        dbcsr_conversion size (32 is the current default)
     949              : !> \param matrix dbcsr matrix to be created
     950              : !> \param template template dbcsr matrix giving its mp_env
     951              : !> \param n global col size of output matrix
     952              : !> \param sym ...
     953              : ! **************************************************************************************************
     954       504160 :    SUBROUTINE cp_dbcsr_m_by_n_from_row_template(matrix, template, n, sym)
     955              :       TYPE(dbcsr_type), INTENT(INOUT)                    :: matrix, template
     956              :       INTEGER                                            :: n
     957              :       CHARACTER, OPTIONAL                                :: sym
     958              : 
     959              :       CHARACTER                                          :: mysym
     960              :       INTEGER                                            :: npcols
     961       126040 :       INTEGER, DIMENSION(:), POINTER                     :: col_blk_size, col_dist, row_blk_size, &
     962       126040 :                                                             row_dist
     963              :       TYPE(dbcsr_distribution_type)                      :: dist_m_n, tmpl_dist
     964              : 
     965       252080 :       mysym = dbcsr_get_matrix_type(template)
     966       126040 :       IF (PRESENT(sym)) mysym = sym
     967              : 
     968       126040 :       CALL dbcsr_get_info(template, distribution=tmpl_dist)
     969              :       CALL dbcsr_distribution_get(tmpl_dist, &
     970              :                                   npcols=npcols, &
     971       126040 :                                   row_dist=row_dist)
     972              : 
     973       126040 :       NULLIFY (col_dist, col_blk_size)
     974       126040 :       CALL create_bl_distribution(col_dist, col_blk_size, n, npcols)
     975              :       CALL dbcsr_distribution_new(dist_m_n, template=tmpl_dist, &
     976       126040 :                                   row_dist=row_dist, col_dist=col_dist)
     977              : 
     978       126040 :       CALL dbcsr_get_info(template, row_blk_size=row_blk_size)
     979       126040 :       CALL dbcsr_create(matrix, "m_n_template", dist_m_n, mysym, row_blk_size, col_blk_size)
     980              : 
     981       126040 :       DEALLOCATE (col_dist, col_blk_size)
     982       126040 :       CALL dbcsr_distribution_release(dist_m_n)
     983              : 
     984       126040 :    END SUBROUTINE cp_dbcsr_m_by_n_from_row_template
     985              : 
     986              : ! **************************************************************************************************
     987              : !> \brief Distributes elements into blocks and into bins
     988              : !>
     989              : !> \param[out] block_distribution       block distribution to bins
     990              : !> \param[out] block_size       sizes of blocks
     991              : !> \param[in] nelements number of elements to bin
     992              : !> \param[in] nbins             number of bins
     993              : !> \par Term clarification
     994              : !>      An example: blocks are atom blocks and bins are process rows/columns.
     995              : ! **************************************************************************************************
     996      1066829 :    SUBROUTINE create_bl_distribution(block_distribution, &
     997              :                                      block_size, nelements, nbins)
     998              :       INTEGER, DIMENSION(:), INTENT(OUT), POINTER        :: block_distribution, block_size
     999              :       INTEGER, INTENT(IN)                                :: nelements, nbins
    1000              : 
    1001              :       CHARACTER(len=*), PARAMETER :: routineN = 'create_bl_distribution', &
    1002              :          routineP = moduleN//':'//routineN
    1003              : 
    1004              :       INTEGER                                            :: bin, blk_layer, element_stack, els, &
    1005              :                                                             estimated_blocks, max_blocks_per_bin, &
    1006              :                                                             nblks, nblocks, stat
    1007      1066829 :       INTEGER, DIMENSION(:), POINTER                     :: blk_dist, blk_sizes
    1008              : 
    1009              : !   ---------------------------------------------------------------------------
    1010              : 
    1011      1066829 :       NULLIFY (block_distribution)
    1012      1066829 :       NULLIFY (block_size)
    1013              :       ! Define the sizes on which we build the distribution.
    1014      1066829 :       IF (nelements .GT. 0) THEN
    1015              : 
    1016      1059649 :          nblocks = CEILING(REAL(nelements, KIND=dp)/REAL(max_elements_per_block, KIND=dp))
    1017      1059649 :          max_blocks_per_bin = CEILING(REAL(nblocks, KIND=dp)/REAL(nbins, KIND=dp))
    1018              : 
    1019              :          IF (debug_mod) THEN
    1020              :             WRITE (*, '(1X,A,1X,A,I7,A,I7,A)') routineP, "For", nelements, &
    1021              :                " elements and", nbins, " bins"
    1022              :             WRITE (*, '(1X,A,1X,A,I7,A)') routineP, "There are", &
    1023              :                max_elements_per_block, " max elements per block"
    1024              :             WRITE (*, '(1X,A,1X,A,I7,A)') routineP, "There are", &
    1025              :                nblocks, " blocks"
    1026              :             WRITE (*, '(1X,A,1X,A,I7,A)') routineP, "There are", &
    1027              :                max_blocks_per_bin, " max blocks/bin"
    1028              :          END IF
    1029              : 
    1030      1059649 :          estimated_blocks = max_blocks_per_bin*nbins
    1031      3178947 :          ALLOCATE (blk_dist(estimated_blocks), stat=stat)
    1032      1059649 :          IF (stat /= 0) &
    1033            0 :             CPABORT("blk_dist")
    1034      2119298 :          ALLOCATE (blk_sizes(estimated_blocks), stat=stat)
    1035              :          IF (stat /= 0) &
    1036            0 :             CPABORT("blk_sizes")
    1037      1059649 :          element_stack = 0
    1038      1059649 :          nblks = 0
    1039      2180964 :          DO blk_layer = 1, max_blocks_per_bin
    1040      3425329 :             DO bin = 0, nbins - 1
    1041      1244365 :                els = MIN(max_elements_per_block, nelements - element_stack)
    1042      2365680 :                IF (els .GT. 0) THEN
    1043      1132675 :                   element_stack = element_stack + els
    1044      1132675 :                   nblks = nblks + 1
    1045      1132675 :                   blk_dist(nblks) = bin
    1046      1132675 :                   blk_sizes(nblks) = els
    1047              :                   IF (debug_mod) WRITE (*, '(1X,A,I5,A,I5,A,I5)') routineP//" Assigning", &
    1048              :                      els, " elements as block", nblks, " to bin", bin
    1049              :                END IF
    1050              :             END DO
    1051              :          END DO
    1052              :          ! Create the output arrays.
    1053      1059649 :          IF (nblks .EQ. estimated_blocks) THEN
    1054       947959 :             block_distribution => blk_dist
    1055       947959 :             block_size => blk_sizes
    1056              :          ELSE
    1057       335070 :             ALLOCATE (block_distribution(nblks), stat=stat)
    1058              :             IF (stat /= 0) &
    1059            0 :                CPABORT("blk_dist")
    1060       452280 :             block_distribution(:) = blk_dist(1:nblks)
    1061       111690 :             DEALLOCATE (blk_dist)
    1062       223380 :             ALLOCATE (block_size(nblks), stat=stat)
    1063              :             IF (stat /= 0) &
    1064            0 :                CPABORT("blk_sizes")
    1065       452280 :             block_size(:) = blk_sizes(1:nblks)
    1066       111690 :             DEALLOCATE (blk_sizes)
    1067              :          END IF
    1068              :       ELSE
    1069         7180 :          ALLOCATE (block_distribution(0), stat=stat)
    1070              :          IF (stat /= 0) &
    1071            0 :             CPABORT("blk_dist")
    1072         7180 :          ALLOCATE (block_size(0), stat=stat)
    1073              :          IF (stat /= 0) &
    1074            0 :             CPABORT("blk_sizes")
    1075              :       END IF
    1076              : 1579  FORMAT(I5, 1X, I5, 1X, I5, 1X, I5, 1X, I5, 1X, I5, 1X, I5, 1X, I5, 1X, I5, 1X, I5)
    1077              :       IF (debug_mod) THEN
    1078              :          WRITE (*, '(1X,A,A)') routineP//" Distribution"
    1079              :          WRITE (*, 1579) block_distribution(:)
    1080              :          WRITE (*, '(1X,A,A)') routineP//" Sizes"
    1081              :          WRITE (*, 1579) block_size(:)
    1082              :       END IF
    1083      1066829 :    END SUBROUTINE create_bl_distribution
    1084              : 
    1085              : ! **************************************************************************************************
    1086              : !> \brief Creates a new distribution for the right matrix in a matrix
    1087              : !>        multiplication with unrotated grid.
    1088              : !> \param[out] dist_right     new distribution for the right matrix
    1089              : !> \param[in] dist_left       the distribution of the left matrix
    1090              : !> \param[in] ncolumns        number of columns in right matrix
    1091              : !> \param[out] right_col_blk_sizes      sizes of blocks in the created column
    1092              : !> \par The new row distribution for the right matrix is the same as the row
    1093              : !>      distribution of the left matrix, while the column distribution is
    1094              : !>      created so that it is appropriate to the parallel environment.
    1095              : ! **************************************************************************************************
    1096       455279 :    SUBROUTINE dbcsr_create_dist_r_unrot(dist_right, dist_left, ncolumns, &
    1097              :                                         right_col_blk_sizes)
    1098              :       TYPE(dbcsr_distribution_type), INTENT(OUT)         :: dist_right
    1099              :       TYPE(dbcsr_distribution_type), INTENT(IN)          :: dist_left
    1100              :       INTEGER, INTENT(IN)                                :: ncolumns
    1101              :       INTEGER, DIMENSION(:), INTENT(OUT), POINTER        :: right_col_blk_sizes
    1102              : 
    1103              :       INTEGER                                            :: multiplicity, ncols, nimages, npcols, &
    1104              :                                                             nprows
    1105              :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: tmp_images
    1106       455279 :       INTEGER, DIMENSION(:), POINTER                     :: old_col_dist, right_col_dist, &
    1107       455279 :                                                             right_row_dist
    1108              : 
    1109              :       CALL dbcsr_distribution_get(dist_left, &
    1110              :                                   ncols=ncols, &
    1111              :                                   col_dist=old_col_dist, &
    1112              :                                   nprows=nprows, &
    1113       455279 :                                   npcols=npcols)
    1114              : 
    1115              :       ! Create the column distribution
    1116       455279 :       CALL create_bl_distribution(right_col_dist, right_col_blk_sizes, ncolumns, npcols)
    1117              :       ! Create an even row distribution.
    1118      1821116 :       ALLOCATE (right_row_dist(ncols), tmp_images(ncols))
    1119       455279 :       nimages = lcm(nprows, npcols)/nprows
    1120       455279 :       multiplicity = nprows/gcd(nprows, npcols)
    1121       455279 :       CALL rebin_distribution(right_row_dist, tmp_images, old_col_dist, nprows, multiplicity, nimages)
    1122              : 
    1123              :       CALL dbcsr_distribution_new(dist_right, &
    1124              :                                   template=dist_left, &
    1125              :                                   row_dist=right_row_dist, &
    1126              :                                   col_dist=right_col_dist, &
    1127              :                                   !row_cluster=dummy,&
    1128              :                                   !col_cluster=dummy,&
    1129       455279 :                                   reuse_arrays=.TRUE.)
    1130       455279 :       DEALLOCATE (tmp_images)
    1131       455279 :    END SUBROUTINE dbcsr_create_dist_r_unrot
    1132              : 
    1133              : ! **************************************************************************************************
    1134              : !> \brief Makes new distribution with decimation and multiplicity
    1135              : !> \param[out] new_bins      new real distribution
    1136              : !> \param[out] images        new image distribution
    1137              : !> \param[in] source_bins    Basis for the new distribution and images
    1138              : !> \param[in] nbins          number of bins in the new real distribution
    1139              : !> \param[in] multiplicity   multiplicity
    1140              : !> \param[in] nimages        number of images in the new distribution
    1141              : !> \par Definition of multiplicity and nimages
    1142              : !>      Multiplicity and decimation (number of images) are used to
    1143              : !>      match process grid coordinates on non-square process
    1144              : !>      grids. Given source_nbins and target_nbins, their relation is
    1145              : !>                source_nbins * target_multiplicity
    1146              : !>              = target_nbins * target_nimages.
    1147              : !>      It is best when both multiplicity and nimages are small. To
    1148              : !>      get these two factors, then, one can use the following formulas:
    1149              : !>          nimages      = lcm(source_nbins, target_nbins) / target_nbins
    1150              : !>          multiplicity = target_nbins / gcd(source_nbins, target_nbins)
    1151              : !>      from the target's point of view (nimages = target_nimages).
    1152              : !> \par Mapping
    1153              : !>      The new distribution comprises of real bins and images within
    1154              : !>      bins. These can be view as target_nbins*nimages virtual
    1155              : !>      columns. These same virtual columns are also
    1156              : !>      source_nbins*multiplicity in number. Therefore these virtual
    1157              : !>      columns are mapped from source_nbins*multiplicity onto
    1158              : !>      target_bins*nimages (each target bin has nimages images):
    1159              : !>      Source 4: |1 2 3|4 5 6|7 8 9|A B C| (4*3)
    1160              : !>      Target 6: |1 2|3 4|5 6|7 8|9 A|B C| (6*2)
    1161              : !>      multiplicity=3, nimages=2, 12 virtual columns (1-C).
    1162              : !>      Source bin elements are evenly mapped into one of multiplicity
    1163              : !>      virtual columns. Other (non-even, block-size aware) mappings
    1164              : !>      could be better.
    1165              : ! **************************************************************************************************
    1166       455279 :    SUBROUTINE rebin_distribution(new_bins, images, source_bins, &
    1167              :                                  nbins, multiplicity, nimages)
    1168              :       INTEGER, DIMENSION(:), INTENT(OUT)                 :: new_bins, images
    1169              :       INTEGER, DIMENSION(:), INTENT(IN)                  :: source_bins
    1170              :       INTEGER, INTENT(IN)                                :: nbins, multiplicity, nimages
    1171              : 
    1172              :       INTEGER                                            :: bin, i, old_nbins, virtual_bin
    1173       455279 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: bin_multiplier
    1174              : 
    1175              : !   ---------------------------------------------------------------------------
    1176              : 
    1177       455279 :       IF (MOD(nbins*nimages, multiplicity) .NE. 0) THEN
    1178            0 :          CPWARN("mulitplicity is not divisor of new process grid coordinate")
    1179              :       END IF
    1180       455279 :       old_nbins = (nbins*nimages)/multiplicity
    1181      1365837 :       ALLOCATE (bin_multiplier(0:old_nbins - 1))
    1182       910558 :       bin_multiplier(:) = 0
    1183      2167131 :       DO i = 1, SIZE(new_bins)
    1184      1711852 :          IF (i .LE. SIZE(source_bins)) THEN
    1185      1711852 :             bin = source_bins(i)
    1186              :          ELSE
    1187              :             ! Fill remainder with a cyclic distribution
    1188            0 :             bin = MOD(i, old_nbins)
    1189              :          END IF
    1190      1711852 :          virtual_bin = bin*multiplicity + bin_multiplier(bin)
    1191      1711852 :          new_bins(i) = virtual_bin/nimages
    1192      1711852 :          images(i) = 1 + MOD(virtual_bin, nimages)
    1193      1711852 :          bin_multiplier(bin) = bin_multiplier(bin) + 1
    1194      2167131 :          IF (bin_multiplier(bin) .GE. multiplicity) THEN
    1195       732550 :             bin_multiplier(bin) = 0
    1196              :          END IF
    1197              :       END DO
    1198       455279 :    END SUBROUTINE rebin_distribution
    1199              : 
    1200              : ! **************************************************************************************************
    1201              : !> \brief Creates a block-cyclic compatible distribution
    1202              : !>
    1203              : !>        All blocks in a dimension, except for possibly the last
    1204              : !>        block, have the same size.
    1205              : !> \param[out] dist           the elemental distribution
    1206              : !> \param[in] nrows           number of full rows
    1207              : !> \param[in] ncolumns        number of full columns
    1208              : !> \param[in] nrow_block      size of row blocks
    1209              : !> \param[in] ncol_block      size of column blocks
    1210              : !> \param group_handle ...
    1211              : !> \param pgrid ...
    1212              : !> \param[out] row_blk_sizes  row block sizes
    1213              : !> \param[out] col_blk_sizes  column block sizes
    1214              : ! **************************************************************************************************
    1215      2129443 :    SUBROUTINE dbcsr_create_dist_block_cyclic(dist, nrows, ncolumns, &
    1216              :                                              nrow_block, ncol_block, group_handle, pgrid, row_blk_sizes, col_blk_sizes)
    1217              :       TYPE(dbcsr_distribution_type), INTENT(OUT)         :: dist
    1218              :       INTEGER, INTENT(IN)                                :: nrows, ncolumns, nrow_block, ncol_block, &
    1219              :                                                             group_handle
    1220              :       INTEGER, DIMENSION(:, :), POINTER                  :: pgrid
    1221              :       INTEGER, DIMENSION(:), INTENT(OUT), POINTER        :: row_blk_sizes, col_blk_sizes
    1222              : 
    1223              :       CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_create_dist_block_cyclic'
    1224              : 
    1225              :       INTEGER                                            :: nblkcols, nblkrows, npcols, nprows, &
    1226              :                                                             pdim, sz
    1227      2129443 :       INTEGER, DIMENSION(:), POINTER                     :: cd_data, rd_data
    1228              : 
    1229              :       ! Row sizes
    1230      2129443 :       IF (nrow_block .EQ. 0) THEN
    1231              :          nblkrows = 0
    1232              :          sz = 0
    1233              :       ELSE
    1234      2129443 :          nblkrows = nrows/nrow_block
    1235      2129443 :          sz = MOD(nrows, nrow_block)
    1236              :       END IF
    1237      2129443 :       IF (sz .GT. 0) nblkrows = nblkrows + 1
    1238     10645583 :       ALLOCATE (row_blk_sizes(nblkrows), rd_data(nblkrows))
    1239      8435421 :       row_blk_sizes = nrow_block
    1240      2129443 :       IF (sz .NE. 0) row_blk_sizes(nblkrows) = sz
    1241              : 
    1242              :       ! Column sizes
    1243      2129443 :       IF (ncol_block .EQ. 0) THEN
    1244              :          nblkcols = 0
    1245              :          sz = 0
    1246              :       ELSE
    1247      2129443 :          nblkcols = ncolumns/ncol_block
    1248      2129443 :          sz = MOD(ncolumns, ncol_block)
    1249              :       END IF
    1250      2129443 :       IF (sz .GT. 0) nblkcols = nblkcols + 1
    1251     10642115 :       ALLOCATE (col_blk_sizes(nblkcols), cd_data(nblkcols))
    1252      6376771 :       col_blk_sizes = ncol_block
    1253      2129443 :       IF (sz .NE. 0) col_blk_sizes(nblkcols) = sz
    1254              :       !
    1255              :       IF (debug_mod) THEN
    1256              :          WRITE (*, *) routineN//" nrows,nrow_block,nblkrows=", &
    1257              :             nrows, nrow_block, nblkrows
    1258              :          WRITE (*, *) routineN//" ncols,ncol_block,nblkcols=", &
    1259              :             ncolumns, ncol_block, nblkcols
    1260              :       END IF
    1261              :       ! Calculate process row distribution
    1262      2129443 :       nprows = SIZE(pgrid, 1)
    1263      6315350 :       DO pdim = 0, MIN(nprows - 1, nblkrows - 1)
    1264     12621328 :          rd_data(1 + pdim:nblkrows:nprows) = pdim
    1265              :       END DO
    1266              :       ! Calculate process column distribution
    1267      2129443 :       npcols = SIZE(pgrid, 2)
    1268      4257186 :       DO pdim = 0, MIN(npcols - 1, nblkcols - 1)
    1269      8504514 :          cd_data(1 + pdim:nblkcols:npcols) = pdim
    1270              :       END DO
    1271              :       !
    1272              :       IF (debug_mod) THEN
    1273              :          WRITE (*, *) routineN//" row_dist", &
    1274              :             rd_data
    1275              :          WRITE (*, *) routineN//" col_dist", &
    1276              :             cd_data
    1277              :       END IF
    1278              :       !
    1279              :       CALL dbcsr_distribution_new(dist, &
    1280              :                                   group=group_handle, pgrid=pgrid, &
    1281              :                                   row_dist=rd_data, &
    1282              :                                   col_dist=cd_data, &
    1283      2129443 :                                   reuse_arrays=.TRUE.)
    1284              : 
    1285      2129443 :    END SUBROUTINE dbcsr_create_dist_block_cyclic
    1286              : 
    1287              : ! **************************************************************************************************
    1288              : !> \brief   Allocate and initialize a real matrix 1-dimensional set.
    1289              : !> \param[in,out] matrix_set  Set containing the DBCSR matrices
    1290              : !> \param[in] nmatrix         Size of set
    1291              : !> \par History
    1292              : !>      2009-08-17 Adapted from sparse_matrix_type for DBCSR
    1293              : ! **************************************************************************************************
    1294       181489 :    SUBROUTINE allocate_dbcsr_matrix_set_1d(matrix_set, nmatrix)
    1295              :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_set
    1296              :       INTEGER, INTENT(IN)                                :: nmatrix
    1297              : 
    1298              :       INTEGER                                            :: imatrix
    1299              : 
    1300       181489 :       IF (ASSOCIATED(matrix_set)) CALL dbcsr_deallocate_matrix_set(matrix_set)
    1301       859400 :       ALLOCATE (matrix_set(nmatrix))
    1302       496422 :       DO imatrix = 1, nmatrix
    1303       496422 :          NULLIFY (matrix_set(imatrix)%matrix)
    1304              :       END DO
    1305       181489 :    END SUBROUTINE allocate_dbcsr_matrix_set_1d
    1306              : 
    1307              : ! **************************************************************************************************
    1308              : !> \brief   Allocate and initialize a real matrix 2-dimensional set.
    1309              : !> \param[in,out] matrix_set  Set containing the DBCSR matrix pointer type
    1310              : !> \param[in] nmatrix         Size of set
    1311              : !> \param mmatrix ...
    1312              : !> \par History
    1313              : !>      2009-08-17 Adapted from sparse_matrix_type for DBCSR
    1314              : ! **************************************************************************************************
    1315       150854 :    SUBROUTINE allocate_dbcsr_matrix_set_2d(matrix_set, nmatrix, mmatrix)
    1316              :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_set
    1317              :       INTEGER, INTENT(IN)                                :: nmatrix, mmatrix
    1318              : 
    1319              :       INTEGER                                            :: imatrix, jmatrix
    1320              : 
    1321       150854 :       IF (ASSOCIATED(matrix_set)) CALL dbcsr_deallocate_matrix_set(matrix_set)
    1322      2227498 :       ALLOCATE (matrix_set(nmatrix, mmatrix))
    1323       907044 :       DO jmatrix = 1, mmatrix
    1324      1774936 :          DO imatrix = 1, nmatrix
    1325      1624082 :             NULLIFY (matrix_set(imatrix, jmatrix)%matrix)
    1326              :          END DO
    1327              :       END DO
    1328       150854 :    END SUBROUTINE allocate_dbcsr_matrix_set_2d
    1329              : 
    1330              : ! **************************************************************************************************
    1331              : !> \brief   Allocate and initialize a real matrix 3-dimensional set.
    1332              : !> \param[in,out] matrix_set  Set containing the DBCSR matrix pointer type
    1333              : !> \param[in] nmatrix         Size of set
    1334              : !> \param mmatrix ...
    1335              : !> \param pmatrix ...
    1336              : !> \par History
    1337              : !>      2009-08-17 Adapted from sparse_matrix_type for DBCSR
    1338              : ! **************************************************************************************************
    1339            0 :    SUBROUTINE allocate_dbcsr_matrix_set_3d(matrix_set, nmatrix, mmatrix, pmatrix)
    1340              :       TYPE(dbcsr_p_type), DIMENSION(:, :, :), POINTER    :: matrix_set
    1341              :       INTEGER, INTENT(IN)                                :: nmatrix, mmatrix, pmatrix
    1342              : 
    1343              :       INTEGER                                            :: imatrix, jmatrix, kmatrix
    1344              : 
    1345            0 :       IF (ASSOCIATED(matrix_set)) CALL dbcsr_deallocate_matrix_set(matrix_set)
    1346            0 :       ALLOCATE (matrix_set(nmatrix, mmatrix, pmatrix))
    1347            0 :       DO kmatrix = 1, pmatrix
    1348            0 :          DO jmatrix = 1, mmatrix
    1349            0 :             DO imatrix = 1, nmatrix
    1350            0 :                NULLIFY (matrix_set(imatrix, jmatrix, kmatrix)%matrix)
    1351              :             END DO
    1352              :          END DO
    1353              :       END DO
    1354            0 :    END SUBROUTINE allocate_dbcsr_matrix_set_3d
    1355              : 
    1356              : ! **************************************************************************************************
    1357              : !> \brief   Allocate and initialize a real matrix 4-dimensional set.
    1358              : !> \param[in,out] matrix_set  Set containing the DBCSR matrix pointer type
    1359              : !> \param[in] nmatrix         Size of set
    1360              : !> \param mmatrix ...
    1361              : !> \param pmatrix ...
    1362              : !> \param qmatrix ...
    1363              : !> \par History
    1364              : !>      2009-08-17 Adapted from sparse_matrix_type for DBCSR
    1365              : ! **************************************************************************************************
    1366            0 :    SUBROUTINE allocate_dbcsr_matrix_set_4d(matrix_set, nmatrix, mmatrix, pmatrix, qmatrix)
    1367              :       TYPE(dbcsr_p_type), DIMENSION(:, :, :, :), POINTER :: matrix_set
    1368              :       INTEGER, INTENT(IN)                                :: nmatrix, mmatrix, pmatrix, qmatrix
    1369              : 
    1370              :       INTEGER                                            :: imatrix, jmatrix, kmatrix, lmatrix
    1371              : 
    1372            0 :       IF (ASSOCIATED(matrix_set)) CALL dbcsr_deallocate_matrix_set(matrix_set)
    1373            0 :       ALLOCATE (matrix_set(nmatrix, mmatrix, pmatrix, qmatrix))
    1374            0 :       DO lmatrix = 1, qmatrix
    1375            0 :       DO kmatrix = 1, pmatrix
    1376            0 :          DO jmatrix = 1, mmatrix
    1377            0 :             DO imatrix = 1, nmatrix
    1378            0 :                NULLIFY (matrix_set(imatrix, jmatrix, kmatrix, lmatrix)%matrix)
    1379              :             END DO
    1380              :          END DO
    1381              :       END DO
    1382              :       END DO
    1383            0 :    END SUBROUTINE allocate_dbcsr_matrix_set_4d
    1384              : 
    1385              : ! **************************************************************************************************
    1386              : !> \brief   Allocate and initialize a real matrix 5-dimensional set.
    1387              : !> \param[in,out] matrix_set  Set containing the DBCSR matrix pointer type
    1388              : !> \param[in] nmatrix         Size of set
    1389              : !> \param mmatrix ...
    1390              : !> \param pmatrix ...
    1391              : !> \param qmatrix ...
    1392              : !> \param smatrix ...
    1393              : !> \par History
    1394              : !>      2009-08-17 Adapted from sparse_matrix_type for DBCSR
    1395              : ! **************************************************************************************************
    1396            0 :    SUBROUTINE allocate_dbcsr_matrix_set_5d(matrix_set, nmatrix, mmatrix, pmatrix, qmatrix, smatrix)
    1397              :       TYPE(dbcsr_p_type), DIMENSION(:, :, :, :, :), &
    1398              :          POINTER                                         :: matrix_set
    1399              :       INTEGER, INTENT(IN)                                :: nmatrix, mmatrix, pmatrix, qmatrix, &
    1400              :                                                             smatrix
    1401              : 
    1402              :       INTEGER                                            :: hmatrix, imatrix, jmatrix, kmatrix, &
    1403              :                                                             lmatrix
    1404              : 
    1405            0 :       IF (ASSOCIATED(matrix_set)) CALL dbcsr_deallocate_matrix_set(matrix_set)
    1406            0 :       ALLOCATE (matrix_set(nmatrix, mmatrix, pmatrix, qmatrix, smatrix))
    1407            0 :       DO hmatrix = 1, smatrix
    1408            0 :       DO lmatrix = 1, qmatrix
    1409            0 :       DO kmatrix = 1, pmatrix
    1410            0 :          DO jmatrix = 1, mmatrix
    1411            0 :             DO imatrix = 1, nmatrix
    1412            0 :                NULLIFY (matrix_set(imatrix, jmatrix, kmatrix, lmatrix, hmatrix)%matrix)
    1413              :             END DO
    1414              :          END DO
    1415              :       END DO
    1416              :       END DO
    1417              :       END DO
    1418            0 :    END SUBROUTINE allocate_dbcsr_matrix_set_5d
    1419              : 
    1420              :    ! **************************************************************************************************
    1421              : !> \brief Deallocate a real matrix set and release all of the member matrices.
    1422              : !> \param[in,out] matrix_set  Set containing the DBCSR matrix pointer type
    1423              : !> \par History
    1424              : !>      2009-08-17 Adapted from sparse_matrix_type for DBCSR
    1425              : ! **************************************************************************************************
    1426       180650 :    SUBROUTINE deallocate_dbcsr_matrix_set_1d(matrix_set)
    1427              : 
    1428              :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_set
    1429              : 
    1430              :       INTEGER                                            :: imatrix
    1431              : 
    1432       180650 :       IF (ASSOCIATED(matrix_set)) THEN
    1433       495234 :          DO imatrix = 1, SIZE(matrix_set)
    1434       495234 :             CALL dbcsr_deallocate_matrix(matrix_set(imatrix)%matrix)
    1435              :          END DO
    1436       179158 :          DEALLOCATE (matrix_set)
    1437              :       END IF
    1438              : 
    1439       180650 :    END SUBROUTINE deallocate_dbcsr_matrix_set_1d
    1440              : 
    1441              : ! **************************************************************************************************
    1442              : !> \brief Deallocate a real matrix set and release all of the member matrices.
    1443              : !> \param[in,out] matrix_set  Set containing the DBCSR matrix pointer type
    1444              : !> \par History
    1445              : !>      2009-08-17 Adapted from sparse_matrix_type for DBCSR
    1446              : ! **************************************************************************************************
    1447       155275 :    SUBROUTINE deallocate_dbcsr_matrix_set_2d(matrix_set)
    1448              : 
    1449              :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_set
    1450              : 
    1451              :       INTEGER                                            :: imatrix, jmatrix
    1452              : 
    1453       155275 :       IF (ASSOCIATED(matrix_set)) THEN
    1454       915032 :          DO jmatrix = 1, SIZE(matrix_set, 2)
    1455      1790631 :             DO imatrix = 1, SIZE(matrix_set, 1)
    1456      1638358 :                CALL dbcsr_deallocate_matrix(matrix_set(imatrix, jmatrix)%matrix)
    1457              :             END DO
    1458              :          END DO
    1459       152273 :          DEALLOCATE (matrix_set)
    1460              :       END IF
    1461       155275 :    END SUBROUTINE deallocate_dbcsr_matrix_set_2d
    1462              : 
    1463              : ! **************************************************************************************************
    1464              : !> \brief Deallocate a real matrix set and release all of the member matrices.
    1465              : !> \param[in,out] matrix_set  Set containing the DBCSR matrix pointer type
    1466              : !> \par History
    1467              : !>      2009-08-17 Adapted from sparse_matrix_type for DBCSR
    1468              : ! **************************************************************************************************
    1469            0 :    SUBROUTINE deallocate_dbcsr_matrix_set_3d(matrix_set)
    1470              : 
    1471              :       TYPE(dbcsr_p_type), DIMENSION(:, :, :), POINTER    :: matrix_set
    1472              : 
    1473              :       INTEGER                                            :: imatrix, jmatrix, kmatrix
    1474              : 
    1475            0 :       IF (ASSOCIATED(matrix_set)) THEN
    1476            0 :          DO kmatrix = 1, SIZE(matrix_set, 3)
    1477            0 :             DO jmatrix = 1, SIZE(matrix_set, 2)
    1478            0 :                DO imatrix = 1, SIZE(matrix_set, 1)
    1479            0 :                   CALL dbcsr_deallocate_matrix(matrix_set(imatrix, jmatrix, kmatrix)%matrix)
    1480              :                END DO
    1481              :             END DO
    1482              :          END DO
    1483            0 :          DEALLOCATE (matrix_set)
    1484              :       END IF
    1485            0 :    END SUBROUTINE deallocate_dbcsr_matrix_set_3d
    1486              : 
    1487              : ! **************************************************************************************************
    1488              : !> \brief Deallocate a real matrix set and release all of the member matrices.
    1489              : !> \param[in,out] matrix_set  Set containing the DBCSR matrix pointer type
    1490              : !> \par History
    1491              : !>      2009-08-17 Adapted from sparse_matrix_type for DBCSR
    1492              : ! **************************************************************************************************
    1493            0 :    SUBROUTINE deallocate_dbcsr_matrix_set_4d(matrix_set)
    1494              : 
    1495              :       TYPE(dbcsr_p_type), DIMENSION(:, :, :, :), POINTER :: matrix_set
    1496              : 
    1497              :       INTEGER                                            :: imatrix, jmatrix, kmatrix, lmatrix
    1498              : 
    1499            0 :       IF (ASSOCIATED(matrix_set)) THEN
    1500            0 :          DO lmatrix = 1, SIZE(matrix_set, 4)
    1501            0 :          DO kmatrix = 1, SIZE(matrix_set, 3)
    1502            0 :             DO jmatrix = 1, SIZE(matrix_set, 2)
    1503            0 :                DO imatrix = 1, SIZE(matrix_set, 1)
    1504            0 :                   CALL dbcsr_deallocate_matrix(matrix_set(imatrix, jmatrix, kmatrix, lmatrix)%matrix)
    1505              :                END DO
    1506              :             END DO
    1507              :          END DO
    1508              :          END DO
    1509            0 :          DEALLOCATE (matrix_set)
    1510              :       END IF
    1511            0 :    END SUBROUTINE deallocate_dbcsr_matrix_set_4d
    1512              : 
    1513              : ! **************************************************************************************************
    1514              : !> \brief Deallocate a real matrix set and release all of the member matrices.
    1515              : !> \param[in,out] matrix_set  Set containing the DBCSR matrix pointer type
    1516              : !> \par History
    1517              : !>      2009-08-17 Adapted from sparse_matrix_type for DBCSR
    1518              : ! **************************************************************************************************
    1519            0 :    SUBROUTINE deallocate_dbcsr_matrix_set_5d(matrix_set)
    1520              : 
    1521              :       TYPE(dbcsr_p_type), DIMENSION(:, :, :, :, :), &
    1522              :          POINTER                                         :: matrix_set
    1523              : 
    1524              :       INTEGER                                            :: hmatrix, imatrix, jmatrix, kmatrix, &
    1525              :                                                             lmatrix
    1526              : 
    1527            0 :       IF (ASSOCIATED(matrix_set)) THEN
    1528            0 :          DO hmatrix = 1, SIZE(matrix_set, 5)
    1529            0 :             DO lmatrix = 1, SIZE(matrix_set, 4)
    1530            0 :             DO kmatrix = 1, SIZE(matrix_set, 3)
    1531            0 :                DO jmatrix = 1, SIZE(matrix_set, 2)
    1532            0 :                   DO imatrix = 1, SIZE(matrix_set, 1)
    1533            0 :                      CALL dbcsr_deallocate_matrix(matrix_set(imatrix, jmatrix, kmatrix, lmatrix, hmatrix)%matrix)
    1534              :                   END DO
    1535              :                END DO
    1536              :             END DO
    1537              :             END DO
    1538              :          END DO
    1539            0 :          DEALLOCATE (matrix_set)
    1540              :       END IF
    1541            0 :    END SUBROUTINE deallocate_dbcsr_matrix_set_5d
    1542              : 
    1543              : END MODULE cp_dbcsr_operations
        

Generated by: LCOV version 2.0-1