LCOV - code coverage report
Current view: top level - shared/common/src/28_numeric_noabirule - m_abi_linalg.F90 (source / functions) Coverage Total Hit
Test: coverage.info Lines: 71.5 % 158 113
Test Date: 2026-11-08 08:00:28 Functions: 60.0 % 5 3
Branches: 46.0 % 252 116

             Branch data     Line data    Source code
       1                 :             : !!****m* ABINIT/m_abi_linalg
       2                 :             : !! NAME
       3                 :             : !!  m_abi_linalg
       4                 :             : !!
       5                 :             : !! FUNCTION
       6                 :             : !!  management of Linear Algebra wrappers routines
       7                 :             : !!  with support of different external library (scalapack, elpa, plasma, magma, ... )
       8                 :             : !!
       9                 :             : !! COPYRIGHT
      10                 :             : !!  Copyright (C) 2012-2026 ABINIT group (LNguyen,FDahm,MT)
      11                 :             : !!  This file is distributed under the terms of the
      12                 :             : !!  GNU General Public License, see ~abinit/COPYING
      13                 :             : !!  or http://www.gnu.org/copyleft/gpl.txt .
      14                 :             : !!
      15                 :             : !! SOURCE
      16                 :             : 
      17                 :             : #if defined HAVE_CONFIG_H
      18                 :             : #include "config.h"
      19                 :             : #endif
      20                 :             : 
      21                 :             : #include "abi_common.h"
      22                 :             : 
      23                 :             : module m_abi_linalg
      24                 :             : 
      25                 :             :  use, intrinsic :: iso_c_binding
      26                 :             :  USE_MPI
      27                 :             :  use defs_basis
      28                 :             :  use m_errors
      29                 :             :  use m_abicore
      30                 :             :  use m_xmpi
      31                 :             :  use m_xomp
      32                 :             :  use m_gputk
      33                 :             :  use m_slk
      34                 :             : 
      35                 :             : !#ifdef HAVE_LINALG_ELPA
      36                 :             : ! use m_elpa
      37                 :             : !#endif
      38                 :             : #ifdef HAVE_LINALG_PLASMA
      39                 :             :  use plasma, except_dp => dp, except_sp => sp
      40                 :             : #endif
      41                 :             : 
      42                 :             : #if defined HAVE_GPU
      43                 :             :  use m_gpu_toolbox
      44                 :             : #endif
      45                 :             : 
      46                 :             : #if defined HAVE_YAKL
      47                 :             :  use gator_mod, only: gator_allocate, gator_deallocate
      48                 :             : #endif
      49                 :             : 
      50                 :             : #if defined HAVE_MPI1
      51                 :             :  include 'mpif.h'
      52                 :             : #endif
      53                 :             : 
      54                 :             :  use m_time,  only : timab
      55                 :             :  use m_fstrings, only : sjoin, itoa
      56                 :             : 
      57                 :             :  implicit none
      58                 :             : 
      59                 :             :  private
      60                 :             : !!***
      61                 :             : 
      62                 :             : !This flag is ON if abi_linalg functions are in use (if this module is in use)
      63                 :             :  logical :: abi_linalg_in_use=.true.
      64                 :             : 
      65                 :             : !These flags enable the different versions in the BLAS/LAPACK wrappers
      66                 :             :  logical,private,save :: ABI_LINALG_SCALAPACK_ISON=.False.
      67                 :             :  logical,private,save :: ABI_LINALG_MAGMA_ISON=.False.
      68                 :             :  logical,private,save :: ABI_LINALG_PLASMA_ISON=.False.
      69                 :             : 
      70                 :             : !Working arrays for eigen problem
      71                 :             : 
      72                 :             :  logical,save :: lapack_single_precision=.false.
      73                 :             :  logical,save :: lapack_double_precision=.false.
      74                 :             :  logical,save :: lapack_full_storage    =.false.
      75                 :             :  logical,save :: lapack_packed_storage  =.false.
      76                 :             :  logical,save :: lapack_divide_conquer  =.false.
      77                 :             : 
      78                 :             :  integer,save :: eigen_s_maxsize=0
      79                 :             :  integer,save :: eigen_d_maxsize=0
      80                 :             :  integer,save :: eigen_c_maxsize=0
      81                 :             :  integer,save :: eigen_z_maxsize=0
      82                 :             :  integer,save :: eigen_s_lwork=0
      83                 :             :  integer,save :: eigen_d_lwork=0
      84                 :             :  integer,save :: eigen_c_lwork=0
      85                 :             :  integer,save :: eigen_z_lwork=0
      86                 :             :  integer,save :: eigen_c_lrwork=0
      87                 :             :  integer,save :: eigen_z_lrwork=0
      88                 :             :  integer,save :: eigen_liwork=0
      89                 :             : 
      90                 :             :  ! MG FIXME: I really do not understand why we should use global variables to wrap Scalapack routines.
      91                 :             :  !
      92                 :             :  ! 1) The procedures are not thread-safe.
      93                 :             :  ! 2) The procedures cannot be reused with different dimensions without dellocating these globals
      94                 :             : 
      95                 :             :  integer,save,target,allocatable :: eigen_iwork(:)
      96                 :             :  real(sp),save,target,allocatable :: eigen_c_rwork(:)
      97                 :             :  real(dp),save,target,allocatable :: eigen_z_rwork(:)
      98                 :             :  real(sp),save,target,allocatable :: eigen_s_work(:)
      99                 :             :  real(dp),save,target,allocatable :: eigen_d_work(:)
     100                 :             :  complex(sp),save,target,allocatable :: eigen_c_work(:)
     101                 :             :  complex(dp),save,target,allocatable :: eigen_z_work(:)
     102                 :             : 
     103                 :             :  integer,save,private :: slk_minsize=1
     104                 :             :  integer,save,public :: slk_communicator=xmpi_comm_null
     105                 :             :  integer,save,public :: slk_complement_communicator=xmpi_comm_null
     106                 :             : #ifdef HAVE_LINALG_SCALAPACK
     107                 :             :  type(slk_processor_t),save,public :: slk_processor
     108                 :             : #endif
     109                 :             : 
     110                 :             : !Plasma can be activated via command line
     111                 :             : !Use XPLASMA_ISON flag and modifiy it with linalg_allow_plasma
     112                 :             :  logical,private,save :: XPLASMA_ISON=.false.
     113                 :             :  public :: linalg_allow_plasma
     114                 :             : #ifdef HAVE_LINALG_PLASMA
     115                 :             :  type(c_ptr) :: plasma_work
     116                 :             : #endif
     117                 :             : 
     118                 :             :  integer, save, private     :: abi_linalg_gpu_mode = ABI_GPU_DISABLED
     119                 :             : 
     120                 :             : #ifdef HAVE_GPU
     121                 :             :  integer,                       allocatable,save,private,target :: i_work(:)
     122                 :             :  real(kind=c_double),           allocatable,save,private,target :: r_work(:)
     123                 :             :  complex(kind=c_double_complex),allocatable,save,private,target :: c_work(:)
     124                 :             :  type(c_ptr),save,private :: gpu_work
     125                 :             : 
     126                 :             :  !FIXME *_managed arrays are only used with YAKL, in place of previous ones
     127                 :             :  integer(kind=c_int32_t),        contiguous, pointer,save,private :: i_work_managed(:) => null()
     128                 :             :  real(kind=c_double),            contiguous, pointer,save,private :: r_work_managed(:) => null()
     129                 :             :  complex(kind=c_double_complex), contiguous, pointer,save,private :: c_work_managed(:) => null()
     130                 :             : 
     131                 :             :  integer, save, private :: i_work_len = 0
     132                 :             :  integer, save, private :: r_work_len = 0
     133                 :             :  integer, save, private :: c_work_len = 0
     134                 :             :  integer(c_size_t), save, private :: gpu_work_len = 0
     135                 :             : #endif
     136                 :             : 
     137                 :             : !----------------------------------------------------------------------
     138                 :             : !!***
     139                 :             : 
     140                 :             :  !Procedures ------------------------------------
     141                 :             :  public :: abi_linalg_init          ! Initialization routine
     142                 :             :  public :: abi_linalg_finalize      ! CleanuUp routine
     143                 :             :  public :: abi_linalg_work_allocate ! Allocate work arrays
     144                 :             :  !----------------------------------------------------------------------
     145                 :             : 
     146                 :             : !BLAS INTERFACE
     147                 :             :  !public :: abi_zgemm
     148                 :             :  public :: abi_xgemm
     149                 :             : 
     150                 :             :  interface abi_xgemm
     151                 :             :     module procedure abi_zgemm_2d
     152                 :             :     module procedure abi_zgemm_3d
     153                 :             :     module procedure abi_zgemm_2dd
     154                 :             :     module procedure abi_d2zgemm
     155                 :             :     module procedure abi_d2zgemm_2d
     156                 :             :     module procedure abi_d2zgemm_233
     157                 :             :     module procedure abi_d2zgemm_313
     158                 :             :     module procedure abi_d2zgemm_331
     159                 :             :     module procedure abi_d2zgemm_333
     160                 :             :     module procedure abi_d2zgemm_334
     161                 :             :  end interface abi_xgemm
     162                 :             : 
     163                 :             :  interface abi_gpu_xgemm
     164                 :             :     module procedure abi_gpu_xgemm_cptr
     165                 :             :     module procedure abi_gpu_xgemm_d
     166                 :             :     module procedure abi_gpu_xgemm_z
     167                 :             :     module procedure abi_gpu_xgemm_2d
     168                 :             :     module procedure abi_gpu_xgemm_2z
     169                 :             :  end interface abi_gpu_xgemm
     170                 :             : 
     171                 :             :  public :: abi_gpu_xgemm_d
     172                 :             : 
     173                 :             :  interface abi_gpu_xgemm_strided
     174                 :             :     module procedure abi_gpu_xgemm_strided_cptr
     175                 :             :     module procedure abi_gpu_xgemm_strided_d
     176                 :             :     module procedure abi_gpu_xgemm_strided_z
     177                 :             :     module procedure abi_gpu_xgemm_strided_2d
     178                 :             :     module procedure abi_gpu_xgemm_strided_2z
     179                 :             :  end interface abi_gpu_xgemm_strided
     180                 :             : 
     181                 :             :  interface abi_gpu_xsymm
     182                 :             :     module procedure abi_gpu_xsymm_cptr
     183                 :             :     module procedure abi_gpu_xsymm_d
     184                 :             :     module procedure abi_gpu_xsymm_z
     185                 :             :     module procedure abi_gpu_xsymm_2d
     186                 :             :     module procedure abi_gpu_xsymm_2z
     187                 :             :  end interface abi_gpu_xsymm
     188                 :             : 
     189                 :             :  interface abi_gpu_zhemm
     190                 :             :     module procedure abi_gpu_zhemm_cptr
     191                 :             :     module procedure abi_gpu_zhemm_d
     192                 :             :     module procedure abi_gpu_zhemm_z
     193                 :             :     module procedure abi_gpu_zhemm_2d
     194                 :             :     module procedure abi_gpu_zhemm_2z
     195                 :             :  end interface abi_gpu_zhemm
     196                 :             : 
     197                 :             :  interface abi_gpu_xscal
     198                 :             :     module procedure abi_gpu_xscal_cptr
     199                 :             :     module procedure abi_gpu_xscal_d
     200                 :             :     module procedure abi_gpu_xscal_z
     201                 :             :     module procedure abi_gpu_xscal_2d
     202                 :             :     module procedure abi_gpu_xscal_2z
     203                 :             :  end interface abi_gpu_xscal
     204                 :             : 
     205                 :             :  interface abi_gpu_xdot
     206                 :             :     module procedure abi_gpu_xdot_cptr
     207                 :             :     module procedure abi_gpu_xdot_d
     208                 :             :     module procedure abi_gpu_xdot_z
     209                 :             :  end interface abi_gpu_xdot
     210                 :             : 
     211                 :             :  interface abi_gpu_xaxpy
     212                 :             :     module procedure abi_gpu_xaxpy_cptr
     213                 :             :     module procedure abi_gpu_xaxpy_d
     214                 :             :     module procedure abi_gpu_xaxpy_z
     215                 :             :     module procedure abi_gpu_xaxpy_2d
     216                 :             :     module procedure abi_gpu_xaxpy_2z
     217                 :             :  end interface abi_gpu_xaxpy
     218                 :             : 
     219                 :             :  interface abi_gpu_xheevd
     220                 :             :     module procedure abi_gpu_xheevd_cptr
     221                 :             :     module procedure abi_gpu_xheevd_d
     222                 :             :     module procedure abi_gpu_xheevd_z
     223                 :             :     module procedure abi_gpu_xheevd_2d
     224                 :             :     module procedure abi_gpu_xheevd_2z
     225                 :             :  end interface abi_gpu_xheevd
     226                 :             : 
     227                 :             :  interface abi_gpu_xhegvd
     228                 :             :     module procedure abi_gpu_xhegvd_cptr
     229                 :             :     module procedure abi_gpu_xhegvd_d
     230                 :             :     module procedure abi_gpu_xhegvd_z
     231                 :             :     module procedure abi_gpu_xhegvd_2d
     232                 :             :     module procedure abi_gpu_xhegvd_2z
     233                 :             :  end interface abi_gpu_xhegvd
     234                 :             : 
     235                 :             :  interface abi_gpu_xtrsm
     236                 :             :     module procedure abi_gpu_xtrsm_cptr
     237                 :             :     module procedure abi_gpu_xtrsm_d
     238                 :             :     module procedure abi_gpu_xtrsm_z
     239                 :             :     module procedure abi_gpu_xtrsm_2d
     240                 :             :     module procedure abi_gpu_xtrsm_2z
     241                 :             :  end interface abi_gpu_xtrsm
     242                 :             : 
     243                 :             :  interface abi_gpu_xpotrf
     244                 :             :     module procedure abi_gpu_xpotrf_cptr
     245                 :             :     module procedure abi_gpu_xpotrf_d
     246                 :             :     module procedure abi_gpu_xpotrf_z
     247                 :             :     module procedure abi_gpu_xpotrf_2d
     248                 :             :     module procedure abi_gpu_xpotrf_2z
     249                 :             :  end interface abi_gpu_xpotrf
     250                 :             : 
     251                 :             :  interface abi_gpu_xcopy
     252                 :             :     module procedure abi_gpu_xcopy_cptr
     253                 :             :     module procedure abi_gpu_xcopy_d
     254                 :             :     module procedure abi_gpu_xcopy_z
     255                 :             :     module procedure abi_gpu_xcopy_2d
     256                 :             :     module procedure abi_gpu_xcopy_2z
     257                 :             :  end interface abi_gpu_xcopy
     258                 :             : 
     259                 :             :  interface abi_gpu_work_resize
     260                 :             :     module procedure abi_gpu_work_resizeI
     261                 :             :     module procedure abi_gpu_work_resizeR
     262                 :             :     module procedure abi_gpu_work_resizeC
     263                 :             :  end interface abi_gpu_work_resize
     264                 :             : 
     265                 :             :  public :: abi_zgemm
     266                 :             :  public :: abi_zgemm_2d
     267                 :             :  public :: abi_zgemm_2dd
     268                 :             :  public :: abi_zgemm_2r
     269                 :             :  interface abi_zgemm  ! No x_cplx stuff here!
     270                 :             :     module procedure abi_zgemm_2d
     271                 :             :     module procedure abi_zgemm_3d
     272                 :             :     module procedure abi_zgemm_2r
     273                 :             :  end interface abi_zgemm
     274                 :             : 
     275                 :             :  !----------------------------------------------------------------------
     276                 :             :  public :: abi_xscal
     277                 :             :  interface abi_xscal
     278                 :             :     module procedure abi_d2zscal
     279                 :             :     module procedure abi_d2zscal_3d
     280                 :             :     module procedure abi_d2zscal_4d
     281                 :             :     module procedure abi_d2zscal_5d
     282                 :             :     module procedure abi_d2zscal_7d
     283                 :             :     module procedure abi_dscal
     284                 :             :     module procedure abi_dscal_2d
     285                 :             :     module procedure abi_zscal
     286                 :             :     module procedure abi_zscal_2d
     287                 :             :     module procedure abi_zscal_3d
     288                 :             :  end interface abi_xscal
     289                 :             :  !----------------------------------------------------------------------
     290                 :             :  public :: abi_xaxpy
     291                 :             :  interface abi_xaxpy
     292                 :             :     module procedure abi_daxpy
     293                 :             :     module procedure abi_daxpy_2d
     294                 :             :     module procedure abi_d2zaxpy
     295                 :             :     module procedure abi_d2zaxpy_2d
     296                 :             :     module procedure abi_d2zaxpy_3d
     297                 :             :     module procedure abi_d2zaxpy_5d
     298                 :             :     module procedure abi_zaxpy
     299                 :             :     module procedure abi_zaxpy_2d
     300                 :             :     module procedure abi_zaxpy_3d
     301                 :             :  end interface abi_xaxpy
     302                 :             :  !----------------------------------------------------------------------
     303                 :             :  public :: abi_xcopy
     304                 :             :  interface abi_xcopy
     305                 :             :     module procedure abi_zcopy
     306                 :             :     module procedure abi_zcopy_1d
     307                 :             :     module procedure abi_zcopy_2d
     308                 :             :     module procedure abi_dcopy
     309                 :             :     module procedure abi_dcopy_1d
     310                 :             :     module procedure abi_dcopy_2d     ! FIXME To be removed. One can pass the base adress of the array!
     311                 :             :     module procedure abi_dcopy_0d_1d
     312                 :             :     module procedure abi_dcopy_1d_0d
     313                 :             :     module procedure abi_d2zcopy_2d  ! FIXME To be removed. One can pass the base adress of the array!
     314                 :             :     module procedure abi_z2dcopy_2d  ! FIXME To be removed. One can pass the base adress of the array!
     315                 :             :  end interface abi_xcopy
     316                 :             : 
     317                 :             :  !----------------------------------------------------------------------
     318                 :             :  public :: abi_xtrsm
     319                 :             :  interface abi_xtrsm
     320                 :             :     module procedure abi_ztrsm
     321                 :             :     module procedure abi_dtrsm
     322                 :             :     module procedure abi_d2ztrsm
     323                 :             :     !module procedure abi_d2ztrsm_3d
     324                 :             :  end interface abi_xtrsm
     325                 :             : 
     326                 :             :  public :: abi_d2ztrsm_3d ! Used in bestwfk TODO to be Removed
     327                 :             :  !----------------------------------------------------------------------
     328                 :             : 
     329                 :             : !LAPACK INTERFACE
     330                 :             :  public :: abi_xheev
     331                 :             :  interface abi_xheev
     332                 :             :     module procedure abi_dheev
     333                 :             :     module procedure abi_cheev
     334                 :             :     module procedure abi_zheev
     335                 :             :  end interface
     336                 :             :  !----------------------------------------------------------------------
     337                 :             :  public :: abi_xheevd
     338                 :             :  interface abi_xheevd
     339                 :             :     module procedure abi_d2zheevd
     340                 :             :     module procedure abi_zheevd_2d
     341                 :             :  end interface
     342                 :             :  !----------------------------------------------------------------------
     343                 :             :  public :: abi_xhegv
     344                 :             :  interface abi_xhegv
     345                 :             :     module procedure abi_dhegv
     346                 :             :     module procedure abi_chegv
     347                 :             :     module procedure abi_zhegv
     348                 :             :  end interface
     349                 :             :  !----------------------------------------------------------------------
     350                 :             :  public :: abi_xhegvd
     351                 :             :  interface abi_xhegvd
     352                 :             :     module procedure abi_d2zhegvd
     353                 :             :     module procedure abi_zhegvd_2d
     354                 :             :  end interface
     355                 :             :  !----------------------------------------------------------------------
     356                 :             :  public :: abi_xhpev
     357                 :             :  interface abi_xhpev
     358                 :             :     module procedure abi_dhpev
     359                 :             :     module procedure abi_chpev
     360                 :             :     module procedure abi_zhpev
     361                 :             :  end interface
     362                 :             :  !----------------------------------------------------------------------
     363                 :             :  public :: abi_xhpgv
     364                 :             :  interface abi_xhpgv
     365                 :             :     module procedure abi_dhpgv
     366                 :             :     module procedure abi_chpgv
     367                 :             :     module procedure abi_zhpgv
     368                 :             :  end interface
     369                 :             :  !----------------------------------------------------------------------
     370                 :             :  public :: abi_xpotrf
     371                 :             :  interface abi_xpotrf
     372                 :             :     module procedure abi_dpotrf
     373                 :             :     module procedure abi_d2zpotrf
     374                 :             :     module procedure abi_d2zpotrf_3d
     375                 :             :     module procedure abi_zpotrf_2d
     376                 :             :     module procedure abi_zpotrf
     377                 :             :  end interface
     378                 :             :  !----------------------------------------------------------------------
     379                 :             :  public :: abi_xorthonormalize
     380                 :             :  interface abi_xorthonormalize
     381                 :             :     module procedure xorthonormalize
     382                 :             :     module procedure zorthonormalize
     383                 :             :  end interface
     384                 :             : 
     385                 :             :  ! This version operates on arrays in which the real and the imaginary part
     386                 :             :  ! are packed together (real parts first, them imaginary parts), used when gamma-point and istwfk=2
     387                 :             :  public :: ortho_reim
     388                 :             :  !----------------------------------------------------------------------
     389                 :             : 
     390                 :             : #ifndef HAVE_GPU
     391                 :             :  !dummy routines replace gpu helper routines
     392                 :             :  public :: gpu_linalg_init
     393                 :             :  public :: gpu_linalg_shutdown
     394                 :             :  public :: gpu_xgemm
     395                 :             :  public :: gpu_xtrsm
     396                 :             :  public :: gpu_xaxpy
     397                 :             :  public :: gpu_xcopy
     398                 :             :  public :: gpu_xscal
     399                 :             :  public :: gpu_xdot
     400                 :             :  public :: gpu_xsygvd
     401                 :             :  public :: gpu_xsygvd_bufferSize
     402                 :             : #endif
     403                 :             : 
     404                 :             :  public :: gpu_xorthonormalize
     405                 :             : 
     406                 :             :  public :: abi_gpu_xgemm
     407                 :             :  public :: abi_gpu_xgemm_strided
     408                 :             :  public :: abi_gpu_xsymm
     409                 :             :  public :: abi_gpu_zhemm
     410                 :             :  public :: abi_gpu_xscal
     411                 :             :  public :: abi_gpu_xdot
     412                 :             :  public :: abi_gpu_xaxpy
     413                 :             :  public :: abi_gpu_xcopy
     414                 :             :  public :: abi_gpu_xtrsm
     415                 :             :  public :: abi_gpu_xhegvd
     416                 :             :  public :: abi_gpu_xheevd
     417                 :             :  public :: abi_gpu_xpotrf
     418                 :             : 
     419                 :             :  logical,external :: LSAME
     420                 :             : 
     421                 :             :  ! Timab slots, used if we want to profile BLAS calls
     422                 :             :  ! Fine-grained profiling, must be enabled with the CPP option DEV_LINALG_TIMING
     423                 :             :  ! For the time being, I use the same slots employed in lobpcgwf although
     424                 :             :  ! one should define specialized entries.If the index of the slot 0, no profiling is done.
     425                 :             : 
     426                 :             :  integer,parameter,private :: TIMAB_XCOPY=584
     427                 :             :  integer,parameter,private :: TIMAB_XGEMM=532
     428                 :             :  integer,parameter,private :: TIMAB_XORTHO=535
     429                 :             :  integer,parameter,private :: TIMAB_XEIGEN=587
     430                 :             :  integer,parameter,private :: TIMAB_XPRECO=536
     431                 :             :  integer,parameter,private :: TIMAB_WFCOPY=584
     432                 :             :  integer,parameter,private :: TIMAB_XTRSM=535
     433                 :             : 
     434                 :             : ! Define this variable to activate timing routines
     435                 :             : !#define DEV_LINALG_TIMING 1
     436                 :             : 
     437                 :             : ! Support for [Z,C]GEMM3M routines
     438                 :             :  logical,save,private :: XGEMM3M_ISON = .False.
     439                 :             :  !logical,save,private :: XGEMM3M_ISON = .True.
     440                 :             :  ! True if [Z,C]GEMM3M can be used (can be set with linalg_allow_gemm3m)
     441                 :             : 
     442                 :             :  public :: linalg_allow_gemm3m
     443                 :             : 
     444                 :             :  ! Thresholds for the activation of [Z,C]GEMM3M
     445                 :             :  integer,parameter,private :: ZGEMM3M_LIMIT = 325000
     446                 :             :  integer,parameter,private :: CGEMM3M_LIMIT = 200000
     447                 :             : 
     448                 :             : ! Handy macros
     449                 :             : #ifdef HAVE_LINALG_GEMM3M
     450                 :             : #define _ZGEMM3M ZGEMM3M
     451                 :             : #define _CGEMM3M CGEMM3M
     452                 :             : #else
     453                 :             : #define _ZGEMM3M ZGEMM
     454                 :             : #define _CGEMM3M CGEMM
     455                 :             : #endif
     456                 :             : 
     457                 :             : 
     458                 :             : CONTAINS  !===========================================================
     459                 :             : !!***
     460                 :             : 
     461                 :             : !!****f* m_abi_linalg/abi_linalg_init
     462                 :             : !! NAME
     463                 :             : !! abi_linalg_init
     464                 :             : !!
     465                 :             : !! FUNCTION
     466                 :             : !! Initalization of linear algebra environnement
     467                 :             : !!
     468                 :             : !! INPUTS
     469                 :             : !! max_eigen_pb_size= max. size of eigenproblem during calculation
     470                 :             : !! optdriver= type of calculation (ground-state, response function, GW, ...)
     471                 :             : !! wfoptalg= wave functions optimization algorithm (CG, LOBPCG, CHEBFI, ...)
     472                 :             : !! paral_kgb= 1 if (k,g,b) parallelism is on
     473                 :             : !! gpu_option = GPU implementation to use, i.e. cuda, openMP, ... (0=not using GPU)
     474                 :             : !! use_slk= 1 if use of Scalapack is on
     475                 :             : !! np_slk= max. number of processes to be used in Scalapack calls
     476                 :             : !! comm_scalapack= global communicator to be used in case of Scalapack
     477                 :             : !!
     478                 :             : !! SOURCE
     479                 :             : !!
     480                 :             : 
     481                 :        5461 :  subroutine abi_linalg_init(max_eigen_pb_size,optdriver,wfoptalg,paral_kgb,&
     482                 :             : &                           gpu_option,use_slk,np_slk,comm_scalapack)
     483                 :             : 
     484                 :             : !Arguments ------------------------------------
     485                 :             :  integer,intent(in) :: max_eigen_pb_size
     486                 :             :  integer,intent(in) :: optdriver,wfoptalg,paral_kgb
     487                 :             :  integer,intent(in) :: comm_scalapack,np_slk
     488                 :             :  integer,intent(in) :: gpu_option,use_slk
     489                 :             : 
     490                 :             : !Local variables ------------------------------
     491                 :             :  integer :: max_eigen_pb_size_eff=0
     492                 :             :  logical :: need_work_space=.true.
     493                 :             : #ifdef HAVE_LINALG_SCALAPACK
     494                 :             :  integer :: abi_info1,rank,commsize,commcart,sizecart(2)
     495                 :             :  logical :: reorder,periodic(2),keepdim(2)
     496                 :             : #endif
     497                 :             : #ifdef HAVE_LINALG_PLASMA
     498                 :             :  integer :: abi_info2,core_id,rank
     499                 :             :  integer :: num_cores=0,num_cores_node=0
     500                 :             :  integer,allocatable :: affinity(:)
     501                 :             : #endif
     502                 :             : !******************************************************************
     503                 :             : 
     504                 :             : !Use only abi_linalg in case of GS calculations
     505                 :             :  !abi_linalg_in_use=(optdriver==RUNL_GSTATE.or.optdriver==RUNL_GWLS.or.optdriver==RUNL_RESPFN)
     506                 :             :  !abi_linalg_in_use= any(optdriver == [RUNL_GSTATE, RUNL_GWLS, RUNL_RESPFN, RUNL_GWR, RUNL_EPH])
     507   [ +  +  +  + ]:        9222 :  abi_linalg_in_use= any(optdriver == [RUNL_GSTATE, RUNL_GWLS, RUNL_RESPFN, RUNL_EPH])
     508                 :             : 
     509                 :        5461 :  max_eigen_pb_size_eff=0
     510                 :        5461 :  lapack_single_precision=.false.
     511                 :        5461 :  lapack_double_precision=.false.
     512                 :        5461 :  lapack_full_storage    =.false.
     513                 :        5461 :  lapack_packed_storage  =.false.
     514                 :        5461 :  lapack_divide_conquer  =.false.
     515                 :        5461 :  eigen_s_maxsize=0 ; eigen_d_maxsize=0
     516                 :        5461 :  eigen_c_maxsize=0 ; eigen_z_maxsize=0
     517                 :        5461 :  eigen_s_lwork=0   ; eigen_d_lwork=0
     518                 :        5461 :  eigen_c_lwork=0   ; eigen_z_lwork=0
     519                 :        5461 :  eigen_c_lrwork=0  ; eigen_z_lrwork=0
     520                 :        5461 :  eigen_liwork=0
     521                 :        5461 :  ABI_LINALG_SCALAPACK_ISON=.False.
     522                 :        5461 :  ABI_LINALG_MAGMA_ISON=.False.
     523                 :        5461 :  ABI_LINALG_PLASMA_ISON=.False.
     524                 :        5461 :  slk_communicator=xmpi_comm_null
     525                 :        5461 :  slk_complement_communicator=xmpi_comm_null
     526                 :        5461 :  slk_minsize=1
     527                 :             : 
     528                 :             : !Exit here if we don't use this abi_linalg module
     529         [ +  + ]:        5461 :  if (.not.abi_linalg_in_use) return
     530                 :             : 
     531                 :             : !Set Lapack parameters
     532                 :        4965 :  max_eigen_pb_size_eff=max_eigen_pb_size
     533   [ +  +  -  + ]:        4965 :  if (wfoptalg==4.or.wfoptalg==14.or.gpu_option/=ABI_GPU_DISABLED) max_eigen_pb_size_eff=3*max_eigen_pb_size_eff
     534   [ +  +  -  + ]:        4965 :  lapack_full_storage=(wfoptalg==4.or.wfoptalg==14.or.gpu_option/=ABI_GPU_DISABLED)
     535                 :        4965 :  lapack_packed_storage=.true.
     536                 :        4965 :  lapack_single_precision=.false.
     537                 :        4965 :  lapack_double_precision=.true.
     538                 :        4965 :  lapack_divide_conquer=.false.
     539                 :             : 
     540                 :             : !Set maximum sizes
     541                 :        4965 :  eigen_s_maxsize = max_eigen_pb_size_eff
     542                 :        4965 :  eigen_d_maxsize = max_eigen_pb_size_eff
     543                 :        4965 :  eigen_c_maxsize = max_eigen_pb_size_eff
     544                 :        4965 :  eigen_z_maxsize = max_eigen_pb_size_eff
     545                 :        4965 :  need_work_space=.true.
     546                 :             : 
     547                 :             : #if defined HAVE_LINALG_SCALAPACK && defined HAVE_MPI
     548                 :             :  if ((paral_kgb==1.or.use_slk==1).and.np_slk>0) then
     549                 :             :    rank=xmpi_comm_rank(comm_scalapack)
     550                 :             :    ! We create slk_communicator using a cartesian grid, and store its complement
     551                 :             :    commsize = MIN(np_slk, xmpi_comm_size(comm_scalapack))
     552                 :             :    sizecart = (/commsize, xmpi_comm_size(comm_scalapack)/commsize/)
     553                 :             :    periodic = (/.true.,.true./) ; reorder = .false.
     554                 :             :    call MPI_CART_CREATE(comm_scalapack,2,sizecart,periodic,reorder,commcart,abi_info1)
     555                 :             :    keepdim = (/.true., .false./)
     556                 :             :    call MPI_CART_SUB(commcart, keepdim, slk_communicator,abi_info1)
     557                 :             :    keepdim = (/.false., .true./)
     558                 :             :    call MPI_CART_SUB(commcart, keepdim, slk_complement_communicator,abi_info1)
     559                 :             :    call slk_processor%init(slk_communicator)
     560                 :             :    slk_minsize=maxval(slk_processor%grid%dims(1:2))
     561                 :             :    need_work_space=(use_slk/=1) ! In this case we never use the work arrays
     562                 :             :    ABI_LINALG_SCALAPACK_ISON = .true.
     563                 :             :  end if
     564                 :             : #else
     565                 :             :  ABI_UNUSED(comm_scalapack)
     566                 :             :  ABI_UNUSED(paral_kgb)
     567                 :             :  ABI_UNUSED(use_slk)
     568                 :             :  ABI_UNUSED(np_slk)
     569                 :             : #endif
     570                 :             : 
     571                 :             : !#ifdef HAVE_LINALG_ELPA
     572                 :             : ! call elpa_func_init()
     573                 :             : !#endif
     574                 :             : 
     575                 :             : #ifdef HAVE_LINALG_PLASMA
     576                 :             : !Plasma Initialization
     577                 :             : !Because use of hybrid use of mpi+openmp+plasma,
     578                 :             : !we need to set manually the thread bindings policy
     579                 :             : !to avoid conflicts between mpi process due to plasma
     580                 :             :  if (XPLASMA_ISON) then
     581                 :             :    num_cores=xomp_get_max_threads()
     582                 :             :    num_cores_node=xomp_get_num_cores_node()
     583                 :             :    rank=xmpi_comm_rank(xmpi_world)
     584                 :             :    if (num_cores_node == 0) then ! This means that OMP is not enabled.
     585                 :             :      num_cores_node = 1
     586                 :             :      ABI_WARNING("You are using PLASMA but OpenMP is not enabled in Abinit!")
     587                 :             :    end if
     588                 :             :    ABI_MALLOC(affinity,(num_cores))
     589                 :             :    do core_id =1,num_cores
     590                 :             :      affinity(core_id) = MOD(rank*num_cores + (core_id-1), num_cores_node)
     591                 :             :    end do
     592                 :             :    call PLASMA_Init_Affinity(num_cores,affinity(1),abi_info2)
     593                 :             :    ABI_FREE(affinity)
     594                 :             :    ABI_LINALG_PLASMA_ISON = .True.
     595                 :             :    lapack_divide_conquer=.true.
     596                 :             :  end if
     597                 :             : #endif
     598                 :             : 
     599                 :             : #ifdef HAVE_LINALG_MAGMA
     600                 :             : #ifdef HAVE_LINALG_MAGMA_15
     601                 :             :  call magmaf_init()
     602                 :             :  ABI_LINALG_MAGMA_ISON = .true.
     603                 :             :  lapack_divide_conquer=.true.
     604                 :             : #endif
     605                 :             : #endif
     606                 :             : 
     607                 :             : #ifdef HAVE_GPU
     608                 :             : !Cublas initialization
     609                 :             :  if (gpu_option/=ABI_GPU_DISABLED) call gpu_linalg_init()
     610                 :             :  abi_linalg_gpu_mode = gpu_option !FIXME Add a check for this
     611                 :             : #endif
     612                 :             : 
     613                 :        4965 :  if (need_work_space) call abi_linalg_work_allocate()
     614                 :             : 
     615                 :             :  end subroutine abi_linalg_init
     616                 :             : !!***
     617                 :             : 
     618                 :             : !!****f* m_abi_linalg/abi_linalg_work_allocate
     619                 :             : !! NAME
     620                 :             : !! abi_linalg_work_allocate
     621                 :             : !!
     622                 :             : !! FUNCTION
     623                 :             : !!
     624                 :             : !! INPUTS
     625                 :             : !!
     626                 :             : !! SOURCE
     627                 :             : !!
     628                 :        4965 :  subroutine abi_linalg_work_allocate()
     629                 :             : 
     630                 :             : !Arguments ------------------------------------
     631                 :             : 
     632                 :             : !Local variables ------------------------------
     633                 :             : #ifdef HAVE_LINALG_MAGMA
     634                 :             :  integer :: nb, magmaf_get_ssytrd_nb, magmaf_get_dsytrd_nb, magmaf_get_chetrd_nb, magmaf_get_zhetrd_nb
     635                 :             : #endif
     636                 :             : !******************************************************************
     637                 :             : 
     638                 :             : !Single precision WORK
     639                 :        4965 :  eigen_s_lwork = 0
     640         [ +  - ]:        4965 :  if (eigen_s_maxsize>0) then
     641         [ -  + ]:        4965 :    if (lapack_single_precision) then
     642         [ #  # ]:           0 :      if (lapack_full_storage) then
     643                 :           0 :        eigen_s_lwork = max(eigen_s_lwork,3*eigen_s_maxsize-1) ! SSYEV, SSYGV
     644                 :             :      end if
     645         [ #  # ]:           0 :      if (lapack_packed_storage) then
     646                 :           0 :        eigen_s_lwork = max(eigen_s_lwork,3*eigen_s_maxsize) ! SSPEV[D], SSPGV[D]
     647                 :             :      end if
     648         [ #  # ]:           0 :      if (lapack_divide_conquer) then
     649                 :           0 :        eigen_s_lwork = max(eigen_s_lwork,1+6*eigen_s_maxsize+2*eigen_s_maxsize**2) ! SSYEVD, SSYGVD
     650                 :             :      end if
     651                 :             :      if (ABI_LINALG_MAGMA_ISON) then
     652                 :             :        if (lapack_full_storage.and.lapack_divide_conquer) then
     653                 :             : #if defined HAVE_LINALG_MAGMA
     654                 :             :          nb=magmaf_get_ssytrd_nb(eigen_s_maxsize)
     655                 :             :          eigen_s_lwork = max(eigen_s_lwork,eigen_s_maxsize*(nb+2)) ! MAGMAF_SSYEVD, MAGMAF_SSYGVD
     656                 :             : #endif
     657                 :             :        end if
     658                 :             :      end if
     659         [ #  # ]:           0 :      if (ABI_LINALG_PLASMA_ISON) then
     660   [ #  #  #  # ]:           0 :        if (lapack_full_storage.and.lapack_divide_conquer) then
     661                 :           0 :          eigen_s_lwork = max(eigen_s_lwork,eigen_s_maxsize**2) ! PLASMA_SSYEV
     662                 :             :        end if
     663                 :             :      end if
     664                 :             :    end if
     665                 :             :  end if
     666         [ +  + ]:        4965 :  ABI_SFREE(eigen_s_work)
     667   [ +  -  -  -  :       14895 :  ABI_MALLOC(eigen_s_work,(eigen_s_lwork))
          +  -  -  +  -  
             +  -  +  -  
                      + ]
     668                 :             : 
     669                 :             : !Double precision WORK
     670                 :        4965 :  eigen_d_lwork = 0
     671         [ +  - ]:        4965 :  if (eigen_d_maxsize>0) then
     672         [ +  - ]:        4965 :    if (lapack_double_precision) then
     673         [ +  + ]:        4965 :      if (lapack_full_storage) then
     674                 :          25 :        eigen_d_lwork = max(eigen_d_lwork,3*eigen_d_maxsize-1) ! DSYEV, DSYGV
     675                 :             :      end if
     676         [ +  - ]:        4965 :      if (lapack_packed_storage) then
     677                 :        4965 :        eigen_d_lwork = max(eigen_d_lwork,3*eigen_d_maxsize) ! DSPEV[D], DSPGV[D]
     678                 :             :      end if
     679         [ -  + ]:        4965 :      if (lapack_divide_conquer) then
     680                 :           0 :        eigen_d_lwork = max(eigen_d_lwork,1+6*eigen_d_maxsize+2*eigen_d_maxsize**2) ! DSYEVD, DSYGVD
     681                 :             :      end if
     682                 :             :      if (ABI_LINALG_MAGMA_ISON) then
     683                 :             :        if (lapack_full_storage.and.lapack_divide_conquer) then
     684                 :             : #if defined HAVE_LINALG_MAGMA
     685                 :             :          nb=magmaf_get_dsytrd_nb(eigen_d_maxsize)
     686                 :             :          eigen_d_lwork = max(eigen_d_lwork,eigen_d_maxsize*(nb+2)) ! MAGMAF_DSYEVD, MAGMAF_DSYGVD
     687                 :             : #endif
     688                 :             :        end if
     689                 :             :      end if
     690         [ -  + ]:        4965 :      if (ABI_LINALG_PLASMA_ISON) then
     691   [ #  #  #  # ]:           0 :        if (lapack_full_storage.and.lapack_divide_conquer) then
     692                 :           0 :          eigen_d_lwork = max(eigen_d_lwork,eigen_d_maxsize**2) ! PLASMA_DSYEV
     693                 :             :        end if
     694                 :             :      end if
     695                 :             :    end if
     696                 :             :  end if
     697         [ +  + ]:        4965 :  ABI_SFREE(eigen_d_work)
     698   [ -  +  +  -  :       14895 :  ABI_MALLOC(eigen_d_work,(eigen_d_lwork))
          +  -  +  -  -  
             +  -  +  -  
                      + ]
     699                 :             : 
     700                 :             : !Single complex WORK
     701                 :        4965 :  eigen_c_lwork = 0
     702         [ +  - ]:        4965 :  if (eigen_c_maxsize>0) then
     703         [ -  + ]:        4965 :    if (lapack_single_precision) then
     704         [ #  # ]:           0 :      if (lapack_full_storage) then
     705                 :           0 :        eigen_c_lwork = max(eigen_c_lwork,2*eigen_c_maxsize-1) ! CHEEV, CHEGV
     706                 :             :      end if
     707         [ #  # ]:           0 :      if (lapack_packed_storage) then
     708                 :           0 :        eigen_c_lwork = max(eigen_c_lwork,2*eigen_c_maxsize) ! CHPEV[D], CHPGV[D]
     709                 :             :      end if
     710         [ #  # ]:           0 :      if (lapack_divide_conquer) then
     711                 :           0 :        eigen_c_lwork = max(eigen_c_lwork,2*eigen_c_maxsize+eigen_c_maxsize**2) ! CHEEVD, CHEGVD
     712                 :             :      end if
     713                 :             :      if (ABI_LINALG_MAGMA_ISON) then
     714                 :             :        if (lapack_full_storage.and.lapack_divide_conquer) then
     715                 :             : #if defined HAVE_LINALG_MAGMA
     716                 :             :          nb=magmaf_get_chetrd_nb(eigen_c_maxsize)
     717                 :             :          eigen_c_lwork = max(eigen_c_lwork,eigen_c_maxsize*(nb+1)) ! MAGMAF_CHEEVD, MAGMAF_CHEGVD
     718                 :             : #endif
     719                 :             :        end if
     720                 :             :      end if
     721         [ #  # ]:           0 :      if (ABI_LINALG_PLASMA_ISON) then
     722   [ #  #  #  # ]:           0 :        if (lapack_full_storage.and.lapack_divide_conquer) then
     723                 :           0 :          eigen_c_lwork = max(eigen_c_lwork,eigen_c_maxsize**2) ! PLASMA_CHEEV
     724                 :             :        end if
     725                 :             :      end if
     726                 :             :    end if
     727                 :             :  end if
     728         [ +  + ]:        4965 :  ABI_SFREE(eigen_c_work)
     729   [ +  -  -  -  :       14895 :  ABI_MALLOC(eigen_c_work,(eigen_c_lwork))
          +  -  -  +  -  
             +  -  +  -  
                      + ]
     730                 :             : 
     731                 :             : !Double complex WORK
     732                 :        4965 :  eigen_z_lwork = 0
     733         [ +  - ]:        4965 :  if (eigen_z_maxsize>0) then
     734         [ +  - ]:        4965 :    if (lapack_double_precision) then
     735         [ +  + ]:        4965 :      if (lapack_full_storage) then
     736                 :          25 :        eigen_z_lwork = max(eigen_z_lwork,2*eigen_z_maxsize-1) ! ZHEEV, ZHEGV
     737                 :             :      end if
     738         [ +  - ]:        4965 :      if (lapack_packed_storage) then
     739                 :        4965 :        eigen_z_lwork = max(eigen_z_lwork,2*eigen_z_maxsize) ! ZHPEV[D], ZHPGV[D]
     740                 :             :      end if
     741         [ -  + ]:        4965 :      if (lapack_divide_conquer) then
     742                 :           0 :        eigen_z_lwork = max(eigen_z_lwork,2*eigen_z_maxsize+eigen_z_maxsize**2) ! ZHEEVD, ZHEGVD
     743                 :             :      end if
     744                 :             :      if (ABI_LINALG_MAGMA_ISON) then
     745                 :             :        if (lapack_full_storage.and.lapack_divide_conquer) then
     746                 :             : #if defined HAVE_LINALG_MAGMA
     747                 :             :          nb=magmaf_get_zhetrd_nb(eigen_z_maxsize)
     748                 :             :          eigen_z_lwork = max(eigen_z_lwork,eigen_z_maxsize*(nb+1)) ! MAGMAF_ZHEEVD, MAGMAF_ZHEGVD
     749                 :             : #endif
     750                 :             :        end if
     751                 :             :      end if
     752         [ -  + ]:        4965 :      if (ABI_LINALG_PLASMA_ISON) then
     753   [ #  #  #  # ]:           0 :        if (lapack_full_storage.and.lapack_divide_conquer) then
     754                 :           0 :          eigen_z_lwork = max(eigen_z_lwork,eigen_z_maxsize**2) ! PLASMA_ZHEEV
     755                 :             :        end if
     756                 :             :      end if
     757                 :             :    end if
     758                 :             :  end if
     759         [ +  + ]:        4965 :  ABI_SFREE(eigen_z_work)
     760   [ -  +  +  -  :       14895 :  ABI_MALLOC(eigen_z_work,(eigen_z_lwork))
          +  -  +  -  -  
             +  -  +  -  
                      + ]
     761                 :             : 
     762                 :             : !Single precision RWORK
     763                 :        4965 :  eigen_c_lrwork = 0
     764         [ +  - ]:        4965 :  if (eigen_c_maxsize>0) then
     765         [ -  + ]:        4965 :    if (lapack_single_precision) then
     766   [ #  #  #  # ]:           0 :      if (lapack_full_storage.or.lapack_packed_storage) then
     767                 :           0 :        eigen_c_lrwork = max(eigen_c_lrwork,3*eigen_c_maxsize-2) ! CHEEV, CHEGV, CHPEV, CHPGV
     768                 :             :      end if
     769         [ #  # ]:           0 :      if (lapack_divide_conquer) then
     770                 :           0 :        eigen_c_lrwork = max(eigen_c_lrwork,1+5*eigen_c_maxsize+2*eigen_c_maxsize**2) ! CHEEVD, CHEGVD, CHPEVD, CHPGVD
     771                 :             :      end if
     772                 :             :    end if
     773                 :             :  end if
     774         [ +  + ]:        4965 :  ABI_SFREE(eigen_c_rwork)
     775   [ +  -  -  -  :       14895 :  ABI_MALLOC(eigen_c_rwork,(eigen_c_lrwork))
          +  -  -  +  -  
             +  -  +  -  
                      + ]
     776                 :             : 
     777                 :             : !Double precision RWORK
     778                 :        4965 :  eigen_z_lrwork = 0
     779         [ +  - ]:        4965 :  if (eigen_z_maxsize>0) then
     780         [ +  - ]:        4965 :    if (lapack_double_precision) then
     781   [ +  +  +  - ]:        4965 :      if (lapack_full_storage.or.lapack_packed_storage) then
     782                 :        4965 :        eigen_z_lrwork = max(eigen_z_lrwork,3*eigen_z_maxsize-2) ! ZHEEV, ZHEGV, ZHPEV, ZHPGV
     783                 :             :      end if
     784         [ -  + ]:        4965 :      if (lapack_divide_conquer) then
     785                 :           0 :        eigen_z_lrwork = max(eigen_z_lrwork,1+5*eigen_z_maxsize+2*eigen_z_maxsize**2) ! ZHEEVD, ZHEGVD, ZHPEVD, ZHPGVD
     786                 :             :      end if
     787                 :             :    end if
     788                 :             :  end if
     789         [ +  + ]:        4965 :  ABI_SFREE(eigen_z_rwork)
     790   [ -  +  +  -  :       14895 :  ABI_MALLOC(eigen_z_rwork,(eigen_z_lrwork))
          +  -  +  -  -  
             +  -  +  -  
                      + ]
     791                 :             : 
     792                 :             : !Integer IWORK
     793                 :        4965 :  eigen_liwork = 0
     794         [ -  + ]:        4965 :  if (lapack_divide_conquer) then
     795         [ #  # ]:           0 :    if (lapack_single_precision) then
     796         [ #  # ]:           0 :      if (eigen_s_maxsize>0) eigen_liwork = max(eigen_liwork,3+5*eigen_s_maxsize)
     797         [ #  # ]:           0 :      if (eigen_c_maxsize>0) eigen_liwork = max(eigen_liwork,3+5*eigen_c_maxsize)
     798                 :             :    end if
     799         [ #  # ]:           0 :    if (lapack_double_precision) then
     800         [ #  # ]:           0 :      if (eigen_d_maxsize>0) eigen_liwork = max(eigen_liwork,3+5*eigen_d_maxsize)
     801         [ #  # ]:           0 :      if (eigen_z_maxsize>0) eigen_liwork = max(eigen_liwork,3+5*eigen_z_maxsize)
     802                 :             :    end if
     803                 :             :  end if
     804         [ +  + ]:        4965 :  ABI_SFREE(eigen_iwork)
     805   [ +  -  -  -  :       14895 :  ABI_MALLOC(eigen_iwork,(eigen_liwork))
          +  -  -  +  -  
             +  -  +  -  
                      + ]
     806                 :             : 
     807                 :        4965 :  end subroutine abi_linalg_work_allocate
     808                 :             : !!***
     809                 :             : 
     810                 :             : !!****f* m_abi_linalg/abi_linalg_finalize
     811                 :             : !! NAME
     812                 :             : !! abi_linalg_finalize
     813                 :             : !!
     814                 :             : !! FUNCTION
     815                 :             : !!
     816                 :             : !! INPUTS
     817                 :             : !!
     818                 :             : !! SOURCE
     819                 :             : !!
     820                 :        5460 :  subroutine abi_linalg_finalize(gpu_option)
     821                 :             : 
     822                 :             : !Arguments ------------------------------------
     823                 :             :  integer, intent(in) :: gpu_option
     824                 :             : !Local variables ------------------------------
     825                 :             : #ifdef HAVE_LINALG_PLASMA
     826                 :             :  integer :: info
     827                 :             : #endif
     828                 :             : !******************************************************************
     829                 :             : 
     830         [ +  + ]:        5460 :  if (.not.abi_linalg_in_use) return
     831                 :             : 
     832                 :        4965 :  eigen_s_maxsize = 0
     833                 :        4965 :  eigen_d_maxsize = 0
     834                 :        4965 :  eigen_c_maxsize = 0
     835                 :        4965 :  eigen_z_maxsize = 0
     836                 :        4965 :  eigen_s_lwork   = 0
     837                 :        4965 :  eigen_d_lwork   = 0
     838                 :        4965 :  eigen_c_lwork   = 0
     839                 :        4965 :  eigen_z_lwork   = 0
     840                 :        4965 :  eigen_c_lrwork  = 0
     841                 :        4965 :  eigen_z_lrwork  = 0
     842                 :        4965 :  eigen_liwork    = 0
     843                 :             : 
     844                 :        4965 :  lapack_full_storage=.False.
     845                 :        4965 :  lapack_packed_storage=.False.
     846                 :        4965 :  lapack_single_precision=.False.
     847                 :        4965 :  lapack_double_precision=.False.
     848                 :        4965 :  lapack_divide_conquer=.false.
     849                 :             : 
     850                 :             : #ifdef HAVE_LINALG_SCALAPACK
     851                 :             :  if (ABI_LINALG_SCALAPACK_ISON) then
     852                 :             :    call slk_processor%free()
     853                 :             :    call xmpi_comm_free(slk_communicator)
     854                 :             :    call xmpi_comm_free(slk_complement_communicator)
     855                 :             :    slk_communicator=xmpi_comm_null
     856                 :             :    slk_complement_communicator=xmpi_comm_null
     857                 :             :    slk_minsize=1
     858                 :             :  end if
     859                 :             : #endif
     860                 :             : 
     861                 :             : !#ifdef HAVE_LINALG_ELPA
     862                 :             : ! call elpa_func_uninit()
     863                 :             : !#endif
     864                 :             : 
     865                 :             : #ifdef HAVE_LINALG_PLASMA
     866                 :             :    call PLASMA_Finalize(info)
     867                 :             :    ABI_LINALG_PLASMA_ISON=.False.
     868                 :             :  end if
     869                 :             : #endif
     870                 :             : 
     871                 :             : #ifdef HAVE_LINALG_MAGMA
     872                 :             : #ifdef HAVE_LINALG_MAGMA_15
     873                 :             :  if (ABI_LINALG_MAGMA_ISON) then
     874                 :             :    call magmaf_finalize()
     875                 :             :  end if
     876                 :             : #endif
     877                 :             : #endif
     878                 :             : 
     879                 :             : #ifdef HAVE_GPU
     880                 :             :  if (gpu_option/=ABI_GPU_DISABLED) then
     881                 :             :    call abi_gpu_work_finalize()
     882                 :             :    call gpu_linalg_shutdown()
     883                 :             :  end if
     884                 :             :  abi_linalg_gpu_mode = ABI_GPU_DISABLED
     885                 :             : #else
     886                 :             :  ABI_UNUSED(gpu_option)
     887                 :             : #endif
     888                 :             : 
     889                 :             : !Memory freeing
     890         [ +  + ]:        4965 :  ABI_SFREE(eigen_s_work)
     891         [ +  + ]:        4965 :  ABI_SFREE(eigen_d_work)
     892         [ +  + ]:        4965 :  ABI_SFREE(eigen_c_work)
     893         [ +  + ]:        4965 :  ABI_SFREE(eigen_z_work)
     894         [ +  + ]:        4965 :  ABI_SFREE(eigen_c_rwork)
     895         [ +  + ]:        4965 :  ABI_SFREE(eigen_z_rwork)
     896         [ +  + ]:        4965 :  ABI_SFREE(eigen_iwork)
     897                 :             : 
     898                 :             :  end subroutine abi_linalg_finalize
     899                 :             : !!***
     900                 :             : 
     901                 :             : !----------------------------------------------------------------------
     902                 :             : 
     903                 :             : !!****f* m_abi_linalg/linalg_allow_gemm3m
     904                 :             : !! NAME
     905                 :             : !!
     906                 :             : !! FUNCTION
     907                 :             : !!  Programmatic interface to enable the use of [Z,C]GEMM3M calls
     908                 :             : !!
     909                 :             : !! SOURCE
     910                 :             : 
     911                 :           0 : subroutine linalg_allow_gemm3m(bool, write_msg)
     912                 :             : 
     913                 :             : !Arguments ------------------------------------
     914                 :             :  logical,intent(in) :: bool, write_msg
     915                 :             : ! *************************************************************************
     916                 :             : 
     917                 :           0 :  XGEMM3M_ISON = bool
     918         [ #  # ]:           0 :  if (write_msg) then
     919                 :             : #ifdef HAVE_LINALG_GEMM3M
     920         [ #  # ]:           0 :    if (bool) then
     921                 :           0 :      ABI_COMMENT("Activating ZGEMM3M version instead of ZGEMM")
     922                 :             :    else
     923                 :           0 :      ABI_COMMENT("Using ZGEMM instead of ZGEMM3M")
     924                 :             :    end if
     925                 :             : #else
     926                 :             :    if (bool) then
     927                 :             :      ABI_WARNING("Cannot activate ZGEMM3M as HAVE_LINALG_GEMM3M is not defined!")
     928                 :             :    end if
     929                 :             : #endif
     930                 :             :  endif
     931                 :             : 
     932                 :           0 : end subroutine linalg_allow_gemm3m
     933                 :             : !!***
     934                 :             : 
     935                 :             : !----------------------------------------------------------------------
     936                 :             : 
     937                 :             : !!****f* m_abi_linalg/use_zgemm3m
     938                 :             : !! NAME
     939                 :             : !!  use_zgemm3m
     940                 :             : !!
     941                 :             : !! FUNCTION
     942                 :             : !!  Enable the use of ZGEMM3M
     943                 :             : !!
     944                 :             : !! NOTES
     945                 :             : !!  The CGEMM3M and ZGEMM3M routines use an algorithm requiring 3 real matrix
     946                 :             : !!  multiplications and 5 real matrix additions to compute the complex matrix
     947                 :             : !!  product; CGEMM(3S) and ZGEMM(3S) use 4 real matrix multiplications and 2
     948                 :             : !!  real matrix additions. Because the matrix multiplication time is usually
     949                 :             : !!  the limiting performance factor in these routines, CGEMM3M and ZGEMM3M
     950                 :             : !!  may run up to 33 percent faster than CGEMM and ZGEMM.  Because of other
     951                 :             : !!  overhead associated with the 3M routines, however, these performance
     952                 :             : !!  improvements may not always be realized.  For example, on one processor
     953                 :             : !!  the 3M routines will generally run more slowly than the standard complex
     954                 :             : !!  matrix multiplication routines when m * n * k < FACTOR, where m, n, and k
     955                 :             : !!  are the input matrix dimensions and FACTOR is approximately 200000 for
     956                 :             : !!  CGEMM3M and 325000 for ZGEMM3M.
     957                 :             : !!  from: http://techpubs.sgi.com/library/tpl/cgi-bin/getdoc.cgi?coll=0650&db=man&raw=1&fname=/usr/share/catman/p_man/cat3/SCSL/ZGEMM3M.z
     958                 :             : !!
     959                 :             : !! SOURCE
     960                 :             : 
     961                 :   178125528 : pure logical function use_zgemm3m(m, n, k)
     962                 :             : 
     963                 :             : !Arguments ------------------------------------
     964                 :             :  integer,intent(in) :: m,n,k
     965                 :             : ! *************************************************************************
     966                 :             : 
     967                 :   178125528 :  use_zgemm3m = .False.
     968   [ -  +  -  -  :   178125528 :  if (XGEMM3M_ISON) use_zgemm3m = ((m * n * k) > ZGEMM3M_LIMIT)
             -  +  -  - ]
     969                 :             :  !if (XGEMM3M_ISON) use_zgemm3m = .True.
     970                 :             : 
     971                 :             : #ifndef HAVE_LINALG_GEMM3M
     972                 :             :  use_zgemm3m = .False.
     973                 :             : #endif
     974                 :             : 
     975                 :             : end function use_zgemm3m
     976                 :             : !!***
     977                 :             : 
     978                 :             : !----------------------------------------------------------------------
     979                 :             : 
     980                 :             : !!****f* m_abi_linalg/use_cgemm3m
     981                 :             : !! NAME
     982                 :             : !!  use_cgemm3m
     983                 :             : !!
     984                 :             : !! FUNCTION
     985                 :             : !!  Enable the use of CGEMM3M
     986                 :             : !!
     987                 :             : !! NOTES
     988                 :             : !!  See use_zgemm3m
     989                 :             : !!
     990                 :             : !! SOURCE
     991                 :             : 
     992                 :             : pure logical function use_cgemm3m(m, n, k)
     993                 :             : 
     994                 :             : !Arguments ------------------------------------
     995                 :             :  integer,intent(in) :: m,n,k
     996                 :             : ! *************************************************************************
     997                 :             : 
     998                 :             :  use_cgemm3m = .False.
     999                 :             :  if (XGEMM3M_ISON) use_cgemm3m = ((m * n * k) > CGEMM3M_LIMIT)
    1000                 :             : #ifndef HAVE_LINALG_GEMM3M
    1001                 :             :  use_cgemm3m = .False.
    1002                 :             : #endif
    1003                 :             : 
    1004                 :             : end function use_cgemm3m
    1005                 :             : !!***
    1006                 :             : 
    1007                 :             : !----------------------------------------------------------------------
    1008                 :             : 
    1009                 :             : !!****f* m_abi_linalg/linalg_allow_plasma
    1010                 :             : !! NAME
    1011                 :             : !!
    1012                 :             : !! FUNCTION
    1013                 :             : !!  Programmatic interface to enable the use of PLASMA
    1014                 :             : !!  False to disable PLASMA version.
    1015                 :             : !!
    1016                 :             : !! SOURCE
    1017                 :             : 
    1018                 :           0 : subroutine linalg_allow_plasma(bool)
    1019                 :             : 
    1020                 :             : !Arguments ------------------------------------
    1021                 :             :  logical,intent(in) :: bool
    1022                 :             : ! *************************************************************************
    1023                 :             : 
    1024                 :             :  XPLASMA_ISON = bool
    1025                 :             : #ifndef HAVE_LINALG_PLASMA
    1026                 :             :  ! Just to be on the safe-side.
    1027                 :             :  ! I have to use a weird set of branches to make abirules happy in the BLAS/LAPACK
    1028                 :             :  ! wrappers, and one cannot set XPLASMA_MODE to .True. if PLASMA is not available.
    1029                 :           0 :  XPLASMA_ISON = .False.
    1030                 :             : #endif
    1031                 :             : 
    1032                 :           0 : end subroutine linalg_allow_plasma
    1033                 :             : !!***
    1034                 :             : 
    1035                 :             : #ifdef HAVE_LINALG_PLASMA
    1036                 :             : 
    1037                 :             : !!****f* m_abi_linalg/uplo_plasma
    1038                 :             : !! NAME
    1039                 :             : !!
    1040                 :             : !! FUNCTION
    1041                 :             : !!  Convert uplo character to PLASMA integer
    1042                 :             : !!
    1043                 :             : !! SOURCE
    1044                 :             : 
    1045                 :             : integer function uplo_plasma(uplo)
    1046                 :             : 
    1047                 :             : !Arguments ------------------------------------
    1048                 :             :  character(len=1),intent(in) :: uplo
    1049                 :             : ! *************************************************************************
    1050                 :             : 
    1051                 :             :  if (LSAME(uplo,'U')) then
    1052                 :             :     uplo_plasma = PlasmaUpper
    1053                 :             :  else
    1054                 :             :     uplo_plasma = PlasmaLower
    1055                 :             :  end if
    1056                 :             : 
    1057                 :             : end function uplo_plasma
    1058                 :             : !!***
    1059                 :             : 
    1060                 :             : !----------------------------------------------------------------------
    1061                 :             : 
    1062                 :             : !!****f* m_abi_linalg/trans_plasma
    1063                 :             : !! NAME
    1064                 :             : !!
    1065                 :             : !! FUNCTION
    1066                 :             : !!  Convert trans character to PLASMA integer
    1067                 :             : !!
    1068                 :             : !! SOURCE
    1069                 :             : 
    1070                 :             : integer function trans_plasma(trans)
    1071                 :             : 
    1072                 :             : !Arguments ------------------------------------
    1073                 :             :  character(len=1),intent(in) :: trans
    1074                 :             : ! *************************************************************************
    1075                 :             : 
    1076                 :             :  if (LSAME(trans,'C')) then
    1077                 :             :    trans_plasma = PlasmaConjTrans
    1078                 :             :  else if (LSAME(trans,'T')) then
    1079                 :             :    trans_plasma = PlasmaTrans
    1080                 :             :  else
    1081                 :             :    trans_plasma = PlasmaNoTrans
    1082                 :             :  end if
    1083                 :             : 
    1084                 :             : end function trans_plasma
    1085                 :             : !!***
    1086                 :             : 
    1087                 :             : !----------------------------------------------------------------------
    1088                 :             : 
    1089                 :             : !!****f* m_abi_linalg/side_plasma
    1090                 :             : !! NAME
    1091                 :             : !!
    1092                 :             : !! FUNCTION
    1093                 :             : !!  Convert side character to PLASMA integer
    1094                 :             : !!
    1095                 :             : !! SOURCE
    1096                 :             : 
    1097                 :             : integer function side_plasma(side)
    1098                 :             : 
    1099                 :             : !Arguments ------------------------------------
    1100                 :             :  character(len=1),intent(in) :: side
    1101                 :             : ! *************************************************************************
    1102                 :             : 
    1103                 :             :  if(LSAME(side,'L')) then
    1104                 :             :     side_plasma = PlasmaLeft
    1105                 :             :  else
    1106                 :             :     side_plasma = PlasmaRight
    1107                 :             :  end if
    1108                 :             : 
    1109                 :             : end function side_plasma
    1110                 :             : !!***
    1111                 :             : 
    1112                 :             : !----------------------------------------------------------------------
    1113                 :             : 
    1114                 :             : !!****f* m_abi_linalg/diag_plasma
    1115                 :             : !! NAME
    1116                 :             : !!
    1117                 :             : !! FUNCTION
    1118                 :             : !!  Convert diag character to PLASMA integer
    1119                 :             : !!
    1120                 :             : !! SOURCE
    1121                 :             : 
    1122                 :             : integer function diag_plasma(diag)
    1123                 :             : 
    1124                 :             : !Arguments ------------------------------------
    1125                 :             :  character(len=1),intent(in) :: diag
    1126                 :             : ! *************************************************************************
    1127                 :             : 
    1128                 :             :  if (LSAME(diag,'U')) then
    1129                 :             :    diag_plasma = PlasmaUnit
    1130                 :             :  else
    1131                 :             :    diag_plasma = PlasmaNonUnit
    1132                 :             :  end if
    1133                 :             : 
    1134                 :             : end function diag_plasma
    1135                 :             : !!***
    1136                 :             : 
    1137                 :             : !----------------------------------------------------------------------
    1138                 :             : 
    1139                 :             : !!****f* m_abi_linalg/jobz_plasma
    1140                 :             : !! NAME
    1141                 :             : !!
    1142                 :             : !! FUNCTION
    1143                 :             : !!  Convert jobz character to PLASMA integer
    1144                 :             : !!
    1145                 :             : !! SOURCE
    1146                 :             : 
    1147                 :             : integer function jobz_plasma(jobz)
    1148                 :             : 
    1149                 :             : !Arguments ------------------------------------
    1150                 :             :  character(len=1),intent(in) :: jobz
    1151                 :             : ! *************************************************************************
    1152                 :             : 
    1153                 :             :  if (LSAME(jobz,'N')) then
    1154                 :             :    jobz_plasma = PlasmaNoVec
    1155                 :             :  else
    1156                 :             :    jobz_plasma = PlasmaVec
    1157                 :             :  end if
    1158                 :             : 
    1159                 :             : end function jobz_plasma
    1160                 :             : !!***
    1161                 :             : 
    1162                 :             : #endif
    1163                 :             : 
    1164                 :             : ! Include files providing wrappers for some of the most commonly used BLAS & LAPACK routines
    1165                 :             : 
    1166                 :             : ! *********************************************************************!
    1167                 :             : ! ******************* BLAS_XGEMM interface ****************************!
    1168                 :             : ! *********************************************************************!
    1169                 :             : #include "abi_xgemm.f90"
    1170                 :             : 
    1171                 :             : ! *********************************************************************!
    1172                 :             : ! ******************** BLAS_XSCAL interface ***************************!
    1173                 :             : ! *********************************************************************!
    1174                 :             : #include "abi_xscal.f90"
    1175                 :             : 
    1176                 :             : ! *********************************************************************!
    1177                 :             : ! ******************** BLAS_XAXPY interface ***************************!
    1178                 :             : ! *********************************************************************!
    1179                 :             : #include "abi_xaxpy.f90"
    1180                 :             : 
    1181                 :             : ! *********************************************************************!
    1182                 :             : ! ******************** BLAS_XCOPY interface ***************************!
    1183                 :             : ! *********************************************************************!
    1184                 :             : #include "abi_xcopy.f90"
    1185                 :             : 
    1186                 :             : ! *********************************************************************!
    1187                 :             : ! ******************** BLAS_XTRSM interface ***************************!
    1188                 :             : ! *********************************************************************!
    1189                 :             : #include "abi_xtrsm.f90"
    1190                 :             : 
    1191                 :             : ! *********************************************************************!
    1192                 :             : ! ********************* LAPACK XHEEV  interface ***********************!
    1193                 :             : ! *********************************************************************!
    1194                 :             : #include "abi_xheev.f90"
    1195                 :             : 
    1196                 :             : ! *********************************************************************!
    1197                 :             : ! ********************* LAPACK XHEGV  interface ***********************!
    1198                 :             : ! *********************************************************************!
    1199                 :             : #include "abi_xhegv.f90"
    1200                 :             : 
    1201                 :             : ! *********************************************************************!
    1202                 :             : ! ********************* LAPACK XHPEV  interface ***********************!
    1203                 :             : ! *********************************************************************!
    1204                 :             : #include "abi_xhpev.f90"
    1205                 :             : 
    1206                 :             : ! *********************************************************************!
    1207                 :             : ! ********************* LAPACK XHPGV  interface ***********************!
    1208                 :             : ! *********************************************************************!
    1209                 :             : #include "abi_xhpgv.f90"
    1210                 :             : 
    1211                 :             : ! *********************************************************************!
    1212                 :             : ! ****************** LAPACK XPOTRF  interface *************************!
    1213                 :             : ! *********************************************************************!
    1214                 :             : #include "abi_xpotrf.f90"
    1215                 :             : 
    1216                 :             : ! *********************************************************************!
    1217                 :             : ! ************  ABINIT Orthonormalization interface *******************!
    1218                 :             : ! *********************************************************************!
    1219                 :             : #include "abi_xorthonormalize.f90"
    1220                 :             : 
    1221                 :             : ! *********************************************************************!
    1222                 :             : ! ******************    GPU LINALG interface  *************************!
    1223                 :             : ! *********************************************************************!
    1224                 :             : #include "abi_gpu_linalg.f90"
    1225                 :             : 
    1226                 :             : !----------------------------------------------------------------------
    1227                 :             : 
    1228                 :             : end module m_abi_linalg
    1229                 :             : !!***
        

Generated by: LCOV version 2.3-1