Classic Newton implementation attempt. Not functional yet.

src/AbstractSystems.f90

@@ -0,0 +1,64 @@
+module LightKrylov_AbstractSystems
+    use LightKrylov_AbstractVectors
+    use LightKrylov_AbstractLinops
+    implicit none
+    private
+
+    character*128, parameter :: this_module = 'LightKrylov_AbstractSystems'
+
+    !> General abstract type for general systems.
+    type, abstract, public :: abstract_dynamical_system
+    end type abstract_dynamical_system
+
+    !-----------------------------------------------------------
+    !-----     ABSTRACT GENERAL SYSTEM TYPE DEFINITION     -----
+    !-----------------------------------------------------------
+
+    !> Abstract Jacobian linop.
+    type, abstract, extends(abstract_linop_rdp), public :: abstract_jacobian_linop_rdp
+        class(abstract_vector_rdp), allocatable :: X
+    contains
+        private
+        procedure(abstract_jac_matvec_rdp), pass(self), deferred, public :: matvec
+        procedure(abstract_jac_matvec_rdp), pass(self), deferred, public :: rmatvec
+    end type
+
+    abstract interface
+        subroutine abstract_jac_matvec_rdp(self, vec_in, vec_out)
+            !! Interface for the matrix-vector product.
+            use lightkrylov_AbstractVectors
+            import abstract_jacobian_linop_rdp
+            class(abstract_jacobian_linop_rdp) , intent(in)  :: self
+            !! Linear operator \(\mathbf{A}\).
+            class(abstract_vector_rdp), intent(in)  :: vec_in
+            !! Vector to be multiplied by \(\mathbf{A}\).
+            class(abstract_vector_rdp), intent(out) :: vec_out
+            !! Result of the matrix-vector product.
+        end subroutine abstract_jac_matvec_rdp
+    end interface
+
+    !> Abstract continuous system.
+    type, abstract, extends(abstract_dynamical_system), public :: abstract_system_rdp
+        class(abstract_jacobian_linop_rdp), allocatable :: jacobian
+        !! System Jacobian \( \left. \frac{\partial \mathbf{F}}{\partial \mathbf{X}} \right|_{X^*} \).
+    contains
+        private
+        procedure(abstract_eval_rdp), pass(self), deferred, public :: eval
+        !! Procedure to evaluate the system response \( \mathbf{Y} = \mathbf{F}(\mathbf{X}) \).
+    end type
+
+    abstract interface
+        subroutine abstract_eval_rdp(self, vec_in, vec_out)
+            !! Interface for the evaluation of the system response.
+            use LightKrylov_AbstractVectors
+            import abstract_system_rdp
+            class(abstract_system_rdp), intent(in)  :: self
+            !! System
+            class(abstract_vector_rdp), intent(in)  :: vec_in
+            !! Base state
+            class(abstract_vector_rdp), intent(out) :: vec_out
+            !! System response
+        end subroutine abstract_eval_rdp
+    end interface
+
+end module LightKrylov_AbstractSystems

src/LightKrylov.f90

@@ -7,6 +7,8 @@ module LightKrylov
     use LightKrylov_AbstractVectors
     ! --> Definitions of the abstract linear operators.
     use LightKrylov_AbstractLinops
+    ! --> Definitions of the abstract dynamical systems.
+    use LightKrylov_AbstractSystems
     ! --> Standard Krylov techniques.
     use LightKrylov_BaseKrylov
     ! --> Iterative solvers.
@@ -68,6 +70,10 @@ module LightKrylov
     public :: scaled_linop_cdp
     public :: axpby_linop_cdp
     public :: abstract_hermitian_linop_cdp
+
+    ! AbstractSystems exports.
+    public :: abstract_system_rdp
+    public :: abstract_jacobian_linop_rdp
 
     ! BaseKrylov exports.
     public :: qr

src/LightKrylov.fypp

@@ -9,6 +9,8 @@ module LightKrylov
     use LightKrylov_AbstractVectors
     ! --> Definitions of the abstract linear operators.
     use LightKrylov_AbstractLinops
+    ! --> Definitions of the abstract dynamical systems.
+    use LightKrylov_AbstractSystems
     ! --> Standard Krylov techniques.
     use LightKrylov_BaseKrylov
     ! --> Iterative solvers.
@@ -57,6 +59,10 @@ module LightKrylov
     public :: abstract_hermitian_linop_${type[0]}$${kind}$
     #:endif
     #:endfor
+
+    ! AbstractSystems exports.
+    public :: abstract_system_rdp
+    public :: abstract_jacobian_linop_rdp
 
     ! BaseKrylov exports.
     public :: qr

src/NewtonKrylov.f90

@@ -0,0 +1,91 @@
+module LightKrylov_NewtonKrylov
+    use stdlib_optval, only: optval
+    use LightKrylov_Constants
+    use LightKrylov_Logger
+    use LightKrylov_AbstractVectors
+    use LightKrylov_AbstractLinops
+    use LightKrylov_AbstractSystems
+    use LightKrylov_IterativeSolvers
+
+    implicit none
+    private
+
+    character*128, parameter :: this_module = 'LightKrylov_NewtonKrylov'
+
+    public :: newton
+
+    interface newton
+       module procedure newton_classic_rdp
+    end interface
+
+contains
+
+    subroutine newton_classic_rdp(sys, X0, tol, verb, info)
+       !! Classic no-frills implementation of the Newton-Krylov root-finding algorithm
+       class(abstract_system_rdp), intent(inout)  :: sys
+       !! Dynamical system for which we wish to compute a fixed point
+       class(abstract_vector_rdp), intent(inout)  :: X0
+       !! Initial guess for the fixed point, will be overwritten with solution
+       real(dp),                   intent(in)     :: tol
+       !! Solver tolerance
+       logical, optional,          intent(in)     :: verb
+       !! verbosity toggle
+       integer,                    intent(out)    :: info
+       !! Information flag
+
+       !--------------------------------------
+       !-----     Internal variables     -----
+       !--------------------------------------
+       ! residual vector
+       class(abstract_vector_rdp), allocatable :: residual, X, dX
+       real(dp) :: rnorm
+       logical :: converged, verb_
+       integer :: i, maxiter
+
+       ! Optional parameters
+       verb_   = optval(verb, .false.)
+
+       ! Initialisation
+       maxiter = 100
+       converged = .false.
+       allocate(residual, source=X0); call residual%zero()
+       allocate(dX,       source=X0); call dX%zero()
+       allocate(X,        source=X0);
+
+       ! Newton iteration
+       newton: do i = 1, maxiter
+
+          call sys%eval(X0, X)
+          rnorm = X%norm()
+          if (verb) write(*,*) "Iteration", i, ": Residual norm = ", rnorm
+
+          ! Check for convergence.
+          if (rnorm < tol) then
+            if (verb) write(*,*) "Convergence achieved."
+            converged = .true.
+            exit newton
+          end if
+
+          ! Define the Jacobian
+          call sys%jacobian%X%axpby(zero_rdp, X, one_rdp)
+
+          ! Solve the linear system using GMRES.
+          call residual%chsgn()
+          call gmres(sys%jacobian, residual, dX, info)
+          call check_info(info, 'gmres', module=this_module, procedure='newton_classic_rdp')
+
+          ! Update the solution and overwrite X0
+          call X0%axpby(zero_rdp, X, one_rdp)
+          call X0%add(dX)
+
+       enddo newton
+
+       if (.not.converged) then
+         if (verb) write(*,*) 'Newton iteration did not converge within', maxiter, 'steps.'
+          info = -1
+       endif
+
+       return
+    end subroutine newton_classic_rdp
+
+end module LightKrylov_NewtonKrylov

src/TestUtils.f90

@@ -142,6 +142,40 @@ module LightKrylov_TestUtils
     end type
 
 
+    ! ROESSLER SYSTEM
+
+    ! Mesh related parameters.
+    real(dp), parameter :: a = 0.2
+    real(dp), parameter :: b = 0.2
+    real(dp), parameter :: c = 5.7
+
+    type, extends(abstract_vector_rdp), public :: state_vector
+       real(dp) :: x = 0.0_dp
+       real(dp) :: y = 0.0_dp
+       real(dp) :: z = 0.0_dp
+    contains
+       private
+       procedure, pass(self), public :: zero
+       procedure, pass(self), public :: dot
+       procedure, pass(self), public :: scal
+       procedure, pass(self), public :: axpby
+       procedure, pass(self), public :: rand
+       procedure, pass(self), public :: get_size
+    end type state_vector
+
+    type, extends(abstract_system_rdp), public :: roessler
+    contains
+       private
+       procedure, pass(self), public :: eval => eval_roessler
+    end type roessler
+
+    type, extends(abstract_jacobian_linop_rdp), public :: jacobian
+    contains
+       private
+       procedure, pass(self), public :: matvec => lin_roessler
+       procedure, pass(self), public :: rmatvec => adj_lin_roessler
+    end type jacobian
+
     interface get_data
         module procedure get_data_vec_rsp
         module procedure get_data_vec_basis_rsp
@@ -1065,5 +1099,136 @@ subroutine get_err_str_dp(msg, info, err)
       msg = indent // info // value_str // achar(10)
 
     end subroutine get_err_str_dp
+
+    !-----------------------------------------------
+    !-----     ROESSLER SYSTEM TYPE DEFINITION     -----
+    !-----------------------------------------------
+
+    subroutine zero(self)
+      class(state_vector), intent(inout) :: self
+      self%x = 0.0_dp
+      self%y = 0.0_dp
+      self%z = 0.0_dp
+      return
+    end subroutine zero
+
+    real(dp) function dot(self, vec) result(alpha)
+      class(state_vector)       , intent(in) :: self
+      class(abstract_vector_rdp), intent(in) :: vec
+      select type(vec)
+      type is(state_vector)
+         alpha = self%x*vec%x + self%y*vec%y + self%z*vec%z
+      end select
+      return
+    end function dot
+
+    subroutine scal(self, alpha)
+      class(state_vector), intent(inout) :: self
+      real(dp)           , intent(in)    :: alpha
+      self%x = self%x * alpha
+      self%y = self%y * alpha
+      self%z = self%z * alpha
+      return
+    end subroutine scal
+
+    subroutine axpby(self, alpha, vec, beta)
+      class(state_vector)       , intent(inout) :: self
+      class(abstract_vector_rdp), intent(in)    :: vec
+      real(dp)                  , intent(in)    :: alpha, beta
+      select type(vec)
+      type is(state_vector)
+         self%x = alpha*self%x + beta*vec%x
+         self%y = alpha*self%y + beta*vec%y
+         self%z = alpha*self%z + beta*vec%z
+      end select
+      return
+    end subroutine axpby
+
+    integer function get_size(self) result(N)
+      class(state_vector), intent(in) :: self
+      N = 3
+      return
+    end function get_size
+
+    subroutine rand(self, ifnorm)
+      class(state_vector), intent(inout) :: self
+      logical, optional,   intent(in)    :: ifnorm
+      logical :: normalized
+      real(dp) :: mu, var
+      real(dp) :: alpha
+
+      mu = 0.0_dp
+      var = 1.0_dp
+      self%x = normal(mu, var)
+      self%y = normal(mu, var)
+      self%z = normal(mu, var)
+
+      normalized = optval(ifnorm, .false.)
+      if (normalized) then
+         alpha = self%norm()
+         call self%scal(1.0_dp/alpha)
+      endif
+      return
+    end subroutine rand
+
+    subroutine eval_roessler(self, vec_in, vec_out)
+      class(roessler),            intent(in)  :: self
+      class(abstract_vector_rdp), intent(in)  :: vec_in
+      class(abstract_vector_rdp), intent(out) :: vec_out
+
+      select type(vec_in)
+        type is(state_vector)
+            select type(vec_out)
+            type is(state_vector)
+
+            vec_out%x = -vec_in%y - vec_in%z
+            vec_out%y = vec_in%x + a * vec_in%y
+            vec_out%z = b + vec_in%z * (vec_in%x - c)
+
+            end select
+        end select
+
+      return
+    end subroutine eval_roessler
+
+    subroutine lin_roessler(self, vec_in, vec_out)
+      class(jacobian),            intent(in)  :: self
+      class(abstract_vector_rdp), intent(in)  :: vec_in
+      class(abstract_vector_rdp), intent(out) :: vec_out
+
+      select type(vec_in)
+      type is(state_vector)
+         select type(vec_out)
+         type is(state_vector)
+
+            vec_out%x = -vec_in%y - vec_in%z
+            vec_out%y =  vec_in%x + a*vec_in%y
+            vec_out%z =  vec_in%x*self%X%z + vec_in%z*(self%X%x - c)
+
+         end select
+      end select
+
+      return
+    end subroutine lin_roessler
+
+    subroutine adj_lin_roessler(self, vec_in, vec_out)
+      class(jacobian),            intent(in)  :: self
+      class(abstract_vector_rdp), intent(in)  :: vec_in
+      class(abstract_vector_rdp), intent(out) :: vec_out
+
+      select type(vec_in)
+      type is(state_vector)
+         select type(vec_out)
+         type is(state_vector)
+
+            vec_out%x =  vec_in%y + vec_in%z*self%X%z
+            vec_out%y = -vec_in%x + a*vec_in%y
+            vec_out%z = -vec_in%x + vec_in%z*(self%X%x - c)
+
+         end select
+      end select
+
+      return
+    end subroutine adj_lin_roessler
 
 end module LightKrylov_TestUtils