1! 2! Copyright (C) 1996-2016 The SIESTA group 3! This file is distributed under the terms of the 4! GNU General Public License: see COPYING in the top directory 5! or http://www.gnu.org/copyleft/gpl.txt. 6! See Docs/Contributors.txt for a list of contributors. 7! 8! This code segment has been fully created by: 9! Nick Papior Andersen, 2013, nickpapior@gmail.com 10! Please conctact the author, prior to re-using this code. 11! * It has been heavily inspired by the original authors of the 12! Transiesta code (hence the references here are still remaining) * 13 14module m_ts_mumps_init 15 16 use precision, only : dp 17 18 implicit none 19 20#ifdef SIESTA__MUMPS 21 22 integer, public, save :: MUMPS_mem = 20 23 integer, public, save :: MUMPS_ordering = 7 24 integer, public, save :: MUMPS_block = -8 ! blocking factor 25 26 public :: read_ts_mumps 27 public :: mum_err 28 public :: init_MUMPS 29 public :: analyze_MUMPS 30 public :: prep_LHS 31 public :: prep_RHS_Eq 32 public :: prep_RHS_nEq 33 public :: insert_Self_Energies 34 35 private 36 37contains 38 39 subroutine read_ts_mumps( ) 40 41 use fdf, only : fdf_get, leqi 42 character(len=200) :: chars 43 44 MUMPS_mem = fdf_get('TS.MUMPS.Mem',20) 45 MUMPS_block = fdf_get('TS.MUMPS.BlockingFactor',112) 46 chars = fdf_get('TS.MUMPS.Ordering','auto') 47 if ( leqi(chars,'auto') ) then 48 MUMPS_ordering = 7 49 else if ( leqi(chars,'amd') ) then 50 MUMPS_ordering = 0 51 else if ( leqi(chars,'amf') ) then 52 MUMPS_ordering = 2 53 else if ( leqi(chars,'scotch') ) then 54 MUMPS_ordering = 3 55 else if ( leqi(chars,'pord') ) then 56 MUMPS_ordering = 4 57 else if ( leqi(chars,'metis') ) then 58 MUMPS_ordering = 5 59 else if ( leqi(chars,'qamd') ) then 60 MUMPS_ordering = 6 61 else 62 call die('Unknown MUMPS ordering.') 63 end if 64 65 66 end subroutine read_ts_mumps 67 68 subroutine init_MUMPS(mum,ID) 69#ifdef MPI 70 use mpi_siesta 71#endif 72 include 'zmumps_struc.h' 73 type(zMUMPS_STRUC), intent(inout) :: mum 74 integer, intent(in) :: ID 75 character(len=25) :: file 76 77 ! define processors and MATRIX-type 78 mum%SYM = 0 ! non-symmetric 79 mum%PAR = 1 ! sequential 80 mum%COMM = MPI_COMM_SELF ! only communicate with it-self 81 82 mum%JOB = -1 ! initialise 83 call zMUMPS(mum) 84 call mum_err(mum, 'MUMPS initialization had an error.') 85 86 ! print out control (verbosity) 87 mum%ICNTL(4) = 2 88 89 ! Setup the control parameters 90 mum%ICNTL(5) = 0 ! assembled format 91 ! The ordering of the matrix 92 mum%ICNTL(7) = MUMPS_ordering 93 94 ! We request a sparse right-hand side: 95 ! 20 == 1, we need not initialize RHS_SPARSE! 96 mum%ICNTL(20) = 1 97 mum%ICNTL(21) = 0 98 99 ! Allow memory increase handled by the user 100 mum%ICNTL(14) = MUMPS_Mem 101 102 ! Sets the blocking factor (opt for change in future versions) 103 ! Current MUMPS version: 4.10.0 104 mum%ICNTL(27) = MUMPS_block 105 106 ! Request specific elements of the inverse matrix: 107 ! 30 == 1 we MUST allocate it, but need not initialize it 108 mum%ICNTL(30) = 1 109 ! this however, posses some other problems. 110 111 ! For each processor, add their own output files. 112 ! As they become quite big we make them rewritten 113 ! in every SCF 114 write(file,'(a,i0,a)') 'TS_MUMPS_',ID,'.dat' 115 call io_assign(mum%ICNTL(1)) 116 mum%ICNTL(2) = mum%ICNTL(1) 117 mum%ICNTL(3) = mum%ICNTL(1) 118 open(mum%ICNTL(1),file=trim(file),status='replace', & 119 action='write',form='formatted') 120 121 end subroutine init_MUMPS 122 123 subroutine analyze_MUMPS(mum) 124 include 'zmumps_struc.h' 125 type(zMUMPS_STRUC), intent(inout) :: mum 126 integer :: iu 127 128 ! We force the analysis to be done without 129 ! knowing the values of A 130 ! otherwise A contains "spurious" values 131 mum%ICNTL(6) = 1 132 133 ! analyse the MUMPS solver, this will determine 134 ! factorization strategy 135 mum%JOB = 1 136 call zMUMPS(mum) 137 call mum_err(mum, 'MUMPS analysis step had an error.') 138 139 ! Write out estimated memory requirements 140 iu = mum%ICNTL(1) 141 write(iu,'(/,a)')'### Memory estimation from ANALYSIS step...' 142 write(iu,'(a,i0,a)')'### Minimum memory required for calculation: ', & 143 mum%INFO(15),' MB' 144 write(iu,'(a,i0,a)')'### MUMPS is allocating: ', & 145 nint(mum%INFO(15)*real(100+mum%ICNTL(14),dp)/100._dp),' MB' 146 write(iu,'(a,/)')"### MUMPS memory can be altered using TS.MUMPS.Mem." 147 148 end subroutine analyze_MUMPS 149 150 subroutine prep_LHS(IsVolt, mum,N_Elec,Elecs) 151#ifdef MPI 152 use mpi_siesta, only : MPI_Comm_Self 153#endif 154 use class_OrbitalDistribution 155 use class_Sparsity 156 use m_ts_elec_se, only: UC_minimum_worksize 157 use m_ts_sparse, only : ts_sp_uc 158 use m_ts_electype 159 use m_ts_method, only : orb_offset 160 use create_Sparsity_Union, only: crtSparsity_Union 161 include 'zmumps_struc.h' 162 163 logical, intent(in) :: IsVolt 164 type(zMUMPS_STRUC), intent(inout) :: mum 165 integer, intent(in) :: N_Elec 166 type(Elec), intent(inout) :: Elecs(N_Elec) 167 168 type(OrbitalDistribution) :: dit 169 type(Sparsity) :: tmpSp1, tmpSp2 170 integer :: iEl, idx, no, io, ind, nr, ioff 171 integer, pointer :: l_ptr(:),l_ncol(:), l_col(:) 172 173#ifdef MPI 174 call newDistribution(nrows_g(ts_sp_uc),MPI_Comm_Self,dit, & 175 name='MUMPS UC distribution') 176#else 177 call newDistribution(nrows_g(ts_sp_uc),-1,dit, & 178 name='MUMPS UC distribution') 179#endif 180 181 ! This works as creating a new sparsity deletes the previous 182 ! and as it is referenced several times it will not be actually 183 ! deleted... 184 tmpSp2 = ts_sp_uc 185 do iEl = 1 , N_Elec 186 187 idx = Elecs(iEl)%idx_o 188 no = TotUsedOrbs(Elecs(iEl)) 189 190 ! we first create the super-set sparsity 191 tmpSp1 = tmpSp2 192 call crtSparsity_Union(dit,tmpSp1, & 193 idx,idx,no,no, tmpSp2) 194 end do 195 call delete(tmpSp1) 196 call delete(dit) 197 198 ! Create the index 199 call attach(tmpSp2,list_ptr=l_ptr, & 200 n_col=l_ncol,list_col=l_col,nrows_g=nr, & 201 nnzs=mum%NZ) 202 call UC_minimum_worksize(IsVolt, N_Elec, Elecs, io) 203 io = max(mum%NZ, io) 204 205 ! Allocate LHS, first we need to 206 ! calculate the actual size of the matrix 207 ! We know that it must be the Hamiltonian sparsity 208 ! pattern UNION DENSE-electrodes! 209 call memory('A', 'I', mum%NZ * 2, 'prep_LHS') 210 allocate( mum%IRN( mum%NZ ) ) 211 allocate( mum%JCN( mum%NZ ) ) 212 call memory('A', 'Z', io, 'prep_LHS') 213 allocate( mum%A( io ) ) 214 215!$OMP parallel do default(shared), private(io,ioff,ind) 216 do io = 1 , nr 217 218 if ( l_ncol(io) /= 0 ) then 219 220 ioff = io - orb_offset(io) 221 222 do ind = l_ptr(io) + 1 , l_ptr(io) + l_ncol(io) 223 224 mum%IRN(ind) = l_col(ind) - orb_offset(l_col(ind)) 225 mum%JCN(ind) = ioff 226 227 end do 228 229 end if 230 end do 231!$OMP end parallel do 232 233 call delete(tmpSp2) 234 235 end subroutine prep_LHS 236 237 subroutine allocate_mum(mum,nsize,N_Elec,Elecs,GF) 238 use m_ts_electype 239 include 'zmumps_struc.h' 240 type(zMUMPS_STRUC), intent(inout) :: mum 241 integer, intent(in) :: nsize, N_Elec 242 type(Elec), intent(inout) :: Elecs(N_Elec) 243 complex(dp), pointer :: Gf(:) 244 245 integer :: no, iEl, io 246 247 ! We allocate for the equilibrium GF 248 mum%NZ_RHS = nsize ! number of non-zero RHS elements 249 250 if ( associated(mum%IRHS_PTR) ) then 251 call memory('D','I',size(mum%IRHS_PTR), 'allocate_mum') 252 deallocate(mum%IRHS_PTR) 253 call memory('D','I',size(mum%IRHS_SPARSE), 'allocate_mum') 254 deallocate(mum%IRHS_SPARSE) 255 nullify(mum%IRHS_PTR,mum%IRHS_SPARSE,mum%RHS_SPARSE) 256 end if 257 call memory('A','I',mum%NRHS+1, 'allocate_mum') 258 allocate( mum%IRHS_PTR(mum%NRHS+1) ) 259 call memory('A','I',mum%NZ_RHS, 'allocate_mum') 260 allocate( mum%IRHS_SPARSE(mum%NZ_RHS) ) 261 262 no = 0 263 do iEl = 1 , N_Elec 264 io = TotUsedOrbs(Elecs(iEl)) 265 no = no + io ** 2 266 end do 267 ! Allocate maximum space available 268 if ( associated(Gf) ) then 269 call memory('D','Z',size(Gf), 'allocate_mum') 270 deallocate(Gf) 271 nullify(Gf) 272 end if 273 ! for large electrodes and not so large device 274 ! we need to allocate more space 275 allocate( Gf(max(no,mum%NZ_RHS)) ) 276 call memory('A','Z',size(Gf), 'allocate_mum') 277 mum%RHS_SPARSE => Gf(1:mum%NZ_RHS) 278 279 no = 0 280 do iEl = 1 , N_Elec 281 282 ! This seems stupid, however, we never use the Sigma and 283 ! GF at the same time. Hence it will be safe 284 ! to have them point to the same array. 285 ! When the UC_expansion_Sigma_GammaT is called 286 ! first the Sigma is assigned and then 287 ! it is required that prepare_GF_inv is called 288 ! immediately (which it is) 289 ! Hence the GF must NOT be used in between these two calls! 290 io = TotUsedOrbs(Elecs(iEl)) ** 2 291 Elecs(iEl)%Sigma => Gf(no+1:no+io) 292 no = no + io 293 294 end do 295 296 end subroutine allocate_mum 297 298 subroutine prep_RHS_Eq(mum,no_u_TS,nzs,N_Elec,Elecs,GF) 299 use m_ts_electype 300 use m_ts_sparse, only : tsup_sp_uc 301 use m_ts_method, only : orb_offset, orb_type, TYP_BUFFER 302 use class_Sparsity 303 include 'zmumps_struc.h' 304 type(zMUMPS_STRUC), intent(inout) :: mum 305 integer, intent(in) :: no_u_TS, nzs, N_Elec 306 type(Elec), intent(inout) :: Elecs(N_Elec) 307 complex(dp), pointer :: Gf(:) 308 integer :: io, j, nr, ind 309 integer, pointer :: l_ptr(:), l_ncol(:), l_col(:) 310 311 ! Create the index 312 call attach(tsup_sp_uc,list_ptr=l_ptr, & 313 n_col=l_ncol,list_col=l_col,nrows_g=nr) 314 315 call allocate_mum(mum,nzs,N_Elec,Elecs,GF) 316 317 ! TODO, this requires that the sparsity pattern is symmetric 318 ! Which it always is! 319 mum%IRHS_PTR(:) = 1 320 321!$OMP parallel do default(shared), private(io,j,ind) 322 do io = 1 , nr 323 if ( orb_type(io) /= TYP_BUFFER ) then 324 mum%IRHS_PTR(io-orb_offset(io)) = l_ptr(io) + 1 325 if ( l_ncol(io) /= 0 ) then ! no entries 326 ! Create the row-index 327 do j = 1 , l_ncol(io) 328 ind = l_ptr(io) + j 329 mum%IRHS_SPARSE(ind) = l_col(ind) - orb_offset(l_col(ind)) 330 end do 331 332 end if 333 end if 334 335 end do 336!$OMP end parallel do 337 338 mum%IRHS_PTR(no_u_TS+1) = nzs + 1 339 340 end subroutine prep_RHS_Eq 341 342 subroutine prep_RHS_nEq(mum,no_u_TS,N_Elec, Elecs,Gf) 343 use m_ts_electype 344 use m_ts_method, only : ts2s_orb 345 use class_Sparsity 346 include 'zmumps_struc.h' 347 type(zMUMPS_STRUC), intent(inout) :: mum 348 349 integer, intent(in) :: no_u_TS, N_Elec 350 type(Elec), intent(inout) :: Elecs(N_Elec) 351 complex(dp), pointer :: Gf(:) 352 integer :: iEl, no, i, j, io, ind 353 354 ! We only need a partial size of the Green function 355 no = sum(TotUsedOrbs(Elecs)) 356 357 call allocate_mum(mum,no*no_u_TS,N_Elec,Elecs,GF) 358 359 ind = 0 360 do i = 1 , no_u_TS 361 mum%IRHS_PTR(i) = ind + 1 362 ! get correct siesta-orbital 363 io = ts2s_orb(i) 364 iElec: do iEl = 1 , N_Elec 365 if ( .not. OrbInElec(Elecs(iEl),io) ) cycle 366 ! Create the row-index 367 do j = 1 , no_u_TS 368 ind = ind + 1 369 mum%IRHS_SPARSE(ind) = j 370 end do 371 exit iElec 372 end do iElec 373 end do 374 mum%IRHS_PTR(no_u_TS+1) = no*no_u_TS + 1 375 376 if ( ind /= mum%NZ_RHS ) then 377 call die('Error in sparsity pattern of non-equilibrium') 378 end if 379 380 end subroutine prep_RHS_nEq 381 382 subroutine insert_Self_Energies(mum, El) 383 use m_ts_electype 384 use m_ts_method, only : orb_offset, ts2s_orb 385 include 'zmumps_struc.h' 386 type(zMUMPS_STRUC), intent(inout) :: mum 387 type(Elec), intent(in) :: El 388 389 integer :: off, ii, no, ind, iso, jso 390 391 no = TotUsedOrbs(El) 392 off = El%idx_o - 1 393 394 if ( El%Bulk ) then 395!$OMP do private(ind,jso,iso,ii) 396 do ind = 1 , mum%NZ 397 jso = ts2s_orb(mum%JCN(ind)) 398 if ( OrbInElec(El,jso) ) then 399 iso = ts2s_orb(mum%IRN(ind)) 400 if ( OrbInElec(El,iso) ) then 401 ii = (jso - El%idx_o) * no + iso - off 402 mum%A(ind) = El%Sigma(ii) 403 end if 404 end if 405 end do 406!$OMP end do nowait 407 else 408!$OMP do private(ind,jso,iso,ii) 409 do ind = 1 , mum%NZ 410 jso = ts2s_orb(mum%JCN(ind)) 411 if ( OrbInElec(El,jso) ) then 412 iso = ts2s_orb(mum%IRN(ind)) 413 if ( OrbInElec(El,iso) ) then 414 ii = (jso - El%idx_o) * no + iso - off 415 mum%A(ind) = mum%A(ind) - El%Sigma(ii) 416 end if 417 end if 418 end do 419!$OMP end do nowait 420 end if 421 422 end subroutine insert_Self_Energies 423 424 425 subroutine mum_err(mum,m) 426 include 'zmumps_struc.h' 427 type(zMUMPS_STRUC), intent(inout) :: mum 428 character(len=*), intent(in) :: m 429 430 ! We here check the error message that MUMPS 431 ! has given 432 if ( mum%INFOG(1) == 0 .and. mum%INFO(1) == 0 ) return ! no error 433 434 ! write out the message 435 call msg(m) 436 437 select case ( mum%INFOG(1) ) 438 case ( -5 ) 439 call die('Analysis step could not allocate space. & 440 &Do you have enough memory on your machine?') 441 case ( -8 , -14 ) 442 call msg('Integer work-array too small, increase & 443 &TS.MUMPS.Mem.') 444 write(*,'(a,i0)') 'Additional elements are needed: ',mum%INFOG(2) 445 write(0,'(a,i0)') 'Additional elements are needed: ',mum%INFOG(2) 446 case ( -9 ) 447 call msg('Work-array S too small, increase & 448 &TS.MUMPS.Mem.') 449 write(*,'(a,i0)') 'Additional elements are needed: ',-mum%INFOG(2)*1e6 450 write(0,'(a,i0)') 'Additional elements are needed: ',-mum%INFOG(2)*1e6 451 case ( -17, -20 ) 452 call msg('Work-array S too small, increase & 453 &TS.MUMPS.Mem.') 454 end select 455 456 call die('Check the MUMPS documentation for & 457 &the error message as well as the output.') 458 459 contains 460 461 subroutine msg(m) 462 character(len=*), intent(in) :: m 463 write(*,'(a)') m 464 write(0,'(a)') m 465 end subroutine msg 466 467 end subroutine mum_err 468 469 470#endif 471end module m_ts_mumps_init 472