Commit af3efe33 authored by Alberto Garcia's avatar Alberto Garcia
Browse files

Wrap PEXSI code within preprocessor blocks

The mechanism is similar to that used for TranSiesta.

(Also wrapped ELPA code within MPI preprocessor blocks)

parent a074e723
------------------------------------------------------------
July 30, 2016 A. Garcia trunk-537--pexsi-v0.8-spin-11
Wrap PEXSI code within preprocessor blocks
The mechanism is similar to that used for TranSiesta.
(Also wrapped ELPA code within MPI preprocessor blocks)
------------------------------------------------------------
July 29, 2016 A. Garcia trunk-537--pexsi-v0.8-spin-10
Remove old pexsi code
......
......@@ -410,17 +410,6 @@ $(ALL_OBJS): $(XC) $(COMP_LIBS)
# Libraries that might need to be compiled
# PEXSI dummy routines
libdummy_pexsi.a:
@(cd dummy_pexsi ; $(MAKE) \
"VPATH=$(VPATH)/dummy_pexsi" \
"CC=$(CC)" \
"CFLAGS=$(CFLAGS)" libdummy_pexsi.a)
libdummy_raw_pexsi.a:
@(cd dummy_pexsi ; $(MAKE) \
"VPATH=$(VPATH)/dummy_pexsi" \
"CC=$(CC)" \
"CFLAGS=$(CFLAGS)" libdummy_raw_pexsi.a)
#
# Linear algebra routines
#
......
......@@ -91,9 +91,11 @@ module class_Distribution
type(Distribution_) :: spdata
integer :: ierr
#ifdef MPI
if (spdata%group /= MPI_GROUP_NULL) then
call MPI_Group_free(spdata%group,ierr)
endif
#endif
if (allocated(spdata%ranks_in_ref_comm)) then
deallocate(spdata%ranks_in_ref_comm)
endif
......@@ -202,8 +204,12 @@ module class_Distribution
if (Node >= obj%Nodes) then
nl = 0
else
#ifdef MPI
nl = numroc(nels,obj%blocksize,Node, &
obj%isrcproc,obj%nodes)
#else
nl = nels
#endif
endif
case (TYPE_PEXSI)
remainder = nels - obj%blocksize * obj%nodes
......@@ -236,8 +242,12 @@ module class_Distribution
if (Node >= obj%Nodes) then
ig = 0
else
#ifdef MPI
ig = indxl2g(il,obj%blocksize,Node, &
obj%isrcproc,obj%nodes)
#else
ig = il
#endif
endif
case (TYPE_PEXSI)
if (Node >= obj%Nodes) then
......@@ -268,16 +278,24 @@ module class_Distribution
if (Node >= obj%nodes) then
il = 0
else
#ifdef MPI
il = indxg2l(ig,obj%blocksize,Node, &
obj%isrcproc,obj%nodes)
#else
il = ig
#endif
endif
else
! Assume that we only want a non-zero value if the orb
! is local to this node
owner = node_handling_element_(this,ig)
if (owner == obj%node) then
#ifdef MPI
il = indxg2l(ig,obj%blocksize,myrank, &
obj%isrcproc,obj%nodes)
#else
il = ig
#endif
else
il = 0
endif
......@@ -317,9 +335,12 @@ module class_Distribution
select case(obj%dist_type)
case (TYPE_BLOCK_CYCLIC)
#ifdef MPI
proc = indxg2p(ig,obj%blocksize,dummy_Node, &
obj%isrcproc,obj%nodes)
#else
proc = 0
#endif
case (TYPE_PEXSI)
! Assume bs=2, norbs=13, nodes=5
......
......@@ -34,10 +34,13 @@
use parallel, only: IONode
use parallel, only: SIESTA_worker
use m_compute_ebs_shift, only: compute_ebs_shift
#ifdef PEXSI
use m_pexsi_solver, only: pexsi_solver
#endif
use m_hsx, only : write_hs_formatted
c$$$ use m_matrixwrite
#ifdef MPI
use mpi_siesta
#endif
#ifdef CDF
use iodmhs_netcdf, only: write_dmh_netcdf
#endif
......@@ -105,8 +108,8 @@ c$$$ use m_matrixwrite
endif
endif
#ifdef PEXSI
if (isolve .eq. SOLVE_PEXSI) then
! Prototype interface
call pexsi_solver(iscf, no_u, no_l, nspin,
$ maxnh, numh, listhptr, listh,
$ H, S, qtot, Dscf, Escf,
......@@ -114,7 +117,7 @@ c$$$ use m_matrixwrite
if (ionode) write (6,"(/a,f14.6)") 'Entropy/k:', Entropy
endif
if (.not. SIESTA_worker) RETURN
#endif
! Here we decide if we want to calculate one or more SCF steps by
! diagonalization before proceeding with the OMM routine
CallDiagon=.false.
......
......@@ -6,7 +6,7 @@
! distributed along with the original code in the file "COPYING".
module ELPA1
#ifdef MPI
! Version 1.1.2, 2011-02-21
implicit none
......@@ -3804,14 +3804,14 @@ subroutine hh_transform_complex(alpha, xnorm_sq, xf, tau)
end subroutine
! --------------------------------------------------------------------------------------------------
#endif
end module ELPA1
! --------------------------------------------------------------------------------------------------
! Please note that the following routines are outside of the module ELPA1
! so that they can be used with real or complex data
! --------------------------------------------------------------------------------------------------
#ifdef MPI
subroutine elpa_transpose_vectors(vmat_s,ld_s,comm_s,vmat_t,ld_t,comm_t,nvs,nvr,nvc,nblk)
!-------------------------------------------------------------------------------
......@@ -4014,5 +4014,5 @@ subroutine elpa_reduce_add_vectors(vmat_s,ld_s,comm_s,vmat_t,ld_t,comm_t,nvr,nvc
deallocate(aux2)
end subroutine
#endif
!-------------------------------------------------------------------------------
......@@ -6,7 +6,7 @@
! distributed along with the original code in the file "COPYING".
module ELPA2
#ifdef MPI
! Version 1.1.2, 2011-02-21
USE ELPA1
......@@ -3408,5 +3408,5 @@ subroutine divide_band(nblocks_total, n_pes, block_limits)
endif
end subroutine
#endif
end module ELPA2
! For consistency, include the actual interface file in the library directory
#ifdef PEXSI
!
include "f_interface.f90"
#endif
module m_pexsi
#ifdef PEXSI
use precision, only: dp
use iso_c_binding
......@@ -67,5 +68,5 @@ module m_pexsi
ncol = ncol - 1
enddo
end subroutine get_row_col
#endif
end module m_pexsi
! Tangled code
module m_pexsi_dos
#ifdef PEXSI
use precision, only : dp
public :: pexsi_dos
......@@ -463,5 +464,6 @@ character(len=*), intent(in) :: str
end subroutine check_info
end subroutine pexsi_dos
#endif
end module m_pexsi_dos
! End of tangled code
! Tangled code
module m_pexsi_solver
use precision, only : dp
implicit none
public :: pexsi_solver
real(dp), save :: prevDmax ! For communication of max diff in DM in scf loop
! used in the heuristics for N_el tolerance
public :: prevDmax
#ifdef PEXSI
public :: pexsi_solver
CONTAINS
......@@ -1229,5 +1230,6 @@ character(len=*), intent(in) :: str
end subroutine check_info
end subroutine pexsi_solver
#endif
end module m_pexsi_solver
! End of tangled code
This diff is collapsed.
! --- Tangled code
module m_redist_spmatrix
#ifdef PEXSI
implicit none
type, public :: comm_t
integer :: src, dst, i1, i2, nitems
......@@ -542,5 +543,6 @@ CONTAINS
deallocate(local_reqR, local_reqS, statuses)
end subroutine do_transfers_dp
#endif
end module m_redist_spmatrix
! --- End of tangled code
subroutine memory_all(str,comm)
use m_rusage, only : rss_max
#ifdef MPI
use mpi
#endif
character(len=*), intent(in) :: str
integer, intent(in) :: comm
......@@ -19,9 +20,14 @@ subroutine memory_all(str,comm)
my_mem = rss_max()
#ifdef MPI
call MPI_Reduce(my_mem,max_mem,1,MPI_Real,MPI_max,0,comm,MPIerror)
call MPI_Reduce(my_mem,min_mem,1,MPI_Real,MPI_min,0,comm,MPIerror)
#else
max_mem = my_mem
min_mem = my_mem
#endif
if (myrank == 0) then
write(6,"(a,2f12.2)") " &m -- Peak memory (Mb) " // trim(str) // " (max,min): ", max_mem, min_mem
endif
......
......@@ -21,6 +21,7 @@ This is the main structure of the code.
#+begin_src f90 :noweb-ref code-structure
module m_pexsi_dos
#ifdef PEXSI
use precision, only : dp
public :: pexsi_dos
......@@ -48,6 +49,7 @@ This is the main structure of the code.
<<support-routines>>
end subroutine pexsi_dos
#endif
end module m_pexsi_dos
#+end_src
......
......@@ -10,12 +10,14 @@ Module structure:
#+BEGIN_SRC f90 :noweb-ref module-structure
MODULE m_pexsi_local_dos
#ifdef PEXSI
private
public :: pexsi_local_dos
CONTAINS
<<siesta-side-parent-routine>>
<<pexsi-ldos-routine>>
#endif
end module m_pexsi_local_dos
#+END_SRC
......@@ -102,9 +104,6 @@ This is the main structure of the LDOS routine.
<<routine-header>>
<<routine-variables>>
! -------- for serial compilation
#ifndef MPI
call die("PEXSI LDOS needs MPI")
#else
<<define-communicators>>
<<re-distribute-matrices>>
<<set-options>>
......@@ -114,7 +113,6 @@ This is the main structure of the LDOS routine.
<<compute-ldos>>
<<copy-to-siesta-side>>
<<clean-up>>
#endif
CONTAINS
......@@ -144,19 +142,17 @@ end subroutine get_LDOS_SI
*** Used modules
#+BEGIN_SRC f90 :noweb-ref used-modules
use precision, only : dp
use fdf
use units, only: eV, pi
use sys, only: die
use m_mpi_utils, only: broadcast
use parallel, only : SIESTA_worker, BlockSize
use parallel, only : SIESTA_Group, SIESTA_Comm
use m_redist_spmatrix, only: aux_matrix, redistribute_spmatrix
use class_Distribution
use alloc, only: re_alloc, de_alloc
#ifdef MPI
use precision, only : dp
use fdf
use units, only: eV, pi
use sys, only: die
use m_mpi_utils, only: broadcast
use parallel, only : SIESTA_worker, BlockSize
use parallel, only : SIESTA_Group, SIESTA_Comm
use m_redist_spmatrix, only: aux_matrix, redistribute_spmatrix
use class_Distribution
use alloc, only: re_alloc, de_alloc
use mpi_siesta
#endif
use iso_c_binding
use f_ppexsi_interface, only: f_ppexsi_options
use f_ppexsi_interface, only: f_ppexsi_plan_finalize
......
......@@ -15,8 +15,9 @@ factorization information.
** TODO Maybe integrate the plan handling in the =pexsi_solver= module.
#+BEGIN_SRC f90 :noweb yes :tangle m_pexsi.f90
#+BEGIN_SRC f90 :noweb yes :tangle m_pexsi.F90
module m_pexsi
#ifdef PEXSI
use precision, only: dp
use iso_c_binding
......@@ -32,7 +33,7 @@ factorization information.
<<initialize-scf-loop>>
<<finalize-scf-loop>>
<<get-optimal-row-col>>
#endif
end module m_pexsi
#+END_SRC
......@@ -110,15 +111,16 @@ This is the main structure of the code.
#+begin_src f90 :noweb-ref code-structure
module m_pexsi_solver
use precision, only : dp
implicit none
public :: pexsi_solver
real(dp), save :: prevDmax ! For communication of max diff in DM in scf loop
! used in the heuristics for N_el tolerance
public :: prevDmax
#ifdef PEXSI
public :: pexsi_solver
CONTAINS
......@@ -146,6 +148,7 @@ CONTAINS
<<support-routines>>
end subroutine pexsi_solver
#endif
end module m_pexsi_solver
#+end_src
......
......@@ -43,12 +43,12 @@
use parallel, only : Node, Nodes, BlockSize
#ifdef MPI
use mpi_siesta, only : mpi_comm_world
use elpa1, only: get_elpa_row_col_comms
use elpa2, only: solve_evp_real_2stage
#endif
use alloc
use sys, only : die
use elpa1, only: get_elpa_row_col_comms
use elpa2, only: solve_evp_real_2stage
implicit none
......
......@@ -695,11 +695,16 @@ subroutine read_options( na, ns, nspin )
write(6,3) 'redata: Method of Calculation', 'Orbital Minimization Method'
endif
else if (leqi(method,"pexsi")) then
isolve = SOLVE_PEXSI
#ifdef PEXSI
isolve = SOLVE_PEXSI
if (ionode) then
write(6,'(a,4x,a)') 'redata: Method of Calculation = ', &
'PEXSI'
endif
#else
call die("PEXSI solver is not compiled in. Use -DPEXSI")
#endif
#ifdef TRANSIESTA
else if (leqi(method,'transi') .or. leqi(method,'transiesta') ) then
isolve = SOLVE_TRANSI
......
......@@ -25,12 +25,14 @@ non-disjoint groups.
The general structure is
#+BEGIN_SRC f90 :noweb-ref code-structure
module m_redist_spmatrix
#ifdef PEXSI
implicit none
<<module-type-declarations>>
public :: redistribute_spmatrix
CONTAINS
<<routine-redist-spmatrix>>
<<transfers>>
#endif
end module m_redist_spmatrix
#+END_SRC
......
......@@ -17,7 +17,9 @@
use parallel, only: SIESTA_worker ! whether part of Siesta's communicator
use parallel, only: ionode
#ifdef MPI
use mpi_siesta, only: true_mpi_comm_world
#endif
use m_mpi_utils, only: broadcast
#ifdef TRACING_SOLVEONLY
......@@ -58,18 +60,20 @@
C Begin of coordinate relaxation iteration
relaxd = .false.
#ifdef PEXSI
! Broadcast relevant things for program logic
! These were set in siesta_options, called only by "SIESTA_workers".
call broadcast(inicoor,comm=true_MPI_Comm_World)
call broadcast(fincoor,comm=true_MPI_Comm_World)
#endif
istep = inicoor
DO WHILE ((istep.le.fincoor) .AND. (.not. relaxd))
call siesta_forces( istep )
if (SIESTA_worker) call siesta_move( istep, relaxd )
#ifdef PEXSI
call broadcast(relaxd,comm=true_MPI_Comm_World)
#endif
if (.not. relaxd) then
istep = istep + 1
endif
......
......@@ -21,8 +21,10 @@
use m_ksv
USE m_projected_DOS, only: projected_DOS
USE m_local_DOS, only: local_DOS
#ifdef PEXSI
USE m_pexsi_local_DOS, only: pexsi_local_DOS
USE m_pexsi_dos, only: pexsi_dos
#endif
USE siesta_options
use units, only: Debye, eV
use sparse_matrices, only: maxnh, listh, listhptr, numh
......@@ -124,11 +126,13 @@
endif
endif
#ifdef PEXSI
if (fdf_get("PEXSI.DOS",.false.)) then
call pexsi_dos(no_u, no_l, nspin,
$ maxnh, numh, listhptr, listh, H, S, qtot, ef)
endif
#endif
! section done by Siesta subset of nodes
if (SIESTA_worker) then
......@@ -546,11 +550,11 @@ C Find local density of states
! in 'siesta_move'
endif ! SIESTA_worker
#ifdef PEXSI
if (fdf_get("PEXSI.LocalDOS",.false.)) then
call pexsi_local_DOS()
endif
#endif
if (SIESTA_worker) call timer("Analysis",2)
!------------------------------------------------------------------------- END
......
......@@ -14,7 +14,9 @@
CONTAINS
subroutine siesta_forces(istep)
#ifdef MPI
use mpi_siesta
#endif
use precision, only: dp
use siesta_cml
#ifdef SIESTA__FLOOK
......@@ -50,9 +52,7 @@
use siesta_master, only: siesta_server ! Is siesta a server?
use m_save_density_matrix, only: save_density_matrix
use m_iodm_old, only: write_spmatrix
#ifdef MPI
use mpi_siesta, only: MPI_Barrier, MPI_Comm_World
#endif
!
use units, only: eV, Ang
use sparse_matrices, only: H, Hold, Dold, Dscf
use m_pexsi_solver, only: prevDmax
......@@ -68,7 +68,9 @@
use m_compute_energies, only: compute_energies
use m_mpi_utils, only: broadcast
use fdf
#ifdef PEXSI
use m_pexsi, only: pexsi_finalize_scfloop
#endif
#ifdef TRANSIESTA
use m_ts_options, only : N_Elec
......@@ -113,10 +115,12 @@
call write_debug( ' PRE siesta_forces' )
#endif
#ifdef PEXSI
! Broadcast relevant things for program logic
! These were set in read_options, called only by "SIESTA_workers".
call broadcast(nscf,comm=true_MPI_Comm_World)
#endif
if (SIESTA_worker) then
! Initialization tasks for a given geometry
call state_init( istep )
......@@ -368,17 +372,19 @@
endif ! SIESTA_worker
#ifdef PEXSI
call broadcast(last_step,comm=true_MPI_Comm_World)
call broadcast(iscf,comm=true_MPI_Comm_World)
#endif
if ( last_step ) exit
end do
#ifdef PEXSI
if (isolve == SOLVE_PEXSI) then
call pexsi_finalize_scfloop()
endif
#endif
if (.not. SIESTA_worker) RETURN
......
trunk-537--pexsi-v0.8-spin-10
trunk-537--pexsi-v0.8-spin-11
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment