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