Line data    Source code

       1             : module m_apply_inverse_olap
       2             :    use m_glob_tofrom_loc
       3             :    USE m_types_mpimat
       4             : contains
       5          36 :    subroutine apply_inverse_olaps(mpdata, atoms, cell, hybdat, fmpi, sym, ikpt, coulomb)
       6             :       USE m_olap, ONLY: olap_pw
       7             :       USE m_types
       8             :       use m_judft
       9             :       implicit none
      10             :       type(t_mpdata), intent(in)  :: mpdata
      11             :       type(t_atoms), intent(in)   :: atoms
      12             :       type(t_cell), intent(in)    :: cell
      13             :       type(t_hybdat), intent(in)  :: hybdat
      14             :       type(t_mpi), intent(in)     :: fmpi
      15             :       type(t_sym), intent(in)     :: sym
      16             :       class(t_mat), intent(inout) :: coulomb
      17             :       integer, intent(in)         :: ikpt
      18             : 
      19          36 :       type(t_mat)               :: olap
      20          36 :       class(t_mat), allocatable :: coul_submtx
      21          36 :       type(t_mpimat)            :: olap_mpi
      22             : 
      23             :       integer         :: nbasm, loc_size, i, j, i_loc, ierr, pe_i, pe_j, pe_recv, pe_send, recv_loc, send_loc, j_loc
      24             :       complex         :: cdum
      25             : 
      26          36 :       call timestart("solve olap linear eq. sys")
      27          36 :       nbasm = hybdat%nbasp + mpdata%n_g(ikpt)
      28          36 :       CALL olap%alloc(.false., mpdata%n_g(ikpt), mpdata%n_g(ikpt), 0.0)
      29             :       !calculate IR overlap-matrix
      30       24652 :       CALL olap_pw(olap, mpdata%g(:, mpdata%gptm_ptr(:mpdata%n_g(ikpt), ikpt)), mpdata%n_g(ikpt), atoms, cell, fmpi)
      31             : 
      32             :       ! perform O^-1 * coulhlp%data_r(hybdat%nbasp + 1:, :) = x
      33             :       ! rewritten as O * x = C
      34             : 
      35          36 :       loc_size = 0
      36       15392 :       do i = 1, nbasm
      37       15356 :          call glob_to_loc(fmpi, i, pe_i, i_loc)
      38       15392 :          if (fmpi%n_rank == pe_i) loc_size = loc_size + 1
      39             :       end do
      40             : 
      41          36 :       call timestart("copy in 1")
      42          36 :       allocate(t_mat::coul_submtx)
      43          36 :       call coul_submtx%alloc(.false., mpdata%n_g(ikpt), loc_size)
      44     1217516 :       coul_submtx%data_c(:, :) = coulomb%data_c(hybdat%nbasp + 1:, :)
      45          36 :       call timestop("copy in 1")
      46             : 
      47             :       !$acc data copyin(olap, olap%data_r, olap%data_c, coul_submtx) copy(coul_submtx%data_r, coul_submtx%data_c)
      48          36 :          call olap%linear_problem(coul_submtx)
      49             :       !$acc end data
      50          36 :       call timestart("copy out 1")
      51     1217516 :       coulomb%data_c(hybdat%nbasp + 1:, :) = coul_submtx%data_c
      52          36 :       call coul_submtx%free()
      53          36 :       deallocate(coul_submtx)
      54          36 :       call timestop("copy out 1")
      55             : 
      56             : 
      57             :       ! perform  coulomb%data_r(hybdat%nbasp + 1:, :) * O^-1  = X
      58             :       ! rewritten as O^T * x^T = C^T
      59          36 :       call copy_in_2(fmpi, sym, mpdata, hybdat, coulomb, ikpt, coul_submtx)
      60             : 
      61             :       ! reload O, since the solver destroys it.
      62       24652 :       CALL olap_pw(olap, mpdata%g(:, mpdata%gptm_ptr(:mpdata%n_g(ikpt), ikpt)), mpdata%n_g(ikpt), atoms, cell, fmpi)
      63             :       ! Notice O = O^T since it's symmetric
      64             : 
      65             :       SELECT TYPE(coul_submtx)
      66             :       CLASS is (t_mat)
      67             :          !$acc data copyin(olap, olap%data_r, olap%data_c, coul_submtx) copy(coul_submtx%data_r, coul_submtx%data_c)
      68           0 :             call olap%linear_problem(coul_submtx)
      69             :          !$acc end data
      70           0 :          call olap%free()
      71             :       class is (t_mpimat)
      72          36 :          call olap_mpi%init(coul_submtx,  olap%matsize1, olap%matsize2)
      73          36 :          call olap_mpi%from_non_dist(olap)
      74          36 :          call olap_mpi%linear_problem(coul_submtx)
      75          36 :          call olap_mpi%free()
      76             :       end select
      77             : 
      78          36 :       call copy_out_2(fmpi, sym, mpdata, hybdat, ikpt, coul_submtx, coulomb)
      79          36 :       deallocate(coul_submtx)
      80          36 :       call timestop("solve olap linear eq. sys")
      81          36 :    end subroutine apply_inverse_olaps
      82             : 
      83          36 :    subroutine copy_in_2(fmpi, sym, mpdata, hybdat, coulomb, ikpt, coul_submtx)
      84             :       USE m_types
      85             :       implicit none 
      86             :       type(t_mpi), intent(in)      :: fmpi 
      87             :       integer, intent(in)          :: ikpt
      88             :       type(t_sym), intent(in)      :: sym
      89             :       type(t_mpdata), intent(in)   :: mpdata 
      90             :       type(t_hybdat), intent(in)   :: hybdat
      91             :       class(t_mat), intent(in)     :: coulomb 
      92             :       class(t_mat), intent(inout), allocatable  :: coul_submtx
      93             : 
      94             :       integer :: i, j, ierr, i_loc, j_loc, pe_i, pe_j
      95             :       complex :: cdum
      96             : 
      97          36 :       call timestart("copy in 2")
      98             : 
      99             :       SELECT TYPE(coulomb)
     100             :       CLASS is (t_mat)
     101           0 :          allocate(t_mat::coul_submtx)
     102           0 :          call coul_submtx%alloc(.false., mpdata%n_g(ikpt), mpdata%n_g(ikpt))
     103           0 :          do j = 1, mpdata%n_g(ikpt)
     104           0 :             do i = 1, mpdata%n_g(ikpt)
     105           0 :                coul_submtx%data_c(j, i) = conjg(coulomb%data_c(hybdat%nbasp+i, hybdat%nbasp + j))
     106             :             enddo 
     107             :          enddo
     108             :       class is (t_mpimat)
     109             : #ifdef CPP_SCALAPACK
     110          36 :          allocate(t_mpimat::coul_submtx)
     111          36 :          call coul_submtx%init(.False., mpdata%n_g(ikpt), mpdata%n_g(ikpt), fmpi%sub_comm, .True.)
     112             :          select type(coul_submtx)
     113             :          class is (t_mpimat)
     114             :             ! copy bottom right corner of coulomb to coul_submtx
     115             :             !call pzgemr2d(m,              n,               a,                ia,           ja,             desca, 
     116             :             call pzgemr2d(mpdata%n_g(ikpt),mpdata%n_g(ikpt),coulomb%data_c, hybdat%nbasp+1, hybdat%nbasp+1, coulomb%blacsdata%blacs_desc,&
     117             :             !             b, ib, jb,             descb, ictxt)
     118          36 :                         coul_submtx%data_c, 1, 1, coul_submtx%blacsdata%blacs_desc, coulomb%blacsdata%blacs_desc(2))
     119          36 :             call coul_submtx%transpose()
     120             :          class default
     121           0 :             call judft_error("coul_submtx should also be mpimat")
     122             :          end select
     123             : #endif
     124             :       END SELECT
     125          36 :       call timestop("copy in 2")
     126          36 :    end subroutine copy_in_2
     127             : 
     128          36 :    subroutine copy_out_2(fmpi, sym, mpdata, hybdat, ikpt, coul_submtx, coulomb)
     129             :       USE m_types
     130             :       implicit none 
     131             :       type(t_mpi), intent(in)      :: fmpi 
     132             :       integer, intent(in)          :: ikpt
     133             :       type(t_sym), intent(in)      :: sym
     134             :       type(t_mpdata), intent(in)   :: mpdata 
     135             :       type(t_hybdat), intent(in)   :: hybdat
     136             :       class(t_mat), intent(inout)  :: coulomb 
     137             :       class(t_mat), intent(inout)  :: coul_submtx
     138             : 
     139             :       integer :: i, j
     140             : 
     141          36 :       call timestart("copy out 2")
     142             : 
     143             :       SELECT TYPE(coulomb)
     144             :       CLASS is (t_mat)
     145           0 :          do j = 1, mpdata%n_g(ikpt)
     146           0 :             do i = 1, mpdata%n_g(ikpt)
     147           0 :                coulomb%data_c(hybdat%nbasp+i, hybdat%nbasp + j) = conjg(coul_submtx%data_c(j, i))
     148             :             enddo 
     149             :          enddo
     150             :       class is (t_mpimat)
     151             : #ifdef CPP_SCALAPACK
     152             :          select type(coul_submtx)
     153             :          class is (t_mpimat)
     154          36 :             call coul_submtx%transpose()
     155             :             ! copy coul_submtx to bottom right corner of coulomb
     156             :             !call pzgemr2d(m,              n,               a,                  ia, ja,             desca, 
     157             :             call pzgemr2d(mpdata%n_g(ikpt),mpdata%n_g(ikpt),coul_submtx%data_c, 1, 1, coul_submtx%blacsdata%blacs_desc,&
     158             :             !             b,             ib,            jb,             descb, ictxt)
     159          36 :                         coulomb%data_c, hybdat%nbasp+1, hybdat%nbasp+1, coulomb%blacsdata%blacs_desc, coulomb%blacsdata%blacs_desc(2))
     160             :          class default
     161           0 :             call judft_error("coul_submtx should also be mpimat")
     162             :          end select
     163             : #endif
     164             :       end select
     165          36 :       call timestop("copy out 2")
     166          36 :    end subroutine copy_out_2
     167         144 : end module m_apply_inverse_olap