Merge pull request #25 from NLESC-JCER/opt

implementing deflation method

Author: felipeZ

CHANGELOG.md

@@ -1,4 +1,10 @@
 
+# Version 0.0.5
+
+### New
+
+ * Select the initial orthonormal basis set based on the lowest diagonal elements of the matrix
+
 # Version 0.0.4
 
 ## changed

src/array_utils.f90

@@ -2,13 +2,14 @@ module array_utils
 
   use numeric_kinds, only: dp
   use lapack_wrapper, only: lapack_generalized_eigensolver, lapack_matmul, lapack_matrix_vector, &
-       lapack_qr, lapack_solver
+       lapack_qr, lapack_solver, lapack_sort
   implicit none
 
   !> \private
   private
   !> \public
-  public :: concatenate, diagonal,eye, generate_diagonal_dominant, norm
+  public :: concatenate, diagonal,eye, generate_diagonal_dominant, norm, &
+       generate_preconditioner
 
 contains
 
@@ -33,7 +34,7 @@ pure function eye(m, n, alpha)
     do i=1, m
        do j=1, n
           if (i /= j) then
-             eye(i, j) = 0
+             eye(i, j) = 0.d0
           else
              eye(i, i) = x
           end if
@@ -132,4 +133,49 @@ function diagonal(matrix)
 
   end function diagonal  
 
+  function generate_preconditioner(diag, dim_sub) result(precond)
+    !> \brief generates a diagonal preconditioner for `matrix`.
+    !> \return diagonal matrix
+
+    ! input variable
+    real(dp), dimension(:), intent(inout) :: diag
+    integer, intent(in) :: dim_sub
+
+    ! local variables
+    real(dp), dimension(size(diag), dim_sub) :: precond
+    integer, dimension(size(diag)) :: keys
+    integer :: i, k
+    
+    ! sort diagonal
+    keys = lapack_sort('I', diag)
+    ! Fill matrix with zeros
+    precond = 0.0_dp
+
+    ! Add one depending on the order of the matrix diagonal
+    do i=1, dim_sub
+       k = search_key(keys, i)
+       precond(k, i) = 1.d0
+    end do
+    
+  end function generate_preconditioner
+
+  function search_key(keys, i) result(k)
+    !> \brief Search for a given index `i` in a vector `keys`
+    !> \param keys Vector of index
+    !> \param i Index to search for
+    !> \return index of i inside keys
+
+    integer, dimension(:), intent(in) :: keys
+    integer, intent(in) :: i
+    integer :: j, k
+    
+    do j=1,size(keys)
+       if (keys(j) == i) then
+          k = j
+          exit
+       end if
+    end do
+
+  end function search_key
+    
 end module array_utils

src/davidson.f90

@@ -1,12 +1,20 @@
 !> \namespace Davidson eigensolver
 !> \author Felipe Zapata
+!> The current implementation uses a general  davidson algorithm, meaning
+!> that it compute all the eigenvalues simultaneusly using a variable size block approach.
+!> The family of Davidson algorithm only differ in the way that the correction
+!> vector is computed.
+!> Computed pairs of eigenvalues/eigenvectors are deflated using algorithm
+!> described at: https://doi.org/10.1023/A:101919970
+
+
 module davidson_dense
   !> Submodule containing the implementation of the Davidson diagonalization method
   !> for dense matrices
   use numeric_kinds, only: dp
   use lapack_wrapper, only: lapack_generalized_eigensolver, lapack_matmul, lapack_matrix_vector, &
-       lapack_qr, lapack_solver
-  use array_utils, only: concatenate, eye, norm
+       lapack_qr, lapack_solver, lapack_sort
+  use array_utils, only: concatenate, diagonal, eye, generate_preconditioner, norm
 
   implicit none
 
@@ -41,11 +49,8 @@ end function compute_correction_generalized_dense
 
     subroutine generalized_eigensolver_dense(mtx, eigenvalues, ritz_vectors, lowest, method, max_iters, &
         tolerance, iters, max_dim_sub, stx)
-    !> The current implementation uses a general  davidson algorithm, meaning
-    !> that it compute all the eigenvalues simultaneusly using a block approach.
-    !> The family of Davidson algorithm only differ in the way that the correction
-    !> vector is computed.
-    
+     !> Implementation storing in memory the initial densed matrix mtx.
+      
     !> \param[in] mtx: Matrix to diagonalize
     !> \param[in, opt] Optional matrix to solve the general eigenvalue problem:
     !> \f$ mtx \lambda = V stx \lambda \f$
@@ -78,6 +83,7 @@ subroutine generalized_eigensolver_dense(mtx, eigenvalues, ritz_vectors, lowest,
 
     !local variables
     integer :: i, j, dim_sub, max_dim
+    integer :: n_converged ! Number of converged eigenvalue/eigenvector pairs
 
     ! Basis of subspace of approximants
     real(dp), dimension(size(mtx, 1)) :: guess, rs
@@ -87,12 +93,22 @@ subroutine generalized_eigensolver_dense(mtx, eigenvalues, ritz_vectors, lowest,
     real(dp), dimension(:), allocatable :: eigenvalues_sub
     real(dp), dimension(:, :), allocatable :: correction, eigenvectors_sub, mtx_proj, stx_proj, V
 
+    ! Diagonal matrix
+    real(dp), dimension(size(mtx, 1)) :: d
+    
     ! generalize problem
     logical :: gev 
 
+    ! indices of the eigenvalues/eigenvectors pair that have not converged
+    logical, dimension(lowest) :: has_converged
+
     ! Iteration subpsace dimension
     dim_sub = lowest * 2
 
+    ! Initial number of converged eigenvalue/eigenvector pairs
+    n_converged = 0
+    has_converged = .False.
+    
     ! maximum dimension of the basis for the subspace
     if (present(max_dim_sub)) then
        max_dim  = max_dim_sub
@@ -104,8 +120,10 @@ subroutine generalized_eigensolver_dense(mtx, eigenvalues, ritz_vectors, lowest,
     gev = present(stx)
 
     ! 1. Variables initialization
-    V = eye(size(ritz_vectors, 1), dim_sub) ! Initial orthonormal basis
-
+    ! Select the initial ortogonal subspace based on lowest elements
+    ! of the diagonal of the matrix
+    d = diagonal(mtx)
+    V = generate_preconditioner(d, dim_sub)
 
    ! 2. Generate subpace matrix problem by projecting into V
    mtx_proj = lapack_matmul('T', 'N', V, lapack_matmul('N', 'N', mtx, V))
@@ -145,9 +163,16 @@ subroutine generalized_eigensolver_dense(mtx, eigenvalues, ritz_vectors, lowest,
           end if
           rs = lapack_matrix_vector('N', mtx, ritz_vectors(:, j)) - guess
           errors(j) = norm(rs)
+          ! Check which eigenvalues has converged
+          if (errors(j) < tolerance) then
+             has_converged(j) = .true.
+          end if
        end do
-
-       if (all(errors < tolerance)) then
+       
+       ! Count converged pairs of eigenvalues/eigenvectors
+       n_converged = n_converged + count(errors < tolerance)
+       
+       if (all(has_converged)) then
           iters = i
           exit
        end if
@@ -197,7 +222,8 @@ subroutine generalized_eigensolver_dense(mtx, eigenvalues, ritz_vectors, lowest,
     end do outer_loop
 
     !  8. Check convergence
-    if (i > max_iters / dim_sub) then
+    if (i > max_iters) then
+       iters = i
        print *, "Warning: Algorithm did not converge!!"
     end if
 
@@ -265,7 +291,7 @@ module davidson_free
   use numeric_kinds, only: dp
   use lapack_wrapper, only: lapack_generalized_eigensolver, lapack_matmul, lapack_matrix_vector, &
        lapack_qr, lapack_solver
-  use array_utils, only: concatenate, eye, norm
+  use array_utils, only: concatenate, generate_preconditioner, norm
   use davidson_dense, only: generalized_eigensolver_dense
   implicit none
 
@@ -343,7 +369,7 @@ end function fun_stx_gemv
 
     ! ! Basis of subspace of approximants
     real(dp), dimension(size(ritz_vectors, 1),1) :: guess, rs
-    real(dp), dimension(size(ritz_vectors, 1)  ) :: diag_mtx, diag_stx
+    real(dp), dimension(size(ritz_vectors, 1)  ) :: diag_mtx, diag_stx, copy_d
     real(dp), dimension(lowest):: errors
 
     ! ! Working arrays
@@ -368,7 +394,10 @@ end function fun_stx_gemv
     diag_stx = extract_diagonal_free(fun_stx_gemv,dim_mtx)
 
     ! 1. Variables initialization
-    V = eye(dim_mtx, dim_sub) ! Initial orthonormal basis
+    ! Select the initial ortogonal subspace based on lowest elements
+    ! of the diagonal of the matrix
+    copy_d = diag_mtx
+    V = generate_preconditioner(copy_d, dim_sub) ! Initial orthonormal basis
 
     ! 2. Generate subspace matrix problem by projecting into V
     mtx_proj = lapack_matmul('T', 'N', V, fun_mtx_gemv(V))
@@ -757,7 +786,7 @@ function compute_DPR_generalized_dense(mtx, V, eigenvalues, eigenvectors, stx) r
     real(dp), dimension(:, :), intent(in) :: mtx, V, eigenvectors
     real(dp), dimension(:, :), intent(in), optional ::  stx 
     real(dp), dimension(size(mtx, 1), size(V, 2)) :: correction
-
+    
     ! local variables
     integer :: ii,j, m
     real(dp), dimension(size(mtx, 1), size(mtx, 2)) :: diag, arr
@@ -780,12 +809,12 @@ function compute_DPR_generalized_dense(mtx, V, eigenvalues, eigenvectors, stx) r
        correction(:, j) = lapack_matrix_vector('N', arr, vec) 
 
        do ii=1,size(correction,1)
-        if (gev) then
-          correction(ii, j) = correction(ii, j) / (eigenvalues(j) * stx(ii,ii) - mtx(ii, ii))
-        else 
-          correction(ii, j) = correction(ii, j) / (eigenvalues(j)  - mtx(ii, ii))
-        end if
-       end do
+          if (gev) then
+             correction(ii, j) = correction(ii, j) / (eigenvalues(j) * stx(ii,ii) - mtx(ii, ii))
+           else
+              correction(ii, j) = correction(ii, j) / (eigenvalues(j)  - mtx(ii, ii))
+           endif
+        end do
     end do
 
   end function compute_DPR_generalized_dense

src/lapack_wrapper.f90

@@ -7,7 +7,7 @@ module lapack_wrapper
   private
   !> \public
   public :: lapack_generalized_eigensolver, lapack_matmul, lapack_matrix_vector, &
-       lapack_qr, lapack_solver
+       lapack_qr, lapack_solver, lapack_sort
 
 contains
 
@@ -280,6 +280,35 @@ function lapack_matrix_vector(transA, mtx, vector, alpha) result(rs)
 
   end function lapack_matrix_vector
 
+
+  function lapack_sort(id, vector) result(keys)
+    !> \brief sort a vector of douboles
+    !> \param vector Array to sort
+    !> \param keys Indices of the sorted array
+    !> \param id Sorted either `I` increasing or `D` decreasing
+    real(dp), dimension(:), intent(inout) :: vector
+    character(len=1), intent(in) :: id
+    integer, dimension(size(vector)) :: keys
+
+    ! local variables
+    real(dp), dimension(size(vector)) :: xs
+    integer :: i, j, info
+    xs = vector
+    
+    call DLASRT(id, size(vector), vector, info)
+    call check_lapack_call(info, "DLASRT")
+    
+    do i=1,size(vector)
+       do j=1, size(vector)
+          if (abs(vector(j) - xs(i)) < 1e-16) then
+             keys(i) = j
+          end if
+       end do
+    end do
+    
+  end function lapack_sort
+  
+  
   subroutine check_lapack_call(info, name)
     !> Check if a subroutine finishes sucessfully
     !> \param info: Termination signal

src/main.f90

@@ -37,7 +37,7 @@ program main
 
   implicit none
 
-  integer, parameter :: dim = 50
+  integer, parameter :: dim = 100
   real(dp), dimension(3) :: eigenvalues_DPR, eigenvalues_GJD
   real(dp), dimension(dim, 3) :: eigenvectors_DPR, eigenvectors_GJD
   real(dp), dimension(dim, dim) :: mtx
@@ -46,27 +46,29 @@ program main
   real(dp), dimension(dim) :: xs
   integer :: iter_i, j
 
-  mtx = generate_diagonal_dominant(dim, 1d-4)
-  stx = generate_diagonal_dominant(dim, 1d-4, 1d0)
-
-  call generalized_eigensolver(mtx, eigenvalues_GJD, eigenvectors_GJD, 3, "GJD", 1000, 1d-8, iter_i, 10, stx)
-  call generalized_eigensolver(mtx, eigenvalues_DPR, eigenvectors_DPR, 3, "DPR", 1000, 1d-8, iter_i, 10, stx)
+  mtx = generate_diagonal_dominant(dim, 1d-3)
+  stx = generate_diagonal_dominant(dim, 1d-3, 1d0)
 
+  call generalized_eigensolver(mtx, eigenvalues_GJD, eigenvectors_GJD, 3, "GJD", 100, 1d-5, iter_i, 10, stx)
+  print *, "GJD algorithm converged in: ", iter_i, " iterations!"
+  call generalized_eigensolver(mtx, eigenvalues_DPR, eigenvectors_DPR, 3, "DPR", 100, 1d-5, iter_i, 10, stx)
+  print *, "DPR algorithm converged in: ", iter_i, " iterations!"
+  
   print *, "Test 1"
   test_norm_eigenvalues = norm(eigenvalues_GJD - eigenvalues_DPR)
   print *, "Check that eigenvalues norm computed by different methods are the same: ", test_norm_eigenvalues < 1e-6
 
   print *, "Test 2"
-  print *, "Check that eigenvalue equation:  H V = l S V "
+  print *, "Check that eigenvalue equation:  H V = l S V  holds!"
   print *, "DPR method:"
   do j=1,3
      xs = matmul(mtx, eigenvectors_DPR(:, j)) - (eigenvalues_DPR(j) * matmul(stx, eigenvectors_DPR(:, j)))
-     print *, "eigenvalue ", j, ": ", eigenvalues_DPR(j), " : ", norm(xs) , " : ", norm(xs)< 1.d-7
+     print *, "eigenvalue ", j, ": ", eigenvalues_DPR(j), "||Error||: ", norm(xs)
   end do
   print *, "GJD method:"
   do j=1,3
      xs = matmul(mtx, eigenvectors_GJD(:, j)) - (eigenvalues_GJD(j) * matmul( stx, eigenvectors_GJD(:, j)))
-     print *, "eigenvalue ", j, ": ",eigenvalues_GJD(j), " : ",  norm(xs), " : ", norm(xs)< 1.d-7
+     print *, "eigenvalue ", j, ": ",eigenvalues_GJD(j), "||Error||: ", norm(xs)
   end do  
 
 end program main

src/tests/CMakeLists.txt

@@ -17,7 +17,8 @@ foreach(PROG ${SOURCE_DAVIDSON})
 endforeach(PROG)
 # set(NUMPY_TEST_SRC ${CMAKE_SOURCE_DIR}/src/davidson.f90 ${CMAKE_SOURCE_DIR}/src/lapack_wrapper.f90 ${CMAKE_SOURCE_DIR}/src/numeric_kinds.f90)
 
-set(test_cases test_call_lapack test_dense_numpy test_free_numpy test_dense_properties test_free_properties)
+set(test_cases test_call_lapack test_dense_numpy test_free_numpy test_dense_properties
+  test_free_properties test_reorder)
 
 foreach(PROG ${test_cases})
   add_executable(${PROG}
@@ -50,6 +51,12 @@ add_test(
   WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/bin
   )
 
+add_test(
+  NAME test_reorder
+  COMMAND $<TARGET_FILE:test_reorder>
+  WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/bin  
+  )
+
 # Tests requiring python, numpy and scipy
 if(${_numpy_status} EQUAL 0)
   add_test(