diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index f05eba6..0b639a6 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,7 +1,7 @@ # create binary -set(UTILS numeric_kinds.f90 lapack_wrapper.f90) -set(SOURCES main.f90 davidson.f90 benchmark.f90 ${UTILS}) +set(UTILS array_utils.f90 numeric_kinds.f90 lapack_wrapper.f90) +set(SOURCES main.f90 davidson.f90 ${UTILS}) message (STATUS "SOURCES: " ${SOURCES}) add_executable(main ${SOURCES}) diff --git a/src/array_utils.f90 b/src/array_utils.f90 new file mode 100644 index 0000000..771651c --- /dev/null +++ b/src/array_utils.f90 @@ -0,0 +1,135 @@ +module array_utils + + 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 :: concatenate, diagonal,eye, generate_diagonal_dominant, norm + +contains + + pure function eye(m, n, alpha) + !> Create a matrix with ones in the diagonal and zero everywhere else + !> \param m: number of rows + !> \param n: number of colums + !> \param alpha: optional diagonal value + !> \return matrix of size n x m + integer, intent(in) :: n, m + real(dp), dimension(m, n) :: eye + real(dp), intent(in), optional :: alpha + + !local variable + integer :: i, j + real(dp) :: x + + ! check optional values + x = 1.d0 + if (present(alpha)) x = alpha + + do i=1, m + do j=1, n + if (i /= j) then + eye(i, j) = 0 + else + eye(i, i) = x + end if + end do + end do + + end function eye + + pure function norm(vector) + !> compute the norm-2 of a vector + real(dp), dimension(:), intent(in) :: vector + real(dp) :: norm + + norm = sqrt(sum(vector ** 2.d0)) + + end function norm + + subroutine concatenate(arr, brr) + + !> Concatenate two matrices + !> \param arr: first array + !> \param brr: second array + !> \return arr concatenate brr (overwrites arr) + + real(dp), dimension(:, :), intent(inout), allocatable :: arr + real(dp), dimension(:, :), intent(in) :: brr + real(dp), dimension(:, :), allocatable :: tmp_array + integer :: new_dim, dim_cols, dim_rows + + ! dimension + dim_rows = size(arr, 1) + dim_cols = size(arr, 2) + + ! Number of columns of the new matrix + new_dim = dim_cols + size(brr, 2) + + ! move to temporal array + allocate(tmp_array(dim_rows, new_dim)) + tmp_array(:, :dim_cols) = arr + + ! Move to new expanded matrix + deallocate(arr) + call move_alloc(tmp_array, arr) + + arr(:, dim_cols + 1:) = brr + + end subroutine concatenate + + function generate_diagonal_dominant(m, sparsity, diag_val) result(arr) + !> Generate a diagonal dominant square matrix of dimension m + !> \param m dimension of the matrix + !> \param sparsity magnitude order of the off-diagonal values + + integer, intent(in) :: m ! size of the square matrix + real(dp), optional :: diag_val + integer :: i, j + real(dp) :: sparsity + real(dp), dimension(m, m) :: arr + call random_number(arr) + + arr = arr * sparsity + do j=1, m + do i=1, m + if (i > j) then + arr(i, j) = arr(j, i) + else if(i == j) then + if (present(diag_val))then + arr(i,i) = diag_val + else + arr(i, i) = i + end if + end if + end do + end do + + end function generate_diagonal_dominant + + function diagonal(matrix) + !> return the diagonal of a matrix + real(dp), dimension(:, :), intent(in) :: matrix + real(dp), dimension(size(matrix, 1)) :: diagonal + + ! local variables + integer :: i, j, m + + ! dimension of the matrix + m = size(matrix, 1) + + do i=1,m + do j=1,m + if (i == j) then + diagonal(i) = matrix(i, j) + end if + end do + end do + + end function diagonal + +end module array_utils diff --git a/src/benchmark.f90 b/src/benchmark.f90 deleted file mode 100644 index 4cd38a0..0000000 --- a/src/benchmark.f90 +++ /dev/null @@ -1,94 +0,0 @@ -module benchmark - - use numeric_kinds, only: dp - use Davidson, only: generalized_eigensolver, generate_diagonal_dominant - - implicit none - -contains - - subroutine compute_benchmark(dims, lowest, sparsity, times, iters) - !> Benchmark the Davidson algorithm - !> \param dims: Vector of integer with the dimension of the matrices to test - !> \param lowest: number of eigenvalues/eigenvectors to compute - !> \param sparsity: magnitud of the off-diagonal terms - !> \return times: resulting times - !> \param number of iterations - - integer, dimension(:), intent(in) :: dims - integer, intent(in) :: lowest - integer, dimension(size(dims), 2), intent(out) :: iters - real(dp), dimension(size(dims), 2), intent(out) :: times - real(dp) :: sparsity - - ! local variable - integer :: i, iter_i - real(dp), dimension(:), allocatable :: eigenvalues - real(dp), dimension(:, :), allocatable :: mtx, eigenvectors - real(dp) :: dt - - - do i=1, size(dims) - call benchmark_method(mtx, eigenvalues, eigenvectors, "DPR", dims(i), lowest, sparsity, dt, iter_i) - print *, "cycles: ", iter_i - iters(i, 1) = iter_i - times(i, 1) = dt - - call benchmark_method(mtx, eigenvalues, eigenvectors, "GJD", dims(i), lowest, sparsity, dt, iter_i) - print *, "cycles: ", iter_i - iters(i, 2) = iter_i - times(i, 2) = dt - - end do - - call free_space(mtx, eigenvalues, eigenvectors) - - end subroutine compute_benchmark - - subroutine benchmark_method(mtx, eigenvalues, eigenvectors, method, dim, lowest, sparsity, dt, iter_i) - !> Benchmark different correction methods - real(dp), dimension(:), allocatable :: eigenvalues - real(dp), dimension(:, :), allocatable :: mtx, eigenvectors - real(dp), intent(out) :: dt - integer, intent(in) :: dim, lowest - integer, intent(out) :: iter_i - character(len=*), intent(in) :: method - real(dp) :: t1, t2, sparsity - - call free_space(mtx, eigenvalues, eigenvectors) - ! generate test matrix - mtx = generate_diagonal_dominant(dim, sparsity) - ! diagonalize - allocate(eigenvalues(lowest)) - allocate(eigenvectors(dim, lowest)) - print *, "Davidson method: ", method, " dimension: ", dim - call cpu_time(t1) - call generalized_eigensolver(mtx, eigenvalues, eigenvectors, 3, method, 1000, 1d-8, iter_i) - call cpu_time(t2) - dt = t2 - t1 - print *, "time: ", dt - - - end subroutine benchmark_method - - subroutine free_space(mtx, es, vs) - !> Deallocate variables - real(dp), dimension(:, :), allocatable, intent(inout) :: mtx, vs - real(dp), dimension(:), allocatable, intent(inout) :: es - - if (allocated(mtx)) then - deallocate(mtx) - end if - - if (allocated(vs)) then - deallocate(vs) - end if - - if (allocated(es)) then - deallocate(es) - end if - - end subroutine free_space - - -end module benchmark diff --git a/src/benchmark_free.f90 b/src/benchmark_free.f90 index cb7607d..be33490 100644 --- a/src/benchmark_free.f90 +++ b/src/benchmark_free.f90 @@ -51,7 +51,8 @@ end module matrix_free program main use numeric_kinds, only: dp - use davidson, only: generalized_eigensolver, generate_diagonal_dominant, norm + use davidson, only: generalized_eigensolver + use array_utils, only: generate_diagonal_dominant, norm use matrix_free, only: compute_matrix_on_the_fly, compute_stx_on_the_fly implicit none diff --git a/src/davidson.f90 b/src/davidson.f90 index ed1d436..785ce14 100644 --- a/src/davidson.f90 +++ b/src/davidson.f90 @@ -1,84 +1,16 @@ !> \namespace Davidson eigensolver !> \author Felipe Zapata -module davidson - +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 - implicit none - - !> \private - private - !> \public - public :: eye, generalized_eigensolver, generate_diagonal_dominant, norm + use array_utils, only: concatenate, eye, norm - interface generalized_eigensolver - !> \brief Solve a (general) eigenvalue problem using different types of Davidson algorithms. - - !> \param[in] mtx: Matrix to diagonalize - !> \param[in, opt] stx: Optional matrix for 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 - module procedure generalized_eigensolver_dense - module procedure generalized_eigensolver_free - - end interface generalized_eigensolver - 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. - - !> \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 - - 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 @@ -89,7 +21,7 @@ module function compute_correction_generalized_dense(mtx, V, eigenvalues, eigenv !> \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 @@ -97,22 +29,21 @@ module function compute_correction_generalized_dense(mtx, V, eigenvalues, eigenv 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 - !> the general eigenvalue problem + contains + + 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] fun_mtx: Function to compute the Matrix to diagonalize - !> \param[in, opt] fun_stx: function to compute the optional general eigenvalue problem. + !> \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 @@ -126,9 +57,12 @@ subroutine generalized_eigensolver_free(fun_mtx, eigenvalues, ritz_vectors, lowe !> \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 @@ -137,85 +71,77 @@ subroutine generalized_eigensolver_free(fun_mtx, eigenvalues, ritz_vectors, lowe 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) - !> \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 :: dim_mtx, dim_sub, max_dim, i, j + integer :: i, j, dim_sub, max_dim - ! ! Basis of subspace of approximants - real(dp), dimension(size(ritz_vectors, 1)) :: guess, rs + ! Basis of subspace of approximants + real(dp), dimension(size(mtx, 1)) :: guess, rs real(dp), dimension(lowest):: errors - - ! ! Working arrays + + ! 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 - - ! dimension of the matrix - dim_mtx = size(ritz_vectors, 1) - + + ! generalied problem + gev = present(stx) + ! 1. Variables initialization - V = eye(dim_mtx, dim_sub) ! Initial orthonormal basis - - ! 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 + 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))) - - call lapack_generalized_eigensolver(mtx_proj, eigenvalues_sub, eigenvectors_sub, stx_proj) - + + + 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 - guess = eigenvalues_sub(j) * free_matrix_vector(fun_stx, ritz_vectors(:, j), dim_mtx) - rs = free_matrix_vector(fun_mtx, ritz_vectors(:, j), dim_mtx) - guess + 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 @@ -223,43 +149,53 @@ 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) - + 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_free(fun_mtx, V, mtx_proj) - call update_projection_free(fun_stx, V, stx_proj) - + 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, free_matmul(fun_mtx, V)) - stx_proj = lapack_matmul('T', 'N', V, free_matmul(fun_stx, 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 + 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) @@ -267,37 +203,270 @@ 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) + if (gev) then + call check_deallocate_matrix(stx_proj) + endif - end subroutine generalized_eigensolver_free + 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:)) - 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 + ! Move to new expanded matrix + deallocate(A_proj) + call move_alloc(tmp_array, A_proj) + + end subroutine update_projection_dense + + subroutine check_deallocate_matrix(mtx) + !> deallocate a matrix if allocated + real(dp), dimension(:, :), allocatable, intent(inout) :: mtx + + if (allocated(mtx)) then + deallocate(mtx) + end if + + end subroutine check_deallocate_matrix + +end module davidson_dense + + +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 davidson_dense, only: generalized_eigensolver_dense + implicit none + + !> \private + private + !> \public + public :: generalized_eigensolver_free + +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 + !> the general eigenvalue problem + !> 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] 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 + !> \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(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 + + ! Function to compute the target matrix on the fly + 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 :: dim_mtx, dim_sub, max_dim, i, j + + ! ! 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 + 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 + ! 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))) + + 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 + guess = eigenvalues_sub(j) * free_matrix_vector(fun_stx, ritz_vectors(:, j), dim_mtx) + 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 + 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(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 + eigenvalues = eigenvalues_sub(:lowest) + + ! 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 @@ -517,354 +686,58 @@ subroutine check_deallocate_matrix(mtx) if (allocated(mtx)) then deallocate(mtx) end if - - end subroutine check_deallocate_matrix - - pure function eye(m, n, alpha) - !> Create a matrix with ones in the diagonal and zero everywhere else - !> \param m: number of rows - !> \param n: number of colums - !> \param alpha: optional diagonal value - !> \return matrix of size n x m - integer, intent(in) :: n, m - real(dp), dimension(m, n) :: eye - real(dp), intent(in), optional :: alpha - - !local variable - integer :: i, j - real(dp) :: x - - ! check optional values - x = 1.d0 - if (present(alpha)) x = alpha - - do i=1, m - do j=1, n - if (i /= j) then - eye(i, j) = 0 - else - eye(i, i) = x - end if - end do - end do - - end function eye - - pure function norm(vector) - !> compute the norm-2 of a vector - real(dp), dimension(:), intent(in) :: vector - real(dp) :: norm - - norm = sqrt(sum(vector ** 2.d0)) - - end function norm - - subroutine concatenate(arr, brr) - - !> Concatenate two matrices - !> \param arr: first array - !> \param brr: second array - !> \return arr concatenate brr (overwrites arr) - - real(dp), dimension(:, :), intent(inout), allocatable :: arr - real(dp), dimension(:, :), intent(in) :: brr - real(dp), dimension(:, :), allocatable :: tmp_array - integer :: new_dim, dim_cols, dim_rows - - ! dimension - dim_rows = size(arr, 1) - dim_cols = size(arr, 2) - - ! Number of columns of the new matrix - new_dim = dim_cols + size(brr, 2) - - ! move to temporal array - allocate(tmp_array(dim_rows, new_dim)) - tmp_array(:, :dim_cols) = arr - - ! Move to new expanded matrix - deallocate(arr) - call move_alloc(tmp_array, arr) - - arr(:, dim_cols + 1:) = brr - - end subroutine concatenate - function generate_diagonal_dominant(m, sparsity, diag_val) result(arr) - !> Generate a diagonal dominant square matrix of dimension m - !> \param m dimension of the matrix - !> \param sparsity magnitude order of the off-diagonal values - - integer, intent(in) :: m ! size of the square matrix - real(dp), optional :: diag_val - integer :: i, j - real(dp) :: sparsity - real(dp), dimension(m, m) :: arr - call random_number(arr) - - arr = arr * sparsity - do j=1, m - do i=1, m - if (i > j) then - arr(i, j) = arr(j, i) - else if(i == j) then - if (present(diag_val))then - arr(i,i) = diag_val - else - arr(i, i) = i - end if - end if - end do - end do - - end function generate_diagonal_dominant - - function diagonal(matrix) - !> return the diagonal of a matrix - real(dp), dimension(:, :), intent(in) :: matrix - real(dp), dimension(size(matrix, 1)) :: diagonal - - ! local variables - integer :: i, j, m + end subroutine check_deallocate_matrix - ! dimension of the matrix - m = size(matrix, 1) - - do i=1,m - do j=1,m - if (i == j) then - diagonal(i) = matrix(i, j) - end if - end do - end do - end function diagonal +end module davidson_free - -end module davidson +module davidson -submodule (davidson) 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 + use davidson_dense, only: generalized_eigensolver_dense + use davidson_free, only: generalized_eigensolver_free implicit none - contains + !> \private + private + !> \public + public :: generalized_eigensolver - 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 + interface generalized_eigensolver + !> \brief Solve a (general) eigenvalue problem using different types of Davidson algorithms. - 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 + !> \param[in] mtx: Matrix to diagonalize + !> \param[in, opt] stx: Optional matrix for 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` - ! local variables - integer :: nvec, old_dim + procedure generalized_eigensolver_dense + procedure generalized_eigensolver_free + + end interface generalized_eigensolver - ! dimension of the matrices - nvec = size(V,2) - old_dim = size(A_proj,1) +end module davidson - ! 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 (davidson_dense) correction_methods_generalized_dense !> submodule containing the implementations of different kind !> algorithms to compute the correction vectors for the Davidson's diagonalization diff --git a/src/main.f90 b/src/main.f90 index fe6b103..dbec38e 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -3,45 +3,6 @@ module test_utils implicit none contains - function diagonal(matrix) - !> return the diagonal of a matrix - real(dp), dimension(:, :), intent(in) :: matrix - real(dp), dimension(size(matrix, 1)) :: diagonal - - ! local variables - integer :: i, j, m - - ! dimension of the matrix - m = size(matrix, 1) - - do i=1,m - do j=1,m - if (i == j) then - diagonal(i) = matrix(i, j) - end if - end do - end do - - end function diagonal - - function read_matrix(path_file, dim) result(mtx) - !> read a row-major square matrix from a file - !> \param path_file: path to the file - !> \param dim: dimension of the square matrix - !> \return matrix - - character(len=*), intent(in) :: path_file - integer, intent(in) :: dim - real(dp), dimension(dim, dim) :: mtx - integer :: i - - open(unit=3541, file=path_file, status="OLD") - do i=1,dim - read(3541, *) mtx(i, :) - end do - close(3541) - - end function read_matrix subroutine write_vector(path_file, vector) !> Write vector to path_file @@ -69,27 +30,24 @@ end module test_utils program main use numeric_kinds, only: dp - use davidson, only: generalized_eigensolver, norm, generate_diagonal_dominant + use davidson, only: generalized_eigensolver use lapack_wrapper, only: lapack_generalized_eigensolver - use test_utils, only: diagonal , read_matrix, write_vector, cast_to_double - use benchmark, only: compute_benchmark + use array_utils, only: norm, generate_diagonal_dominant + use test_utils, only: write_vector, cast_to_double implicit none + integer, parameter :: dim = 50 real(dp), dimension(3) :: eigenvalues_DPR, eigenvalues_GJD - real(dp), dimension(50, 3) :: eigenvectors_DPR, eigenvectors_GJD - real(dp), dimension(50, 50) :: mtx - real(dp), dimension(50, 50) :: stx - real(dp) :: test_norm_eigenvalues, sparsity - real(dp), dimension(50) :: xs - real(dp), dimension(5, 2) :: times - integer, dimension(5, 2) :: iters - integer, dimension(5) :: dims + real(dp), dimension(dim, 3) :: eigenvectors_DPR, eigenvectors_GJD + real(dp), dimension(dim, dim) :: mtx + real(dp), dimension(dim, dim) :: stx + real(dp) :: test_norm_eigenvalues + real(dp), dimension(dim) :: xs integer :: iter_i, j - character(len=20) :: arg1 - mtx = generate_diagonal_dominant(50, 1d-4) - stx = generate_diagonal_dominant(50, 1d-4, 1d0) + 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) @@ -109,22 +67,6 @@ program main 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 - end do - - ! Run benchmark - call get_command_argument(1, arg1) - if (arg1 == "benchmark") then - print *, "Running Benchmark! " - dims = [10, 50, 100, 500, 1000] - sparsity = 1d-3 - call compute_benchmark(dims, 3, sparsity, times, iters) - - call write_vector("times_DPR.txt", times(:, 1)) - call write_vector("times_GJD.txt", times(:, 2)) - - call write_vector("cycles_DPR.txt", cast_to_double(iters(:, 1))) - call write_vector("cycles_GJD.txt", cast_to_double(iters(:, 2))) - end if - + end do end program main diff --git a/src/tests/CMakeLists.txt b/src/tests/CMakeLists.txt index 0ba11da..3d36e19 100644 --- a/src/tests/CMakeLists.txt +++ b/src/tests/CMakeLists.txt @@ -11,7 +11,11 @@ execute_process( OUTPUT_STRIP_TRAILING_WHITESPACE ) -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(SOURCE_DAVIDSON array_utils davidson lapack_wrapper numeric_kinds) +foreach(PROG ${SOURCE_DAVIDSON}) + list(APPEND NUMPY_TEST_SRC ${CMAKE_SOURCE_DIR}/src/${PROG}.f90) +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) diff --git a/src/tests/test_call_lapack.f90 b/src/tests/test_call_lapack.f90 index 6f17c26..746c0d2 100644 --- a/src/tests/test_call_lapack.f90 +++ b/src/tests/test_call_lapack.f90 @@ -1,7 +1,7 @@ program main use numeric_kinds, only: dp + use array_utils, only: generate_diagonal_dominant use test_utils, only: write_matrix, write_vector - use davidson, only: generate_diagonal_dominant use lapack_wrapper, only: lapack_generalized_eigensolver, lapack_qr implicit none diff --git a/src/tests/test_dense_numpy.f90 b/src/tests/test_dense_numpy.f90 index b8d56cd..59098be 100644 --- a/src/tests/test_dense_numpy.f90 +++ b/src/tests/test_dense_numpy.f90 @@ -1,6 +1,7 @@ program main use numeric_kinds, only: dp - use davidson, only: generalized_eigensolver, generate_diagonal_dominant + use davidson, only: generalized_eigensolver + use array_utils, only: generate_diagonal_dominant use test_utils, only: write_matrix, write_vector implicit none diff --git a/src/tests/test_dense_properties.f90 b/src/tests/test_dense_properties.f90 index cfd3d3a..236159d 100644 --- a/src/tests/test_dense_properties.f90 +++ b/src/tests/test_dense_properties.f90 @@ -1,8 +1,8 @@ program main use numeric_kinds, only: dp - use davidson, only: generalized_eigensolver, norm, generate_diagonal_dominant - use test_utils, only: diagonal + use davidson, only: generalized_eigensolver + use array_utils, only: diagonal, norm, generate_diagonal_dominant implicit none diff --git a/src/tests/test_free_numpy.f90 b/src/tests/test_free_numpy.f90 index 999eea4..6420e97 100644 --- a/src/tests/test_free_numpy.f90 +++ b/src/tests/test_free_numpy.f90 @@ -1,6 +1,7 @@ program main use numeric_kinds, only: dp - use davidson, only: generalized_eigensolver, generate_diagonal_dominant + use davidson, only: generalized_eigensolver + use array_utils, only: generate_diagonal_dominant use test_utils, only: compute_matrix_on_the_fly, compute_stx_on_the_fly, write_matrix, write_vector implicit none diff --git a/src/tests/test_free_properties.f90 b/src/tests/test_free_properties.f90 index 4dc7b5f..5e85d0d 100644 --- a/src/tests/test_free_properties.f90 +++ b/src/tests/test_free_properties.f90 @@ -1,8 +1,9 @@ program main use numeric_kinds, only: dp - 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 + use davidson, only: generalized_eigensolver + use array_utils, only: diagonal, norm, generate_diagonal_dominant + use test_utils, only: compute_matrix_on_the_fly, compute_stx_on_the_fly, write_matrix implicit none diff --git a/src/tests/test_utils.f90 b/src/tests/test_utils.f90 index 0816458..791f0c7 100644 --- a/src/tests/test_utils.f90 +++ b/src/tests/test_utils.f90 @@ -1,7 +1,7 @@ module test_utils use numeric_kinds, only: dp - use davidson, only: generate_diagonal_dominant + use array_utils, only: generate_diagonal_dominant implicit none contains @@ -85,30 +85,7 @@ function compute_stx_on_the_fly(i, dim) result (vector) vector(i) = 1d0 end function compute_stx_on_the_fly - - - function diagonal(matrix) - !> return the diagonal of a matrix - real(dp), dimension(:, :), intent(in) :: matrix - real(dp), dimension(size(matrix, 1)) :: diagonal - - ! local variables - integer :: i, j, m - - ! dimension of the matrix - m = size(matrix, 1) - - do i=1,m - do j=1,m - if (i == j) then - diagonal(i) = matrix(i, j) - end if - end do - end do - - end function diagonal - function read_matrix(path_file, dim) result(mtx) !> read a row-major square matrix from a file !> \param path_file: path to the file