diff --git a/example/ginzburg_landau/main.f90 b/example/ginzburg_landau/main.f90 index 7b57b2d..9435d37 100644 --- a/example/ginzburg_landau/main.f90 +++ b/example/ginzburg_landau/main.f90 @@ -42,6 +42,9 @@ program demo !----- INITIALIZATION ----- !---------------------------------- + !> Set up logging + call logger_setup() + !> Initialize physical parameters. call initialize_parameters() diff --git a/fpm.toml b/fpm.toml index 1e44bc6..ddfec6d 100644 --- a/fpm.toml +++ b/fpm.toml @@ -11,6 +11,7 @@ auto-executables = true auto-tests = true auto-examples = true module-naming = false +external-modules = "mpi_f08" [install] library = true @@ -22,6 +23,7 @@ source-form = "free" [dependencies] stdlib = "*" +FACE = {git="https://github.com/szaghi/FACE.git"} [dev-dependencies] test-drive.git = "https://github.com/nekStab/test-drive.git" diff --git a/src/BaseKrylov.f90 b/src/BaseKrylov.f90 index 2eb6459..f1aff72 100644 --- a/src/BaseKrylov.f90 +++ b/src/BaseKrylov.f90 @@ -2,7 +2,7 @@ module lightkrylov_BaseKrylov use iso_fortran_env use stdlib_optval, only: optval use stdlib_linalg, only: eye - use LightKrylov_constants + use LightKrylov_Constants use LightKrylov_Logger use LightKrylov_Utils use LightKrylov_AbstractVectors @@ -11,12 +11,13 @@ module lightkrylov_BaseKrylov private character(len=128), parameter :: this_module = 'LightKrylov_BaseKrylov' - public :: qr public :: apply_permutation_matrix public :: apply_inverse_permutation_matrix public :: arnoldi public :: initialize_krylov_subspace + public :: is_orthonormal + public :: orthonormalize_basis public :: orthogonalize_against_basis public :: lanczos_bidiagonalization public :: lanczos_tridiagonalization @@ -79,6 +80,20 @@ module lightkrylov_BaseKrylov module procedure initialize_krylov_subspace_cdp end interface + interface is_orthonormal + module procedure is_orthonormal_rsp + module procedure is_orthonormal_rdp + module procedure is_orthonormal_csp + module procedure is_orthonormal_cdp + end interface + + interface orthonormalize_basis + module procedure orthonormalize_basis_rsp + module procedure orthonormalize_basis_rdp + module procedure orthonormalize_basis_csp + module procedure orthonormalize_basis_cdp + end interface + interface orthogonalize_against_basis module procedure orthogonalize_vector_against_basis_rsp module procedure orthogonalize_basis_against_basis_rsp @@ -145,15 +160,54 @@ end function eigvals_select_dp !----- UTILITY FUNCTIONS ----- !------------------------------------- + logical function is_orthonormal_rsp(X) result(ortho) + class(abstract_vector_rsp), intent(in) :: X(:) + real(sp), dimension(size(X), size(X)) :: G + ortho = .true. + call innerprod(G, X, X) + if (norm2(abs(G - eye(size(X)))) > rtol_sp) then + ! The basis is not orthonormal. Cannot orthonormalize. + ortho = .false. + end if + end function + logical function is_orthonormal_rdp(X) result(ortho) + class(abstract_vector_rdp), intent(in) :: X(:) + real(dp), dimension(size(X), size(X)) :: G + ortho = .true. + call innerprod(G, X, X) + if (norm2(abs(G - eye(size(X)))) > rtol_sp) then + ! The basis is not orthonormal. Cannot orthonormalize. + ortho = .false. + end if + end function + logical function is_orthonormal_csp(X) result(ortho) + class(abstract_vector_csp), intent(in) :: X(:) + complex(sp), dimension(size(X), size(X)) :: G + ortho = .true. + call innerprod(G, X, X) + if (norm2(abs(G - eye(size(X)))) > rtol_sp) then + ! The basis is not orthonormal. Cannot orthonormalize. + ortho = .false. + end if + end function + logical function is_orthonormal_cdp(X) result(ortho) + class(abstract_vector_cdp), intent(in) :: X(:) + complex(dp), dimension(size(X), size(X)) :: G + ortho = .true. + call innerprod(G, X, X) + if (norm2(abs(G - eye(size(X)))) > rtol_sp) then + ! The basis is not orthonormal. Cannot orthonormalize. + ortho = .false. + end if + end function + subroutine initialize_krylov_subspace_rsp(X, X0) class(abstract_vector_rsp), intent(inout) :: X(:) class(abstract_vector_rsp), optional, intent(in) :: X0(:) ! Internal variables. class(abstract_vector_rsp), allocatable :: Q(:) - real(sp), allocatable :: R(:, :) - integer, allocatable :: perm(:) - integer :: i, p, info + integer :: p ! Zero-out X. call zero_basis(X) @@ -161,28 +215,34 @@ subroutine initialize_krylov_subspace_rsp(X, X0) ! Deals with optional args. if(present(X0)) then p = size(X0) - ! Initialize. - allocate(Q(1:p), source=X0(1:p)) - do i = 1, p - call Q(i)%zero() - call Q(i)%add(X0(i)) - enddo - + allocate(Q(p), source=X0) ! Orthonormalize. - allocate(R(p, p)) ; R = zero_rsp - allocate(perm(p)) ; perm = 0 - call qr(Q, R, perm, info) - do i = 1, p - call X(i)%add(Q(i)) - enddo + call orthonormalize_basis(Q) + call copy_basis(X(:p), Q) endif return end subroutine initialize_krylov_subspace_rsp + subroutine orthonormalize_basis_rsp(X) + !! Orthonormalizes the `abstract_vector` basis `X` + class(abstract_vector_rsp), intent(inout) :: X(:) + !! Input `abstract_vector` basis to orthogonalize against + + ! internals + real(sp) :: R(size(X),size(X)) + integer :: info + + ! internals + call qr(X, R, info) + call check_info(info, 'qr', module=this_module, procedure='orthonormalize_basis_rsp') + + return + end subroutine orthonormalize_basis_rsp + subroutine orthogonalize_vector_against_basis_rsp(y, X, info, if_chk_orthonormal, beta) - !! Orthonormalizes the `abstract_vector` `y` against a basis `X` of `abstract_vector`. + !! Orthogonalizes the `abstract_vector` `y` against a basis `X` of `abstract_vector`. class(abstract_vector_rsp), intent(inout) :: y !! Input `abstract_vector` to orthogonalize class(abstract_vector_rsp), intent(in) :: X(:) @@ -239,7 +299,7 @@ subroutine orthogonalize_vector_against_basis_rsp(y, X, info, if_chk_orthonormal end subroutine orthogonalize_vector_against_basis_rsp subroutine orthogonalize_basis_against_basis_rsp(Y, X, info, if_chk_orthonormal, beta) - !! Orthonormalizes the `abstract_vector` basis `Y` against a basis `X` of `abstract_vector`. + !! Orthogonalizes the `abstract_vector` basis `Y` against a basis `X` of `abstract_vector`. class(abstract_vector_rsp), intent(inout) :: Y(:) !! Input `abstract_vector` basis to orthogonalize class(abstract_vector_rsp), intent(in) :: X(:) @@ -326,11 +386,12 @@ subroutine DGS_vector_against_basis_rsp(y, X, info, if_chk_orthonormal, beta) ! Orthogonalize vector y w.r.t. to Krylov basis X in two passes of GS. ! first pass call orthogonalize_against_basis(y, X, info, if_chk_orthonormal=.false., beta=proj_coefficients) - call check_info(info, 'orthogonalize_against_basis', module=this_module, procedure='DGS_vector_against_basis, first pass') + call check_info(info, 'orthogonalize_against_basis', module=this_module, & + & procedure='DGS_vector_against_basis_rsp, pass 1') ! second pass call orthogonalize_against_basis(y, X, info, if_chk_orthonormal=.false., beta=wrk) - call check_info(info, 'orthogonalize_against_basis', module=this_module, procedure='DGS_vector_against_basis_rsp, second& - & pass') + call check_info(info, 'orthogonalize_against_basis', module=this_module, & + & procedure='DGS_vector_against_basis_rsp, pass 2') ! combine passes proj_coefficients = proj_coefficients + wrk @@ -391,9 +452,7 @@ subroutine initialize_krylov_subspace_rdp(X, X0) ! Internal variables. class(abstract_vector_rdp), allocatable :: Q(:) - real(dp), allocatable :: R(:, :) - integer, allocatable :: perm(:) - integer :: i, p, info + integer :: p ! Zero-out X. call zero_basis(X) @@ -401,28 +460,34 @@ subroutine initialize_krylov_subspace_rdp(X, X0) ! Deals with optional args. if(present(X0)) then p = size(X0) - ! Initialize. - allocate(Q(1:p), source=X0(1:p)) - do i = 1, p - call Q(i)%zero() - call Q(i)%add(X0(i)) - enddo - + allocate(Q(p), source=X0) ! Orthonormalize. - allocate(R(p, p)) ; R = zero_rdp - allocate(perm(p)) ; perm = 0 - call qr(Q, R, perm, info) - do i = 1, p - call X(i)%add(Q(i)) - enddo + call orthonormalize_basis(Q) + call copy_basis(X(:p), Q) endif return end subroutine initialize_krylov_subspace_rdp + subroutine orthonormalize_basis_rdp(X) + !! Orthonormalizes the `abstract_vector` basis `X` + class(abstract_vector_rdp), intent(inout) :: X(:) + !! Input `abstract_vector` basis to orthogonalize against + + ! internals + real(dp) :: R(size(X),size(X)) + integer :: info + + ! internals + call qr(X, R, info) + call check_info(info, 'qr', module=this_module, procedure='orthonormalize_basis_rdp') + + return + end subroutine orthonormalize_basis_rdp + subroutine orthogonalize_vector_against_basis_rdp(y, X, info, if_chk_orthonormal, beta) - !! Orthonormalizes the `abstract_vector` `y` against a basis `X` of `abstract_vector`. + !! Orthogonalizes the `abstract_vector` `y` against a basis `X` of `abstract_vector`. class(abstract_vector_rdp), intent(inout) :: y !! Input `abstract_vector` to orthogonalize class(abstract_vector_rdp), intent(in) :: X(:) @@ -479,7 +544,7 @@ subroutine orthogonalize_vector_against_basis_rdp(y, X, info, if_chk_orthonormal end subroutine orthogonalize_vector_against_basis_rdp subroutine orthogonalize_basis_against_basis_rdp(Y, X, info, if_chk_orthonormal, beta) - !! Orthonormalizes the `abstract_vector` basis `Y` against a basis `X` of `abstract_vector`. + !! Orthogonalizes the `abstract_vector` basis `Y` against a basis `X` of `abstract_vector`. class(abstract_vector_rdp), intent(inout) :: Y(:) !! Input `abstract_vector` basis to orthogonalize class(abstract_vector_rdp), intent(in) :: X(:) @@ -566,11 +631,12 @@ subroutine DGS_vector_against_basis_rdp(y, X, info, if_chk_orthonormal, beta) ! Orthogonalize vector y w.r.t. to Krylov basis X in two passes of GS. ! first pass call orthogonalize_against_basis(y, X, info, if_chk_orthonormal=.false., beta=proj_coefficients) - call check_info(info, 'orthogonalize_against_basis', module=this_module, procedure='DGS_vector_against_basis, first pass') + call check_info(info, 'orthogonalize_against_basis', module=this_module, & + & procedure='DGS_vector_against_basis_rdp, pass 1') ! second pass call orthogonalize_against_basis(y, X, info, if_chk_orthonormal=.false., beta=wrk) - call check_info(info, 'orthogonalize_against_basis', module=this_module, procedure='DGS_vector_against_basis_rdp, second& - & pass') + call check_info(info, 'orthogonalize_against_basis', module=this_module, & + & procedure='DGS_vector_against_basis_rdp, pass 2') ! combine passes proj_coefficients = proj_coefficients + wrk @@ -631,9 +697,7 @@ subroutine initialize_krylov_subspace_csp(X, X0) ! Internal variables. class(abstract_vector_csp), allocatable :: Q(:) - complex(sp), allocatable :: R(:, :) - integer, allocatable :: perm(:) - integer :: i, p, info + integer :: p ! Zero-out X. call zero_basis(X) @@ -641,28 +705,34 @@ subroutine initialize_krylov_subspace_csp(X, X0) ! Deals with optional args. if(present(X0)) then p = size(X0) - ! Initialize. - allocate(Q(1:p), source=X0(1:p)) - do i = 1, p - call Q(i)%zero() - call Q(i)%add(X0(i)) - enddo - + allocate(Q(p), source=X0) ! Orthonormalize. - allocate(R(p, p)) ; R = zero_rsp - allocate(perm(p)) ; perm = 0 - call qr(Q, R, perm, info) - do i = 1, p - call X(i)%add(Q(i)) - enddo + call orthonormalize_basis(Q) + call copy_basis(X(:p), Q) endif return end subroutine initialize_krylov_subspace_csp + subroutine orthonormalize_basis_csp(X) + !! Orthonormalizes the `abstract_vector` basis `X` + class(abstract_vector_csp), intent(inout) :: X(:) + !! Input `abstract_vector` basis to orthogonalize against + + ! internals + complex(sp) :: R(size(X),size(X)) + integer :: info + + ! internals + call qr(X, R, info) + call check_info(info, 'qr', module=this_module, procedure='orthonormalize_basis_csp') + + return + end subroutine orthonormalize_basis_csp + subroutine orthogonalize_vector_against_basis_csp(y, X, info, if_chk_orthonormal, beta) - !! Orthonormalizes the `abstract_vector` `y` against a basis `X` of `abstract_vector`. + !! Orthogonalizes the `abstract_vector` `y` against a basis `X` of `abstract_vector`. class(abstract_vector_csp), intent(inout) :: y !! Input `abstract_vector` to orthogonalize class(abstract_vector_csp), intent(in) :: X(:) @@ -719,7 +789,7 @@ subroutine orthogonalize_vector_against_basis_csp(y, X, info, if_chk_orthonormal end subroutine orthogonalize_vector_against_basis_csp subroutine orthogonalize_basis_against_basis_csp(Y, X, info, if_chk_orthonormal, beta) - !! Orthonormalizes the `abstract_vector` basis `Y` against a basis `X` of `abstract_vector`. + !! Orthogonalizes the `abstract_vector` basis `Y` against a basis `X` of `abstract_vector`. class(abstract_vector_csp), intent(inout) :: Y(:) !! Input `abstract_vector` basis to orthogonalize class(abstract_vector_csp), intent(in) :: X(:) @@ -806,11 +876,12 @@ subroutine DGS_vector_against_basis_csp(y, X, info, if_chk_orthonormal, beta) ! Orthogonalize vector y w.r.t. to Krylov basis X in two passes of GS. ! first pass call orthogonalize_against_basis(y, X, info, if_chk_orthonormal=.false., beta=proj_coefficients) - call check_info(info, 'orthogonalize_against_basis', module=this_module, procedure='DGS_vector_against_basis, first pass') + call check_info(info, 'orthogonalize_against_basis', module=this_module, & + & procedure='DGS_vector_against_basis_csp, pass 1') ! second pass call orthogonalize_against_basis(y, X, info, if_chk_orthonormal=.false., beta=wrk) - call check_info(info, 'orthogonalize_against_basis', module=this_module, procedure='DGS_vector_against_basis_csp, second& - & pass') + call check_info(info, 'orthogonalize_against_basis', module=this_module, & + & procedure='DGS_vector_against_basis_csp, pass 2') ! combine passes proj_coefficients = proj_coefficients + wrk @@ -871,9 +942,7 @@ subroutine initialize_krylov_subspace_cdp(X, X0) ! Internal variables. class(abstract_vector_cdp), allocatable :: Q(:) - complex(dp), allocatable :: R(:, :) - integer, allocatable :: perm(:) - integer :: i, p, info + integer :: p ! Zero-out X. call zero_basis(X) @@ -881,28 +950,34 @@ subroutine initialize_krylov_subspace_cdp(X, X0) ! Deals with optional args. if(present(X0)) then p = size(X0) - ! Initialize. - allocate(Q(1:p), source=X0(1:p)) - do i = 1, p - call Q(i)%zero() - call Q(i)%add(X0(i)) - enddo - + allocate(Q(p), source=X0) ! Orthonormalize. - allocate(R(p, p)) ; R = zero_rdp - allocate(perm(p)) ; perm = 0 - call qr(Q, R, perm, info) - do i = 1, p - call X(i)%add(Q(i)) - enddo + call orthonormalize_basis(Q) + call copy_basis(X(:p), Q) endif return end subroutine initialize_krylov_subspace_cdp + subroutine orthonormalize_basis_cdp(X) + !! Orthonormalizes the `abstract_vector` basis `X` + class(abstract_vector_cdp), intent(inout) :: X(:) + !! Input `abstract_vector` basis to orthogonalize against + + ! internals + complex(dp) :: R(size(X),size(X)) + integer :: info + + ! internals + call qr(X, R, info) + call check_info(info, 'qr', module=this_module, procedure='orthonormalize_basis_cdp') + + return + end subroutine orthonormalize_basis_cdp + subroutine orthogonalize_vector_against_basis_cdp(y, X, info, if_chk_orthonormal, beta) - !! Orthonormalizes the `abstract_vector` `y` against a basis `X` of `abstract_vector`. + !! Orthogonalizes the `abstract_vector` `y` against a basis `X` of `abstract_vector`. class(abstract_vector_cdp), intent(inout) :: y !! Input `abstract_vector` to orthogonalize class(abstract_vector_cdp), intent(in) :: X(:) @@ -959,7 +1034,7 @@ subroutine orthogonalize_vector_against_basis_cdp(y, X, info, if_chk_orthonormal end subroutine orthogonalize_vector_against_basis_cdp subroutine orthogonalize_basis_against_basis_cdp(Y, X, info, if_chk_orthonormal, beta) - !! Orthonormalizes the `abstract_vector` basis `Y` against a basis `X` of `abstract_vector`. + !! Orthogonalizes the `abstract_vector` basis `Y` against a basis `X` of `abstract_vector`. class(abstract_vector_cdp), intent(inout) :: Y(:) !! Input `abstract_vector` basis to orthogonalize class(abstract_vector_cdp), intent(in) :: X(:) @@ -1046,11 +1121,12 @@ subroutine DGS_vector_against_basis_cdp(y, X, info, if_chk_orthonormal, beta) ! Orthogonalize vector y w.r.t. to Krylov basis X in two passes of GS. ! first pass call orthogonalize_against_basis(y, X, info, if_chk_orthonormal=.false., beta=proj_coefficients) - call check_info(info, 'orthogonalize_against_basis', module=this_module, procedure='DGS_vector_against_basis, first pass') + call check_info(info, 'orthogonalize_against_basis', module=this_module, & + & procedure='DGS_vector_against_basis_cdp, pass 1') ! second pass call orthogonalize_against_basis(y, X, info, if_chk_orthonormal=.false., beta=wrk) - call check_info(info, 'orthogonalize_against_basis', module=this_module, procedure='DGS_vector_against_basis_cdp, second& - & pass') + call check_info(info, 'orthogonalize_against_basis', module=this_module, & + & procedure='DGS_vector_against_basis_cdp, pass 2') ! combine passes proj_coefficients = proj_coefficients + wrk diff --git a/src/BaseKrylov.fypp b/src/BaseKrylov.fypp index 12e2b3b..9d520ac 100644 --- a/src/BaseKrylov.fypp +++ b/src/BaseKrylov.fypp @@ -4,7 +4,7 @@ module lightkrylov_BaseKrylov use iso_fortran_env use stdlib_optval, only: optval use stdlib_linalg, only: eye - use LightKrylov_constants + use LightKrylov_Constants use LightKrylov_Logger use LightKrylov_Utils use LightKrylov_AbstractVectors @@ -13,12 +13,13 @@ module lightkrylov_BaseKrylov private character(len=128), parameter :: this_module = 'LightKrylov_BaseKrylov' - public :: qr public :: apply_permutation_matrix public :: apply_inverse_permutation_matrix public :: arnoldi public :: initialize_krylov_subspace + public :: is_orthonormal + public :: orthonormalize_basis public :: orthogonalize_against_basis public :: lanczos_bidiagonalization public :: lanczos_tridiagonalization @@ -67,6 +68,18 @@ module lightkrylov_BaseKrylov #:endfor end interface + interface is_orthonormal + #:for kind, type in RC_KINDS_TYPES + module procedure is_orthonormal_${type[0]}$${kind}$ + #:endfor + end interface + + interface orthonormalize_basis + #:for kind, type in RC_KINDS_TYPES + module procedure orthonormalize_basis_${type[0]}$${kind}$ + #:endfor + end interface + interface orthogonalize_against_basis #:for kind, type in RC_KINDS_TYPES module procedure orthogonalize_vector_against_basis_${type[0]}$${kind}$ @@ -119,6 +132,19 @@ contains !----- UTILITY FUNCTIONS ----- !------------------------------------- + #:for kind, type in RC_KINDS_TYPES + logical function is_orthonormal_${type[0]}$${kind}$(X) result(ortho) + class(abstract_vector_${type[0]}$${kind}$), intent(in) :: X(:) + ${type}$, dimension(size(X), size(X)) :: G + ortho = .true. + call innerprod(G, X, X) + if (norm2(abs(G - eye(size(X)))) > rtol_sp) then + ! The basis is not orthonormal. Cannot orthonormalize. + ortho = .false. + end if + end function + #:endfor + #:for kind, type in RC_KINDS_TYPES subroutine initialize_krylov_subspace_${type[0]}$${kind}$(X, X0) class(abstract_vector_${type[0]}$${kind}$), intent(inout) :: X(:) @@ -126,9 +152,7 @@ contains ! Internal variables. class(abstract_vector_${type[0]}$${kind}$), allocatable :: Q(:) - ${type}$, allocatable :: R(:, :) - integer, allocatable :: perm(:) - integer :: i, p, info + integer :: p ! Zero-out X. call zero_basis(X) @@ -136,28 +160,34 @@ contains ! Deals with optional args. if(present(X0)) then p = size(X0) - ! Initialize. - allocate(Q(1:p), source=X0(1:p)) - do i = 1, p - call Q(i)%zero() - call Q(i)%add(X0(i)) - enddo - + allocate(Q(p), source=X0) ! Orthonormalize. - allocate(R(p, p)) ; R = zero_r${kind}$ - allocate(perm(p)) ; perm = 0 - call qr(Q, R, perm, info) - do i = 1, p - call X(i)%add(Q(i)) - enddo + call orthonormalize_basis(Q) + call copy_basis(X(:p), Q) endif return end subroutine initialize_krylov_subspace_${type[0]}$${kind}$ + subroutine orthonormalize_basis_${type[0]}$${kind}$(X) + !! Orthonormalizes the `abstract_vector` basis `X` + class(abstract_vector_${type[0]}$${kind}$), intent(inout) :: X(:) + !! Input `abstract_vector` basis to orthogonalize against + + ! internals + ${type}$ :: R(size(X),size(X)) + integer :: info + + ! internals + call qr(X, R, info) + call check_info(info, 'qr', module=this_module, procedure='orthonormalize_basis_${type[0]}$${kind}$') + + return + end subroutine orthonormalize_basis_${type[0]}$${kind}$ + subroutine orthogonalize_vector_against_basis_${type[0]}$${kind}$(y, X, info, if_chk_orthonormal, beta) - !! Orthonormalizes the `abstract_vector` `y` against a basis `X` of `abstract_vector`. + !! Orthogonalizes the `abstract_vector` `y` against a basis `X` of `abstract_vector`. class(abstract_vector_${type[0]}$${kind}$), intent(inout) :: y !! Input `abstract_vector` to orthogonalize class(abstract_vector_${type[0]}$${kind}$), intent(in) :: X(:) @@ -214,7 +244,7 @@ contains end subroutine orthogonalize_vector_against_basis_${type[0]}$${kind}$ subroutine orthogonalize_basis_against_basis_${type[0]}$${kind}$(Y, X, info, if_chk_orthonormal, beta) - !! Orthonormalizes the `abstract_vector` basis `Y` against a basis `X` of `abstract_vector`. + !! Orthogonalizes the `abstract_vector` basis `Y` against a basis `X` of `abstract_vector`. class(abstract_vector_${type[0]}$${kind}$), intent(inout) :: Y(:) !! Input `abstract_vector` basis to orthogonalize class(abstract_vector_${type[0]}$${kind}$), intent(in) :: X(:) @@ -301,10 +331,12 @@ contains ! Orthogonalize vector y w.r.t. to Krylov basis X in two passes of GS. ! first pass call orthogonalize_against_basis(y, X, info, if_chk_orthonormal=.false., beta=proj_coefficients) - call check_info(info, 'orthogonalize_against_basis', module=this_module, procedure='DGS_vector_against_basis, first pass') + call check_info(info, 'orthogonalize_against_basis', module=this_module, & + & procedure='DGS_vector_against_basis_${type[0]}$${kind}$, pass 1') ! second pass call orthogonalize_against_basis(y, X, info, if_chk_orthonormal=.false., beta=wrk) - call check_info(info, 'orthogonalize_against_basis', module=this_module, procedure='DGS_vector_against_basis_${type[0]}$${kind}$, second pass') + call check_info(info, 'orthogonalize_against_basis', module=this_module, & + & procedure='DGS_vector_against_basis_${type[0]}$${kind}$, pass 2') ! combine passes proj_coefficients = proj_coefficients + wrk diff --git a/src/Constants.f90 b/src/Constants.f90 index b722bc4..b6d5da9 100644 --- a/src/Constants.f90 +++ b/src/Constants.f90 @@ -1,30 +1,77 @@ -module lightkrylov_constants - implicit none - private - - integer , parameter, public :: sp = selected_real_kind(6, 37) - !! Definition of the single precision data type. - real(sp), parameter, public :: atol_sp = 10.0_sp ** -precision(1.0_sp) - !! Definition of the absolute tolerance for single precision computations. - real(sp), parameter, public :: rtol_sp = sqrt(atol_sp) - !! Definition of the relative tolerance for single precision computations. - - integer , parameter, public :: dp = selected_real_kind(15, 307) - !! Definition of the double precision data type. - real(dp), parameter, public :: atol_dp = 10.0_dp ** -precision(1.0_dp) - !! Definition of the absolute tolerance for double precision computations. - real(dp), parameter, public :: rtol_dp = sqrt(atol_dp) - !! Definition of the relative tolerance for double precision computations. - - real(sp), parameter, public :: one_rsp = 1.0_sp - real(sp), parameter, public :: zero_rsp = 0.0_sp - real(dp), parameter, public :: one_rdp = 1.0_dp - real(dp), parameter, public :: zero_rdp = 0.0_dp - complex(sp), parameter, public :: one_csp = cmplx(1.0_sp, 0.0_sp, kind=sp) - complex(sp), parameter, public :: one_im_csp = cmplx(0.0_sp, 1.0_sp, kind=sp) - complex(sp), parameter, public :: zero_csp = cmplx(0.0_sp, 0.0_sp, kind=sp) - complex(dp), parameter, public :: one_cdp = cmplx(1.0_dp, 0.0_dp, kind=dp) - complex(sp), parameter, public :: one_im_cdp = cmplx(0.0_dp, 1.0_dp, kind=dp) - complex(dp), parameter, public :: zero_cdp = cmplx(0.0_dp, 0.0_dp, kind=dp) - -end module lightkrylov_constants +module LightKrylov_Constants + implicit none + private + + integer, private :: nid = 0 + integer, private :: comm_size = 1 + integer, private :: nio = 0 + + integer , parameter, public :: sp = selected_real_kind(6, 37) + !! Definition of the single precision data type. + real(sp), parameter, public :: atol_sp = 10.0_sp ** -precision(1.0_sp) + !! Definition of the absolute tolerance for single precision computations. + real(sp), parameter, public :: rtol_sp = sqrt(atol_sp) + !! Definition of the relative tolerance for single precision computations. + + integer , parameter, public :: dp = selected_real_kind(15, 307) + !! Definition of the double precision data type. + real(dp), parameter, public :: atol_dp = 10.0_dp ** -precision(1.0_dp) + !! Definition of the absolute tolerance for double precision computations. + real(dp), parameter, public :: rtol_dp = sqrt(atol_dp) + !! Definition of the relative tolerance for double precision computations. + + real(sp), parameter, public :: one_rsp = 1.0_sp + real(sp), parameter, public :: zero_rsp = 0.0_sp + real(dp), parameter, public :: one_rdp = 1.0_dp + real(dp), parameter, public :: zero_rdp = 0.0_dp + complex(sp), parameter, public :: one_csp = cmplx(1.0_sp, 0.0_sp, kind=sp) + complex(sp), parameter, public :: one_im_csp = cmplx(0.0_sp, 1.0_sp, kind=sp) + complex(sp), parameter, public :: zero_csp = cmplx(0.0_sp, 0.0_sp, kind=sp) + complex(dp), parameter, public :: one_cdp = cmplx(1.0_dp, 0.0_dp, kind=dp) + complex(sp), parameter, public :: one_im_cdp = cmplx(0.0_dp, 1.0_dp, kind=dp) + complex(dp), parameter, public :: zero_cdp = cmplx(0.0_dp, 0.0_dp, kind=dp) + + ! Getter/setter routines + public :: set_comm_size + public :: set_rank + public :: set_io_rank + public :: get_rank + public :: get_comm_size + public :: io_rank + +contains + + subroutine set_rank(rank) + integer :: rank + nid = rank + end subroutine set_rank + + subroutine set_comm_size(c_size) + integer :: c_size + comm_size = c_size + end subroutine set_comm_size + + subroutine set_io_rank(rk) + integer, intent(in) :: rk + if (rk > comm_size .or. rk < 0) then + if (io_rank()) print *, 'Invalid I/O rank specified!', rk + else + nio = rk + if (io_rank()) print *, 'I/O rank --> rank ', nio + end if + end + + integer function get_rank() result(rank) + rank = nid + end function get_rank + + integer function get_comm_size() result(c_size) + c_size = comm_size + end function get_comm_size + + logical function io_rank() result(is_io) + is_io = .false. + if (nid == nio) is_io = .true. + end function io_rank + +end module LightKrylov_Constants diff --git a/src/ExpmLib.f90 b/src/ExpmLib.f90 index 84258e9..83fd583 100644 --- a/src/ExpmLib.f90 +++ b/src/ExpmLib.f90 @@ -8,17 +8,16 @@ module lightkrylov_expmlib use stdlib_linalg, only: eye ! LightKrylov. - use lightkrylov_constants + use LightKrylov_Constants use LightKrylov_Logger - use lightkrylov_utils - use lightkrylov_AbstractVectors - use lightkrylov_AbstractLinops - use lightkrylov_BaseKrylov + use LightKrylov_Utils + use LightKrylov_AbstractVectors + use LightKrylov_AbstractLinops + use LightKrylov_BaseKrylov implicit none character(len=128), parameter, private :: this_module = 'LightKrylov_ExpmLib' - public :: abstract_exptA_rsp public :: abstract_exptA_rdp public :: abstract_exptA_csp @@ -163,10 +162,10 @@ subroutine expm_rsp(E, A, order) allocate(wrk(n)) ! Compute the L-infinity norm. - a_norm = norml_rsp(A) + a_norm = norml(A) ! Determine scaling factor for the matrix. - ee = int(log2_rsp(a_norm)) + 1 + ee = int(log2(a_norm)) + 1 s = max(0, ee+1) ! Scale the input matrix & initialize polynomial. @@ -531,10 +530,10 @@ subroutine expm_rdp(E, A, order) allocate(wrk(n)) ! Compute the L-infinity norm. - a_norm = norml_rdp(A) + a_norm = norml(A) ! Determine scaling factor for the matrix. - ee = int(log2_rdp(a_norm)) + 1 + ee = int(log2(a_norm)) + 1 s = max(0, ee+1) ! Scale the input matrix & initialize polynomial. @@ -899,10 +898,10 @@ subroutine expm_csp(E, A, order) allocate(wrk(n)) ! Compute the L-infinity norm. - a_norm = norml_csp(A) + a_norm = norml(A) ! Determine scaling factor for the matrix. - ee = int(log2_rsp(a_norm)) + 1 + ee = int(log2(a_norm)) + 1 s = max(0, ee+1) ! Scale the input matrix & initialize polynomial. @@ -1267,10 +1266,10 @@ subroutine expm_cdp(E, A, order) allocate(wrk(n)) ! Compute the L-infinity norm. - a_norm = norml_cdp(A) + a_norm = norml(A) ! Determine scaling factor for the matrix. - ee = int(log2_rdp(a_norm)) + 1 + ee = int(log2(a_norm)) + 1 s = max(0, ee+1) ! Scale the input matrix & initialize polynomial. diff --git a/src/ExpmLib.fypp b/src/ExpmLib.fypp index 234c819..99fc241 100644 --- a/src/ExpmLib.fypp +++ b/src/ExpmLib.fypp @@ -10,17 +10,16 @@ module lightkrylov_expmlib use stdlib_linalg, only: eye ! LightKrylov. - use lightkrylov_constants + use LightKrylov_Constants use LightKrylov_Logger - use lightkrylov_utils - use lightkrylov_AbstractVectors - use lightkrylov_AbstractLinops - use lightkrylov_BaseKrylov + use LightKrylov_Utils + use LightKrylov_AbstractVectors + use LightKrylov_AbstractLinops + use LightKrylov_BaseKrylov implicit none character(len=128), parameter, private :: this_module = 'LightKrylov_ExpmLib' - #:for kind, type in RC_KINDS_TYPES public :: abstract_exptA_${type[0]}$${kind}$ #:endfor @@ -104,10 +103,10 @@ contains allocate(wrk(n)) ! Compute the L-infinity norm. - a_norm = norml_${type[0]}$${kind}$(A) + a_norm = norml(A) ! Determine scaling factor for the matrix. - ee = int(log2_r${kind}$(a_norm)) + 1 + ee = int(log2(a_norm)) + 1 s = max(0, ee+1) ! Scale the input matrix & initialize polynomial. diff --git a/src/IterativeSolvers.f90 b/src/IterativeSolvers.f90 index 192ec8a..eb5c14b 100644 --- a/src/IterativeSolvers.f90 +++ b/src/IterativeSolvers.f90 @@ -8,12 +8,12 @@ module lightkrylov_IterativeSolvers use stdlib_linalg, only: lstsq, svd use stdlib_stats, only: median - use lightkrylov_constants + use LightKrylov_Constants Use LightKrylov_Logger - use lightkrylov_Utils - use lightkrylov_AbstractVectors - use lightkrylov_AbstractLinops - use lightkrylov_BaseKrylov + use LightKrylov_Utils + use LightKrylov_AbstractVectors + use LightKrylov_AbstractLinops + use LightKrylov_BaseKrylov implicit none private @@ -363,6 +363,8 @@ subroutine eigs_rsp(A, X, eigvals, residuals, info, kdim, select, tolerance, ver integer :: indices(kdim_) real(sp) :: abs_eigvals(kdim_) + ! Re-compute eigenvalues and eigenvectors. + k = min(k, kdim_) ; call eig(H(:k, :k), eigvecs_wrk(:k, :k), eigvals_wrk(:k)) ! Sort eigenvalues. abs_eigvals = abs(eigvals_wrk) ; call sort_index(abs_eigvals, indices, reverse=.true.) eigvals_wrk = eigvals_wrk(indices) ; eigvecs_wrk = eigvecs_wrk(:, indices) @@ -373,7 +375,6 @@ subroutine eigs_rsp(A, X, eigvals, residuals, info, kdim, select, tolerance, ver end block ! Construct eigenvectors. - k = min(k, kdim_) do i = 1, nev call X(i)%zero() do j = 1, k @@ -493,6 +494,8 @@ subroutine eigs_rdp(A, X, eigvals, residuals, info, kdim, select, tolerance, ver integer :: indices(kdim_) real(dp) :: abs_eigvals(kdim_) + ! Re-compute eigenvalues and eigenvectors. + k = min(k, kdim_) ; call eig(H(:k, :k), eigvecs_wrk(:k, :k), eigvals_wrk(:k)) ! Sort eigenvalues. abs_eigvals = abs(eigvals_wrk) ; call sort_index(abs_eigvals, indices, reverse=.true.) eigvals_wrk = eigvals_wrk(indices) ; eigvecs_wrk = eigvecs_wrk(:, indices) @@ -503,7 +506,6 @@ subroutine eigs_rdp(A, X, eigvals, residuals, info, kdim, select, tolerance, ver end block ! Construct eigenvectors. - k = min(k, kdim_) do i = 1, nev call X(i)%zero() do j = 1, k @@ -613,6 +615,8 @@ subroutine eigs_csp(A, X, eigvals, residuals, info, kdim, select, tolerance, ver integer :: indices(kdim_) real(sp) :: abs_eigvals(kdim_) + ! Re-compute eigenvalues and eigenvectors. + k = min(k, kdim_) ; call eig(H(:k, :k), eigvecs_wrk(:k, :k), eigvals_wrk(:k)) ! Sort eigenvalues. abs_eigvals = abs(eigvals_wrk) ; call sort_index(abs_eigvals, indices, reverse=.true.) eigvals_wrk = eigvals_wrk(indices) ; eigvecs_wrk = eigvecs_wrk(:, indices) @@ -623,7 +627,6 @@ subroutine eigs_csp(A, X, eigvals, residuals, info, kdim, select, tolerance, ver end block ! Construct eigenvectors. - k = min(k, kdim_) do i = 1, nev call X(i)%zero() do j = 1, k @@ -733,6 +736,8 @@ subroutine eigs_cdp(A, X, eigvals, residuals, info, kdim, select, tolerance, ver integer :: indices(kdim_) real(dp) :: abs_eigvals(kdim_) + ! Re-compute eigenvalues and eigenvectors. + k = min(k, kdim_) ; call eig(H(:k, :k), eigvecs_wrk(:k, :k), eigvals_wrk(:k)) ! Sort eigenvalues. abs_eigvals = abs(eigvals_wrk) ; call sort_index(abs_eigvals, indices, reverse=.true.) eigvals_wrk = eigvals_wrk(indices) ; eigvecs_wrk = eigvecs_wrk(:, indices) @@ -743,7 +748,6 @@ subroutine eigs_cdp(A, X, eigvals, residuals, info, kdim, select, tolerance, ver end block ! Construct eigenvectors. - k = min(k, kdim_) do i = 1, nev call X(i)%zero() do j = 1, k @@ -1242,7 +1246,7 @@ subroutine svds_rsp(A, U, S, V, residuals, info, kdim, tolerance, verbosity) S = svdvals_wrk(:nsv) ; residuals = residuals_wrk(:nsv) ! Singular vectors. - k = min(k, kdim_) + k = min(k, kdim_) ; info = k do i = 1, nsv call U(i)%zero() ; call V(i)%zero() do j = 1, k @@ -1336,7 +1340,7 @@ subroutine svds_rdp(A, U, S, V, residuals, info, kdim, tolerance, verbosity) S = svdvals_wrk(:nsv) ; residuals = residuals_wrk(:nsv) ! Singular vectors. - k = min(k, kdim_) + k = min(k, kdim_) ; info = k do i = 1, nsv call U(i)%zero() ; call V(i)%zero() do j = 1, k @@ -1430,7 +1434,7 @@ subroutine svds_csp(A, U, S, V, residuals, info, kdim, tolerance, verbosity) S = svdvals_wrk(:nsv) ; residuals = residuals_wrk(:nsv) ! Singular vectors. - k = min(k, kdim_) + k = min(k, kdim_) ; info = k do i = 1, nsv call U(i)%zero() ; call V(i)%zero() do j = 1, k @@ -1524,7 +1528,7 @@ subroutine svds_cdp(A, U, S, V, residuals, info, kdim, tolerance, verbosity) S = svdvals_wrk(:nsv) ; residuals = residuals_wrk(:nsv) ! Singular vectors. - k = min(k, kdim_) + k = min(k, kdim_) ; info = k do i = 1, nsv call U(i)%zero() ; call V(i)%zero() do j = 1, k @@ -1580,7 +1584,7 @@ subroutine gmres_rsp(A, b, x, info, preconditioner, options, transpose) class(abstract_precond_rsp), allocatable :: precond ! Miscellaneous. - integer :: i, k + integer :: i, k, niter real(sp), allocatable :: alpha(:) class(abstract_vector_rsp), allocatable :: dx, wrk @@ -1617,7 +1621,7 @@ subroutine gmres_rsp(A, b, x, info, preconditioner, options, transpose) allocate(alpha(kdim)) ; alpha = 0.0_sp allocate(e(kdim+1)) ; e = 0.0_sp - info = 0 + info = 0 ; niter = 0 ! Initial Krylov vector. if (x%norm() > 0) then @@ -1666,7 +1670,7 @@ subroutine gmres_rsp(A, b, x, info, preconditioner, options, transpose) beta = norm2(abs(e(:k+1) - matmul(H(:k+1, :k), y(:k)))) ! Current number of iterations performed. - info = info + 1 + niter = niter + 1 ! Check convergence. if (abs(beta) <= tol) then @@ -1694,6 +1698,9 @@ subroutine gmres_rsp(A, b, x, info, preconditioner, options, transpose) enddo gmres_iter + ! Returns the number of iterations. + info = niter + return end subroutine gmres_rsp @@ -1736,7 +1743,7 @@ subroutine gmres_rdp(A, b, x, info, preconditioner, options, transpose) class(abstract_precond_rdp), allocatable :: precond ! Miscellaneous. - integer :: i, k + integer :: i, k, niter real(dp), allocatable :: alpha(:) class(abstract_vector_rdp), allocatable :: dx, wrk @@ -1773,7 +1780,7 @@ subroutine gmres_rdp(A, b, x, info, preconditioner, options, transpose) allocate(alpha(kdim)) ; alpha = 0.0_dp allocate(e(kdim+1)) ; e = 0.0_dp - info = 0 + info = 0 ; niter = 0 ! Initial Krylov vector. if (x%norm() > 0) then @@ -1822,7 +1829,7 @@ subroutine gmres_rdp(A, b, x, info, preconditioner, options, transpose) beta = norm2(abs(e(:k+1) - matmul(H(:k+1, :k), y(:k)))) ! Current number of iterations performed. - info = info + 1 + niter = niter + 1 ! Check convergence. if (abs(beta) <= tol) then @@ -1850,6 +1857,9 @@ subroutine gmres_rdp(A, b, x, info, preconditioner, options, transpose) enddo gmres_iter + ! Returns the number of iterations. + info = niter + return end subroutine gmres_rdp @@ -1892,7 +1902,7 @@ subroutine gmres_csp(A, b, x, info, preconditioner, options, transpose) class(abstract_precond_csp), allocatable :: precond ! Miscellaneous. - integer :: i, k + integer :: i, k, niter complex(sp), allocatable :: alpha(:) class(abstract_vector_csp), allocatable :: dx, wrk @@ -1929,7 +1939,7 @@ subroutine gmres_csp(A, b, x, info, preconditioner, options, transpose) allocate(alpha(kdim)) ; alpha = 0.0_sp allocate(e(kdim+1)) ; e = 0.0_sp - info = 0 + info = 0 ; niter = 0 ! Initial Krylov vector. if (x%norm() > 0) then @@ -1978,7 +1988,7 @@ subroutine gmres_csp(A, b, x, info, preconditioner, options, transpose) beta = norm2(abs(e(:k+1) - matmul(H(:k+1, :k), y(:k)))) ! Current number of iterations performed. - info = info + 1 + niter = niter + 1 ! Check convergence. if (abs(beta) <= tol) then @@ -2006,6 +2016,9 @@ subroutine gmres_csp(A, b, x, info, preconditioner, options, transpose) enddo gmres_iter + ! Returns the number of iterations. + info = niter + return end subroutine gmres_csp @@ -2048,7 +2061,7 @@ subroutine gmres_cdp(A, b, x, info, preconditioner, options, transpose) class(abstract_precond_cdp), allocatable :: precond ! Miscellaneous. - integer :: i, k + integer :: i, k, niter complex(dp), allocatable :: alpha(:) class(abstract_vector_cdp), allocatable :: dx, wrk @@ -2085,7 +2098,7 @@ subroutine gmres_cdp(A, b, x, info, preconditioner, options, transpose) allocate(alpha(kdim)) ; alpha = 0.0_dp allocate(e(kdim+1)) ; e = 0.0_dp - info = 0 + info = 0 ; niter = 0 ! Initial Krylov vector. if (x%norm() > 0) then @@ -2134,7 +2147,7 @@ subroutine gmres_cdp(A, b, x, info, preconditioner, options, transpose) beta = norm2(abs(e(:k+1) - matmul(H(:k+1, :k), y(:k)))) ! Current number of iterations performed. - info = info + 1 + niter = niter + 1 ! Check convergence. if (abs(beta) <= tol) then @@ -2162,6 +2175,9 @@ subroutine gmres_cdp(A, b, x, info, preconditioner, options, transpose) enddo gmres_iter + ! Returns the number of iterations. + info = niter + return end subroutine gmres_cdp diff --git a/src/IterativeSolvers.fypp b/src/IterativeSolvers.fypp index 7653c5f..4268aea 100644 --- a/src/IterativeSolvers.fypp +++ b/src/IterativeSolvers.fypp @@ -10,12 +10,12 @@ module lightkrylov_IterativeSolvers use stdlib_linalg, only: lstsq, svd use stdlib_stats, only: median - use lightkrylov_constants + use LightKrylov_Constants Use LightKrylov_Logger - use lightkrylov_Utils - use lightkrylov_AbstractVectors - use lightkrylov_AbstractLinops - use lightkrylov_BaseKrylov + use LightKrylov_Utils + use LightKrylov_AbstractVectors + use LightKrylov_AbstractLinops + use LightKrylov_BaseKrylov implicit none private @@ -264,6 +264,8 @@ contains integer :: indices(kdim_) real(${kind}$) :: abs_eigvals(kdim_) + ! Re-compute eigenvalues and eigenvectors. + k = min(k, kdim_) ; call eig(H(:k, :k), eigvecs_wrk(:k, :k), eigvals_wrk(:k)) ! Sort eigenvalues. abs_eigvals = abs(eigvals_wrk) ; call sort_index(abs_eigvals, indices, reverse=.true.) eigvals_wrk = eigvals_wrk(indices) ; eigvecs_wrk = eigvecs_wrk(:, indices) @@ -274,7 +276,6 @@ contains end block ! Construct eigenvectors. - k = min(k, kdim_) do i = 1, nev call X(i)%zero() do j = 1, k @@ -488,7 +489,7 @@ contains S = svdvals_wrk(:nsv) ; residuals = residuals_wrk(:nsv) ! Singular vectors. - k = min(k, kdim_) + k = min(k, kdim_) ; info = k do i = 1, nsv call U(i)%zero() ; call V(i)%zero() do j = 1, k @@ -546,7 +547,7 @@ contains class(abstract_precond_${type[0]}$${kind}$), allocatable :: precond ! Miscellaneous. - integer :: i, k + integer :: i, k, niter ${type}$, allocatable :: alpha(:) class(abstract_vector_${type[0]}$${kind}$), allocatable :: dx, wrk @@ -583,7 +584,7 @@ contains allocate(alpha(kdim)) ; alpha = 0.0_${kind}$ allocate(e(kdim+1)) ; e = 0.0_${kind}$ - info = 0 + info = 0 ; niter = 0 ! Initial Krylov vector. if (x%norm() > 0) then @@ -632,7 +633,7 @@ contains beta = norm2(abs(e(:k+1) - matmul(H(:k+1, :k), y(:k)))) ! Current number of iterations performed. - info = info + 1 + niter = niter + 1 ! Check convergence. if (abs(beta) <= tol) then @@ -660,6 +661,9 @@ contains enddo gmres_iter + ! Returns the number of iterations. + info = niter + return end subroutine gmres_${type[0]}$${kind}$ diff --git a/src/LightKrylov.f90 b/src/LightKrylov.f90 index b634ebd..c6055b8 100644 --- a/src/LightKrylov.f90 +++ b/src/LightKrylov.f90 @@ -36,7 +36,11 @@ module LightKrylov public :: abstract_vector_rdp public :: abstract_vector_csp public :: abstract_vector_cdp + public :: innerprod + public :: linear_combination + public :: axpby_basis public :: zero_basis + public :: copy_basis ! AbstractLinops exports. public :: abstract_linop @@ -70,6 +74,8 @@ module LightKrylov public :: apply_permutation_matrix, apply_inverse_permutation_matrix public :: arnoldi public :: initialize_krylov_subspace + public :: orthogonalize_against_basis + public :: orthonormalize_basis public :: lanczos_bidiagonalization public :: lanczos_tridiagonalization public :: krylov_schur @@ -110,10 +116,10 @@ subroutine greetings() write (*, *) " |___/ |___/ " write (*, *) - write (*, *) "Developped by: Jean-Christophe Loiseau" - write (*, *) " J. Simon Kern" - write (*, *) " Arts & Métiers Institute of Technology, 2023." - write (*, *) " jean-christophe.loiseau@ensam.eu" + write (*, *) "Developed by: Jean-Christophe Loiseau" + write (*, *) " J. Simon Kern" + write (*, *) " Arts & Métiers Institute of Technology, 2023." + write (*, *) " jean-christophe.loiseau@ensam.eu" write (*, *) write (*, *) "Version -- beta 0.1.0" diff --git a/src/LightKrylov.fypp b/src/LightKrylov.fypp index a1c8f9c..e5968b1 100644 --- a/src/LightKrylov.fypp +++ b/src/LightKrylov.fypp @@ -37,7 +37,11 @@ module LightKrylov #:for kind, type in RC_KINDS_TYPES public :: abstract_vector_${type[0]}$${kind}$ #:endfor + public :: innerprod + public :: linear_combination + public :: axpby_basis public :: zero_basis + public :: copy_basis ! AbstractLinops exports. public :: abstract_linop @@ -59,6 +63,8 @@ module LightKrylov public :: apply_permutation_matrix, apply_inverse_permutation_matrix public :: arnoldi public :: initialize_krylov_subspace + public :: orthogonalize_against_basis + public :: orthonormalize_basis public :: lanczos_bidiagonalization public :: lanczos_tridiagonalization public :: krylov_schur @@ -97,10 +103,10 @@ contains write (*, *) " |___/ |___/ " write (*, *) - write (*, *) "Developped by: Jean-Christophe Loiseau" - write (*, *) " J. Simon Kern" - write (*, *) " Arts & Métiers Institute of Technology, 2023." - write (*, *) " jean-christophe.loiseau@ensam.eu" + write (*, *) "Developed by: Jean-Christophe Loiseau" + write (*, *) " J. Simon Kern" + write (*, *) " Arts & Métiers Institute of Technology, 2023." + write (*, *) " jean-christophe.loiseau@ensam.eu" write (*, *) write (*, *) "Version -- beta 0.1.0" diff --git a/src/Logger.f90 b/src/Logger.f90 index 2419200..60e6ec7 100644 --- a/src/Logger.f90 +++ b/src/Logger.f90 @@ -1,24 +1,135 @@ module LightKrylov_Logger +#ifdef MPI + use mpi_f08 +#endif + ! Fortran Standard Library use stdlib_optval, only : optval - use stdlib_logger use stdlib_logger, only: logger => global_logger use stdlib_ascii, only : to_lower use stdlib_strings, only : chomp, replace_all ! Testdrive use testdrive, only: error_type, test_failed + ! LightKrylov + use LightKrylov_Constants + implicit none private - logical, parameter :: exit_on_error = .true. - logical, parameter :: exit_on_test_error = .true. + character(len=128), parameter, private :: this_module = 'LightKrylov_Logger' + + logical, parameter, private :: exit_on_error = .true. + logical, parameter, private :: exit_on_test_error = .true. public :: stop_error public :: check_info public :: check_test public :: logger + public :: logger_setup + + ! MPI subroutines + public :: comm_setup + public :: comm_close + contains + subroutine logger_setup(logfile, nio, log_level, log_stdout, log_timestamp) + !! Wrapper to set up MPI if needed and initialize log files + character(len=*), optional, intent(in) :: logfile + !! name of the dedicated LightKrylov logfile + integer, optional, intent(in) :: nio + !! I/O rank for logging + integer, optional, intent(in) :: log_level + !! set logging level + !! 0 : all_level + !! 10 : debug_level + !! 20 : information_level + !! 30 : warning_level + !! 40 : error_level + !! 100 : none_level + logical, optional, intent(in) :: log_stdout + !! duplicate log messages to stdout? + logical, optional, intent(in) :: log_timestamp + !! add timestamp to log messages + + ! internals + character(len=:), allocatable :: logfile_ + integer :: nio_ + integer :: log_level_ + logical :: log_stdout_ + logical :: log_timestamp_ + ! misc + integer :: stat, iunit + + logfile_ = optval(logfile, 'lightkrylov.log') + nio_ = optval(nio, 0) + log_level_ = optval(log_level, 20) + log_level_ = max(0, min(log_level_, 100)) + log_stdout_ = optval(log_stdout, .true.) + log_timestamp_ = optval(log_timestamp, .true.) + + ! set log level + call logger%configure(level=log_level_, time_stamp=log_timestamp_) + + ! set up LightKrylov log file + call logger%add_log_file(logfile_, unit=iunit, stat=stat) + if (stat /= 0) call stop_error('Unable to open logfile '//trim(logfile_)//'.', module=this_module, procedure='logger_setup') + + ! Set up comms + call comm_setup() + + ! Set I/O rank + if (nio_/=0) call set_io_rank(nio_) + + ! log to stdout + if (log_stdout_) then + call logger%add_log_unit(6, stat=stat) + if (stat /= 0) call stop_error('Unable to add stdout to logger.', module=this_module, procedure='logger_setup') + end if + return + end subroutine logger_setup + + subroutine comm_setup() + ! internal + integer :: ierr, nid, comm_size + logical :: mpi_is_initialized + character(len=128) :: msg +#ifdef MPI + ! check if MPI has already been initialized and if not, initialize + call MPI_Initialized(mpi_is_initialized, ierr) + if (.not. mpi_is_initialized) then + call logger%log_debug('Set up parallel run with MPI.', module='LightKrylov', procedure='comm_setup') + call MPI_Init(ierr) + if (ierr /= MPI_SUCCESS) call stop_error("Error initializing MPI", module='LightKrylov',procedure='mpi_init') + else + call logger%log_debug('MPI already initialized.', module='LightKrylov', procedure='comm_setup') + end if + call MPI_Comm_rank(MPI_COMM_WORLD, nid, ierr); call set_rank(nid) + call MPI_Comm_size(MPI_COMM_WORLD, comm_size, ierr); call set_comm_size(comm_size) + write(msg, '(A,I4,A,I4)') 'rank', nid, ', comm_size = ', comm_size + call logger%log_debug(trim(msg), module='LightKrylov', procedure='comm_setup') +#else + write(msg, *) 'Setup serial run' + call set_rank(0) + call set_comm_size(1) + call logger%log_debug(trim(msg), module='LightKrylov', procedure='comm_setup') +#endif + return + end subroutine comm_setup + + subroutine comm_close() + integer :: ierr +#ifdef MPI + character(len=128) :: msg + ! Finalize MPI + call MPI_Finalize(ierr) + if (ierr /= MPI_SUCCESS) call stop_error("Error finalizing MPI", module='LightKrylov',procedure='comm_close') +#else + ierr = 0 +#endif + return + end subroutine comm_close + subroutine stop_error(msg, module, procedure) character(len=*), intent(in) :: msg !! The name of the procedure in which the call happens @@ -26,7 +137,7 @@ subroutine stop_error(msg, module, procedure) !! The name of the module in which the call happens character(len=*), optional, intent(in) :: procedure !! The name of the procedure in which the call happens - call check_info(-1, origin='', module=module, procedure=procedure, info_msg=msg) + call check_info(-1, origin="STOP_ERROR", module=module, procedure=procedure, info_msg=msg) return end subroutine stop_error @@ -337,6 +448,12 @@ subroutine check_info(info, origin, module, procedure, info_msg) ierr = -1 end if ! + ! stop error + ! + else if (trim(origin) == 'STOP_ERROR') then + call logger%log_error(trim(origin), module=module, procedure=procedure, stat=info, errmsg=trim(str)) + ierr = -1 + ! ! Default ! else diff --git a/src/TestTypes.f90 b/src/TestUtils.f90 similarity index 53% rename from src/TestTypes.f90 rename to src/TestUtils.f90 index 0ef52a6..a19781c 100644 --- a/src/TestTypes.f90 +++ b/src/TestUtils.f90 @@ -1,18 +1,25 @@ -module LightKrylov_TestTypes - ! Frortran Standard Library +module LightKrylov_TestUtils + ! Fortran Standard Library use stdlib_stats_distribution_normal, only: normal => rvs_normal use stdlib_optval, only: optval + use stdlib_linalg, only: eye ! LightKrylov use LightKrylov + use LightKrylov_Constants implicit none private - character(len=128), parameter, private :: this_module = 'LightKrylov_TestTypes' + character(len=128), parameter, private :: this_module = 'LightKrylov_TestUtils' integer, parameter, public :: test_size = 128 + public :: get_data + public :: put_data + public :: init_rand + public :: get_err_str + !----------------------------------------------- !----- TEST VECTOR TYPE DEFINITION ----- !----------------------------------------------- @@ -21,7 +28,7 @@ module LightKrylov_TestTypes real(sp), dimension(test_size) :: data = 0.0_sp contains private - procedure, pass(self), public :: zero => zero_rsp + procedure, pass(self), public :: zero => init_zero_rsp procedure, pass(self), public :: dot => dot_rsp procedure, pass(self), public :: scal => scal_rsp procedure, pass(self), public :: axpby => axpby_rsp @@ -33,7 +40,7 @@ module LightKrylov_TestTypes real(dp), dimension(test_size) :: data = 0.0_dp contains private - procedure, pass(self), public :: zero => zero_rdp + procedure, pass(self), public :: zero => init_zero_rdp procedure, pass(self), public :: dot => dot_rdp procedure, pass(self), public :: scal => scal_rdp procedure, pass(self), public :: axpby => axpby_rdp @@ -45,7 +52,7 @@ module LightKrylov_TestTypes complex(sp), dimension(test_size) :: data = 0.0_sp contains private - procedure, pass(self), public :: zero => zero_csp + procedure, pass(self), public :: zero => init_zero_csp procedure, pass(self), public :: dot => dot_csp procedure, pass(self), public :: scal => scal_csp procedure, pass(self), public :: axpby => axpby_csp @@ -57,7 +64,7 @@ module LightKrylov_TestTypes complex(dp), dimension(test_size) :: data = 0.0_dp contains private - procedure, pass(self), public :: zero => zero_cdp + procedure, pass(self), public :: zero => init_zero_cdp procedure, pass(self), public :: dot => dot_cdp procedure, pass(self), public :: scal => scal_cdp procedure, pass(self), public :: axpby => axpby_cdp @@ -134,17 +141,72 @@ module LightKrylov_TestTypes procedure, pass(self), public :: rmatvec => hermitian_matvec_cdp end type + + interface get_data + module procedure get_data_vec_rsp + module procedure get_data_vec_basis_rsp + module procedure get_data_linop_rsp + module procedure get_data_vec_rdp + module procedure get_data_vec_basis_rdp + module procedure get_data_linop_rdp + module procedure get_data_vec_csp + module procedure get_data_vec_basis_csp + module procedure get_data_linop_csp + module procedure get_data_vec_cdp + module procedure get_data_vec_basis_cdp + module procedure get_data_linop_cdp + end interface + + interface put_data + module procedure put_data_vec_rsp + module procedure put_data_vec_basis_rsp + module procedure put_data_linop_rsp + module procedure put_data_vec_rdp + module procedure put_data_vec_basis_rdp + module procedure put_data_linop_rdp + module procedure put_data_vec_csp + module procedure put_data_vec_basis_csp + module procedure put_data_linop_csp + module procedure put_data_vec_cdp + module procedure put_data_vec_basis_cdp + module procedure put_data_linop_cdp + end interface + + interface init_rand + module procedure init_rand_vec_rsp + module procedure init_rand_basis_rsp + module procedure init_rand_linop_rsp + module procedure init_rand_spd_linop_rsp + module procedure init_rand_vec_rdp + module procedure init_rand_basis_rdp + module procedure init_rand_linop_rdp + module procedure init_rand_spd_linop_rdp + module procedure init_rand_vec_csp + module procedure init_rand_basis_csp + module procedure init_rand_linop_csp + module procedure init_rand_hermitian_linop_csp + module procedure init_rand_vec_cdp + module procedure init_rand_basis_cdp + module procedure init_rand_linop_cdp + module procedure init_rand_hermitian_linop_cdp + end interface + + interface get_err_str + module procedure get_err_str_sp + module procedure get_err_str_dp + end interface + contains !---------------------------------------------------------- !----- TYPE-BOUND PROCEDURES FOR TEST VECTORS ----- !---------------------------------------------------------- - subroutine zero_rsp(self) + subroutine init_zero_rsp(self) class(vector_rsp), intent(inout) :: self self%data = 0.0_sp return - end subroutine zero_rsp + end subroutine init_zero_rsp function dot_rsp(self, vec) result(alpha) class(vector_rsp), intent(in) :: self @@ -200,11 +262,11 @@ subroutine rand_rsp(self, ifnorm) endif end subroutine rand_rsp - subroutine zero_rdp(self) + subroutine init_zero_rdp(self) class(vector_rdp), intent(inout) :: self self%data = 0.0_dp return - end subroutine zero_rdp + end subroutine init_zero_rdp function dot_rdp(self, vec) result(alpha) class(vector_rdp), intent(in) :: self @@ -260,11 +322,11 @@ subroutine rand_rdp(self, ifnorm) endif end subroutine rand_rdp - subroutine zero_csp(self) + subroutine init_zero_csp(self) class(vector_csp), intent(inout) :: self self%data = 0.0_sp return - end subroutine zero_csp + end subroutine init_zero_csp function dot_csp(self, vec) result(alpha) class(vector_csp), intent(in) :: self @@ -320,11 +382,11 @@ subroutine rand_csp(self, ifnorm) endif end subroutine rand_csp - subroutine zero_cdp(self) + subroutine init_zero_cdp(self) class(vector_cdp), intent(inout) :: self self%data = 0.0_dp return - end subroutine zero_cdp + end subroutine init_zero_cdp function dot_cdp(self, vec) result(alpha) class(vector_cdp), intent(in) :: self @@ -601,5 +663,407 @@ subroutine hermitian_matvec_cdp(self, vec_in, vec_out) return end subroutine hermitian_matvec_cdp + + !---------------------------------------------------- + !----- EXTRACT DATA FROM ABSTRACT TYPES ----- + !---------------------------------------------------- + + subroutine get_data_vec_rsp(vec_out, vec_in) + real(sp), intent(out) :: vec_out(:) + type(vector_rsp), intent(in) :: vec_in + vec_out = vec_in%data + return + end subroutine get_data_vec_rsp + + subroutine get_data_vec_basis_rsp(basis_out, basis_in) + real(sp), intent(out) :: basis_out(:, :) + type(vector_rsp), intent(in) :: basis_in(:) + ! Internal variables. + integer :: k + do k = 1, size(basis_in) + basis_out(:, k) = basis_in(k)%data + enddo + return + end subroutine get_data_vec_basis_rsp + + subroutine get_data_linop_rsp(mat_out, linop_in) + real(sp), intent(out) :: mat_out(:, :) + type(linop_rsp), intent(in) :: linop_in + mat_out = linop_in%data + return + end subroutine get_data_linop_rsp + + subroutine get_data_vec_rdp(vec_out, vec_in) + real(dp), intent(out) :: vec_out(:) + type(vector_rdp), intent(in) :: vec_in + vec_out = vec_in%data + return + end subroutine get_data_vec_rdp + + subroutine get_data_vec_basis_rdp(basis_out, basis_in) + real(dp), intent(out) :: basis_out(:, :) + type(vector_rdp), intent(in) :: basis_in(:) + ! Internal variables. + integer :: k + do k = 1, size(basis_in) + basis_out(:, k) = basis_in(k)%data + enddo + return + end subroutine get_data_vec_basis_rdp + + subroutine get_data_linop_rdp(mat_out, linop_in) + real(dp), intent(out) :: mat_out(:, :) + type(linop_rdp), intent(in) :: linop_in + mat_out = linop_in%data + return + end subroutine get_data_linop_rdp + + subroutine get_data_vec_csp(vec_out, vec_in) + complex(sp), intent(out) :: vec_out(:) + type(vector_csp), intent(in) :: vec_in + vec_out = vec_in%data + return + end subroutine get_data_vec_csp + + subroutine get_data_vec_basis_csp(basis_out, basis_in) + complex(sp), intent(out) :: basis_out(:, :) + type(vector_csp), intent(in) :: basis_in(:) + ! Internal variables. + integer :: k + do k = 1, size(basis_in) + basis_out(:, k) = basis_in(k)%data + enddo + return + end subroutine get_data_vec_basis_csp + + subroutine get_data_linop_csp(mat_out, linop_in) + complex(sp), intent(out) :: mat_out(:, :) + type(linop_csp), intent(in) :: linop_in + mat_out = linop_in%data + return + end subroutine get_data_linop_csp + + subroutine get_data_vec_cdp(vec_out, vec_in) + complex(dp), intent(out) :: vec_out(:) + type(vector_cdp), intent(in) :: vec_in + vec_out = vec_in%data + return + end subroutine get_data_vec_cdp + + subroutine get_data_vec_basis_cdp(basis_out, basis_in) + complex(dp), intent(out) :: basis_out(:, :) + type(vector_cdp), intent(in) :: basis_in(:) + ! Internal variables. + integer :: k + do k = 1, size(basis_in) + basis_out(:, k) = basis_in(k)%data + enddo + return + end subroutine get_data_vec_basis_cdp + + subroutine get_data_linop_cdp(mat_out, linop_in) + complex(dp), intent(out) :: mat_out(:, :) + type(linop_cdp), intent(in) :: linop_in + mat_out = linop_in%data + return + end subroutine get_data_linop_cdp + + + !---------------------------------------------- + !----- PUT DATA TO ABSTRACT TYPES ----- + !---------------------------------------------- + + subroutine put_data_vec_rsp(vec_out, vec_in) + type(vector_rsp), intent(out) :: vec_out + real(sp), intent(in) :: vec_in + vec_out%data = vec_in + return + end subroutine put_data_vec_rsp + + subroutine put_data_vec_basis_rsp(basis_out, basis_in) + type(vector_rsp), intent(out) :: basis_out(:) + real(sp), intent(in) :: basis_in(:, :) + ! Internal variables. + integer :: k + do k = 1, size(basis_out) + basis_out(k)%data = basis_in(:, k) + enddo + return + end subroutine put_data_vec_basis_rsp + + subroutine put_data_linop_rsp(linop_out, mat_in) + type(linop_rsp), intent(out) :: linop_out + real(sp), intent(in) :: mat_in(:, :) + ! Internal variables. + linop_out%data = mat_in + return + end subroutine put_data_linop_rsp + + subroutine put_data_vec_rdp(vec_out, vec_in) + type(vector_rdp), intent(out) :: vec_out + real(dp), intent(in) :: vec_in + vec_out%data = vec_in + return + end subroutine put_data_vec_rdp + + subroutine put_data_vec_basis_rdp(basis_out, basis_in) + type(vector_rdp), intent(out) :: basis_out(:) + real(dp), intent(in) :: basis_in(:, :) + ! Internal variables. + integer :: k + do k = 1, size(basis_out) + basis_out(k)%data = basis_in(:, k) + enddo + return + end subroutine put_data_vec_basis_rdp + + subroutine put_data_linop_rdp(linop_out, mat_in) + type(linop_rdp), intent(out) :: linop_out + real(dp), intent(in) :: mat_in(:, :) + ! Internal variables. + linop_out%data = mat_in + return + end subroutine put_data_linop_rdp + + subroutine put_data_vec_csp(vec_out, vec_in) + type(vector_csp), intent(out) :: vec_out + complex(sp), intent(in) :: vec_in + vec_out%data = vec_in + return + end subroutine put_data_vec_csp + + subroutine put_data_vec_basis_csp(basis_out, basis_in) + type(vector_csp), intent(out) :: basis_out(:) + complex(sp), intent(in) :: basis_in(:, :) + ! Internal variables. + integer :: k + do k = 1, size(basis_out) + basis_out(k)%data = basis_in(:, k) + enddo + return + end subroutine put_data_vec_basis_csp + + subroutine put_data_linop_csp(linop_out, mat_in) + type(linop_csp), intent(out) :: linop_out + complex(sp), intent(in) :: mat_in(:, :) + ! Internal variables. + linop_out%data = mat_in + return + end subroutine put_data_linop_csp + + subroutine put_data_vec_cdp(vec_out, vec_in) + type(vector_cdp), intent(out) :: vec_out + complex(dp), intent(in) :: vec_in + vec_out%data = vec_in + return + end subroutine put_data_vec_cdp + + subroutine put_data_vec_basis_cdp(basis_out, basis_in) + type(vector_cdp), intent(out) :: basis_out(:) + complex(dp), intent(in) :: basis_in(:, :) + ! Internal variables. + integer :: k + do k = 1, size(basis_out) + basis_out(k)%data = basis_in(:, k) + enddo + return + end subroutine put_data_vec_basis_cdp + + subroutine put_data_linop_cdp(linop_out, mat_in) + type(linop_cdp), intent(out) :: linop_out + complex(dp), intent(in) :: mat_in(:, :) + ! Internal variables. + linop_out%data = mat_in + return + end subroutine put_data_linop_cdp + + + !-------------------------------------------------------------- + !----- INITIALIZE ABSTRACT TYPES WITH RANDOM DATA ----- + !-------------------------------------------------------------- + + subroutine init_rand_vec_rsp(x) + type(vector_rsp), intent(inout) :: x + call x%rand() + return + end subroutine init_rand_vec_rsp + + subroutine init_rand_basis_rsp(X) + type(vector_rsp), intent(inout) :: X(:) + integer :: i + do i = 1, size(X) + call X(i)%rand() + enddo + return + end subroutine init_rand_basis_rsp + + subroutine init_rand_linop_rsp(linop) + type(linop_rsp), intent(inout) :: linop + real(sp), allocatable :: mu(:, :), var(:, :) + allocate(mu(test_size, test_size)) ; mu = 0.0_sp + allocate(var(test_size, test_size)) + var = 1.0_sp + linop%data = normal(mu, var) + return + end subroutine init_rand_linop_rsp + + subroutine init_rand_spd_linop_rsp(linop) + type(spd_linop_rsp), intent(inout) :: linop + real(sp), allocatable :: mu(:, :), var(:, :) + real(sp), allocatable :: data(:, :) + allocate(mu(test_size, test_size)) ; mu = zero_rsp + allocate(var(test_size, test_size)) ; var = one_rsp + + data = normal(mu, var) + linop%data = matmul(data, transpose(data))/test_size + 0.01*eye(test_size) + + return + end subroutine init_rand_spd_linop_rsp + + subroutine init_rand_vec_rdp(x) + type(vector_rdp), intent(inout) :: x + call x%rand() + return + end subroutine init_rand_vec_rdp + + subroutine init_rand_basis_rdp(X) + type(vector_rdp), intent(inout) :: X(:) + integer :: i + do i = 1, size(X) + call X(i)%rand() + enddo + return + end subroutine init_rand_basis_rdp + + subroutine init_rand_linop_rdp(linop) + type(linop_rdp), intent(inout) :: linop + real(dp), allocatable :: mu(:, :), var(:, :) + allocate(mu(test_size, test_size)) ; mu = 0.0_dp + allocate(var(test_size, test_size)) + var = 1.0_dp + linop%data = normal(mu, var) + return + end subroutine init_rand_linop_rdp + + subroutine init_rand_spd_linop_rdp(linop) + type(spd_linop_rdp), intent(inout) :: linop + real(dp), allocatable :: mu(:, :), var(:, :) + real(dp), allocatable :: data(:, :) + allocate(mu(test_size, test_size)) ; mu = zero_rdp + allocate(var(test_size, test_size)) ; var = one_rdp + + data = normal(mu, var) + linop%data = matmul(data, transpose(data))/test_size + 0.01*eye(test_size) + + return + end subroutine init_rand_spd_linop_rdp + + subroutine init_rand_vec_csp(x) + type(vector_csp), intent(inout) :: x + call x%rand() + return + end subroutine init_rand_vec_csp + + subroutine init_rand_basis_csp(X) + type(vector_csp), intent(inout) :: X(:) + integer :: i + do i = 1, size(X) + call X(i)%rand() + enddo + return + end subroutine init_rand_basis_csp + + subroutine init_rand_linop_csp(linop) + type(linop_csp), intent(inout) :: linop + complex(sp), allocatable :: mu(:, :), var(:, :) + allocate(mu(test_size, test_size)) ; mu = 0.0_sp + allocate(var(test_size, test_size)) + var = cmplx(1.0_sp, 1.0_sp, kind=sp) + linop%data = normal(mu, var) + return + end subroutine init_rand_linop_csp + + subroutine init_rand_hermitian_linop_csp(linop) + type(hermitian_linop_csp), intent(inout) :: linop + complex(sp), allocatable :: data(:, :) + complex(sp), allocatable :: mu(:, :), var(:, :) + + allocate(mu(test_size, test_size)) ; mu = 0.0_sp + allocate(var(test_size, test_size)) ; var = cmplx(1.0_sp, 1.0_sp, kind=sp) + + data = normal(mu, var) + data = matmul(data, transpose(conjg(data)))/test_size + 0.01*eye(test_size) + linop%data = data + + return + end subroutine init_rand_hermitian_linop_csp + + subroutine init_rand_vec_cdp(x) + type(vector_cdp), intent(inout) :: x + call x%rand() + return + end subroutine init_rand_vec_cdp + + subroutine init_rand_basis_cdp(X) + type(vector_cdp), intent(inout) :: X(:) + integer :: i + do i = 1, size(X) + call X(i)%rand() + enddo + return + end subroutine init_rand_basis_cdp + + subroutine init_rand_linop_cdp(linop) + type(linop_cdp), intent(inout) :: linop + complex(dp), allocatable :: mu(:, :), var(:, :) + allocate(mu(test_size, test_size)) ; mu = 0.0_dp + allocate(var(test_size, test_size)) + var = cmplx(1.0_dp, 1.0_dp, kind=dp) + linop%data = normal(mu, var) + return + end subroutine init_rand_linop_cdp + + subroutine init_rand_hermitian_linop_cdp(linop) + type(hermitian_linop_cdp), intent(inout) :: linop + complex(dp), allocatable :: data(:, :) + complex(dp), allocatable :: mu(:, :), var(:, :) + + allocate(mu(test_size, test_size)) ; mu = 0.0_dp + allocate(var(test_size, test_size)) ; var = cmplx(1.0_dp, 1.0_dp, kind=dp) + + data = normal(mu, var) + data = matmul(data, transpose(conjg(data)))/test_size + 0.01*eye(test_size) + linop%data = data + + return + end subroutine init_rand_hermitian_linop_cdp + + + subroutine get_err_str_sp(msg, info, err) + character(len=*), intent(inout) :: msg + character(len=*), intent(in) :: info + real(sp) :: err + + ! internals + character*8 :: value_str + character(len=*), parameter :: indent = repeat(" ", 4) + + write(value_str, '(E8.2)') err + msg = indent // info // value_str // achar(10) + + end subroutine get_err_str_sp + subroutine get_err_str_dp(msg, info, err) + character(len=*), intent(inout) :: msg + character(len=*), intent(in) :: info + real(dp) :: err + + ! internals + character*8 :: value_str + character(len=*), parameter :: indent = repeat(" ", 4) + + write(value_str, '(E8.2)') err + msg = indent // info // value_str // achar(10) + + end subroutine get_err_str_dp -end module LightKrylov_TestTypes +end module LightKrylov_TestUtils diff --git a/src/TestTypes.fypp b/src/TestUtils.fypp similarity index 54% rename from src/TestTypes.fypp rename to src/TestUtils.fypp index 227df9b..94661f6 100644 --- a/src/TestTypes.fypp +++ b/src/TestUtils.fypp @@ -1,20 +1,27 @@ #:include "../include/common.fypp" #:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES -module LightKrylov_TestTypes - ! Frortran Standard Library +module LightKrylov_TestUtils + ! Fortran Standard Library use stdlib_stats_distribution_normal, only: normal => rvs_normal use stdlib_optval, only: optval + use stdlib_linalg, only: eye ! LightKrylov use LightKrylov + use LightKrylov_Constants implicit none private - character(len=128), parameter, private :: this_module = 'LightKrylov_TestTypes' + character(len=128), parameter, private :: this_module = 'LightKrylov_TestUtils' integer, parameter, public :: test_size = 128 + public :: get_data + public :: put_data + public :: init_rand + public :: get_err_str + !----------------------------------------------- !----- TEST VECTOR TYPE DEFINITION ----- !----------------------------------------------- @@ -24,7 +31,7 @@ module LightKrylov_TestTypes ${type}$, dimension(test_size) :: data = 0.0_${kind}$ contains private - procedure, pass(self), public :: zero => zero_${type[0]}$${kind}$ + procedure, pass(self), public :: zero => init_zero_${type[0]}$${kind}$ procedure, pass(self), public :: dot => dot_${type[0]}$${kind}$ procedure, pass(self), public :: scal => scal_${type[0]}$${kind}$ procedure, pass(self), public :: axpby => axpby_${type[0]}$${kind}$ @@ -74,6 +81,42 @@ module LightKrylov_TestTypes #:endif #:endfor + + interface get_data + #:for kind, type in RC_KINDS_TYPES + module procedure get_data_vec_${type[0]}$${kind}$ + module procedure get_data_vec_basis_${type[0]}$${kind}$ + module procedure get_data_linop_${type[0]}$${kind}$ + #:endfor + end interface + + interface put_data + #:for kind, type in RC_KINDS_TYPES + module procedure put_data_vec_${type[0]}$${kind}$ + module procedure put_data_vec_basis_${type[0]}$${kind}$ + module procedure put_data_linop_${type[0]}$${kind}$ + #:endfor + end interface + + interface init_rand + #:for kind, type in RC_KINDS_TYPES + module procedure init_rand_vec_${type[0]}$${kind}$ + module procedure init_rand_basis_${type[0]}$${kind}$ + module procedure init_rand_linop_${type[0]}$${kind}$ + #:if type[0] == "r" + module procedure init_rand_spd_linop_${type[0]}$${kind}$ + #:else + module procedure init_rand_hermitian_linop_${type[0]}$${kind}$ + #:endif + #:endfor + end interface + + interface get_err_str + #:for kind, type in REAL_KINDS_TYPES + module procedure get_err_str_${kind}$ + #:endfor + end interface + contains !---------------------------------------------------------- @@ -81,11 +124,11 @@ contains !---------------------------------------------------------- #:for kind, type in RC_KINDS_TYPES - subroutine zero_${type[0]}$${kind}$(self) + subroutine init_zero_${type[0]}$${kind}$(self) class(vector_${type[0]}$${kind}$), intent(inout) :: self self%data = 0.0_${kind}$ return - end subroutine zero_${type[0]}$${kind}$ + end subroutine init_zero_${type[0]}$${kind}$ function dot_${type[0]}$${kind}$(self, vec) result(alpha) class(vector_${type[0]}$${kind}$), intent(in) :: self @@ -231,5 +274,152 @@ contains #:endif #:endfor + + !---------------------------------------------------- + !----- EXTRACT DATA FROM ABSTRACT TYPES ----- + !---------------------------------------------------- + + #:for kind, type in RC_KINDS_TYPES + subroutine get_data_vec_${type[0]}$${kind}$(vec_out, vec_in) + ${type}$, intent(out) :: vec_out(:) + type(vector_${type[0]}$${kind}$), intent(in) :: vec_in + vec_out = vec_in%data + return + end subroutine get_data_vec_${type[0]}$${kind}$ + + subroutine get_data_vec_basis_${type[0]}$${kind}$(basis_out, basis_in) + ${type}$, intent(out) :: basis_out(:, :) + type(vector_${type[0]}$${kind}$), intent(in) :: basis_in(:) + ! Internal variables. + integer :: k + do k = 1, size(basis_in) + basis_out(:, k) = basis_in(k)%data + enddo + return + end subroutine get_data_vec_basis_${type[0]}$${kind}$ + + subroutine get_data_linop_${type[0]}$${kind}$(mat_out, linop_in) + ${type}$, intent(out) :: mat_out(:, :) + type(linop_${type[0]}$${kind}$), intent(in) :: linop_in + mat_out = linop_in%data + return + end subroutine get_data_linop_${type[0]}$${kind}$ + + #:endfor + + !---------------------------------------------- + !----- PUT DATA TO ABSTRACT TYPES ----- + !---------------------------------------------- + + #:for kind, type in RC_KINDS_TYPES + subroutine put_data_vec_${type[0]}$${kind}$(vec_out, vec_in) + type(vector_${type[0]}$${kind}$), intent(out) :: vec_out + ${type}$, intent(in) :: vec_in + vec_out%data = vec_in + return + end subroutine put_data_vec_${type[0]}$${kind}$ + + subroutine put_data_vec_basis_${type[0]}$${kind}$(basis_out, basis_in) + type(vector_${type[0]}$${kind}$), intent(out) :: basis_out(:) + ${type}$, intent(in) :: basis_in(:, :) + ! Internal variables. + integer :: k + do k = 1, size(basis_out) + basis_out(k)%data = basis_in(:, k) + enddo + return + end subroutine put_data_vec_basis_${type[0]}$${kind}$ + + subroutine put_data_linop_${type[0]}$${kind}$(linop_out, mat_in) + type(linop_${type[0]}$${kind}$), intent(out) :: linop_out + ${type}$, intent(in) :: mat_in(:, :) + ! Internal variables. + linop_out%data = mat_in + return + end subroutine put_data_linop_${type[0]}$${kind}$ + + #:endfor + + !-------------------------------------------------------------- + !----- INITIALIZE ABSTRACT TYPES WITH RANDOM DATA ----- + !-------------------------------------------------------------- + + #:for kind, type in RC_KINDS_TYPES + subroutine init_rand_vec_${type[0]}$${kind}$(x) + type(vector_${type[0]}$${kind}$), intent(inout) :: x + call x%rand() + return + end subroutine init_rand_vec_${type[0]}$${kind}$ + + subroutine init_rand_basis_${type[0]}$${kind}$(X) + type(vector_${type[0]}$${kind}$), intent(inout) :: X(:) + integer :: i + do i = 1, size(X) + call X(i)%rand() + enddo + return + end subroutine init_rand_basis_${type[0]}$${kind}$ + + subroutine init_rand_linop_${type[0]}$${kind}$(linop) + type(linop_${type[0]}$${kind}$), intent(inout) :: linop + ${type}$, allocatable :: mu(:, :), var(:, :) + allocate(mu(test_size, test_size)) ; mu = 0.0_${kind}$ + allocate(var(test_size, test_size)) + #:if type[0] == "r" + var = 1.0_${kind}$ + #:else + var = cmplx(1.0_${kind}$, 1.0_${kind}$, kind=${kind}$) + #:endif + linop%data = normal(mu, var) + return + end subroutine init_rand_linop_${type[0]}$${kind}$ + + #:if type[0] == "r" + subroutine init_rand_spd_linop_${type[0]}$${kind}$(linop) + type(spd_linop_${type[0]}$${kind}$), intent(inout) :: linop + ${type}$, allocatable :: mu(:, :), var(:, :) + ${type}$, allocatable :: data(:, :) + allocate(mu(test_size, test_size)) ; mu = zero_${type[0]}$${kind}$ + allocate(var(test_size, test_size)) ; var = one_${type[0]}$${kind}$ + + data = normal(mu, var) + linop%data = matmul(data, transpose(data))/test_size + 0.01*eye(test_size) + + return + end subroutine init_rand_spd_linop_${type[0]}$${kind}$ + #:else + subroutine init_rand_hermitian_linop_${type[0]}$${kind}$(linop) + type(hermitian_linop_${type[0]}$${kind}$), intent(inout) :: linop + ${type}$, allocatable :: data(:, :) + ${type}$, allocatable :: mu(:, :), var(:, :) + + allocate(mu(test_size, test_size)) ; mu = 0.0_${kind}$ + allocate(var(test_size, test_size)) ; var = cmplx(1.0_${kind}$, 1.0_${kind}$, kind=${kind}$) + + data = normal(mu, var) + data = matmul(data, transpose(conjg(data)))/test_size + 0.01*eye(test_size) + linop%data = data + + return + end subroutine init_rand_hermitian_linop_${type[0]}$${kind}$ + #:endif + + #:endfor + + #:for kind, type in REAL_KINDS_TYPES + subroutine get_err_str_${kind}$(msg, info, err) + character(len=*), intent(inout) :: msg + character(len=*), intent(in) :: info + real(${kind}$) :: err + + ! internals + character*8 :: value_str + character(len=*), parameter :: indent = repeat(" ", 4) + + write(value_str, '(E8.2)') err + msg = indent // info // value_str // achar(10) + + end subroutine get_err_str_${kind}$ + #:endfor -end module LightKrylov_TestTypes +end module LightKrylov_TestUtils diff --git a/src/Utils.f90 b/src/Utils.f90 index ea031f9..7c25b84 100644 --- a/src/Utils.f90 +++ b/src/Utils.f90 @@ -3,8 +3,7 @@ module lightkrylov_utils !----- Standard Fortran Library ----- !-------------------------------------------- use iso_fortran_env, only: output_unit - use LightKrylov_Logger - use stdlib_linalg, only: is_hermitian, is_symmetric, diag + use stdlib_linalg, only: is_hermitian, is_symmetric, diag, svd ! Matrix inversion. use stdlib_linalg_lapack, only: getrf, getri ! Eigenvalue problem (general + symmetric). @@ -16,14 +15,15 @@ module lightkrylov_utils !----- LightKrylov ----- !------------------------------- ! Various constants. + use LightKrylov_Logger use LightKrylov_Constants implicit none private - character*128, parameter :: this_module = 'LightKrylov_Utils' + character(len=128), parameter :: this_module = 'LightKrylov_Utils' - public :: assert_shape + public :: assert_shape, norml, log2 ! Compute B = inv(A) in-place for dense matrices. public :: inv ! Compute AX = XD for general dense matrices. @@ -32,18 +32,12 @@ module lightkrylov_utils public :: eigh ! Compute matrix sqrt of input symmetric/hermitian positive definite matrix A public :: sqrtm + public :: sqrtm_eig ! Compute AX = XS where S is in Schur form. public :: schur ! Re-orders the Schur factorization of A. public :: ordschur - public :: log2_rsp - public :: norml_rsp - public :: log2_rdp - public :: norml_rdp - public :: norml_csp - public :: norml_cdp - interface assert_shape module procedure assert_shape_vector_rsp module procedure assert_shape_matrix_rsp @@ -55,6 +49,18 @@ module lightkrylov_utils module procedure assert_shape_matrix_cdp end interface + interface norml + module procedure norml_rsp + module procedure norml_rdp + module procedure norml_csp + module procedure norml_cdp + end interface + + interface log2 + module procedure log2_rsp + module procedure log2_rdp + end interface + interface inv module procedure inv_rsp module procedure inv_rdp @@ -97,6 +103,13 @@ module lightkrylov_utils module procedure sqrtm_cdp end interface + interface sqrtm_eig + module procedure sqrtm_eig_rsp + module procedure sqrtm_eig_rdp + module procedure sqrtm_eig_csp + module procedure sqrtm_eig_cdp + end interface + !------------------------------------------------ !----- OPTS TYPE FOR LINEAR SOLVERS ----- !------------------------------------------------ @@ -176,7 +189,7 @@ subroutine assert_shape_vector_rsp(v, size, routine, matname) !! Name of the asserted vector. ! internals - character*128 :: msg + character(len=256) :: msg if(any(shape(v) /= size)) then write(msg, *) "In routine "//routine//" vector "//matname//" has illegal length ", shape(v), & @@ -198,7 +211,7 @@ subroutine assert_shape_matrix_rsp(A, size, routine, matname) !! Name of the asserted matrix. ! internals - character*128 :: msg + character(len=256) :: msg if(any(shape(A) /= size)) then write(msg, *) "In routine "//routine//" matrix "//matname//" has illegal shape ", shape(A), & @@ -219,7 +232,7 @@ subroutine assert_shape_vector_rdp(v, size, routine, matname) !! Name of the asserted vector. ! internals - character*128 :: msg + character(len=256) :: msg if(any(shape(v) /= size)) then write(msg, *) "In routine "//routine//" vector "//matname//" has illegal length ", shape(v), & @@ -241,7 +254,7 @@ subroutine assert_shape_matrix_rdp(A, size, routine, matname) !! Name of the asserted matrix. ! internals - character*128 :: msg + character(len=256) :: msg if(any(shape(A) /= size)) then write(msg, *) "In routine "//routine//" matrix "//matname//" has illegal shape ", shape(A), & @@ -262,7 +275,7 @@ subroutine assert_shape_vector_csp(v, size, routine, matname) !! Name of the asserted vector. ! internals - character*128 :: msg + character(len=256) :: msg if(any(shape(v) /= size)) then write(msg, *) "In routine "//routine//" vector "//matname//" has illegal length ", shape(v), & @@ -284,7 +297,7 @@ subroutine assert_shape_matrix_csp(A, size, routine, matname) !! Name of the asserted matrix. ! internals - character*128 :: msg + character(len=256) :: msg if(any(shape(A) /= size)) then write(msg, *) "In routine "//routine//" matrix "//matname//" has illegal shape ", shape(A), & @@ -305,7 +318,7 @@ subroutine assert_shape_vector_cdp(v, size, routine, matname) !! Name of the asserted vector. ! internals - character*128 :: msg + character(len=256) :: msg if(any(shape(v) /= size)) then write(msg, *) "In routine "//routine//" vector "//matname//" has illegal length ", shape(v), & @@ -327,7 +340,7 @@ subroutine assert_shape_matrix_cdp(A, size, routine, matname) !! Name of the asserted matrix. ! internals - character*128 :: msg + character(len=256) :: msg if(any(shape(A) /= size)) then write(msg, *) "In routine "//routine//" matrix "//matname//" has illegal shape ", shape(A), & @@ -487,6 +500,54 @@ subroutine ordschur_rsp(T, Q, selected) end subroutine ordschur_rsp subroutine sqrtm_rsp(X, sqrtmX, info) + !! Matrix-valued sqrt function for dense symmetric/hermitian positive (semi-)definite matrices + real(sp), intent(inout) :: X(:,:) + !! Matrix of which to compute the sqrt + real(sp), intent(out) :: sqrtmX(size(X,1),size(X,1)) + !! Return matrix + integer, intent(out) :: info + !! Information flag + + ! internal + real(sp) :: S(size(X,1)) + real(sp) :: U(size(X,1), size(X,1)), VT(size(X,1), size(X,1)) + integer :: i + real(sp) :: symmetry_error + character(len=128) :: msg + + info = 0 + + ! Check if the matrix is symmetric + symmetry_error = 0.5*maxval(X - transpose(X)) + if (symmetry_error > rtol_sp) then + write(msg,*) "Input matrix is not symmetric. 0.5*max(X-X.T) = ", & + & symmetry_error, ", tol = ", rtol_sp + call stop_error(msg, module=this_module, procedure='sqrtm_rsp') + else if (symmetry_error > 10*atol_sp) then + write(msg,*) "Input matrix is not exactly symmetric. 0.5*max(X-X.T) = ", symmetry_error + call logger%log_warning(trim(msg), module=this_module, procedure='sqrtm_rsp') + end if + + ! Perform svd + call svd(X, S, U, VT) + + ! Check if the matrix is positive definite (up to tol) + do i = 1, size(S) + if (S(i) .gt. 10*atol_sp ) then + S(i) = sqrt(S(i)) + else + S(i) = zero_rsp + info = 1 + end if + end do + + ! Reconstruct the square root matrix + sqrtmX = matmul(U, matmul(diag(S), VT)) + + return + end subroutine + + subroutine sqrtm_eig_rsp(X, sqrtmX, info) !! Matrix-valued sqrt function for dense symmetric/hermitian positive (semi-)definite matrices real(sp), intent(in) :: X(:,:) !! Matrix of which to compute the sqrt @@ -499,7 +560,7 @@ subroutine sqrtm_rsp(X, sqrtmX, info) real(sp) :: lambda(size(X,1)) real(sp) :: V(size(X,1), size(X,1)) integer :: i - character*128 :: msg + character(len=128) :: msg info = 0 @@ -532,7 +593,6 @@ subroutine sqrtm_rsp(X, sqrtmX, info) return end subroutine - subroutine inv_rdp(A) !! In-place inversion of A using LAPACK. real(dp), intent(inout) :: A(:, :) @@ -679,6 +739,54 @@ subroutine ordschur_rdp(T, Q, selected) end subroutine ordschur_rdp subroutine sqrtm_rdp(X, sqrtmX, info) + !! Matrix-valued sqrt function for dense symmetric/hermitian positive (semi-)definite matrices + real(dp), intent(inout) :: X(:,:) + !! Matrix of which to compute the sqrt + real(dp), intent(out) :: sqrtmX(size(X,1),size(X,1)) + !! Return matrix + integer, intent(out) :: info + !! Information flag + + ! internal + real(dp) :: S(size(X,1)) + real(dp) :: U(size(X,1), size(X,1)), VT(size(X,1), size(X,1)) + integer :: i + real(dp) :: symmetry_error + character(len=128) :: msg + + info = 0 + + ! Check if the matrix is symmetric + symmetry_error = 0.5*maxval(X - transpose(X)) + if (symmetry_error > rtol_dp) then + write(msg,*) "Input matrix is not symmetric. 0.5*max(X-X.T) = ", & + & symmetry_error, ", tol = ", rtol_dp + call stop_error(msg, module=this_module, procedure='sqrtm_rdp') + else if (symmetry_error > 10*atol_dp) then + write(msg,*) "Input matrix is not exactly symmetric. 0.5*max(X-X.T) = ", symmetry_error + call logger%log_warning(trim(msg), module=this_module, procedure='sqrtm_rdp') + end if + + ! Perform svd + call svd(X, S, U, VT) + + ! Check if the matrix is positive definite (up to tol) + do i = 1, size(S) + if (S(i) .gt. 10*atol_dp ) then + S(i) = sqrt(S(i)) + else + S(i) = zero_rdp + info = 1 + end if + end do + + ! Reconstruct the square root matrix + sqrtmX = matmul(U, matmul(diag(S), VT)) + + return + end subroutine + + subroutine sqrtm_eig_rdp(X, sqrtmX, info) !! Matrix-valued sqrt function for dense symmetric/hermitian positive (semi-)definite matrices real(dp), intent(in) :: X(:,:) !! Matrix of which to compute the sqrt @@ -691,7 +799,7 @@ subroutine sqrtm_rdp(X, sqrtmX, info) real(dp) :: lambda(size(X,1)) real(dp) :: V(size(X,1), size(X,1)) integer :: i - character*128 :: msg + character(len=128) :: msg info = 0 @@ -724,7 +832,6 @@ subroutine sqrtm_rdp(X, sqrtmX, info) return end subroutine - subroutine inv_csp(A) !! In-place inversion of A using LAPACK. complex(sp), intent(inout) :: A(:, :) @@ -866,6 +973,54 @@ subroutine ordschur_csp(T, Q, selected) end subroutine ordschur_csp subroutine sqrtm_csp(X, sqrtmX, info) + !! Matrix-valued sqrt function for dense symmetric/hermitian positive (semi-)definite matrices + complex(sp), intent(inout) :: X(:,:) + !! Matrix of which to compute the sqrt + complex(sp), intent(out) :: sqrtmX(size(X,1),size(X,1)) + !! Return matrix + integer, intent(out) :: info + !! Information flag + + ! internal + real(sp) :: S(size(X,1)) + complex(sp) :: U(size(X,1), size(X,1)), VT(size(X,1), size(X,1)) + integer :: i + real(sp) :: symmetry_error + character(len=128) :: msg + + info = 0 + + ! Check if the matrix is hermitian + symmetry_error = 0.5*maxval(abs(X - conjg(transpose(X)))) + if (symmetry_error > rtol_sp) then + write(msg,*) "Input matrix is not hermitian. 0.5*max(abs(X-X.H)) = ", & + & symmetry_error, ", tol = ", rtol_sp + call stop_error(msg, module=this_module, procedure='sqrtm_csp') + else if (symmetry_error > 10*atol_sp) then + write(msg,*) "Input matrix is not exactly hermitian. 0.5*max(X-X.T) = ", symmetry_error + call logger%log_warning(trim(msg), module=this_module, procedure='sqrtm_csp') + end if + + ! Perform svd + call svd(X, S, U, VT) + + ! Check if the matrix is positive definite (up to tol) + do i = 1, size(S) + if (S(i) .gt. 10*atol_sp ) then + S(i) = sqrt(S(i)) + else + S(i) = zero_rsp + info = 1 + end if + end do + + ! Reconstruct the square root matrix + sqrtmX = matmul(U, matmul(diag(S), VT)) + + return + end subroutine + + subroutine sqrtm_eig_csp(X, sqrtmX, info) !! Matrix-valued sqrt function for dense symmetric/hermitian positive (semi-)definite matrices complex(sp), intent(in) :: X(:,:) !! Matrix of which to compute the sqrt @@ -878,7 +1033,7 @@ subroutine sqrtm_csp(X, sqrtmX, info) real(sp) :: lambda(size(X,1)) complex(sp) :: V(size(X,1), size(X,1)) integer :: i - character*128 :: msg + character(len=128) :: msg info = 0 @@ -911,7 +1066,6 @@ subroutine sqrtm_csp(X, sqrtmX, info) return end subroutine - subroutine inv_cdp(A) !! In-place inversion of A using LAPACK. complex(dp), intent(inout) :: A(:, :) @@ -1053,6 +1207,54 @@ subroutine ordschur_cdp(T, Q, selected) end subroutine ordschur_cdp subroutine sqrtm_cdp(X, sqrtmX, info) + !! Matrix-valued sqrt function for dense symmetric/hermitian positive (semi-)definite matrices + complex(dp), intent(inout) :: X(:,:) + !! Matrix of which to compute the sqrt + complex(dp), intent(out) :: sqrtmX(size(X,1),size(X,1)) + !! Return matrix + integer, intent(out) :: info + !! Information flag + + ! internal + real(dp) :: S(size(X,1)) + complex(dp) :: U(size(X,1), size(X,1)), VT(size(X,1), size(X,1)) + integer :: i + real(dp) :: symmetry_error + character(len=128) :: msg + + info = 0 + + ! Check if the matrix is hermitian + symmetry_error = 0.5*maxval(abs(X - conjg(transpose(X)))) + if (symmetry_error > rtol_dp) then + write(msg,*) "Input matrix is not hermitian. 0.5*max(abs(X-X.H)) = ", & + & symmetry_error, ", tol = ", rtol_dp + call stop_error(msg, module=this_module, procedure='sqrtm_cdp') + else if (symmetry_error > 10*atol_dp) then + write(msg,*) "Input matrix is not exactly hermitian. 0.5*max(X-X.T) = ", symmetry_error + call logger%log_warning(trim(msg), module=this_module, procedure='sqrtm_cdp') + end if + + ! Perform svd + call svd(X, S, U, VT) + + ! Check if the matrix is positive definite (up to tol) + do i = 1, size(S) + if (S(i) .gt. 10*atol_dp ) then + S(i) = sqrt(S(i)) + else + S(i) = zero_rdp + info = 1 + end if + end do + + ! Reconstruct the square root matrix + sqrtmX = matmul(U, matmul(diag(S), VT)) + + return + end subroutine + + subroutine sqrtm_eig_cdp(X, sqrtmX, info) !! Matrix-valued sqrt function for dense symmetric/hermitian positive (semi-)definite matrices complex(dp), intent(in) :: X(:,:) !! Matrix of which to compute the sqrt @@ -1065,7 +1267,7 @@ subroutine sqrtm_cdp(X, sqrtmX, info) real(dp) :: lambda(size(X,1)) complex(dp) :: V(size(X,1), size(X,1)) integer :: i - character*128 :: msg + character(len=128) :: msg info = 0 @@ -1099,7 +1301,6 @@ subroutine sqrtm_cdp(X, sqrtmX, info) return end subroutine - !--------------------------------- !----- MISCELLANEOUS ----- !--------------------------------- diff --git a/src/Utils.fypp b/src/Utils.fypp index 1653990..2f2145a 100644 --- a/src/Utils.fypp +++ b/src/Utils.fypp @@ -5,8 +5,7 @@ module lightkrylov_utils !----- Standard Fortran Library ----- !-------------------------------------------- use iso_fortran_env, only: output_unit - use LightKrylov_Logger - use stdlib_linalg, only: is_hermitian, is_symmetric, diag + use stdlib_linalg, only: is_hermitian, is_symmetric, diag, svd ! Matrix inversion. use stdlib_linalg_lapack, only: getrf, getri ! Eigenvalue problem (general + symmetric). @@ -18,14 +17,15 @@ module lightkrylov_utils !----- LightKrylov ----- !------------------------------- ! Various constants. + use LightKrylov_Logger use LightKrylov_Constants implicit none private - character*128, parameter :: this_module = 'LightKrylov_Utils' + character(len=128), parameter :: this_module = 'LightKrylov_Utils' - public :: assert_shape + public :: assert_shape, norml, log2 ! Compute B = inv(A) in-place for dense matrices. public :: inv ! Compute AX = XD for general dense matrices. @@ -34,18 +34,12 @@ module lightkrylov_utils public :: eigh ! Compute matrix sqrt of input symmetric/hermitian positive definite matrix A public :: sqrtm + public :: sqrtm_eig ! Compute AX = XS where S is in Schur form. public :: schur ! Re-orders the Schur factorization of A. public :: ordschur - #:for kind, type in RC_KINDS_TYPES - #:if type[0] == "r" - public :: log2_${type[0]}$${kind}$ - #:endif - public :: norml_${type[0]}$${kind}$ - #:endfor - interface assert_shape #:for kind, type in RC_KINDS_TYPES module procedure assert_shape_vector_${type[0]}$${kind}$ @@ -53,6 +47,18 @@ module lightkrylov_utils #:endfor end interface + interface norml + #:for kind, type in RC_KINDS_TYPES + module procedure norml_${type[0]}$${kind}$ + #:endfor + end interface + + interface log2 + #:for kind, type in REAL_KINDS_TYPES + module procedure log2_${type[0]}$${kind}$ + #:endfor + end interface + interface inv #:for kind, type in RC_KINDS_TYPES module procedure inv_${type[0]}$${kind}$ @@ -89,6 +95,12 @@ module lightkrylov_utils #:endfor end interface + interface sqrtm_eig + #:for kind, type in RC_KINDS_TYPES + module procedure sqrtm_eig_${type[0]}$${kind}$ + #:endfor + end interface + !------------------------------------------------ !----- OPTS TYPE FOR LINEAR SOLVERS ----- !------------------------------------------------ @@ -145,7 +157,7 @@ contains !! Name of the asserted vector. ! internals - character*128 :: msg + character(len=256) :: msg if(any(shape(v) /= size)) then write(msg, *) "In routine "//routine//" vector "//matname//" has illegal length ", shape(v), & @@ -167,7 +179,7 @@ contains !! Name of the asserted matrix. ! internals - character*128 :: msg + character(len=256) :: msg if(any(shape(A) /= size)) then write(msg, *) "In routine "//routine//" matrix "//matname//" has illegal shape ", shape(A), & @@ -389,6 +401,67 @@ contains end subroutine ordschur_${type[0]}$${kind}$ subroutine sqrtm_${type[0]}$${kind}$(X, sqrtmX, info) + !! Matrix-valued sqrt function for dense symmetric/hermitian positive (semi-)definite matrices + ${type}$, intent(inout) :: X(:,:) + !! Matrix of which to compute the sqrt + ${type}$, intent(out) :: sqrtmX(size(X,1),size(X,1)) + !! Return matrix + integer, intent(out) :: info + !! Information flag + + ! internal + real(${kind}$) :: S(size(X,1)) + ${type}$ :: U(size(X,1), size(X,1)), VT(size(X,1), size(X,1)) + integer :: i + real(${kind}$) :: symmetry_error + character(len=128) :: msg + + info = 0 + + #:if type[0] == "r" + ! Check if the matrix is symmetric + symmetry_error = 0.5*maxval(X - transpose(X)) + if (symmetry_error > rtol_${kind}$) then + write(msg,*) "Input matrix is not symmetric. 0.5*max(X-X.T) = ", & + & symmetry_error, ", tol = ", rtol_${kind}$ + call stop_error(msg, module=this_module, procedure='sqrtm_${type[0]}$${kind}$') + else if (symmetry_error > 10*atol_${kind}$) then + write(msg,*) "Input matrix is not exactly symmetric. 0.5*max(X-X.T) = ", symmetry_error + call logger%log_warning(trim(msg), module=this_module, procedure='sqrtm_${type[0]}$${kind}$') + end if + #:else + ! Check if the matrix is hermitian + symmetry_error = 0.5*maxval(abs(X - conjg(transpose(X)))) + if (symmetry_error > rtol_${kind}$) then + write(msg,*) "Input matrix is not hermitian. 0.5*max(abs(X-X.H)) = ", & + & symmetry_error, ", tol = ", rtol_${kind}$ + call stop_error(msg, module=this_module, procedure='sqrtm_${type[0]}$${kind}$') + else if (symmetry_error > 10*atol_${kind}$) then + write(msg,*) "Input matrix is not exactly hermitian. 0.5*max(X-X.T) = ", symmetry_error + call logger%log_warning(trim(msg), module=this_module, procedure='sqrtm_${type[0]}$${kind}$') + end if + #:endif + + ! Perform svd + call svd(X, S, U, VT) + + ! Check if the matrix is positive definite (up to tol) + do i = 1, size(S) + if (S(i) .gt. 10*atol_${kind}$ ) then + S(i) = sqrt(S(i)) + else + S(i) = zero_r${kind}$ + info = 1 + end if + end do + + ! Reconstruct the square root matrix + sqrtmX = matmul(U, matmul(diag(S), VT)) + + return + end subroutine + + subroutine sqrtm_eig_${type[0]}$${kind}$(X, sqrtmX, info) !! Matrix-valued sqrt function for dense symmetric/hermitian positive (semi-)definite matrices ${type}$, intent(in) :: X(:,:) !! Matrix of which to compute the sqrt @@ -401,7 +474,7 @@ contains real(${kind}$) :: lambda(size(X,1)) ${type}$ :: V(size(X,1), size(X,1)) integer :: i - character*128 :: msg + character(len=128) :: msg info = 0 @@ -446,7 +519,6 @@ contains return end subroutine - #:endfor !--------------------------------- diff --git a/test/TestExpmlib.f90 b/test/TestExpmlib.f90 index 4d02947..597d806 100644 --- a/test/TestExpmlib.f90 +++ b/test/TestExpmlib.f90 @@ -4,7 +4,7 @@ module TestExpmlib use stdlib_math, only: is_close, all_close use stdlib_linalg, only: eye, diag use stdlib_io_npy, only: save_npy - ! Testdrive + ! Testdrive use testdrive, only: new_unittest, unittest_type, error_type, check ! LightKrylov use LightKrylov @@ -12,8 +12,7 @@ module TestExpmlib use LightKrylov_Logger use LightKrylov_Utils, only : eig, sqrtm ! Test Utilities - use LightKrylov_TestTypes - use TestUtils + use LightKrylov_TestUtils implicit none @@ -80,7 +79,7 @@ subroutine test_dense_expm_rsp(error) err = maxval(abs(E-Eref)) call get_err_str(msg, "max err: ", err) call check(error, err < rtol_sp) - call check_test(error, 'test_dense_expm_rsp', eq='maxval(abs(E-Eref)) < rtol', context=msg) + call check_test(error, 'test_dense_expm_rsp', eq='maxval(abs(E-Eref))', context=msg) return end subroutine test_dense_expm_rsp @@ -122,7 +121,7 @@ subroutine test_kexptA_rsp(error) call get_err_str(msg, "max err: ", err) call check(error, err < rtol_sp) call check_test(error, 'test_kexptA_rsp', & - & info='Comparison with matrix exponential', context=msg) + & eq='Comparison with dense matrix exponential', context=msg) return end subroutine test_kexptA_rsp @@ -210,7 +209,7 @@ subroutine test_block_kexptA_rsp(error) call check(error, errv < rtol_sp) call check_test(error, 'test_block_kexptA_rsp', & - & info='Comparison with matrix exponential', context=msg) + & eq='Comparison with dense matrix exponential', context=msg) return end subroutine test_block_kexptA_rsp @@ -258,7 +257,7 @@ subroutine test_dense_expm_rdp(error) err = maxval(abs(E-Eref)) call get_err_str(msg, "max err: ", err) call check(error, err < rtol_dp) - call check_test(error, 'test_dense_expm_rdp', eq='maxval(abs(E-Eref)) < rtol', context=msg) + call check_test(error, 'test_dense_expm_rdp', eq='maxval(abs(E-Eref))', context=msg) return end subroutine test_dense_expm_rdp @@ -300,7 +299,7 @@ subroutine test_kexptA_rdp(error) call get_err_str(msg, "max err: ", err) call check(error, err < rtol_dp) call check_test(error, 'test_kexptA_rdp', & - & info='Comparison with matrix exponential', context=msg) + & eq='Comparison with dense matrix exponential', context=msg) return end subroutine test_kexptA_rdp @@ -388,7 +387,7 @@ subroutine test_block_kexptA_rdp(error) call check(error, errv < rtol_dp) call check_test(error, 'test_block_kexptA_rdp', & - & info='Comparison with matrix exponential', context=msg) + & eq='Comparison with dense matrix exponential', context=msg) return end subroutine test_block_kexptA_rdp @@ -436,7 +435,7 @@ subroutine test_dense_expm_csp(error) err = maxval(abs(E-Eref)) call get_err_str(msg, "max err: ", err) call check(error, err < rtol_sp) - call check_test(error, 'test_dense_expm_csp', eq='maxval(abs(E-Eref)) < rtol', context=msg) + call check_test(error, 'test_dense_expm_csp', eq='maxval(abs(E-Eref))', context=msg) return end subroutine test_dense_expm_csp @@ -478,7 +477,7 @@ subroutine test_kexptA_csp(error) call get_err_str(msg, "max err: ", err) call check(error, err < rtol_sp) call check_test(error, 'test_kexptA_csp', & - & info='Comparison with matrix exponential', context=msg) + & eq='Comparison with dense matrix exponential', context=msg) return end subroutine test_kexptA_csp @@ -566,7 +565,7 @@ subroutine test_block_kexptA_csp(error) call check(error, errv < rtol_sp) call check_test(error, 'test_block_kexptA_csp', & - & info='Comparison with matrix exponential', context=msg) + & eq='Comparison with dense matrix exponential', context=msg) return end subroutine test_block_kexptA_csp @@ -614,7 +613,7 @@ subroutine test_dense_expm_cdp(error) err = maxval(abs(E-Eref)) call get_err_str(msg, "max err: ", err) call check(error, err < rtol_dp) - call check_test(error, 'test_dense_expm_cdp', eq='maxval(abs(E-Eref)) < rtol', context=msg) + call check_test(error, 'test_dense_expm_cdp', eq='maxval(abs(E-Eref))', context=msg) return end subroutine test_dense_expm_cdp @@ -656,7 +655,7 @@ subroutine test_kexptA_cdp(error) call get_err_str(msg, "max err: ", err) call check(error, err < rtol_dp) call check_test(error, 'test_kexptA_cdp', & - & info='Comparison with matrix exponential', context=msg) + & eq='Comparison with dense matrix exponential', context=msg) return end subroutine test_kexptA_cdp @@ -744,7 +743,7 @@ subroutine test_block_kexptA_cdp(error) call check(error, errv < rtol_dp) call check_test(error, 'test_block_kexptA_cdp', & - & info='Comparison with matrix exponential', context=msg) + & eq='Comparison with dense matrix exponential', context=msg) return end subroutine test_block_kexptA_cdp @@ -790,8 +789,6 @@ subroutine test_dense_sqrtm_pos_def_rsp(error) end do ! reconstruct matrix A = matmul(sqrtmA, matmul(diag(abs(lambda)), transpose(sqrtmA))) - ! ensure it is exactly symmetric/hermitian - A = 0.5_sp*(A + transpose(A)) ! compute matrix square root call sqrtm(A, sqrtmA, info) @@ -812,7 +809,7 @@ subroutine test_dense_sqrtm_pos_semi_def_rsp(error) !> Error type to be returned. type(error_type), allocatable, intent(out) :: error !> Problem dimension. - integer, parameter :: n = 5 + integer, parameter :: n = 25 !> Test matrix. real(sp) :: A(n, n) real(sp) :: sqrtmA(n, n) @@ -832,8 +829,6 @@ subroutine test_dense_sqrtm_pos_semi_def_rsp(error) lambda(n) = zero_rsp ! reconstruct matrix A = matmul(sqrtmA, matmul(diag(abs(lambda)), transpose(sqrtmA))) - ! ensure it is exactly symmetric/hermitian - A = 0.5_sp*(A + transpose(A)) ! compute matrix square root call sqrtm(A, sqrtmA, info) @@ -884,8 +879,6 @@ subroutine test_dense_sqrtm_pos_def_rdp(error) end do ! reconstruct matrix A = matmul(sqrtmA, matmul(diag(abs(lambda)), transpose(sqrtmA))) - ! ensure it is exactly symmetric/hermitian - A = 0.5_dp*(A + transpose(A)) ! compute matrix square root call sqrtm(A, sqrtmA, info) @@ -906,7 +899,7 @@ subroutine test_dense_sqrtm_pos_semi_def_rdp(error) !> Error type to be returned. type(error_type), allocatable, intent(out) :: error !> Problem dimension. - integer, parameter :: n = 5 + integer, parameter :: n = 25 !> Test matrix. real(dp) :: A(n, n) real(dp) :: sqrtmA(n, n) @@ -926,8 +919,6 @@ subroutine test_dense_sqrtm_pos_semi_def_rdp(error) lambda(n) = zero_rdp ! reconstruct matrix A = matmul(sqrtmA, matmul(diag(abs(lambda)), transpose(sqrtmA))) - ! ensure it is exactly symmetric/hermitian - A = 0.5_dp*(A + transpose(A)) ! compute matrix square root call sqrtm(A, sqrtmA, info) @@ -979,8 +970,6 @@ subroutine test_dense_sqrtm_pos_def_csp(error) end do ! reconstruct matrix A = matmul(sqrtmA, matmul(diag(abs(lambda)), conjg(transpose(sqrtmA)))) - ! ensure it is exactly symmetric/hermitian - A = 0.5_sp*(A + conjg(transpose(A))) ! compute matrix square root call sqrtm(A, sqrtmA, info) @@ -1001,7 +990,7 @@ subroutine test_dense_sqrtm_pos_semi_def_csp(error) !> Error type to be returned. type(error_type), allocatable, intent(out) :: error !> Problem dimension. - integer, parameter :: n = 5 + integer, parameter :: n = 25 !> Test matrix. complex(sp) :: A(n, n) complex(sp) :: sqrtmA(n, n) @@ -1022,8 +1011,6 @@ subroutine test_dense_sqrtm_pos_semi_def_csp(error) lambda(n) = zero_rsp ! reconstruct matrix A = matmul(sqrtmA, matmul(diag(abs(lambda)), conjg(transpose(sqrtmA)))) - ! ensure it is exactly symmetric/hermitian - A = 0.5_sp*(A + conjg(transpose(A))) ! compute matrix square root call sqrtm(A, sqrtmA, info) @@ -1075,8 +1062,6 @@ subroutine test_dense_sqrtm_pos_def_cdp(error) end do ! reconstruct matrix A = matmul(sqrtmA, matmul(diag(abs(lambda)), conjg(transpose(sqrtmA)))) - ! ensure it is exactly symmetric/hermitian - A = 0.5_dp*(A + conjg(transpose(A))) ! compute matrix square root call sqrtm(A, sqrtmA, info) @@ -1097,7 +1082,7 @@ subroutine test_dense_sqrtm_pos_semi_def_cdp(error) !> Error type to be returned. type(error_type), allocatable, intent(out) :: error !> Problem dimension. - integer, parameter :: n = 5 + integer, parameter :: n = 25 !> Test matrix. complex(dp) :: A(n, n) complex(dp) :: sqrtmA(n, n) @@ -1118,8 +1103,6 @@ subroutine test_dense_sqrtm_pos_semi_def_cdp(error) lambda(n) = zero_rdp ! reconstruct matrix A = matmul(sqrtmA, matmul(diag(abs(lambda)), conjg(transpose(sqrtmA)))) - ! ensure it is exactly symmetric/hermitian - A = 0.5_dp*(A + conjg(transpose(A))) ! compute matrix square root call sqrtm(A, sqrtmA, info) diff --git a/test/TestExpmlib.fypp b/test/TestExpmlib.fypp index d37831b..97c6de0 100644 --- a/test/TestExpmlib.fypp +++ b/test/TestExpmlib.fypp @@ -6,7 +6,7 @@ module TestExpmlib use stdlib_math, only: is_close, all_close use stdlib_linalg, only: eye, diag use stdlib_io_npy, only: save_npy - ! Testdrive + ! Testdrive use testdrive, only: new_unittest, unittest_type, error_type, check ! LightKrylov use LightKrylov @@ -14,8 +14,7 @@ module TestExpmlib use LightKrylov_Logger use LightKrylov_Utils, only : eig, sqrtm ! Test Utilities - use LightKrylov_TestTypes - use TestUtils + use LightKrylov_TestUtils implicit none @@ -81,7 +80,7 @@ contains err = maxval(abs(E-Eref)) call get_err_str(msg, "max err: ", err) call check(error, err < rtol_${kind}$) - call check_test(error, 'test_dense_expm_${type[0]}$${kind}$', eq='maxval(abs(E-Eref)) < rtol', context=msg) + call check_test(error, 'test_dense_expm_${type[0]}$${kind}$', eq='maxval(abs(E-Eref))', context=msg) return end subroutine test_dense_expm_${type[0]}$${kind}$ @@ -123,7 +122,7 @@ contains call get_err_str(msg, "max err: ", err) call check(error, err < rtol_${kind}$) call check_test(error, 'test_kexptA_${type[0]}$${kind}$', & - & info='Comparison with matrix exponential', context=msg) + & eq='Comparison with dense matrix exponential', context=msg) return end subroutine test_kexptA_${type[0]}$${kind}$ @@ -211,7 +210,7 @@ contains call check(error, errv < rtol_${kind}$) call check_test(error, 'test_block_kexptA_${type[0]}$${kind}$', & - & info='Comparison with matrix exponential', context=msg) + & eq='Comparison with dense matrix exponential', context=msg) return end subroutine test_block_kexptA_${type[0]}$${kind}$ @@ -269,12 +268,8 @@ contains ! reconstruct matrix #:if type[0] == "r" A = matmul(sqrtmA, matmul(diag(abs(lambda)), transpose(sqrtmA))) - ! ensure it is exactly symmetric/hermitian - A = 0.5_${kind}$*(A + transpose(A)) #:else A = matmul(sqrtmA, matmul(diag(abs(lambda)), conjg(transpose(sqrtmA)))) - ! ensure it is exactly symmetric/hermitian - A = 0.5_${kind}$*(A + conjg(transpose(A))) #:endif ! compute matrix square root @@ -302,7 +297,7 @@ contains !> Error type to be returned. type(error_type), allocatable, intent(out) :: error !> Problem dimension. - integer, parameter :: n = 5 + integer, parameter :: n = 25 !> Test matrix. ${type}$ :: A(n, n) ${type}$ :: sqrtmA(n, n) @@ -332,12 +327,8 @@ contains ! reconstruct matrix #:if type[0] == "r" A = matmul(sqrtmA, matmul(diag(abs(lambda)), transpose(sqrtmA))) - ! ensure it is exactly symmetric/hermitian - A = 0.5_${kind}$*(A + transpose(A)) #:else A = matmul(sqrtmA, matmul(diag(abs(lambda)), conjg(transpose(sqrtmA)))) - ! ensure it is exactly symmetric/hermitian - A = 0.5_${kind}$*(A + conjg(transpose(A))) #:endif ! compute matrix square root diff --git a/test/TestIterativeSolvers.f90 b/test/TestIterativeSolvers.f90 index b81799e..0ee7db6 100644 --- a/test/TestIterativeSolvers.f90 +++ b/test/TestIterativeSolvers.f90 @@ -13,8 +13,7 @@ module TestIterativeSolvers use LightKrylov_Logger use LightKrylov_AbstractVectors ! Test Utilities - use LightKrylov_TestTypes - use TestUtils + use LightKrylov_TestUtils implicit none diff --git a/test/TestIterativeSolvers.fypp b/test/TestIterativeSolvers.fypp index 82dac9e..578b3c0 100644 --- a/test/TestIterativeSolvers.fypp +++ b/test/TestIterativeSolvers.fypp @@ -15,8 +15,7 @@ module TestIterativeSolvers use LightKrylov_Logger use LightKrylov_AbstractVectors ! Test Utilities - use LightKrylov_TestTypes - use TestUtils + use LightKrylov_TestUtils implicit none diff --git a/test/TestKrylov.f90 b/test/TestKrylov.f90 index 2951025..d988199 100644 --- a/test/TestKrylov.f90 +++ b/test/TestKrylov.f90 @@ -12,8 +12,7 @@ module TestKrylov use LightKrylov_Logger use LightKrylov_AbstractVectors ! Test Utilities - use LightKrylov_TestTypes - use TestUtils + use LightKrylov_TestUtils implicit none diff --git a/test/TestKrylov.fypp b/test/TestKrylov.fypp index 8076d0f..b1e453f 100644 --- a/test/TestKrylov.fypp +++ b/test/TestKrylov.fypp @@ -14,8 +14,7 @@ module TestKrylov use LightKrylov_Logger use LightKrylov_AbstractVectors ! Test Utilities - use LightKrylov_TestTypes - use TestUtils + use LightKrylov_TestUtils implicit none diff --git a/test/TestLinops.f90 b/test/TestLinops.f90 index 040e0e3..849e409 100644 --- a/test/TestLinops.f90 +++ b/test/TestLinops.f90 @@ -7,7 +7,8 @@ module TestLinops use LightKrylov use LightKrylov_Constants use LightKrylov_Logger - use LightKrylov_TestTypes + ! Test Utilities + use LightKrylov_TestUtils implicit none diff --git a/test/TestLinops.fypp b/test/TestLinops.fypp index c04ce7b..7182c3d 100644 --- a/test/TestLinops.fypp +++ b/test/TestLinops.fypp @@ -9,7 +9,8 @@ module TestLinops use LightKrylov use LightKrylov_Constants use LightKrylov_Logger - use LightKrylov_TestTypes + ! Test Utilities + use LightKrylov_TestUtils implicit none diff --git a/test/TestUtils.f90 b/test/TestUtils.f90 deleted file mode 100644 index 7832fd1..0000000 --- a/test/TestUtils.f90 +++ /dev/null @@ -1,478 +0,0 @@ -module TestUtils - use stdlib_io_npy, only: save_npy - use stdlib_linalg, only: eye, diag - use stdlib_stats_distribution_normal, only: normal => rvs_normal - use LightKrylov - use LightKrylov_Constants - use LightKrylov_TestTypes - - implicit none - - private - - character(len=128), parameter, private :: this_module = 'LightKrylov_TestUtils' - - public :: get_data - public :: put_data - public :: init_rand - public :: get_err_str - - interface get_data - module procedure get_data_vec_rsp - module procedure get_data_vec_basis_rsp - module procedure get_data_linop_rsp - module procedure get_data_vec_rdp - module procedure get_data_vec_basis_rdp - module procedure get_data_linop_rdp - module procedure get_data_vec_csp - module procedure get_data_vec_basis_csp - module procedure get_data_linop_csp - module procedure get_data_vec_cdp - module procedure get_data_vec_basis_cdp - module procedure get_data_linop_cdp - end interface - - interface put_data - module procedure put_data_vec_rsp - module procedure put_data_vec_basis_rsp - module procedure put_data_linop_rsp - module procedure put_data_vec_rdp - module procedure put_data_vec_basis_rdp - module procedure put_data_linop_rdp - module procedure put_data_vec_csp - module procedure put_data_vec_basis_csp - module procedure put_data_linop_csp - module procedure put_data_vec_cdp - module procedure put_data_vec_basis_cdp - module procedure put_data_linop_cdp - end interface - - interface init_rand - module procedure init_rand_vec_rsp - module procedure init_rand_basis_rsp - module procedure init_rand_linop_rsp - module procedure init_rand_spd_linop_rsp - module procedure init_rand_vec_rdp - module procedure init_rand_basis_rdp - module procedure init_rand_linop_rdp - module procedure init_rand_spd_linop_rdp - module procedure init_rand_vec_csp - module procedure init_rand_basis_csp - module procedure init_rand_linop_csp - module procedure init_rand_hermitian_linop_csp - module procedure init_rand_vec_cdp - module procedure init_rand_basis_cdp - module procedure init_rand_linop_cdp - module procedure init_rand_hermitian_linop_cdp - end interface - - interface get_err_str - module procedure get_err_str_sp - module procedure get_err_str_dp - end interface - -contains - - !---------------------------------------------------- - !----- EXTRACT DATA FROM ABSTRACT TYPES ----- - !---------------------------------------------------- - - subroutine get_data_vec_rsp(vec_out, vec_in) - real(sp), intent(out) :: vec_out(:) - type(vector_rsp), intent(in) :: vec_in - vec_out = vec_in%data - return - end subroutine get_data_vec_rsp - - subroutine get_data_vec_basis_rsp(basis_out, basis_in) - real(sp), intent(out) :: basis_out(:, :) - type(vector_rsp), intent(in) :: basis_in(:) - ! Internal variables. - integer :: k - do k = 1, size(basis_in) - basis_out(:, k) = basis_in(k)%data - enddo - return - end subroutine get_data_vec_basis_rsp - - subroutine get_data_linop_rsp(mat_out, linop_in) - real(sp), intent(out) :: mat_out(:, :) - type(linop_rsp), intent(in) :: linop_in - mat_out = linop_in%data - return - end subroutine get_data_linop_rsp - - subroutine get_data_vec_rdp(vec_out, vec_in) - real(dp), intent(out) :: vec_out(:) - type(vector_rdp), intent(in) :: vec_in - vec_out = vec_in%data - return - end subroutine get_data_vec_rdp - - subroutine get_data_vec_basis_rdp(basis_out, basis_in) - real(dp), intent(out) :: basis_out(:, :) - type(vector_rdp), intent(in) :: basis_in(:) - ! Internal variables. - integer :: k - do k = 1, size(basis_in) - basis_out(:, k) = basis_in(k)%data - enddo - return - end subroutine get_data_vec_basis_rdp - - subroutine get_data_linop_rdp(mat_out, linop_in) - real(dp), intent(out) :: mat_out(:, :) - type(linop_rdp), intent(in) :: linop_in - mat_out = linop_in%data - return - end subroutine get_data_linop_rdp - - subroutine get_data_vec_csp(vec_out, vec_in) - complex(sp), intent(out) :: vec_out(:) - type(vector_csp), intent(in) :: vec_in - vec_out = vec_in%data - return - end subroutine get_data_vec_csp - - subroutine get_data_vec_basis_csp(basis_out, basis_in) - complex(sp), intent(out) :: basis_out(:, :) - type(vector_csp), intent(in) :: basis_in(:) - ! Internal variables. - integer :: k - do k = 1, size(basis_in) - basis_out(:, k) = basis_in(k)%data - enddo - return - end subroutine get_data_vec_basis_csp - - subroutine get_data_linop_csp(mat_out, linop_in) - complex(sp), intent(out) :: mat_out(:, :) - type(linop_csp), intent(in) :: linop_in - mat_out = linop_in%data - return - end subroutine get_data_linop_csp - - subroutine get_data_vec_cdp(vec_out, vec_in) - complex(dp), intent(out) :: vec_out(:) - type(vector_cdp), intent(in) :: vec_in - vec_out = vec_in%data - return - end subroutine get_data_vec_cdp - - subroutine get_data_vec_basis_cdp(basis_out, basis_in) - complex(dp), intent(out) :: basis_out(:, :) - type(vector_cdp), intent(in) :: basis_in(:) - ! Internal variables. - integer :: k - do k = 1, size(basis_in) - basis_out(:, k) = basis_in(k)%data - enddo - return - end subroutine get_data_vec_basis_cdp - - subroutine get_data_linop_cdp(mat_out, linop_in) - complex(dp), intent(out) :: mat_out(:, :) - type(linop_cdp), intent(in) :: linop_in - mat_out = linop_in%data - return - end subroutine get_data_linop_cdp - - - !---------------------------------------------- - !----- PUT DATA TO ABSTRACT TYPES ----- - !---------------------------------------------- - - subroutine put_data_vec_rsp(vec_out, vec_in) - type(vector_rsp), intent(out) :: vec_out - real(sp), intent(in) :: vec_in - vec_out%data = vec_in - return - end subroutine put_data_vec_rsp - - subroutine put_data_vec_basis_rsp(basis_out, basis_in) - type(vector_rsp), intent(out) :: basis_out(:) - real(sp), intent(in) :: basis_in(:, :) - ! Internal variables. - integer :: k - do k = 1, size(basis_out) - basis_out(k)%data = basis_in(:, k) - enddo - return - end subroutine put_data_vec_basis_rsp - - subroutine put_data_linop_rsp(linop_out, mat_in) - type(linop_rsp), intent(out) :: linop_out - real(sp), intent(in) :: mat_in(:, :) - ! Internal variables. - linop_out%data = mat_in - return - end subroutine put_data_linop_rsp - - subroutine put_data_vec_rdp(vec_out, vec_in) - type(vector_rdp), intent(out) :: vec_out - real(dp), intent(in) :: vec_in - vec_out%data = vec_in - return - end subroutine put_data_vec_rdp - - subroutine put_data_vec_basis_rdp(basis_out, basis_in) - type(vector_rdp), intent(out) :: basis_out(:) - real(dp), intent(in) :: basis_in(:, :) - ! Internal variables. - integer :: k - do k = 1, size(basis_out) - basis_out(k)%data = basis_in(:, k) - enddo - return - end subroutine put_data_vec_basis_rdp - - subroutine put_data_linop_rdp(linop_out, mat_in) - type(linop_rdp), intent(out) :: linop_out - real(dp), intent(in) :: mat_in(:, :) - ! Internal variables. - linop_out%data = mat_in - return - end subroutine put_data_linop_rdp - - subroutine put_data_vec_csp(vec_out, vec_in) - type(vector_csp), intent(out) :: vec_out - complex(sp), intent(in) :: vec_in - vec_out%data = vec_in - return - end subroutine put_data_vec_csp - - subroutine put_data_vec_basis_csp(basis_out, basis_in) - type(vector_csp), intent(out) :: basis_out(:) - complex(sp), intent(in) :: basis_in(:, :) - ! Internal variables. - integer :: k - do k = 1, size(basis_out) - basis_out(k)%data = basis_in(:, k) - enddo - return - end subroutine put_data_vec_basis_csp - - subroutine put_data_linop_csp(linop_out, mat_in) - type(linop_csp), intent(out) :: linop_out - complex(sp), intent(in) :: mat_in(:, :) - ! Internal variables. - linop_out%data = mat_in - return - end subroutine put_data_linop_csp - - subroutine put_data_vec_cdp(vec_out, vec_in) - type(vector_cdp), intent(out) :: vec_out - complex(dp), intent(in) :: vec_in - vec_out%data = vec_in - return - end subroutine put_data_vec_cdp - - subroutine put_data_vec_basis_cdp(basis_out, basis_in) - type(vector_cdp), intent(out) :: basis_out(:) - complex(dp), intent(in) :: basis_in(:, :) - ! Internal variables. - integer :: k - do k = 1, size(basis_out) - basis_out(k)%data = basis_in(:, k) - enddo - return - end subroutine put_data_vec_basis_cdp - - subroutine put_data_linop_cdp(linop_out, mat_in) - type(linop_cdp), intent(out) :: linop_out - complex(dp), intent(in) :: mat_in(:, :) - ! Internal variables. - linop_out%data = mat_in - return - end subroutine put_data_linop_cdp - - - !-------------------------------------------------------------- - !----- INITIALIZE ABSTRACT TYPES WITH RANDOM DATA ----- - !-------------------------------------------------------------- - - subroutine init_rand_vec_rsp(x) - type(vector_rsp), intent(inout) :: x - call x%rand() - return - end subroutine init_rand_vec_rsp - - subroutine init_rand_basis_rsp(X) - type(vector_rsp), intent(inout) :: X(:) - integer :: i - do i = 1, size(X) - call X(i)%rand() - enddo - return - end subroutine init_rand_basis_rsp - - subroutine init_rand_linop_rsp(linop) - type(linop_rsp), intent(inout) :: linop - real(sp), allocatable :: mu(:, :), var(:, :) - allocate(mu(test_size, test_size)) ; mu = 0.0_sp - allocate(var(test_size, test_size)) - var = 1.0_sp - linop%data = normal(mu, var) - return - end subroutine init_rand_linop_rsp - - subroutine init_rand_spd_linop_rsp(linop) - type(spd_linop_rsp), intent(inout) :: linop - real(sp), allocatable :: mu(:, :), var(:, :) - real(sp), allocatable :: data(:, :) - allocate(mu(test_size, test_size)) ; mu = zero_rsp - allocate(var(test_size, test_size)) ; var = one_rsp - - data = normal(mu, var) - linop%data = matmul(data, transpose(data))/test_size + 0.01*eye(test_size) - - return - end subroutine init_rand_spd_linop_rsp - - subroutine init_rand_vec_rdp(x) - type(vector_rdp), intent(inout) :: x - call x%rand() - return - end subroutine init_rand_vec_rdp - - subroutine init_rand_basis_rdp(X) - type(vector_rdp), intent(inout) :: X(:) - integer :: i - do i = 1, size(X) - call X(i)%rand() - enddo - return - end subroutine init_rand_basis_rdp - - subroutine init_rand_linop_rdp(linop) - type(linop_rdp), intent(inout) :: linop - real(dp), allocatable :: mu(:, :), var(:, :) - allocate(mu(test_size, test_size)) ; mu = 0.0_dp - allocate(var(test_size, test_size)) - var = 1.0_dp - linop%data = normal(mu, var) - return - end subroutine init_rand_linop_rdp - - subroutine init_rand_spd_linop_rdp(linop) - type(spd_linop_rdp), intent(inout) :: linop - real(dp), allocatable :: mu(:, :), var(:, :) - real(dp), allocatable :: data(:, :) - allocate(mu(test_size, test_size)) ; mu = zero_rdp - allocate(var(test_size, test_size)) ; var = one_rdp - - data = normal(mu, var) - linop%data = matmul(data, transpose(data))/test_size + 0.01*eye(test_size) - - return - end subroutine init_rand_spd_linop_rdp - - subroutine init_rand_vec_csp(x) - type(vector_csp), intent(inout) :: x - call x%rand() - return - end subroutine init_rand_vec_csp - - subroutine init_rand_basis_csp(X) - type(vector_csp), intent(inout) :: X(:) - integer :: i - do i = 1, size(X) - call X(i)%rand() - enddo - return - end subroutine init_rand_basis_csp - - subroutine init_rand_linop_csp(linop) - type(linop_csp), intent(inout) :: linop - complex(sp), allocatable :: mu(:, :), var(:, :) - allocate(mu(test_size, test_size)) ; mu = 0.0_sp - allocate(var(test_size, test_size)) - var = cmplx(1.0_sp, 1.0_sp, kind=sp) - linop%data = normal(mu, var) - return - end subroutine init_rand_linop_csp - - subroutine init_rand_hermitian_linop_csp(linop) - type(hermitian_linop_csp), intent(inout) :: linop - complex(sp), allocatable :: data(:, :) - complex(sp), allocatable :: mu(:, :), var(:, :) - - allocate(mu(test_size, test_size)) ; mu = 0.0_sp - allocate(var(test_size, test_size)) ; var = cmplx(1.0_sp, 1.0_sp, kind=sp) - - data = normal(mu, var) - data = matmul(data, transpose(conjg(data)))/test_size + 0.01*eye(test_size) - linop%data = data - - return - end subroutine init_rand_hermitian_linop_csp - - subroutine init_rand_vec_cdp(x) - type(vector_cdp), intent(inout) :: x - call x%rand() - return - end subroutine init_rand_vec_cdp - - subroutine init_rand_basis_cdp(X) - type(vector_cdp), intent(inout) :: X(:) - integer :: i - do i = 1, size(X) - call X(i)%rand() - enddo - return - end subroutine init_rand_basis_cdp - - subroutine init_rand_linop_cdp(linop) - type(linop_cdp), intent(inout) :: linop - complex(dp), allocatable :: mu(:, :), var(:, :) - allocate(mu(test_size, test_size)) ; mu = 0.0_dp - allocate(var(test_size, test_size)) - var = cmplx(1.0_dp, 1.0_dp, kind=dp) - linop%data = normal(mu, var) - return - end subroutine init_rand_linop_cdp - - subroutine init_rand_hermitian_linop_cdp(linop) - type(hermitian_linop_cdp), intent(inout) :: linop - complex(dp), allocatable :: data(:, :) - complex(dp), allocatable :: mu(:, :), var(:, :) - - allocate(mu(test_size, test_size)) ; mu = 0.0_dp - allocate(var(test_size, test_size)) ; var = cmplx(1.0_dp, 1.0_dp, kind=dp) - - data = normal(mu, var) - data = matmul(data, transpose(conjg(data)))/test_size + 0.01*eye(test_size) - linop%data = data - - return - end subroutine init_rand_hermitian_linop_cdp - - - subroutine get_err_str_sp(msg, info, err) - character(len=*), intent(inout) :: msg - character(len=*), intent(in) :: info - real(sp) :: err - - ! internals - character*8 :: value_str - character(len=*), parameter :: indent = repeat(" ", 4) - - write(value_str, '(E8.2)') err - msg = indent // info // value_str // achar(10) - - end subroutine get_err_str_sp - subroutine get_err_str_dp(msg, info, err) - character(len=*), intent(inout) :: msg - character(len=*), intent(in) :: info - real(dp) :: err - - ! internals - character*8 :: value_str - character(len=*), parameter :: indent = repeat(" ", 4) - - write(value_str, '(E8.2)') err - msg = indent // info // value_str // achar(10) - - end subroutine get_err_str_dp - -end module diff --git a/test/TestUtils.fypp b/test/TestUtils.fypp deleted file mode 100644 index 98145d7..0000000 --- a/test/TestUtils.fypp +++ /dev/null @@ -1,206 +0,0 @@ -#:include "../include/common.fypp" -#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES -module TestUtils - use stdlib_io_npy, only: save_npy - use stdlib_linalg, only: eye, diag - use stdlib_stats_distribution_normal, only: normal => rvs_normal - use LightKrylov - use LightKrylov_Constants - use LightKrylov_TestTypes - - implicit none - - private - - character(len=128), parameter, private :: this_module = 'LightKrylov_TestUtils' - - public :: get_data - public :: put_data - public :: init_rand - public :: get_err_str - - interface get_data - #:for kind, type in RC_KINDS_TYPES - module procedure get_data_vec_${type[0]}$${kind}$ - module procedure get_data_vec_basis_${type[0]}$${kind}$ - module procedure get_data_linop_${type[0]}$${kind}$ - #:endfor - end interface - - interface put_data - #:for kind, type in RC_KINDS_TYPES - module procedure put_data_vec_${type[0]}$${kind}$ - module procedure put_data_vec_basis_${type[0]}$${kind}$ - module procedure put_data_linop_${type[0]}$${kind}$ - #:endfor - end interface - - interface init_rand - #:for kind, type in RC_KINDS_TYPES - module procedure init_rand_vec_${type[0]}$${kind}$ - module procedure init_rand_basis_${type[0]}$${kind}$ - module procedure init_rand_linop_${type[0]}$${kind}$ - #:if type[0] == "r" - module procedure init_rand_spd_linop_${type[0]}$${kind}$ - #:else - module procedure init_rand_hermitian_linop_${type[0]}$${kind}$ - #:endif - #:endfor - end interface - - interface get_err_str - #:for kind, type in REAL_KINDS_TYPES - module procedure get_err_str_${kind}$ - #:endfor - end interface - -contains - - !---------------------------------------------------- - !----- EXTRACT DATA FROM ABSTRACT TYPES ----- - !---------------------------------------------------- - - #:for kind, type in RC_KINDS_TYPES - subroutine get_data_vec_${type[0]}$${kind}$(vec_out, vec_in) - ${type}$, intent(out) :: vec_out(:) - type(vector_${type[0]}$${kind}$), intent(in) :: vec_in - vec_out = vec_in%data - return - end subroutine get_data_vec_${type[0]}$${kind}$ - - subroutine get_data_vec_basis_${type[0]}$${kind}$(basis_out, basis_in) - ${type}$, intent(out) :: basis_out(:, :) - type(vector_${type[0]}$${kind}$), intent(in) :: basis_in(:) - ! Internal variables. - integer :: k - do k = 1, size(basis_in) - basis_out(:, k) = basis_in(k)%data - enddo - return - end subroutine get_data_vec_basis_${type[0]}$${kind}$ - - subroutine get_data_linop_${type[0]}$${kind}$(mat_out, linop_in) - ${type}$, intent(out) :: mat_out(:, :) - type(linop_${type[0]}$${kind}$), intent(in) :: linop_in - mat_out = linop_in%data - return - end subroutine get_data_linop_${type[0]}$${kind}$ - - #:endfor - - !---------------------------------------------- - !----- PUT DATA TO ABSTRACT TYPES ----- - !---------------------------------------------- - - #:for kind, type in RC_KINDS_TYPES - subroutine put_data_vec_${type[0]}$${kind}$(vec_out, vec_in) - type(vector_${type[0]}$${kind}$), intent(out) :: vec_out - ${type}$, intent(in) :: vec_in - vec_out%data = vec_in - return - end subroutine put_data_vec_${type[0]}$${kind}$ - - subroutine put_data_vec_basis_${type[0]}$${kind}$(basis_out, basis_in) - type(vector_${type[0]}$${kind}$), intent(out) :: basis_out(:) - ${type}$, intent(in) :: basis_in(:, :) - ! Internal variables. - integer :: k - do k = 1, size(basis_out) - basis_out(k)%data = basis_in(:, k) - enddo - return - end subroutine put_data_vec_basis_${type[0]}$${kind}$ - - subroutine put_data_linop_${type[0]}$${kind}$(linop_out, mat_in) - type(linop_${type[0]}$${kind}$), intent(out) :: linop_out - ${type}$, intent(in) :: mat_in(:, :) - ! Internal variables. - linop_out%data = mat_in - return - end subroutine put_data_linop_${type[0]}$${kind}$ - - #:endfor - - !-------------------------------------------------------------- - !----- INITIALIZE ABSTRACT TYPES WITH RANDOM DATA ----- - !-------------------------------------------------------------- - - #:for kind, type in RC_KINDS_TYPES - subroutine init_rand_vec_${type[0]}$${kind}$(x) - type(vector_${type[0]}$${kind}$), intent(inout) :: x - call x%rand() - return - end subroutine init_rand_vec_${type[0]}$${kind}$ - - subroutine init_rand_basis_${type[0]}$${kind}$(X) - type(vector_${type[0]}$${kind}$), intent(inout) :: X(:) - integer :: i - do i = 1, size(X) - call X(i)%rand() - enddo - return - end subroutine init_rand_basis_${type[0]}$${kind}$ - - subroutine init_rand_linop_${type[0]}$${kind}$(linop) - type(linop_${type[0]}$${kind}$), intent(inout) :: linop - ${type}$, allocatable :: mu(:, :), var(:, :) - allocate(mu(test_size, test_size)) ; mu = 0.0_${kind}$ - allocate(var(test_size, test_size)) - #:if type[0] == "r" - var = 1.0_${kind}$ - #:else - var = cmplx(1.0_${kind}$, 1.0_${kind}$, kind=${kind}$) - #:endif - linop%data = normal(mu, var) - return - end subroutine init_rand_linop_${type[0]}$${kind}$ - - #:if type[0] == "r" - subroutine init_rand_spd_linop_${type[0]}$${kind}$(linop) - type(spd_linop_${type[0]}$${kind}$), intent(inout) :: linop - ${type}$, allocatable :: mu(:, :), var(:, :) - ${type}$, allocatable :: data(:, :) - allocate(mu(test_size, test_size)) ; mu = zero_${type[0]}$${kind}$ - allocate(var(test_size, test_size)) ; var = one_${type[0]}$${kind}$ - - data = normal(mu, var) - linop%data = matmul(data, transpose(data))/test_size + 0.01*eye(test_size) - - return - end subroutine init_rand_spd_linop_${type[0]}$${kind}$ - #:else - subroutine init_rand_hermitian_linop_${type[0]}$${kind}$(linop) - type(hermitian_linop_${type[0]}$${kind}$), intent(inout) :: linop - ${type}$, allocatable :: data(:, :) - ${type}$, allocatable :: mu(:, :), var(:, :) - - allocate(mu(test_size, test_size)) ; mu = 0.0_${kind}$ - allocate(var(test_size, test_size)) ; var = cmplx(1.0_${kind}$, 1.0_${kind}$, kind=${kind}$) - - data = normal(mu, var) - data = matmul(data, transpose(conjg(data)))/test_size + 0.01*eye(test_size) - linop%data = data - - return - end subroutine init_rand_hermitian_linop_${type[0]}$${kind}$ - #:endif - - #:endfor - - #:for kind, type in REAL_KINDS_TYPES - subroutine get_err_str_${kind}$(msg, info, err) - character(len=*), intent(inout) :: msg - character(len=*), intent(in) :: info - real(${kind}$) :: err - - ! internals - character*8 :: value_str - character(len=*), parameter :: indent = repeat(" ", 4) - - write(value_str, '(E8.2)') err - msg = indent // info // value_str // achar(10) - - end subroutine get_err_str_${kind}$ - #:endfor - -end module diff --git a/test/TestVectors.f90 b/test/TestVectors.f90 index ce440a7..2471a33 100644 --- a/test/TestVectors.f90 +++ b/test/TestVectors.f90 @@ -8,7 +8,8 @@ module TestVectors ! LightKrylov use LightKrylov use LightKrylov_Logger - use LightKrylov_TestTypes + ! Test Utilities + use LightKrylov_TestUtils implicit none diff --git a/test/TestVectors.fypp b/test/TestVectors.fypp index d104ccf..2884439 100644 --- a/test/TestVectors.fypp +++ b/test/TestVectors.fypp @@ -10,7 +10,8 @@ module TestVectors ! LightKrylov use LightKrylov use LightKrylov_Logger - use LightKrylov_TestTypes + ! Test Utilities + use LightKrylov_TestUtils implicit none