From a6ef4367e715554f0561c4e2a1f47568b7edc745 Mon Sep 17 00:00:00 2001 From: felipe Date: Fri, 22 Feb 2019 14:17:27 +0100 Subject: [PATCH] split implementation into smaller modules #19 --- src/CMakeLists.txt | 5 +- src/davidson.f90 | 1130 +++++++++++---------------- src/lapack_wrapper.f90 | 298 +++++++ src/main.f90 | 3 +- src/tests/CMakeLists.txt | 2 +- src/tests/test_call_lapack.f90 | 3 +- src/tests/test_dense_properties.f90 | 2 +- src/tests/test_free_properties.f90 | 2 +- 8 files changed, 746 insertions(+), 699 deletions(-) create mode 100644 src/lapack_wrapper.f90 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 83642e3..f05eba6 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,6 +1,7 @@ # create binary -set(SOURCES main.f90 numeric_kinds.f90 davidson.f90 benchmark.f90) +set(UTILS numeric_kinds.f90 lapack_wrapper.f90) +set(SOURCES main.f90 davidson.f90 benchmark.f90 ${UTILS}) message (STATUS "SOURCES: " ${SOURCES}) add_executable(main ${SOURCES}) @@ -29,7 +30,7 @@ target_compile_options(main target_link_libraries(main PRIVATE ${LINEAR_ALGEBRA} ${OpenMP_Fortran_FLAGS}) # Benchmark free matrix -add_executable(benchmark_free benchmark_free.f90 davidson.f90 numeric_kinds.f90) +add_executable(benchmark_free benchmark_free.f90 davidson.f90 ${UTILS}) target_compile_options(main PRIVATE diff --git a/src/davidson.f90 b/src/davidson.f90 index 0ad5420..ed1d436 100644 --- a/src/davidson.f90 +++ b/src/davidson.f90 @@ -3,36 +3,15 @@ module davidson use numeric_kinds, only: dp + use lapack_wrapper, only: lapack_generalized_eigensolver, lapack_matmul, lapack_matrix_vector, & + lapack_qr, lapack_solver implicit none !> \private private !> \public - public :: eye, generalized_eigensolver, generate_diagonal_dominant, norm, lapack_generalized_eigensolver, lapack_qr + public :: eye, generalized_eigensolver, generate_diagonal_dominant, norm - interface - - module function compute_correction_generalized_dense(mtx, V, eigenvalues, eigenvectors, method, stx) & - result(correction) - !> compute the correction vector using a given `method` for the Davidson algorithm - !> See correction_methods submodule for the implementations - !> \param[in] mtx: Original matrix - !> \param[in] stx: Matrix to compute the general eigenvalue problem - !> \param[in] V: Basis of the iteration subspace - !> \param[in] eigenvalues: of the reduce problem - !> \param[in] eigenvectors: of the reduce problem - !> \param[in] method: name of the method to compute the correction - - real(dp), dimension(:), intent(in) :: eigenvalues - real(dp), dimension(:, :), intent(in) :: mtx, V, eigenvectors - real(dp), dimension(:, :), intent(in), optional :: stx - character(len=*), optional, intent(in) :: method - real(dp), dimension(size(mtx, 1), size(V, 2)) :: correction - - end function compute_correction_generalized_dense - - end interface - interface generalized_eigensolver !> \brief Solve a (general) eigenvalue problem using different types of Davidson algorithms. @@ -56,191 +35,73 @@ end function compute_correction_generalized_dense module procedure generalized_eigensolver_dense module procedure generalized_eigensolver_free - - end interface generalized_eigensolver -contains + end interface generalized_eigensolver - 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. - - !> \param[in] mtx: Matrix to diagonalize - !> \param[in, opt] Optional matrix to solve the general eigenvalue problem: - !> \f$ mtx \lambda = V stx \lambda \f$ - !> \param[out] eigenvalues Computed eigenvalues - !> \param[out] ritz_vectors approximation to the eigenvectors - !> \param[in] lowest Number of lowest eigenvalues/ritz_vectors to compute - !> \param[in] method Method to compute the correction vector. Available - !> methods are, - !> DPR: Diagonal-Preconditioned-Residue - !> GJD: Generalized Jacobi Davidson - !> \param[in] max_iters: Maximum number of iterations - !> \param[in] tolerance norm-2 error of the eigenvalues - !> \param[in] method: Method to compute the correction vectors - !> \param[in, opt] max_dim_sub: maximum dimension of the subspace search - !> \param[out] iters: Number of iterations until convergence - !> \return eigenvalues and ritz_vectors of the matrix `mtx` - - implicit none - ! input/output variable - integer, intent(in) :: lowest - real(dp), dimension(:, :), intent(in) :: mtx - real(dp), dimension(:, :), intent(in), optional :: stx - real(dp), dimension(lowest), intent(out) :: eigenvalues - real(dp), dimension(:, :), intent(out) :: ritz_vectors - integer, intent(in) :: max_iters - integer, intent(in), optional :: max_dim_sub - real(dp), intent(in) :: tolerance - character(len=*), intent(in) :: method - integer, intent(out) :: iters - - !local variables - integer :: i, j, dim_sub, max_dim - - ! Basis of subspace of approximants - real(dp), dimension(size(mtx, 1)) :: guess, rs - real(dp), dimension(lowest):: errors - - ! Working arrays - real(dp), dimension(:), allocatable :: eigenvalues_sub - real(dp), dimension(:, :), allocatable :: correction, eigenvectors_sub, mtx_proj, stx_proj, V - - ! generalize problem - logical :: gev - - ! Iteration subpsace dimension - dim_sub = lowest * 2 - - ! maximum dimension of the basis for the subspace - if (present(max_dim_sub)) then - max_dim = max_dim_sub - else - max_dim = lowest * 10 - endif - - ! generalied problem - gev = present(stx) - - ! 1. Variables initialization - V = eye(size(ritz_vectors, 1), dim_sub) ! Initial orthonormal basis - - - ! 2. Generate subpace matrix problem by projecting into V - mtx_proj = lapack_matmul('T', 'N', V, lapack_matmul('N', 'N', mtx, V)) - - if(gev) then - stx_proj = lapack_matmul('T', 'N', V, lapack_matmul('N', 'N', stx, V)) - end if - - ! ! Outer loop block Davidson schema - outer_loop: do i=1, max_iters - - ! 3. compute the eigenvalues and their corresponding ritz_vectors - ! for the projected matrix using lapack - call check_deallocate_matrix(eigenvectors_sub) - - if (allocated(eigenvalues_sub)) then - deallocate(eigenvalues_sub) - end if - - allocate(eigenvalues_sub(size(mtx_proj, 1))) - allocate(eigenvectors_sub(size(mtx_proj, 1), size(mtx_proj, 2))) - - - if (gev) then - call lapack_generalized_eigensolver(mtx_proj, eigenvalues_sub, eigenvectors_sub, stx_proj) - else - call lapack_generalized_eigensolver(mtx_proj, eigenvalues_sub, eigenvectors_sub) - end if - - ! 4. Check for convergence - ritz_vectors = lapack_matmul('N', 'N', V, eigenvectors_sub(:, :lowest)) - do j=1,lowest - if(gev) then - guess = eigenvalues_sub(j) * lapack_matrix_vector('N',stx,ritz_vectors(:, j)) - else - guess = eigenvalues_sub(j) * ritz_vectors(:, j) - end if - rs = lapack_matrix_vector('N', mtx, ritz_vectors(:, j)) - guess - errors(j) = norm(rs) - end do - - if (all(errors < tolerance)) then - iters = i - exit - end if + interface + + module 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. - ! 5. Add the correction vectors to the current basis - if (size(V, 2) <= max_dim) then - - ! append correction to the current basis - call check_deallocate_matrix(correction) - allocate(correction(size(mtx, 1), size(V, 2))) - - if(gev) then - correction = compute_correction_generalized_dense(mtx, V, eigenvalues_sub, eigenvectors_sub, method, stx) - else - correction = compute_correction_generalized_dense(mtx, V, eigenvalues_sub, eigenvectors_sub, method) - end if - - - ! 6. Increase Basis size - call concatenate(V, correction) + !> \param[in] mtx: Matrix to diagonalize + !> \param[in, opt] Optional matrix to solve the general eigenvalue problem: + !> \f$ mtx \lambda = V stx \lambda \f$ + !> \param[out] eigenvalues Computed eigenvalues + !> \param[out] ritz_vectors approximation to the eigenvectors + !> \param[in] lowest Number of lowest eigenvalues/ritz_vectors to compute + !> \param[in] method Method to compute the correction vector. Available + !> methods are, + !> DPR: Diagonal-Preconditioned-Residue + !> GJD: Generalized Jacobi Davidson + !> \param[in] max_iters: Maximum number of iterations + !> \param[in] tolerance norm-2 error of the eigenvalues + !> \param[in] method: Method to compute the correction vectors + !> \param[in, opt] max_dim_sub: maximum dimension of the subspace search + !> \param[out] iters: Number of iterations until convergence + !> \return eigenvalues and ritz_vectors of the matrix `mtx` + use numeric_kinds, only: dp + implicit none + ! input/output variable + integer, intent(in) :: lowest + real(dp), dimension(:, :), intent(in) :: mtx + real(dp), dimension(:, :), intent(in), optional :: stx + real(dp), dimension(lowest), intent(out) :: eigenvalues + real(dp), dimension(:, :), intent(out) :: ritz_vectors + integer, intent(in) :: max_iters + integer, intent(in), optional :: max_dim_sub + real(dp), intent(in) :: tolerance + character(len=*), intent(in) :: method + integer, intent(out) :: iters - ! 7. Orthogonalize basis - call lapack_qr(V) - - - ! 8. Update the the projection - call update_projection_dense(mtx, V, mtx_proj) - if (gev) then - call update_projection_dense(stx, V, stx_proj) - end if - - else - - ! 6. Otherwise reduce the basis of the subspace to the current correction - V = lapack_matmul('N', 'N', V, eigenvectors_sub(:, :dim_sub)) - - ! we refresh the projected matrices - mtx_proj = lapack_matmul('T', 'N', V, lapack_matmul('N', 'N', mtx, V)) - - if(gev) then - stx_proj = lapack_matmul('T', 'N', V, lapack_matmul('N', 'N', stx, V)) - end if - - end if - - - end do outer_loop - - ! 8. Check convergence - if (i > max_iters / dim_sub) then - print *, "Warning: Algorithm did not converge!!" - end if - - - ! Select the lowest eigenvalues and their corresponding ritz_vectors - ! They are sort in increasing order - eigenvalues = eigenvalues_sub(:lowest) - - ! Free memory - call check_deallocate_matrix(correction) - deallocate(eigenvalues_sub, eigenvectors_sub, V, mtx_proj) - - ! free optional matrix - if (gev) then - call check_deallocate_matrix(stx_proj) - endif - - end subroutine generalized_eigensolver_dense - + end subroutine generalized_eigensolver_dense + module function compute_correction_generalized_dense(mtx, V, eigenvalues, eigenvectors, method, stx) & + result(correction) + !> compute the correction vector using a given `method` for the Davidson algorithm + !> See correction_methods submodule for the implementations + !> \param[in] mtx: Original matrix + !> \param[in] stx: Matrix to compute the general eigenvalue problem + !> \param[in] V: Basis of the iteration subspace + !> \param[in] eigenvalues: of the reduce problem + !> \param[in] eigenvectors: of the reduce problem + !> \param[in] method: name of the method to compute the correction + + real(dp), dimension(:), intent(in) :: eigenvalues + real(dp), dimension(:, :), intent(in) :: mtx, V, eigenvectors + real(dp), dimension(:, :), intent(in), optional :: stx + character(len=*), optional, intent(in) :: method + real(dp), dimension(size(mtx, 1), size(V, 2)) :: correction + + end function compute_correction_generalized_dense + + end interface +contains + subroutine generalized_eigensolver_free(fun_mtx, eigenvalues, ritz_vectors, lowest, method, max_iters, & tolerance, iters, max_dim_sub, fun_stx) !> \brief use a pair of functions fun_mtx and fun_stx to compute on the fly the matrices to solve @@ -249,7 +110,7 @@ subroutine generalized_eigensolver_free(fun_mtx, eigenvalues, ritz_vectors, lowe !> 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. - + !> \param[in] fun_mtx: Function to compute the Matrix to diagonalize !> \param[in, opt] fun_stx: function to compute the optional general eigenvalue problem. !> \param[out] eigenvalues Computed eigenvalues @@ -275,7 +136,7 @@ subroutine generalized_eigensolver_free(fun_mtx, eigenvalues, ritz_vectors, lowe real(dp), intent(in) :: tolerance character(len=*), intent(in) :: method integer, intent(out) :: iters - + ! Function to compute the target matrix on the fly interface function fun_mtx(i, dim) result(vec) @@ -288,7 +149,7 @@ function fun_mtx(i, dim) result(vec) real(dp), dimension(dim) :: vec end function fun_mtx - + function fun_stx(i, dim) result(vec) !> \brief Fucntion to compute the optional stx matrix on the fly !> \param[in] i column/row to compute from stx @@ -297,7 +158,7 @@ function fun_stx(i, dim) result(vec) integer, intent(in) :: i integer, intent(in) :: dim real(dp), dimension(dim) :: vec - + end function fun_stx end interface @@ -307,31 +168,31 @@ end function fun_stx ! ! Basis of subspace of approximants real(dp), dimension(size(ritz_vectors, 1)) :: guess, rs real(dp), dimension(lowest):: errors - + ! ! Working arrays real(dp), dimension(:), allocatable :: eigenvalues_sub real(dp), dimension(:, :), allocatable :: correction, eigenvectors_sub, mtx_proj, stx_proj, V - + ! Iteration subpsace dimension dim_sub = lowest * 2 - + ! maximum dimension of the basis for the subspace if (present(max_dim_sub)) then max_dim = max_dim_sub else max_dim = lowest * 10 endif - + ! dimension of the matrix dim_mtx = size(ritz_vectors, 1) ! 1. Variables initialization V = eye(dim_mtx, dim_sub) ! Initial orthonormal basis - - ! 2. Generate subspace matrix problem by projecting into V + + ! 2. Generate subspace matrix problem by projecting into V mtx_proj = lapack_matmul('T', 'N', V, free_matmul(fun_mtx, V)) stx_proj = lapack_matmul('T', 'N', V, free_matmul(fun_stx, V)) - + ! Outer loop block Davidson schema outer_loop: do i=1, max_iters ! 3. compute the eigenvalues and their corresponding ritz_vectors @@ -341,12 +202,12 @@ end function fun_stx if (allocated(eigenvalues_sub)) then deallocate(eigenvalues_sub) end if - + allocate(eigenvalues_sub(size(mtx_proj, 1))) allocate(eigenvectors_sub(size(mtx_proj, 1), size(mtx_proj, 2))) - + call lapack_generalized_eigensolver(mtx_proj, eigenvalues_sub, eigenvectors_sub, stx_proj) - + ! 4. Check for convergence ritz_vectors = lapack_matmul('N', 'N', V, eigenvectors_sub(:, :lowest)) do j=1,lowest @@ -354,7 +215,7 @@ end function fun_stx rs = free_matrix_vector(fun_mtx, ritz_vectors(:, j), dim_mtx) - guess errors(j) = norm(rs) end do - + if (all(errors < tolerance)) then iters = i exit @@ -362,42 +223,42 @@ end function fun_stx ! 5. Add the correction vectors to the current basis if (size(V, 2) <= max_dim) then - + ! append correction to the current basis call check_deallocate_matrix(correction) allocate(correction(size(ritz_vectors, 1), size(V, 2))) - + correction = compute_DPR_free(fun_mtx, fun_stx, V, eigenvalues_sub, eigenvectors_sub) - + ! 6. Increase Basis size call concatenate(V, correction) - + ! 7. Orthogonalize basis call lapack_qr(V) - + ! 8. Update the the projection call update_projection_free(fun_mtx, V, mtx_proj) call update_projection_free(fun_stx, V, stx_proj) - + else - + ! 6. Otherwise reduce the basis of the subspace to the current correction V = lapack_matmul('N', 'N', V, eigenvectors_sub(:, :dim_sub)) - + ! we refresh the projected matrices mtx_proj = lapack_matmul('T', 'N', V, free_matmul(fun_mtx, V)) stx_proj = lapack_matmul('T', 'N', V, free_matmul(fun_stx, V)) - + end if - + end do outer_loop - + ! 8. Check convergence if (i > max_iters / dim_sub) then print *, "Warning: Algorithm did not converge!!" end if - + ! Select the lowest eigenvalues and their corresponding ritz_vectors ! They are sort in increasing order @@ -406,157 +267,124 @@ end function fun_stx ! Free memory call check_deallocate_matrix(correction) deallocate(eigenvalues_sub, eigenvectors_sub, V, mtx_proj) - + ! free optional matrix call check_deallocate_matrix(stx_proj) end subroutine generalized_eigensolver_free - - function compute_DPR_free(fun_mtx, fun_stx, V, eigenvalues, eigenvectors) result(correction) - !> compute the correction vector using the DPR method for a matrix free diagonalization - !> See correction_methods submodule for the implementations - !> \param[in] fun_mtx: function to compute matrix - !> \param[in] fun_stx: function to compute the matrix for the generalized case - !> \param[in] V: Basis of the iteration subspace - !> \param[in] eigenvalues: of the reduce problem - !> \param[in] eigenvectors: of the reduce problem - !> \return correction matrix - - real(dp), dimension(:), intent(in) :: eigenvalues - real(dp), dimension(:, :), intent(in) :: V, eigenvectors - real(dp), dimension(size(V, 1), size(V, 2)) :: correction - - interface - function fun_mtx(i, dim) result(vec) - !> \brief Function to compute the optional mtx on the fly - !> \param[in] i column/row to compute from mtx - !> \param vec column/row from mtx - use numeric_kinds, only: dp - integer, intent(in) :: i - integer, intent(in) :: dim - real(dp), dimension(dim) :: vec - - end function fun_mtx - - function fun_stx(i, dim) result(vec) - !> \brief Fucntion to compute the optional stx matrix on the fly - !> \param[in] i column/row to compute from stx - !> \param vec column/row from stx - use numeric_kinds, only: dp - integer, intent(in) :: i - integer, intent(in) :: dim - real(dp), dimension(dim) :: vec - - end function fun_stx - - end interface - - ! local variables - integer :: ii, j - real(dp), dimension(size(V, 1)) :: diag_mtx, diag_stx, rs, vector - - do j=1, size(V, 2) - vector = lapack_matrix_vector('N', V, eigenvectors(:, j)) - call compute_error(fun_mtx, fun_stx, eigenvalues(j), vector, rs, diag_mtx, diag_stx) - correction(:, j) = rs - do ii=1,size(correction,1) - correction(ii, j) = correction(ii, j) / (eigenvalues(j) * diag_stx(ii) - diag_mtx(ii)) - end do - end do - - end function compute_DPR_free - - subroutine compute_error(fun_mtx, fun_stx, eigenval, vector, rs, diag_mtx, diag_stx) - !> \brief compute the multiplication: (mtx - eigenval) * vector and - !> return the diagonal of the matrix - !> \param[in] fun_mtx: function to compute the matrix - !> \param[in] fun_stx: function to compute the matrix for the generalized case - !> \param[in] eigenval - !> \param[in] vector - !> \param[out] rs correction vector - !> \param[out] diagonal elements of the matrix - !> \param[out] diagonal elements of the matrix to compute the generalized problem - real(dp), dimension(:), intent(in) :: vector - real(dp), dimension(:), intent(out) :: diag_mtx, diag_stx, rs - real(dp) :: eigenval - - ! local variables - integer :: i - real(dp), dimension(size(vector)) ::xs, ys - - interface - function fun_mtx(i, dim) result(vec) - !> \brief Function to compute the optional mtx on the fly - !> \param[in] i column/row to compute from mtx - !> \param[in] j optional second index - !> \param vec column/row from mtx - use numeric_kinds, only: dp - integer, intent(in) :: i - integer, intent(in) :: dim - real(dp), dimension(dim) :: vec + + function compute_DPR_free(fun_mtx, fun_stx, V, eigenvalues, eigenvectors) result(correction) + !> compute the correction vector using the DPR method for a matrix free diagonalization + !> See correction_methods submodule for the implementations + !> \param[in] fun_mtx: function to compute matrix + !> \param[in] fun_stx: function to compute the matrix for the generalized case + !> \param[in] V: Basis of the iteration subspace + !> \param[in] eigenvalues: of the reduce problem + !> \param[in] eigenvectors: of the reduce problem + !> \return correction matrix + + real(dp), dimension(:), intent(in) :: eigenvalues + real(dp), dimension(:, :), intent(in) :: V, eigenvectors + real(dp), dimension(size(V, 1), size(V, 2)) :: correction + + interface + function fun_mtx(i, dim) result(vec) + !> \brief Function to compute the optional mtx on the fly + !> \param[in] i column/row to compute from mtx + !> \param vec column/row from mtx + use numeric_kinds, only: dp + integer, intent(in) :: i + integer, intent(in) :: dim + real(dp), dimension(dim) :: vec + + end function fun_mtx - end function fun_mtx - - function fun_stx(i, dim) result(vec) - !> \brief Fucntion to compute the optional stx matrix on the fly - !> \param[in] i column/row to compute from stx - !> \param vec column/row from stx - use numeric_kinds, only: dp - integer, intent(in) :: i - integer, intent(in) :: dim - real(dp), dimension(dim) :: vec - - end function fun_stx - - end interface - - !$OMP PARALLEL DO & - !$OMP PRIVATE(i, xs, ys) - do i = 1, size(vector) - xs = fun_mtx(i, size(vector)) - ys = fun_stx(i, size(vector)) - diag_mtx(i) = xs(i) - diag_stx(i) = ys(i) - xs = xs - eigenval * ys - rs(i) = dot_product(xs, vector) - end do - !$OMP END PARALLEL DO - - end subroutine compute_error - + function fun_stx(i, dim) result(vec) + !> \brief Fucntion to compute the optional stx matrix on the fly + !> \param[in] i column/row to compute from stx + !> \param vec column/row from stx + use numeric_kinds, only: dp + integer, intent(in) :: i + integer, intent(in) :: dim + real(dp), dimension(dim) :: vec + + end function fun_stx + + end interface + + ! local variables + integer :: ii, j + real(dp), dimension(size(V, 1)) :: diag_mtx, diag_stx, rs, vector + + do j=1, size(V, 2) + vector = lapack_matrix_vector('N', V, eigenvectors(:, j)) + call compute_error(fun_mtx, fun_stx, eigenvalues(j), vector, rs, diag_mtx, diag_stx) + correction(:, j) = rs + do ii=1,size(correction,1) + correction(ii, j) = correction(ii, j) / (eigenvalues(j) * diag_stx(ii) - diag_mtx(ii)) + end do + end do + + end function compute_DPR_free + + subroutine compute_error(fun_mtx, fun_stx, eigenval, vector, rs, diag_mtx, diag_stx) + !> \brief compute the multiplication: (mtx - eigenval) * vector and + !> return the diagonal of the matrix + !> \param[in] fun_mtx: function to compute the matrix + !> \param[in] fun_stx: function to compute the matrix for the generalized case + !> \param[in] eigenval + !> \param[in] vector + !> \param[out] rs correction vector + !> \param[out] diagonal elements of the matrix + !> \param[out] diagonal elements of the matrix to compute the generalized problem + real(dp), dimension(:), intent(in) :: vector + real(dp), dimension(:), intent(out) :: diag_mtx, diag_stx, rs + real(dp) :: eigenval + + ! local variables + integer :: i + real(dp), dimension(size(vector)) ::xs, ys + + interface + function fun_mtx(i, dim) result(vec) + !> \brief Function to compute the optional mtx on the fly + !> \param[in] i column/row to compute from mtx + !> \param[in] j optional second index + !> \param vec column/row from mtx + use numeric_kinds, only: dp + integer, intent(in) :: i + integer, intent(in) :: dim + real(dp), dimension(dim) :: vec + + end function fun_mtx + + function fun_stx(i, dim) result(vec) + !> \brief Fucntion to compute the optional stx matrix on the fly + !> \param[in] i column/row to compute from stx + !> \param vec column/row from stx + use numeric_kinds, only: dp + integer, intent(in) :: i + integer, intent(in) :: dim + real(dp), dimension(dim) :: vec + + end function fun_stx + + end interface + + !$OMP PARALLEL DO & + !$OMP PRIVATE(i, xs, ys) + do i = 1, size(vector) + xs = fun_mtx(i, size(vector)) + ys = fun_stx(i, size(vector)) + diag_mtx(i) = xs(i) + diag_stx(i) = ys(i) + xs = xs - eigenval * ys + rs(i) = dot_product(xs, vector) + end do + !$OMP END PARALLEL DO + + end subroutine compute_error - subroutine update_projection_dense(A, V, A_proj) - !> \brief update the projected matrices - !> \param A: full matrix - !> \param V: projector - !> \param A_proj: projected matrix - - implicit none - real(dp), dimension(:, :), intent(in) :: A - real(dp), dimension(:, :), intent(in) :: V - real(dp), dimension(:, :), intent(inout), allocatable :: A_proj - real(dp), dimension(:, :), allocatable :: tmp_array - - ! local variables - integer :: nvec, old_dim - - ! dimension of the matrices - nvec = size(V,2) - old_dim = size(A_proj,1) - - ! move to temporal array - allocate(tmp_array(nvec, nvec)) - tmp_array(:old_dim, :old_dim) = A_proj - tmp_array(:,old_dim+1:) = lapack_matmul('T', 'N', V, lapack_matmul('N', 'N', A, V(:, old_dim+1:))) - tmp_array( old_dim+1:,:old_dim ) = transpose(tmp_array(:old_dim, old_dim+1:)) - - ! Move to new expanded matrix - deallocate(A_proj) - call move_alloc(tmp_array, A_proj) - - - end subroutine update_projection_dense - subroutine update_projection_free(fun_A, V, A_proj) !> \brief update the projected matrices !> \param A: full matrix @@ -579,262 +407,26 @@ function fun_A(i, dim) result(vec) real(dp), dimension(dim) :: vec end function fun_A - end interface - - ! local variables - integer :: nvec, old_dim - - ! dimension of the matrices - nvec = size(V,2) - old_dim = size(A_proj,1) - - ! move to temporal array - allocate(tmp_array(nvec, nvec)) - tmp_array(:old_dim, :old_dim) = A_proj - tmp_array(:,old_dim+1:) = lapack_matmul('T', 'N', V, free_matmul(fun_A, V(:, old_dim+1:))) - tmp_array( old_dim+1:,:old_dim ) = transpose(tmp_array(:old_dim, old_dim+1:)) - - ! Move to new expanded matrix - deallocate(A_proj) - call move_alloc(tmp_array, A_proj) - - end subroutine update_projection_free - - - subroutine lapack_generalized_eigensolver(mtx, eigenvalues, eigenvectors, stx) - !> Call the DSYGV subroutine lapack to compute ALL the eigenvalues - !> and corresponding eigenvectors of mtx - !> \param mtx: Matrix to diaogonalize - !> \param stx: Overlap Matrix to diaogonalize - !> \param eigenvalues: lowest eigenvalues - !> \param eigenvectors: corresponding eigenvectors - !> \return eigenvalues/eigenvectors - - ! input/output - implicit none - real(dp), dimension(:, :), intent(in) :: mtx - real(dp), dimension(:, :), intent(in), optional :: stx - real(dp), dimension(size(mtx, 1)), intent(inout) :: eigenvalues - real(dp), dimension(size(mtx, 1), size(mtx, 2)), intent(inout) :: eigenvectors - - real(dp), dimension(:, :), allocatable :: mtx_copy - real(dp), dimension(:, :), allocatable :: stx_copy - - - ! Local variables - integer :: dim, info, lwork, itype = 1 - logical :: gev - - ! ALL the eigenvalues of the subpace (re, im) - real(dp), dimension(size(mtx, 1)) :: eigenvalues_work - real(dp), dimension(:), allocatable :: work ! workspace, see lapack documentation - - ! ! dimension of the guess space - dim = size(mtx, 1) - gev = present(stx) - - ! local copy of the matrices - allocate(mtx_copy(dim,dim)) - mtx_copy = mtx - - if (gev) then - allocate(stx_copy(dim,dim)) - stx_copy = stx - end if - - ! Query size of the optimal workspace - allocate(work(1)) - - if (gev) then - call DSYGV(itype,"V", "U", dim, mtx_copy, dim, stx_copy, dim, eigenvalues_work, work, -1, info) - call check_lapack_call(info, "DSYGV") - else - call DSYEV("V", "U", dim, mtx_copy, dim, eigenvalues_work, work, -1, info) - call check_lapack_call(info, "DSYEV") - end if - - ! Allocate memory for the workspace - lwork = max(1, int(work(1))) - deallocate(work) - allocate(work(lwork)) - - ! Compute Eigenvalues - if (gev) then - call DSYGV(itype,"V", "U", dim, mtx_copy, dim, stx_copy, dim, eigenvalues_work, work, lwork, info) - call check_lapack_call(info, "DSYGV") - else - call DSYEV("V", "U", dim, mtx_copy, dim, eigenvalues_work, work, lwork, info) - call check_lapack_call(info, "DSYEV") - end if - - ! Sort the eigenvalues and eigenvectors of the basis - eigenvalues = eigenvalues_work - eigenvectors = mtx_copy - - ! release memory - deallocate(work) - deallocate(mtx_copy) - if (gev) then - deallocate(stx_copy) - end if - - end subroutine lapack_generalized_eigensolver - - - subroutine lapack_qr(basis) - !> Orthoghonalize the basis using the QR factorization. - !> QR factorization of the M-by-N (M>N) matrx A=Q*R in the form where - !> Q is square M-by-M matrix and R is an upper triangular M-by-N matrix. - !> The equality A=Q*R can be re-written also as a product Q1*R1 where Q1 - !> is a rectangular M-by-N submatrix of the matrix Q and R1 is M-by-M - !> submatrix of the R. Let us note that columns of Q1 are orthonormal - !> (they are orthogonal to each other and have norms equal to 1). - !> The equality A=Q1*R1 can be treated as every column of A is a linear - !> combination of Q1 columns, i.e. they span the same linear space. - !> In other words, columns of Q1 is the result of ortogonalization of columns A. - !> DGEQRF does not not compute Q directly, DORGQR must be call subsequently. - - !> \param basis - !> \return orthogonal basis - - implicit none - real(dp), dimension(:, :), intent(inout) :: basis - real(dp), dimension(:), allocatable :: work ! workspace, see lapack documentation - real(dp), dimension(size(basis, 2)) :: tau ! see DGEQRF documentation - integer :: info, lwork, m, n - - ! Matrix shape - m = size(basis, 1) - n = size(basis, 2) - - ! 1. Call the QR decomposition - ! 1.1 Query size of the workspace (Check lapack documentation) - allocate(work(1)) - call DGEQRF(m, n, basis, m, tau, work, -1, info) - call check_lapack_call(info, "DGEQRF") - - ! 1.2 Allocate memory for the workspace - lwork = max(1, int(work(1))) - deallocate(work) - allocate(work(lwork)) - - ! 1.3 Call QR factorization - call DGEQRF(m, n, basis, m, tau, work, lwork, info) - call check_lapack_call(info, "DGEQRF") - deallocate(work) - - ! 2. Generates an orthonormal matrix - ! 2.1 Query size of the workspace (Check lapack documentation) - allocate(work(1)) - call DORGQR(m, n, min(m, n), basis, m, tau, work, -1, info) - call check_lapack_call(info, "DORGQR") - - ! 2.2 Allocate memory fo the workspace - lwork = max(1, int(work(1))) - deallocate(work) - allocate(work(lwork)) - - ! 2.3 compute the matrix Q - call DORGQR(m, n, min(m, n), basis, m, tau, work, lwork, info) - call check_lapack_call(info, "DORGQR") - - ! release memory - deallocate(work) - - end subroutine lapack_qr - - subroutine lapack_solver(arr, brr) - !> Call lapack DSYSV subroutine to solve a AX=B Linear system - !> \param arr: matrix with the coefficients of the linear system - !> \param brr: Vector with the constant terms - !> \returns: Solution vector X (overwriten brr) - - implicit none - - real(dp), dimension(:, :), intent(inout) :: arr, brr - - ! local variables - real(dp), dimension(:), allocatable :: work - integer :: n, info, lwork - integer, dimension(size(arr, 1)) :: ipiv - - n = size(arr, 1) - - ! query spacework size - allocate(work(1)) - call DSYSV("U", n, 1, arr, n, ipiv, brr, n, work, -1, info) - call check_lapack_call(info, "DSYSV") - - ! Allocate memory fo the workspace - lwork = max(1, int(work(1))) - deallocate(work) - allocate(work(lwork)) - - ! run linear solver - call DSYSV("U", n, 1, arr, n, ipiv, brr, n, work, lwork, info) - ! If the diagonalization fails due to a singular value try to recover - ! by replacing the 0 value with a tiny one - if (info > 0) then - arr(info, info) = tiny(arr(1, 1)) - call DSYSV("U", n, 1, arr, n, ipiv, brr, n, work, lwork, info) - call check_lapack_call(info, "DSYSV") - end if - - deallocate(work) - - end subroutine lapack_solver - - function lapack_matmul(transA, transB, arr, brr, alpha) result (mtx) - !> perform the matrix multiplication alpha * arr ^ (transA) * brr ^ (transB) - !> see Lapack DGEMM for further details - !> \param transA: 'T' transpose A, 'N' do not tranpose - !> \param transB: 'T' transpose B, 'N' do not tranpose - !> \param arr: first matrix to multiply - !> \param brr: second matrix - !> \param alpha: optional scalar number - !> \return matrix multiplication - - implicit none - - character(len=1), intent(in) :: transA, transB - real(dp), dimension(:, :), intent(in) :: arr, brr - real(dp), optional, intent(in) :: alpha - real(dp), dimension(:, :), allocatable :: mtx - - ! local variables - real(dp) :: x - integer :: m, n, k, lda, ldb - x = 1.d0 - - ! check optional variable - if (present(alpha)) x=alpha - - if (transA == 'T') then - k = size(arr, 1) - m = size(arr, 2) - lda = k - else - k = size(arr, 2) - m = size(arr, 1) - lda = m - end if + end interface - if (transB == 'T') then - n = size(brr, 1) - ldb = n - else - n = size(brr, 2) - ldb = k - end if + ! local variables + integer :: nvec, old_dim - ! resulting array - allocate(mtx(m, n)) - mtx = 0.d0 + ! dimension of the matrices + nvec = size(V,2) + old_dim = size(A_proj,1) - call DGEMM(transA, transB, m, n, k, x, arr, lda, brr, ldb, 0.d0, mtx, m) + ! move to temporal array + allocate(tmp_array(nvec, nvec)) + tmp_array(:old_dim, :old_dim) = A_proj + tmp_array(:,old_dim+1:) = lapack_matmul('T', 'N', V, free_matmul(fun_A, V(:, old_dim+1:))) + tmp_array( old_dim+1:,:old_dim ) = transpose(tmp_array(:old_dim, old_dim+1:)) + ! Move to new expanded matrix + deallocate(A_proj) + call move_alloc(tmp_array, A_proj) - end function lapack_matmul + end subroutine update_projection_free function free_matmul(fun, array) result (mtx) !> \brief perform a matrix-matrix multiplication by generating a matrix on the fly using `fun` @@ -881,44 +473,6 @@ end function fun end function free_matmul - function lapack_matrix_vector(transA, mtx, vector, alpha) result(rs) - !> perform the Matrix vector multiplication alpha * mtx ^ (transA) * vector - !> see DGEMV for details - !> \param transA: 'T' transpose A; 'N' do not transpose - !> \param mtx: matrix to multiply - !> \param vector: vector to multiply - !> \param alpha: optional scalar value - !> \return resulting vector - - implicit none - - character(len=1), intent(in) :: transA - real(dp), dimension(:, :), intent(in) :: mtx - real(dp), dimension(:), intent(in) :: vector - real(dp), optional, intent(in) :: alpha - real(dp), dimension(:), allocatable :: rs - - ! local variable - integer :: m, n - real(dp) :: scalar - scalar = 1.d0 - - ! check optional variable - if (present(alpha)) scalar=alpha - - ! number of row of mtx - m = size(mtx, 1) - n = size(mtx, 2) - - allocate(rs(m)) - rs = 0.d0 - - call DGEMV(transA, m, n, scalar, mtx, m, vector, 1, 0.d0, rs, 1) - - - end function lapack_matrix_vector - - function free_matrix_vector(fun, vector, dim_result) result(rs) !> \brief perform a matrix vector multiplcation computing the matrix on the fly !> \param[in] fun Function to compute a matrix on the fly @@ -956,7 +510,6 @@ end function fun end function free_matrix_vector - subroutine check_deallocate_matrix(mtx) !> deallocate a matrix if allocated real(dp), dimension(:, :), allocatable, intent(inout) :: mtx @@ -966,21 +519,6 @@ subroutine check_deallocate_matrix(mtx) end if end subroutine check_deallocate_matrix - - subroutine check_lapack_call(info, name) - !> Check if a subroutine finishes sucessfully - !> \param info: Termination signal - !> \param name: Name of the subroutine - integer :: info - character(len=*), intent(in) :: name - - if (info /= 0) then - print *, "call to subroutine: ", name, " has failed!" - print *, "info: ", info - error stop - end if - - end subroutine check_lapack_call pure function eye(m, n, alpha) !> Create a matrix with ones in the diagonal and zero everywhere else @@ -1106,6 +644,226 @@ end function diagonal end module davidson +submodule (davidson) davidson_dense + !> Submodule containing the implementation of the Davidson diagonalization method + !> for dense matrices + use lapack_wrapper, only: lapack_generalized_eigensolver, lapack_matmul, lapack_matrix_vector, & + lapack_qr, lapack_solver + + implicit none + + contains + + module 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. + + !> \param[in] mtx: Matrix to diagonalize + !> \param[in, opt] Optional matrix to solve the general eigenvalue problem: + !> \f$ mtx \lambda = V stx \lambda \f$ + !> \param[out] eigenvalues Computed eigenvalues + !> \param[out] ritz_vectors approximation to the eigenvectors + !> \param[in] lowest Number of lowest eigenvalues/ritz_vectors to compute + !> \param[in] method Method to compute the correction vector. Available + !> methods are, + !> DPR: Diagonal-Preconditioned-Residue + !> GJD: Generalized Jacobi Davidson + !> \param[in] max_iters: Maximum number of iterations + !> \param[in] tolerance norm-2 error of the eigenvalues + !> \param[in] method: Method to compute the correction vectors + !> \param[in, opt] max_dim_sub: maximum dimension of the subspace search + !> \param[out] iters: Number of iterations until convergence + !> \return eigenvalues and ritz_vectors of the matrix `mtx` + + implicit none + ! input/output variable + integer, intent(in) :: lowest + real(dp), dimension(:, :), intent(in) :: mtx + real(dp), dimension(:, :), intent(in), optional :: stx + real(dp), dimension(lowest), intent(out) :: eigenvalues + real(dp), dimension(:, :), intent(out) :: ritz_vectors + integer, intent(in) :: max_iters + integer, intent(in), optional :: max_dim_sub + real(dp), intent(in) :: tolerance + character(len=*), intent(in) :: method + integer, intent(out) :: iters + + !local variables + integer :: i, j, dim_sub, max_dim + + ! Basis of subspace of approximants + real(dp), dimension(size(mtx, 1)) :: guess, rs + real(dp), dimension(lowest):: errors + + ! Working arrays + real(dp), dimension(:), allocatable :: eigenvalues_sub + real(dp), dimension(:, :), allocatable :: correction, eigenvectors_sub, mtx_proj, stx_proj, V + + ! generalize problem + logical :: gev + + ! Iteration subpsace dimension + dim_sub = lowest * 2 + + ! maximum dimension of the basis for the subspace + if (present(max_dim_sub)) then + max_dim = max_dim_sub + else + max_dim = lowest * 10 + endif + + ! generalied problem + gev = present(stx) + + ! 1. Variables initialization + V = eye(size(ritz_vectors, 1), dim_sub) ! Initial orthonormal basis + + + ! 2. Generate subpace matrix problem by projecting into V + mtx_proj = lapack_matmul('T', 'N', V, lapack_matmul('N', 'N', mtx, V)) + + if(gev) then + stx_proj = lapack_matmul('T', 'N', V, lapack_matmul('N', 'N', stx, V)) + end if + + ! ! Outer loop block Davidson schema + outer_loop: do i=1, max_iters + + ! 3. compute the eigenvalues and their corresponding ritz_vectors + ! for the projected matrix using lapack + call check_deallocate_matrix(eigenvectors_sub) + + if (allocated(eigenvalues_sub)) then + deallocate(eigenvalues_sub) + end if + + allocate(eigenvalues_sub(size(mtx_proj, 1))) + allocate(eigenvectors_sub(size(mtx_proj, 1), size(mtx_proj, 2))) + + + if (gev) then + call lapack_generalized_eigensolver(mtx_proj, eigenvalues_sub, eigenvectors_sub, stx_proj) + else + call lapack_generalized_eigensolver(mtx_proj, eigenvalues_sub, eigenvectors_sub) + end if + + ! 4. Check for convergence + ritz_vectors = lapack_matmul('N', 'N', V, eigenvectors_sub(:, :lowest)) + do j=1,lowest + if(gev) then + guess = eigenvalues_sub(j) * lapack_matrix_vector('N',stx,ritz_vectors(:, j)) + else + guess = eigenvalues_sub(j) * ritz_vectors(:, j) + end if + rs = lapack_matrix_vector('N', mtx, ritz_vectors(:, j)) - guess + errors(j) = norm(rs) + end do + + if (all(errors < tolerance)) then + iters = i + exit + end if + + ! 5. Add the correction vectors to the current basis + if (size(V, 2) <= max_dim) then + + ! append correction to the current basis + call check_deallocate_matrix(correction) + allocate(correction(size(mtx, 1), size(V, 2))) + + if(gev) then + correction = compute_correction_generalized_dense(mtx, V, eigenvalues_sub, eigenvectors_sub, method, stx) + else + correction = compute_correction_generalized_dense(mtx, V, eigenvalues_sub, eigenvectors_sub, method) + end if + + + ! 6. Increase Basis size + call concatenate(V, correction) + + ! 7. Orthogonalize basis + call lapack_qr(V) + + + ! 8. Update the the projection + call update_projection_dense(mtx, V, mtx_proj) + if (gev) then + call update_projection_dense(stx, V, stx_proj) + end if + + else + + ! 6. Otherwise reduce the basis of the subspace to the current correction + V = lapack_matmul('N', 'N', V, eigenvectors_sub(:, :dim_sub)) + + ! we refresh the projected matrices + mtx_proj = lapack_matmul('T', 'N', V, lapack_matmul('N', 'N', mtx, V)) + + if(gev) then + stx_proj = lapack_matmul('T', 'N', V, lapack_matmul('N', 'N', stx, V)) + end if + + end if + + + end do outer_loop + + ! 8. Check convergence + if (i > max_iters / dim_sub) then + print *, "Warning: Algorithm did not converge!!" + end if + + ! Select the lowest eigenvalues and their corresponding ritz_vectors + ! They are sort in increasing order + eigenvalues = eigenvalues_sub(:lowest) + + ! Free memory + call check_deallocate_matrix(correction) + deallocate(eigenvalues_sub, eigenvectors_sub, V, mtx_proj) + + ! free optional matrix + if (gev) then + call check_deallocate_matrix(stx_proj) + endif + + end subroutine generalized_eigensolver_dense + + subroutine update_projection_dense(A, V, A_proj) + !> \brief update the projected matrices + !> \param A: full matrix + !> \param V: projector + !> \param A_proj: projected matrix + + implicit none + real(dp), dimension(:, :), intent(in) :: A + real(dp), dimension(:, :), intent(in) :: V + real(dp), dimension(:, :), intent(inout), allocatable :: A_proj + real(dp), dimension(:, :), allocatable :: tmp_array + + ! local variables + integer :: nvec, old_dim + + ! dimension of the matrices + nvec = size(V,2) + old_dim = size(A_proj,1) + + ! move to temporal array + allocate(tmp_array(nvec, nvec)) + tmp_array(:old_dim, :old_dim) = A_proj + tmp_array(:,old_dim+1:) = lapack_matmul('T', 'N', V, lapack_matmul('N', 'N', A, V(:, old_dim+1:))) + tmp_array( old_dim+1:,:old_dim ) = transpose(tmp_array(:old_dim, old_dim+1:)) + + ! Move to new expanded matrix + deallocate(A_proj) + call move_alloc(tmp_array, A_proj) + + end subroutine update_projection_dense + +end submodule davidson_dense + submodule (davidson) correction_methods_generalized_dense !> submodule containing the implementations of different kind !> algorithms to compute the correction vectors for the Davidson's diagonalization @@ -1244,15 +1002,3 @@ function substract_from_diagonal(mtx, alpha) result(arr) end function substract_from_diagonal end submodule correction_methods_generalized_dense - - -! submodule (davidson) correction_methods_free -! !> \brief submodule to compute the correction vectors for the matrix free algorithm - -! implicit none - -! contains - - -! end submodule correction_methods_free - diff --git a/src/lapack_wrapper.f90 b/src/lapack_wrapper.f90 new file mode 100644 index 0000000..9015fd4 --- /dev/null +++ b/src/lapack_wrapper.f90 @@ -0,0 +1,298 @@ +module lapack_wrapper + + use numeric_kinds, only: dp + implicit none + + !> \private + private + !> \public + public :: lapack_generalized_eigensolver, lapack_matmul, lapack_matrix_vector, & + lapack_qr, lapack_solver + +contains + + subroutine lapack_generalized_eigensolver(mtx, eigenvalues, eigenvectors, stx) + !> Call the DSYGV subroutine lapack to compute ALL the eigenvalues + !> and corresponding eigenvectors of mtx + !> \param mtx: Matrix to diaogonalize + !> \param stx: Overlap Matrix to diaogonalize + !> \param eigenvalues: lowest eigenvalues + !> \param eigenvectors: corresponding eigenvectors + !> \return eigenvalues/eigenvectors + + ! input/output + implicit none + real(dp), dimension(:, :), intent(in) :: mtx + real(dp), dimension(:, :), intent(in), optional :: stx + real(dp), dimension(size(mtx, 1)), intent(inout) :: eigenvalues + real(dp), dimension(size(mtx, 1), size(mtx, 2)), intent(inout) :: eigenvectors + + real(dp), dimension(:, :), allocatable :: mtx_copy + real(dp), dimension(:, :), allocatable :: stx_copy + + + ! Local variables + integer :: dim, info, lwork, itype = 1 + logical :: gev + + ! ALL the eigenvalues of the subpace (re, im) + real(dp), dimension(size(mtx, 1)) :: eigenvalues_work + real(dp), dimension(:), allocatable :: work ! workspace, see lapack documentation + + ! ! dimension of the guess space + dim = size(mtx, 1) + gev = present(stx) + + ! local copy of the matrices + allocate(mtx_copy(dim,dim)) + mtx_copy = mtx + + if (gev) then + allocate(stx_copy(dim,dim)) + stx_copy = stx + end if + + ! Query size of the optimal workspace + allocate(work(1)) + + if (gev) then + call DSYGV(itype,"V", "U", dim, mtx_copy, dim, stx_copy, dim, eigenvalues_work, work, -1, info) + call check_lapack_call(info, "DSYGV") + else + call DSYEV("V", "U", dim, mtx_copy, dim, eigenvalues_work, work, -1, info) + call check_lapack_call(info, "DSYEV") + end if + + ! Allocate memory for the workspace + lwork = max(1, int(work(1))) + deallocate(work) + allocate(work(lwork)) + + ! Compute Eigenvalues + if (gev) then + call DSYGV(itype,"V", "U", dim, mtx_copy, dim, stx_copy, dim, eigenvalues_work, work, lwork, info) + call check_lapack_call(info, "DSYGV") + else + call DSYEV("V", "U", dim, mtx_copy, dim, eigenvalues_work, work, lwork, info) + call check_lapack_call(info, "DSYEV") + end if + + ! Sort the eigenvalues and eigenvectors of the basis + eigenvalues = eigenvalues_work + eigenvectors = mtx_copy + + ! release memory + deallocate(work) + deallocate(mtx_copy) + if (gev) then + deallocate(stx_copy) + end if + + end subroutine lapack_generalized_eigensolver + + subroutine lapack_qr(basis) + !> Orthoghonalize the basis using the QR factorization. + !> QR factorization of the M-by-N (M>N) matrx A=Q*R in the form where + !> Q is square M-by-M matrix and R is an upper triangular M-by-N matrix. + !> The equality A=Q*R can be re-written also as a product Q1*R1 where Q1 + !> is a rectangular M-by-N submatrix of the matrix Q and R1 is M-by-M + !> submatrix of the R. Let us note that columns of Q1 are orthonormal + !> (they are orthogonal to each other and have norms equal to 1). + !> The equality A=Q1*R1 can be treated as every column of A is a linear + !> combination of Q1 columns, i.e. they span the same linear space. + !> In other words, columns of Q1 is the result of ortogonalization of columns A. + !> DGEQRF does not not compute Q directly, DORGQR must be call subsequently. + + !> \param basis + !> \return orthogonal basis + + implicit none + real(dp), dimension(:, :), intent(inout) :: basis + real(dp), dimension(:), allocatable :: work ! workspace, see lapack documentation + real(dp), dimension(size(basis, 2)) :: tau ! see DGEQRF documentation + integer :: info, lwork, m, n + + ! Matrix shape + m = size(basis, 1) + n = size(basis, 2) + + ! 1. Call the QR decomposition + ! 1.1 Query size of the workspace (Check lapack documentation) + allocate(work(1)) + call DGEQRF(m, n, basis, m, tau, work, -1, info) + call check_lapack_call(info, "DGEQRF") + + ! 1.2 Allocate memory for the workspace + lwork = max(1, int(work(1))) + deallocate(work) + allocate(work(lwork)) + + ! 1.3 Call QR factorization + call DGEQRF(m, n, basis, m, tau, work, lwork, info) + call check_lapack_call(info, "DGEQRF") + deallocate(work) + + ! 2. Generates an orthonormal matrix + ! 2.1 Query size of the workspace (Check lapack documentation) + allocate(work(1)) + call DORGQR(m, n, min(m, n), basis, m, tau, work, -1, info) + call check_lapack_call(info, "DORGQR") + + ! 2.2 Allocate memory fo the workspace + lwork = max(1, int(work(1))) + deallocate(work) + allocate(work(lwork)) + + ! 2.3 compute the matrix Q + call DORGQR(m, n, min(m, n), basis, m, tau, work, lwork, info) + call check_lapack_call(info, "DORGQR") + + ! release memory + deallocate(work) + + end subroutine lapack_qr + + subroutine lapack_solver(arr, brr) + !> Call lapack DSYSV subroutine to solve a AX=B Linear system + !> \param arr: matrix with the coefficients of the linear system + !> \param brr: Vector with the constant terms + !> \returns: Solution vector X (overwriten brr) + + implicit none + + real(dp), dimension(:, :), intent(inout) :: arr, brr + + ! local variables + real(dp), dimension(:), allocatable :: work + integer :: n, info, lwork + integer, dimension(size(arr, 1)) :: ipiv + + n = size(arr, 1) + + ! query spacework size + allocate(work(1)) + call DSYSV("U", n, 1, arr, n, ipiv, brr, n, work, -1, info) + call check_lapack_call(info, "DSYSV") + + ! Allocate memory fo the workspace + lwork = max(1, int(work(1))) + deallocate(work) + allocate(work(lwork)) + + ! run linear solver + call DSYSV("U", n, 1, arr, n, ipiv, brr, n, work, lwork, info) + ! If the diagonalization fails due to a singular value try to recover + ! by replacing the 0 value with a tiny one + if (info > 0) then + arr(info, info) = tiny(arr(1, 1)) + call DSYSV("U", n, 1, arr, n, ipiv, brr, n, work, lwork, info) + call check_lapack_call(info, "DSYSV") + end if + + deallocate(work) + + end subroutine lapack_solver + + function lapack_matmul(transA, transB, arr, brr, alpha) result (mtx) + !> perform the matrix multiplication alpha * arr ^ (transA) * brr ^ (transB) + !> see Lapack DGEMM for further details + !> \param transA: 'T' transpose A, 'N' do not tranpose + !> \param transB: 'T' transpose B, 'N' do not tranpose + !> \param arr: first matrix to multiply + !> \param brr: second matrix + !> \param alpha: optional scalar number + !> \return matrix multiplication + + implicit none + + character(len=1), intent(in) :: transA, transB + real(dp), dimension(:, :), intent(in) :: arr, brr + real(dp), optional, intent(in) :: alpha + real(dp), dimension(:, :), allocatable :: mtx + + ! local variables + real(dp) :: x + integer :: m, n, k, lda, ldb + x = 1.d0 + + ! check optional variable + if (present(alpha)) x=alpha + + if (transA == 'T') then + k = size(arr, 1) + m = size(arr, 2) + lda = k + else + k = size(arr, 2) + m = size(arr, 1) + lda = m + end if + + if (transB == 'T') then + n = size(brr, 1) + ldb = n + else + n = size(brr, 2) + ldb = k + end if + + ! resulting array + allocate(mtx(m, n)) + mtx = 0.d0 + + call DGEMM(transA, transB, m, n, k, x, arr, lda, brr, ldb, 0.d0, mtx, m) + + end function lapack_matmul + + function lapack_matrix_vector(transA, mtx, vector, alpha) result(rs) + !> perform the Matrix vector multiplication alpha * mtx ^ (transA) * vector + !> see DGEMV for details + !> \param transA: 'T' transpose A; 'N' do not transpose + !> \param mtx: matrix to multiply + !> \param vector: vector to multiply + !> \param alpha: optional scalar value + !> \return resulting vector + + implicit none + + character(len=1), intent(in) :: transA + real(dp), dimension(:, :), intent(in) :: mtx + real(dp), dimension(:), intent(in) :: vector + real(dp), optional, intent(in) :: alpha + real(dp), dimension(:), allocatable :: rs + + ! local variable + integer :: m, n + real(dp) :: scalar + scalar = 1.d0 + + ! check optional variable + if (present(alpha)) scalar=alpha + + ! number of row of mtx + m = size(mtx, 1) + n = size(mtx, 2) + + allocate(rs(m)) + rs = 0.d0 + + call DGEMV(transA, m, n, scalar, mtx, m, vector, 1, 0.d0, rs, 1) + + end function lapack_matrix_vector + + subroutine check_lapack_call(info, name) + !> Check if a subroutine finishes sucessfully + !> \param info: Termination signal + !> \param name: Name of the subroutine + integer :: info + character(len=*), intent(in) :: name + + if (info /= 0) then + print *, "call to subroutine: ", name, " has failed!" + print *, "info: ", info + error stop + end if + + end subroutine check_lapack_call + +end module lapack_wrapper diff --git a/src/main.f90 b/src/main.f90 index fbcdb11..fe6b103 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -69,7 +69,8 @@ end module test_utils program main use numeric_kinds, only: dp - use davidson, only: generalized_eigensolver, norm, lapack_generalized_eigensolver, generate_diagonal_dominant + use davidson, only: generalized_eigensolver, norm, generate_diagonal_dominant + use lapack_wrapper, only: lapack_generalized_eigensolver use test_utils, only: diagonal , read_matrix, write_vector, cast_to_double use benchmark, only: compute_benchmark diff --git a/src/tests/CMakeLists.txt b/src/tests/CMakeLists.txt index 256d9b1..0ba11da 100644 --- a/src/tests/CMakeLists.txt +++ b/src/tests/CMakeLists.txt @@ -11,7 +11,7 @@ execute_process( OUTPUT_STRIP_TRAILING_WHITESPACE ) -set(NUMPY_TEST_SRC ${CMAKE_SOURCE_DIR}/src/davidson.f90 ${CMAKE_SOURCE_DIR}/src/numeric_kinds.f90) +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) diff --git a/src/tests/test_call_lapack.f90 b/src/tests/test_call_lapack.f90 index d68114e..6f17c26 100644 --- a/src/tests/test_call_lapack.f90 +++ b/src/tests/test_call_lapack.f90 @@ -1,7 +1,8 @@ program main use numeric_kinds, only: dp use test_utils, only: write_matrix, write_vector - use davidson, only: generate_diagonal_dominant, lapack_generalized_eigensolver, lapack_qr + use davidson, only: generate_diagonal_dominant + use lapack_wrapper, only: lapack_generalized_eigensolver, lapack_qr implicit none diff --git a/src/tests/test_dense_properties.f90 b/src/tests/test_dense_properties.f90 index a895298..cfd3d3a 100644 --- a/src/tests/test_dense_properties.f90 +++ b/src/tests/test_dense_properties.f90 @@ -1,7 +1,7 @@ program main use numeric_kinds, only: dp - use davidson, only: generalized_eigensolver, norm, lapack_generalized_eigensolver, generate_diagonal_dominant + use davidson, only: generalized_eigensolver, norm, generate_diagonal_dominant use test_utils, only: diagonal implicit none diff --git a/src/tests/test_free_properties.f90 b/src/tests/test_free_properties.f90 index 79a1360..4dc7b5f 100644 --- a/src/tests/test_free_properties.f90 +++ b/src/tests/test_free_properties.f90 @@ -1,7 +1,7 @@ program main use numeric_kinds, only: dp - use davidson, only: generalized_eigensolver, norm, lapack_generalized_eigensolver, generate_diagonal_dominant + use davidson, only: generalized_eigensolver, norm, generate_diagonal_dominant use test_utils, only: compute_matrix_on_the_fly, compute_stx_on_the_fly, diagonal, write_matrix implicit none