diff --git a/src/ExpmLib.f90 b/src/ExpmLib.f90 index 1acd025..6598254 100644 --- a/src/ExpmLib.f90 +++ b/src/ExpmLib.f90 @@ -8,12 +8,12 @@ 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 @@ -25,6 +25,7 @@ module lightkrylov_expmlib public :: abstract_exptA_cdp public :: expm public :: kexpm + public :: kexpm_var_dt public :: k_exptA abstract interface @@ -132,6 +133,17 @@ end subroutine abstract_exptA_cdp module procedure kexpm_mat_cdp end interface + interface kexpm_var_dt + module procedure kexpm_vec_variable_dt_rsp + module procedure kexpm_mat_variable_dt_rsp + module procedure kexpm_vec_variable_dt_rdp + module procedure kexpm_mat_variable_dt_rdp + module procedure kexpm_vec_variable_dt_csp + module procedure kexpm_mat_variable_dt_csp + module procedure kexpm_vec_variable_dt_cdp + module procedure kexpm_mat_variable_dt_cdp + end interface + interface k_exptA module procedure k_exptA_rsp module procedure k_exptA_rdp @@ -211,7 +223,276 @@ subroutine expm_rsp(E, A, order) return end subroutine expm_rsp - subroutine kexpm_vec_rsp(c, A, b, tau, dt, tol, info, trans, verbosity, kdim) + subroutine kexpm_vec_rsp(c, A, b, tau, tol, info, trans, verbosity, kdim) + class(abstract_vector_rsp), intent(out) :: c + !! Best approximation of \( \exp(\tau \mathbf{A}) \mathbf{b} \) in the computed Krylov subspace + class(abstract_linop_rsp), intent(in) :: A + !! Linear operator to be exponentiated. + class(abstract_vector_rsp), intent(in) :: b + !! Input vector on which to apply \( \exp(\tau \mathbf{A}) \). + real(sp), intent(in) :: tau + !! Time horizon for the exponentiation. + real(sp), intent(in) :: tol + !! Solution tolerance based on error estimates. + integer, intent(out) :: info + !! Information flag. + logical, optional, intent(in) :: trans + !! Use transpose? + logical, optional, intent(in) :: verbosity + !! Verbosity control. + integer, optional, intent(in) :: kdim + !! Maximum size of the Krylov subspace. + + ! ----- Internal variables ----- + integer, parameter :: kmax = 100 + integer :: k, km, kp, nk + ! Arnoldi factorization. + class(abstract_vector_rsp), allocatable :: X(:) + real(sp), allocatable :: H(:, :) + ! Normaliztion & temporary arrays. + real(sp), allocatable :: E(:, :) + class(abstract_vector_rsp), allocatable :: Xwrk + real(sp) :: err_est, beta + character(len=512) :: msg + + ! Optional arguments. + logical :: transpose + logical :: verbose + integer :: nsteps + + ! Deals with optional args. + transpose = optval(trans, .false.) + verbose = optval(verbosity, .false.) + nsteps = optval(kdim, kmax) + nk = nsteps + + info = 0 + + ! Allocate arrays. + allocate(X(nk+1), source=b) ; allocate(H(nk+1, nk+1)) + allocate(E(nk+1, nk+1)) ; allocate(Xwrk, source=b) + + ! Normalize input vector and initialize Krylov subspace. + beta = b%norm() + if (beta == 0.0_sp) then + ! Input is zero => Output is zero. + call c%zero() + err_est = 0.0_sp + kp = 1 + else + call initialize_krylov_subspace(X) + call X(1)%add(b) ; call X(1)%scal(one_rsp / beta) + H = 0.0_sp + + expm_arnoldi: do k = 1, nk + km = k-1 ; kp = k+1 + + ! Reset work arrays. + E = 0.0_sp + + ! Compute k-th step Arnoldi factorization. + call arnoldi(A, X, H, info, kstart=k, kend=k, transpose=transpose) + call check_info(info, 'arnoldi', module=this_module, procedure='kexpm_vec_rsp') + + ! Compute approximation. + if (info == k) then + ! Arnoldi breakdown, do not consider extended matrix. + kp = k + info = -2 + endif + + ! Compute the (dense) matrix exponential of the extended Hessenberg matrix. + call expm(E(:kp, :kp), tau*H(:kp, :kp)) + + ! Project back into original space. + call linear_combination(Xwrk, X(:kp), E(:kp, 1)) + call c%axpby(zero_rsp, Xwrk, one_rsp*beta) + + ! Cheap error esimate (this actually is the magnitude of the included correction + ! and thus is very conservative). + if (info == k) then + ! Approximation is exact. + err_est = 0.0_sp + else + err_est = abs(E(kp, 1) * beta) + endif + + ! Check convergence. + if (err_est <= tol) exit expm_arnoldi + + enddo expm_arnoldi + endif + + if (err_est <= tol) then + info = kp + write(msg, *) new_line('A'),' Arnoldi approxmation of the action of the exp. propagator converged!' + write(msg, '(A,4X,A,I3)') trim(msg)//new_line('A'), 'n° of vectors:', kp + write(msg, '(A,4X,A,E10.4)') trim(msg)//new_line('A'), 'estimated error: ||err_est||_2 = ', err_est + write(msg, '(A,4X,A,E10.4)') trim(msg)//new_line('A'), 'desired tolerance: tol = ', tol + else + info = -1 + write(msg, *) new_line('A'),' Arnoldi approxmation of the action of the exp. propagator did not converge!' + write(msg, '(A,4X,A,I3)') trim(msg)//new_line('A'), 'Maximum n° of vectors reached: ', nk + 1 + write(msg, '(A,4X,A,E10.4)') trim(msg)//new_line('A'), 'Estimated error : ||err_est||_2 = ', err_est + write(msg, '(A,4X,A,E10.4)') trim(msg)//new_line('A'), 'Desired tolerance : tol = ', tol + endif + if (verbose) call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_rsp') + + return + end subroutine kexpm_vec_rsp + + subroutine kexpm_mat_rsp(C, A, B, tau, tol, info, trans, verbosity, kdim) + class(abstract_vector_rsp), intent(out) :: C(:) + !! Best Krylov approximation of \( \mathbf{C} = \exp(\tau \mathbf{A}) \mathbf{B} \). + class(abstract_linop_rsp), intent(in) :: A + !! Linear operator to be exponentiated. + class(abstract_vector_rsp), intent(in) :: B(:) + !! Input matrix on which to apply \( \exp(\tau \mathbf{A}) \). + real(sp), intent(in) :: tau + !! Time horizon for the exponentiation. + real(sp), intent(in) :: tol + !! Solution toleance based on error estimates. + integer, intent(out) :: info + !! Information flag. + logical, optional, intent(in) :: trans + !! Use transpose ? + logical, optional, intent(in) :: verbosity + !! Verbosity control. + integer, optional, intent(in) :: kdim + !! Maximum size of the Krylov subspace. + + ! ----- Internal variables ----- + integer, parameter :: kmax = 100 + integer :: i, j, k, p, kpm, kp, kpp, nk + ! Block-Arnoldi factorization. + class(abstract_vector_rsp), allocatable :: X(:) + real(sp), allocatable :: H(:, :) + ! Normalization & temporary arrays. + real(sp), allocatable :: R(:, :), E(:, :), em(:, :) + integer, allocatable :: perm(:), ptrans(:) + class(abstract_vector_rsp), allocatable :: Xwrk(:), Cwrk(:) + real(sp) :: err_est + character(len=512) :: msg + + ! Optional arguments. + logical :: transpose + logical :: verbose + integer :: nsteps + + ! Determine block size. + p = size(B) + + ! Deals with the optional args. + transpose = optval(trans, .false.) + verbose = optval(verbosity, .false.) + nsteps = optval(kdim, kmax) + nk = nsteps*p + + info = 0 + + ! Allocate arrays. + allocate(R(p, p)) ; allocate(perm(p)) ; allocate(ptrans(p)) + allocate(X(p*(nk+1)), source=B(1)) ; allocate(H(p*(nk+1), p*(nk+1))) + allocate(E(p*(nk+1), p*(nk+1))) ; allocate(em(p, p)) + + ! Scratch arrays. + allocate(Xwrk(p), source=B) ; allocate(Cwrk(p), source=B(1)) + + ! Normalize input matrix and initialize Krylov subspace. + R = 0.0_sp + call qr(Xwrk, R, perm, info) ; call apply_inverse_permutation_matrix(R, perm) + + if (norm2(abs(R)) == 0.0_sp) then + ! Input matrix is zero. + do i = 1, size(C) + call C(i)%zero() + enddo + err_est = 0.0_sp ; k = 0 ; kpp = p + else + call initialize_krylov_subspace(X, Xwrk) ; H = 0.0_sp + + expm_arnoldi: do k = 1, nk + ! Set counters. + kpm = (k-1)*p ; kp = kpm + p ; kpp = kp + p + + ! Reset working arrays. + E = 0.0_sp + do i = 1, size(Xwrk) + call Xwrk(i)%zero() + enddo + + ! Compute the k-th step of the Arnoldi factorization. + call arnoldi(A, X, H, info, kstart=k, kend=k, transpose=transpose, blksize=p) + call check_info(info, 'arnoldi', module=this_module, procedure='kexpm_mat_rsp') + + + if (info == kp) then + ! Arnoldi breakdown. Do not consider extended matrix. + kpp = kp + info = -2 + endif + + ! Compute the (dense) matrix exponential of the extended Hessenberg matrix. + call expm(E(:kpp, :kpp), tau*H(:kpp, :kpp)) + + ! Project back to original space. + do i = 1, size(Xwrk) + call Xwrk(i)%zero() + do j = 1, kpp + call Xwrk(i)%axpby(one_rsp, X(j), E(j, i)) + enddo + enddo + + do i = 1, p + call C(i)%zero() + do j = 1, p + call C(i)%axpby(one_rsp, Xwrk(j), R(j, i)) + enddo + enddo + + ! Cheap error estimate. + if (info == kp) then + ! Approximation is exact. + err_est = 0.0_sp + else + em = matmul(E(kp+1:kpp, :p), R(:p, :p)) + err_est = norm2(abs(em)) + endif + + if (err_est <= tol) exit expm_arnoldi + + enddo expm_arnoldi + endif + + if (err_est .le. tol) then + info = kpp + if (p.eq.1) then + write(msg, *) new_line('A'),' Arnoldi approxmation of the action of the exp. propagator converged!' + else + write(msg, *) new_line('A'),' Block Arnoldi approxmation of the action of the exp. propagator converged!' + endif + write(msg, '(A,2(A,I3))') trim(msg)//new_line('A'), 'n° of vectors:', k+1, ' per input vector, total:', kpp + write(msg, '(A,4X,A,E10.4)') trim(msg)//new_line('A'), 'estimated error: ||err_est||_2 = ', err_est + write(msg, '(A,4X,A,E10.4)') trim(msg)//new_line('A'), 'desired tolerance: tol = ', tol + else + info = -1 + if (p.eq.1) then + write(msg, *) 'Arnoldi approxmation of the action of the exp. propagator did not converge!' + else + write(msg, *) 'Block Arnoldi approxmation of the action of the exp. propagator did not converge!' + endif + write(msg, '(A,4X,2(A,I3))') trim(msg)//new_line('A'), 'maximum n° of vectors reached: ', nsteps+1, & + & ' per input vector, total:', kpp + write(msg, '(A,4X,A,E10.4)') trim(msg)//new_line('A'), 'estimated error: ||err_est||_2 = ', err_est + write(msg, '(A,4X,A,E10.4)') trim(msg)//new_line('A'), 'desired tolerance: tol = ', tol + endif + + if (verbose) call logger%log_message(trim(msg), module=this_module, procedure='kexpm_mat_rsp') + + return + end subroutine kexpm_mat_rsp + + subroutine kexpm_vec_variable_dt_rsp(c, A, b, tau, dt, tol, info, trans, verbosity, kdim) class(abstract_vector_rsp), intent(out) :: c !! Best approximation of \( \exp(\tau \mathbf{A}) \mathbf{b} \) in the computed Krylov subspace class(abstract_linop_rsp), intent(in) :: A @@ -235,14 +516,14 @@ subroutine kexpm_vec_rsp(c, A, b, tau, dt, tol, info, trans, verbosity, kdim) ! ----- Internal variables ----- integer, parameter :: kmax = 100 - integer :: k, km, kp, nk, istep + integer :: k, km, kp, nk, tstep, rstep ! Arnoldi factorization. class(abstract_vector_rsp), allocatable :: X(:) real(sp), allocatable :: H(:, :) ! Normaliztion & temporary arrays. real(sp), allocatable :: E(:, :) class(abstract_vector_rsp), allocatable :: b_ - real(sp) :: err_est, beta, Tend, dt_ + real(sp) :: err_est, beta, Tend, dt_, tol_, scale character(len=128) :: msg ! Optional arguments. @@ -257,6 +538,7 @@ subroutine kexpm_vec_rsp(c, A, b, tau, dt, tol, info, trans, verbosity, kdim) nk = nsteps info = 0 + scale = 1.5_sp ! Allocate arrays. allocate(X(nk+1), b_, source=b) @@ -268,11 +550,11 @@ subroutine kexpm_vec_rsp(c, A, b, tau, dt, tol, info, trans, verbosity, kdim) ! Input is zero => Output is zero. call c%zero() else - istep = 1 + tstep = 1 Tend = 0.0_sp if (verbose) then - write(msg, '(2(A,E12.6))') 'tau = ', tau, ', dt_0 = ', dt - call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_rsp') + write(msg, '(3(A,E12.6))') 'tau = ', tau, ', dt_0 = ', dt, ', tol = ', tol + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_variable_dt_rsp') end if do while ( Tend < tau ) @@ -282,19 +564,20 @@ subroutine kexpm_vec_rsp(c, A, b, tau, dt, tol, info, trans, verbosity, kdim) H = 0.0_sp ! choose time step - dt = min(tau, 2*dt) ! Reuse step size of last call (and try to increase) + dt = min(tau, scale*dt) ! Reuse step size of last call (and try to increase) if (dt > tau-Tend) then dt_ = tau-Tend ! If a small final step is needed to reach Tend, ignore it on restart if (verbose) then write(msg, '(A,E12.6)') 'Filler step, set dt_ = ', dt_ - call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_rsp') + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_variable_dt_rsp') end if else dt_ = dt end if + tol_ = dt_/tau*tol if (verbose) then - write(msg, '(A,I3,3(A,E10.4))') 'step ', istep, ': T = ', Tend, ', dt = ', dt, ', dt_ = ', dt_ - call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_rsp') + write(msg, '(A,I3,4(A,E8.2))') 'step ', tstep, ': T= ', Tend, ', dt= ', dt, ', dt_= ', dt_, ', tol_= ', tol_ + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_variable_dt_rsp') end if expm_arnoldi: do k = 1, nk @@ -305,7 +588,7 @@ subroutine kexpm_vec_rsp(c, A, b, tau, dt, tol, info, trans, verbosity, kdim) ! Compute k-th step Arnoldi factorization. call arnoldi(A, X, H, info, kstart=k, kend=k, transpose=transpose) - call check_info(info, 'arnoldi', module=this_module, procedure='kexpm_vec_rsp') + call check_info(info, 'arnoldi', module=this_module, procedure='kexpm_vec_variable_dt_rsp') ! Compute approximation. if (info == k) then @@ -327,24 +610,29 @@ subroutine kexpm_vec_rsp(c, A, b, tau, dt, tol, info, trans, verbosity, kdim) endif ! Check convergence. - if (err_est <= tol) exit expm_arnoldi + if (err_est <= tol_) exit expm_arnoldi enddo expm_arnoldi if (verbose) then - write(msg, '(30X,A,E12.6)') 'err_est = ', err_est - call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_rsp') + write(msg, '(34X,A,E12.6)') 'err_est = ', err_est + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_variable_dt_rsp') end if - if (err_est > tol) then ! decrease dt_ until error estimate is below tolerance + if (err_est > tol_) then ! decrease dt_ until error estimate is below tolerance ! NOTE: We do not need to recompute X,H! - do while (err_est > tol) + rstep = 1 + do while (err_est > tol_) + + !hfact = (tol/err_est)**(one_rsp/(kdim+1)) + !dt_ = 0.5_sp*sfact_sp*hfact*dt_ + dt_ = dt_/scale ! this is a quick hack, should choose step better - dt_ = dt_/2 ! this is a quick hack, should choose step better + tol_ = dt_/tau*tol ! update accepted tol based on new time step if (dt_ < atol_sp) then write(msg, *) 'dt < atol_sp. Cannot compute action of the matrix exponential!' - call stop_error(trim(msg), module=this_module, procedure='kexpm_vec_rsp') + call stop_error(trim(msg), module=this_module, procedure='kexpm_vec_variable_dt_rsp') end if ! Compute the (dense) matrix exponential of the extended Hessenberg matrix. @@ -354,15 +642,17 @@ subroutine kexpm_vec_rsp(c, A, b, tau, dt, tol, info, trans, verbosity, kdim) err_est = abs(E(kp, 1) * beta) if (verbose) then - write(msg, '(2(A,E12.6))') ' reduce dt_ to ', dt_, ', err_est = ', err_est - call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_rsp') + write(msg, '(6X,3(A,E10.4))') 'v dt_ = ', dt_, ', err_est = ', err_est,', tol_ = ', tol_ + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_variable_dt_rsp') end if + rstep = rstep + 1 + end do ! save dt for reuse dt = dt_ - + end if ! Project back into original space. @@ -371,23 +661,28 @@ subroutine kexpm_vec_rsp(c, A, b, tau, dt, tol, info, trans, verbosity, kdim) call linear_combination(xwrk, X(:kp), E(:kp, 1)) call b_%axpby(zero_rsp, xwrk, one_rsp*beta) end block - + ! move forward in time - Tend = Tend + dt + Tend = Tend + dt_ + + tstep = tstep + 1 - istep = istep + 1 + if (verbose) then + write(msg, '(A,E10.4,2(A,I3),A,E10.4)') ' post: T = ', Tend, ', k = ', k, '/', kdim+1, ', dt => ', dt + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_variable_dt_rsp') + end if end do ! integration ! copy result into output vector call c%axpby(zero_rsp, b_, one_rsp) - endif + endif return - end subroutine kexpm_vec_rsp + end subroutine kexpm_vec_variable_dt_rsp - subroutine kexpm_mat_rsp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) + subroutine kexpm_mat_variable_dt_rsp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) class(abstract_vector_rsp), intent(out) :: C(:) !! Best Krylov approximation of \( \mathbf{C} = \exp(\tau \mathbf{A}) \mathbf{B} \). class(abstract_linop_rsp), intent(in) :: A @@ -411,7 +706,7 @@ subroutine kexpm_mat_rsp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) ! ----- Internal variables ----- integer, parameter :: kmax = 100 - integer :: i, j, k, p, kpm, kp, kpp, nk, istep + integer :: i, j, k, p, kpm, kp, kpp, nk, tstep ! Block-Arnoldi factorization. class(abstract_vector_rsp), allocatable :: X(:) real(sp), allocatable :: H(:, :) @@ -419,7 +714,7 @@ subroutine kexpm_mat_rsp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) real(sp), allocatable :: R(:, :), E(:, :), em(:, :) integer, allocatable :: perm(:), ptrans(:) class(abstract_vector_rsp), allocatable :: Xwrk(:), Cwrk(:), B_(:) - real(sp) :: err_est, Tend, dt_ + real(sp) :: err_est, Tend, dt_, tol_, scale character(len=128) :: msg ! Optional arguments. @@ -437,6 +732,7 @@ subroutine kexpm_mat_rsp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) nk = nsteps*p info = 0 + scale = 1.5 ! Allocate arrays. allocate(perm(p), ptrans(p)) @@ -450,7 +746,7 @@ subroutine kexpm_mat_rsp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) ! Input matrix is zero. call zero_basis(C) else - istep = 1 + tstep = 1 Tend = 0.0_sp do while ( Tend < tau ) @@ -461,11 +757,20 @@ subroutine kexpm_mat_rsp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) H = 0.0_sp ! choose time step - dt = min(tau, 2*dt) ! Reuse step size of last call (and try to increase) - dt_ = max(tau-Tend, dt) ! If a small final step is needed to reach Tend, ignore it on restart + dt = min(tau, scale*dt) ! Reuse step size of last call (and try to increase) + if (dt > tau-Tend) then + dt_ = tau-Tend ! If a small final step is needed to reach Tend, ignore it on restart + if (verbose) then + write(msg, '(A,E12.6)') 'Filler step, set dt_ = ', dt_ + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_mat_variable_dt_rsp') + end if + else + dt_ = dt + end if + tol_ = dt_/tau*tol if (verbose) then - write(msg, '(A,I3,4(A,E8.2))') 'step ', istep, ': tau=', tau, ', T=', Tend, ', dt=', dt, ', dt_=', dt_ - call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_rsp') + write(msg, '(A,I3,4(A,E8.2))') 'step ', tstep, ': T= ', Tend, ', dt= ', dt, ', dt_= ', dt_, ', tol_= ', tol_ + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_mat_variable_dt_rsp') end if expm_arnoldi: do k = 1, nk @@ -477,7 +782,7 @@ subroutine kexpm_mat_rsp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) ! Compute the k-th step of the Arnoldi factorization. call arnoldi(A, X, H, info, kstart=k, kend=k, transpose=transpose, blksize=p) - call check_info(info, 'arnoldi', module=this_module, procedure='kexpm_mat_rsp') + call check_info(info, 'arnoldi', module=this_module, procedure='kexpm_mat_variable_dt_rsp') if (info == kp) then ! Arnoldi breakdown. Do not consider extended matrix. @@ -497,24 +802,28 @@ subroutine kexpm_mat_rsp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) err_est = norm2(abs(em)) endif - if (err_est <= tol) exit expm_arnoldi + if (err_est <= tol_) exit expm_arnoldi enddo expm_arnoldi if (verbose) then - write(msg, '(A,I3,A,I3,A,E8.2)') 'step ', istep, ': k=', k, ', err_est=', err_est - call logger%log_debug(trim(msg), module=this_module, procedure='kexpm_vec_rsp') + write(msg, '(34X,A,E10.4)') 'err_est = ', err_est + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_mat_variable_dt_rsp') end if - if (err_est > tol) then ! decrease dt_ until error estimate is below tolerance + if (err_est > tol_) then ! decrease dt_ until error estimate is below tolerance ! NOTE: We do not need to recompute X,H! - do while (err_est > tol) + do while (err_est > tol_) + + !hfact = (tol/err_est)**(one_rsp/kdim) + !dt_ = sfact_sp*hfact*dt_ + dt_ = dt_/scale ! this is a quick hack, should choose step better - dt_ = dt_/2 ! this is a quick hack, should choose step better + tol_ = dt_/tau*tol ! update accepted tol based on new time step if (dt_ < atol_sp) then write(msg, *) 'dt < atol_sp. Cannot compute action of the matrix exponential!' - call stop_error(trim(msg), module=this_module, procedure='kexpm_vec_rsp') + call stop_error(trim(msg), module=this_module, procedure='kexpm_mat_variable_dt_rsp') end if ! Compute the (dense) matrix exponential of the extended Hessenberg matrix. @@ -525,8 +834,8 @@ subroutine kexpm_mat_rsp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) err_est = norm2(abs(em)) if (verbose) then - write(msg, '(2(A,E8.2))') ' : dt_ =', dt_, ', err_est=', err_est - call logger%log_debug(trim(msg), module=this_module, procedure='kexpm_vec_rsp') + write(msg, '(6X,3(A,E12.6))') 'v dt_ = ', dt_, ', err_est = ', err_est, ', tol_ = ', tol_ + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_mat_variable_dt_rsp') end if end do @@ -543,9 +852,14 @@ subroutine kexpm_mat_rsp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) end block ! move forward in time - Tend = Tend + dt + Tend = Tend + dt_ + + tstep = tstep + 1 - istep = istep + 1 + if (verbose) then + write(msg, '(A,E10.4,2(A,I3),A,E10.4)') ' post: T = ', Tend, ', k = ', k, '/', kdim+1, ', dt => ', dt + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_mat_variable_dt_rsp') + end if end do ! integration @@ -554,7 +868,7 @@ subroutine kexpm_mat_rsp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) endif return - end subroutine kexpm_mat_rsp + end subroutine kexpm_mat_variable_dt_rsp subroutine k_exptA_rsp(vec_out, A, vec_in, tau, dt, info, trans) class(abstract_vector_rsp), intent(out) :: vec_out @@ -581,7 +895,7 @@ subroutine k_exptA_rsp(vec_out, A, vec_in, tau, dt, info, trans) kdim = 30 verbose = .false. - call kexpm(vec_out, A, vec_in, tau, tol, dt, info, trans=trans, verbosity=verbose, kdim=kdim) + call kexpm_var_dt(vec_out, A, vec_in, tau, tol, dt, info, trans=trans, verbosity=verbose, kdim=kdim) call check_info(info, 'kexpm', module=this_module, procedure='k_exptA_rsp') return @@ -653,7 +967,276 @@ subroutine expm_rdp(E, A, order) return end subroutine expm_rdp - subroutine kexpm_vec_rdp(c, A, b, tau, dt, tol, info, trans, verbosity, kdim) + subroutine kexpm_vec_rdp(c, A, b, tau, tol, info, trans, verbosity, kdim) + class(abstract_vector_rdp), intent(out) :: c + !! Best approximation of \( \exp(\tau \mathbf{A}) \mathbf{b} \) in the computed Krylov subspace + class(abstract_linop_rdp), intent(in) :: A + !! Linear operator to be exponentiated. + class(abstract_vector_rdp), intent(in) :: b + !! Input vector on which to apply \( \exp(\tau \mathbf{A}) \). + real(dp), intent(in) :: tau + !! Time horizon for the exponentiation. + real(dp), intent(in) :: tol + !! Solution tolerance based on error estimates. + integer, intent(out) :: info + !! Information flag. + logical, optional, intent(in) :: trans + !! Use transpose? + logical, optional, intent(in) :: verbosity + !! Verbosity control. + integer, optional, intent(in) :: kdim + !! Maximum size of the Krylov subspace. + + ! ----- Internal variables ----- + integer, parameter :: kmax = 100 + integer :: k, km, kp, nk + ! Arnoldi factorization. + class(abstract_vector_rdp), allocatable :: X(:) + real(dp), allocatable :: H(:, :) + ! Normaliztion & temporary arrays. + real(dp), allocatable :: E(:, :) + class(abstract_vector_rdp), allocatable :: Xwrk + real(dp) :: err_est, beta + character(len=512) :: msg + + ! Optional arguments. + logical :: transpose + logical :: verbose + integer :: nsteps + + ! Deals with optional args. + transpose = optval(trans, .false.) + verbose = optval(verbosity, .false.) + nsteps = optval(kdim, kmax) + nk = nsteps + + info = 0 + + ! Allocate arrays. + allocate(X(nk+1), source=b) ; allocate(H(nk+1, nk+1)) + allocate(E(nk+1, nk+1)) ; allocate(Xwrk, source=b) + + ! Normalize input vector and initialize Krylov subspace. + beta = b%norm() + if (beta == 0.0_dp) then + ! Input is zero => Output is zero. + call c%zero() + err_est = 0.0_dp + kp = 1 + else + call initialize_krylov_subspace(X) + call X(1)%add(b) ; call X(1)%scal(one_rdp / beta) + H = 0.0_dp + + expm_arnoldi: do k = 1, nk + km = k-1 ; kp = k+1 + + ! Reset work arrays. + E = 0.0_dp + + ! Compute k-th step Arnoldi factorization. + call arnoldi(A, X, H, info, kstart=k, kend=k, transpose=transpose) + call check_info(info, 'arnoldi', module=this_module, procedure='kexpm_vec_rdp') + + ! Compute approximation. + if (info == k) then + ! Arnoldi breakdown, do not consider extended matrix. + kp = k + info = -2 + endif + + ! Compute the (dense) matrix exponential of the extended Hessenberg matrix. + call expm(E(:kp, :kp), tau*H(:kp, :kp)) + + ! Project back into original space. + call linear_combination(Xwrk, X(:kp), E(:kp, 1)) + call c%axpby(zero_rdp, Xwrk, one_rdp*beta) + + ! Cheap error esimate (this actually is the magnitude of the included correction + ! and thus is very conservative). + if (info == k) then + ! Approximation is exact. + err_est = 0.0_dp + else + err_est = abs(E(kp, 1) * beta) + endif + + ! Check convergence. + if (err_est <= tol) exit expm_arnoldi + + enddo expm_arnoldi + endif + + if (err_est <= tol) then + info = kp + write(msg, *) new_line('A'),' Arnoldi approxmation of the action of the exp. propagator converged!' + write(msg, '(A,4X,A,I3)') trim(msg)//new_line('A'), 'n° of vectors:', kp + write(msg, '(A,4X,A,E10.4)') trim(msg)//new_line('A'), 'estimated error: ||err_est||_2 = ', err_est + write(msg, '(A,4X,A,E10.4)') trim(msg)//new_line('A'), 'desired tolerance: tol = ', tol + else + info = -1 + write(msg, *) new_line('A'),' Arnoldi approxmation of the action of the exp. propagator did not converge!' + write(msg, '(A,4X,A,I3)') trim(msg)//new_line('A'), 'Maximum n° of vectors reached: ', nk + 1 + write(msg, '(A,4X,A,E10.4)') trim(msg)//new_line('A'), 'Estimated error : ||err_est||_2 = ', err_est + write(msg, '(A,4X,A,E10.4)') trim(msg)//new_line('A'), 'Desired tolerance : tol = ', tol + endif + if (verbose) call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_rdp') + + return + end subroutine kexpm_vec_rdp + + subroutine kexpm_mat_rdp(C, A, B, tau, tol, info, trans, verbosity, kdim) + class(abstract_vector_rdp), intent(out) :: C(:) + !! Best Krylov approximation of \( \mathbf{C} = \exp(\tau \mathbf{A}) \mathbf{B} \). + class(abstract_linop_rdp), intent(in) :: A + !! Linear operator to be exponentiated. + class(abstract_vector_rdp), intent(in) :: B(:) + !! Input matrix on which to apply \( \exp(\tau \mathbf{A}) \). + real(dp), intent(in) :: tau + !! Time horizon for the exponentiation. + real(dp), intent(in) :: tol + !! Solution toleance based on error estimates. + integer, intent(out) :: info + !! Information flag. + logical, optional, intent(in) :: trans + !! Use transpose ? + logical, optional, intent(in) :: verbosity + !! Verbosity control. + integer, optional, intent(in) :: kdim + !! Maximum size of the Krylov subspace. + + ! ----- Internal variables ----- + integer, parameter :: kmax = 100 + integer :: i, j, k, p, kpm, kp, kpp, nk + ! Block-Arnoldi factorization. + class(abstract_vector_rdp), allocatable :: X(:) + real(dp), allocatable :: H(:, :) + ! Normalization & temporary arrays. + real(dp), allocatable :: R(:, :), E(:, :), em(:, :) + integer, allocatable :: perm(:), ptrans(:) + class(abstract_vector_rdp), allocatable :: Xwrk(:), Cwrk(:) + real(dp) :: err_est + character(len=512) :: msg + + ! Optional arguments. + logical :: transpose + logical :: verbose + integer :: nsteps + + ! Determine block size. + p = size(B) + + ! Deals with the optional args. + transpose = optval(trans, .false.) + verbose = optval(verbosity, .false.) + nsteps = optval(kdim, kmax) + nk = nsteps*p + + info = 0 + + ! Allocate arrays. + allocate(R(p, p)) ; allocate(perm(p)) ; allocate(ptrans(p)) + allocate(X(p*(nk+1)), source=B(1)) ; allocate(H(p*(nk+1), p*(nk+1))) + allocate(E(p*(nk+1), p*(nk+1))) ; allocate(em(p, p)) + + ! Scratch arrays. + allocate(Xwrk(p), source=B) ; allocate(Cwrk(p), source=B(1)) + + ! Normalize input matrix and initialize Krylov subspace. + R = 0.0_dp + call qr(Xwrk, R, perm, info) ; call apply_inverse_permutation_matrix(R, perm) + + if (norm2(abs(R)) == 0.0_dp) then + ! Input matrix is zero. + do i = 1, size(C) + call C(i)%zero() + enddo + err_est = 0.0_dp ; k = 0 ; kpp = p + else + call initialize_krylov_subspace(X, Xwrk) ; H = 0.0_dp + + expm_arnoldi: do k = 1, nk + ! Set counters. + kpm = (k-1)*p ; kp = kpm + p ; kpp = kp + p + + ! Reset working arrays. + E = 0.0_dp + do i = 1, size(Xwrk) + call Xwrk(i)%zero() + enddo + + ! Compute the k-th step of the Arnoldi factorization. + call arnoldi(A, X, H, info, kstart=k, kend=k, transpose=transpose, blksize=p) + call check_info(info, 'arnoldi', module=this_module, procedure='kexpm_mat_rdp') + + + if (info == kp) then + ! Arnoldi breakdown. Do not consider extended matrix. + kpp = kp + info = -2 + endif + + ! Compute the (dense) matrix exponential of the extended Hessenberg matrix. + call expm(E(:kpp, :kpp), tau*H(:kpp, :kpp)) + + ! Project back to original space. + do i = 1, size(Xwrk) + call Xwrk(i)%zero() + do j = 1, kpp + call Xwrk(i)%axpby(one_rdp, X(j), E(j, i)) + enddo + enddo + + do i = 1, p + call C(i)%zero() + do j = 1, p + call C(i)%axpby(one_rdp, Xwrk(j), R(j, i)) + enddo + enddo + + ! Cheap error estimate. + if (info == kp) then + ! Approximation is exact. + err_est = 0.0_dp + else + em = matmul(E(kp+1:kpp, :p), R(:p, :p)) + err_est = norm2(abs(em)) + endif + + if (err_est <= tol) exit expm_arnoldi + + enddo expm_arnoldi + endif + + if (err_est .le. tol) then + info = kpp + if (p.eq.1) then + write(msg, *) new_line('A'),' Arnoldi approxmation of the action of the exp. propagator converged!' + else + write(msg, *) new_line('A'),' Block Arnoldi approxmation of the action of the exp. propagator converged!' + endif + write(msg, '(A,2(A,I3))') trim(msg)//new_line('A'), 'n° of vectors:', k+1, ' per input vector, total:', kpp + write(msg, '(A,4X,A,E10.4)') trim(msg)//new_line('A'), 'estimated error: ||err_est||_2 = ', err_est + write(msg, '(A,4X,A,E10.4)') trim(msg)//new_line('A'), 'desired tolerance: tol = ', tol + else + info = -1 + if (p.eq.1) then + write(msg, *) 'Arnoldi approxmation of the action of the exp. propagator did not converge!' + else + write(msg, *) 'Block Arnoldi approxmation of the action of the exp. propagator did not converge!' + endif + write(msg, '(A,4X,2(A,I3))') trim(msg)//new_line('A'), 'maximum n° of vectors reached: ', nsteps+1, & + & ' per input vector, total:', kpp + write(msg, '(A,4X,A,E10.4)') trim(msg)//new_line('A'), 'estimated error: ||err_est||_2 = ', err_est + write(msg, '(A,4X,A,E10.4)') trim(msg)//new_line('A'), 'desired tolerance: tol = ', tol + endif + + if (verbose) call logger%log_message(trim(msg), module=this_module, procedure='kexpm_mat_rdp') + + return + end subroutine kexpm_mat_rdp + + subroutine kexpm_vec_variable_dt_rdp(c, A, b, tau, dt, tol, info, trans, verbosity, kdim) class(abstract_vector_rdp), intent(out) :: c !! Best approximation of \( \exp(\tau \mathbf{A}) \mathbf{b} \) in the computed Krylov subspace class(abstract_linop_rdp), intent(in) :: A @@ -677,14 +1260,14 @@ subroutine kexpm_vec_rdp(c, A, b, tau, dt, tol, info, trans, verbosity, kdim) ! ----- Internal variables ----- integer, parameter :: kmax = 100 - integer :: k, km, kp, nk, istep + integer :: k, km, kp, nk, tstep, rstep ! Arnoldi factorization. class(abstract_vector_rdp), allocatable :: X(:) real(dp), allocatable :: H(:, :) ! Normaliztion & temporary arrays. real(dp), allocatable :: E(:, :) class(abstract_vector_rdp), allocatable :: b_ - real(dp) :: err_est, beta, Tend, dt_ + real(dp) :: err_est, beta, Tend, dt_, tol_, scale character(len=128) :: msg ! Optional arguments. @@ -699,6 +1282,7 @@ subroutine kexpm_vec_rdp(c, A, b, tau, dt, tol, info, trans, verbosity, kdim) nk = nsteps info = 0 + scale = 1.5_dp ! Allocate arrays. allocate(X(nk+1), b_, source=b) @@ -710,11 +1294,11 @@ subroutine kexpm_vec_rdp(c, A, b, tau, dt, tol, info, trans, verbosity, kdim) ! Input is zero => Output is zero. call c%zero() else - istep = 1 + tstep = 1 Tend = 0.0_dp if (verbose) then - write(msg, '(2(A,E12.6))') 'tau = ', tau, ', dt_0 = ', dt - call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_rdp') + write(msg, '(3(A,E12.6))') 'tau = ', tau, ', dt_0 = ', dt, ', tol = ', tol + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_variable_dt_rdp') end if do while ( Tend < tau ) @@ -724,19 +1308,20 @@ subroutine kexpm_vec_rdp(c, A, b, tau, dt, tol, info, trans, verbosity, kdim) H = 0.0_dp ! choose time step - dt = min(tau, 2*dt) ! Reuse step size of last call (and try to increase) + dt = min(tau, scale*dt) ! Reuse step size of last call (and try to increase) if (dt > tau-Tend) then dt_ = tau-Tend ! If a small final step is needed to reach Tend, ignore it on restart if (verbose) then write(msg, '(A,E12.6)') 'Filler step, set dt_ = ', dt_ - call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_rdp') + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_variable_dt_rdp') end if else dt_ = dt end if + tol_ = dt_/tau*tol if (verbose) then - write(msg, '(A,I3,3(A,E10.4))') 'step ', istep, ': T = ', Tend, ', dt = ', dt, ', dt_ = ', dt_ - call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_rdp') + write(msg, '(A,I3,4(A,E8.2))') 'step ', tstep, ': T= ', Tend, ', dt= ', dt, ', dt_= ', dt_, ', tol_= ', tol_ + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_variable_dt_rdp') end if expm_arnoldi: do k = 1, nk @@ -747,7 +1332,7 @@ subroutine kexpm_vec_rdp(c, A, b, tau, dt, tol, info, trans, verbosity, kdim) ! Compute k-th step Arnoldi factorization. call arnoldi(A, X, H, info, kstart=k, kend=k, transpose=transpose) - call check_info(info, 'arnoldi', module=this_module, procedure='kexpm_vec_rdp') + call check_info(info, 'arnoldi', module=this_module, procedure='kexpm_vec_variable_dt_rdp') ! Compute approximation. if (info == k) then @@ -769,24 +1354,29 @@ subroutine kexpm_vec_rdp(c, A, b, tau, dt, tol, info, trans, verbosity, kdim) endif ! Check convergence. - if (err_est <= tol) exit expm_arnoldi + if (err_est <= tol_) exit expm_arnoldi enddo expm_arnoldi if (verbose) then - write(msg, '(30X,A,E12.6)') 'err_est = ', err_est - call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_rdp') + write(msg, '(34X,A,E12.6)') 'err_est = ', err_est + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_variable_dt_rdp') end if - if (err_est > tol) then ! decrease dt_ until error estimate is below tolerance + if (err_est > tol_) then ! decrease dt_ until error estimate is below tolerance ! NOTE: We do not need to recompute X,H! - do while (err_est > tol) + rstep = 1 + do while (err_est > tol_) - dt_ = dt_/2 ! this is a quick hack, should choose step better + !hfact = (tol/err_est)**(one_rdp/(kdim+1)) + !dt_ = 0.5_dp*sfact_dp*hfact*dt_ + dt_ = dt_/scale ! this is a quick hack, should choose step better + + tol_ = dt_/tau*tol ! update accepted tol based on new time step if (dt_ < atol_dp) then write(msg, *) 'dt < atol_dp. Cannot compute action of the matrix exponential!' - call stop_error(trim(msg), module=this_module, procedure='kexpm_vec_rdp') + call stop_error(trim(msg), module=this_module, procedure='kexpm_vec_variable_dt_rdp') end if ! Compute the (dense) matrix exponential of the extended Hessenberg matrix. @@ -796,15 +1386,17 @@ subroutine kexpm_vec_rdp(c, A, b, tau, dt, tol, info, trans, verbosity, kdim) err_est = abs(E(kp, 1) * beta) if (verbose) then - write(msg, '(2(A,E12.6))') ' reduce dt_ to ', dt_, ', err_est = ', err_est - call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_rdp') + write(msg, '(6X,3(A,E10.4))') 'v dt_ = ', dt_, ', err_est = ', err_est,', tol_ = ', tol_ + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_variable_dt_rdp') end if + rstep = rstep + 1 + end do ! save dt for reuse dt = dt_ - + end if ! Project back into original space. @@ -813,23 +1405,28 @@ subroutine kexpm_vec_rdp(c, A, b, tau, dt, tol, info, trans, verbosity, kdim) call linear_combination(xwrk, X(:kp), E(:kp, 1)) call b_%axpby(zero_rdp, xwrk, one_rdp*beta) end block - + ! move forward in time - Tend = Tend + dt + Tend = Tend + dt_ - istep = istep + 1 + tstep = tstep + 1 + + if (verbose) then + write(msg, '(A,E10.4,2(A,I3),A,E10.4)') ' post: T = ', Tend, ', k = ', k, '/', kdim+1, ', dt => ', dt + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_variable_dt_rdp') + end if end do ! integration ! copy result into output vector call c%axpby(zero_rdp, b_, one_rdp) - endif + endif return - end subroutine kexpm_vec_rdp + end subroutine kexpm_vec_variable_dt_rdp - subroutine kexpm_mat_rdp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) + subroutine kexpm_mat_variable_dt_rdp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) class(abstract_vector_rdp), intent(out) :: C(:) !! Best Krylov approximation of \( \mathbf{C} = \exp(\tau \mathbf{A}) \mathbf{B} \). class(abstract_linop_rdp), intent(in) :: A @@ -853,7 +1450,7 @@ subroutine kexpm_mat_rdp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) ! ----- Internal variables ----- integer, parameter :: kmax = 100 - integer :: i, j, k, p, kpm, kp, kpp, nk, istep + integer :: i, j, k, p, kpm, kp, kpp, nk, tstep ! Block-Arnoldi factorization. class(abstract_vector_rdp), allocatable :: X(:) real(dp), allocatable :: H(:, :) @@ -861,7 +1458,7 @@ subroutine kexpm_mat_rdp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) real(dp), allocatable :: R(:, :), E(:, :), em(:, :) integer, allocatable :: perm(:), ptrans(:) class(abstract_vector_rdp), allocatable :: Xwrk(:), Cwrk(:), B_(:) - real(dp) :: err_est, Tend, dt_ + real(dp) :: err_est, Tend, dt_, tol_, scale character(len=128) :: msg ! Optional arguments. @@ -879,6 +1476,7 @@ subroutine kexpm_mat_rdp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) nk = nsteps*p info = 0 + scale = 1.5 ! Allocate arrays. allocate(perm(p), ptrans(p)) @@ -892,7 +1490,7 @@ subroutine kexpm_mat_rdp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) ! Input matrix is zero. call zero_basis(C) else - istep = 1 + tstep = 1 Tend = 0.0_dp do while ( Tend < tau ) @@ -903,11 +1501,20 @@ subroutine kexpm_mat_rdp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) H = 0.0_dp ! choose time step - dt = min(tau, 2*dt) ! Reuse step size of last call (and try to increase) - dt_ = max(tau-Tend, dt) ! If a small final step is needed to reach Tend, ignore it on restart + dt = min(tau, scale*dt) ! Reuse step size of last call (and try to increase) + if (dt > tau-Tend) then + dt_ = tau-Tend ! If a small final step is needed to reach Tend, ignore it on restart + if (verbose) then + write(msg, '(A,E12.6)') 'Filler step, set dt_ = ', dt_ + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_mat_variable_dt_rdp') + end if + else + dt_ = dt + end if + tol_ = dt_/tau*tol if (verbose) then - write(msg, '(A,I3,4(A,E8.2))') 'step ', istep, ': tau=', tau, ', T=', Tend, ', dt=', dt, ', dt_=', dt_ - call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_rdp') + write(msg, '(A,I3,4(A,E8.2))') 'step ', tstep, ': T= ', Tend, ', dt= ', dt, ', dt_= ', dt_, ', tol_= ', tol_ + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_mat_variable_dt_rdp') end if expm_arnoldi: do k = 1, nk @@ -919,7 +1526,7 @@ subroutine kexpm_mat_rdp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) ! Compute the k-th step of the Arnoldi factorization. call arnoldi(A, X, H, info, kstart=k, kend=k, transpose=transpose, blksize=p) - call check_info(info, 'arnoldi', module=this_module, procedure='kexpm_mat_rdp') + call check_info(info, 'arnoldi', module=this_module, procedure='kexpm_mat_variable_dt_rdp') if (info == kp) then ! Arnoldi breakdown. Do not consider extended matrix. @@ -939,24 +1546,28 @@ subroutine kexpm_mat_rdp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) err_est = norm2(abs(em)) endif - if (err_est <= tol) exit expm_arnoldi + if (err_est <= tol_) exit expm_arnoldi enddo expm_arnoldi if (verbose) then - write(msg, '(A,I3,A,I3,A,E8.2)') 'step ', istep, ': k=', k, ', err_est=', err_est - call logger%log_debug(trim(msg), module=this_module, procedure='kexpm_vec_rdp') + write(msg, '(34X,A,E10.4)') 'err_est = ', err_est + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_mat_variable_dt_rdp') end if - if (err_est > tol) then ! decrease dt_ until error estimate is below tolerance + if (err_est > tol_) then ! decrease dt_ until error estimate is below tolerance ! NOTE: We do not need to recompute X,H! - do while (err_est > tol) + do while (err_est > tol_) - dt_ = dt_/2 ! this is a quick hack, should choose step better + !hfact = (tol/err_est)**(one_rdp/kdim) + !dt_ = sfact_dp*hfact*dt_ + dt_ = dt_/scale ! this is a quick hack, should choose step better + + tol_ = dt_/tau*tol ! update accepted tol based on new time step if (dt_ < atol_dp) then write(msg, *) 'dt < atol_dp. Cannot compute action of the matrix exponential!' - call stop_error(trim(msg), module=this_module, procedure='kexpm_vec_rdp') + call stop_error(trim(msg), module=this_module, procedure='kexpm_mat_variable_dt_rdp') end if ! Compute the (dense) matrix exponential of the extended Hessenberg matrix. @@ -967,8 +1578,8 @@ subroutine kexpm_mat_rdp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) err_est = norm2(abs(em)) if (verbose) then - write(msg, '(2(A,E8.2))') ' : dt_ =', dt_, ', err_est=', err_est - call logger%log_debug(trim(msg), module=this_module, procedure='kexpm_vec_rdp') + write(msg, '(6X,3(A,E12.6))') 'v dt_ = ', dt_, ', err_est = ', err_est, ', tol_ = ', tol_ + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_mat_variable_dt_rdp') end if end do @@ -985,9 +1596,14 @@ subroutine kexpm_mat_rdp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) end block ! move forward in time - Tend = Tend + dt + Tend = Tend + dt_ + + tstep = tstep + 1 - istep = istep + 1 + if (verbose) then + write(msg, '(A,E10.4,2(A,I3),A,E10.4)') ' post: T = ', Tend, ', k = ', k, '/', kdim+1, ', dt => ', dt + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_mat_variable_dt_rdp') + end if end do ! integration @@ -996,7 +1612,7 @@ subroutine kexpm_mat_rdp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) endif return - end subroutine kexpm_mat_rdp + end subroutine kexpm_mat_variable_dt_rdp subroutine k_exptA_rdp(vec_out, A, vec_in, tau, dt, info, trans) class(abstract_vector_rdp), intent(out) :: vec_out @@ -1023,79 +1639,348 @@ subroutine k_exptA_rdp(vec_out, A, vec_in, tau, dt, info, trans) kdim = 30 verbose = .false. - call kexpm(vec_out, A, vec_in, tau, tol, dt, info, trans=trans, verbosity=verbose, kdim=kdim) + call kexpm_var_dt(vec_out, A, vec_in, tau, tol, dt, info, trans=trans, verbosity=verbose, kdim=kdim) call check_info(info, 'kexpm', module=this_module, procedure='k_exptA_rdp') return end subroutine k_exptA_rdp - subroutine expm_csp(E, A, order) - complex(sp), intent(in) :: A(:, :) - !! Matrix to be exponentiated. - complex(sp), intent(out) :: E(:, :) - !! Output matrix E = exp(tA). - integer, optional :: order - !! Order of the Pade approximation. + subroutine expm_csp(E, A, order) + complex(sp), intent(in) :: A(:, :) + !! Matrix to be exponentiated. + complex(sp), intent(out) :: E(:, :) + !! Output matrix E = exp(tA). + integer, optional :: order + !! Order of the Pade approximation. + + !----- Internal variables ----- + complex(sp), allocatable :: A2(:, :), Q(:, :), X(:, :), invQ(:, :), wrk(:) + real(sp) :: a_norm, c + integer :: n, ee, k, s + logical :: p + integer :: p_order + + ! Deal with optional args. + p_order = optval(order, 10) + + n = size(A, 1) + + ! Allocate arrays. + allocate(A2(n, n)) ; allocate(X(n, n)) + allocate(Q(n, n)) ; allocate(invQ(n, n)) + allocate(wrk(n)) + + ! Compute the L-infinity norm. + a_norm = norml_csp(A) + + ! Determine scaling factor for the matrix. + ee = int(log2_rsp(a_norm)) + 1 + s = max(0, ee+1) + + ! Scale the input matrix & initialize polynomial. + A2 = A / 2.0_sp**s + X = A2 + + ! Initialize P & Q and add first step. + c = 0.5_sp + E = eye(n) ; E = E + c*A2 + + Q = eye(n) ; Q = Q - c*A2 + + ! Iteratively compute the Pade approximation. + p = .true. + do k = 2, p_order + c = c*(p_order - k + 1) / (k * (2*p_order - k + 1)) + X = matmul(A2, X) + E = E + c*X + if (p) then + Q = Q + c*X + else + Q = Q - c*X + endif + p = .not. p + enddo + + invQ = Q ; call inv(invQ) + E = matmul(invQ, E) + + do k = 1, s + E = matmul(E, E) + enddo + + return + end subroutine expm_csp + + subroutine kexpm_vec_csp(c, A, b, tau, tol, info, trans, verbosity, kdim) + class(abstract_vector_csp), intent(out) :: c + !! Best approximation of \( \exp(\tau \mathbf{A}) \mathbf{b} \) in the computed Krylov subspace + class(abstract_linop_csp), intent(in) :: A + !! Linear operator to be exponentiated. + class(abstract_vector_csp), intent(in) :: b + !! Input vector on which to apply \( \exp(\tau \mathbf{A}) \). + real(sp), intent(in) :: tau + !! Time horizon for the exponentiation. + real(sp), intent(in) :: tol + !! Solution tolerance based on error estimates. + integer, intent(out) :: info + !! Information flag. + logical, optional, intent(in) :: trans + !! Use transpose? + logical, optional, intent(in) :: verbosity + !! Verbosity control. + integer, optional, intent(in) :: kdim + !! Maximum size of the Krylov subspace. + + ! ----- Internal variables ----- + integer, parameter :: kmax = 100 + integer :: k, km, kp, nk + ! Arnoldi factorization. + class(abstract_vector_csp), allocatable :: X(:) + complex(sp), allocatable :: H(:, :) + ! Normaliztion & temporary arrays. + complex(sp), allocatable :: E(:, :) + class(abstract_vector_csp), allocatable :: Xwrk + real(sp) :: err_est, beta + character(len=512) :: msg + + ! Optional arguments. + logical :: transpose + logical :: verbose + integer :: nsteps + + ! Deals with optional args. + transpose = optval(trans, .false.) + verbose = optval(verbosity, .false.) + nsteps = optval(kdim, kmax) + nk = nsteps + + info = 0 + + ! Allocate arrays. + allocate(X(nk+1), source=b) ; allocate(H(nk+1, nk+1)) + allocate(E(nk+1, nk+1)) ; allocate(Xwrk, source=b) + + ! Normalize input vector and initialize Krylov subspace. + beta = b%norm() + if (beta == 0.0_sp) then + ! Input is zero => Output is zero. + call c%zero() + err_est = 0.0_sp + kp = 1 + else + call initialize_krylov_subspace(X) + call X(1)%add(b) ; call X(1)%scal(one_csp / beta) + H = 0.0_sp + + expm_arnoldi: do k = 1, nk + km = k-1 ; kp = k+1 + + ! Reset work arrays. + E = 0.0_sp + + ! Compute k-th step Arnoldi factorization. + call arnoldi(A, X, H, info, kstart=k, kend=k, transpose=transpose) + call check_info(info, 'arnoldi', module=this_module, procedure='kexpm_vec_csp') + + ! Compute approximation. + if (info == k) then + ! Arnoldi breakdown, do not consider extended matrix. + kp = k + info = -2 + endif + + ! Compute the (dense) matrix exponential of the extended Hessenberg matrix. + call expm(E(:kp, :kp), tau*H(:kp, :kp)) + + ! Project back into original space. + call linear_combination(Xwrk, X(:kp), E(:kp, 1)) + call c%axpby(zero_csp, Xwrk, one_csp*beta) + + ! Cheap error esimate (this actually is the magnitude of the included correction + ! and thus is very conservative). + if (info == k) then + ! Approximation is exact. + err_est = 0.0_sp + else + err_est = abs(E(kp, 1) * beta) + endif + + ! Check convergence. + if (err_est <= tol) exit expm_arnoldi + + enddo expm_arnoldi + endif + + if (err_est <= tol) then + info = kp + write(msg, *) new_line('A'),' Arnoldi approxmation of the action of the exp. propagator converged!' + write(msg, '(A,4X,A,I3)') trim(msg)//new_line('A'), 'n° of vectors:', kp + write(msg, '(A,4X,A,E10.4)') trim(msg)//new_line('A'), 'estimated error: ||err_est||_2 = ', err_est + write(msg, '(A,4X,A,E10.4)') trim(msg)//new_line('A'), 'desired tolerance: tol = ', tol + else + info = -1 + write(msg, *) new_line('A'),' Arnoldi approxmation of the action of the exp. propagator did not converge!' + write(msg, '(A,4X,A,I3)') trim(msg)//new_line('A'), 'Maximum n° of vectors reached: ', nk + 1 + write(msg, '(A,4X,A,E10.4)') trim(msg)//new_line('A'), 'Estimated error : ||err_est||_2 = ', err_est + write(msg, '(A,4X,A,E10.4)') trim(msg)//new_line('A'), 'Desired tolerance : tol = ', tol + endif + if (verbose) call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_csp') + + return + end subroutine kexpm_vec_csp + + subroutine kexpm_mat_csp(C, A, B, tau, tol, info, trans, verbosity, kdim) + class(abstract_vector_csp), intent(out) :: C(:) + !! Best Krylov approximation of \( \mathbf{C} = \exp(\tau \mathbf{A}) \mathbf{B} \). + class(abstract_linop_csp), intent(in) :: A + !! Linear operator to be exponentiated. + class(abstract_vector_csp), intent(in) :: B(:) + !! Input matrix on which to apply \( \exp(\tau \mathbf{A}) \). + real(sp), intent(in) :: tau + !! Time horizon for the exponentiation. + real(sp), intent(in) :: tol + !! Solution toleance based on error estimates. + integer, intent(out) :: info + !! Information flag. + logical, optional, intent(in) :: trans + !! Use transpose ? + logical, optional, intent(in) :: verbosity + !! Verbosity control. + integer, optional, intent(in) :: kdim + !! Maximum size of the Krylov subspace. + + ! ----- Internal variables ----- + integer, parameter :: kmax = 100 + integer :: i, j, k, p, kpm, kp, kpp, nk + ! Block-Arnoldi factorization. + class(abstract_vector_csp), allocatable :: X(:) + complex(sp), allocatable :: H(:, :) + ! Normalization & temporary arrays. + complex(sp), allocatable :: R(:, :), E(:, :), em(:, :) + integer, allocatable :: perm(:), ptrans(:) + class(abstract_vector_csp), allocatable :: Xwrk(:), Cwrk(:) + real(sp) :: err_est + character(len=512) :: msg - !----- Internal variables ----- - complex(sp), allocatable :: A2(:, :), Q(:, :), X(:, :), invQ(:, :), wrk(:) - real(sp) :: a_norm, c - integer :: n, ee, k, s - logical :: p - integer :: p_order + ! Optional arguments. + logical :: transpose + logical :: verbose + integer :: nsteps - ! Deal with optional args. - p_order = optval(order, 10) + ! Determine block size. + p = size(B) - n = size(A, 1) + ! Deals with the optional args. + transpose = optval(trans, .false.) + verbose = optval(verbosity, .false.) + nsteps = optval(kdim, kmax) + nk = nsteps*p + + info = 0 ! Allocate arrays. - allocate(A2(n, n)) ; allocate(X(n, n)) - allocate(Q(n, n)) ; allocate(invQ(n, n)) - allocate(wrk(n)) - - ! Compute the L-infinity norm. - a_norm = norml_csp(A) - - ! Determine scaling factor for the matrix. - ee = int(log2_rsp(a_norm)) + 1 - s = max(0, ee+1) + allocate(R(p, p)) ; allocate(perm(p)) ; allocate(ptrans(p)) + allocate(X(p*(nk+1)), source=B(1)) ; allocate(H(p*(nk+1), p*(nk+1))) + allocate(E(p*(nk+1), p*(nk+1))) ; allocate(em(p, p)) - ! Scale the input matrix & initialize polynomial. - A2 = A / 2.0_sp**s - X = A2 + ! Scratch arrays. + allocate(Xwrk(p), source=B) ; allocate(Cwrk(p), source=B(1)) - ! Initialize P & Q and add first step. - c = 0.5_sp - E = eye(n) ; E = E + c*A2 + ! Normalize input matrix and initialize Krylov subspace. + R = 0.0_sp + call qr(Xwrk, R, perm, info) ; call apply_inverse_permutation_matrix(R, perm) - Q = eye(n) ; Q = Q - c*A2 + if (norm2(abs(R)) == 0.0_sp) then + ! Input matrix is zero. + do i = 1, size(C) + call C(i)%zero() + enddo + err_est = 0.0_sp ; k = 0 ; kpp = p + else + call initialize_krylov_subspace(X, Xwrk) ; H = 0.0_sp + + expm_arnoldi: do k = 1, nk + ! Set counters. + kpm = (k-1)*p ; kp = kpm + p ; kpp = kp + p + + ! Reset working arrays. + E = 0.0_sp + do i = 1, size(Xwrk) + call Xwrk(i)%zero() + enddo + + ! Compute the k-th step of the Arnoldi factorization. + call arnoldi(A, X, H, info, kstart=k, kend=k, transpose=transpose, blksize=p) + call check_info(info, 'arnoldi', module=this_module, procedure='kexpm_mat_csp') + + + if (info == kp) then + ! Arnoldi breakdown. Do not consider extended matrix. + kpp = kp + info = -2 + endif + + ! Compute the (dense) matrix exponential of the extended Hessenberg matrix. + call expm(E(:kpp, :kpp), tau*H(:kpp, :kpp)) + + ! Project back to original space. + do i = 1, size(Xwrk) + call Xwrk(i)%zero() + do j = 1, kpp + call Xwrk(i)%axpby(one_csp, X(j), E(j, i)) + enddo + enddo + + do i = 1, p + call C(i)%zero() + do j = 1, p + call C(i)%axpby(one_csp, Xwrk(j), R(j, i)) + enddo + enddo + + ! Cheap error estimate. + if (info == kp) then + ! Approximation is exact. + err_est = 0.0_sp + else + em = matmul(E(kp+1:kpp, :p), R(:p, :p)) + err_est = norm2(abs(em)) + endif + + if (err_est <= tol) exit expm_arnoldi + + enddo expm_arnoldi + endif - ! Iteratively compute the Pade approximation. - p = .true. - do k = 2, p_order - c = c*(p_order - k + 1) / (k * (2*p_order - k + 1)) - X = matmul(A2, X) - E = E + c*X - if (p) then - Q = Q + c*X + if (err_est .le. tol) then + info = kpp + if (p.eq.1) then + write(msg, *) new_line('A'),' Arnoldi approxmation of the action of the exp. propagator converged!' else - Q = Q - c*X + write(msg, *) new_line('A'),' Block Arnoldi approxmation of the action of the exp. propagator converged!' + endif + write(msg, '(A,2(A,I3))') trim(msg)//new_line('A'), 'n° of vectors:', k+1, ' per input vector, total:', kpp + write(msg, '(A,4X,A,E10.4)') trim(msg)//new_line('A'), 'estimated error: ||err_est||_2 = ', err_est + write(msg, '(A,4X,A,E10.4)') trim(msg)//new_line('A'), 'desired tolerance: tol = ', tol + else + info = -1 + if (p.eq.1) then + write(msg, *) 'Arnoldi approxmation of the action of the exp. propagator did not converge!' + else + write(msg, *) 'Block Arnoldi approxmation of the action of the exp. propagator did not converge!' endif - p = .not. p - enddo - - invQ = Q ; call inv(invQ) - E = matmul(invQ, E) + write(msg, '(A,4X,2(A,I3))') trim(msg)//new_line('A'), 'maximum n° of vectors reached: ', nsteps+1, & + & ' per input vector, total:', kpp + write(msg, '(A,4X,A,E10.4)') trim(msg)//new_line('A'), 'estimated error: ||err_est||_2 = ', err_est + write(msg, '(A,4X,A,E10.4)') trim(msg)//new_line('A'), 'desired tolerance: tol = ', tol + endif - do k = 1, s - E = matmul(E, E) - enddo + if (verbose) call logger%log_message(trim(msg), module=this_module, procedure='kexpm_mat_csp') return - end subroutine expm_csp + end subroutine kexpm_mat_csp - subroutine kexpm_vec_csp(c, A, b, tau, dt, tol, info, trans, verbosity, kdim) + subroutine kexpm_vec_variable_dt_csp(c, A, b, tau, dt, tol, info, trans, verbosity, kdim) class(abstract_vector_csp), intent(out) :: c !! Best approximation of \( \exp(\tau \mathbf{A}) \mathbf{b} \) in the computed Krylov subspace class(abstract_linop_csp), intent(in) :: A @@ -1119,14 +2004,14 @@ subroutine kexpm_vec_csp(c, A, b, tau, dt, tol, info, trans, verbosity, kdim) ! ----- Internal variables ----- integer, parameter :: kmax = 100 - integer :: k, km, kp, nk, istep + integer :: k, km, kp, nk, tstep, rstep ! Arnoldi factorization. class(abstract_vector_csp), allocatable :: X(:) complex(sp), allocatable :: H(:, :) ! Normaliztion & temporary arrays. complex(sp), allocatable :: E(:, :) class(abstract_vector_csp), allocatable :: b_ - real(sp) :: err_est, beta, Tend, dt_ + real(sp) :: err_est, beta, Tend, dt_, tol_, scale character(len=128) :: msg ! Optional arguments. @@ -1141,6 +2026,7 @@ subroutine kexpm_vec_csp(c, A, b, tau, dt, tol, info, trans, verbosity, kdim) nk = nsteps info = 0 + scale = 1.5_sp ! Allocate arrays. allocate(X(nk+1), b_, source=b) @@ -1152,11 +2038,11 @@ subroutine kexpm_vec_csp(c, A, b, tau, dt, tol, info, trans, verbosity, kdim) ! Input is zero => Output is zero. call c%zero() else - istep = 1 + tstep = 1 Tend = 0.0_sp if (verbose) then - write(msg, '(2(A,E12.6))') 'tau = ', tau, ', dt_0 = ', dt - call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_csp') + write(msg, '(3(A,E12.6))') 'tau = ', tau, ', dt_0 = ', dt, ', tol = ', tol + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_variable_dt_csp') end if do while ( Tend < tau ) @@ -1166,19 +2052,20 @@ subroutine kexpm_vec_csp(c, A, b, tau, dt, tol, info, trans, verbosity, kdim) H = 0.0_sp ! choose time step - dt = min(tau, 2*dt) ! Reuse step size of last call (and try to increase) + dt = min(tau, scale*dt) ! Reuse step size of last call (and try to increase) if (dt > tau-Tend) then dt_ = tau-Tend ! If a small final step is needed to reach Tend, ignore it on restart if (verbose) then write(msg, '(A,E12.6)') 'Filler step, set dt_ = ', dt_ - call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_csp') + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_variable_dt_csp') end if else dt_ = dt end if + tol_ = dt_/tau*tol if (verbose) then - write(msg, '(A,I3,3(A,E10.4))') 'step ', istep, ': T = ', Tend, ', dt = ', dt, ', dt_ = ', dt_ - call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_csp') + write(msg, '(A,I3,4(A,E8.2))') 'step ', tstep, ': T= ', Tend, ', dt= ', dt, ', dt_= ', dt_, ', tol_= ', tol_ + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_variable_dt_csp') end if expm_arnoldi: do k = 1, nk @@ -1189,7 +2076,7 @@ subroutine kexpm_vec_csp(c, A, b, tau, dt, tol, info, trans, verbosity, kdim) ! Compute k-th step Arnoldi factorization. call arnoldi(A, X, H, info, kstart=k, kend=k, transpose=transpose) - call check_info(info, 'arnoldi', module=this_module, procedure='kexpm_vec_csp') + call check_info(info, 'arnoldi', module=this_module, procedure='kexpm_vec_variable_dt_csp') ! Compute approximation. if (info == k) then @@ -1211,24 +2098,29 @@ subroutine kexpm_vec_csp(c, A, b, tau, dt, tol, info, trans, verbosity, kdim) endif ! Check convergence. - if (err_est <= tol) exit expm_arnoldi + if (err_est <= tol_) exit expm_arnoldi enddo expm_arnoldi if (verbose) then - write(msg, '(30X,A,E12.6)') 'err_est = ', err_est - call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_csp') + write(msg, '(34X,A,E12.6)') 'err_est = ', err_est + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_variable_dt_csp') end if - if (err_est > tol) then ! decrease dt_ until error estimate is below tolerance + if (err_est > tol_) then ! decrease dt_ until error estimate is below tolerance ! NOTE: We do not need to recompute X,H! - do while (err_est > tol) + rstep = 1 + do while (err_est > tol_) - dt_ = dt_/2 ! this is a quick hack, should choose step better + !hfact = (tol/err_est)**(one_csp/(kdim+1)) + !dt_ = 0.5_sp*sfact_sp*hfact*dt_ + dt_ = dt_/scale ! this is a quick hack, should choose step better + + tol_ = dt_/tau*tol ! update accepted tol based on new time step if (dt_ < atol_sp) then write(msg, *) 'dt < atol_sp. Cannot compute action of the matrix exponential!' - call stop_error(trim(msg), module=this_module, procedure='kexpm_vec_csp') + call stop_error(trim(msg), module=this_module, procedure='kexpm_vec_variable_dt_csp') end if ! Compute the (dense) matrix exponential of the extended Hessenberg matrix. @@ -1238,15 +2130,17 @@ subroutine kexpm_vec_csp(c, A, b, tau, dt, tol, info, trans, verbosity, kdim) err_est = abs(E(kp, 1) * beta) if (verbose) then - write(msg, '(2(A,E12.6))') ' reduce dt_ to ', dt_, ', err_est = ', err_est - call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_csp') + write(msg, '(6X,3(A,E10.4))') 'v dt_ = ', dt_, ', err_est = ', err_est,', tol_ = ', tol_ + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_variable_dt_csp') end if + rstep = rstep + 1 + end do ! save dt for reuse dt = dt_ - + end if ! Project back into original space. @@ -1255,23 +2149,28 @@ subroutine kexpm_vec_csp(c, A, b, tau, dt, tol, info, trans, verbosity, kdim) call linear_combination(xwrk, X(:kp), E(:kp, 1)) call b_%axpby(zero_csp, xwrk, one_csp*beta) end block - + ! move forward in time - Tend = Tend + dt + Tend = Tend + dt_ + + tstep = tstep + 1 - istep = istep + 1 + if (verbose) then + write(msg, '(A,E10.4,2(A,I3),A,E10.4)') ' post: T = ', Tend, ', k = ', k, '/', kdim+1, ', dt => ', dt + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_variable_dt_csp') + end if end do ! integration ! copy result into output vector call c%axpby(zero_csp, b_, one_csp) - endif + endif return - end subroutine kexpm_vec_csp + end subroutine kexpm_vec_variable_dt_csp - subroutine kexpm_mat_csp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) + subroutine kexpm_mat_variable_dt_csp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) class(abstract_vector_csp), intent(out) :: C(:) !! Best Krylov approximation of \( \mathbf{C} = \exp(\tau \mathbf{A}) \mathbf{B} \). class(abstract_linop_csp), intent(in) :: A @@ -1295,7 +2194,7 @@ subroutine kexpm_mat_csp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) ! ----- Internal variables ----- integer, parameter :: kmax = 100 - integer :: i, j, k, p, kpm, kp, kpp, nk, istep + integer :: i, j, k, p, kpm, kp, kpp, nk, tstep ! Block-Arnoldi factorization. class(abstract_vector_csp), allocatable :: X(:) complex(sp), allocatable :: H(:, :) @@ -1303,7 +2202,7 @@ subroutine kexpm_mat_csp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) complex(sp), allocatable :: R(:, :), E(:, :), em(:, :) integer, allocatable :: perm(:), ptrans(:) class(abstract_vector_csp), allocatable :: Xwrk(:), Cwrk(:), B_(:) - real(sp) :: err_est, Tend, dt_ + real(sp) :: err_est, Tend, dt_, tol_, scale character(len=128) :: msg ! Optional arguments. @@ -1321,6 +2220,7 @@ subroutine kexpm_mat_csp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) nk = nsteps*p info = 0 + scale = 1.5 ! Allocate arrays. allocate(perm(p), ptrans(p)) @@ -1334,7 +2234,7 @@ subroutine kexpm_mat_csp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) ! Input matrix is zero. call zero_basis(C) else - istep = 1 + tstep = 1 Tend = 0.0_sp do while ( Tend < tau ) @@ -1345,11 +2245,20 @@ subroutine kexpm_mat_csp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) H = 0.0_sp ! choose time step - dt = min(tau, 2*dt) ! Reuse step size of last call (and try to increase) - dt_ = max(tau-Tend, dt) ! If a small final step is needed to reach Tend, ignore it on restart + dt = min(tau, scale*dt) ! Reuse step size of last call (and try to increase) + if (dt > tau-Tend) then + dt_ = tau-Tend ! If a small final step is needed to reach Tend, ignore it on restart + if (verbose) then + write(msg, '(A,E12.6)') 'Filler step, set dt_ = ', dt_ + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_mat_variable_dt_csp') + end if + else + dt_ = dt + end if + tol_ = dt_/tau*tol if (verbose) then - write(msg, '(A,I3,4(A,E8.2))') 'step ', istep, ': tau=', tau, ', T=', Tend, ', dt=', dt, ', dt_=', dt_ - call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_csp') + write(msg, '(A,I3,4(A,E8.2))') 'step ', tstep, ': T= ', Tend, ', dt= ', dt, ', dt_= ', dt_, ', tol_= ', tol_ + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_mat_variable_dt_csp') end if expm_arnoldi: do k = 1, nk @@ -1361,7 +2270,7 @@ subroutine kexpm_mat_csp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) ! Compute the k-th step of the Arnoldi factorization. call arnoldi(A, X, H, info, kstart=k, kend=k, transpose=transpose, blksize=p) - call check_info(info, 'arnoldi', module=this_module, procedure='kexpm_mat_csp') + call check_info(info, 'arnoldi', module=this_module, procedure='kexpm_mat_variable_dt_csp') if (info == kp) then ! Arnoldi breakdown. Do not consider extended matrix. @@ -1381,24 +2290,28 @@ subroutine kexpm_mat_csp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) err_est = norm2(abs(em)) endif - if (err_est <= tol) exit expm_arnoldi + if (err_est <= tol_) exit expm_arnoldi enddo expm_arnoldi if (verbose) then - write(msg, '(A,I3,A,I3,A,E8.2)') 'step ', istep, ': k=', k, ', err_est=', err_est - call logger%log_debug(trim(msg), module=this_module, procedure='kexpm_vec_csp') + write(msg, '(34X,A,E10.4)') 'err_est = ', err_est + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_mat_variable_dt_csp') end if - if (err_est > tol) then ! decrease dt_ until error estimate is below tolerance + if (err_est > tol_) then ! decrease dt_ until error estimate is below tolerance ! NOTE: We do not need to recompute X,H! - do while (err_est > tol) + do while (err_est > tol_) - dt_ = dt_/2 ! this is a quick hack, should choose step better + !hfact = (tol/err_est)**(one_csp/kdim) + !dt_ = sfact_sp*hfact*dt_ + dt_ = dt_/scale ! this is a quick hack, should choose step better + + tol_ = dt_/tau*tol ! update accepted tol based on new time step if (dt_ < atol_sp) then write(msg, *) 'dt < atol_sp. Cannot compute action of the matrix exponential!' - call stop_error(trim(msg), module=this_module, procedure='kexpm_vec_csp') + call stop_error(trim(msg), module=this_module, procedure='kexpm_mat_variable_dt_csp') end if ! Compute the (dense) matrix exponential of the extended Hessenberg matrix. @@ -1409,8 +2322,8 @@ subroutine kexpm_mat_csp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) err_est = norm2(abs(em)) if (verbose) then - write(msg, '(2(A,E8.2))') ' : dt_ =', dt_, ', err_est=', err_est - call logger%log_debug(trim(msg), module=this_module, procedure='kexpm_vec_csp') + write(msg, '(6X,3(A,E12.6))') 'v dt_ = ', dt_, ', err_est = ', err_est, ', tol_ = ', tol_ + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_mat_variable_dt_csp') end if end do @@ -1427,9 +2340,14 @@ subroutine kexpm_mat_csp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) end block ! move forward in time - Tend = Tend + dt + Tend = Tend + dt_ + + tstep = tstep + 1 - istep = istep + 1 + if (verbose) then + write(msg, '(A,E10.4,2(A,I3),A,E10.4)') ' post: T = ', Tend, ', k = ', k, '/', kdim+1, ', dt => ', dt + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_mat_variable_dt_csp') + end if end do ! integration @@ -1438,7 +2356,7 @@ subroutine kexpm_mat_csp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) endif return - end subroutine kexpm_mat_csp + end subroutine kexpm_mat_variable_dt_csp subroutine k_exptA_csp(vec_out, A, vec_in, tau, dt, info, trans) class(abstract_vector_csp), intent(out) :: vec_out @@ -1465,7 +2383,7 @@ subroutine k_exptA_csp(vec_out, A, vec_in, tau, dt, info, trans) kdim = 30 verbose = .false. - call kexpm(vec_out, A, vec_in, tau, tol, dt, info, trans=trans, verbosity=verbose, kdim=kdim) + call kexpm_var_dt(vec_out, A, vec_in, tau, tol, dt, info, trans=trans, verbosity=verbose, kdim=kdim) call check_info(info, 'kexpm', module=this_module, procedure='k_exptA_csp') return @@ -1537,7 +2455,276 @@ subroutine expm_cdp(E, A, order) return end subroutine expm_cdp - subroutine kexpm_vec_cdp(c, A, b, tau, dt, tol, info, trans, verbosity, kdim) + subroutine kexpm_vec_cdp(c, A, b, tau, tol, info, trans, verbosity, kdim) + class(abstract_vector_cdp), intent(out) :: c + !! Best approximation of \( \exp(\tau \mathbf{A}) \mathbf{b} \) in the computed Krylov subspace + class(abstract_linop_cdp), intent(in) :: A + !! Linear operator to be exponentiated. + class(abstract_vector_cdp), intent(in) :: b + !! Input vector on which to apply \( \exp(\tau \mathbf{A}) \). + real(dp), intent(in) :: tau + !! Time horizon for the exponentiation. + real(dp), intent(in) :: tol + !! Solution tolerance based on error estimates. + integer, intent(out) :: info + !! Information flag. + logical, optional, intent(in) :: trans + !! Use transpose? + logical, optional, intent(in) :: verbosity + !! Verbosity control. + integer, optional, intent(in) :: kdim + !! Maximum size of the Krylov subspace. + + ! ----- Internal variables ----- + integer, parameter :: kmax = 100 + integer :: k, km, kp, nk + ! Arnoldi factorization. + class(abstract_vector_cdp), allocatable :: X(:) + complex(dp), allocatable :: H(:, :) + ! Normaliztion & temporary arrays. + complex(dp), allocatable :: E(:, :) + class(abstract_vector_cdp), allocatable :: Xwrk + real(dp) :: err_est, beta + character(len=512) :: msg + + ! Optional arguments. + logical :: transpose + logical :: verbose + integer :: nsteps + + ! Deals with optional args. + transpose = optval(trans, .false.) + verbose = optval(verbosity, .false.) + nsteps = optval(kdim, kmax) + nk = nsteps + + info = 0 + + ! Allocate arrays. + allocate(X(nk+1), source=b) ; allocate(H(nk+1, nk+1)) + allocate(E(nk+1, nk+1)) ; allocate(Xwrk, source=b) + + ! Normalize input vector and initialize Krylov subspace. + beta = b%norm() + if (beta == 0.0_dp) then + ! Input is zero => Output is zero. + call c%zero() + err_est = 0.0_dp + kp = 1 + else + call initialize_krylov_subspace(X) + call X(1)%add(b) ; call X(1)%scal(one_cdp / beta) + H = 0.0_dp + + expm_arnoldi: do k = 1, nk + km = k-1 ; kp = k+1 + + ! Reset work arrays. + E = 0.0_dp + + ! Compute k-th step Arnoldi factorization. + call arnoldi(A, X, H, info, kstart=k, kend=k, transpose=transpose) + call check_info(info, 'arnoldi', module=this_module, procedure='kexpm_vec_cdp') + + ! Compute approximation. + if (info == k) then + ! Arnoldi breakdown, do not consider extended matrix. + kp = k + info = -2 + endif + + ! Compute the (dense) matrix exponential of the extended Hessenberg matrix. + call expm(E(:kp, :kp), tau*H(:kp, :kp)) + + ! Project back into original space. + call linear_combination(Xwrk, X(:kp), E(:kp, 1)) + call c%axpby(zero_cdp, Xwrk, one_cdp*beta) + + ! Cheap error esimate (this actually is the magnitude of the included correction + ! and thus is very conservative). + if (info == k) then + ! Approximation is exact. + err_est = 0.0_dp + else + err_est = abs(E(kp, 1) * beta) + endif + + ! Check convergence. + if (err_est <= tol) exit expm_arnoldi + + enddo expm_arnoldi + endif + + if (err_est <= tol) then + info = kp + write(msg, *) new_line('A'),' Arnoldi approxmation of the action of the exp. propagator converged!' + write(msg, '(A,4X,A,I3)') trim(msg)//new_line('A'), 'n° of vectors:', kp + write(msg, '(A,4X,A,E10.4)') trim(msg)//new_line('A'), 'estimated error: ||err_est||_2 = ', err_est + write(msg, '(A,4X,A,E10.4)') trim(msg)//new_line('A'), 'desired tolerance: tol = ', tol + else + info = -1 + write(msg, *) new_line('A'),' Arnoldi approxmation of the action of the exp. propagator did not converge!' + write(msg, '(A,4X,A,I3)') trim(msg)//new_line('A'), 'Maximum n° of vectors reached: ', nk + 1 + write(msg, '(A,4X,A,E10.4)') trim(msg)//new_line('A'), 'Estimated error : ||err_est||_2 = ', err_est + write(msg, '(A,4X,A,E10.4)') trim(msg)//new_line('A'), 'Desired tolerance : tol = ', tol + endif + if (verbose) call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_cdp') + + return + end subroutine kexpm_vec_cdp + + subroutine kexpm_mat_cdp(C, A, B, tau, tol, info, trans, verbosity, kdim) + class(abstract_vector_cdp), intent(out) :: C(:) + !! Best Krylov approximation of \( \mathbf{C} = \exp(\tau \mathbf{A}) \mathbf{B} \). + class(abstract_linop_cdp), intent(in) :: A + !! Linear operator to be exponentiated. + class(abstract_vector_cdp), intent(in) :: B(:) + !! Input matrix on which to apply \( \exp(\tau \mathbf{A}) \). + real(dp), intent(in) :: tau + !! Time horizon for the exponentiation. + real(dp), intent(in) :: tol + !! Solution toleance based on error estimates. + integer, intent(out) :: info + !! Information flag. + logical, optional, intent(in) :: trans + !! Use transpose ? + logical, optional, intent(in) :: verbosity + !! Verbosity control. + integer, optional, intent(in) :: kdim + !! Maximum size of the Krylov subspace. + + ! ----- Internal variables ----- + integer, parameter :: kmax = 100 + integer :: i, j, k, p, kpm, kp, kpp, nk + ! Block-Arnoldi factorization. + class(abstract_vector_cdp), allocatable :: X(:) + complex(dp), allocatable :: H(:, :) + ! Normalization & temporary arrays. + complex(dp), allocatable :: R(:, :), E(:, :), em(:, :) + integer, allocatable :: perm(:), ptrans(:) + class(abstract_vector_cdp), allocatable :: Xwrk(:), Cwrk(:) + real(dp) :: err_est + character(len=512) :: msg + + ! Optional arguments. + logical :: transpose + logical :: verbose + integer :: nsteps + + ! Determine block size. + p = size(B) + + ! Deals with the optional args. + transpose = optval(trans, .false.) + verbose = optval(verbosity, .false.) + nsteps = optval(kdim, kmax) + nk = nsteps*p + + info = 0 + + ! Allocate arrays. + allocate(R(p, p)) ; allocate(perm(p)) ; allocate(ptrans(p)) + allocate(X(p*(nk+1)), source=B(1)) ; allocate(H(p*(nk+1), p*(nk+1))) + allocate(E(p*(nk+1), p*(nk+1))) ; allocate(em(p, p)) + + ! Scratch arrays. + allocate(Xwrk(p), source=B) ; allocate(Cwrk(p), source=B(1)) + + ! Normalize input matrix and initialize Krylov subspace. + R = 0.0_dp + call qr(Xwrk, R, perm, info) ; call apply_inverse_permutation_matrix(R, perm) + + if (norm2(abs(R)) == 0.0_dp) then + ! Input matrix is zero. + do i = 1, size(C) + call C(i)%zero() + enddo + err_est = 0.0_dp ; k = 0 ; kpp = p + else + call initialize_krylov_subspace(X, Xwrk) ; H = 0.0_dp + + expm_arnoldi: do k = 1, nk + ! Set counters. + kpm = (k-1)*p ; kp = kpm + p ; kpp = kp + p + + ! Reset working arrays. + E = 0.0_dp + do i = 1, size(Xwrk) + call Xwrk(i)%zero() + enddo + + ! Compute the k-th step of the Arnoldi factorization. + call arnoldi(A, X, H, info, kstart=k, kend=k, transpose=transpose, blksize=p) + call check_info(info, 'arnoldi', module=this_module, procedure='kexpm_mat_cdp') + + + if (info == kp) then + ! Arnoldi breakdown. Do not consider extended matrix. + kpp = kp + info = -2 + endif + + ! Compute the (dense) matrix exponential of the extended Hessenberg matrix. + call expm(E(:kpp, :kpp), tau*H(:kpp, :kpp)) + + ! Project back to original space. + do i = 1, size(Xwrk) + call Xwrk(i)%zero() + do j = 1, kpp + call Xwrk(i)%axpby(one_cdp, X(j), E(j, i)) + enddo + enddo + + do i = 1, p + call C(i)%zero() + do j = 1, p + call C(i)%axpby(one_cdp, Xwrk(j), R(j, i)) + enddo + enddo + + ! Cheap error estimate. + if (info == kp) then + ! Approximation is exact. + err_est = 0.0_dp + else + em = matmul(E(kp+1:kpp, :p), R(:p, :p)) + err_est = norm2(abs(em)) + endif + + if (err_est <= tol) exit expm_arnoldi + + enddo expm_arnoldi + endif + + if (err_est .le. tol) then + info = kpp + if (p.eq.1) then + write(msg, *) new_line('A'),' Arnoldi approxmation of the action of the exp. propagator converged!' + else + write(msg, *) new_line('A'),' Block Arnoldi approxmation of the action of the exp. propagator converged!' + endif + write(msg, '(A,2(A,I3))') trim(msg)//new_line('A'), 'n° of vectors:', k+1, ' per input vector, total:', kpp + write(msg, '(A,4X,A,E10.4)') trim(msg)//new_line('A'), 'estimated error: ||err_est||_2 = ', err_est + write(msg, '(A,4X,A,E10.4)') trim(msg)//new_line('A'), 'desired tolerance: tol = ', tol + else + info = -1 + if (p.eq.1) then + write(msg, *) 'Arnoldi approxmation of the action of the exp. propagator did not converge!' + else + write(msg, *) 'Block Arnoldi approxmation of the action of the exp. propagator did not converge!' + endif + write(msg, '(A,4X,2(A,I3))') trim(msg)//new_line('A'), 'maximum n° of vectors reached: ', nsteps+1, & + & ' per input vector, total:', kpp + write(msg, '(A,4X,A,E10.4)') trim(msg)//new_line('A'), 'estimated error: ||err_est||_2 = ', err_est + write(msg, '(A,4X,A,E10.4)') trim(msg)//new_line('A'), 'desired tolerance: tol = ', tol + endif + + if (verbose) call logger%log_message(trim(msg), module=this_module, procedure='kexpm_mat_cdp') + + return + end subroutine kexpm_mat_cdp + + subroutine kexpm_vec_variable_dt_cdp(c, A, b, tau, dt, tol, info, trans, verbosity, kdim) class(abstract_vector_cdp), intent(out) :: c !! Best approximation of \( \exp(\tau \mathbf{A}) \mathbf{b} \) in the computed Krylov subspace class(abstract_linop_cdp), intent(in) :: A @@ -1561,14 +2748,14 @@ subroutine kexpm_vec_cdp(c, A, b, tau, dt, tol, info, trans, verbosity, kdim) ! ----- Internal variables ----- integer, parameter :: kmax = 100 - integer :: k, km, kp, nk, istep + integer :: k, km, kp, nk, tstep, rstep ! Arnoldi factorization. class(abstract_vector_cdp), allocatable :: X(:) complex(dp), allocatable :: H(:, :) ! Normaliztion & temporary arrays. complex(dp), allocatable :: E(:, :) class(abstract_vector_cdp), allocatable :: b_ - real(dp) :: err_est, beta, Tend, dt_ + real(dp) :: err_est, beta, Tend, dt_, tol_, scale character(len=128) :: msg ! Optional arguments. @@ -1583,6 +2770,7 @@ subroutine kexpm_vec_cdp(c, A, b, tau, dt, tol, info, trans, verbosity, kdim) nk = nsteps info = 0 + scale = 1.5_dp ! Allocate arrays. allocate(X(nk+1), b_, source=b) @@ -1594,11 +2782,11 @@ subroutine kexpm_vec_cdp(c, A, b, tau, dt, tol, info, trans, verbosity, kdim) ! Input is zero => Output is zero. call c%zero() else - istep = 1 + tstep = 1 Tend = 0.0_dp if (verbose) then - write(msg, '(2(A,E12.6))') 'tau = ', tau, ', dt_0 = ', dt - call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_cdp') + write(msg, '(3(A,E12.6))') 'tau = ', tau, ', dt_0 = ', dt, ', tol = ', tol + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_variable_dt_cdp') end if do while ( Tend < tau ) @@ -1608,19 +2796,20 @@ subroutine kexpm_vec_cdp(c, A, b, tau, dt, tol, info, trans, verbosity, kdim) H = 0.0_dp ! choose time step - dt = min(tau, 2*dt) ! Reuse step size of last call (and try to increase) + dt = min(tau, scale*dt) ! Reuse step size of last call (and try to increase) if (dt > tau-Tend) then dt_ = tau-Tend ! If a small final step is needed to reach Tend, ignore it on restart if (verbose) then write(msg, '(A,E12.6)') 'Filler step, set dt_ = ', dt_ - call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_cdp') + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_variable_dt_cdp') end if else dt_ = dt end if + tol_ = dt_/tau*tol if (verbose) then - write(msg, '(A,I3,3(A,E10.4))') 'step ', istep, ': T = ', Tend, ', dt = ', dt, ', dt_ = ', dt_ - call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_cdp') + write(msg, '(A,I3,4(A,E8.2))') 'step ', tstep, ': T= ', Tend, ', dt= ', dt, ', dt_= ', dt_, ', tol_= ', tol_ + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_variable_dt_cdp') end if expm_arnoldi: do k = 1, nk @@ -1631,7 +2820,7 @@ subroutine kexpm_vec_cdp(c, A, b, tau, dt, tol, info, trans, verbosity, kdim) ! Compute k-th step Arnoldi factorization. call arnoldi(A, X, H, info, kstart=k, kend=k, transpose=transpose) - call check_info(info, 'arnoldi', module=this_module, procedure='kexpm_vec_cdp') + call check_info(info, 'arnoldi', module=this_module, procedure='kexpm_vec_variable_dt_cdp') ! Compute approximation. if (info == k) then @@ -1653,24 +2842,29 @@ subroutine kexpm_vec_cdp(c, A, b, tau, dt, tol, info, trans, verbosity, kdim) endif ! Check convergence. - if (err_est <= tol) exit expm_arnoldi + if (err_est <= tol_) exit expm_arnoldi enddo expm_arnoldi if (verbose) then - write(msg, '(30X,A,E12.6)') 'err_est = ', err_est - call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_cdp') + write(msg, '(34X,A,E12.6)') 'err_est = ', err_est + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_variable_dt_cdp') end if - if (err_est > tol) then ! decrease dt_ until error estimate is below tolerance + if (err_est > tol_) then ! decrease dt_ until error estimate is below tolerance ! NOTE: We do not need to recompute X,H! - do while (err_est > tol) + rstep = 1 + do while (err_est > tol_) + + !hfact = (tol/err_est)**(one_cdp/(kdim+1)) + !dt_ = 0.5_dp*sfact_dp*hfact*dt_ + dt_ = dt_/scale ! this is a quick hack, should choose step better - dt_ = dt_/2 ! this is a quick hack, should choose step better + tol_ = dt_/tau*tol ! update accepted tol based on new time step if (dt_ < atol_dp) then write(msg, *) 'dt < atol_dp. Cannot compute action of the matrix exponential!' - call stop_error(trim(msg), module=this_module, procedure='kexpm_vec_cdp') + call stop_error(trim(msg), module=this_module, procedure='kexpm_vec_variable_dt_cdp') end if ! Compute the (dense) matrix exponential of the extended Hessenberg matrix. @@ -1680,15 +2874,17 @@ subroutine kexpm_vec_cdp(c, A, b, tau, dt, tol, info, trans, verbosity, kdim) err_est = abs(E(kp, 1) * beta) if (verbose) then - write(msg, '(2(A,E12.6))') ' reduce dt_ to ', dt_, ', err_est = ', err_est - call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_cdp') + write(msg, '(6X,3(A,E10.4))') 'v dt_ = ', dt_, ', err_est = ', err_est,', tol_ = ', tol_ + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_variable_dt_cdp') end if + rstep = rstep + 1 + end do ! save dt for reuse dt = dt_ - + end if ! Project back into original space. @@ -1697,23 +2893,28 @@ subroutine kexpm_vec_cdp(c, A, b, tau, dt, tol, info, trans, verbosity, kdim) call linear_combination(xwrk, X(:kp), E(:kp, 1)) call b_%axpby(zero_cdp, xwrk, one_cdp*beta) end block - + ! move forward in time - Tend = Tend + dt + Tend = Tend + dt_ - istep = istep + 1 + tstep = tstep + 1 + + if (verbose) then + write(msg, '(A,E10.4,2(A,I3),A,E10.4)') ' post: T = ', Tend, ', k = ', k, '/', kdim+1, ', dt => ', dt + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_variable_dt_cdp') + end if end do ! integration ! copy result into output vector call c%axpby(zero_cdp, b_, one_cdp) - endif + endif return - end subroutine kexpm_vec_cdp + end subroutine kexpm_vec_variable_dt_cdp - subroutine kexpm_mat_cdp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) + subroutine kexpm_mat_variable_dt_cdp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) class(abstract_vector_cdp), intent(out) :: C(:) !! Best Krylov approximation of \( \mathbf{C} = \exp(\tau \mathbf{A}) \mathbf{B} \). class(abstract_linop_cdp), intent(in) :: A @@ -1737,7 +2938,7 @@ subroutine kexpm_mat_cdp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) ! ----- Internal variables ----- integer, parameter :: kmax = 100 - integer :: i, j, k, p, kpm, kp, kpp, nk, istep + integer :: i, j, k, p, kpm, kp, kpp, nk, tstep ! Block-Arnoldi factorization. class(abstract_vector_cdp), allocatable :: X(:) complex(dp), allocatable :: H(:, :) @@ -1745,7 +2946,7 @@ subroutine kexpm_mat_cdp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) complex(dp), allocatable :: R(:, :), E(:, :), em(:, :) integer, allocatable :: perm(:), ptrans(:) class(abstract_vector_cdp), allocatable :: Xwrk(:), Cwrk(:), B_(:) - real(dp) :: err_est, Tend, dt_ + real(dp) :: err_est, Tend, dt_, tol_, scale character(len=128) :: msg ! Optional arguments. @@ -1763,6 +2964,7 @@ subroutine kexpm_mat_cdp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) nk = nsteps*p info = 0 + scale = 1.5 ! Allocate arrays. allocate(perm(p), ptrans(p)) @@ -1776,7 +2978,7 @@ subroutine kexpm_mat_cdp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) ! Input matrix is zero. call zero_basis(C) else - istep = 1 + tstep = 1 Tend = 0.0_dp do while ( Tend < tau ) @@ -1787,11 +2989,20 @@ subroutine kexpm_mat_cdp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) H = 0.0_dp ! choose time step - dt = min(tau, 2*dt) ! Reuse step size of last call (and try to increase) - dt_ = max(tau-Tend, dt) ! If a small final step is needed to reach Tend, ignore it on restart + dt = min(tau, scale*dt) ! Reuse step size of last call (and try to increase) + if (dt > tau-Tend) then + dt_ = tau-Tend ! If a small final step is needed to reach Tend, ignore it on restart + if (verbose) then + write(msg, '(A,E12.6)') 'Filler step, set dt_ = ', dt_ + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_mat_variable_dt_cdp') + end if + else + dt_ = dt + end if + tol_ = dt_/tau*tol if (verbose) then - write(msg, '(A,I3,4(A,E8.2))') 'step ', istep, ': tau=', tau, ', T=', Tend, ', dt=', dt, ', dt_=', dt_ - call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_cdp') + write(msg, '(A,I3,4(A,E8.2))') 'step ', tstep, ': T= ', Tend, ', dt= ', dt, ', dt_= ', dt_, ', tol_= ', tol_ + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_mat_variable_dt_cdp') end if expm_arnoldi: do k = 1, nk @@ -1803,7 +3014,7 @@ subroutine kexpm_mat_cdp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) ! Compute the k-th step of the Arnoldi factorization. call arnoldi(A, X, H, info, kstart=k, kend=k, transpose=transpose, blksize=p) - call check_info(info, 'arnoldi', module=this_module, procedure='kexpm_mat_cdp') + call check_info(info, 'arnoldi', module=this_module, procedure='kexpm_mat_variable_dt_cdp') if (info == kp) then ! Arnoldi breakdown. Do not consider extended matrix. @@ -1823,24 +3034,28 @@ subroutine kexpm_mat_cdp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) err_est = norm2(abs(em)) endif - if (err_est <= tol) exit expm_arnoldi + if (err_est <= tol_) exit expm_arnoldi enddo expm_arnoldi if (verbose) then - write(msg, '(A,I3,A,I3,A,E8.2)') 'step ', istep, ': k=', k, ', err_est=', err_est - call logger%log_debug(trim(msg), module=this_module, procedure='kexpm_vec_cdp') + write(msg, '(34X,A,E10.4)') 'err_est = ', err_est + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_mat_variable_dt_cdp') end if - if (err_est > tol) then ! decrease dt_ until error estimate is below tolerance + if (err_est > tol_) then ! decrease dt_ until error estimate is below tolerance ! NOTE: We do not need to recompute X,H! - do while (err_est > tol) + do while (err_est > tol_) + + !hfact = (tol/err_est)**(one_cdp/kdim) + !dt_ = sfact_dp*hfact*dt_ + dt_ = dt_/scale ! this is a quick hack, should choose step better - dt_ = dt_/2 ! this is a quick hack, should choose step better + tol_ = dt_/tau*tol ! update accepted tol based on new time step if (dt_ < atol_dp) then write(msg, *) 'dt < atol_dp. Cannot compute action of the matrix exponential!' - call stop_error(trim(msg), module=this_module, procedure='kexpm_vec_cdp') + call stop_error(trim(msg), module=this_module, procedure='kexpm_mat_variable_dt_cdp') end if ! Compute the (dense) matrix exponential of the extended Hessenberg matrix. @@ -1851,8 +3066,8 @@ subroutine kexpm_mat_cdp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) err_est = norm2(abs(em)) if (verbose) then - write(msg, '(2(A,E8.2))') ' : dt_ =', dt_, ', err_est=', err_est - call logger%log_debug(trim(msg), module=this_module, procedure='kexpm_vec_cdp') + write(msg, '(6X,3(A,E12.6))') 'v dt_ = ', dt_, ', err_est = ', err_est, ', tol_ = ', tol_ + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_mat_variable_dt_cdp') end if end do @@ -1869,9 +3084,14 @@ subroutine kexpm_mat_cdp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) end block ! move forward in time - Tend = Tend + dt + Tend = Tend + dt_ - istep = istep + 1 + tstep = tstep + 1 + + if (verbose) then + write(msg, '(A,E10.4,2(A,I3),A,E10.4)') ' post: T = ', Tend, ', k = ', k, '/', kdim+1, ', dt => ', dt + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_mat_variable_dt_cdp') + end if end do ! integration @@ -1880,7 +3100,7 @@ subroutine kexpm_mat_cdp(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) endif return - end subroutine kexpm_mat_cdp + end subroutine kexpm_mat_variable_dt_cdp subroutine k_exptA_cdp(vec_out, A, vec_in, tau, dt, info, trans) class(abstract_vector_cdp), intent(out) :: vec_out @@ -1907,7 +3127,7 @@ subroutine k_exptA_cdp(vec_out, A, vec_in, tau, dt, info, trans) kdim = 30 verbose = .false. - call kexpm(vec_out, A, vec_in, tau, tol, dt, info, trans=trans, verbosity=verbose, kdim=kdim) + call kexpm_var_dt(vec_out, A, vec_in, tau, tol, dt, info, trans=trans, verbosity=verbose, kdim=kdim) call check_info(info, 'kexpm', module=this_module, procedure='k_exptA_cdp') return diff --git a/src/ExpmLib.fypp b/src/ExpmLib.fypp index b76d9a2..a6613b6 100644 --- a/src/ExpmLib.fypp +++ b/src/ExpmLib.fypp @@ -10,12 +10,12 @@ 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 @@ -26,6 +26,7 @@ module lightkrylov_expmlib #:endfor public :: expm public :: kexpm + public :: kexpm_var_dt public :: k_exptA abstract interface @@ -67,6 +68,13 @@ module lightkrylov_expmlib #:endfor end interface + interface kexpm_var_dt + #:for kind, type in RC_KINDS_TYPES + module procedure kexpm_vec_variable_dt_${type[0]}$${kind}$ + module procedure kexpm_mat_variable_dt_${type[0]}$${kind}$ + #:endfor + end interface + interface k_exptA #:for kind, type in RC_KINDS_TYPES module procedure k_exptA_${type[0]}$${kind}$ @@ -146,7 +154,276 @@ contains return end subroutine expm_${type[0]}$${kind}$ - subroutine kexpm_vec_${type[0]}$${kind}$(c, A, b, tau, dt, tol, info, trans, verbosity, kdim) + subroutine kexpm_vec_${type[0]}$${kind}$(c, A, b, tau, tol, info, trans, verbosity, kdim) + class(abstract_vector_${type[0]}$${kind}$), intent(out) :: c + !! Best approximation of \( \exp(\tau \mathbf{A}) \mathbf{b} \) in the computed Krylov subspace + class(abstract_linop_${type[0]}$${kind}$), intent(in) :: A + !! Linear operator to be exponentiated. + class(abstract_vector_${type[0]}$${kind}$), intent(in) :: b + !! Input vector on which to apply \( \exp(\tau \mathbf{A}) \). + real(${kind}$), intent(in) :: tau + !! Time horizon for the exponentiation. + real(${kind}$), intent(in) :: tol + !! Solution tolerance based on error estimates. + integer, intent(out) :: info + !! Information flag. + logical, optional, intent(in) :: trans + !! Use transpose? + logical, optional, intent(in) :: verbosity + !! Verbosity control. + integer, optional, intent(in) :: kdim + !! Maximum size of the Krylov subspace. + + ! ----- Internal variables ----- + integer, parameter :: kmax = 100 + integer :: k, km, kp, nk + ! Arnoldi factorization. + class(abstract_vector_${type[0]}$${kind}$), allocatable :: X(:) + ${type}$, allocatable :: H(:, :) + ! Normaliztion & temporary arrays. + ${type}$, allocatable :: E(:, :) + class(abstract_vector_${type[0]}$${kind}$), allocatable :: Xwrk + real(${kind}$) :: err_est, beta + character(len=512) :: msg + + ! Optional arguments. + logical :: transpose + logical :: verbose + integer :: nsteps + + ! Deals with optional args. + transpose = optval(trans, .false.) + verbose = optval(verbosity, .false.) + nsteps = optval(kdim, kmax) + nk = nsteps + + info = 0 + + ! Allocate arrays. + allocate(X(nk+1), source=b) ; allocate(H(nk+1, nk+1)) + allocate(E(nk+1, nk+1)) ; allocate(Xwrk, source=b) + + ! Normalize input vector and initialize Krylov subspace. + beta = b%norm() + if (beta == 0.0_${kind}$) then + ! Input is zero => Output is zero. + call c%zero() + err_est = 0.0_${kind}$ + kp = 1 + else + call initialize_krylov_subspace(X) + call X(1)%add(b) ; call X(1)%scal(one_${type[0]}$${kind}$ / beta) + H = 0.0_${kind}$ + + expm_arnoldi: do k = 1, nk + km = k-1 ; kp = k+1 + + ! Reset work arrays. + E = 0.0_${kind}$ + + ! Compute k-th step Arnoldi factorization. + call arnoldi(A, X, H, info, kstart=k, kend=k, transpose=transpose) + call check_info(info, 'arnoldi', module=this_module, procedure='kexpm_vec_${type[0]}$${kind}$') + + ! Compute approximation. + if (info == k) then + ! Arnoldi breakdown, do not consider extended matrix. + kp = k + info = -2 + endif + + ! Compute the (dense) matrix exponential of the extended Hessenberg matrix. + call expm(E(:kp, :kp), tau*H(:kp, :kp)) + + ! Project back into original space. + call linear_combination(Xwrk, X(:kp), E(:kp, 1)) + call c%axpby(zero_${type[0]}$${kind}$, Xwrk, one_${type[0]}$${kind}$*beta) + + ! Cheap error esimate (this actually is the magnitude of the included correction + ! and thus is very conservative). + if (info == k) then + ! Approximation is exact. + err_est = 0.0_${kind}$ + else + err_est = abs(E(kp, 1) * beta) + endif + + ! Check convergence. + if (err_est <= tol) exit expm_arnoldi + + enddo expm_arnoldi + endif + + if (err_est <= tol) then + info = kp + write(msg, *) new_line('A'),' Arnoldi approxmation of the action of the exp. propagator converged!' + write(msg, '(A,4X,A,I3)') trim(msg)//new_line('A'), 'n° of vectors:', kp + write(msg, '(A,4X,A,E10.4)') trim(msg)//new_line('A'), 'estimated error: ||err_est||_2 = ', err_est + write(msg, '(A,4X,A,E10.4)') trim(msg)//new_line('A'), 'desired tolerance: tol = ', tol + else + info = -1 + write(msg, *) new_line('A'),' Arnoldi approxmation of the action of the exp. propagator did not converge!' + write(msg, '(A,4X,A,I3)') trim(msg)//new_line('A'), 'Maximum n° of vectors reached: ', nk + 1 + write(msg, '(A,4X,A,E10.4)') trim(msg)//new_line('A'), 'Estimated error : ||err_est||_2 = ', err_est + write(msg, '(A,4X,A,E10.4)') trim(msg)//new_line('A'), 'Desired tolerance : tol = ', tol + endif + if (verbose) call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_${type[0]}$${kind}$') + + return + end subroutine kexpm_vec_${type[0]}$${kind}$ + + subroutine kexpm_mat_${type[0]}$${kind}$(C, A, B, tau, tol, info, trans, verbosity, kdim) + class(abstract_vector_${type[0]}$${kind}$), intent(out) :: C(:) + !! Best Krylov approximation of \( \mathbf{C} = \exp(\tau \mathbf{A}) \mathbf{B} \). + class(abstract_linop_${type[0]}$${kind}$), intent(in) :: A + !! Linear operator to be exponentiated. + class(abstract_vector_${type[0]}$${kind}$), intent(in) :: B(:) + !! Input matrix on which to apply \( \exp(\tau \mathbf{A}) \). + real(${kind}$), intent(in) :: tau + !! Time horizon for the exponentiation. + real(${kind}$), intent(in) :: tol + !! Solution toleance based on error estimates. + integer, intent(out) :: info + !! Information flag. + logical, optional, intent(in) :: trans + !! Use transpose ? + logical, optional, intent(in) :: verbosity + !! Verbosity control. + integer, optional, intent(in) :: kdim + !! Maximum size of the Krylov subspace. + + ! ----- Internal variables ----- + integer, parameter :: kmax = 100 + integer :: i, j, k, p, kpm, kp, kpp, nk + ! Block-Arnoldi factorization. + class(abstract_vector_${type[0]}$${kind}$), allocatable :: X(:) + ${type}$, allocatable :: H(:, :) + ! Normalization & temporary arrays. + ${type}$, allocatable :: R(:, :), E(:, :), em(:, :) + integer, allocatable :: perm(:), ptrans(:) + class(abstract_vector_${type[0]}$${kind}$), allocatable :: Xwrk(:), Cwrk(:) + real(${kind}$) :: err_est + character(len=512) :: msg + + ! Optional arguments. + logical :: transpose + logical :: verbose + integer :: nsteps + + ! Determine block size. + p = size(B) + + ! Deals with the optional args. + transpose = optval(trans, .false.) + verbose = optval(verbosity, .false.) + nsteps = optval(kdim, kmax) + nk = nsteps*p + + info = 0 + + ! Allocate arrays. + allocate(R(p, p)) ; allocate(perm(p)) ; allocate(ptrans(p)) + allocate(X(p*(nk+1)), source=B(1)) ; allocate(H(p*(nk+1), p*(nk+1))) + allocate(E(p*(nk+1), p*(nk+1))) ; allocate(em(p, p)) + + ! Scratch arrays. + allocate(Xwrk(p), source=B) ; allocate(Cwrk(p), source=B(1)) + + ! Normalize input matrix and initialize Krylov subspace. + R = 0.0_${kind}$ + call qr(Xwrk, R, perm, info) ; call apply_inverse_permutation_matrix(R, perm) + + if (norm2(abs(R)) == 0.0_${kind}$) then + ! Input matrix is zero. + do i = 1, size(C) + call C(i)%zero() + enddo + err_est = 0.0_${kind}$ ; k = 0 ; kpp = p + else + call initialize_krylov_subspace(X, Xwrk) ; H = 0.0_${kind}$ + + expm_arnoldi: do k = 1, nk + ! Set counters. + kpm = (k-1)*p ; kp = kpm + p ; kpp = kp + p + + ! Reset working arrays. + E = 0.0_${kind}$ + do i = 1, size(Xwrk) + call Xwrk(i)%zero() + enddo + + ! Compute the k-th step of the Arnoldi factorization. + call arnoldi(A, X, H, info, kstart=k, kend=k, transpose=transpose, blksize=p) + call check_info(info, 'arnoldi', module=this_module, procedure='kexpm_mat_${type[0]}$${kind}$') + + + if (info == kp) then + ! Arnoldi breakdown. Do not consider extended matrix. + kpp = kp + info = -2 + endif + + ! Compute the (dense) matrix exponential of the extended Hessenberg matrix. + call expm(E(:kpp, :kpp), tau*H(:kpp, :kpp)) + + ! Project back to original space. + do i = 1, size(Xwrk) + call Xwrk(i)%zero() + do j = 1, kpp + call Xwrk(i)%axpby(one_${type[0]}$${kind}$, X(j), E(j, i)) + enddo + enddo + + do i = 1, p + call C(i)%zero() + do j = 1, p + call C(i)%axpby(one_${type[0]}$${kind}$, Xwrk(j), R(j, i)) + enddo + enddo + + ! Cheap error estimate. + if (info == kp) then + ! Approximation is exact. + err_est = 0.0_${kind}$ + else + em = matmul(E(kp+1:kpp, :p), R(:p, :p)) + err_est = norm2(abs(em)) + endif + + if (err_est <= tol) exit expm_arnoldi + + enddo expm_arnoldi + endif + + if (err_est .le. tol) then + info = kpp + if (p.eq.1) then + write(msg, *) new_line('A'),' Arnoldi approxmation of the action of the exp. propagator converged!' + else + write(msg, *) new_line('A'),' Block Arnoldi approxmation of the action of the exp. propagator converged!' + endif + write(msg, '(A,2(A,I3))') trim(msg)//new_line('A'), 'n° of vectors:', k+1, ' per input vector, total:', kpp + write(msg, '(A,4X,A,E10.4)') trim(msg)//new_line('A'), 'estimated error: ||err_est||_2 = ', err_est + write(msg, '(A,4X,A,E10.4)') trim(msg)//new_line('A'), 'desired tolerance: tol = ', tol + else + info = -1 + if (p.eq.1) then + write(msg, *) 'Arnoldi approxmation of the action of the exp. propagator did not converge!' + else + write(msg, *) 'Block Arnoldi approxmation of the action of the exp. propagator did not converge!' + endif + write(msg, '(A,4X,2(A,I3))') trim(msg)//new_line('A'), 'maximum n° of vectors reached: ', nsteps+1, & + & ' per input vector, total:', kpp + write(msg, '(A,4X,A,E10.4)') trim(msg)//new_line('A'), 'estimated error: ||err_est||_2 = ', err_est + write(msg, '(A,4X,A,E10.4)') trim(msg)//new_line('A'), 'desired tolerance: tol = ', tol + endif + + if (verbose) call logger%log_message(trim(msg), module=this_module, procedure='kexpm_mat_${type[0]}$${kind}$') + + return + end subroutine kexpm_mat_${type[0]}$${kind}$ + + subroutine kexpm_vec_variable_dt_${type[0]}$${kind}$(c, A, b, tau, dt, tol, info, trans, verbosity, kdim) class(abstract_vector_${type[0]}$${kind}$), intent(out) :: c !! Best approximation of \( \exp(\tau \mathbf{A}) \mathbf{b} \) in the computed Krylov subspace class(abstract_linop_${type[0]}$${kind}$), intent(in) :: A @@ -170,14 +447,14 @@ contains ! ----- Internal variables ----- integer, parameter :: kmax = 100 - integer :: k, km, kp, nk, istep + integer :: k, km, kp, nk, tstep, rstep ! Arnoldi factorization. class(abstract_vector_${type[0]}$${kind}$), allocatable :: X(:) ${type}$, allocatable :: H(:, :) ! Normaliztion & temporary arrays. ${type}$, allocatable :: E(:, :) class(abstract_vector_${type[0]}$${kind}$), allocatable :: b_ - real(${kind}$) :: err_est, beta, Tend, dt_ + real(${kind}$) :: err_est, beta, Tend, dt_, tol_, scale character(len=128) :: msg ! Optional arguments. @@ -192,6 +469,7 @@ contains nk = nsteps info = 0 + scale = 1.5_${kind}$ ! Allocate arrays. allocate(X(nk+1), b_, source=b) @@ -203,11 +481,11 @@ contains ! Input is zero => Output is zero. call c%zero() else - istep = 1 + tstep = 1 Tend = 0.0_${kind}$ if (verbose) then - write(msg, '(2(A,E12.6))') 'tau = ', tau, ', dt_0 = ', dt - call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_${type[0]}$${kind}$') + write(msg, '(3(A,E12.6))') 'tau = ', tau, ', dt_0 = ', dt, ', tol = ', tol + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_variable_dt_${type[0]}$${kind}$') end if do while ( Tend < tau ) @@ -217,19 +495,20 @@ contains H = 0.0_${kind}$ ! choose time step - dt = min(tau, 2*dt) ! Reuse step size of last call (and try to increase) + dt = min(tau, scale*dt) ! Reuse step size of last call (and try to increase) if (dt > tau-Tend) then dt_ = tau-Tend ! If a small final step is needed to reach Tend, ignore it on restart if (verbose) then write(msg, '(A,E12.6)') 'Filler step, set dt_ = ', dt_ - call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_${type[0]}$${kind}$') + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_variable_dt_${type[0]}$${kind}$') end if else dt_ = dt end if + tol_ = dt_/tau*tol if (verbose) then - write(msg, '(A,I3,3(A,E10.4))') 'step ', istep, ': T = ', Tend, ', dt = ', dt, ', dt_ = ', dt_ - call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_${type[0]}$${kind}$') + write(msg, '(A,I3,4(A,E8.2))') 'step ', tstep, ': T= ', Tend, ', dt= ', dt, ', dt_= ', dt_, ', tol_= ', tol_ + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_variable_dt_${type[0]}$${kind}$') end if expm_arnoldi: do k = 1, nk @@ -240,7 +519,7 @@ contains ! Compute k-th step Arnoldi factorization. call arnoldi(A, X, H, info, kstart=k, kend=k, transpose=transpose) - call check_info(info, 'arnoldi', module=this_module, procedure='kexpm_vec_${type[0]}$${kind}$') + call check_info(info, 'arnoldi', module=this_module, procedure='kexpm_vec_variable_dt_${type[0]}$${kind}$') ! Compute approximation. if (info == k) then @@ -262,24 +541,29 @@ contains endif ! Check convergence. - if (err_est <= tol) exit expm_arnoldi + if (err_est <= tol_) exit expm_arnoldi enddo expm_arnoldi if (verbose) then - write(msg, '(30X,A,E12.6)') 'err_est = ', err_est - call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_${type[0]}$${kind}$') + write(msg, '(34X,A,E12.6)') 'err_est = ', err_est + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_variable_dt_${type[0]}$${kind}$') end if - if (err_est > tol) then ! decrease dt_ until error estimate is below tolerance + if (err_est > tol_) then ! decrease dt_ until error estimate is below tolerance ! NOTE: We do not need to recompute X,H! - do while (err_est > tol) + rstep = 1 + do while (err_est > tol_) + + !hfact = (tol/err_est)**(one_${type[0]}$${kind}$/(kdim+1)) + !dt_ = 0.5_${kind}$*sfact_${kind}$*hfact*dt_ + dt_ = dt_/scale ! this is a quick hack, should choose step better - dt_ = dt_/2 ! this is a quick hack, should choose step better + tol_ = dt_/tau*tol ! update accepted tol based on new time step if (dt_ < atol_${kind}$) then write(msg, *) 'dt < atol_${kind}$. Cannot compute action of the matrix exponential!' - call stop_error(trim(msg), module=this_module, procedure='kexpm_vec_${type[0]}$${kind}$') + call stop_error(trim(msg), module=this_module, procedure='kexpm_vec_variable_dt_${type[0]}$${kind}$') end if ! Compute the (dense) matrix exponential of the extended Hessenberg matrix. @@ -289,15 +573,17 @@ contains err_est = abs(E(kp, 1) * beta) if (verbose) then - write(msg, '(2(A,E12.6))') ' reduce dt_ to ', dt_, ', err_est = ', err_est - call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_${type[0]}$${kind}$') + write(msg, '(6X,3(A,E10.4))') 'v dt_ = ', dt_, ', err_est = ', err_est,', tol_ = ', tol_ + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_variable_dt_${type[0]}$${kind}$') end if + rstep = rstep + 1 + end do ! save dt for reuse dt = dt_ - + end if ! Project back into original space. @@ -306,23 +592,28 @@ contains call linear_combination(xwrk, X(:kp), E(:kp, 1)) call b_%axpby(zero_${type[0]}$${kind}$, xwrk, one_${type[0]}$${kind}$*beta) end block - + ! move forward in time - Tend = Tend + dt + Tend = Tend + dt_ - istep = istep + 1 + tstep = tstep + 1 + + if (verbose) then + write(msg, '(A,E10.4,2(A,I3),A,E10.4)') ' post: T = ', Tend, ', k = ', k, '/', kdim+1, ', dt => ', dt + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_variable_dt_${type[0]}$${kind}$') + end if end do ! integration ! copy result into output vector call c%axpby(zero_${type[0]}$${kind}$, b_, one_${type[0]}$${kind}$) - endif + endif return - end subroutine kexpm_vec_${type[0]}$${kind}$ + end subroutine kexpm_vec_variable_dt_${type[0]}$${kind}$ - subroutine kexpm_mat_${type[0]}$${kind}$(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) + subroutine kexpm_mat_variable_dt_${type[0]}$${kind}$(C, A, B, tau, dt, tol, info, trans, verbosity, kdim) class(abstract_vector_${type[0]}$${kind}$), intent(out) :: C(:) !! Best Krylov approximation of \( \mathbf{C} = \exp(\tau \mathbf{A}) \mathbf{B} \). class(abstract_linop_${type[0]}$${kind}$), intent(in) :: A @@ -346,7 +637,7 @@ contains ! ----- Internal variables ----- integer, parameter :: kmax = 100 - integer :: i, j, k, p, kpm, kp, kpp, nk, istep + integer :: i, j, k, p, kpm, kp, kpp, nk, tstep ! Block-Arnoldi factorization. class(abstract_vector_${type[0]}$${kind}$), allocatable :: X(:) ${type}$, allocatable :: H(:, :) @@ -354,7 +645,7 @@ contains ${type}$, allocatable :: R(:, :), E(:, :), em(:, :) integer, allocatable :: perm(:), ptrans(:) class(abstract_vector_${type[0]}$${kind}$), allocatable :: Xwrk(:), Cwrk(:), B_(:) - real(${kind}$) :: err_est, Tend, dt_ + real(${kind}$) :: err_est, Tend, dt_, tol_, scale character(len=128) :: msg ! Optional arguments. @@ -372,6 +663,7 @@ contains nk = nsteps*p info = 0 + scale = 1.5 ! Allocate arrays. allocate(perm(p), ptrans(p)) @@ -385,7 +677,7 @@ contains ! Input matrix is zero. call zero_basis(C) else - istep = 1 + tstep = 1 Tend = 0.0_${kind}$ do while ( Tend < tau ) @@ -396,11 +688,20 @@ contains H = 0.0_${kind}$ ! choose time step - dt = min(tau, 2*dt) ! Reuse step size of last call (and try to increase) - dt_ = max(tau-Tend, dt) ! If a small final step is needed to reach Tend, ignore it on restart + dt = min(tau, scale*dt) ! Reuse step size of last call (and try to increase) + if (dt > tau-Tend) then + dt_ = tau-Tend ! If a small final step is needed to reach Tend, ignore it on restart + if (verbose) then + write(msg, '(A,E12.6)') 'Filler step, set dt_ = ', dt_ + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_mat_variable_dt_${type[0]}$${kind}$') + end if + else + dt_ = dt + end if + tol_ = dt_/tau*tol if (verbose) then - write(msg, '(A,I3,4(A,E8.2))') 'step ', istep, ': tau=', tau, ', T=', Tend, ', dt=', dt, ', dt_=', dt_ - call logger%log_message(trim(msg), module=this_module, procedure='kexpm_vec_${type[0]}$${kind}$') + write(msg, '(A,I3,4(A,E8.2))') 'step ', tstep, ': T= ', Tend, ', dt= ', dt, ', dt_= ', dt_, ', tol_= ', tol_ + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_mat_variable_dt_${type[0]}$${kind}$') end if expm_arnoldi: do k = 1, nk @@ -412,7 +713,7 @@ contains ! Compute the k-th step of the Arnoldi factorization. call arnoldi(A, X, H, info, kstart=k, kend=k, transpose=transpose, blksize=p) - call check_info(info, 'arnoldi', module=this_module, procedure='kexpm_mat_${type[0]}$${kind}$') + call check_info(info, 'arnoldi', module=this_module, procedure='kexpm_mat_variable_dt_${type[0]}$${kind}$') if (info == kp) then ! Arnoldi breakdown. Do not consider extended matrix. @@ -432,24 +733,28 @@ contains err_est = norm2(abs(em)) endif - if (err_est <= tol) exit expm_arnoldi + if (err_est <= tol_) exit expm_arnoldi enddo expm_arnoldi if (verbose) then - write(msg, '(A,I3,A,I3,A,E8.2)') 'step ', istep, ': k=', k, ', err_est=', err_est - call logger%log_debug(trim(msg), module=this_module, procedure='kexpm_vec_${type[0]}$${kind}$') + write(msg, '(34X,A,E10.4)') 'err_est = ', err_est + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_mat_variable_dt_${type[0]}$${kind}$') end if - if (err_est > tol) then ! decrease dt_ until error estimate is below tolerance + if (err_est > tol_) then ! decrease dt_ until error estimate is below tolerance ! NOTE: We do not need to recompute X,H! - do while (err_est > tol) + do while (err_est > tol_) + + !hfact = (tol/err_est)**(one_${type[0]}$${kind}$/kdim) + !dt_ = sfact_${kind}$*hfact*dt_ + dt_ = dt_/scale ! this is a quick hack, should choose step better - dt_ = dt_/2 ! this is a quick hack, should choose step better + tol_ = dt_/tau*tol ! update accepted tol based on new time step if (dt_ < atol_${kind}$) then write(msg, *) 'dt < atol_${kind}$. Cannot compute action of the matrix exponential!' - call stop_error(trim(msg), module=this_module, procedure='kexpm_vec_${type[0]}$${kind}$') + call stop_error(trim(msg), module=this_module, procedure='kexpm_mat_variable_dt_${type[0]}$${kind}$') end if ! Compute the (dense) matrix exponential of the extended Hessenberg matrix. @@ -460,8 +765,8 @@ contains err_est = norm2(abs(em)) if (verbose) then - write(msg, '(2(A,E8.2))') ' : dt_ =', dt_, ', err_est=', err_est - call logger%log_debug(trim(msg), module=this_module, procedure='kexpm_vec_${type[0]}$${kind}$') + write(msg, '(6X,3(A,E12.6))') 'v dt_ = ', dt_, ', err_est = ', err_est, ', tol_ = ', tol_ + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_mat_variable_dt_${type[0]}$${kind}$') end if end do @@ -478,9 +783,14 @@ contains end block ! move forward in time - Tend = Tend + dt + Tend = Tend + dt_ - istep = istep + 1 + tstep = tstep + 1 + + if (verbose) then + write(msg, '(A,E10.4,2(A,I3),A,E10.4)') ' post: T = ', Tend, ', k = ', k, '/', kdim+1, ', dt => ', dt + call logger%log_message(trim(msg), module=this_module, procedure='kexpm_mat_variable_dt_${type[0]}$${kind}$') + end if end do ! integration @@ -489,7 +799,7 @@ contains endif return - end subroutine kexpm_mat_${type[0]}$${kind}$ + end subroutine kexpm_mat_variable_dt_${type[0]}$${kind}$ subroutine k_exptA_${type[0]}$${kind}$(vec_out, A, vec_in, tau, dt, info, trans) class(abstract_vector_${type[0]}$${kind}$), intent(out) :: vec_out @@ -516,7 +826,7 @@ contains kdim = 30 verbose = .false. - call kexpm(vec_out, A, vec_in, tau, tol, dt, info, trans=trans, verbosity=verbose, kdim=kdim) + call kexpm_var_dt(vec_out, A, vec_in, tau, tol, dt, info, trans=trans, verbosity=verbose, kdim=kdim) call check_info(info, 'kexpm', module=this_module, procedure='k_exptA_${type[0]}$${kind}$') return diff --git a/test/TestExpmlib.f90 b/test/TestExpmlib.f90 index 015c2b2..e28076a 100644 --- a/test/TestExpmlib.f90 +++ b/test/TestExpmlib.f90 @@ -78,7 +78,7 @@ subroutine test_dense_expm_rsp(error) ! Check result. err = maxval(abs(E-Eref)) - call get_err_str(msg, "max err: ", err) + call get_err_str(msg, "max err: ", err, rtol_sp) call check(error, err < rtol_sp) call check_test(error, 'test_dense_expm_rsp', eq='maxval(abs(E-Eref))', context=msg) @@ -99,7 +99,7 @@ subroutine test_kexptA_rsp(error) real(sp), allocatable :: E(:, :) real(sp), parameter :: tau = 0.1_sp integer, parameter :: nkmax = 50 - real(sp) :: err, dt + real(sp) :: err integer :: info character(len=256) :: msg @@ -116,13 +116,12 @@ subroutine test_kexptA_rsp(error) Xref%data = matmul(E, Q%data) ! Krylov exponential. - dt = tau - call kexpm(Xkryl, A, Q, tau, dt, rtol_sp, info, verbosity=verb, kdim=nkmax) + call kexpm(Xkryl, A, Q, tau, rtol_sp, info, verbosity=verb, kdim=nkmax) call check_info(info, 'kexpm', module=this_module, procedure='test_kexptA_rsp') ! Check result. call Xkryl%sub(Xref) ; err = Xkryl%norm() - call get_err_str(msg, "max err: ", err) + call get_err_str(msg, "max err: ", err, rtol_sp) call check(error, err < rtol_sp) call check_test(error, 'test_kexptA_rsp', & & eq='Comparison with dense matrix exponential', context=msg) @@ -142,13 +141,13 @@ subroutine test_kexptA_variable_dt_rsp(error) ! ----- Internal variables ----- real(sp), allocatable :: E(:, :) - real(sp), parameter :: tau = 0.1_sp + real(sp), parameter :: tau = 0.5_sp integer, parameter :: nkmax = 10 real(sp) :: err, dt integer :: info character(len=256) :: msg - logical, parameter :: verb = .true. + logical, parameter :: verb = .false. ! Initialize data. A = linop_rsp() ; call init_rand(A) @@ -162,15 +161,15 @@ subroutine test_kexptA_variable_dt_rsp(error) ! Krylov exponential. dt = tau - call kexpm(Xkryl, A, Q, tau, dt, rtol_sp, info, verbosity=verb, kdim=nkmax) - call check_info(info, 'kexpm', module=this_module, procedure='test_kexptA_rsp') + call kexpm_var_dt(Xkryl, A, Q, tau, dt, rtol_sp, info, verbosity=verb, kdim=nkmax) + call check_info(info, 'kexpm_var_dt', module=this_module, procedure='test_kexptA_variable_dt_rsp') ! Check result. call Xkryl%sub(Xref) ; err = Xkryl%norm() - call get_err_str(msg, "max err: ", err) + call get_err_str(msg, "max err: ", err, rtol_sp) call check(error, err < rtol_sp) - call check_test(error, 'test_kexptA_rsp', & - & info='Comparison with matrix exponential', context=msg) + call check_test(error, 'test_kexptA_variable_dt_rsp', & + & eq='Comparison with dense matrix exponential', context=msg) return end subroutine test_kexptA_variable_dt_rsp @@ -201,7 +200,7 @@ subroutine test_block_kexptA_rsp(error) ! Misc. integer :: i, j real(sp), allocatable :: Xdata(:, :), Qdata(:, :) - real(sp) :: errv, dt + real(sp) :: errv real(sp) :: err(p, p) character(len=256) :: msg @@ -227,14 +226,12 @@ subroutine test_block_kexptA_rsp(error) ! Compute Krylov matrix exponential using sequential arnoldi method for each input column do i = 1,p - dt = tau - call kexpm(C(i), A, B(i), tau, dt, tol, info, verbosity=verb, kdim=nkmax) + call kexpm(C(i), A, B(i), tau, tol, info, verbosity=verb, kdim=nkmax) call check_info(info, 'kexpm', module=this_module, procedure='test_block_kexptA_rsp, 1') end do ! Compute Krylov matrix exponential using block-arnoldi method - dt = tau - call kexpm(Cblk, A, B, tau, dt, tol, info, verbosity=verb, kdim=nkmax) + call kexpm(Cblk, A, B, tau, tol, info, verbosity=verb, kdim=nkmax) call check_info(info, 'kexpm', module=this_module, procedure='test_block_kexptA_rsp, 2') do i = 1, p @@ -250,7 +247,7 @@ subroutine test_block_kexptA_rsp(error) enddo enddo errv = norm2(abs(err)) - !call get_err_str(msg, "max err: ", errv) + !call get_err_str(msg, "max err: ", errv, rtol_sp) do i = 1, size(Cblk) do j = 1, size(Cblk) @@ -258,7 +255,7 @@ subroutine test_block_kexptA_rsp(error) enddo enddo errv = norm2(abs(err)) - call get_err_str(msg, "max err: ", errv) + call get_err_str(msg, "max err: ", errv, rtol_sp) call check(error, errv < rtol_sp) call check_test(error, 'test_block_kexptA_rsp', & @@ -288,7 +285,7 @@ subroutine test_block_kexptA_variable_dt_rsp(error) ! Test parameters integer , parameter :: nkmax = 10 integer , parameter :: p = 3 - real(sp), parameter :: tau = 0.1_sp + real(sp), parameter :: tau = 0.5_sp real(sp), parameter :: tol = rtol_sp ! Misc. integer :: i, j @@ -297,7 +294,7 @@ subroutine test_block_kexptA_variable_dt_rsp(error) real(sp) :: err(p, p) character(len=256) :: msg - logical, parameter :: verb = .true. + logical, parameter :: verb = .false. allocate(Adata(kdim, kdim)) ; Adata = 0.0_sp allocate(Edata(kdim, kdim)) ; Edata = 0.0_sp @@ -318,16 +315,16 @@ subroutine test_block_kexptA_variable_dt_rsp(error) call expm(Edata, tau*Adata) ; Xdata = matmul(Edata,Qdata) ; call put_data(Cref, Xdata) ! Compute Krylov matrix exponential using sequential arnoldi method for each input column + dt = tau do i = 1,p - dt = tau - call kexpm(C(i), A, B(i), tau, dt, tol, info, verbosity=verb, kdim=nkmax) - call check_info(info, 'kexpm', module=this_module, procedure='test_block_kexptA_rsp, 1') + call kexpm_var_dt(C(i), A, B(i), tau, dt, tol, info, verbosity=verb, kdim=nkmax) + call check_info(info, 'kexpm_var_dt', module=this_module, procedure='test_block_kexptA_variable_dt_rsp, 1') end do ! Compute Krylov matrix exponential using block-arnoldi method dt = tau - call kexpm(Cblk, A, B, tau, dt, tol, info, verbosity=verb, kdim=nkmax) - call check_info(info, 'kexpm', module=this_module, procedure='test_block_kexptA_rsp, 2') + call kexpm_var_dt(Cblk, A, B, tau, dt, tol, info, verbosity=verb, kdim=nkmax) + call check_info(info, 'kexpm_var_dt', module=this_module, procedure='test_block_kexptA_variable_dt_rsp, 2') do i = 1, p write(output_unit, *) C(i)%norm(), Cblk(i)%norm() @@ -342,7 +339,7 @@ subroutine test_block_kexptA_variable_dt_rsp(error) enddo enddo errv = norm2(abs(err)) - !call get_err_str(msg, "max err: ", errv) + !call get_err_str(msg, "max err: ", errv, rtol_sp) do i = 1, size(Cblk) do j = 1, size(Cblk) @@ -350,11 +347,11 @@ subroutine test_block_kexptA_variable_dt_rsp(error) enddo enddo errv = norm2(abs(err)) - call get_err_str(msg, "max err: ", errv) + call get_err_str(msg, "max err: ", errv, rtol_sp) call check(error, errv < rtol_sp) - call check_test(error, 'test_block_kexptA_rsp', & - & info='Comparison with matrix exponential', context=msg) + call check_test(error, 'test_block_kexptA_variable_dt_rsp', & + & eq='Comparison with matrix exponential', context=msg) return end subroutine test_block_kexptA_variable_dt_rsp @@ -401,7 +398,7 @@ subroutine test_dense_expm_rdp(error) ! Check result. err = maxval(abs(E-Eref)) - call get_err_str(msg, "max err: ", err) + call get_err_str(msg, "max err: ", err, rtol_dp) call check(error, err < rtol_dp) call check_test(error, 'test_dense_expm_rdp', eq='maxval(abs(E-Eref))', context=msg) @@ -422,7 +419,7 @@ subroutine test_kexptA_rdp(error) real(dp), allocatable :: E(:, :) real(dp), parameter :: tau = 0.1_dp integer, parameter :: nkmax = 50 - real(dp) :: err, dt + real(dp) :: err integer :: info character(len=256) :: msg @@ -439,13 +436,12 @@ subroutine test_kexptA_rdp(error) Xref%data = matmul(E, Q%data) ! Krylov exponential. - dt = tau - call kexpm(Xkryl, A, Q, tau, dt, rtol_dp, info, verbosity=verb, kdim=nkmax) + call kexpm(Xkryl, A, Q, tau, rtol_dp, info, verbosity=verb, kdim=nkmax) call check_info(info, 'kexpm', module=this_module, procedure='test_kexptA_rdp') ! Check result. call Xkryl%sub(Xref) ; err = Xkryl%norm() - call get_err_str(msg, "max err: ", err) + call get_err_str(msg, "max err: ", err, rtol_dp) call check(error, err < rtol_dp) call check_test(error, 'test_kexptA_rdp', & & eq='Comparison with dense matrix exponential', context=msg) @@ -465,13 +461,13 @@ subroutine test_kexptA_variable_dt_rdp(error) ! ----- Internal variables ----- real(dp), allocatable :: E(:, :) - real(dp), parameter :: tau = 0.1_dp + real(dp), parameter :: tau = 0.5_dp integer, parameter :: nkmax = 10 real(dp) :: err, dt integer :: info character(len=256) :: msg - logical, parameter :: verb = .true. + logical, parameter :: verb = .false. ! Initialize data. A = linop_rdp() ; call init_rand(A) @@ -485,15 +481,15 @@ subroutine test_kexptA_variable_dt_rdp(error) ! Krylov exponential. dt = tau - call kexpm(Xkryl, A, Q, tau, dt, rtol_dp, info, verbosity=verb, kdim=nkmax) - call check_info(info, 'kexpm', module=this_module, procedure='test_kexptA_rdp') + call kexpm_var_dt(Xkryl, A, Q, tau, dt, rtol_dp, info, verbosity=verb, kdim=nkmax) + call check_info(info, 'kexpm_var_dt', module=this_module, procedure='test_kexptA_variable_dt_rdp') ! Check result. call Xkryl%sub(Xref) ; err = Xkryl%norm() - call get_err_str(msg, "max err: ", err) + call get_err_str(msg, "max err: ", err, rtol_dp) call check(error, err < rtol_dp) - call check_test(error, 'test_kexptA_rdp', & - & info='Comparison with matrix exponential', context=msg) + call check_test(error, 'test_kexptA_variable_dt_rdp', & + & eq='Comparison with dense matrix exponential', context=msg) return end subroutine test_kexptA_variable_dt_rdp @@ -524,7 +520,7 @@ subroutine test_block_kexptA_rdp(error) ! Misc. integer :: i, j real(dp), allocatable :: Xdata(:, :), Qdata(:, :) - real(dp) :: errv, dt + real(dp) :: errv real(dp) :: err(p, p) character(len=256) :: msg @@ -550,14 +546,12 @@ subroutine test_block_kexptA_rdp(error) ! Compute Krylov matrix exponential using sequential arnoldi method for each input column do i = 1,p - dt = tau - call kexpm(C(i), A, B(i), tau, dt, tol, info, verbosity=verb, kdim=nkmax) + call kexpm(C(i), A, B(i), tau, tol, info, verbosity=verb, kdim=nkmax) call check_info(info, 'kexpm', module=this_module, procedure='test_block_kexptA_rdp, 1') end do ! Compute Krylov matrix exponential using block-arnoldi method - dt = tau - call kexpm(Cblk, A, B, tau, dt, tol, info, verbosity=verb, kdim=nkmax) + call kexpm(Cblk, A, B, tau, tol, info, verbosity=verb, kdim=nkmax) call check_info(info, 'kexpm', module=this_module, procedure='test_block_kexptA_rdp, 2') do i = 1, p @@ -573,7 +567,7 @@ subroutine test_block_kexptA_rdp(error) enddo enddo errv = norm2(abs(err)) - !call get_err_str(msg, "max err: ", errv) + !call get_err_str(msg, "max err: ", errv, rtol_dp) do i = 1, size(Cblk) do j = 1, size(Cblk) @@ -581,7 +575,7 @@ subroutine test_block_kexptA_rdp(error) enddo enddo errv = norm2(abs(err)) - call get_err_str(msg, "max err: ", errv) + call get_err_str(msg, "max err: ", errv, rtol_dp) call check(error, errv < rtol_dp) call check_test(error, 'test_block_kexptA_rdp', & @@ -611,7 +605,7 @@ subroutine test_block_kexptA_variable_dt_rdp(error) ! Test parameters integer , parameter :: nkmax = 10 integer , parameter :: p = 3 - real(dp), parameter :: tau = 0.1_dp + real(dp), parameter :: tau = 0.5_dp real(dp), parameter :: tol = rtol_dp ! Misc. integer :: i, j @@ -620,7 +614,7 @@ subroutine test_block_kexptA_variable_dt_rdp(error) real(dp) :: err(p, p) character(len=256) :: msg - logical, parameter :: verb = .true. + logical, parameter :: verb = .false. allocate(Adata(kdim, kdim)) ; Adata = 0.0_dp allocate(Edata(kdim, kdim)) ; Edata = 0.0_dp @@ -641,16 +635,16 @@ subroutine test_block_kexptA_variable_dt_rdp(error) call expm(Edata, tau*Adata) ; Xdata = matmul(Edata,Qdata) ; call put_data(Cref, Xdata) ! Compute Krylov matrix exponential using sequential arnoldi method for each input column + dt = tau do i = 1,p - dt = tau - call kexpm(C(i), A, B(i), tau, dt, tol, info, verbosity=verb, kdim=nkmax) - call check_info(info, 'kexpm', module=this_module, procedure='test_block_kexptA_rdp, 1') + call kexpm_var_dt(C(i), A, B(i), tau, dt, tol, info, verbosity=verb, kdim=nkmax) + call check_info(info, 'kexpm_var_dt', module=this_module, procedure='test_block_kexptA_variable_dt_rdp, 1') end do ! Compute Krylov matrix exponential using block-arnoldi method dt = tau - call kexpm(Cblk, A, B, tau, dt, tol, info, verbosity=verb, kdim=nkmax) - call check_info(info, 'kexpm', module=this_module, procedure='test_block_kexptA_rdp, 2') + call kexpm_var_dt(Cblk, A, B, tau, dt, tol, info, verbosity=verb, kdim=nkmax) + call check_info(info, 'kexpm_var_dt', module=this_module, procedure='test_block_kexptA_variable_dt_rdp, 2') do i = 1, p write(output_unit, *) C(i)%norm(), Cblk(i)%norm() @@ -665,7 +659,7 @@ subroutine test_block_kexptA_variable_dt_rdp(error) enddo enddo errv = norm2(abs(err)) - !call get_err_str(msg, "max err: ", errv) + !call get_err_str(msg, "max err: ", errv, rtol_dp) do i = 1, size(Cblk) do j = 1, size(Cblk) @@ -673,11 +667,11 @@ subroutine test_block_kexptA_variable_dt_rdp(error) enddo enddo errv = norm2(abs(err)) - call get_err_str(msg, "max err: ", errv) + call get_err_str(msg, "max err: ", errv, rtol_dp) call check(error, errv < rtol_dp) - call check_test(error, 'test_block_kexptA_rdp', & - & info='Comparison with matrix exponential', context=msg) + call check_test(error, 'test_block_kexptA_variable_dt_rdp', & + & eq='Comparison with matrix exponential', context=msg) return end subroutine test_block_kexptA_variable_dt_rdp @@ -724,7 +718,7 @@ subroutine test_dense_expm_csp(error) ! Check result. err = maxval(abs(E-Eref)) - call get_err_str(msg, "max err: ", err) + call get_err_str(msg, "max err: ", err, rtol_sp) call check(error, err < rtol_sp) call check_test(error, 'test_dense_expm_csp', eq='maxval(abs(E-Eref))', context=msg) @@ -745,7 +739,7 @@ subroutine test_kexptA_csp(error) complex(sp), allocatable :: E(:, :) real(sp), parameter :: tau = 0.1_sp integer, parameter :: nkmax = 50 - real(sp) :: err, dt + real(sp) :: err integer :: info character(len=256) :: msg @@ -762,13 +756,12 @@ subroutine test_kexptA_csp(error) Xref%data = matmul(E, Q%data) ! Krylov exponential. - dt = tau - call kexpm(Xkryl, A, Q, tau, dt, rtol_sp, info, verbosity=verb, kdim=nkmax) + call kexpm(Xkryl, A, Q, tau, rtol_sp, info, verbosity=verb, kdim=nkmax) call check_info(info, 'kexpm', module=this_module, procedure='test_kexptA_csp') ! Check result. call Xkryl%sub(Xref) ; err = Xkryl%norm() - call get_err_str(msg, "max err: ", err) + call get_err_str(msg, "max err: ", err, rtol_sp) call check(error, err < rtol_sp) call check_test(error, 'test_kexptA_csp', & & eq='Comparison with dense matrix exponential', context=msg) @@ -788,13 +781,13 @@ subroutine test_kexptA_variable_dt_csp(error) ! ----- Internal variables ----- complex(sp), allocatable :: E(:, :) - real(sp), parameter :: tau = 0.1_sp + real(sp), parameter :: tau = 0.5_sp integer, parameter :: nkmax = 10 real(sp) :: err, dt integer :: info character(len=256) :: msg - logical, parameter :: verb = .true. + logical, parameter :: verb = .false. ! Initialize data. A = linop_csp() ; call init_rand(A) @@ -808,15 +801,15 @@ subroutine test_kexptA_variable_dt_csp(error) ! Krylov exponential. dt = tau - call kexpm(Xkryl, A, Q, tau, dt, rtol_sp, info, verbosity=verb, kdim=nkmax) - call check_info(info, 'kexpm', module=this_module, procedure='test_kexptA_csp') + call kexpm_var_dt(Xkryl, A, Q, tau, dt, rtol_sp, info, verbosity=verb, kdim=nkmax) + call check_info(info, 'kexpm_var_dt', module=this_module, procedure='test_kexptA_variable_dt_csp') ! Check result. call Xkryl%sub(Xref) ; err = Xkryl%norm() - call get_err_str(msg, "max err: ", err) + call get_err_str(msg, "max err: ", err, rtol_sp) call check(error, err < rtol_sp) - call check_test(error, 'test_kexptA_csp', & - & info='Comparison with matrix exponential', context=msg) + call check_test(error, 'test_kexptA_variable_dt_csp', & + & eq='Comparison with dense matrix exponential', context=msg) return end subroutine test_kexptA_variable_dt_csp @@ -847,7 +840,7 @@ subroutine test_block_kexptA_csp(error) ! Misc. integer :: i, j complex(sp), allocatable :: Xdata(:, :), Qdata(:, :) - real(sp) :: errv, dt + real(sp) :: errv complex(sp) :: err(p, p) character(len=256) :: msg @@ -873,14 +866,12 @@ subroutine test_block_kexptA_csp(error) ! Compute Krylov matrix exponential using sequential arnoldi method for each input column do i = 1,p - dt = tau - call kexpm(C(i), A, B(i), tau, dt, tol, info, verbosity=verb, kdim=nkmax) + call kexpm(C(i), A, B(i), tau, tol, info, verbosity=verb, kdim=nkmax) call check_info(info, 'kexpm', module=this_module, procedure='test_block_kexptA_csp, 1') end do ! Compute Krylov matrix exponential using block-arnoldi method - dt = tau - call kexpm(Cblk, A, B, tau, dt, tol, info, verbosity=verb, kdim=nkmax) + call kexpm(Cblk, A, B, tau, tol, info, verbosity=verb, kdim=nkmax) call check_info(info, 'kexpm', module=this_module, procedure='test_block_kexptA_csp, 2') do i = 1, p @@ -896,7 +887,7 @@ subroutine test_block_kexptA_csp(error) enddo enddo errv = norm2(abs(err)) - !call get_err_str(msg, "max err: ", errv) + !call get_err_str(msg, "max err: ", errv, rtol_sp) do i = 1, size(Cblk) do j = 1, size(Cblk) @@ -904,7 +895,7 @@ subroutine test_block_kexptA_csp(error) enddo enddo errv = norm2(abs(err)) - call get_err_str(msg, "max err: ", errv) + call get_err_str(msg, "max err: ", errv, rtol_sp) call check(error, errv < rtol_sp) call check_test(error, 'test_block_kexptA_csp', & @@ -934,7 +925,7 @@ subroutine test_block_kexptA_variable_dt_csp(error) ! Test parameters integer , parameter :: nkmax = 10 integer , parameter :: p = 3 - real(sp), parameter :: tau = 0.1_sp + real(sp), parameter :: tau = 0.5_sp real(sp), parameter :: tol = rtol_sp ! Misc. integer :: i, j @@ -943,7 +934,7 @@ subroutine test_block_kexptA_variable_dt_csp(error) complex(sp) :: err(p, p) character(len=256) :: msg - logical, parameter :: verb = .true. + logical, parameter :: verb = .false. allocate(Adata(kdim, kdim)) ; Adata = 0.0_sp allocate(Edata(kdim, kdim)) ; Edata = 0.0_sp @@ -964,16 +955,16 @@ subroutine test_block_kexptA_variable_dt_csp(error) call expm(Edata, tau*Adata) ; Xdata = matmul(Edata,Qdata) ; call put_data(Cref, Xdata) ! Compute Krylov matrix exponential using sequential arnoldi method for each input column + dt = tau do i = 1,p - dt = tau - call kexpm(C(i), A, B(i), tau, dt, tol, info, verbosity=verb, kdim=nkmax) - call check_info(info, 'kexpm', module=this_module, procedure='test_block_kexptA_csp, 1') + call kexpm_var_dt(C(i), A, B(i), tau, dt, tol, info, verbosity=verb, kdim=nkmax) + call check_info(info, 'kexpm_var_dt', module=this_module, procedure='test_block_kexptA_variable_dt_csp, 1') end do ! Compute Krylov matrix exponential using block-arnoldi method dt = tau - call kexpm(Cblk, A, B, tau, dt, tol, info, verbosity=verb, kdim=nkmax) - call check_info(info, 'kexpm', module=this_module, procedure='test_block_kexptA_csp, 2') + call kexpm_var_dt(Cblk, A, B, tau, dt, tol, info, verbosity=verb, kdim=nkmax) + call check_info(info, 'kexpm_var_dt', module=this_module, procedure='test_block_kexptA_variable_dt_csp, 2') do i = 1, p write(output_unit, *) C(i)%norm(), Cblk(i)%norm() @@ -988,7 +979,7 @@ subroutine test_block_kexptA_variable_dt_csp(error) enddo enddo errv = norm2(abs(err)) - !call get_err_str(msg, "max err: ", errv) + !call get_err_str(msg, "max err: ", errv, rtol_sp) do i = 1, size(Cblk) do j = 1, size(Cblk) @@ -996,11 +987,11 @@ subroutine test_block_kexptA_variable_dt_csp(error) enddo enddo errv = norm2(abs(err)) - call get_err_str(msg, "max err: ", errv) + call get_err_str(msg, "max err: ", errv, rtol_sp) call check(error, errv < rtol_sp) - call check_test(error, 'test_block_kexptA_csp', & - & info='Comparison with matrix exponential', context=msg) + call check_test(error, 'test_block_kexptA_variable_dt_csp', & + & eq='Comparison with matrix exponential', context=msg) return end subroutine test_block_kexptA_variable_dt_csp @@ -1047,7 +1038,7 @@ subroutine test_dense_expm_cdp(error) ! Check result. err = maxval(abs(E-Eref)) - call get_err_str(msg, "max err: ", err) + call get_err_str(msg, "max err: ", err, rtol_dp) call check(error, err < rtol_dp) call check_test(error, 'test_dense_expm_cdp', eq='maxval(abs(E-Eref))', context=msg) @@ -1068,7 +1059,7 @@ subroutine test_kexptA_cdp(error) complex(dp), allocatable :: E(:, :) real(dp), parameter :: tau = 0.1_dp integer, parameter :: nkmax = 50 - real(dp) :: err, dt + real(dp) :: err integer :: info character(len=256) :: msg @@ -1085,13 +1076,12 @@ subroutine test_kexptA_cdp(error) Xref%data = matmul(E, Q%data) ! Krylov exponential. - dt = tau - call kexpm(Xkryl, A, Q, tau, dt, rtol_dp, info, verbosity=verb, kdim=nkmax) + call kexpm(Xkryl, A, Q, tau, rtol_dp, info, verbosity=verb, kdim=nkmax) call check_info(info, 'kexpm', module=this_module, procedure='test_kexptA_cdp') ! Check result. call Xkryl%sub(Xref) ; err = Xkryl%norm() - call get_err_str(msg, "max err: ", err) + call get_err_str(msg, "max err: ", err, rtol_dp) call check(error, err < rtol_dp) call check_test(error, 'test_kexptA_cdp', & & eq='Comparison with dense matrix exponential', context=msg) @@ -1111,13 +1101,13 @@ subroutine test_kexptA_variable_dt_cdp(error) ! ----- Internal variables ----- complex(dp), allocatable :: E(:, :) - real(dp), parameter :: tau = 0.1_dp + real(dp), parameter :: tau = 0.5_dp integer, parameter :: nkmax = 10 real(dp) :: err, dt integer :: info character(len=256) :: msg - logical, parameter :: verb = .true. + logical, parameter :: verb = .false. ! Initialize data. A = linop_cdp() ; call init_rand(A) @@ -1131,15 +1121,15 @@ subroutine test_kexptA_variable_dt_cdp(error) ! Krylov exponential. dt = tau - call kexpm(Xkryl, A, Q, tau, dt, rtol_dp, info, verbosity=verb, kdim=nkmax) - call check_info(info, 'kexpm', module=this_module, procedure='test_kexptA_cdp') + call kexpm_var_dt(Xkryl, A, Q, tau, dt, rtol_dp, info, verbosity=verb, kdim=nkmax) + call check_info(info, 'kexpm_var_dt', module=this_module, procedure='test_kexptA_variable_dt_cdp') ! Check result. call Xkryl%sub(Xref) ; err = Xkryl%norm() - call get_err_str(msg, "max err: ", err) + call get_err_str(msg, "max err: ", err, rtol_dp) call check(error, err < rtol_dp) - call check_test(error, 'test_kexptA_cdp', & - & info='Comparison with matrix exponential', context=msg) + call check_test(error, 'test_kexptA_variable_dt_cdp', & + & eq='Comparison with dense matrix exponential', context=msg) return end subroutine test_kexptA_variable_dt_cdp @@ -1170,7 +1160,7 @@ subroutine test_block_kexptA_cdp(error) ! Misc. integer :: i, j complex(dp), allocatable :: Xdata(:, :), Qdata(:, :) - real(dp) :: errv, dt + real(dp) :: errv complex(dp) :: err(p, p) character(len=256) :: msg @@ -1196,14 +1186,12 @@ subroutine test_block_kexptA_cdp(error) ! Compute Krylov matrix exponential using sequential arnoldi method for each input column do i = 1,p - dt = tau - call kexpm(C(i), A, B(i), tau, dt, tol, info, verbosity=verb, kdim=nkmax) + call kexpm(C(i), A, B(i), tau, tol, info, verbosity=verb, kdim=nkmax) call check_info(info, 'kexpm', module=this_module, procedure='test_block_kexptA_cdp, 1') end do ! Compute Krylov matrix exponential using block-arnoldi method - dt = tau - call kexpm(Cblk, A, B, tau, dt, tol, info, verbosity=verb, kdim=nkmax) + call kexpm(Cblk, A, B, tau, tol, info, verbosity=verb, kdim=nkmax) call check_info(info, 'kexpm', module=this_module, procedure='test_block_kexptA_cdp, 2') do i = 1, p @@ -1219,7 +1207,7 @@ subroutine test_block_kexptA_cdp(error) enddo enddo errv = norm2(abs(err)) - !call get_err_str(msg, "max err: ", errv) + !call get_err_str(msg, "max err: ", errv, rtol_dp) do i = 1, size(Cblk) do j = 1, size(Cblk) @@ -1227,7 +1215,7 @@ subroutine test_block_kexptA_cdp(error) enddo enddo errv = norm2(abs(err)) - call get_err_str(msg, "max err: ", errv) + call get_err_str(msg, "max err: ", errv, rtol_dp) call check(error, errv < rtol_dp) call check_test(error, 'test_block_kexptA_cdp', & @@ -1257,7 +1245,7 @@ subroutine test_block_kexptA_variable_dt_cdp(error) ! Test parameters integer , parameter :: nkmax = 10 integer , parameter :: p = 3 - real(dp), parameter :: tau = 0.1_dp + real(dp), parameter :: tau = 0.5_dp real(dp), parameter :: tol = rtol_dp ! Misc. integer :: i, j @@ -1266,7 +1254,7 @@ subroutine test_block_kexptA_variable_dt_cdp(error) complex(dp) :: err(p, p) character(len=256) :: msg - logical, parameter :: verb = .true. + logical, parameter :: verb = .false. allocate(Adata(kdim, kdim)) ; Adata = 0.0_dp allocate(Edata(kdim, kdim)) ; Edata = 0.0_dp @@ -1287,16 +1275,16 @@ subroutine test_block_kexptA_variable_dt_cdp(error) call expm(Edata, tau*Adata) ; Xdata = matmul(Edata,Qdata) ; call put_data(Cref, Xdata) ! Compute Krylov matrix exponential using sequential arnoldi method for each input column + dt = tau do i = 1,p - dt = tau - call kexpm(C(i), A, B(i), tau, dt, tol, info, verbosity=verb, kdim=nkmax) - call check_info(info, 'kexpm', module=this_module, procedure='test_block_kexptA_cdp, 1') + call kexpm_var_dt(C(i), A, B(i), tau, dt, tol, info, verbosity=verb, kdim=nkmax) + call check_info(info, 'kexpm_var_dt', module=this_module, procedure='test_block_kexptA_variable_dt_cdp, 1') end do ! Compute Krylov matrix exponential using block-arnoldi method dt = tau - call kexpm(Cblk, A, B, tau, dt, tol, info, verbosity=verb, kdim=nkmax) - call check_info(info, 'kexpm', module=this_module, procedure='test_block_kexptA_cdp, 2') + call kexpm_var_dt(Cblk, A, B, tau, dt, tol, info, verbosity=verb, kdim=nkmax) + call check_info(info, 'kexpm_var_dt', module=this_module, procedure='test_block_kexptA_variable_dt_cdp, 2') do i = 1, p write(output_unit, *) C(i)%norm(), Cblk(i)%norm() @@ -1311,7 +1299,7 @@ subroutine test_block_kexptA_variable_dt_cdp(error) enddo enddo errv = norm2(abs(err)) - !call get_err_str(msg, "max err: ", errv) + !call get_err_str(msg, "max err: ", errv, rtol_dp) do i = 1, size(Cblk) do j = 1, size(Cblk) @@ -1319,11 +1307,11 @@ subroutine test_block_kexptA_variable_dt_cdp(error) enddo enddo errv = norm2(abs(err)) - call get_err_str(msg, "max err: ", errv) + call get_err_str(msg, "max err: ", errv, rtol_dp) call check(error, errv < rtol_dp) - call check_test(error, 'test_block_kexptA_cdp', & - & info='Comparison with matrix exponential', context=msg) + call check_test(error, 'test_block_kexptA_variable_dt_cdp', & + & eq='Comparison with matrix exponential', context=msg) return end subroutine test_block_kexptA_variable_dt_cdp @@ -1377,7 +1365,7 @@ subroutine test_dense_sqrtm_pos_def_rsp(error) call check_info(info, 'sqrtm', module=this_module, procedure='test_dense_sqrtm_pos_def_rsp') err = maxval(matmul(sqrtmA, sqrtmA) - A) - call get_err_str(msg, "max err: ", err) + call get_err_str(msg, "max err: ", err, rtol_sp) call check(error, err < rtol_sp) call check_test(error, 'test_dense_sqrtm_pos_def_rsp', eq='sqrt(A)**2 = A', context=msg) @@ -1419,7 +1407,7 @@ subroutine test_dense_sqrtm_pos_semi_def_rsp(error) call check_info(info, 'sqrtm', module=this_module, procedure='test_dense_sqrtm_pos_semi_def_rsp') err = maxval(matmul(sqrtmA, sqrtmA) - A) - call get_err_str(msg, "max err: ", err) + call get_err_str(msg, "max err: ", err, rtol_sp) call check(error, err < rtol_sp) call check_test(error, 'test_dense_sqrtm_pos_semi_def_rsp', eq='sqrt(A)**2 = A', context=trim(msg)) @@ -1471,7 +1459,7 @@ subroutine test_dense_sqrtm_pos_def_rdp(error) call check_info(info, 'sqrtm', module=this_module, procedure='test_dense_sqrtm_pos_def_rdp') err = maxval(matmul(sqrtmA, sqrtmA) - A) - call get_err_str(msg, "max err: ", err) + call get_err_str(msg, "max err: ", err, rtol_dp) call check(error, err < rtol_dp) call check_test(error, 'test_dense_sqrtm_pos_def_rdp', eq='sqrt(A)**2 = A', context=msg) @@ -1513,7 +1501,7 @@ subroutine test_dense_sqrtm_pos_semi_def_rdp(error) call check_info(info, 'sqrtm', module=this_module, procedure='test_dense_sqrtm_pos_semi_def_rdp') err = maxval(matmul(sqrtmA, sqrtmA) - A) - call get_err_str(msg, "max err: ", err) + call get_err_str(msg, "max err: ", err, rtol_dp) call check(error, err < rtol_dp) call check_test(error, 'test_dense_sqrtm_pos_semi_def_rdp', eq='sqrt(A)**2 = A', context=trim(msg)) @@ -1566,7 +1554,7 @@ subroutine test_dense_sqrtm_pos_def_csp(error) call check_info(info, 'sqrtm', module=this_module, procedure='test_dense_sqrtm_pos_def_csp') err = maxval(abs(matmul(sqrtmA, sqrtmA) - A)) - call get_err_str(msg, "max err: ", err) + call get_err_str(msg, "max err: ", err, rtol_sp) call check(error, err < rtol_sp) call check_test(error, 'test_dense_sqrtm_pos_def_csp', eq='sqrt(A)**2 = A', context=msg) @@ -1609,7 +1597,7 @@ subroutine test_dense_sqrtm_pos_semi_def_csp(error) call check_info(info, 'sqrtm', module=this_module, procedure='test_dense_sqrtm_pos_semi_def_csp') err = maxval(abs(matmul(sqrtmA, sqrtmA) - A)) - call get_err_str(msg, "max err: ", err) + call get_err_str(msg, "max err: ", err, rtol_sp) call check(error, err < rtol_sp) call check_test(error, 'test_dense_sqrtm_pos_semi_def_csp', eq='sqrt(A)**2 = A', context=trim(msg)) @@ -1662,7 +1650,7 @@ subroutine test_dense_sqrtm_pos_def_cdp(error) call check_info(info, 'sqrtm', module=this_module, procedure='test_dense_sqrtm_pos_def_cdp') err = maxval(abs(matmul(sqrtmA, sqrtmA) - A)) - call get_err_str(msg, "max err: ", err) + call get_err_str(msg, "max err: ", err, rtol_dp) call check(error, err < rtol_dp) call check_test(error, 'test_dense_sqrtm_pos_def_cdp', eq='sqrt(A)**2 = A', context=msg) @@ -1705,7 +1693,7 @@ subroutine test_dense_sqrtm_pos_semi_def_cdp(error) call check_info(info, 'sqrtm', module=this_module, procedure='test_dense_sqrtm_pos_semi_def_cdp') err = maxval(abs(matmul(sqrtmA, sqrtmA) - A)) - call get_err_str(msg, "max err: ", err) + call get_err_str(msg, "max err: ", err, rtol_dp) call check(error, err < rtol_dp) call check_test(error, 'test_dense_sqrtm_pos_semi_def_cdp', eq='sqrt(A)**2 = A', context=trim(msg)) diff --git a/test/TestExpmlib.fypp b/test/TestExpmlib.fypp index 8a5b695..da0b9b9 100644 --- a/test/TestExpmlib.fypp +++ b/test/TestExpmlib.fypp @@ -79,7 +79,7 @@ contains ! Check result. err = maxval(abs(E-Eref)) - call get_err_str(msg, "max err: ", err) + call get_err_str(msg, "max err: ", err, rtol_${kind}$) call check(error, err < rtol_${kind}$) call check_test(error, 'test_dense_expm_${type[0]}$${kind}$', eq='maxval(abs(E-Eref))', context=msg) @@ -100,7 +100,7 @@ contains ${type}$, allocatable :: E(:, :) real(${kind}$), parameter :: tau = 0.1_${kind}$ integer, parameter :: nkmax = 50 - real(${kind}$) :: err, dt + real(${kind}$) :: err integer :: info character(len=256) :: msg @@ -117,13 +117,12 @@ contains Xref%data = matmul(E, Q%data) ! Krylov exponential. - dt = tau - call kexpm(Xkryl, A, Q, tau, dt, rtol_${kind}$, info, verbosity=verb, kdim=nkmax) + call kexpm(Xkryl, A, Q, tau, rtol_${kind}$, info, verbosity=verb, kdim=nkmax) call check_info(info, 'kexpm', module=this_module, procedure='test_kexptA_${type[0]}$${kind}$') ! Check result. call Xkryl%sub(Xref) ; err = Xkryl%norm() - call get_err_str(msg, "max err: ", err) + call get_err_str(msg, "max err: ", err, rtol_${kind}$) call check(error, err < rtol_${kind}$) call check_test(error, 'test_kexptA_${type[0]}$${kind}$', & & eq='Comparison with dense matrix exponential', context=msg) @@ -143,13 +142,13 @@ contains ! ----- Internal variables ----- ${type}$, allocatable :: E(:, :) - real(${kind}$), parameter :: tau = 0.1_${kind}$ + real(${kind}$), parameter :: tau = 0.5_${kind}$ integer, parameter :: nkmax = 10 real(${kind}$) :: err, dt integer :: info character(len=256) :: msg - logical, parameter :: verb = .true. + logical, parameter :: verb = .false. ! Initialize data. A = linop_${type[0]}$${kind}$() ; call init_rand(A) @@ -163,15 +162,15 @@ contains ! Krylov exponential. dt = tau - call kexpm(Xkryl, A, Q, tau, dt, rtol_${kind}$, info, verbosity=verb, kdim=nkmax) - call check_info(info, 'kexpm', module=this_module, procedure='test_kexptA_${type[0]}$${kind}$') + call kexpm_var_dt(Xkryl, A, Q, tau, dt, rtol_${kind}$, info, verbosity=verb, kdim=nkmax) + call check_info(info, 'kexpm_var_dt', module=this_module, procedure='test_kexptA_variable_dt_${type[0]}$${kind}$') ! Check result. call Xkryl%sub(Xref) ; err = Xkryl%norm() - call get_err_str(msg, "max err: ", err) + call get_err_str(msg, "max err: ", err, rtol_${kind}$) call check(error, err < rtol_${kind}$) - call check_test(error, 'test_kexptA_${type[0]}$${kind}$', & - & info='Comparison with matrix exponential', context=msg) + call check_test(error, 'test_kexptA_variable_dt_${type[0]}$${kind}$', & + & eq='Comparison with dense matrix exponential', context=msg) return end subroutine test_kexptA_variable_dt_${type[0]}$${kind}$ @@ -202,7 +201,7 @@ contains ! Misc. integer :: i, j ${type}$, allocatable :: Xdata(:, :), Qdata(:, :) - real(${kind}$) :: errv, dt + real(${kind}$) :: errv ${type}$ :: err(p, p) character(len=256) :: msg @@ -228,14 +227,12 @@ contains ! Compute Krylov matrix exponential using sequential arnoldi method for each input column do i = 1,p - dt = tau - call kexpm(C(i), A, B(i), tau, dt, tol, info, verbosity=verb, kdim=nkmax) + call kexpm(C(i), A, B(i), tau, tol, info, verbosity=verb, kdim=nkmax) call check_info(info, 'kexpm', module=this_module, procedure='test_block_kexptA_${type[0]}$${kind}$, 1') end do ! Compute Krylov matrix exponential using block-arnoldi method - dt = tau - call kexpm(Cblk, A, B, tau, dt, tol, info, verbosity=verb, kdim=nkmax) + call kexpm(Cblk, A, B, tau, tol, info, verbosity=verb, kdim=nkmax) call check_info(info, 'kexpm', module=this_module, procedure='test_block_kexptA_${type[0]}$${kind}$, 2') do i = 1, p @@ -251,7 +248,7 @@ contains enddo enddo errv = norm2(abs(err)) - !call get_err_str(msg, "max err: ", errv) + !call get_err_str(msg, "max err: ", errv, rtol_${kind}$) do i = 1, size(Cblk) do j = 1, size(Cblk) @@ -259,7 +256,7 @@ contains enddo enddo errv = norm2(abs(err)) - call get_err_str(msg, "max err: ", errv) + call get_err_str(msg, "max err: ", errv, rtol_${kind}$) call check(error, errv < rtol_${kind}$) call check_test(error, 'test_block_kexptA_${type[0]}$${kind}$', & @@ -289,7 +286,7 @@ contains ! Test parameters integer , parameter :: nkmax = 10 integer , parameter :: p = 3 - real(${kind}$), parameter :: tau = 0.1_${kind}$ + real(${kind}$), parameter :: tau = 0.5_${kind}$ real(${kind}$), parameter :: tol = rtol_${kind}$ ! Misc. integer :: i, j @@ -298,7 +295,7 @@ contains ${type}$ :: err(p, p) character(len=256) :: msg - logical, parameter :: verb = .true. + logical, parameter :: verb = .false. allocate(Adata(kdim, kdim)) ; Adata = 0.0_${kind}$ allocate(Edata(kdim, kdim)) ; Edata = 0.0_${kind}$ @@ -319,16 +316,16 @@ contains call expm(Edata, tau*Adata) ; Xdata = matmul(Edata,Qdata) ; call put_data(Cref, Xdata) ! Compute Krylov matrix exponential using sequential arnoldi method for each input column + dt = tau do i = 1,p - dt = tau - call kexpm(C(i), A, B(i), tau, dt, tol, info, verbosity=verb, kdim=nkmax) - call check_info(info, 'kexpm', module=this_module, procedure='test_block_kexptA_${type[0]}$${kind}$, 1') + call kexpm_var_dt(C(i), A, B(i), tau, dt, tol, info, verbosity=verb, kdim=nkmax) + call check_info(info, 'kexpm_var_dt', module=this_module, procedure='test_block_kexptA_variable_dt_${type[0]}$${kind}$, 1') end do ! Compute Krylov matrix exponential using block-arnoldi method dt = tau - call kexpm(Cblk, A, B, tau, dt, tol, info, verbosity=verb, kdim=nkmax) - call check_info(info, 'kexpm', module=this_module, procedure='test_block_kexptA_${type[0]}$${kind}$, 2') + call kexpm_var_dt(Cblk, A, B, tau, dt, tol, info, verbosity=verb, kdim=nkmax) + call check_info(info, 'kexpm_var_dt', module=this_module, procedure='test_block_kexptA_variable_dt_${type[0]}$${kind}$, 2') do i = 1, p write(output_unit, *) C(i)%norm(), Cblk(i)%norm() @@ -343,7 +340,7 @@ contains enddo enddo errv = norm2(abs(err)) - !call get_err_str(msg, "max err: ", errv) + !call get_err_str(msg, "max err: ", errv, rtol_${kind}$) do i = 1, size(Cblk) do j = 1, size(Cblk) @@ -351,11 +348,11 @@ contains enddo enddo errv = norm2(abs(err)) - call get_err_str(msg, "max err: ", errv) + call get_err_str(msg, "max err: ", errv, rtol_${kind}$) call check(error, errv < rtol_${kind}$) - call check_test(error, 'test_block_kexptA_${type[0]}$${kind}$', & - & info='Comparison with matrix exponential', context=msg) + call check_test(error, 'test_block_kexptA_variable_dt_${type[0]}$${kind}$', & + & eq='Comparison with matrix exponential', context=msg) return end subroutine test_block_kexptA_variable_dt_${type[0]}$${kind}$ @@ -427,11 +424,11 @@ contains #:if type[0] == "r" err = maxval(matmul(sqrtmA, sqrtmA) - A) - call get_err_str(msg, "max err: ", err) + call get_err_str(msg, "max err: ", err, rtol_${kind}$) call check(error, err < rtol_${kind}$) #:else err = maxval(abs(matmul(sqrtmA, sqrtmA) - A)) - call get_err_str(msg, "max err: ", err) + call get_err_str(msg, "max err: ", err, rtol_${kind}$) call check(error, err < rtol_${kind}$) #:endif call check_test(error, 'test_dense_sqrtm_pos_def_${type[0]}$${kind}$', eq='sqrt(A)**2 = A', context=msg) @@ -490,11 +487,11 @@ contains #:if type[0] == "r" err = maxval(matmul(sqrtmA, sqrtmA) - A) - call get_err_str(msg, "max err: ", err) + call get_err_str(msg, "max err: ", err, rtol_${kind}$) call check(error, err < rtol_${kind}$) #:else err = maxval(abs(matmul(sqrtmA, sqrtmA) - A)) - call get_err_str(msg, "max err: ", err) + call get_err_str(msg, "max err: ", err, rtol_${kind}$) call check(error, err < rtol_${kind}$) #:endif call check_test(error, 'test_dense_sqrtm_pos_semi_def_${type[0]}$${kind}$', eq='sqrt(A)**2 = A', context=trim(msg))