LCOV - code coverage report
Current view: top level - src/fm - cp_fm_types.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:123285f) Lines: 88.2 % 820 723
Test Date: 2026-11-04 00:58:17 Functions: 80.9 % 47 38

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief represent a full matrix distributed on many processors
      10              : !> \par History
      11              : !>      3) separated structure object, removed globenv, renamed to full matrix
      12              : !>         many changes (fawzi 08.2002)
      13              : !> \author Matthias Krack (22.05.2001)
      14              : ! **************************************************************************************************
      15              : MODULE cp_fm_types
      16              :    USE OMP_LIB,                         ONLY: omp_get_max_threads
      17              :    USE cp_blacs_env,                    ONLY: cp_blacs_env_type
      18              :    USE cp_blacs_types,                  ONLY: cp_blacs_type
      19              :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      20              :                                               cp_fm_struct_equivalent,&
      21              :                                               cp_fm_struct_get,&
      22              :                                               cp_fm_struct_release,&
      23              :                                               cp_fm_struct_retain,&
      24              :                                               cp_fm_struct_type,&
      25              :                                               cp_fm_struct_write_info
      26              :    USE kinds,                           ONLY: dp
      27              :    USE message_passing,                 ONLY: cp2k_is_parallel,&
      28              :                                               mp_any_source,&
      29              :                                               mp_para_env_type,&
      30              :                                               mp_proc_null,&
      31              :                                               mp_request_null,&
      32              :                                               mp_request_type,&
      33              :                                               mp_waitall
      34              :    USE parallel_rng_types,              ONLY: UNIFORM,&
      35              :                                               rng_stream_type
      36              : #include "../base/base_uses.f90"
      37              : 
      38              :    IMPLICIT NONE
      39              : 
      40              :    PRIVATE
      41              : 
      42              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_fm_types'
      43              :    LOGICAL, PARAMETER          :: debug_this_module = .TRUE.
      44              :    INTEGER, PARAMETER :: src_tag = 3, dest_tag = 5, send_tag = 7, recv_tag = 11
      45              : 
      46              :    INTEGER, PRIVATE :: cp_fm_mm_type = 1
      47              : 
      48              :    PUBLIC :: cp_fm_type, &
      49              :              cp_fm_p_type, copy_info_type
      50              : 
      51              :    PUBLIC :: cp_fm_add_to_element, &
      52              :              cp_fm_create, &
      53              :              cp_fm_release, &
      54              :              cp_fm_get_info, &
      55              :              cp_fm_set_element, &
      56              :              cp_fm_get_element, &
      57              :              cp_fm_get_diag, & ! get diagonal
      58              :              cp_fm_set_all, & ! set all elements and diagonal
      59              :              cp_fm_set_all_submatrix, & ! set a submatrix to a given value
      60              :              cp_fm_set_submatrix, & ! set a submatrix to given values
      61              :              cp_fm_get_submatrix, & ! get a submatrix of given values
      62              :              cp_fm_init_random, &
      63              :              cp_fm_maxabsval, & ! find the maximum absolute value
      64              :              cp_fm_maxabsrownorm, & ! find the maximum of the sum of the abs of the elements of a row
      65              :              cp_fm_to_fm, & ! copy (parts of) a fm to a fm
      66              :              cp_fm_vectorsnorm, & ! compute the norm of the column-vectors
      67              :              cp_fm_vectorssum, & ! compute the sum of all elements of the column-vectors
      68              :              cp_fm_sumall, & ! sum all elements or a function of them of the matrix
      69              :              cp_fm_to_fm_submat, & ! copy (parts of) a fm to a fm
      70              :              cp_fm_to_fm_triangular, &
      71              :              cp_fm_copy_general, &
      72              :              cp_fm_start_copy_general, &
      73              :              cp_fm_finish_copy_general, &
      74              :              cp_fm_cleanup_copy_general, &
      75              :              cp_fm_write_unformatted, & ! writes a full matrix to an open unit
      76              :              cp_fm_write_formatted, & ! writes a full matrix to an open unit
      77              :              cp_fm_read_unformatted, & ! reads a full matrix from an open unit
      78              :              cp_fm_setup, & ! allows to set flags for fms
      79              :              cp_fm_get_mm_type, &
      80              :              cp_fm_write_info, &
      81              :              cp_fm_to_fm_submat_general ! copy matrix across different contexts
      82              : 
      83              :    PUBLIC :: cp_fm_pilaenv
      84              : 
      85              :    INTERFACE cp_fm_to_fm
      86              :       MODULE PROCEDURE cp_fm_to_fm_matrix, & ! a full matrix
      87              :          cp_fm_to_fm_columns ! just a number of columns
      88              :    END INTERFACE
      89              : 
      90              :    INTERFACE cp_fm_release
      91              :       MODULE PROCEDURE cp_fm_release_aa0, &
      92              :          cp_fm_release_aa1, &
      93              :          cp_fm_release_aa2, &
      94              :          cp_fm_release_aa3, &
      95              :          cp_fm_release_ap1, &
      96              :          cp_fm_release_ap2, &
      97              :          cp_fm_release_pa1, &
      98              :          cp_fm_release_pa2, &
      99              :          cp_fm_release_pa3, &
     100              :          cp_fm_release_pp1, &
     101              :          cp_fm_release_pp2
     102              :    END INTERFACE
     103              : 
     104              : ! **************************************************************************************************
     105              : !> \brief represent a full matrix
     106              : !> \param name the name of the matrix, used for printing
     107              : !> \param matrix_struct structure of this matrix
     108              : !> \param local_data array with the data of the matrix (its contents
     109              : !>        depend on the matrix type used: in parallel runs it will be
     110              : !>        in scalapack format, in sequential, it will simply contain
     111              : !>        the matrix)
     112              : !> \par History
     113              : !>      08.2002 created [fawzi]
     114              : !> \author fawzi
     115              : ! **************************************************************************************************
     116              :    TYPE cp_fm_type
     117              : !    PRIVATE
     118              :       CHARACTER(LEN=60) :: name = ""
     119              :       TYPE(cp_fm_struct_type), POINTER :: matrix_struct => NULL()
     120              :       REAL(KIND=dp), DIMENSION(:, :), POINTER, CONTIGUOUS :: local_data => NULL()
     121              :    END TYPE cp_fm_type
     122              : 
     123              : ! **************************************************************************************************
     124              : !> \brief just to build arrays of pointers to matrices
     125              : !> \param matrix the pointer to the matrix
     126              : !> \par History
     127              : !>      08.2002 created [fawzi]
     128              : !> \author fawzi
     129              : ! **************************************************************************************************
     130              :    TYPE cp_fm_p_type
     131              :       TYPE(cp_fm_type), POINTER :: matrix => NULL()
     132              :    END TYPE cp_fm_p_type
     133              : 
     134              : ! **************************************************************************************************
     135              : !> \brief Stores the state of a copy between cp_fm_start_copy_general
     136              : !>        and cp_fm_finish_copy_general
     137              : !> \par History
     138              : !>      Jan 2017  [Mark T]
     139              : ! **************************************************************************************************
     140              :    TYPE copy_info_type
     141              :       INTEGER :: send_size = -1
     142              :       INTEGER, DIMENSION(2) :: nlocal_recv = -1, nblock_src = -1, src_num_pe = -1 ! 1->row  2->col
     143              :       TYPE(mp_request_type), DIMENSION(:), ALLOCATABLE :: send_request, recv_request
     144              :       INTEGER, DIMENSION(:), ALLOCATABLE   :: recv_disp
     145              :       INTEGER, DIMENSION(:), POINTER       :: recv_col_indices => NULL(), recv_row_indices => NULL()
     146              :       INTEGER, DIMENSION(:, :), ALLOCATABLE :: src_blacs2mpi
     147              :       REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: recv_buf, send_buf
     148              :    END TYPE copy_info_type
     149              : 
     150              : CONTAINS
     151              : 
     152              : ! **************************************************************************************************
     153              : !> \brief creates a new full matrix with the given structure
     154              : !> \param matrix the matrix to be created
     155              : !> \param matrix_struct the structure of matrix
     156              : !> \param name ...
     157              : !> \param nrow ...
     158              : !> \param ncol ...
     159              : !> \param set_zero ...
     160              : !> \par History
     161              : !>      08.2002 created [fawzi]
     162              : !> \author Fawzi Mohamed
     163              : !> \note
     164              : !>      preferred allocation routine
     165              : ! **************************************************************************************************
     166      1609673 :    SUBROUTINE cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
     167              :       TYPE(cp_fm_type), INTENT(OUT)                      :: matrix
     168              :       TYPE(cp_fm_struct_type), INTENT(IN), TARGET        :: matrix_struct
     169              :       CHARACTER(LEN=*), INTENT(in), OPTIONAL             :: name
     170              :       INTEGER, INTENT(IN), OPTIONAL                      :: nrow, ncol
     171              :       LOGICAL, INTENT(in), OPTIONAL                      :: set_zero
     172              : 
     173              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'cp_fm_create'
     174              : 
     175              :       INTEGER                                            :: handle, ncol_global, ncol_local, &
     176              :                                                             nrow_global, nrow_local
     177              :       TYPE(cp_blacs_env_type), POINTER                   :: context
     178              :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct
     179              : 
     180      1609673 :       CALL timeset(routineN, handle)
     181              : 
     182      1609673 :       IF (PRESENT(nrow) .OR. PRESENT(ncol)) THEN
     183        26236 :          CALL cp_fm_struct_get(matrix_struct, nrow_global=nrow_global, ncol_global=ncol_global)
     184        26236 :          IF (PRESENT(nrow)) nrow_global = nrow
     185        26236 :          IF (PRESENT(ncol)) ncol_global = ncol
     186              :          CALL cp_fm_struct_create(fm_struct, template_fmstruct=matrix_struct, &
     187        26236 :                                   nrow_global=nrow_global, ncol_global=ncol_global)
     188              : 
     189        26236 :          context => fm_struct%context
     190        26236 :          matrix%matrix_struct => fm_struct
     191        26236 :          CALL cp_fm_struct_retain(matrix%matrix_struct)
     192              : 
     193        26236 :          nrow_local = fm_struct%local_leading_dimension
     194        26236 :          ncol_local = MAX(1, fm_struct%ncol_locals(context%mepos(2)))
     195              : 
     196        26236 :          CALL cp_fm_struct_release(fm_struct)
     197              :       ELSE
     198              : 
     199      1583437 :          context => matrix_struct%context
     200      1583437 :          matrix%matrix_struct => matrix_struct
     201      1583437 :          CALL cp_fm_struct_retain(matrix%matrix_struct)
     202              : 
     203              :          ! OK, we allocate here at least a 1 x 1 matrix
     204              :          ! this must (and is) compatible with the descinit call
     205              :          ! in cp_fm_struct
     206      1583437 :          nrow_local = matrix_struct%local_leading_dimension
     207      1583437 :          ncol_local = MAX(1, matrix_struct%ncol_locals(context%mepos(2)))
     208              :       END IF
     209              : 
     210      1609673 :       NULLIFY (matrix%local_data)
     211              : 
     212      6438692 :       ALLOCATE (matrix%local_data(nrow_local, ncol_local))
     213              : 
     214      1609673 :       IF (PRESENT(set_zero)) THEN
     215       257445 :          IF (set_zero) THEN
     216     86051847 :             matrix%local_data(1:nrow_local, 1:ncol_local) = 0.0_dp
     217              :          END IF
     218              :       END IF
     219              : 
     220      1609673 :       IF (PRESENT(name)) THEN
     221       666273 :          matrix%name = name
     222              :       ELSE
     223       943400 :          matrix%name = 'full matrix'
     224              :       END IF
     225              : 
     226      1609673 :       CALL timestop(handle)
     227              : 
     228      1609673 :    END SUBROUTINE cp_fm_create
     229              : 
     230              : ! **************************************************************************************************
     231              : !> \brief releases a full matrix
     232              : !> \param matrix the matrix to release
     233              : !> \par History
     234              : !>      08.2002 created [fawzi]
     235              : !> \author Fawzi Mohamed
     236              : ! **************************************************************************************************
     237      1665822 :    SUBROUTINE cp_fm_release_aa0(matrix)
     238              :       TYPE(cp_fm_type), INTENT(INOUT)                    :: matrix
     239              : 
     240      1665822 :       IF (ASSOCIATED(matrix%local_data)) THEN
     241      1608955 :          DEALLOCATE (matrix%local_data)
     242              :          NULLIFY (matrix%local_data)
     243              :       END IF
     244      1665822 :       matrix%name = ""
     245      1665822 :       CALL cp_fm_struct_release(matrix%matrix_struct)
     246              : 
     247      1665822 :    END SUBROUTINE cp_fm_release_aa0
     248              : 
     249              : ! **************************************************************************************************
     250              : !> \brief ...
     251              : !> \param matrices ...
     252              : ! **************************************************************************************************
     253        85314 :    SUBROUTINE cp_fm_release_aa1(matrices)
     254              :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: matrices
     255              : 
     256              :       INTEGER                                            :: i
     257              : 
     258        85314 :       IF (ALLOCATED(matrices)) THEN
     259       225095 :          DO i = 1, SIZE(matrices)
     260       225095 :             CALL cp_fm_release(matrices(i))
     261              :          END DO
     262        81244 :          DEALLOCATE (matrices)
     263              :       END IF
     264        85314 :    END SUBROUTINE cp_fm_release_aa1
     265              : 
     266              : ! **************************************************************************************************
     267              : !> \brief ...
     268              : !> \param matrices ...
     269              : ! **************************************************************************************************
     270         8884 :    SUBROUTINE cp_fm_release_aa2(matrices)
     271              :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :)     :: matrices
     272              : 
     273              :       INTEGER                                            :: i, j
     274              : 
     275         8884 :       IF (ALLOCATED(matrices)) THEN
     276        16160 :          DO i = 1, SIZE(matrices, 1)
     277        44570 :             DO j = 1, SIZE(matrices, 2)
     278        38574 :                CALL cp_fm_release(matrices(i, j))
     279              :             END DO
     280              :          END DO
     281         5996 :          DEALLOCATE (matrices)
     282              :       END IF
     283         8884 :    END SUBROUTINE cp_fm_release_aa2
     284              : 
     285              : ! **************************************************************************************************
     286              : !> \brief ...
     287              : !> \param matrices ...
     288              : ! **************************************************************************************************
     289          142 :    SUBROUTINE cp_fm_release_aa3(matrices)
     290              :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :, :)  :: matrices
     291              : 
     292              :       INTEGER                                            :: i, j, k
     293              : 
     294          142 :       IF (ALLOCATED(matrices)) THEN
     295         2224 :          DO i = 1, SIZE(matrices, 1)
     296         6416 :             DO j = 1, SIZE(matrices, 2)
     297        11146 :                DO k = 1, SIZE(matrices, 3)
     298         9064 :                   CALL cp_fm_release(matrices(i, j, k))
     299              :                END DO
     300              :             END DO
     301              :          END DO
     302          142 :          DEALLOCATE (matrices)
     303              :       END IF
     304          142 :    END SUBROUTINE cp_fm_release_aa3
     305              : 
     306              : ! **************************************************************************************************
     307              : !> \brief ...
     308              : !> \param matrices ...
     309              : ! **************************************************************************************************
     310       217359 :    SUBROUTINE cp_fm_release_pa1(matrices)
     311              :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: matrices
     312              : 
     313              :       INTEGER                                            :: i
     314              : 
     315       217359 :       IF (ASSOCIATED(matrices)) THEN
     316       152188 :          DO i = 1, SIZE(matrices)
     317       152188 :             CALL cp_fm_release(matrices(i))
     318              :          END DO
     319        62209 :          DEALLOCATE (matrices)
     320              :          NULLIFY (matrices)
     321              :       END IF
     322       217359 :    END SUBROUTINE cp_fm_release_pa1
     323              : 
     324              : ! **************************************************************************************************
     325              : !> \brief ...
     326              : !> \param matrices ...
     327              : ! **************************************************************************************************
     328        70187 :    SUBROUTINE cp_fm_release_pa2(matrices)
     329              :       TYPE(cp_fm_type), DIMENSION(:, :), POINTER         :: matrices
     330              : 
     331              :       INTEGER                                            :: i, j
     332              : 
     333        70187 :       IF (ASSOCIATED(matrices)) THEN
     334        87297 :          DO i = 1, SIZE(matrices, 1)
     335       162613 :             DO j = 1, SIZE(matrices, 2)
     336       138216 :                CALL cp_fm_release(matrices(i, j))
     337              :             END DO
     338              :          END DO
     339        24397 :          DEALLOCATE (matrices)
     340              :          NULLIFY (matrices)
     341              :       END IF
     342        70187 :    END SUBROUTINE cp_fm_release_pa2
     343              : 
     344              : ! **************************************************************************************************
     345              : !> \brief ...
     346              : !> \param matrices ...
     347              : ! **************************************************************************************************
     348            0 :    SUBROUTINE cp_fm_release_pa3(matrices)
     349              :       TYPE(cp_fm_type), DIMENSION(:, :, :), POINTER      :: matrices
     350              : 
     351              :       INTEGER                                            :: i, j, k
     352              : 
     353            0 :       IF (ASSOCIATED(matrices)) THEN
     354            0 :          DO i = 1, SIZE(matrices, 1)
     355            0 :             DO j = 1, SIZE(matrices, 2)
     356            0 :                DO k = 1, SIZE(matrices, 3)
     357            0 :                   CALL cp_fm_release(matrices(i, j, k))
     358              :                END DO
     359              :             END DO
     360              :          END DO
     361            0 :          DEALLOCATE (matrices)
     362              :          NULLIFY (matrices)
     363              :       END IF
     364            0 :    END SUBROUTINE cp_fm_release_pa3
     365              : 
     366              : ! **************************************************************************************************
     367              : !> \brief ...
     368              : !> \param matrices ...
     369              : ! **************************************************************************************************
     370            0 :    SUBROUTINE cp_fm_release_ap1(matrices)
     371              :       TYPE(cp_fm_p_type), ALLOCATABLE, DIMENSION(:)      :: matrices
     372              : 
     373              :       INTEGER                                            :: i
     374              : 
     375            0 :       IF (ALLOCATED(matrices)) THEN
     376            0 :          DO i = 1, SIZE(matrices)
     377            0 :             CALL cp_fm_release(matrices(i)%matrix)
     378            0 :             DEALLOCATE (matrices(i)%matrix)
     379              :          END DO
     380            0 :          DEALLOCATE (matrices)
     381              :       END IF
     382            0 :    END SUBROUTINE cp_fm_release_ap1
     383              : 
     384              : ! **************************************************************************************************
     385              : !> \brief ...
     386              : !> \param matrices ...
     387              : ! **************************************************************************************************
     388            0 :    SUBROUTINE cp_fm_release_ap2(matrices)
     389              :       TYPE(cp_fm_p_type), ALLOCATABLE, DIMENSION(:, :)   :: matrices
     390              : 
     391              :       INTEGER                                            :: i, j
     392              : 
     393            0 :       IF (ALLOCATED(matrices)) THEN
     394            0 :          DO i = 1, SIZE(matrices, 1)
     395            0 :             DO j = 1, SIZE(matrices, 2)
     396            0 :                CALL cp_fm_release(matrices(i, j)%matrix)
     397            0 :                DEALLOCATE (matrices(i, j)%matrix)
     398              :             END DO
     399              :          END DO
     400            0 :          DEALLOCATE (matrices)
     401              :       END IF
     402            0 :    END SUBROUTINE cp_fm_release_ap2
     403              : 
     404              : ! **************************************************************************************************
     405              : !> \brief ...
     406              : !> \param matrices ...
     407              : ! **************************************************************************************************
     408            0 :    SUBROUTINE cp_fm_release_pp1(matrices)
     409              :       TYPE(cp_fm_p_type), DIMENSION(:), POINTER          :: matrices
     410              : 
     411              :       INTEGER                                            :: i
     412              : 
     413            0 :       IF (ASSOCIATED(matrices)) THEN
     414            0 :          DO i = 1, SIZE(matrices)
     415            0 :             CALL cp_fm_release(matrices(i)%matrix)
     416            0 :             DEALLOCATE (matrices(i)%matrix)
     417              :          END DO
     418            0 :          DEALLOCATE (matrices)
     419              :          NULLIFY (matrices)
     420              :       END IF
     421            0 :    END SUBROUTINE cp_fm_release_pp1
     422              : 
     423              : ! **************************************************************************************************
     424              : !> \brief ...
     425              : !> \param matrices ...
     426              : ! **************************************************************************************************
     427            0 :    SUBROUTINE cp_fm_release_pp2(matrices)
     428              :       TYPE(cp_fm_p_type), DIMENSION(:, :), POINTER       :: matrices
     429              : 
     430              :       INTEGER                                            :: i, j
     431              : 
     432            0 :       IF (ASSOCIATED(matrices)) THEN
     433            0 :          DO i = 1, SIZE(matrices, 1)
     434            0 :             DO j = 1, SIZE(matrices, 2)
     435            0 :                CALL cp_fm_release(matrices(i, j)%matrix)
     436            0 :                DEALLOCATE (matrices(i, j)%matrix)
     437              :             END DO
     438              :          END DO
     439            0 :          DEALLOCATE (matrices)
     440              :          NULLIFY (matrices)
     441              :       END IF
     442            0 :    END SUBROUTINE cp_fm_release_pp2
     443              : 
     444              : ! **************************************************************************************************
     445              : !> \brief fills a matrix with random numbers
     446              : !> \param matrix : to be initialized
     447              : !> \param ncol : numbers of cols to fill
     448              : !> \param start_col : starting at coll number
     449              : !> \author Joost VandeVondele
     450              : !> \note
     451              : !>      the value of a_ij is independent of the number of cpus
     452              : ! **************************************************************************************************
     453         3262 :    SUBROUTINE cp_fm_init_random(matrix, ncol, start_col)
     454              :       TYPE(cp_fm_type), INTENT(IN)                       :: matrix
     455              :       INTEGER, INTENT(IN), OPTIONAL                      :: ncol, start_col
     456              : 
     457              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'cp_fm_init_random'
     458              : 
     459              :       INTEGER :: handle, icol_global, icol_local, irow_local, my_ncol, my_start_col, ncol_global, &
     460              :          ncol_local, nrow_global, nrow_local
     461         6524 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
     462              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: buff
     463              :       REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
     464         3262 :          POINTER                                         :: local_data
     465              :       REAL(KIND=dp), DIMENSION(3, 2), SAVE :: &
     466              :          seed = RESHAPE([1.0_dp, 2.0_dp, 3.0_dp, 4.0_dp, 5.0_dp, 6.0_dp], [3, 2])
     467              :       TYPE(rng_stream_type)                              :: rng
     468              : 
     469         3262 :       CALL timeset(routineN, handle)
     470              : 
     471              :       ! guarantee same seed over all tasks
     472         3262 :       CALL matrix%matrix_struct%para_env%bcast(seed, 0)
     473              : 
     474              :       rng = rng_stream_type("cp_fm_init_random_stream", distribution_type=UNIFORM, &
     475         3262 :                             extended_precision=.TRUE., seed=seed)
     476              : 
     477              :       CALL cp_fm_get_info(matrix, nrow_global=nrow_global, ncol_global=ncol_global, &
     478              :                           nrow_local=nrow_local, ncol_local=ncol_local, &
     479              :                           local_data=local_data, &
     480         3262 :                           row_indices=row_indices, col_indices=col_indices)
     481              : 
     482         3262 :       my_start_col = 1
     483         3262 :       IF (PRESENT(start_col)) my_start_col = start_col
     484         3262 :       my_ncol = matrix%matrix_struct%ncol_global
     485         3262 :       IF (PRESENT(ncol)) my_ncol = ncol
     486              : 
     487         3262 :       IF (ncol_global < (my_start_col + my_ncol - 1)) THEN
     488            0 :          CPABORT("ncol_global>=(my_start_col+my_ncol-1)")
     489              :       END IF
     490              : 
     491         9786 :       ALLOCATE (buff(nrow_global))
     492              : 
     493              :       ! each global row has its own substream, in order to reach the stream for the local col,
     494              :       ! we just reset to the next substream
     495              :       ! following this, we fill the full buff with random numbers, and pick those we need
     496         3262 :       icol_global = 0
     497        24969 :       DO icol_local = 1, ncol_local
     498        21707 :          CPASSERT(col_indices(icol_local) > icol_global)
     499              :          DO
     500        21707 :             CALL rng%reset_to_next_substream()
     501        21707 :             icol_global = icol_global + 1
     502        21707 :             IF (icol_global == col_indices(icol_local)) EXIT
     503              :          END DO
     504        21707 :          CALL rng%fill(buff)
     505       840665 :          DO irow_local = 1, nrow_local
     506       837403 :             local_data(irow_local, icol_local) = buff(row_indices(irow_local))
     507              :          END DO
     508              :       END DO
     509              : 
     510         3262 :       DEALLOCATE (buff)
     511              : 
     512              :       ! store seed before deletion (unclear if this is the proper seed)
     513              : 
     514              :       ! Note that, the initial state (rng%ig) instead of the current state (rng%cg) is stored in the
     515              :       ! seed variable. As a consequence, each invocation of cp_fm_init_random uses exactly the same
     516              :       ! stream of random numbers. While this seems odd and contrary to the original design,
     517              :       ! it was probably introduced to improve reproducibility.
     518              :       ! See also https://github.com/cp2k/cp2k/pull/506
     519         3262 :       CALL rng%get(ig=seed)
     520              : 
     521         3262 :       CALL timestop(handle)
     522              : 
     523        91336 :    END SUBROUTINE cp_fm_init_random
     524              : 
     525              : ! **************************************************************************************************
     526              : !> \brief set all elements of a matrix to the same value,
     527              : !>      and optionally the diagonal to a different one
     528              : !> \param matrix input matrix
     529              : !> \param alpha scalar used to set all elements of the matrix
     530              : !> \param beta scalar used to set diagonal of the matrix
     531              : !> \note
     532              : !>      can be used to zero a matrix
     533              : !>      can be used to create a unit matrix (I-matrix) alpha=0.0_dp beta=1.0_dp
     534              : ! **************************************************************************************************
     535       389265 :    SUBROUTINE cp_fm_set_all(matrix, alpha, beta)
     536              : 
     537              :       TYPE(cp_fm_type), INTENT(IN)                       :: matrix
     538              :       REAL(KIND=dp), INTENT(IN)                          :: alpha
     539              :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: beta
     540              : 
     541              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'cp_fm_set_all'
     542              : 
     543              :       INTEGER                                            :: handle, i, n
     544              : 
     545       389265 :       CALL timeset(routineN, handle)
     546              : 
     547    225391216 :       matrix%local_data(:, :) = alpha
     548              : 
     549       389265 :       IF (PRESENT(beta)) THEN
     550        67218 :          n = MIN(matrix%matrix_struct%nrow_global, matrix%matrix_struct%ncol_global)
     551       520198 :          DO i = 1, n
     552       520198 :             CALL cp_fm_set_element(matrix, i, i, beta)
     553              :          END DO
     554              :       END IF
     555              : 
     556       389265 :       CALL timestop(handle)
     557              : 
     558       389265 :    END SUBROUTINE cp_fm_set_all
     559              : 
     560              : ! **************************************************************************************************
     561              : !> \brief returns the diagonal elements of a fm
     562              : !> \param matrix ...
     563              : !> \param diag ...
     564              : !> \param n_active optional number of leading diagonal elements to return
     565              : ! **************************************************************************************************
     566        32236 :    SUBROUTINE cp_fm_get_diag(matrix, diag, n_active)
     567              : 
     568              :       TYPE(cp_fm_type), INTENT(IN)                       :: matrix
     569              :       REAL(KIND=dp), DIMENSION(:), INTENT(OUT)           :: diag
     570              :       INTEGER, INTENT(IN), OPTIONAL                      :: n_active
     571              : 
     572              :       INTEGER                                            :: i_global, n_diag, ncol_global, nrow_global
     573              : #if defined(__parallel)
     574              :       INTEGER                                            :: icol_local, ipcol, iprow, irow, irow_local, &
     575              :                                                             mypcol, myprow, npcol, nprow, nrow_local
     576              :       INTEGER, DIMENSION(9)                              :: desca
     577        32236 :       INTEGER, DIMENSION(:), POINTER                     :: row_indices
     578              :       TYPE(cp_blacs_env_type), POINTER                   :: context
     579              : #endif
     580              : 
     581        32236 :       CALL cp_fm_get_info(matrix, nrow_global=nrow_global, ncol_global=ncol_global)
     582              : 
     583        32236 :       n_diag = MIN(nrow_global, ncol_global)
     584        32236 :       IF (PRESENT(n_active)) n_diag = n_active
     585        32236 :       CPASSERT(n_diag >= 0 .AND. n_diag <= MIN(nrow_global, ncol_global))
     586        32236 :       CPASSERT(SIZE(diag) >= n_diag)
     587              : 
     588      2232340 :       diag = 0.0_dp
     589              : #if defined(__parallel)
     590        32236 :       CALL cp_fm_get_info(matrix, nrow_local=nrow_local, row_indices=row_indices)
     591        32236 :       context => matrix%matrix_struct%context
     592        32236 :       myprow = context%mepos(1)
     593        32236 :       mypcol = context%mepos(2)
     594        32236 :       nprow = context%num_pe(1)
     595        32236 :       npcol = context%num_pe(2)
     596              : 
     597       322360 :       desca(:) = matrix%matrix_struct%descriptor(:)
     598              : 
     599      1186125 :       DO irow = 1, nrow_local
     600      1153889 :          i_global = row_indices(irow)
     601      1153889 :          IF (i_global > n_diag) EXIT
     602              :          CALL infog2l(i_global, i_global, desca, nprow, npcol, myprow, mypcol, &
     603      1153889 :                       irow_local, icol_local, iprow, ipcol)
     604      1186125 :          IF ((iprow == myprow) .AND. (ipcol == mypcol)) THEN
     605      1153889 :             diag(i_global) = matrix%local_data(irow_local, icol_local)
     606              :          END IF
     607              :       END DO
     608              : #else
     609              :       DO i_global = 1, n_diag
     610              :          diag(i_global) = matrix%local_data(i_global, i_global)
     611              :       END DO
     612              : #endif
     613      4432444 :       CALL matrix%matrix_struct%para_env%sum(diag)
     614              : 
     615        32236 :    END SUBROUTINE cp_fm_get_diag
     616              : 
     617              : ! **************************************************************************************************
     618              : !> \brief returns an element of a fm
     619              : !>      this value is valid on every cpu
     620              : !>      using this call is expensive
     621              : !> \param matrix the matrix to read
     622              : !> \param irow_global the row
     623              : !> \param icol_global the col
     624              : !> \param alpha the value of matrix(irow_global, icol_global)
     625              : !> \param local true if the element is on this cpu, false otherwise
     626              : !> \note
     627              : !>      - modified semantics. now this function always returns the value
     628              : !>        previously the value was zero on cpus that didn't own the relevant
     629              : !>        part of the matrix (Joost VandeVondele, May 2003)
     630              : !>      - usage of the function should be avoided, as it is likely to rather slow
     631              : !>        using row_indices/col_indices/local_data + some smart scheme normally
     632              : !>        yields a real parallel code
     633              : ! **************************************************************************************************
     634      1097784 :    SUBROUTINE cp_fm_get_element(matrix, irow_global, icol_global, alpha, local)
     635              : 
     636              :       ! arguments
     637              :       TYPE(cp_fm_type), INTENT(IN)          :: matrix
     638              :       REAL(KIND=dp), INTENT(OUT)                     :: alpha
     639              :       INTEGER, INTENT(IN)                       :: icol_global, &
     640              :                                                    irow_global
     641              :       LOGICAL, INTENT(OUT), OPTIONAL            :: local
     642              : 
     643              :       ! locals
     644              : #if defined(__parallel)
     645              :       INTEGER, DIMENSION(9) :: desca
     646              :       TYPE(cp_blacs_env_type), POINTER :: context
     647              :       INTEGER :: icol_local, ipcol, iprow, irow_local, mypcol, myprow, npcol, &
     648              :                  nprow
     649      1097784 :       REAL(KIND=dp), DIMENSION(:, :), POINTER :: a
     650              : #endif
     651              : 
     652              : #if defined(__parallel)
     653      1097784 :       context => matrix%matrix_struct%context
     654      1097784 :       myprow = context%mepos(1)
     655      1097784 :       mypcol = context%mepos(2)
     656      1097784 :       nprow = context%num_pe(1)
     657      1097784 :       npcol = context%num_pe(2)
     658              : 
     659      1097784 :       a => matrix%local_data
     660     10977840 :       desca(:) = matrix%matrix_struct%descriptor(:)
     661              : 
     662              :       CALL infog2l(irow_global, icol_global, desca, nprow, npcol, myprow, mypcol, &
     663      1097784 :                    irow_local, icol_local, iprow, ipcol)
     664              : 
     665      1097784 :       IF ((iprow == myprow) .AND. (ipcol == mypcol)) THEN
     666       548942 :          alpha = a(irow_local, icol_local)
     667       548942 :          CALL context%dgebs2d('A', ' ', 1, 1, alpha, 1)
     668       548942 :          IF (PRESENT(local)) local = .TRUE.
     669              :       ELSE
     670       548842 :          CALL context%dgebr2d('A', ' ', 1, 1, alpha, 1, iprow, ipcol)
     671       548842 :          IF (PRESENT(local)) local = .FALSE.
     672              :       END IF
     673              : 
     674              : #else
     675              :       IF (PRESENT(local)) local = .TRUE.
     676              :       alpha = matrix%local_data(irow_global, icol_global)
     677              : #endif
     678              : 
     679      1097784 :    END SUBROUTINE cp_fm_get_element
     680              : 
     681              : ! **************************************************************************************************
     682              : !> \brief sets an element of a matrix
     683              : !> \param matrix ...
     684              : !> \param irow_global ...
     685              : !> \param icol_global ...
     686              : !> \param alpha ...
     687              : !> \note
     688              : !>      we expect all cpus to have the same arguments in the call to this function
     689              : !>      (otherwise one should use local_data tricks)
     690              : ! **************************************************************************************************
     691       623550 :    SUBROUTINE cp_fm_set_element(matrix, irow_global, icol_global, alpha)
     692              :       TYPE(cp_fm_type), INTENT(IN)          :: matrix
     693              :       INTEGER, INTENT(IN)                      :: irow_global, icol_global
     694              :       REAL(KIND=dp), INTENT(IN)                :: alpha
     695              : 
     696              :       INTEGER                                  :: mypcol, myprow, npcol, nprow
     697              :       TYPE(cp_blacs_env_type), POINTER         :: context
     698              : #if defined(__parallel)
     699              :       INTEGER                                  :: icol_local, ipcol, iprow, &
     700              :                                                   irow_local
     701              :       INTEGER, DIMENSION(9)                    :: desca
     702       623550 :       REAL(KIND=dp), DIMENSION(:, :), POINTER  :: a
     703              : #endif
     704              : 
     705       623550 :       context => matrix%matrix_struct%context
     706       623550 :       myprow = context%mepos(1)
     707       623550 :       mypcol = context%mepos(2)
     708       623550 :       nprow = context%num_pe(1)
     709       623550 :       npcol = context%num_pe(2)
     710              : 
     711              : #if defined(__parallel)
     712              : 
     713       623550 :       a => matrix%local_data
     714              : 
     715      6235500 :       desca(:) = matrix%matrix_struct%descriptor(:)
     716              : 
     717              :       CALL infog2l(irow_global, icol_global, desca, nprow, npcol, myprow, mypcol, &
     718       623550 :                    irow_local, icol_local, iprow, ipcol)
     719              : 
     720       623550 :       IF ((iprow == myprow) .AND. (ipcol == mypcol)) THEN
     721       313680 :          a(irow_local, icol_local) = alpha
     722              :       END IF
     723              : 
     724              : #else
     725              : 
     726              :       matrix%local_data(irow_global, icol_global) = alpha
     727              : 
     728              : #endif
     729       623550 :    END SUBROUTINE cp_fm_set_element
     730              : 
     731              : ! **************************************************************************************************
     732              : !> \brief sets a submatrix of a full matrix
     733              : !>       fm(start_row:start_row+n_rows,start_col:start_col+n_cols)
     734              : !>       = alpha*op(new_values)(1:n_rows,1:n_cols)+ beta
     735              : !>       * fm(start_row:start_row+n_rows,start_col:start_col+n_cols)
     736              : !> \param fm the full to change
     737              : !> \param new_values a replicated full matrix with the new values
     738              : !> \param start_row the starting row of b_matrix (defaults to 1)
     739              : !> \param start_col the starting col of b_matrix (defaults to 1)
     740              : !> \param n_rows the number of row to change in b (defaults to
     741              : !>        size(op(new_values),1))
     742              : !> \param n_cols the number of columns to change in b (defaults to
     743              : !>        size(op(new_values),2))
     744              : !> \param alpha rescaling factor for the new values (defaults to 1.0)
     745              : !> \param beta rescaling factor for the old values (defaults to 0.0)
     746              : !> \param transpose if new_values should be transposed: if true
     747              : !>        op(new_values)=new_values^T, else op(new_values)=new_values
     748              : !>        (defaults to false)
     749              : !> \par History
     750              : !>      07.2002 created borrowing from Joost's blacs_replicated_copy [fawzi]
     751              : !> \author Fawzi Mohamed
     752              : !> \note
     753              : !>      optimized for full column updates and alpha=1.0, beta=0.0
     754              : !>      the new_values need to be valid on all cpus
     755              : ! **************************************************************************************************
     756        90247 :    SUBROUTINE cp_fm_set_submatrix(fm, new_values, start_row, &
     757              :                                   start_col, n_rows, n_cols, alpha, beta, transpose)
     758              :       TYPE(cp_fm_type), INTENT(IN)                       :: fm
     759              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(in)         :: new_values
     760              :       INTEGER, INTENT(in), OPTIONAL                      :: start_row, start_col, n_rows, n_cols
     761              :       REAL(KIND=dp), INTENT(in), OPTIONAL                :: alpha, beta
     762              :       LOGICAL, INTENT(in), OPTIONAL                      :: transpose
     763              : 
     764              :       INTEGER                                            :: i, i0, j, j0, ncol, ncol_block, &
     765              :                                                             ncol_global, ncol_local, nrow, &
     766              :                                                             nrow_block, nrow_global, nrow_local, &
     767              :                                                             this_col, this_row
     768        90247 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
     769              :       LOGICAL                                            :: tr_a
     770              :       REAL(KIND=dp)                                      :: al, be
     771        90247 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: full_block
     772              : 
     773        90247 :       al = 1.0_dp; be = 0.0_dp; i0 = 1; j0 = 1; tr_a = .FALSE.
     774              : 
     775        25048 :       IF (PRESENT(alpha)) al = alpha
     776        90247 :       IF (PRESENT(beta)) be = beta
     777        90247 :       IF (PRESENT(start_row)) i0 = start_row
     778        90247 :       IF (PRESENT(start_col)) j0 = start_col
     779        90247 :       IF (PRESENT(transpose)) tr_a = transpose
     780        16231 :       IF (tr_a) THEN
     781        16139 :          nrow = SIZE(new_values, 2)
     782        16139 :          ncol = SIZE(new_values, 1)
     783              :       ELSE
     784        74108 :          nrow = SIZE(new_values, 1)
     785        74108 :          ncol = SIZE(new_values, 2)
     786              :       END IF
     787        90247 :       IF (PRESENT(n_rows)) nrow = n_rows
     788        90247 :       IF (PRESENT(n_cols)) ncol = n_cols
     789              : 
     790        90247 :       full_block => fm%local_data
     791              : 
     792              :       CALL cp_fm_get_info(matrix=fm, &
     793              :                           nrow_global=nrow_global, ncol_global=ncol_global, &
     794              :                           nrow_block=nrow_block, ncol_block=ncol_block, &
     795              :                           nrow_local=nrow_local, ncol_local=ncol_local, &
     796        90247 :                           row_indices=row_indices, col_indices=col_indices)
     797              : 
     798        90247 :       IF (al == 1.0_dp .AND. be == 0.0_dp) THEN
     799      1435414 :          DO j = 1, ncol_local
     800      1360767 :             this_col = col_indices(j) - j0 + 1
     801      1435414 :             IF (this_col >= 1 .AND. this_col <= ncol) THEN
     802       326615 :                IF (tr_a) THEN
     803         7603 :                   IF (i0 == 1 .AND. nrow_global == nrow) THEN
     804       162336 :                      DO i = 1, nrow_local
     805       162336 :                         full_block(i, j) = new_values(this_col, row_indices(i))
     806              :                      END DO
     807              :                   ELSE
     808        11840 :                      DO i = 1, nrow_local
     809        10968 :                         this_row = row_indices(i) - i0 + 1
     810        11840 :                         IF (this_row >= 1 .AND. this_row <= nrow) THEN
     811         1750 :                            full_block(i, j) = new_values(this_col, this_row)
     812              :                         END IF
     813              :                      END DO
     814              :                   END IF
     815              :                ELSE
     816       319012 :                   IF (i0 == 1 .AND. nrow_global == nrow) THEN
     817      6881523 :                      DO i = 1, nrow_local
     818      6881523 :                         full_block(i, j) = new_values(row_indices(i), this_col)
     819              :                      END DO
     820              :                   ELSE
     821         9070 :                      DO i = 1, nrow_local
     822         8154 :                         this_row = row_indices(i) - i0 + 1
     823         9070 :                         IF (this_row >= 1 .AND. this_row <= nrow) THEN
     824         2703 :                            full_block(i, j) = new_values(this_row, this_col)
     825              :                         END IF
     826              :                      END DO
     827              :                   END IF
     828              :                END IF
     829              :             END IF
     830              :          END DO
     831              :       ELSE
     832       873320 :          DO j = 1, ncol_local
     833       857720 :             this_col = col_indices(j) - j0 + 1
     834       873320 :             IF (this_col >= 1 .AND. this_col <= ncol) THEN
     835       857720 :                IF (tr_a) THEN
     836     88195316 :                   DO i = 1, nrow_local
     837     87373812 :                      this_row = row_indices(i) - i0 + 1
     838     88195316 :                      IF (this_row >= 1 .AND. this_row <= nrow) THEN
     839              :                         full_block(i, j) = al*new_values(this_col, this_row) + &
     840       410752 :                                            be*full_block(i, j)
     841              :                      END IF
     842              :                   END DO
     843              :                ELSE
     844       140364 :                   DO i = 1, nrow_local
     845       104148 :                      this_row = row_indices(i) - i0 + 1
     846       140364 :                      IF (this_row >= 1 .AND. this_row <= nrow) THEN
     847              :                         full_block(i, j) = al*new_values(this_row, this_col) + &
     848       104148 :                                            be*full_block(i, j)
     849              :                      END IF
     850              :                   END DO
     851              :                END IF
     852              :             END IF
     853              :          END DO
     854              :       END IF
     855              : 
     856        90247 :    END SUBROUTINE cp_fm_set_submatrix
     857              : 
     858              : ! **************************************************************************************************
     859              : !> \brief sets a submatrix of a full matrix to a given value
     860              : !>       fm(start_row:start_row+n_rows,start_col:start_col+n_cols) = value
     861              : !> \param fm the full to change
     862              : !> \param new_value  ...
     863              : !> \param start_row the starting row of matrix
     864              : !> \param start_col the starting col of matrix
     865              : !> \param n_rows the number of rows to change
     866              : !> \param n_cols the number of columns to change
     867              : !> \par History
     868              : !>      07.2002 created borrowing from Joost's blacs_replicated_copy [fawzi]
     869              : !>      12.2025 created from cp_fm_set_submatrix
     870              : !> \author JGH
     871              : ! **************************************************************************************************
     872      1048461 :    SUBROUTINE cp_fm_set_all_submatrix(fm, new_value, start_row, start_col, n_rows, n_cols)
     873              :       TYPE(cp_fm_type), INTENT(IN)                       :: fm
     874              :       REAL(KIND=dp), INTENT(in)                          :: new_value
     875              :       INTEGER, INTENT(in)                                :: start_row, start_col, n_rows, n_cols
     876              : 
     877              :       INTEGER                                            :: i, i0, j, j0, ncol_global, ncol_local, &
     878              :                                                             nrow_global, nrow_local, this_col, &
     879              :                                                             this_row
     880      1048461 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
     881      1048461 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: full_block
     882              : 
     883      1048461 :       full_block => fm%local_data
     884              : 
     885              :       CALL cp_fm_get_info(matrix=fm, &
     886              :                           nrow_global=nrow_global, ncol_global=ncol_global, &
     887              :                           nrow_local=nrow_local, ncol_local=ncol_local, &
     888      1048461 :                           row_indices=row_indices, col_indices=col_indices)
     889              : 
     890      1048461 :       i0 = start_row
     891      1048461 :       j0 = start_col
     892     20785456 :       DO j = 1, ncol_local
     893     19736995 :          this_col = col_indices(j) - j0 + 1
     894     20785456 :          IF (this_col >= 1 .AND. this_col <= n_cols) THEN
     895    695992060 :             DO i = 1, nrow_local
     896    676315659 :                this_row = row_indices(i) - i0 + 1
     897    695992060 :                IF (this_row >= 1 .AND. this_row <= n_rows) THEN
     898    669248605 :                   full_block(i, j) = new_value
     899              :                END IF
     900              :             END DO
     901              :          END IF
     902              :       END DO
     903              : 
     904      1048461 :    END SUBROUTINE cp_fm_set_all_submatrix
     905              : 
     906              : ! **************************************************************************************************
     907              : !> \brief gets a submatrix of a full matrix
     908              : !>       op(target_m)(1:n_rows,1:n_cols)
     909              : !>       =fm(start_row:start_row+n_rows,start_col:start_col+n_cols)
     910              : !>      target_m is replicated on all cpus
     911              : !>      using this call is expensive
     912              : !> \param fm the full you want to get the info from
     913              : !> \param target_m a replicated full matrix that will contain the result
     914              : !> \param start_row the starting row of b_matrix (defaults to 1)
     915              : !> \param start_col the starting col of b_matrix (defaults to 1)
     916              : !> \param n_rows the number of row to change in b (defaults to
     917              : !>        size(op(new_values),1))
     918              : !> \param n_cols the number of columns to change in b (defaults to
     919              : !>        size(op(new_values),2))
     920              : !> \param transpose if target_m should be transposed: if true
     921              : !>        op(target_m)=target_m^T, else op(target_m)=target_m
     922              : !>        (defaults to false)
     923              : !> \par History
     924              : !>      07.2002 created borrowing from Joost's blacs_replicated_copy [fawzi]
     925              : !> \author Fawzi Mohamed
     926              : !> \note
     927              : !>      optimized for full column updates. Zeros out a little too much
     928              : !>      of target_m
     929              : !>      the target_m is replicated and valid on all cpus
     930              : ! **************************************************************************************************
     931        88540 :    SUBROUTINE cp_fm_get_submatrix(fm, target_m, start_row, &
     932              :                                   start_col, n_rows, n_cols, transpose)
     933              :       TYPE(cp_fm_type), INTENT(IN)                       :: fm
     934              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(out)        :: target_m
     935              :       INTEGER, INTENT(in), OPTIONAL                      :: start_row, start_col, n_rows, n_cols
     936              :       LOGICAL, INTENT(in), OPTIONAL                      :: transpose
     937              : 
     938              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_get_submatrix'
     939              : 
     940              :       INTEGER                                            :: handle, i, i0, j, j0, ncol, ncol_global, &
     941              :                                                             ncol_local, nrow, nrow_global, &
     942              :                                                             nrow_local, this_col, this_row
     943        88540 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
     944              :       LOGICAL                                            :: tr_a
     945        88540 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: full_block
     946              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     947              : 
     948        88540 :       CALL timeset(routineN, handle)
     949              : 
     950        88540 :       i0 = 1; j0 = 1; tr_a = .FALSE.
     951              : 
     952        88540 :       IF (PRESENT(start_row)) i0 = start_row
     953        88540 :       IF (PRESENT(start_col)) j0 = start_col
     954        88540 :       IF (PRESENT(transpose)) tr_a = transpose
     955         6152 :       IF (tr_a) THEN
     956         2206 :          nrow = SIZE(target_m, 2)
     957         2206 :          ncol = SIZE(target_m, 1)
     958              :       ELSE
     959        86334 :          nrow = SIZE(target_m, 1)
     960        86334 :          ncol = SIZE(target_m, 2)
     961              :       END IF
     962        88540 :       IF (PRESENT(n_rows)) nrow = n_rows
     963        88540 :       IF (PRESENT(n_cols)) ncol = n_cols
     964              : 
     965        88540 :       para_env => fm%matrix_struct%para_env
     966              : 
     967        88540 :       full_block => fm%local_data
     968              : #if defined(__parallel)
     969              :       ! zero-out whole target_m
     970        88540 :       IF (SIZE(target_m, 1)*SIZE(target_m, 2) /= 0) THEN
     971       207522 :          CALL dcopy(SIZE(target_m, 1)*SIZE(target_m, 2), [0.0_dp], 0, target_m, 1)
     972              :       END IF
     973              : #endif
     974              : 
     975              :       CALL cp_fm_get_info(matrix=fm, &
     976              :                           nrow_global=nrow_global, ncol_global=ncol_global, &
     977              :                           nrow_local=nrow_local, ncol_local=ncol_local, &
     978        88540 :                           row_indices=row_indices, col_indices=col_indices)
     979              : 
     980       576062 :       DO j = 1, ncol_local
     981       487522 :          this_col = col_indices(j) - j0 + 1
     982       576062 :          IF (this_col >= 1 .AND. this_col <= ncol) THEN
     983       349318 :             IF (tr_a) THEN
     984         2206 :                IF (i0 == 1 .AND. nrow_global == nrow) THEN
     985        80516 :                   DO i = 1, nrow_local
     986        80516 :                      target_m(this_col, row_indices(i)) = full_block(i, j)
     987              :                   END DO
     988              :                ELSE
     989            0 :                   DO i = 1, nrow_local
     990            0 :                      this_row = row_indices(i) - i0 + 1
     991            0 :                      IF (this_row >= 1 .AND. this_row <= nrow) THEN
     992            0 :                         target_m(this_col, this_row) = full_block(i, j)
     993              :                      END IF
     994              :                   END DO
     995              :                END IF
     996              :             ELSE
     997       347112 :                IF (i0 == 1 .AND. nrow_global == nrow) THEN
     998      6460657 :                   DO i = 1, nrow_local
     999      6460657 :                      target_m(row_indices(i), this_col) = full_block(i, j)
    1000              :                   END DO
    1001              :                ELSE
    1002      1596922 :                   DO i = 1, nrow_local
    1003      1551390 :                      this_row = row_indices(i) - i0 + 1
    1004      1596922 :                      IF (this_row >= 1 .AND. this_row <= nrow) THEN
    1005        55532 :                         target_m(this_row, this_col) = full_block(i, j)
    1006              :                      END IF
    1007              :                   END DO
    1008              :                END IF
    1009              :             END IF
    1010              :          END IF
    1011              :       END DO
    1012              : 
    1013     17754288 :       CALL para_env%sum(target_m)
    1014              : 
    1015        88540 :       CALL timestop(handle)
    1016              : 
    1017        88540 :    END SUBROUTINE cp_fm_get_submatrix
    1018              : 
    1019              : ! **************************************************************************************************
    1020              : !> \brief returns all kind of information about the full matrix
    1021              : !> \param matrix ...
    1022              : !> \param name ...
    1023              : !> \param nrow_global ...
    1024              : !> \param ncol_global ...
    1025              : !> \param nrow_block ...
    1026              : !> \param ncol_block ...
    1027              : !> \param nrow_local ...
    1028              : !> \param ncol_local ...
    1029              : !> \param row_indices ...
    1030              : !> \param col_indices ...
    1031              : !> \param local_data ...
    1032              : !> \param context ...
    1033              : !> \param nrow_locals ...
    1034              : !> \param ncol_locals ...
    1035              : !> \param matrix_struct ...
    1036              : !> \param para_env ...
    1037              : !> \note
    1038              : !>       see also cp_fm_struct for explanation
    1039              : !>       - nrow_local, ncol_local, row_indices, col_indices, local_data are hooks for efficient
    1040              : !>         access to the local blacs block
    1041              : ! **************************************************************************************************
    1042      6778218 :    SUBROUTINE cp_fm_get_info(matrix, name, nrow_global, ncol_global, &
    1043              :                              nrow_block, ncol_block, nrow_local, ncol_local, &
    1044              :                              row_indices, col_indices, local_data, context, &
    1045              :                              nrow_locals, ncol_locals, matrix_struct, para_env)
    1046              : 
    1047              :       TYPE(cp_fm_type), INTENT(IN)                       :: matrix
    1048              :       CHARACTER(LEN=*), INTENT(OUT), OPTIONAL            :: name
    1049              :       INTEGER, INTENT(OUT), OPTIONAL                     :: nrow_global, ncol_global, nrow_block, &
    1050              :                                                             ncol_block, nrow_local, ncol_local
    1051              :       INTEGER, DIMENSION(:), OPTIONAL, POINTER           :: row_indices, col_indices
    1052              :       REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
    1053              :          OPTIONAL, POINTER                               :: local_data
    1054              :       TYPE(cp_blacs_env_type), OPTIONAL, POINTER         :: context
    1055              :       INTEGER, DIMENSION(:), OPTIONAL, POINTER           :: nrow_locals, ncol_locals
    1056              :       TYPE(cp_fm_struct_type), OPTIONAL, POINTER         :: matrix_struct
    1057              :       TYPE(mp_para_env_type), OPTIONAL, POINTER          :: para_env
    1058              : 
    1059            8 :       IF (PRESENT(name)) name = matrix%name
    1060      6778218 :       IF (PRESENT(matrix_struct)) matrix_struct => matrix%matrix_struct
    1061      6778218 :       IF (PRESENT(local_data)) local_data => matrix%local_data ! not hiding things anymore :-(
    1062              : 
    1063              :       CALL cp_fm_struct_get(matrix%matrix_struct, nrow_local=nrow_local, &
    1064              :                             ncol_local=ncol_local, nrow_global=nrow_global, &
    1065              :                             ncol_global=ncol_global, nrow_block=nrow_block, &
    1066              :                             ncol_block=ncol_block, row_indices=row_indices, &
    1067              :                             col_indices=col_indices, nrow_locals=nrow_locals, &
    1068      6778218 :                             ncol_locals=ncol_locals, context=context, para_env=para_env)
    1069              : 
    1070      6778218 :    END SUBROUTINE cp_fm_get_info
    1071              : 
    1072              : ! **************************************************************************************************
    1073              : !> \brief Write nicely formatted info about the FM to the given I/O unit (including the underlying FM struct)
    1074              : !> \param matrix a cp_fm_type instance
    1075              : !> \param io_unit the I/O unit to use for writing
    1076              : ! **************************************************************************************************
    1077            3 :    SUBROUTINE cp_fm_write_info(matrix, io_unit)
    1078              :       TYPE(cp_fm_type), INTENT(IN)                       :: matrix
    1079              :       INTEGER, INTENT(IN)                                :: io_unit
    1080              : 
    1081            3 :       WRITE (io_unit, '(/,A,A12)') "CP_FM | Name:                           ", matrix%name
    1082            3 :       CALL cp_fm_struct_write_info(matrix%matrix_struct, io_unit)
    1083            3 :    END SUBROUTINE cp_fm_write_info
    1084              : 
    1085              : ! **************************************************************************************************
    1086              : !> \brief find the maximum absolute value of the matrix element
    1087              : !>      maxval(abs(matrix))
    1088              : !> \param matrix ...
    1089              : !> \param a_max ...
    1090              : !> \param ir_max ...
    1091              : !> \param ic_max ...
    1092              : ! **************************************************************************************************
    1093       131817 :    SUBROUTINE cp_fm_maxabsval(matrix, a_max, ir_max, ic_max)
    1094              :       TYPE(cp_fm_type), INTENT(IN)                       :: matrix
    1095              :       REAL(KIND=dp), INTENT(OUT)                         :: a_max
    1096              :       INTEGER, INTENT(OUT), OPTIONAL                     :: ir_max, ic_max
    1097              : 
    1098              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'cp_fm_maxabsval'
    1099              : 
    1100              :       INTEGER                                            :: handle, i, ic_max_local, ir_max_local, &
    1101              :                                                             j, mepos, ncol_local, nrow_local, &
    1102              :                                                             num_pe
    1103       131817 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: ic_max_vec, ir_max_vec
    1104       131817 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
    1105              :       REAL(dp)                                           :: my_max
    1106       131817 :       REAL(dp), ALLOCATABLE, DIMENSION(:)                :: a_max_vec
    1107       131817 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: my_block
    1108              : 
    1109       131817 :       CALL timeset(routineN, handle)
    1110              : 
    1111       131817 :       my_block => matrix%local_data
    1112              : 
    1113              :       CALL cp_fm_get_info(matrix, nrow_local=nrow_local, ncol_local=ncol_local, &
    1114       131817 :                           row_indices=row_indices, col_indices=col_indices)
    1115              : 
    1116     72817993 :       a_max = MAXVAL(ABS(my_block(1:nrow_local, 1:ncol_local)))
    1117              : 
    1118       131817 :       IF (PRESENT(ir_max)) THEN
    1119            0 :          num_pe = matrix%matrix_struct%para_env%num_pe
    1120            0 :          mepos = matrix%matrix_struct%para_env%mepos
    1121            0 :          ALLOCATE (ir_max_vec(0:num_pe - 1))
    1122            0 :          ir_max_vec(0:num_pe - 1) = 0
    1123            0 :          ALLOCATE (ic_max_vec(0:num_pe - 1))
    1124            0 :          ic_max_vec(0:num_pe - 1) = 0
    1125            0 :          ALLOCATE (a_max_vec(0:num_pe - 1))
    1126            0 :          a_max_vec(0:num_pe - 1) = 0.0_dp
    1127            0 :          my_max = 0.0_dp
    1128              : 
    1129            0 :          IF ((ncol_local > 0) .AND. (nrow_local > 0)) THEN
    1130            0 :             DO i = 1, ncol_local
    1131            0 :                DO j = 1, nrow_local
    1132            0 :                   IF (ABS(my_block(j, i)) > my_max) THEN
    1133            0 :                      my_max = my_block(j, i)
    1134            0 :                      ir_max_local = j
    1135            0 :                      ic_max_local = i
    1136              :                   END IF
    1137              :                END DO
    1138              :             END DO
    1139              : 
    1140            0 :             a_max_vec(mepos) = my_max
    1141            0 :             ir_max_vec(mepos) = row_indices(ir_max_local)
    1142            0 :             ic_max_vec(mepos) = col_indices(ic_max_local)
    1143              : 
    1144              :          END IF
    1145              : 
    1146            0 :          CALL matrix%matrix_struct%para_env%sum(a_max_vec)
    1147            0 :          CALL matrix%matrix_struct%para_env%sum(ir_max_vec)
    1148            0 :          CALL matrix%matrix_struct%para_env%sum(ic_max_vec)
    1149              : 
    1150            0 :          my_max = 0.0_dp
    1151            0 :          DO i = 0, num_pe - 1
    1152            0 :             IF (a_max_vec(i) > my_max) THEN
    1153            0 :                ir_max = ir_max_vec(i)
    1154            0 :                ic_max = ic_max_vec(i)
    1155              :             END IF
    1156              :          END DO
    1157              : 
    1158            0 :          DEALLOCATE (ir_max_vec, ic_max_vec, a_max_vec)
    1159            0 :          CPASSERT(ic_max > 0)
    1160            0 :          CPASSERT(ir_max > 0)
    1161              : 
    1162              :       END IF
    1163              : 
    1164       131817 :       CALL matrix%matrix_struct%para_env%max(a_max)
    1165              : 
    1166       131817 :       CALL timestop(handle)
    1167              : 
    1168       263634 :    END SUBROUTINE cp_fm_maxabsval
    1169              : 
    1170              : ! **************************************************************************************************
    1171              : !> \brief find the maximum over the rows of the sum of the absolute values of the elements of a given row
    1172              : !>      = || A ||_infinity
    1173              : !> \param matrix ...
    1174              : !> \param a_max ...
    1175              : !> \note
    1176              : !>      for a real symmetric matrix it holds that || A ||_2 = |lambda_max| < || A ||_infinity
    1177              : !>      Hence this can be used to estimate an upper bound for the eigenvalues of a matrix
    1178              : !>      http://mathworld.wolfram.com/MatrixNorm.html
    1179              : !>      (but the bound is not so tight in the general case)
    1180              : ! **************************************************************************************************
    1181         4984 :    SUBROUTINE cp_fm_maxabsrownorm(matrix, a_max)
    1182              :       TYPE(cp_fm_type), INTENT(IN)                       :: matrix
    1183              :       REAL(KIND=dp), INTENT(OUT)                         :: a_max
    1184              : 
    1185              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_maxabsrownorm'
    1186              : 
    1187              :       INTEGER                                            :: handle, i, j, ncol_local, nrow_global, &
    1188              :                                                             nrow_local
    1189         4984 :       INTEGER, DIMENSION(:), POINTER                     :: row_indices
    1190              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: values
    1191         4984 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: my_block
    1192              : 
    1193         4984 :       CALL timeset(routineN, handle)
    1194              : 
    1195         4984 :       my_block => matrix%local_data
    1196              : 
    1197              :       CALL cp_fm_get_info(matrix, row_indices=row_indices, nrow_global=nrow_global, &
    1198         4984 :                           nrow_local=nrow_local, ncol_local=ncol_local)
    1199              : 
    1200              :       ! the efficiency could be improved by making use of the row-col distribution of scalapack
    1201        14952 :       ALLOCATE (values(nrow_global))
    1202         4984 :       values = 0.0_dp
    1203        65802 :       DO j = 1, ncol_local
    1204       534993 :          DO i = 1, nrow_local
    1205       530009 :             values(row_indices(i)) = values(row_indices(i)) + ABS(my_block(i, j))
    1206              :          END DO
    1207              :       END DO
    1208         4984 :       CALL matrix%matrix_struct%para_env%sum(values)
    1209        65802 :       a_max = MAXVAL(values)
    1210         4984 :       DEALLOCATE (values)
    1211              : 
    1212         4984 :       CALL timestop(handle)
    1213         4984 :    END SUBROUTINE cp_fm_maxabsrownorm
    1214              : 
    1215              : ! **************************************************************************************************
    1216              : !> \brief find the norm of each column norm_{j}= sqrt( \sum_{i} A_{ij}*A_{ij} )
    1217              : !> \param matrix ...
    1218              : !> \param norm_array ...
    1219              : !> \param squared return squared norms instead (default: false)
    1220              : ! **************************************************************************************************
    1221         5296 :    SUBROUTINE cp_fm_vectorsnorm(matrix, norm_array, squared)
    1222              :       TYPE(cp_fm_type), INTENT(IN)                       :: matrix
    1223              :       REAL(KIND=dp), DIMENSION(:), INTENT(OUT)           :: norm_array
    1224              :       LOGICAL, INTENT(IN), OPTIONAL                      :: squared
    1225              : 
    1226              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'cp_fm_vectorsnorm'
    1227              : 
    1228              :       INTEGER                                            :: handle, i, j, ncol_local, nrow_local
    1229              : 
    1230              : !$    INTEGER                                            :: nworkers
    1231         5296 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices
    1232              :       LOGICAL                                            :: take_sqrt
    1233              :       REAL(KIND=dp)                                      :: norm_sq
    1234              :       REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
    1235         5296 :          POINTER                                         :: local_data
    1236              : 
    1237         5296 :       CALL timeset(routineN, handle)
    1238              : 
    1239         5296 :       local_data => matrix%local_data
    1240              :       CALL cp_fm_get_info(matrix, col_indices=col_indices, nrow_local=nrow_local, &
    1241         5296 :                           ncol_local=ncol_local)
    1242              : 
    1243        56430 :       norm_array = 0.0_dp
    1244              :       ! Own whole columns and keep at least 8192 elements per worker to amortize team startup.
    1245         5296 : !$    nworkers = MAX(1, MIN(omp_get_max_threads(), ncol_local))
    1246         5296 : !$    nworkers = MAX(1, MIN(nworkers, (nrow_local*ncol_local + 1)/8192))
    1247              : !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(STATIC) &
    1248              : !$OMP SHARED(local_data, col_indices, nrow_local, ncol_local, norm_array) PRIVATE(i, j, norm_sq) &
    1249         5296 : !$OMP NUM_THREADS(nworkers) IF(nworkers > 1)
    1250              :       DO j = 1, ncol_local
    1251              :          norm_sq = 0.0_dp
    1252              : !$OMP SIMD REDUCTION(+:norm_sq)
    1253              :          DO i = 1, nrow_local
    1254              :             norm_sq = norm_sq + local_data(i, j)**2
    1255              :          END DO
    1256              : !$OMP END SIMD
    1257              :          norm_array(col_indices(j)) = norm_sq
    1258              :       END DO
    1259              : !$OMP END PARALLEL DO
    1260       107564 :       CALL matrix%matrix_struct%para_env%sum(norm_array)
    1261         5296 :       take_sqrt = .TRUE.
    1262         5296 :       IF (PRESENT(squared)) take_sqrt = .NOT. squared
    1263        38044 :       IF (take_sqrt) norm_array = SQRT(norm_array)
    1264              : 
    1265         5296 :       CALL timestop(handle)
    1266         5296 :    END SUBROUTINE cp_fm_vectorsnorm
    1267              : 
    1268              : ! **************************************************************************************************
    1269              : !> \brief summing up all the elements along the matrix's i-th index
    1270              : !>        \f$ \mathrm{sum}_{j} = \sum_{i} A_{ij} \f$
    1271              : !>        or
    1272              : !>        \f$ \mathrm{sum}_{i} = \sum_{j} A_{ij} \f$
    1273              : !> \param matrix     an input matrix A
    1274              : !> \param sum_array  sums of elements in each column/row
    1275              : !> \param dir ...
    1276              : !> \note  forked from cp_fm_vectorsnorm() to be used with
    1277              : !>        the maximum overlap method
    1278              : !>        added row variation
    1279              : ! **************************************************************************************************
    1280        11348 :    SUBROUTINE cp_fm_vectorssum(matrix, sum_array, dir)
    1281              :       TYPE(cp_fm_type), INTENT(IN)                       :: matrix
    1282              :       REAL(KIND=dp), DIMENSION(:), INTENT(OUT)           :: sum_array
    1283              :       CHARACTER(LEN=1), INTENT(IN), OPTIONAL             :: dir
    1284              : 
    1285              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'cp_fm_vectorssum'
    1286              : 
    1287              :       INTEGER                                            :: handle, i, j, ncol_local, nrow_local
    1288        11348 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
    1289              :       LOGICAL                                            :: docol
    1290        11348 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: my_block
    1291              : 
    1292        11348 :       CALL timeset(routineN, handle)
    1293              : 
    1294        11348 :       IF (PRESENT(dir)) THEN
    1295        11308 :          IF (dir == 'c' .OR. dir == 'C') THEN
    1296              :             docol = .TRUE.
    1297              :          ELSE IF (dir == 'r' .OR. dir == 'R') THEN
    1298              :             docol = .FALSE.
    1299              :          ELSE
    1300            0 :             CPABORT('Wrong argument DIR')
    1301              :          END IF
    1302              :       ELSE
    1303              :          docol = .TRUE.
    1304              :       END IF
    1305              : 
    1306        11348 :       my_block => matrix%local_data
    1307              : 
    1308              :       CALL cp_fm_get_info(matrix, col_indices=col_indices, row_indices=row_indices, &
    1309        11348 :                           nrow_local=nrow_local, ncol_local=ncol_local)
    1310              : 
    1311              :       ! the efficiency could be improved by making use of the row-col distribution of scalapack
    1312       281172 :       sum_array(:) = 0.0_dp
    1313        11348 :       IF (docol) THEN
    1314          448 :       DO j = 1, ncol_local
    1315         3628 :          DO i = 1, nrow_local
    1316         3588 :             sum_array(col_indices(j)) = sum_array(col_indices(j)) + my_block(i, j)
    1317              :          END DO
    1318              :       END DO
    1319              :       ELSE
    1320       100972 :       DO j = 1, ncol_local
    1321      6517440 :          DO i = 1, nrow_local
    1322      6506132 :             sum_array(row_indices(i)) = sum_array(row_indices(i)) + my_block(i, j)
    1323              :          END DO
    1324              :       END DO
    1325              :       END IF
    1326       550996 :       CALL matrix%matrix_struct%para_env%sum(sum_array)
    1327              : 
    1328        11348 :       CALL timestop(handle)
    1329        11348 :    END SUBROUTINE cp_fm_vectorssum
    1330              : 
    1331              : ! **************************************************************************************************
    1332              : !> \brief summing up all the elements of the matrix
    1333              : !>        \f$ result = \sum_{ij} fun(A_{ij}) \f$
    1334              : !> \param matrix     an input matrix A
    1335              : !> \param fun ...
    1336              : !> \return ...
    1337              : ! **************************************************************************************************
    1338         3234 :    FUNCTION cp_fm_sumall(matrix, fun) RESULT(res)
    1339              :       TYPE(cp_fm_type), INTENT(IN)                       :: matrix
    1340              :       CHARACTER(LEN=3), INTENT(IN), OPTIONAL             :: fun
    1341              :       REAL(KIND=dp)                                      :: res
    1342              : 
    1343              :       CHARACTER(LEN=3)                                   :: funval
    1344         3234 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: data
    1345              : 
    1346         3234 :       IF (PRESENT(fun)) THEN
    1347            0 :          funval = fun
    1348              :       ELSE
    1349         3234 :          funval = "ONE"
    1350              :       END IF
    1351              : 
    1352         3234 :       DATA => matrix%local_data
    1353              : 
    1354              :       SELECT CASE (funval)
    1355              :       CASE ("ONE", "one")
    1356       300613 :          res = SUM(DATA)
    1357              :       CASE ("TWO", "two")
    1358            0 :          res = SUM(DATA**2)
    1359              :       CASE ("ABS", "abs")
    1360            0 :          res = SUM(ABS(DATA))
    1361              :       CASE DEFAULT
    1362         3234 :          CPABORT("Unknown function requested")
    1363              :       END SELECT
    1364              : 
    1365         3234 :       CALL matrix%matrix_struct%para_env%sum(res)
    1366              : 
    1367         3234 :    END FUNCTION cp_fm_sumall
    1368              : 
    1369              : ! **************************************************************************************************
    1370              : !> \brief copy one identically sized matrix in the other
    1371              : !> \param source ...
    1372              : !> \param destination ...
    1373              : !> \note
    1374              : !>      see also cp_fm_to_fm_columns
    1375              : ! **************************************************************************************************
    1376       427551 :    SUBROUTINE cp_fm_to_fm_matrix(source, destination)
    1377              : 
    1378              :       TYPE(cp_fm_type), INTENT(IN)                       :: source, destination
    1379              : 
    1380              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_to_fm_matrix'
    1381              : 
    1382              :       INTEGER                                            :: handle, npcol, nprow
    1383              : 
    1384       427551 :       CALL timeset(routineN, handle)
    1385              : 
    1386       427551 :       nprow = source%matrix_struct%context%num_pe(1)
    1387       427551 :       npcol = source%matrix_struct%context%num_pe(2)
    1388              : 
    1389       427551 :       IF ((.NOT. cp2k_is_parallel) .OR. &
    1390              :           cp_fm_struct_equivalent(source%matrix_struct, &
    1391              :                                   destination%matrix_struct)) THEN
    1392       427551 :          IF (SIZE(source%local_data, 1) /= SIZE(destination%local_data, 1) .OR. &
    1393              :              SIZE(source%local_data, 2) /= SIZE(destination%local_data, 2)) THEN
    1394              :             CALL cp_abort(__LOCATION__, &
    1395              :                           "Cannot copy full matrix <"//TRIM(source%name)// &
    1396              :                           "> to full matrix <"//TRIM(destination%name)// &
    1397            0 :                           ">. The local_data blocks have different sizes.")
    1398              :          END IF
    1399              :          CALL dcopy(SIZE(source%local_data, 1)*SIZE(source%local_data, 2), &
    1400       427551 :                     source%local_data, 1, destination%local_data, 1)
    1401              :       ELSE
    1402            0 :          CPABORT("Data structures of source and target full matrix are not equivalent")
    1403              :       END IF
    1404              : 
    1405       427551 :       CALL timestop(handle)
    1406              : 
    1407       427551 :    END SUBROUTINE cp_fm_to_fm_matrix
    1408              : 
    1409              : ! **************************************************************************************************
    1410              : !> \brief copy just a subset of columns of a fm to a fm
    1411              : !> \param msource ...
    1412              : !> \param mtarget ...
    1413              : !> \param ncol ...
    1414              : !> \param source_start ...
    1415              : !> \param target_start ...
    1416              : ! **************************************************************************************************
    1417       186612 :    SUBROUTINE cp_fm_to_fm_columns(msource, mtarget, ncol, source_start, &
    1418              :                                   target_start)
    1419              : 
    1420              :       TYPE(cp_fm_type), INTENT(IN)          :: msource, mtarget
    1421              :       INTEGER, INTENT(IN)                      :: ncol
    1422              :       INTEGER, INTENT(IN), OPTIONAL            :: source_start, target_start
    1423              : 
    1424              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_to_fm_columns'
    1425              : 
    1426              :       INTEGER                                  :: handle, n, ss, ts
    1427              :       REAL(KIND=dp), DIMENSION(:, :), POINTER  :: a, b
    1428              : #if defined(__parallel)
    1429              :       INTEGER                                  :: i
    1430              :       INTEGER, DIMENSION(9)                    :: desca, descb
    1431              : #endif
    1432              : 
    1433       186612 :       CALL timeset(routineN, handle)
    1434              : 
    1435       186612 :       ss = 1
    1436       186612 :       ts = 1
    1437              : 
    1438       186612 :       IF (PRESENT(source_start)) ss = source_start
    1439       186612 :       IF (PRESENT(target_start)) ts = target_start
    1440              : 
    1441       186612 :       n = msource%matrix_struct%nrow_global
    1442              : 
    1443       186612 :       a => msource%local_data
    1444       186612 :       b => mtarget%local_data
    1445              : 
    1446              : #if defined(__parallel)
    1447      1866120 :       desca(:) = msource%matrix_struct%descriptor(:)
    1448      1866120 :       descb(:) = mtarget%matrix_struct%descriptor(:)
    1449       829734 :       DO i = 0, ncol - 1
    1450       829734 :          CALL pdcopy(n, a, 1, ss + i, desca, 1, b, 1, ts + i, descb, 1)
    1451              :       END DO
    1452              : #else
    1453              :       IF (ss <= SIZE(a, 2) .AND. ts <= SIZE(b, 2)) THEN
    1454              :          CALL dcopy(ncol*n, a(:, ss), 1, b(:, ts), 1)
    1455              :       END IF
    1456              : #endif
    1457              : 
    1458       186612 :       CALL timestop(handle)
    1459              : 
    1460       186612 :    END SUBROUTINE cp_fm_to_fm_columns
    1461              : 
    1462              : ! **************************************************************************************************
    1463              : !> \brief copy just a triangular matrix
    1464              : !> \param msource ...
    1465              : !> \param mtarget ...
    1466              : !> \param uplo ...
    1467              : ! **************************************************************************************************
    1468           58 :    SUBROUTINE cp_fm_to_fm_triangular(msource, mtarget, uplo)
    1469              : 
    1470              :       TYPE(cp_fm_type), INTENT(IN)             :: msource, mtarget
    1471              :       CHARACTER(LEN=1), OPTIONAL, INTENT(IN)   :: uplo
    1472              : 
    1473              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_to_fm_triangular'
    1474              : 
    1475              :       CHARACTER(LEN=1)                         :: myuplo
    1476              :       INTEGER                                  :: handle, ncol, nrow
    1477           58 :       REAL(KIND=dp), DIMENSION(:, :), POINTER  :: a, b
    1478              : #if defined(__parallel)
    1479              :       INTEGER, DIMENSION(9)                    :: desca, descb
    1480              : #endif
    1481              : 
    1482           58 :       CALL timeset(routineN, handle)
    1483              : 
    1484           58 :       myuplo = 'U'
    1485           58 :       IF (PRESENT(uplo)) myuplo = uplo
    1486              : 
    1487           58 :       nrow = msource%matrix_struct%nrow_global
    1488           58 :       ncol = msource%matrix_struct%ncol_global
    1489              : 
    1490           58 :       a => msource%local_data
    1491           58 :       b => mtarget%local_data
    1492              : 
    1493              : #if defined(__parallel)
    1494          580 :       desca(:) = msource%matrix_struct%descriptor(:)
    1495          580 :       descb(:) = mtarget%matrix_struct%descriptor(:)
    1496           58 :       CALL pdlacpy(myuplo, nrow, ncol, a(1, 1), 1, 1, desca, b(1, 1), 1, 1, descb)
    1497              : #else
    1498              :       CALL dlacpy(myuplo, nrow, ncol, a(1, 1), nrow, b(1, 1), nrow)
    1499              : #endif
    1500              : 
    1501           58 :       CALL timestop(handle)
    1502              : 
    1503           58 :    END SUBROUTINE cp_fm_to_fm_triangular
    1504              : 
    1505              : ! **************************************************************************************************
    1506              : !> \brief copy just a part ot the matrix
    1507              : !> \param msource ...
    1508              : !> \param mtarget ...
    1509              : !> \param nrow ...
    1510              : !> \param ncol ...
    1511              : !> \param s_firstrow ...
    1512              : !> \param s_firstcol ...
    1513              : !> \param t_firstrow ...
    1514              : !> \param t_firstcol ...
    1515              : ! **************************************************************************************************
    1516              : 
    1517        16846 :    SUBROUTINE cp_fm_to_fm_submat(msource, mtarget, nrow, ncol, s_firstrow, s_firstcol, t_firstrow, t_firstcol)
    1518              : 
    1519              :       TYPE(cp_fm_type), INTENT(IN)             :: msource, mtarget
    1520              :       INTEGER, INTENT(IN)                      :: nrow, ncol, s_firstrow, &
    1521              :                                                   s_firstcol, t_firstrow, &
    1522              :                                                   t_firstcol
    1523              : 
    1524              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_to_fm_submat'
    1525              : 
    1526              :       INTEGER                                  :: handle, i, na, nb, ss, ts
    1527              :       REAL(KIND=dp), DIMENSION(:, :), POINTER  :: a, b
    1528              : #if defined(__parallel)
    1529              :       INTEGER, DIMENSION(9)                    :: desca, descb
    1530              : #endif
    1531              : 
    1532        16846 :       CALL timeset(routineN, handle)
    1533              : 
    1534        16846 :       a => msource%local_data
    1535        16846 :       b => mtarget%local_data
    1536              : 
    1537        16846 :       na = msource%matrix_struct%nrow_global
    1538        16846 :       nb = mtarget%matrix_struct%nrow_global
    1539              : !    nrow must be <= na and nb
    1540        16846 :       IF (nrow > na) THEN
    1541            0 :          CPABORT("cannot copy because nrow > number of rows of source matrix")
    1542              :       END IF
    1543        16846 :       IF (nrow > nb) THEN
    1544            0 :          CPABORT("cannot copy because nrow > number of rows of target matrix")
    1545              :       END IF
    1546        16846 :       na = msource%matrix_struct%ncol_global
    1547        16846 :       nb = mtarget%matrix_struct%ncol_global
    1548              : !    ncol must be <= na_col and nb_col
    1549        16846 :       IF (ncol > na) THEN
    1550            0 :          CPABORT("cannot copy because nrow > number of rows of source matrix")
    1551              :       END IF
    1552        16846 :       IF (ncol > nb) THEN
    1553            0 :          CPABORT("cannot copy because nrow > number of rows of target matrix")
    1554              :       END IF
    1555              : 
    1556              : #if defined(__parallel)
    1557       168460 :       desca(:) = msource%matrix_struct%descriptor(:)
    1558       168460 :       descb(:) = mtarget%matrix_struct%descriptor(:)
    1559       211064 :       DO i = 0, ncol - 1
    1560       194218 :          ss = s_firstcol + i
    1561       194218 :          ts = t_firstcol + i
    1562       211064 :          CALL pdcopy(nrow, a, s_firstrow, ss, desca, 1, b, t_firstrow, ts, descb, 1)
    1563              :       END DO
    1564              : #else
    1565              :       DO i = 0, ncol - 1
    1566              :          ss = s_firstcol + i
    1567              :          ts = t_firstcol + i
    1568              :          CALL dcopy(nrow, a(s_firstrow:, ss), 1, b(t_firstrow:, ts), 1)
    1569              :       END DO
    1570              : #endif
    1571              : 
    1572        16846 :       CALL timestop(handle)
    1573        16846 :    END SUBROUTINE cp_fm_to_fm_submat
    1574              : 
    1575              : ! **************************************************************************************************
    1576              : !> \brief General copy of a fm matrix to another fm matrix.
    1577              : !>        Uses non-blocking MPI rather than ScaLAPACK.
    1578              : !>
    1579              : !> \param source          input fm matrix
    1580              : !> \param destination     output fm matrix
    1581              : !> \param para_env        parallel environment corresponding to the BLACS env that covers all parts
    1582              : !>                        of the input and output matrices
    1583              : !> \par History
    1584              : !>      31-Jan-2017 : Re-implemented using non-blocking MPI [IainB, MarkT]
    1585              : ! **************************************************************************************************
    1586         5184 :    SUBROUTINE cp_fm_copy_general(source, destination, para_env)
    1587              :       TYPE(cp_fm_type), INTENT(IN)                       :: source, destination
    1588              :       TYPE(mp_para_env_type), INTENT(IN)                 :: para_env
    1589              : 
    1590              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_copy_general'
    1591              : 
    1592              :       INTEGER                                            :: handle
    1593        46656 :       TYPE(copy_info_type)                               :: info
    1594              : 
    1595         5184 :       CALL timeset(routineN, handle)
    1596              : 
    1597         5184 :       CALL cp_fm_start_copy_general(source, destination, para_env, info)
    1598         5184 :       IF (ASSOCIATED(destination%matrix_struct)) THEN
    1599         5170 :          CALL cp_fm_finish_copy_general(destination, info)
    1600              :       END IF
    1601         5184 :       IF (ASSOCIATED(source%matrix_struct)) THEN
    1602         5023 :          CALL cp_fm_cleanup_copy_general(info)
    1603              :       END IF
    1604              : 
    1605         5184 :       CALL timestop(handle)
    1606         5184 :    END SUBROUTINE cp_fm_copy_general
    1607              : 
    1608              : ! **************************************************************************************************
    1609              : !> \brief Initiates the copy operation: get distribution data, post MPI isend and irecvs
    1610              : !> \param source input fm matrix
    1611              : !> \param destination output fm matrix
    1612              : !> \param para_env parallel environment corresponding to the BLACS env that covers all parts
    1613              : !>                        of the input and output matrices
    1614              : !> \param info all of the data that will be needed to complete the copy operation
    1615              : ! **************************************************************************************************
    1616       209400 :    SUBROUTINE cp_fm_start_copy_general(source, destination, para_env, info)
    1617              :       TYPE(cp_fm_type), INTENT(IN)                       :: source, destination
    1618              :       TYPE(mp_para_env_type), INTENT(IN)                 :: para_env
    1619              :       TYPE(copy_info_type), INTENT(OUT)                  :: info
    1620              : 
    1621              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_start_copy_general'
    1622              : 
    1623              :       INTEGER :: dest_p_i, dest_q_j, global_rank, global_size, handle, i, j, k, mpi_rank, &
    1624              :          ncol_block_dest, ncol_block_src, ncol_local_recv, ncol_local_send, ncols, &
    1625              :          nrow_block_dest, nrow_block_src, nrow_local_recv, nrow_local_send, nrows, p, q, &
    1626              :          recv_rank, recv_size, send_rank, send_size
    1627        20940 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: all_ranks, dest2global, dest_p, dest_q, &
    1628        41880 :                                                             recv_count, send_count, send_disp, &
    1629        20940 :                                                             source2global, src_p, src_q
    1630        20940 :       INTEGER, ALLOCATABLE, DIMENSION(:, :)              :: dest_blacs2mpi
    1631              :       INTEGER, DIMENSION(2)                              :: dest_block, dest_block_tmp, dest_num_pe, &
    1632              :                                                             src_block, src_block_tmp, src_num_pe
    1633        41880 :       INTEGER, DIMENSION(:), POINTER                     :: recv_col_indices, recv_row_indices, &
    1634        41880 :                                                             send_col_indices, send_row_indices
    1635              :       TYPE(cp_fm_struct_type), POINTER                   :: recv_dist, send_dist
    1636       293160 :       TYPE(mp_request_type), DIMENSION(6)                :: recv_req, send_req
    1637              : 
    1638        20940 :       CALL timeset(routineN, handle)
    1639              : 
    1640              :       IF (.NOT. cp2k_is_parallel) THEN
    1641              :          ! Just copy all of the matrix data into a 'send buffer', to be unpacked later
    1642              :          nrow_local_send = SIZE(source%local_data, 1)
    1643              :          ncol_local_send = SIZE(source%local_data, 2)
    1644              :          ALLOCATE (info%send_buf(nrow_local_send*ncol_local_send))
    1645              :          k = 0
    1646              :          DO j = 1, ncol_local_send
    1647              :             DO i = 1, nrow_local_send
    1648              :                k = k + 1
    1649              :                info%send_buf(k) = source%local_data(i, j)
    1650              :             END DO
    1651              :          END DO
    1652              :       ELSE
    1653        20940 :          NULLIFY (recv_dist, send_dist)
    1654        20940 :          NULLIFY (recv_col_indices, recv_row_indices, send_col_indices, send_row_indices)
    1655              : 
    1656              :          ! The 'global' communicator contains both the source and destination decompositions
    1657        20940 :          global_size = para_env%num_pe
    1658        20940 :          global_rank = para_env%mepos
    1659              : 
    1660              :          ! The source/send decomposition and destination/recv decompositions may only exist on
    1661              :          ! on a subset of the processes involved in the communication
    1662              :          ! Check if the source and/or destination arguments are .not. ASSOCIATED():
    1663              :          ! if so, skip the send / recv parts (since these processes do not participate in the sending/receiving distribution)
    1664        20940 :          IF (ASSOCIATED(destination%matrix_struct)) THEN
    1665        16140 :             recv_dist => destination%matrix_struct
    1666        16140 :             recv_rank = recv_dist%para_env%mepos
    1667              :          ELSE
    1668         4800 :             recv_rank = mp_proc_null
    1669              :          END IF
    1670              : 
    1671        20940 :          IF (ASSOCIATED(source%matrix_struct)) THEN
    1672        18187 :             send_dist => source%matrix_struct
    1673        18187 :             send_rank = send_dist%para_env%mepos
    1674              :          ELSE
    1675         2753 :             send_rank = mp_proc_null
    1676              :          END IF
    1677              : 
    1678              :          ! Map the rank in the source/dest communicator to the global rank
    1679        62820 :          ALLOCATE (all_ranks(0:global_size - 1))
    1680              : 
    1681        20940 :          CALL para_env%allgather(send_rank, all_ranks)
    1682        20940 :          IF (ASSOCIATED(recv_dist)) THEN
    1683        80700 :             ALLOCATE (source2global(0:COUNT(all_ranks /= mp_proc_null) - 1))
    1684        48420 :             DO i = 0, global_size - 1
    1685        48420 :                IF (all_ranks(i) /= mp_proc_null) THEN
    1686        26774 :                   source2global(all_ranks(i)) = i
    1687              :                END IF
    1688              :             END DO
    1689              :          END IF
    1690              : 
    1691        20940 :          CALL para_env%allgather(recv_rank, all_ranks)
    1692        20940 :          IF (ASSOCIATED(send_dist)) THEN
    1693        90935 :             ALLOCATE (dest2global(0:COUNT(all_ranks /= mp_proc_null) - 1))
    1694        54561 :             DO i = 0, global_size - 1
    1695        54561 :                IF (all_ranks(i) /= mp_proc_null) THEN
    1696        26774 :                   dest2global(all_ranks(i)) = i
    1697              :                END IF
    1698              :             END DO
    1699              :          END IF
    1700        20940 :          DEALLOCATE (all_ranks)
    1701              : 
    1702              :          ! Some data from the two decompositions will be needed by all processes in the global group :
    1703              :          ! process grid shape, block size, and the BLACS-to-MPI mapping
    1704              : 
    1705              :          ! The global root process will receive the data (from the root process in each decomposition)
    1706       146580 :          send_req(:) = mp_request_null
    1707        20940 :          IF (global_rank == 0) THEN
    1708        73290 :             recv_req(:) = mp_request_null
    1709        10470 :             CALL para_env%irecv(src_block, mp_any_source, recv_req(1), tag=src_tag)
    1710        10470 :             CALL para_env%irecv(dest_block, mp_any_source, recv_req(2), tag=dest_tag)
    1711        10470 :             CALL para_env%irecv(src_num_pe, mp_any_source, recv_req(3), tag=src_tag)
    1712        10470 :             CALL para_env%irecv(dest_num_pe, mp_any_source, recv_req(4), tag=dest_tag)
    1713              :          END IF
    1714              : 
    1715        20940 :          IF (ASSOCIATED(send_dist)) THEN
    1716        18187 :             IF ((send_rank == 0)) THEN
    1717              :                ! need to use separate buffers here in case this is actually global rank 0
    1718        31410 :                src_block_tmp = [send_dist%nrow_block, send_dist%ncol_block]
    1719        10470 :                CALL para_env%isend(src_block_tmp, 0, send_req(1), tag=src_tag)
    1720        10470 :                CALL para_env%isend(send_dist%context%num_pe, 0, send_req(2), tag=src_tag)
    1721              :             END IF
    1722              :          END IF
    1723              : 
    1724        20940 :          IF (ASSOCIATED(recv_dist)) THEN
    1725        16140 :             IF ((recv_rank == 0)) THEN
    1726        31410 :                dest_block_tmp = [recv_dist%nrow_block, recv_dist%ncol_block]
    1727        10470 :                CALL para_env%isend(dest_block_tmp, 0, send_req(3), tag=dest_tag)
    1728        10470 :                CALL para_env%isend(recv_dist%context%num_pe, 0, send_req(4), tag=dest_tag)
    1729              :             END IF
    1730              :          END IF
    1731              : 
    1732        20940 :          IF (global_rank == 0) THEN
    1733        10470 :             CALL mp_waitall(recv_req(1:4))
    1734              :             ! Now we know the process decomposition, we can allocate the arrays to hold the blacs2mpi mapping
    1735            0 :             ALLOCATE (info%src_blacs2mpi(0:src_num_pe(1) - 1, 0:src_num_pe(2) - 1), &
    1736              :                       dest_blacs2mpi(0:dest_num_pe(1) - 1, 0:dest_num_pe(2) - 1) &
    1737        73290 :                       )
    1738        10470 :             CALL para_env%irecv(info%src_blacs2mpi, mp_any_source, recv_req(5), tag=src_tag)
    1739        10470 :             CALL para_env%irecv(dest_blacs2mpi, mp_any_source, recv_req(6), tag=dest_tag)
    1740              :          END IF
    1741              : 
    1742        20940 :          IF (ASSOCIATED(send_dist)) THEN
    1743        18187 :             IF ((send_rank == 0)) THEN
    1744        10470 :                CALL para_env%isend(send_dist%context%blacs2mpi(:, :), 0, send_req(5), tag=src_tag)
    1745              :             END IF
    1746              :          END IF
    1747              : 
    1748        20940 :          IF (ASSOCIATED(recv_dist)) THEN
    1749        16140 :             IF ((recv_rank == 0)) THEN
    1750        10470 :                CALL para_env%isend(recv_dist%context%blacs2mpi(:, :), 0, send_req(6), tag=dest_tag)
    1751              :             END IF
    1752              :          END IF
    1753              : 
    1754        20940 :          IF (global_rank == 0) THEN
    1755        10470 :             CALL mp_waitall(recv_req(5:6))
    1756              :          END IF
    1757              : 
    1758              :          ! Finally, broadcast the data to all processes in the global communicator
    1759        20940 :          CALL para_env%bcast(src_block, 0)
    1760        20940 :          CALL para_env%bcast(dest_block, 0)
    1761        20940 :          CALL para_env%bcast(src_num_pe, 0)
    1762        20940 :          CALL para_env%bcast(dest_num_pe, 0)
    1763        62820 :          info%src_num_pe(1:2) = src_num_pe(1:2)
    1764        62820 :          info%nblock_src(1:2) = src_block(1:2)
    1765        20940 :          IF (global_rank /= 0) THEN
    1766            0 :             ALLOCATE (info%src_blacs2mpi(0:src_num_pe(1) - 1, 0:src_num_pe(2) - 1), &
    1767            0 :                       dest_blacs2mpi(0:dest_num_pe(1) - 1, 0:dest_num_pe(2) - 1) &
    1768        73290 :                       )
    1769              :          END IF
    1770        20940 :          CALL para_env%bcast(info%src_blacs2mpi, 0)
    1771        20940 :          CALL para_env%bcast(dest_blacs2mpi, 0)
    1772              : 
    1773        20940 :          recv_size = dest_num_pe(1)*dest_num_pe(2)
    1774        20940 :          send_size = src_num_pe(1)*src_num_pe(2)
    1775        20940 :          info%send_size = send_size
    1776        20940 :          CALL mp_waitall(send_req(:))
    1777              : 
    1778              :          ! Setup is now complete, we can start the actual communication here.
    1779              :          ! The order implemented here is:
    1780              :          !  DEST_1
    1781              :          !      compute recv sizes
    1782              :          !      call irecv
    1783              :          !  SRC_1
    1784              :          !      compute send sizes
    1785              :          !      pack send buffers
    1786              :          !      call isend
    1787              :          !  DEST_2
    1788              :          !      wait for the recvs and unpack buffers (this part eventually will go into another
    1789              :          !         routine to allow comms to run concurrently)
    1790              :          !  SRC_2
    1791              :          !      wait for the sends
    1792              : 
    1793              :          ! DEST_1
    1794        20940 :          IF (ASSOCIATED(recv_dist)) THEN
    1795              :             CALL cp_fm_struct_get(recv_dist, row_indices=recv_row_indices, &
    1796              :                                   col_indices=recv_col_indices &
    1797        16140 :                                   )
    1798        16140 :             info%recv_col_indices => recv_col_indices
    1799        16140 :             info%recv_row_indices => recv_row_indices
    1800        16140 :             nrow_block_src = src_block(1)
    1801        16140 :             ncol_block_src = src_block(2)
    1802       107474 :             ALLOCATE (recv_count(0:send_size - 1), info%recv_disp(0:send_size - 1), info%recv_request(0:send_size - 1))
    1803              : 
    1804              :             ! Determine the recv counts, allocate the receive buffers, call mpi_irecv for all the non-zero sized receives
    1805        16140 :             nrow_local_recv = recv_dist%nrow_locals(recv_dist%context%mepos(1))
    1806        16140 :             ncol_local_recv = recv_dist%ncol_locals(recv_dist%context%mepos(2))
    1807        16140 :             info%nlocal_recv(1) = nrow_local_recv
    1808        16140 :             info%nlocal_recv(2) = ncol_local_recv
    1809              :             ! Initialise src_p, src_q arrays (sized using number of rows/cols in the receiving distribution)
    1810        80700 :             ALLOCATE (src_p(nrow_local_recv), src_q(ncol_local_recv))
    1811       233813 :             DO i = 1, nrow_local_recv
    1812              :                ! For each local row we will receive, we look up its global row (in recv_row_indices),
    1813              :                ! then work out which row block it comes from, and which process row that row block comes from.
    1814       233813 :                src_p(i) = MOD(((recv_row_indices(i) - 1)/nrow_block_src), src_num_pe(1))
    1815              :             END DO
    1816       370332 :             DO j = 1, ncol_local_recv
    1817              :                ! Similarly for the columns
    1818       370332 :                src_q(j) = MOD(((recv_col_indices(j) - 1)/ncol_block_src), src_num_pe(2))
    1819              :             END DO
    1820              :             ! src_p/q now contains the process row/column ID that will send data to that row/column
    1821              : 
    1822        32280 :             DO q = 0, src_num_pe(2) - 1
    1823       370332 :                ncols = COUNT(src_q == q)
    1824        59054 :                DO p = 0, src_num_pe(1) - 1
    1825       450730 :                   nrows = COUNT(src_p == p)
    1826              :                   ! Use the send_dist here as we are looking up the processes where the data comes from
    1827        42914 :                   recv_count(info%src_blacs2mpi(p, q)) = nrows*ncols
    1828              :                END DO
    1829              :             END DO
    1830        16140 :             DEALLOCATE (src_p, src_q)
    1831              : 
    1832              :             ! Use one long buffer (and displacements into that buffer)
    1833              :             !     this prevents the need for a rectangular array where not all elements will be populated
    1834        75194 :             ALLOCATE (info%recv_buf(SUM(recv_count(:))))
    1835        16140 :             info%recv_disp(0) = 0
    1836        26774 :             DO i = 1, send_size - 1
    1837        26774 :                info%recv_disp(i) = info%recv_disp(i - 1) + recv_count(i - 1)
    1838              :             END DO
    1839              : 
    1840              :             ! Issue receive calls on ranks which expect data
    1841        42914 :             DO k = 0, send_size - 1
    1842        42914 :                IF (recv_count(k) > 0) THEN
    1843              :                   CALL para_env%irecv(info%recv_buf(info%recv_disp(k) + 1:info%recv_disp(k) + recv_count(k)), &
    1844        20942 :                                       source2global(k), info%recv_request(k))
    1845              :                END IF
    1846              :             END DO
    1847        16140 :             DEALLOCATE (source2global)
    1848              :          END IF ! ASSOCIATED(recv_dist)
    1849              : 
    1850              :          ! SRC_1
    1851        20940 :          IF (ASSOCIATED(send_dist)) THEN
    1852              :             CALL cp_fm_struct_get(send_dist, row_indices=send_row_indices, &
    1853              :                                   col_indices=send_col_indices &
    1854        18187 :                                   )
    1855        18187 :             nrow_block_dest = dest_block(1)
    1856        18187 :             ncol_block_dest = dest_block(2)
    1857       117709 :             ALLOCATE (send_count(0:recv_size - 1), send_disp(0:recv_size - 1), info%send_request(0:recv_size - 1))
    1858              : 
    1859              :             ! Determine the send counts, allocate the send buffers
    1860        18187 :             nrow_local_send = send_dist%nrow_locals(send_dist%context%mepos(1))
    1861        18187 :             ncol_local_send = send_dist%ncol_locals(send_dist%context%mepos(2))
    1862              : 
    1863              :             ! Initialise dest_p, dest_q arrays (sized nrow_local, ncol_local)
    1864              :             !   i.e. number of rows,cols in the sending distribution
    1865        90935 :             ALLOCATE (dest_p(nrow_local_send), dest_q(ncol_local_send))
    1866              : 
    1867       235860 :             DO i = 1, nrow_local_send
    1868              :                ! Use the send_dist%row_indices() here (we are looping over the local rows we will send)
    1869       235860 :                dest_p(i) = MOD(((send_row_indices(i) - 1)/nrow_block_dest), dest_num_pe(1))
    1870              :             END DO
    1871       430451 :             DO j = 1, ncol_local_send
    1872       430451 :                dest_q(j) = MOD(((send_col_indices(j) - 1)/ncol_block_dest), dest_num_pe(2))
    1873              :             END DO
    1874              :             ! dest_p/q now contain the process row/column ID that will receive data from that row/column
    1875              : 
    1876        36374 :             DO q = 0, dest_num_pe(2) - 1
    1877       430451 :                ncols = COUNT(dest_q == q)
    1878        63148 :                DO p = 0, dest_num_pe(1) - 1
    1879       395740 :                   nrows = COUNT(dest_p == p)
    1880        44961 :                   send_count(dest_blacs2mpi(p, q)) = nrows*ncols
    1881              :                END DO
    1882              :             END DO
    1883        18187 :             DEALLOCATE (dest_p, dest_q)
    1884              : 
    1885              :             ! Allocate the send buffer using send_count -- and calculate the offset into the buffer for each process
    1886        81335 :             ALLOCATE (info%send_buf(SUM(send_count(:))))
    1887        18187 :             send_disp(0) = 0
    1888        26774 :             DO k = 1, recv_size - 1
    1889        26774 :                send_disp(k) = send_disp(k - 1) + send_count(k - 1)
    1890              :             END DO
    1891              : 
    1892              :             ! Loop over the smat, pack the send buffers
    1893        18187 :             send_count(:) = 0
    1894       430451 :             DO j = 1, ncol_local_send
    1895              :                ! Use send_col_indices and row_indices here, as we are looking up the global row/column number of local rows.
    1896       412264 :                dest_q_j = MOD(((send_col_indices(j) - 1)/ncol_block_dest), dest_num_pe(2))
    1897     12265263 :                DO i = 1, nrow_local_send
    1898     11834812 :                   dest_p_i = MOD(((send_row_indices(i) - 1)/nrow_block_dest), dest_num_pe(1))
    1899     11834812 :                   mpi_rank = dest_blacs2mpi(dest_p_i, dest_q_j)
    1900     11834812 :                   send_count(mpi_rank) = send_count(mpi_rank) + 1
    1901     12247076 :                   info%send_buf(send_disp(mpi_rank) + send_count(mpi_rank)) = source%local_data(i, j)
    1902              :                END DO
    1903              :             END DO
    1904              : 
    1905              :             ! For each non-zero send_count, call mpi_isend
    1906        44961 :             DO k = 0, recv_size - 1
    1907        44961 :                IF (send_count(k) > 0) THEN
    1908              :                   CALL para_env%isend(info%send_buf(send_disp(k) + 1:send_disp(k) + send_count(k)), &
    1909        20942 :                                       dest2global(k), info%send_request(k))
    1910              :                END IF
    1911              :             END DO
    1912        18187 :             DEALLOCATE (send_count, send_disp, dest2global)
    1913              :          END IF ! ASSOCIATED(send_dist)
    1914        20940 :          DEALLOCATE (dest_blacs2mpi)
    1915              : 
    1916              :       END IF !IF (.NOT. cp2k_is_parallel)
    1917              : 
    1918        20940 :       CALL timestop(handle)
    1919              : 
    1920        83760 :    END SUBROUTINE cp_fm_start_copy_general
    1921              : 
    1922              : ! **************************************************************************************************
    1923              : !> \brief Completes the copy operation: wait for comms, unpack, clean up MPI state
    1924              : !> \param destination output fm matrix
    1925              : !> \param info all of the data that will be needed to complete the copy operation
    1926              : ! **************************************************************************************************
    1927        16140 :    SUBROUTINE cp_fm_finish_copy_general(destination, info)
    1928              :       TYPE(cp_fm_type), INTENT(IN)                       :: destination
    1929              :       TYPE(copy_info_type), INTENT(INOUT)                :: info
    1930              : 
    1931              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_finish_copy_general'
    1932              : 
    1933              :       INTEGER                                            :: handle, i, j, k, mpi_rank, send_size, &
    1934              :                                                             src_p_i, src_q_j
    1935        16140 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: recv_count
    1936              :       INTEGER, DIMENSION(2)                              :: nblock_src, nlocal_recv, src_num_pe
    1937        16140 :       INTEGER, DIMENSION(:), POINTER                     :: recv_col_indices, recv_row_indices
    1938              : 
    1939        16140 :       CALL timeset(routineN, handle)
    1940              : 
    1941              :       IF (.NOT. cp2k_is_parallel) THEN
    1942              :          ! Now unpack the data from the 'send buffer'
    1943              :          k = 0
    1944              :          DO j = 1, SIZE(destination%local_data, 2)
    1945              :             DO i = 1, SIZE(destination%local_data, 1)
    1946              :                k = k + 1
    1947              :                destination%local_data(i, j) = info%send_buf(k)
    1948              :             END DO
    1949              :          END DO
    1950              :          DEALLOCATE (info%send_buf)
    1951              :       ELSE
    1952              :          ! Set up local variables ...
    1953        16140 :          send_size = info%send_size
    1954        48420 :          nlocal_recv(1:2) = info%nlocal_recv(:)
    1955        48420 :          nblock_src(1:2) = info%nblock_src(:)
    1956        48420 :          src_num_pe(1:2) = info%src_num_pe(:)
    1957        16140 :          recv_col_indices => info%recv_col_indices
    1958        16140 :          recv_row_indices => info%recv_row_indices
    1959              : 
    1960              :          ! ... use the local variables to do the work
    1961              :          ! DEST_2
    1962        16140 :          CALL mp_waitall(info%recv_request(:))
    1963        48420 :          ALLOCATE (recv_count(0:send_size - 1))
    1964              :          ! Loop over the rmat, filling it in with data from the recv buffers
    1965              :          ! (here the block sizes, num_pes refer to the distribution of the source matrix)
    1966        16140 :          recv_count(:) = 0
    1967       370332 :          DO j = 1, nlocal_recv(2)
    1968       354192 :             src_q_j = MOD(((recv_col_indices(j) - 1)/nblock_src(2)), src_num_pe(2))
    1969     12205144 :             DO i = 1, nlocal_recv(1)
    1970     11834812 :                src_p_i = MOD(((recv_row_indices(i) - 1)/nblock_src(1)), src_num_pe(1))
    1971     11834812 :                mpi_rank = info%src_blacs2mpi(src_p_i, src_q_j)
    1972     11834812 :                recv_count(mpi_rank) = recv_count(mpi_rank) + 1
    1973     12189004 :                destination%local_data(i, j) = info%recv_buf(info%recv_disp(mpi_rank) + recv_count(mpi_rank))
    1974              :             END DO
    1975              :          END DO
    1976        16140 :          DEALLOCATE (recv_count, info%recv_disp, info%recv_request, info%recv_buf, info%src_blacs2mpi)
    1977              :          ! Invalidate the stored state
    1978              :          NULLIFY (info%recv_col_indices, &
    1979        16140 :                   info%recv_row_indices)
    1980              : 
    1981              :       END IF
    1982              : 
    1983        16140 :       CALL timestop(handle)
    1984              : 
    1985        16140 :    END SUBROUTINE cp_fm_finish_copy_general
    1986              : 
    1987              : ! **************************************************************************************************
    1988              : !> \brief Completes the copy operation: wait for comms clean up MPI state
    1989              : !> \param info all of the data that will be needed to complete the copy operation
    1990              : ! **************************************************************************************************
    1991        18187 :    SUBROUTINE cp_fm_cleanup_copy_general(info)
    1992              :       TYPE(copy_info_type), INTENT(INOUT)                :: info
    1993              : 
    1994              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_cleanup_copy_general'
    1995              : 
    1996              :       INTEGER                                            :: handle
    1997              : 
    1998        18187 :       CALL timeset(routineN, handle)
    1999              : 
    2000              :       IF (.NOT. cp2k_is_parallel) THEN
    2001              :          ! Don't do anything - no MPI state for the serial case
    2002              :       ELSE
    2003              :          ! SRC_2
    2004              :          ! If this process is also in the destination decomposition, this deallocate
    2005              :          ! Was already done in cp_fm_finish_copy_general
    2006        18187 :          IF (ALLOCATED(info%src_blacs2mpi)) THEN
    2007         4800 :             DEALLOCATE (info%src_blacs2mpi)
    2008              :          END IF
    2009        18187 :          CALL mp_waitall(info%send_request)
    2010        18187 :          DEALLOCATE (info%send_request, info%send_buf)
    2011              : 
    2012              :       END IF
    2013              : 
    2014        18187 :       CALL timestop(handle)
    2015              : 
    2016        18187 :    END SUBROUTINE cp_fm_cleanup_copy_general
    2017              : 
    2018              : ! **************************************************************************************************
    2019              : !> \brief General copy of a submatrix of fm matrix to  a submatrix of another fm matrix.
    2020              : !>        The two matrices can have different contexts.
    2021              : !>
    2022              : !>        Summary of distribution routines for dense matrices
    2023              : !>        The following will copy A(iA:iA+M-1,jA:jA+N-1) to B(iB:iB+M-1,jB:jB+N-1):
    2024              : !>
    2025              : !>        call pdgemr2d(M,N,Aloc,iA,jA,descA,Bloc,iB,jB,descB,context)
    2026              : !>
    2027              : !>        A process that is not a part of the context of A should set descA(2)
    2028              : !>        to -1, and similarly for B.
    2029              : !>
    2030              : !> \param source          input fm matrix
    2031              : !> \param destination     output fm matrix
    2032              : !> \param nrows           number of rows of sub matrix to be copied
    2033              : !> \param ncols           number of cols of sub matrix to be copied
    2034              : !> \param s_firstrow      starting global row index of sub matrix in source
    2035              : !> \param s_firstcol      starting global col index of sub matrix in source
    2036              : !> \param d_firstrow      starting global row index of sub matrix in destination
    2037              : !> \param d_firstcol      starting global col index of sub matrix in destination
    2038              : !> \param global_context  process grid that covers all parts of either A or B.
    2039              : ! **************************************************************************************************
    2040        13789 :    SUBROUTINE cp_fm_to_fm_submat_general(source, &
    2041              :                                          destination, &
    2042              :                                          nrows, &
    2043              :                                          ncols, &
    2044              :                                          s_firstrow, &
    2045              :                                          s_firstcol, &
    2046              :                                          d_firstrow, &
    2047              :                                          d_firstcol, &
    2048              :                                          global_context)
    2049              : 
    2050              :       TYPE(cp_fm_type), INTENT(IN)                       :: source, destination
    2051              :       INTEGER, INTENT(IN)                                :: nrows, ncols, s_firstrow, s_firstcol, &
    2052              :                                                             d_firstrow, d_firstcol
    2053              : 
    2054              :       CLASS(cp_blacs_type), INTENT(IN)        :: global_context
    2055              : 
    2056              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_to_fm_submat_general'
    2057              : 
    2058              :       LOGICAL                                  :: debug
    2059              :       INTEGER                                  :: handle
    2060              : #if defined(__parallel)
    2061              :       INTEGER, DIMENSION(9)                    :: desca, descb
    2062              :       REAL(KIND=dp), DIMENSION(1, 1), TARGET   :: dummy
    2063        13789 :       REAL(KIND=dp), DIMENSION(:, :), POINTER  :: smat, dmat
    2064              : #endif
    2065              : 
    2066        13789 :       CALL timeset(routineN, handle)
    2067              : 
    2068        13789 :       debug = debug_this_module
    2069              : 
    2070              :       IF (.NOT. cp2k_is_parallel) THEN
    2071              :          CALL cp_fm_to_fm_submat(source, &
    2072              :                                  destination, &
    2073              :                                  nrows, &
    2074              :                                  ncols, &
    2075              :                                  s_firstrow, &
    2076              :                                  s_firstcol, &
    2077              :                                  d_firstrow, &
    2078              :                                  d_firstcol)
    2079              :       ELSE
    2080              : #ifdef __parallel
    2081              :          NULLIFY (smat, dmat)
    2082              :          ! check whether source is available on this process
    2083        13789 :          IF (ASSOCIATED(source%matrix_struct)) THEN
    2084       137890 :             desca = source%matrix_struct%descriptor
    2085        13789 :             IF (nrows > source%matrix_struct%nrow_global) THEN
    2086            0 :                CPABORT("nrows is greater than nrow_global of source")
    2087              :             END IF
    2088        13789 :             IF (ncols > source%matrix_struct%ncol_global) THEN
    2089            0 :                CPABORT("ncols is greater than ncol_global of source")
    2090              :             END IF
    2091        13789 :             smat => source%local_data
    2092              :          ELSE
    2093            0 :             desca = -1
    2094            0 :             smat => dummy
    2095              :          END IF
    2096              :          ! check destination is available on this process
    2097        13789 :          IF (ASSOCIATED(destination%matrix_struct)) THEN
    2098       137890 :             descb = destination%matrix_struct%descriptor
    2099        13789 :             IF (nrows > destination%matrix_struct%nrow_global) THEN
    2100            0 :                CPABORT("nrows is greater than nrow_global of destination")
    2101              :             END IF
    2102        13789 :             IF (ncols > destination%matrix_struct%ncol_global) THEN
    2103            0 :                CPABORT("ncols is greater than ncol_global of destination")
    2104              :             END IF
    2105        13789 :             dmat => destination%local_data
    2106              :          ELSE
    2107            0 :             descb = -1
    2108            0 :             dmat => dummy
    2109              :          END IF
    2110              :          ! do copy
    2111              : 
    2112              :          CALL pdgemr2d(nrows, &
    2113              :                        ncols, &
    2114              :                        smat, &
    2115              :                        s_firstrow, &
    2116              :                        s_firstcol, &
    2117              :                        desca, &
    2118              :                        dmat, &
    2119              :                        d_firstrow, &
    2120              :                        d_firstcol, &
    2121              :                        descb, &
    2122        13789 :                        global_context%get_handle())
    2123              : #else
    2124              :          MARK_USED(global_context)
    2125              :          CPABORT("this subroutine only supports SCALAPACK")
    2126              : #endif
    2127              :       END IF
    2128              : 
    2129        13789 :       CALL timestop(handle)
    2130              : 
    2131        13789 :    END SUBROUTINE cp_fm_to_fm_submat_general
    2132              : 
    2133              : ! **************************************************************************************************
    2134              : !> \brief ...
    2135              : !> \param matrix ...
    2136              : !> \param irow_global ...
    2137              : !> \param icol_global ...
    2138              : !> \param alpha ...
    2139              : ! **************************************************************************************************
    2140          240 :    SUBROUTINE cp_fm_add_to_element(matrix, irow_global, icol_global, alpha)
    2141              : 
    2142              :       ! Add alpha to the matrix element specified by the global indices
    2143              :       ! irow_global and icol_global
    2144              : 
    2145              :       ! - Creation (05.05.06,MK)
    2146              : 
    2147              :       TYPE(cp_fm_type), INTENT(IN)          :: matrix
    2148              :       INTEGER, INTENT(IN)                      :: irow_global, icol_global
    2149              :       REAL(KIND=dp), INTENT(IN)                :: alpha
    2150              : 
    2151              :       INTEGER                                  :: mypcol, myprow, npcol, nprow
    2152          240 :       REAL(KIND=dp), DIMENSION(:, :), POINTER  :: a
    2153              :       TYPE(cp_blacs_env_type), POINTER         :: context
    2154              : #if defined(__parallel)
    2155              :       INTEGER                                  :: icol_local, ipcol, iprow, &
    2156              :                                                   irow_local
    2157              :       INTEGER, DIMENSION(9)                    :: desca
    2158              : #endif
    2159              : 
    2160          240 :       context => matrix%matrix_struct%context
    2161              : 
    2162          240 :       myprow = context%mepos(1)
    2163          240 :       mypcol = context%mepos(2)
    2164              : 
    2165          240 :       nprow = context%num_pe(1)
    2166          240 :       npcol = context%num_pe(2)
    2167              : 
    2168          240 :       a => matrix%local_data
    2169              : 
    2170              : #if defined(__parallel)
    2171              : 
    2172         2400 :       desca(:) = matrix%matrix_struct%descriptor(:)
    2173              : 
    2174              :       CALL infog2l(irow_global, icol_global, desca, nprow, npcol, myprow, mypcol, &
    2175          240 :                    irow_local, icol_local, iprow, ipcol)
    2176              : 
    2177          240 :       IF ((iprow == myprow) .AND. (ipcol == mypcol)) THEN
    2178          120 :          a(irow_local, icol_local) = a(irow_local, icol_local) + alpha
    2179              :       END IF
    2180              : 
    2181              : #else
    2182              : 
    2183              :       a(irow_global, icol_global) = a(irow_global, icol_global) + alpha
    2184              : 
    2185              : #endif
    2186              : 
    2187          240 :    END SUBROUTINE cp_fm_add_to_element
    2188              : 
    2189              : ! **************************************************************************************************
    2190              : !> \brief ...
    2191              : !> \param fm ...
    2192              : !> \param unit ...
    2193              : ! **************************************************************************************************
    2194       142516 :    SUBROUTINE cp_fm_write_unformatted(fm, unit)
    2195              :       TYPE(cp_fm_type), INTENT(IN)             :: fm
    2196              :       INTEGER, INTENT(IN)                      :: unit
    2197              : 
    2198              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_write_unformatted'
    2199              : 
    2200              :       INTEGER                                  :: handle, j, max_block, &
    2201              :                                                   ncol_global, nrow_global
    2202              :       TYPE(mp_para_env_type), POINTER          :: para_env
    2203              : #if defined(__parallel)
    2204              :       INTEGER                                  :: i, i_block, icol_local, &
    2205              :                                                   in, info, ipcol, &
    2206              :                                                   iprow, irow_local, &
    2207              :                                                   mepos, &
    2208              :                                                   num_pe, rb, tag
    2209              :       INTEGER, DIMENSION(9)                    :: desc
    2210       142516 :       REAL(KIND=dp), DIMENSION(:), POINTER     :: vecbuf
    2211              :       REAL(KIND=dp), DIMENSION(:, :), POINTER  :: newdat
    2212              :       TYPE(cp_blacs_type) :: ictxt_loc
    2213              :       INTEGER, EXTERNAL                        :: numroc
    2214              : #endif
    2215              : 
    2216       142516 :       CALL timeset(routineN, handle)
    2217              :       CALL cp_fm_get_info(fm, nrow_global=nrow_global, ncol_global=ncol_global, ncol_block=max_block, &
    2218       142516 :                           para_env=para_env)
    2219              : 
    2220              : #if defined(__parallel)
    2221       142516 :       num_pe = para_env%num_pe
    2222       142516 :       mepos = para_env%mepos
    2223       142516 :       rb = nrow_global
    2224       142516 :       tag = 0
    2225              :       ! get a new context
    2226       142516 :       CALL ictxt_loc%gridinit(para_env, 'R', 1, num_pe)
    2227       142516 :       CALL descinit(desc, nrow_global, ncol_global, rb, max_block, 0, 0, ictxt_loc%get_handle(), nrow_global, info)
    2228       142516 :       CPASSERT(info == 0)
    2229              :       ASSOCIATE (nprow => ictxt_loc%num_pe(1), npcol => ictxt_loc%num_pe(2), &
    2230              :                  myprow => ictxt_loc%mepos(1), mypcol => ictxt_loc%mepos(2))
    2231       285032 :          in = numroc(ncol_global, max_block, mypcol, 0, npcol)
    2232              : 
    2233       570064 :          ALLOCATE (newdat(nrow_global, MAX(1, in)))
    2234              : 
    2235              :          ! do the actual scalapack to cols reordering
    2236              :          CALL pdgemr2d(nrow_global, ncol_global, fm%local_data, 1, 1, &
    2237              :                        fm%matrix_struct%descriptor, &
    2238       142516 :                        newdat, 1, 1, desc, ictxt_loc%get_handle())
    2239              : 
    2240       427548 :          ALLOCATE (vecbuf(nrow_global*max_block))
    2241     50236318 :          vecbuf = HUGE(1.0_dp) ! init for valgrind
    2242              : 
    2243       631229 :          DO i = 1, ncol_global, MAX(max_block, 1)
    2244       346197 :             i_block = MIN(max_block, ncol_global - i + 1)
    2245              :             CALL infog2l(1, i, desc, nprow, npcol, myprow, mypcol, &
    2246       346197 :                          irow_local, icol_local, iprow, ipcol)
    2247       346197 :             IF (ipcol == mypcol) THEN
    2248      1425434 :                DO j = 1, i_block
    2249    121104700 :                   vecbuf((j - 1)*nrow_global + 1:nrow_global*j) = newdat(:, icol_local + j - 1)
    2250              :                END DO
    2251              :             END IF
    2252              : 
    2253       346197 :             IF (ipcol == 0) THEN
    2254              :                ! do nothing
    2255              :             ELSE
    2256       137538 :                IF (ipcol == mypcol) THEN
    2257     27137783 :                   CALL para_env%send(vecbuf(:), 0, tag)
    2258              :                END IF
    2259       137538 :                IF (mypcol == 0) THEN
    2260     54206797 :                   CALL para_env%recv(vecbuf(:), ipcol, tag)
    2261              :                END IF
    2262              :             END IF
    2263              : 
    2264       488713 :             IF (unit > 0) THEN
    2265      1413339 :                DO j = 1, i_block
    2266     61147944 :                   WRITE (unit) vecbuf((j - 1)*nrow_global + 1:nrow_global*j)
    2267              :                END DO
    2268              :             END IF
    2269              : 
    2270              :          END DO
    2271              :       END ASSOCIATE
    2272       142516 :       DEALLOCATE (vecbuf)
    2273              : 
    2274       142516 :       CALL ictxt_loc%gridexit()
    2275              : 
    2276       142516 :       DEALLOCATE (newdat)
    2277              : 
    2278              : #else
    2279              : 
    2280              :       IF (unit > 0) THEN
    2281              :          DO j = 1, ncol_global
    2282              :             WRITE (unit) fm%local_data(:, j)
    2283              :          END DO
    2284              :       END IF
    2285              : 
    2286              : #endif
    2287       142516 :       CALL timestop(handle)
    2288              : 
    2289       712580 :    END SUBROUTINE cp_fm_write_unformatted
    2290              : 
    2291              : ! **************************************************************************************************
    2292              : !> \brief Write out a full matrix in plain text.
    2293              : !> \param fm the matrix to be outputted
    2294              : !> \param unit the unit number for I/O
    2295              : !> \param header optional header
    2296              : !> \param value_format ...
    2297              : ! **************************************************************************************************
    2298         3278 :    SUBROUTINE cp_fm_write_formatted(fm, unit, header, value_format)
    2299              :       TYPE(cp_fm_type), INTENT(IN)             :: fm
    2300              :       INTEGER, INTENT(IN)                      :: unit
    2301              :       CHARACTER(LEN=*), INTENT(IN), OPTIONAL   :: header, value_format
    2302              : 
    2303              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_write_formatted'
    2304              : 
    2305              :       CHARACTER(LEN=21)                        :: my_value_format
    2306              :       INTEGER                                  :: handle, i, j, max_block, &
    2307              :                                                   ncol_global, nrow_global
    2308              :       TYPE(mp_para_env_type), POINTER          :: para_env
    2309              : #if defined(__parallel)
    2310              :       INTEGER                                  :: i_block, icol_local, &
    2311              :                                                   in, info, ipcol, &
    2312              :                                                   iprow, irow_local, &
    2313              :                                                   mepos, num_pe, rb, tag, k, &
    2314              :                                                   icol, irow
    2315              :       INTEGER, DIMENSION(9)                    :: desc
    2316         3278 :       REAL(KIND=dp), DIMENSION(:), POINTER     :: vecbuf
    2317         3278 :       REAL(KIND=dp), DIMENSION(:, :), POINTER  :: newdat
    2318              :       TYPE(cp_blacs_type) :: ictxt_loc
    2319              :       INTEGER, EXTERNAL                        :: numroc
    2320              : #endif
    2321              : 
    2322         3278 :       CALL timeset(routineN, handle)
    2323              :       CALL cp_fm_get_info(fm, nrow_global=nrow_global, ncol_global=ncol_global, ncol_block=max_block, &
    2324         3278 :                           para_env=para_env)
    2325              : 
    2326         3278 :       IF (PRESENT(value_format)) THEN
    2327            0 :          CPASSERT(LEN_TRIM(ADJUSTL(value_format)) < 11)
    2328            0 :          my_value_format = "(I10, I10, "//TRIM(ADJUSTL(value_format))//")"
    2329              :       ELSE
    2330         3278 :          my_value_format = "(I10, I10, ES24.12)"
    2331              :       END IF
    2332              : 
    2333         3278 :       IF (unit > 0) THEN
    2334           11 :          IF (PRESENT(header)) WRITE (unit, *) header
    2335           11 :          WRITE (unit, "(A2, A8, A10, A24)") "#", "Row", "Column", ADJUSTL("Value")
    2336              :       END IF
    2337              : 
    2338              : #if defined(__parallel)
    2339         3278 :       num_pe = para_env%num_pe
    2340         3278 :       mepos = para_env%mepos
    2341         3278 :       rb = nrow_global
    2342         3278 :       tag = 0
    2343              :       ! get a new context
    2344         3278 :       CALL ictxt_loc%gridinit(para_env, 'R', 1, num_pe)
    2345         3278 :       CALL descinit(desc, nrow_global, ncol_global, rb, max_block, 0, 0, ictxt_loc%get_handle(), nrow_global, info)
    2346         3278 :       CPASSERT(info == 0)
    2347              :       ASSOCIATE (nprow => ictxt_loc%num_pe(1), npcol => ictxt_loc%num_pe(2), &
    2348              :                  myprow => ictxt_loc%mepos(1), mypcol => ictxt_loc%mepos(2))
    2349         6556 :          in = numroc(ncol_global, max_block, mypcol, 0, npcol)
    2350              : 
    2351        13112 :          ALLOCATE (newdat(nrow_global, MAX(1, in)))
    2352              : 
    2353              :          ! do the actual scalapack to cols reordering
    2354              :          CALL pdgemr2d(nrow_global, ncol_global, fm%local_data, 1, 1, &
    2355              :                        fm%matrix_struct%descriptor, &
    2356         3278 :                        newdat, 1, 1, desc, ictxt_loc%get_handle())
    2357              : 
    2358         9834 :          ALLOCATE (vecbuf(nrow_global*max_block))
    2359       103876 :          vecbuf = HUGE(1.0_dp) ! init for valgrind
    2360         3278 :          irow = 1
    2361         3278 :          icol = 1
    2362              : 
    2363        15250 :          DO i = 1, ncol_global, MAX(max_block, 1)
    2364         8694 :             i_block = MIN(max_block, ncol_global - i + 1)
    2365              :             CALL infog2l(1, i, desc, nprow, npcol, myprow, mypcol, &
    2366         8694 :                          irow_local, icol_local, iprow, ipcol)
    2367         8694 :             IF (ipcol == mypcol) THEN
    2368        17052 :                DO j = 1, i_block
    2369       299892 :                   vecbuf((j - 1)*nrow_global + 1:nrow_global*j) = newdat(:, icol_local + j - 1)
    2370              :                END DO
    2371              :             END IF
    2372              : 
    2373         8694 :             IF (ipcol == 0) THEN
    2374              :                ! do nothing
    2375              :             ELSE
    2376         3430 :                IF (ipcol == mypcol) THEN
    2377        50449 :                   CALL para_env%send(vecbuf(:), 0, tag)
    2378              :                END IF
    2379         3430 :                IF (mypcol == 0) THEN
    2380        99183 :                   CALL para_env%recv(vecbuf(:), ipcol, tag)
    2381              :                END IF
    2382              :             END IF
    2383              : 
    2384        11972 :             IF (unit > 0) THEN
    2385          210 :                DO j = 1, i_block
    2386         6438 :                   DO k = (j - 1)*nrow_global + 1, nrow_global*j
    2387         6228 :                      WRITE (UNIT=unit, FMT=my_value_format) irow, icol, vecbuf(k)
    2388         6228 :                      irow = irow + 1
    2389         6417 :                      IF (irow > nrow_global) THEN
    2390          189 :                         irow = 1
    2391          189 :                         icol = icol + 1
    2392              :                      END IF
    2393              :                   END DO
    2394              :                END DO
    2395              :             END IF
    2396              : 
    2397              :          END DO
    2398              :       END ASSOCIATE
    2399         3278 :       DEALLOCATE (vecbuf)
    2400              : 
    2401         3278 :       CALL ictxt_loc%gridexit()
    2402              : 
    2403         3278 :       DEALLOCATE (newdat)
    2404              : 
    2405              : #else
    2406              : 
    2407              :       IF (unit > 0) THEN
    2408              :          DO j = 1, ncol_global
    2409              :             DO i = 1, nrow_global
    2410              :                WRITE (UNIT=unit, FMT=my_value_format) i, j, fm%local_data(i, j)
    2411              :             END DO
    2412              :          END DO
    2413              :       END IF
    2414              : 
    2415              : #endif
    2416         3278 :       CALL timestop(handle)
    2417              : 
    2418        16390 :    END SUBROUTINE cp_fm_write_formatted
    2419              : 
    2420              : ! **************************************************************************************************
    2421              : !> \brief ...
    2422              : !> \param fm ...
    2423              : !> \param unit ...
    2424              : ! **************************************************************************************************
    2425         5814 :    SUBROUTINE cp_fm_read_unformatted(fm, unit)
    2426              :       TYPE(cp_fm_type), INTENT(INOUT)          :: fm
    2427              :       INTEGER, INTENT(IN)                      :: unit
    2428              : 
    2429              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_read_unformatted'
    2430              : 
    2431              :       INTEGER                                  :: handle, j, max_block, &
    2432              :                                                   ncol_global, nrow_global
    2433              :       TYPE(mp_para_env_type), POINTER          :: para_env
    2434              : #if defined(__parallel)
    2435              :       INTEGER                                  :: k, n_cols
    2436              :       REAL(KIND=dp), DIMENSION(:, :), POINTER  :: vecbuf
    2437              : #endif
    2438              : 
    2439         5814 :       CALL timeset(routineN, handle)
    2440              : 
    2441              :       CALL cp_fm_get_info(fm, nrow_global=nrow_global, ncol_global=ncol_global, ncol_block=max_block, &
    2442         5814 :                           para_env=para_env)
    2443              : 
    2444              : #if defined(__parallel)
    2445              : 
    2446              :       ! the parallel case could be made more efficient (see cp_fm_write_unformatted)
    2447              : 
    2448        23256 :       ALLOCATE (vecbuf(nrow_global, max_block))
    2449              : 
    2450        23474 :       DO j = 1, ncol_global, max_block
    2451              : 
    2452        17660 :          n_cols = MIN(max_block, ncol_global - j + 1)
    2453        17660 :          IF (para_env%mepos == 0) THEN
    2454        40708 :             DO k = 1, n_cols
    2455       626116 :                READ (unit) vecbuf(:, k)
    2456              :             END DO
    2457              :          END IF
    2458      2775924 :          CALL para_env%bcast(vecbuf, 0)
    2459        23474 :          CALL cp_fm_set_submatrix(fm, vecbuf, start_row=1, start_col=j, n_cols=n_cols)
    2460              : 
    2461              :       END DO
    2462              : 
    2463         5814 :       DEALLOCATE (vecbuf)
    2464              : 
    2465              : #else
    2466              : 
    2467              :       DO j = 1, ncol_global
    2468              :          READ (unit) fm%local_data(:, j)
    2469              :       END DO
    2470              : 
    2471              : #endif
    2472              : 
    2473         5814 :       CALL timestop(handle)
    2474              : 
    2475         5814 :    END SUBROUTINE cp_fm_read_unformatted
    2476              : 
    2477              : ! **************************************************************************************************
    2478              : !> \brief ...
    2479              : !> \param mm_type ...
    2480              : ! **************************************************************************************************
    2481        11793 :    SUBROUTINE cp_fm_setup(mm_type)
    2482              :       INTEGER, INTENT(IN)                                :: mm_type
    2483              : 
    2484        11793 :       cp_fm_mm_type = mm_type
    2485        11793 :    END SUBROUTINE cp_fm_setup
    2486              : 
    2487              : ! **************************************************************************************************
    2488              : !> \brief ...
    2489              : !> \return ...
    2490              : ! **************************************************************************************************
    2491      2016402 :    FUNCTION cp_fm_get_mm_type() RESULT(res)
    2492              :       INTEGER                                            :: res
    2493              : 
    2494      2016402 :       res = cp_fm_mm_type
    2495      2016402 :    END FUNCTION cp_fm_get_mm_type
    2496              : 
    2497              : ! **************************************************************************************************
    2498              : !> \brief ...
    2499              : !> \param ictxt ...
    2500              : !> \param prec ...
    2501              : !> \return ...
    2502              : ! **************************************************************************************************
    2503           10 :    FUNCTION cp_fm_pilaenv(ictxt, prec) RESULT(res)
    2504              :       INTEGER                                            :: ictxt
    2505              :       CHARACTER(LEN=1)                                   :: prec
    2506              :       INTEGER                                            :: res
    2507              : #if defined(__parallel)
    2508              :       INTEGER                                            :: pilaenv
    2509           10 :       res = pilaenv(ictxt, prec)
    2510              : #else
    2511              :       MARK_USED(ictxt)
    2512              :       MARK_USED(prec)
    2513              :       res = -1
    2514              : #endif
    2515              : 
    2516           10 :    END FUNCTION cp_fm_pilaenv
    2517              : 
    2518            0 : END MODULE cp_fm_types
        

Generated by: LCOV version 2.0-1