Merge lp:~garsua/siesta/trunk-elsi-dm into lp:~albertog/siesta/trunk-elsi-dm

[Alberto Garcia]

Cosmetic changes on m_elsi_interface.F90
- Removed "use fdf, only: fdf_get" from elsi_real_solver and elsi_complex_solver subroutines
- Unindented lines from 401 to 510 (in elsi_real_solver) and from 1351 to 1460 (in elsi_complex_solver)

Trying to resubmit.

Unmerged revisions

696. By Victor GS on 2019-01-16

Cosmetic changes on m_elsi_interface.F90
- Removed "use fdf, only: fdf_get" from elsi_real_solver and elsi_complex_solver subroutines
- Unindented lines from 401 to 510 (in elsi_real_solver) and from 1351 to 1460 (in elsi_complex_solver)

=== modified file 'Src/m_elsi_interface.F90'
--- Src/m_elsi_interface.F90	2019-01-15 14:46:49 +0000
+++ Src/m_elsi_interface.F90	2019-01-17 09:06:38 +0000
@@ -271,7 +271,6 @@
 subroutine elsi_real_solver(iscf, n_basis, n_basis_l, n_spin, nnz_l, row_ptr, &
   col_idx, qtot, temp, ham, ovlp, dm, edm, ef, ets, Get_EDM_Only)
 
-  use fdf,         only: fdf_get
   use m_mpi_utils, only: globalize_sum
   use parallel,    only: BlockSize
 #ifdef MPI
@@ -399,116 +398,116 @@
 
  endif   ! iscf == 1
 
-    if (n_spin == 1) then
-
-       ! Sparsity pattern
-       call globalize_sum(nnz_l, nnz_g, comm=elsi_global_comm)
-
-       allocate(row_ptr2(n_basis_l+1))
-       row_ptr2(1:n_basis_l) = row_ptr(1:n_basis_l)+1
-       row_ptr2(n_basis_l+1) = nnz_l+1
-
-       call elsi_set_csc(elsi_h, nnz_g, nnz_l, n_basis_l, col_idx, row_ptr2)
-       deallocate(row_ptr2)
-
-       call elsi_set_csc_blk(elsi_h, BlockSize)
-       call elsi_set_mpi(elsi_h, elsi_global_comm)
-
-    else
-
-       ! MPI logic for spin polarization
-
-       ! Re-create numh, as we use it in the transfer
-       allocate(numh(n_basis_l))
-       numh(1) = row_ptr(2)
-       do i = 2, n_basis_l-1
-          numh(i) = row_ptr(i+1)-row_ptr(i)
-       enddo
-       numh(n_basis_l) = nnz_l - row_ptr(n_basis_l)
-
-       ! Split the communicator in spins and get distribution objects
-       ! for the data redistribution needed
-       ! Note that dist_spin is an array
-       call get_spin_comms_and_dists(elsi_global_comm,elsi_global_comm, &
-            blocksize, n_spin, &
-            dist_global,dist_spin, elsi_spatial_comm, elsi_spin_comm)
-
-       ! Find out which spin team we are in, and tag the spin we work on
-       call mpi_comm_rank( elsi_Spin_Comm, spin_rank, ierr )
-       my_spin = spin_rank+1  ! {1,2}
-
-
-       ! This is done serially, each time filling one spin set
-       ! Note that **all processes** need to have the same pkg_global
-
-       do ispin = 1, n_spin
-
-          ! Load pkg_global data package
-          pkg_global%norbs = n_basis
-          pkg_global%no_l  = n_basis_l
-          pkg_global%nnzl  = nnz_l
-          pkg_global%numcols => numh
-          pkg_global%cols    => col_idx
-
-          allocate(pkg_global%vals(2))
-          ! Link the vals items to the appropriate arrays (no extra memory here)
-          pkg_global%vals(1)%data => ovlp(:)
-          ! Note that we *cannot* say  => ham(:,my_spin)
-          ! and avoid the sequential loop, as then half the processors will send
-          ! the information for 'spin up' and the other half the information for 'spin down',
-          ! which is *not* what we want.
-          pkg_global%vals(2)%data => ham(:,ispin)
-
-          call timer("redist_orbs_fwd", 1)
-
-          ! We are doing the transfers sequentially. One spin team is
-          ! 'idle' (in the receiving side) in each pass, as the dist_spin(ispin) distribution
-          ! does not involve them.
-
-          call redistribute_spmatrix(n_basis,pkg_global,dist_global, &
-                                             pkg_spin,dist_spin(ispin),elsi_global_Comm)
-
-          call timer("redist_orbs_fwd", 2)
-
-          if (my_spin == ispin) then  ! Each team gets their own data
-
-             !nrows = pkg_spin%norbs          ! or simply 'norbs'
-             my_no_l = pkg_spin%no_l
-             my_nnz_l    = pkg_spin%nnzl
-             call MPI_AllReduce(my_nnz_l,my_nnz,1,MPI_integer,MPI_sum,elsi_Spatial_Comm,ierr)
-             ! generate off-by-one row pointer
-             call re_alloc(my_row_ptr2,1,my_no_l+1,"my_row_ptr2","elsi_solver")
-             my_row_ptr2(1) = 1
-             do ih = 1,my_no_l
-                my_row_ptr2(ih+1) = my_row_ptr2(ih) + pkg_spin%numcols(ih)
-             enddo
-
-             my_col_idx => pkg_spin%cols
-             my_S => pkg_spin%vals(1)%data
-             my_H => pkg_spin%vals(2)%data
-
-             call re_alloc(my_DM,1,my_nnz_l,"my_DM","elsi_solver")
-             call re_alloc(my_EDM,1,my_nnz_l,"my_EDM","elsi_solver")
-          endif
-
-          ! Clean pkg_global
-          nullify(pkg_global%vals(1)%data)
-          nullify(pkg_global%vals(2)%data)
-          deallocate(pkg_global%vals)
-          nullify(pkg_global%numcols)
-          nullify(pkg_global%cols)
-
-       enddo
-
-       call elsi_set_csc(elsi_h, my_nnz, my_nnz_l, my_no_l, my_col_idx, my_row_ptr2)
-       call de_alloc(my_row_ptr2,"my_row_ptr2","elsi_solver")
-
-       call elsi_set_csc_blk(elsi_h, BlockSize)
-       call elsi_set_spin(elsi_h, n_spin, my_spin)
-       call elsi_set_mpi(elsi_h, elsi_Spatial_comm)
-       call elsi_set_mpi_global(elsi_h, elsi_global_comm)
-
-    endif  ! n_spin
+ if (n_spin == 1) then
+
+    ! Sparsity pattern
+    call globalize_sum(nnz_l, nnz_g, comm=elsi_global_comm)
+
+    allocate(row_ptr2(n_basis_l+1))
+    row_ptr2(1:n_basis_l) = row_ptr(1:n_basis_l)+1
+    row_ptr2(n_basis_l+1) = nnz_l+1
+
+    call elsi_set_csc(elsi_h, nnz_g, nnz_l, n_basis_l, col_idx, row_ptr2)
+    deallocate(row_ptr2)
+
+    call elsi_set_csc_blk(elsi_h, BlockSize)
+    call elsi_set_mpi(elsi_h, elsi_global_comm)
+
+ else
+
+    ! MPI logic for spin polarization
+
+    ! Re-create numh, as we use it in the transfer
+    allocate(numh(n_basis_l))
+    numh(1) = row_ptr(2)
+    do i = 2, n_basis_l-1
+       numh(i) = row_ptr(i+1)-row_ptr(i)
+    enddo
+    numh(n_basis_l) = nnz_l - row_ptr(n_basis_l)
+
+    ! Split the communicator in spins and get distribution objects
+    ! for the data redistribution needed
+    ! Note that dist_spin is an array
+    call get_spin_comms_and_dists(elsi_global_comm,elsi_global_comm, &
+         blocksize, n_spin, &
+         dist_global,dist_spin, elsi_spatial_comm, elsi_spin_comm)
+
+    ! Find out which spin team we are in, and tag the spin we work on
+    call mpi_comm_rank( elsi_Spin_Comm, spin_rank, ierr )
+    my_spin = spin_rank+1  ! {1,2}
+
+
+    ! This is done serially, each time filling one spin set
+    ! Note that **all processes** need to have the same pkg_global
+
+    do ispin = 1, n_spin
+
+       ! Load pkg_global data package
+       pkg_global%norbs = n_basis
+       pkg_global%no_l  = n_basis_l
+       pkg_global%nnzl  = nnz_l
+       pkg_global%numcols => numh
+       pkg_global%cols    => col_idx
+
+       allocate(pkg_global%vals(2))
+       ! Link the vals items to the appropriate arrays (no extra memory here)
+       pkg_global%vals(1)%data => ovlp(:)
+       ! Note that we *cannot* say  => ham(:,my_spin)
+       ! and avoid the sequential loop, as then half the processors will send
+       ! the information for 'spin up' and the other half the information for 'spin down',
+       ! which is *not* what we want.
+       pkg_global%vals(2)%data => ham(:,ispin)
+
+       call timer("redist_orbs_fwd", 1)
+
+       ! We are doing the transfers sequentially. One spin team is
+       ! 'idle' (in the receiving side) in each pass, as the dist_spin(ispin) distribution
+       ! does not involve them.
+
+       call redistribute_spmatrix(n_basis,pkg_global,dist_global, &
+                                          pkg_spin,dist_spin(ispin),elsi_global_Comm)
+
+       call timer("redist_orbs_fwd", 2)
+
+       if (my_spin == ispin) then  ! Each team gets their own data
+
+          !nrows = pkg_spin%norbs          ! or simply 'norbs'
+          my_no_l = pkg_spin%no_l
+          my_nnz_l    = pkg_spin%nnzl
+          call MPI_AllReduce(my_nnz_l,my_nnz,1,MPI_integer,MPI_sum,elsi_Spatial_Comm,ierr)
+          ! generate off-by-one row pointer
+          call re_alloc(my_row_ptr2,1,my_no_l+1,"my_row_ptr2","elsi_solver")
+          my_row_ptr2(1) = 1
+          do ih = 1,my_no_l
+             my_row_ptr2(ih+1) = my_row_ptr2(ih) + pkg_spin%numcols(ih)
+          enddo
+
+          my_col_idx => pkg_spin%cols
+          my_S => pkg_spin%vals(1)%data
+          my_H => pkg_spin%vals(2)%data
+
+          call re_alloc(my_DM,1,my_nnz_l,"my_DM","elsi_solver")
+          call re_alloc(my_EDM,1,my_nnz_l,"my_EDM","elsi_solver")
+       endif
+
+       ! Clean pkg_global
+       nullify(pkg_global%vals(1)%data)
+       nullify(pkg_global%vals(2)%data)
+       deallocate(pkg_global%vals)
+       nullify(pkg_global%numcols)
+       nullify(pkg_global%cols)
+
+    enddo
+
+    call elsi_set_csc(elsi_h, my_nnz, my_nnz_l, my_no_l, my_col_idx, my_row_ptr2)
+    call de_alloc(my_row_ptr2,"my_row_ptr2","elsi_solver")
+
+    call elsi_set_csc_blk(elsi_h, BlockSize)
+    call elsi_set_spin(elsi_h, n_spin, my_spin)
+    call elsi_set_mpi(elsi_h, elsi_Spatial_comm)
+    call elsi_set_mpi_global(elsi_h, elsi_global_comm)
+
+ endif  ! n_spin
 
   call timer("elsi-solver", 1)
 
@@ -1217,7 +1216,6 @@
      col_idx, qtot, temp, ham, ovlp, dm, edm, ef, ets, &
      nkpnt, kpt_n, kpt, weight, kpt_comm, Get_EDM_Only)
 
-  use fdf,         only: fdf_get
   use m_mpi_utils, only: globalize_sum
   use parallel,    only: BlockSize
 #ifdef MPI
@@ -1350,116 +1348,116 @@
 
       !print *, global_rank, "| ", " Entering elsi_complex_solver"
 
-    if (n_spin == 1) then
-
-       ! Sparsity pattern
-       call globalize_sum(nnz_l, nnz_g, comm=kpt_comm)
-
-       allocate(row_ptr2(n_basis_l+1))
-       row_ptr2(1:n_basis_l) = row_ptr(1:n_basis_l)+1
-       row_ptr2(n_basis_l+1) = nnz_l+1
-
-       call elsi_set_csc(elsi_h, nnz_g, nnz_l, n_basis_l, col_idx, row_ptr2)
-       deallocate(row_ptr2)
-
-       call elsi_set_csc_blk(elsi_h, BlockSize)
-       call elsi_set_kpoint(elsi_h, nkpnt, kpt_n, weight)
-       call elsi_set_mpi(elsi_h, kpt_comm)
-       call elsi_set_mpi_global(elsi_h, elsi_global_comm)
-
-    else
-
-       call mpi_comm_rank( elsi_global_Comm, global_rank, ierr )
-
-       ! MPI logic for spin polarization
-
-       ! Split the communicator in spins and get distribution objects
-       ! for the data redistribution needed
-       ! Note that dist_spin is an array
-       call get_spin_comms_and_dists(kpt_comm,kpt_comm, &  !! **** kpt_comm as global?
-            blocksize, n_spin, &
-            dist_global,dist_spin, elsi_spatial_comm, elsi_spin_comm)
-
-       ! Find out which spin team we are in, and tag the spin we work on
-       call mpi_comm_rank( elsi_Spin_Comm, spin_rank, ierr )
-       my_spin = spin_rank+1  ! {1,2}
-
-      !print *, global_rank, "| ", "spin ", my_spin, " After spin splitting"
-
-       ! This is done serially, each time filling one spin set
-       ! Note that **all processes** need to have the same pkg_global
-
-       do ispin = 1, n_spin
-
-          ! Load pkg_global data package
-          pkg_global%norbs = n_basis
-          pkg_global%no_l  = n_basis_l
-          pkg_global%nnzl  = nnz_l
-          pkg_global%numcols => numh
-          pkg_global%cols    => col_idx
-
-          allocate(pkg_global%complex_vals(2))
-          ! Link the vals items to the appropriate arrays (no extra memory here)
-          pkg_global%complex_vals(1)%data => ovlp(:)
-          ! Note that we *cannot* say  => ham(:,my_spin)
-          ! and avoid the sequential loop, as then half the processors will send
-          ! the information for 'spin up' and the other half the information for 'spin down',
-          ! which is *not* what we want.
-          pkg_global%complex_vals(2)%data => ham(:,ispin)
-
-          call timer("redist_orbs_fwd", 1)
-
-          ! We are doing the transfers sequentially. One spin team is
-          ! 'idle' (in the receiving side) in each pass, as the dist_spin(ispin) distribution
-          ! does not involve them.
-
-          call redistribute_spmatrix(n_basis,pkg_global,dist_global, &
+ if (n_spin == 1) then
+
+    ! Sparsity pattern
+    call globalize_sum(nnz_l, nnz_g, comm=kpt_comm)
+
+    allocate(row_ptr2(n_basis_l+1))
+    row_ptr2(1:n_basis_l) = row_ptr(1:n_basis_l)+1
+    row_ptr2(n_basis_l+1) = nnz_l+1
+
+    call elsi_set_csc(elsi_h, nnz_g, nnz_l, n_basis_l, col_idx, row_ptr2)
+    deallocate(row_ptr2)
+
+    call elsi_set_csc_blk(elsi_h, BlockSize)
+    call elsi_set_kpoint(elsi_h, nkpnt, kpt_n, weight)
+    call elsi_set_mpi(elsi_h, kpt_comm)
+    call elsi_set_mpi_global(elsi_h, elsi_global_comm)
+
+ else
+
+    call mpi_comm_rank( elsi_global_Comm, global_rank, ierr )
+
+    ! MPI logic for spin polarization
+
+    ! Split the communicator in spins and get distribution objects
+    ! for the data redistribution needed
+    ! Note that dist_spin is an array
+    call get_spin_comms_and_dists(kpt_comm,kpt_comm, &  !! **** kpt_comm as global?
+         blocksize, n_spin, &
+         dist_global,dist_spin, elsi_spatial_comm, elsi_spin_comm)
+
+    ! Find out which spin team we are in, and tag the spin we work on
+    call mpi_comm_rank( elsi_Spin_Comm, spin_rank, ierr )
+    my_spin = spin_rank+1  ! {1,2}
+
+    !print *, global_rank, "| ", "spin ", my_spin, " After spin splitting"
+
+    ! This is done serially, each time filling one spin set
+    ! Note that **all processes** need to have the same pkg_global
+
+    do ispin = 1, n_spin
+
+       ! Load pkg_global data package
+       pkg_global%norbs = n_basis
+       pkg_global%no_l  = n_basis_l
+       pkg_global%nnzl  = nnz_l
+       pkg_global%numcols => numh
+       pkg_global%cols    => col_idx
+
+       allocate(pkg_global%complex_vals(2))
+       ! Link the vals items to the appropriate arrays (no extra memory here)
+       pkg_global%complex_vals(1)%data => ovlp(:)
+       ! Note that we *cannot* say  => ham(:,my_spin)
+       ! and avoid the sequential loop, as then half the processors will send
+       ! the information for 'spin up' and the other half the information for 'spin down',
+       ! which is *not* what we want.
+       pkg_global%complex_vals(2)%data => ham(:,ispin)
+
+       call timer("redist_orbs_fwd", 1)
+
+       ! We are doing the transfers sequentially. One spin team is
+       ! 'idle' (in the receiving side) in each pass, as the dist_spin(ispin) distribution
+       ! does not involve them.
+
+       call redistribute_spmatrix(n_basis,pkg_global,dist_global, &
                                              pkg_spin,dist_spin(ispin),kpt_Comm)
 
-          call timer("redist_orbs_fwd", 2)
-
-          if (my_spin == ispin) then  ! Each team gets their own data
-
-             !nrows = pkg_spin%norbs          ! or simply 'norbs'
-             my_no_l = pkg_spin%no_l
-             my_nnz_l    = pkg_spin%nnzl
-             call MPI_AllReduce(my_nnz_l,my_nnz,1,MPI_integer,MPI_sum,elsi_Spatial_Comm,ierr)
-             ! generate off-by-one row pointer
-             call re_alloc(my_row_ptr2,1,my_no_l+1,"my_row_ptr2","elsi_solver")
-             my_row_ptr2(1) = 1
-             do ih = 1,my_no_l
-                my_row_ptr2(ih+1) = my_row_ptr2(ih) + pkg_spin%numcols(ih)
-             enddo
-
-             my_col_idx => pkg_spin%cols
-             my_S => pkg_spin%complex_vals(1)%data
-             my_H => pkg_spin%complex_vals(2)%data
-
-             call re_alloc(my_DM,1,my_nnz_l,"my_DM","elsi_solver")
-             call re_alloc(my_EDM,1,my_nnz_l,"my_EDM","elsi_solver")
-          endif
-
-          ! Clean pkg_global
-          nullify(pkg_global%complex_vals(1)%data)
-          nullify(pkg_global%complex_vals(2)%data)
-          deallocate(pkg_global%complex_vals)
-          nullify(pkg_global%numcols)
-          nullify(pkg_global%cols)
-
-       enddo
-
-       !print *, global_rank, "| ", "spin ", my_spin, "Done spin transfers"
-
-       call elsi_set_csc(elsi_h, my_nnz, my_nnz_l, my_no_l, my_col_idx, my_row_ptr2)
-       call de_alloc(my_row_ptr2,"my_row_ptr2","elsi_solver")
-
-       call elsi_set_csc_blk(elsi_h, BlockSize)
-       call elsi_set_spin(elsi_h, n_spin, my_spin)
-       call elsi_set_kpoint(elsi_h, nkpnt, kpt_n, weight)
-       call elsi_set_mpi(elsi_h, elsi_Spatial_comm)
-       call elsi_set_mpi_global(elsi_h, elsi_global_comm)
-
-    endif  ! n_spin
+       call timer("redist_orbs_fwd", 2)
+
+       if (my_spin == ispin) then  ! Each team gets their own data
+
+          !nrows = pkg_spin%norbs          ! or simply 'norbs'
+          my_no_l = pkg_spin%no_l
+          my_nnz_l    = pkg_spin%nnzl
+          call MPI_AllReduce(my_nnz_l,my_nnz,1,MPI_integer,MPI_sum,elsi_Spatial_Comm,ierr)
+          ! generate off-by-one row pointer
+          call re_alloc(my_row_ptr2,1,my_no_l+1,"my_row_ptr2","elsi_solver")
+          my_row_ptr2(1) = 1
+          do ih = 1,my_no_l
+             my_row_ptr2(ih+1) = my_row_ptr2(ih) + pkg_spin%numcols(ih)
+          enddo
+
+          my_col_idx => pkg_spin%cols
+          my_S => pkg_spin%complex_vals(1)%data
+          my_H => pkg_spin%complex_vals(2)%data
+
+          call re_alloc(my_DM,1,my_nnz_l,"my_DM","elsi_solver")
+          call re_alloc(my_EDM,1,my_nnz_l,"my_EDM","elsi_solver")
+       endif
+
+       ! Clean pkg_global
+       nullify(pkg_global%complex_vals(1)%data)
+       nullify(pkg_global%complex_vals(2)%data)
+       deallocate(pkg_global%complex_vals)
+       nullify(pkg_global%numcols)
+       nullify(pkg_global%cols)
+
+    enddo
+
+    !print *, global_rank, "| ", "spin ", my_spin, "Done spin transfers"
+
+    call elsi_set_csc(elsi_h, my_nnz, my_nnz_l, my_no_l, my_col_idx, my_row_ptr2)
+    call de_alloc(my_row_ptr2,"my_row_ptr2","elsi_solver")
+
+    call elsi_set_csc_blk(elsi_h, BlockSize)
+    call elsi_set_spin(elsi_h, n_spin, my_spin)
+    call elsi_set_kpoint(elsi_h, nkpnt, kpt_n, weight)
+    call elsi_set_mpi(elsi_h, elsi_Spatial_comm)
+    call elsi_set_mpi_global(elsi_h, elsi_global_comm)
+
+ endif  ! n_spin
 
   call timer("elsi-solver", 1)