prg_graphsolver_mod.F90

Go to the documentation of this file.

 
module prg_graphsolver_mod
 
  use bml
  use prg_parallel_mod
  use prg_progress_mod
  use prg_partition_mod
  use prg_graph_mod
  use prg_system_mod
  use prg_timer_mod
  use prg_extras_mod
  use prg_densitymatrix_mod
  use prg_subgraphloop_mod
  use prg_genz_mod
 
  implicit none
 
  private  !Everything is private by default
 
  integer, parameter :: dp = kind(1.0d0)
 
  public :: prg_build_densitygp_t0, prg_build_zmatgp
 
contains
 
  subroutine prg_build_densitygp_t0(ham_bml, g_bml, rho_bml, threshold, bndfil, &
       & Ef, nparts, verbose)
 
    integer, parameter ::  dp = kind(1.0d0)
    type(bml_matrix_t), intent(in) :: ham_bml, g_bml
    type(bml_matrix_t), intent(out) :: rho_bml
    real(dp) :: bndfil, threshold
    integer ::  myrank, norbs, nnodes, tnnz, i, vsize(2), inorbs
    integer, intent(inout) :: nparts
    real(dp), intent(inout) :: ef
    character(20) :: bml_type
    integer, allocatable :: xadj(:), adjncy(:), vector(:)
    integer, allocatable  ::  part(:), core_count(:), halo_count(:,:), ch_count(:)
    type(graph_partitioning_t) ::  gpat
    real(dp) :: maxch, pnorm=6, smooth_maxch, sumcubes, mlsi, mlsii
    type(system_type), allocatable    ::  syprt(:)
    integer, optional, intent(in) :: verbose
 
    real(dp) :: error, dnocc_part, dnocc, nocc
 
    ! Initialize progress MPI and get ranks
    if(getnranks() == 0) call prg_progress_init()
    !call prg_progress_init()
    myrank = getmyrank() + 1
    bml_type = bml_get_type(ham_bml)
 
    ! Allocate bml matrices
    norbs = bml_get_n(ham_bml)
    if(.not. bml_allocated(rho_bml))then
      call bml_zero_matrix(bml_type,bml_element_real,dp,norbs,norbs,rho_bml)
    endif
 
    ! Using Metis for graph partitioning
    allocate(xadj(norbs+1)) ! Adjacency in csr format
 
    tnnz = 0
    do i=1,norbs
      tnnz = tnnz + bml_get_row_bandwidth(g_bml,i)
    enddo
 
    allocate(adjncy(tnnz+1))
 
    call bml_adjacency(g_bml, xadj, adjncy, 1)
 
    nnodes = norbs
 
    allocate(part(nnodes))
    allocate(core_count(nparts))
    allocate(ch_count(nparts))
    allocate(halo_count(nparts, nnodes))
 
    call prg_metispartition(gpat, nnodes, nnodes, xadj, adjncy, nparts, part, core_count,&
         ch_count, halo_count, sumcubes, maxch, smooth_maxch, pnorm)
 
    if(present(verbose))then
      if(verbose > 1 .and. myrank == 1)then
        write(*,*)"gpat = ",gpat%pname
        write(*,*)"totalParts = ",gpat%totalParts
        write(*,*)"totalNodes = ",gpat%totalNodes
        write(*,*)"partMin = ",gpat%localPartMin(1)
        write(*,*)"partMax = ",gpat%localPartMax(1)
      endif
    endif
 
    allocate(syprt(gpat%TotalParts))
 
    call prg_wait
 
    if(present(verbose) .and. (verbose > 1 .and. myrank == 1)) mlsi = mls()
 
    dnocc = 0.0_dp
    ! Extract core halo indices from partitions
    do i=1,gpat%TotalParts
      call bml_matrix2submatrix_index(g_bml,&
           gpat%sgraph(i)%nodeInPart,gpat%nnodesInPart(i),&
           gpat%sgraph(i)%core_halo_index, &
           vsize,.true.)
      gpat%sgraph(i)%lsize = vsize(1)
      gpat%sgraph(i)%llsize = vsize(2)
      if(present(verbose))then
        if(verbose > 1 .and. myrank == 1)then
          write(*,*) 'nodeInpart size', size(gpat%sgraph(i)%nodeInPart)
          write(*,*) 'nodeInpart(i)', gpat%nnodesInPart(i)
          write(*,*)"Core and CH sizes:",vsize(2),vsize(1)
        endif
      endif
    enddo
 
    ! Construct DM for every part
    do i = gpat%localPartMin(myrank), gpat%localPartMax(myrank)
      !  do i=1,gpat%TotalParts ! If no MPI
      call bml_matrix2submatrix_index(g_bml,&
           gpat%sgraph(i)%nodeInPart,gpat%nnodesInPart(i),&
           gpat%sgraph(i)%core_halo_index, &
           vsize,.true.)
      gpat%sgraph(i)%lsize = vsize(1)
      gpat%sgraph(i)%llsize = vsize(2)
      inorbs = vsize(1)
      call bml_zero_matrix("dense",bml_element_real,dp,inorbs,inorbs,syprt(i)%estr%ham)
      if(allocated(vector))deallocate(vector)
      allocate(vector(gpat%sgraph(i)%lsize))
      vector(:) = gpat%sgraph(i)%core_halo_index(1:gpat%sgraph(i)%lsize)
      call bml_matrix2submatrix(ham_bml, syprt(i)%estr%ham, vector, &
           & inorbs)
 
      !Computing the density matrix with diagonalization
      mlsii = mls()
      call bml_zero_matrix("dense",bml_element_real,dp,inorbs,inorbs,syprt(i)%estr%rho)
      call prg_build_density_t_fermi(syprt(i)%estr%ham,syprt(i)%estr%rho,bndfil, 0.1_dp, ef, 0, dnocc_part)
      dnocc = dnocc + dnocc_part
 
      !call bml_print_matrix("rho_part_bml",syprt(i)%estr%rho,0,10,0,10)
      if(present(verbose))then
        if(verbose > 1 .and. myrank == 1) write(*,*)"Time for prg_build_density_T0",mls()-mlsii
      endif
      call bml_submatrix2matrix(syprt(i)%estr%rho,rho_bml,vector,gpat%sgraph(i)%lsize,gpat%sgraph(i)%llsize,threshold)
    enddo
 
    ! Collect the small DM matrices into its full system corresponding
    call prg_partordering(gpat)
    call prg_collectmatrixfromparts(gpat, rho_bml)
    call prg_wait
 
    nocc = bml_trace(rho_bml)
    error = abs(nocc-2*bndfil*norbs)
    if (nocc < 2*bndfil*norbs) dnocc = - dnocc
 
    if(present(verbose))then
      if(verbose > 1 .and. myrank == 1) write(6,'(A,10es15.6)') 'Ef, Nocc, delta_nocc, gradient',  &
           ef, nocc, error, dnocc
    endif
 
    if (error > 1.e-6_dp .and. abs(dnocc) > 1.e-6_dp) then
      if ( error > 1.0_dp ) then
        ef = ef - nocc/dnocc
      else
        ef = ef - error / dnocc
      endif
    endif
 
    if(present(verbose))then
      if(verbose > 1 .and. myrank == 1) write(*,*)"Total time for graph solver =",mls()-mlsi
    endif
 
    do i = gpat%localPartMin(myrank), gpat%localPartMax(myrank)
      call bml_deallocate(syprt(i)%estr%ham)
      call bml_deallocate(syprt(i)%estr%rho)
    enddo
    deallocate(syprt)
 
    !call prg_progress_shutdown
    if(getnranks() == 0) call prg_progress_shutdown()
  end subroutine prg_build_densitygp_t0
 
 
  subroutine prg_build_zmatgp(over_bml, g_bml, zmat_bml, threshold, &
       & nparts, verbose)
 
    integer, parameter ::  dp = kind(1.0d0)
    type(bml_matrix_t), intent(in) :: over_bml, g_bml
    type(bml_matrix_t), intent(out) :: zmat_bml
    real(dp) :: threshold
    integer ::  myrank, norbs, nnodes, tnnz, i, vsize(2), inorbs
    integer, intent(inout) :: nparts
    character(20) :: bml_type
    integer, allocatable :: xadj(:), adjncy(:), vector(:)
    integer, allocatable  ::  part(:), core_count(:), halo_count(:,:), ch_count(:)
    type(graph_partitioning_t) ::  gpat
    real(dp) :: maxch, pnorm=6, smooth_maxch, sumcubes, mlsi, mlsii
    type(system_type), allocatable    ::  syprt(:)
    integer, optional, intent(in) :: verbose
 
    ! Initialize progress MPI and get ranks
    if(getnranks() == 0) call prg_progress_init()
    !call prg_progress_init()
    myrank = getmyrank() + 1
 
    bml_type = bml_get_type(over_bml)
 
    ! Allocate bml matrices
    norbs = bml_get_n(over_bml)
    if(.not. bml_allocated(zmat_bml))then
      call bml_zero_matrix(bml_type,bml_element_real,dp,norbs,norbs,zmat_bml)
    endif
 
    ! Using Metis to do graph partitioning
    allocate(xadj(norbs+1)) ! Adjacency in csr format
 
    tnnz = 0
    do i=1,norbs
      tnnz = tnnz + bml_get_row_bandwidth(g_bml,i)
    enddo
 
    allocate(adjncy(tnnz+1))
 
    call bml_adjacency(g_bml, xadj, adjncy, 1)
 
    nnodes = norbs
 
    allocate(part(nnodes))
    allocate(core_count(nparts))
    allocate(ch_count(nparts))
    allocate(halo_count(nparts, nnodes))
 
    call prg_metispartition(gpat, nnodes, nnodes, xadj, adjncy, nparts, part, core_count,&
         ch_count, halo_count, sumcubes, maxch, smooth_maxch, pnorm)
 
    if(present(verbose))then
      if(verbose > 1 .and. myrank == 1)then
        write(*,*)"gpat = ",gpat%pname
        write(*,*)"totalParts = ",gpat%totalParts
        write(*,*)"totalNodes = ",gpat%totalNodes
        write(*,*)"partMin = ",gpat%localPartMin(1)
        write(*,*)"partMax = ",gpat%localPartMax(1)
      endif
    endif
    allocate(syprt(gpat%TotalParts))
 
    call prg_wait
 
    mlsi = mls()
 
    ! Extract core halo indices from partitions
    do i=1,gpat%TotalParts
      call bml_matrix2submatrix_index(g_bml,&
           gpat%sgraph(i)%nodeInPart,gpat%nnodesInPart(i),&
           gpat%sgraph(i)%core_halo_index, &
           vsize,.true.)
      gpat%sgraph(i)%lsize = vsize(1)
      gpat%sgraph(i)%llsize = vsize(2)
      if(present(verbose))then
        if(verbose > 1 .and. myrank == 1)then
          write(*,*) 'nodeInpart size', size(gpat%sgraph(i)%nodeInPart)
          write(*,*) 'nodeInpart(i)', gpat%nnodesInPart(i)
          write(*,*)"Core and CH sizes:",vsize(2),vsize(1)
        endif
      endif
    enddo
 
    ! Construct Zmat for every part
    do i = gpat%localPartMin(myrank), gpat%localPartMax(myrank)
      !  do i=1,gpat%TotalParts ! If no MPI
      call bml_matrix2submatrix_index(g_bml,&
           gpat%sgraph(i)%nodeInPart,gpat%nnodesInPart(i),&
           gpat%sgraph(i)%core_halo_index, &
           vsize,.true.)
      gpat%sgraph(i)%lsize = vsize(1)
      gpat%sgraph(i)%llsize = vsize(2)
      inorbs = vsize(1)
      call bml_zero_matrix("dense",bml_element_real,dp,inorbs,inorbs,syprt(i)%estr%over)
      if(allocated(vector))deallocate(vector)
      allocate(vector(gpat%sgraph(i)%lsize))
      vector(:) = gpat%sgraph(i)%core_halo_index(1:gpat%sgraph(i)%lsize)
      call bml_matrix2submatrix(over_bml, syprt(i)%estr%over, vector, &
           & inorbs)
 
      !Computing the zmatrix with diagonalization
      mlsii = mls()
      call bml_zero_matrix("dense",bml_element_real,dp,inorbs,inorbs,syprt(i)%estr%zmat)
      call prg_buildzdiag(syprt(i)%estr%over,syprt(i)%estr%zmat,threshold,inorbs,"dense")
 
      if(present(verbose))then
        if(verbose > 1 .and. myrank == 1) write(*,*)"Time for prg_build_density_T0",mls()-mlsii
      endif
      call bml_submatrix2matrix(syprt(i)%estr%zmat,zmat_bml,vector,gpat%sgraph(i)%lsize,gpat%sgraph(i)%llsize,threshold)
    enddo
 
    ! Collect the small DM matrices into its full system corresponding
    call prg_partordering(gpat)
    call prg_collectmatrixfromparts(gpat, zmat_bml)
    call prg_wait
 
    if(present(verbose))then
      if(verbose > 1 .and. myrank == 1)write(*,*)"Total time for graph solver =",mls()-mlsi
    endif
 
    do i = gpat%localPartMin(myrank), gpat%localPartMax(myrank)
      call bml_deallocate(syprt(i)%estr%over)
      call bml_deallocate(syprt(i)%estr%zmat)
    enddo
    deallocate(syprt)
 
    !call prg_progress_shutdown
    if(getnranks() == 0) call prg_progress_shutdown()
  end subroutine prg_build_zmatgp
 
end module prg_graphsolver_mod