Merge pull request #8 from nekStab/DLRA

Rank-adaptive Projector-Splitting Integrator for Dynamic Low-Rank Approximation of Lyapunov and Riccati equations

.gitignore

@@ -35,8 +35,13 @@
 doc
 
 # auxiliary files
-example/DLRA_laplacian2D_lti_lyapunov/output.txt
-example/DLRA_laplacian2D_lti_lyapunov/laplacian.gpp
+example/DLRA_laplacian2D_lti_lyapunov/*.txt
+example/DLRA_laplacian2D_lti_lyapunov/*.gpp
+example/DLRA_ginzburg_landau/*.txt
+example/DLRA_ginzburg_landau/*.gpp
+
+# local data
+local/
 
 # data files
 *.npy

example/DLRA_ginzburg_landau/ginzburg_landau_RK_lyapunov.f90

@@ -0,0 +1,268 @@
+module Ginzburg_Landau_RK_Lyapunov
+   ! Standard Library.
+   use stdlib_optval, only : optval
+   ! RKLIB module for time integration.
+   use rklib_module
+   ! LightKrylov for linear algebra.
+   use LightKrylov
+   use LightKrylov, only : wp => dp
+   ! Ginzburg Landau
+   use Ginzburg_Landau_Base
+   use Ginzburg_Landau_Operators
+   implicit none
+
+   private :: this_module
+   character*128, parameter :: this_module = 'Ginzburg_Landau_RK_Lyapunov'
+
+   public  :: GL_mat
+
+   !-------------------------------------------
+   !-----     LIGHTKRYLOV VECTOR TYPE     -----
+   !-------------------------------------------
+
+   type, extends(abstract_vector_rdp), public :: state_matrix
+      real(wp) :: state(N**2) = 0.0_wp
+   contains
+      private
+      procedure, pass(self), public :: zero => matrix_zero
+      procedure, pass(self), public :: dot => matrix_dot
+      procedure, pass(self), public :: scal => matrix_scal
+      procedure, pass(self), public :: axpby => matrix_axpby
+      procedure, pass(self), public :: rand => matrix_rand
+      procedure, pass(self), public :: get_size => matrix_get_size
+   end type state_matrix
+
+   !-------------------------------
+   !-----     RK LYAPUNOV     -----
+   !-------------------------------
+
+   type, extends(abstract_linop_rdp), public :: RK_lyapunov
+      real(wp) :: tau ! Integration time.
+   contains
+      private
+      procedure, pass(self), public :: matvec => direct_solver_lyap
+      procedure, pass(self), public :: rmatvec => adjoint_solver_lyap
+   end type RK_lyapunov
+
+contains
+
+   !-----     TYPE-BOUND PROCEDURE FOR MATRICES     -----
+
+   subroutine matrix_zero(self)
+      class(state_matrix), intent(inout) :: self
+      self%state = 0.0_wp
+      return
+   end subroutine matrix_zero
+
+   real(wp) function matrix_dot(self, vec) result(alpha)
+      class(state_matrix),        intent(in) :: self
+      class(abstract_vector_rdp), intent(in) :: vec
+      select type(vec)
+      type is(state_matrix)
+         alpha = dot_product(self%state, weight_mat*vec%state)
+      end select
+      return
+   end function matrix_dot
+
+   integer function matrix_get_size(self) result(N)
+     class(state_matrix), intent(in) :: self
+     N = N**2
+     return
+   end function matrix_get_size
+
+   subroutine matrix_scal(self, alpha)
+      class(state_matrix), intent(inout) :: self
+      real(wp),            intent(in)    :: alpha
+      self%state = self%state * alpha
+      return
+   end subroutine matrix_scal  
+
+   subroutine matrix_axpby(self, alpha, vec, beta)
+      class(state_matrix),        intent(inout) :: self
+      class(abstract_vector_rdp), intent(in)    :: vec
+      real(wp)         , intent(in)    :: alpha, beta
+      select type(vec)
+      type is(state_matrix)
+         self%state = alpha*self%state + beta*vec%state
+      end select
+      return
+   end subroutine matrix_axpby
+
+   subroutine matrix_rand(self, ifnorm)
+      class(state_matrix), intent(inout) :: self
+      logical, optional,   intent(in)    :: ifnorm
+      ! internals
+      logical :: normalize
+      real(wp) :: alpha
+      normalize = optval(ifnorm, .true.)
+      call random_number(self%state)
+      if (normalize) then
+         alpha = self%norm()
+         call self%scal(1.0/alpha)
+      endif
+      return
+   end subroutine matrix_rand
+
+   subroutine GL_mat(flat_mat_out, flat_mat_in, adjoint, transpose)
+
+      !> State vector.
+      real(wp), dimension(:), intent(in)  :: flat_mat_in
+      !> Time-derivative.
+      real(wp), dimension(:), intent(out) :: flat_mat_out
+      !> Adjoint
+      logical, optional :: adjoint
+      logical           :: adj
+      logical, optional :: transpose
+      logical           :: trans
+
+      !> Internal variables.
+      integer :: j
+      real(wp), dimension(N,N) :: mat, dmat
+
+      !> Deal with optional argument
+      adj   = optval(adjoint,.false.)
+      trans = optval(transpose,.false.)
+
+      !> Sets the internal variables.
+      mat  = reshape(flat_mat_in(1:N**2),(/N, N/))
+      dmat = 0.0_wp
+
+      if (adj) then
+         if (trans) then
+            do j = 1,N
+               call adjoint_GL(mat(:,j), dmat(j,:))
+            end do
+         else
+            do j = 1,N
+               call adjoint_GL(mat(:,j), dmat(:,j))
+            end do
+         end if
+      else
+         if (trans) then
+            do j = 1,N
+               call direct_GL(mat(:,j), dmat(j,:))
+            end do
+         else
+            do j = 1,N
+               call direct_GL(mat(:,j), dmat(:,j))
+            end do
+         end if
+      endif
+
+      !> Reshape for output
+      flat_mat_out = reshape(dmat, shape(flat_mat_in))
+
+      return
+   end subroutine GL_mat
+
+   !--------------------------------------
+   !-----     WRAPPERS FOR RKLIB     -----
+   !--------------------------------------
+
+   subroutine rhs_lyap(me, t, x_flat, f_flat)
+      ! Time-integrator.
+      class(rk_class), intent(inout)             :: me
+      ! Current time.
+      real(wp),        intent(in)                :: t
+      ! State vector.
+      real(wp),        dimension(:), intent(in)  :: x_flat
+      ! Time-derivative.
+      real(wp),        dimension(:), intent(out) :: f_flat
+
+      ! internals
+      real(wp),        dimension(N**2) :: x_tmp, AX_flat, XAH_flat
+
+      f_flat = 0.0_wp; AX_flat = 0.0_wp; XAH_flat = 0.0_wp; x_tmp = 0.0_wp
+      ! A @ X
+      call GL_mat( AX_flat, x_flat, adjoint = .false., transpose = .false.)
+      ! build X.T
+      x_tmp    = reshape(transpose(reshape(x_flat,   (/ N,N /))), shape(x_flat))
+      ! build ( A @ X.T ).T = X @ A.T
+      call GL_mat(XAH_flat, x_tmp,  adjoint = .false., transpose = .true.)
+      ! construct Lyapunov equation
+      f_flat = AX_flat + XAH_flat + BBTW_flat
+
+      return
+   end subroutine rhs_lyap
+
+   subroutine adjoint_rhs_lyap(me, t, x_flat, f_flat)
+      ! Time-integrator.
+      class(rk_class), intent(inout)             :: me
+      ! Current time.
+      real(wp),        intent(in)                :: t
+      ! State vector.
+      real(wp),        dimension(:), intent(in)  :: x_flat
+      ! Time-derivative.
+      real(wp),        dimension(:), intent(out) :: f_flat
+
+      ! internals
+      real(wp),        dimension(N**2) :: x_tmp, AHX_flat, XA_flat
+
+      f_flat = 0.0_wp; AHX_flat = 0.0_wp; XA_flat = 0.0_wp; x_tmp = 0.0_wp
+      ! A.T @ X
+      call GL_mat(AHX_flat, x_flat, adjoint = .true., transpose = .false.)
+      ! build X.T
+      x_tmp   = reshape(transpose(reshape(x_flat,  (/ N,N /))), shape(x_flat))
+      ! build ( A.T @ X.T ).T = X @ A
+      call GL_mat( XA_flat, x_tmp,  adjoint = .true., transpose = .true.)
+      ! construct Lyapunov equation
+      f_flat = AHX_flat + XA_flat + CTCW_flat
+
+      return
+   end subroutine adjoint_rhs_lyap
+
+   !------------------------------------------------------------------------
+   !-----     TYPE-BOUND PROCEDURES FOR THE EXPONENTIAL PROPAGATOR     -----
+   !------------------------------------------------------------------------
+
+   subroutine direct_solver_lyap(self, vec_in, vec_out)
+      ! Linear Operator.
+      class(rk_lyapunov),          intent(in)  :: self
+      ! Input vector.
+      class(abstract_vector_rdp),  intent(in)  :: vec_in
+      ! Output vector.
+      class(abstract_vector_rdp),  intent(out) :: vec_out
+
+      ! Time-integrator.
+      type(rks54_class) :: prop
+      real(kind=wp)     :: dt = 1.0_wp
+
+      select type(vec_in)
+      type is(state_matrix)
+         select type(vec_out)
+         type is(state_matrix)
+            ! Initialize propagator.
+            call prop%initialize(n=N**2, f=rhs_lyap)
+            ! Integrate forward in time.
+            call prop%integrate(0.0_wp, vec_in%state, dt, self%tau, vec_out%state)
+         end select
+      end select
+      return
+   end subroutine direct_solver_lyap
+
+   subroutine adjoint_solver_lyap(self, vec_in, vec_out)
+      ! Linear Operator.
+      class(rk_lyapunov),          intent(in)  :: self
+      ! Input vector.
+      class(abstract_vector_rdp),  intent(in)  :: vec_in
+      ! Output vector.
+      class(abstract_vector_rdp),  intent(out) :: vec_out
+
+      ! Time-integrator.
+      type(rks54_class) :: prop
+      real(kind=wp)     :: dt = 1.0_wp
+
+      select type(vec_in)
+      type is(state_matrix)
+         select type(vec_out)
+         type is(state_matrix)
+            ! Initialize propagator.
+            call prop%initialize(n=N**2, f=adjoint_rhs_lyap)
+            ! Integrate forward in time.
+            call prop%integrate(0.0_wp, vec_in%state, dt, self%tau, vec_out%state)
+         end select
+      end select
+      return
+   end subroutine adjoint_solver_lyap
+
+end module Ginzburg_Landau_RK_Lyapunov