Merge pull request #95 from nekStab/update_eigs

Added default eigenvalue selection routine for `krylov_schur` in `eigs`
nekStab · Jun 14, 2024 · 86ec306 · 86ec306
2 parents 6df6313 + 2f4b248
commit 86ec306
Show file tree

Hide file tree

Showing 15 changed files with 169 additions and 221 deletions.
diff --git a/src/AbstractLinops.f90 b/src/AbstractLinops.f90
@@ -3,8 +3,9 @@ module lightkrylov_AbstractLinops
     use lightkrylov_utils
     use lightkrylov_AbstractVectors
     implicit none
+    private
 
-    character*128, parameter, private :: this_module = 'Lightkrylov_AbstractLinops'
+    character*128, parameter :: this_module = 'Lightkrylov_AbstractLinops'
 
     type, abstract, public :: abstract_linop
     contains
@@ -13,10 +14,6 @@ module lightkrylov_AbstractLinops
         generic, public :: assignment(=) => copy
     end type abstract_linop
 
-
-
-
-
     !------------------------------------------------------------------------------
     !-----     Definition of an abstract real(sp) operator with kind=sp     -----
     !------------------------------------------------------------------------------

diff --git a/src/AbstractLinops.fypp b/src/AbstractLinops.fypp
@@ -5,8 +5,9 @@ module lightkrylov_AbstractLinops
     use lightkrylov_utils
     use lightkrylov_AbstractVectors
     implicit none
+    private
 
-    character*128, parameter, private :: this_module = 'Lightkrylov_AbstractLinops'
+    character*128, parameter :: this_module = 'Lightkrylov_AbstractLinops'
 
     type, abstract, public :: abstract_linop
     contains
@@ -15,10 +16,6 @@ module lightkrylov_AbstractLinops
         generic, public :: assignment(=) => copy
     end type abstract_linop
 
-
-
-
-
     #:for kind, type in RC_KINDS_TYPES
     !------------------------------------------------------------------------------
     !-----     Definition of an abstract ${type}$ operator with kind=${kind}$     -----

diff --git a/src/AbstractVectors.f90 b/src/AbstractVectors.f90
@@ -2,8 +2,9 @@ module lightkrylov_AbstractVectors
     use lightkrylov_constants
     use lightkrylov_utils
     implicit none
+    private
 
-    character*128, parameter, private :: this_module = 'Lightkrylov_AbstractVectors'
+    character*128, parameter :: this_module = 'Lightkrylov_AbstractVectors'
 
     public :: innerprod
     public :: linear_combination

diff --git a/src/AbstractVectors.fypp b/src/AbstractVectors.fypp
@@ -4,8 +4,9 @@ module lightkrylov_AbstractVectors
     use lightkrylov_constants
     use lightkrylov_utils
     implicit none
+    private
 
-    character*128, parameter, private :: this_module = 'Lightkrylov_AbstractVectors'
+    character*128, parameter :: this_module = 'Lightkrylov_AbstractVectors'
 
     public :: innerprod
     public :: linear_combination

diff --git a/src/BaseKrylov.f90 b/src/BaseKrylov.f90
@@ -8,8 +8,9 @@ module lightkrylov_BaseKrylov
     use LightKrylov_AbstractVectors
     use LightKrylov_AbstractLinops
     implicit none
+    private
 
-    character*128, parameter, private :: this_module = 'LightKrylov_BaseKrylov'
+    character*128, parameter :: this_module = 'LightKrylov_BaseKrylov'
 
     public :: qr
     public :: apply_permutation_matrix
@@ -20,6 +21,9 @@ module lightkrylov_BaseKrylov
     public :: lanczos_bidiagonalization
     public :: lanczos_tridiagonalization
     public :: krylov_schur
+    public :: double_gram_schmidt_step
+    public :: eigvals_select_sp
+    public :: eigvals_select_dp
 
     interface qr
         module procedure qr_no_pivoting_rsp
@@ -118,6 +122,23 @@ module lightkrylov_BaseKrylov
         module procedure krylov_schur_cdp
     end interface
 
+    !----------------------------------------------------------
+    !-----     ABSTRACT EIGENVALUE SELECTOR INTERFACE     -----
+    !----------------------------------------------------------
+
+    abstract interface
+        function eigvals_select_sp(lambda) result(selected)
+            import sp
+            complex(sp), intent(in) :: lambda(:)
+            logical                       :: selected(size(lambda))
+        end function eigvals_select_sp
+        function eigvals_select_dp(lambda) result(selected)
+            import dp
+            complex(dp), intent(in) :: lambda(:)
+            logical                       :: selected(size(lambda))
+        end function eigvals_select_dp
+    end interface
+
 contains
 
     !-------------------------------------
@@ -135,9 +156,7 @@ subroutine initialize_krylov_subspace_rsp(X, X0)
         integer :: i, p, info
 
         ! Zero-out X.
-        do i = 1, size(X)
-            call X(i)%zero()
-        enddo
+        call zero_basis(X)
 
         ! Deals with optional args.
         if(present(X0)) then
@@ -151,8 +170,8 @@ subroutine initialize_krylov_subspace_rsp(X, X0)
             enddo
 
             ! Orthonormalize.
-            allocate(R(1:p, 1:p)) ; R = zero_rsp
-            allocate(perm(1:p)) ; perm = 0
+            allocate(R(p, p)) ; R = zero_rsp
+            allocate(perm(p)) ; perm = 0
             call qr(Q, R, perm, info)
             do i = 1, p
                 call X(i)%add(Q(i))
@@ -377,9 +396,7 @@ subroutine initialize_krylov_subspace_rdp(X, X0)
         integer :: i, p, info
 
         ! Zero-out X.
-        do i = 1, size(X)
-            call X(i)%zero()
-        enddo
+        call zero_basis(X)
 
         ! Deals with optional args.
         if(present(X0)) then
@@ -393,8 +410,8 @@ subroutine initialize_krylov_subspace_rdp(X, X0)
             enddo
 
             ! Orthonormalize.
-            allocate(R(1:p, 1:p)) ; R = zero_rdp
-            allocate(perm(1:p)) ; perm = 0
+            allocate(R(p, p)) ; R = zero_rdp
+            allocate(perm(p)) ; perm = 0
             call qr(Q, R, perm, info)
             do i = 1, p
                 call X(i)%add(Q(i))
@@ -619,9 +636,7 @@ subroutine initialize_krylov_subspace_csp(X, X0)
         integer :: i, p, info
 
         ! Zero-out X.
-        do i = 1, size(X)
-            call X(i)%zero()
-        enddo
+        call zero_basis(X)
 
         ! Deals with optional args.
         if(present(X0)) then
@@ -635,8 +650,8 @@ subroutine initialize_krylov_subspace_csp(X, X0)
             enddo
 
             ! Orthonormalize.
-            allocate(R(1:p, 1:p)) ; R = zero_rsp
-            allocate(perm(1:p)) ; perm = 0
+            allocate(R(p, p)) ; R = zero_rsp
+            allocate(perm(p)) ; perm = 0
             call qr(Q, R, perm, info)
             do i = 1, p
                 call X(i)%add(Q(i))
@@ -861,9 +876,7 @@ subroutine initialize_krylov_subspace_cdp(X, X0)
         integer :: i, p, info
 
         ! Zero-out X.
-        do i = 1, size(X)
-            call X(i)%zero()
-        enddo
+        call zero_basis(X)
 
         ! Deals with optional args.
         if(present(X0)) then
@@ -877,8 +890,8 @@ subroutine initialize_krylov_subspace_cdp(X, X0)
             enddo
 
             ! Orthonormalize.
-            allocate(R(1:p, 1:p)) ; R = zero_rdp
-            allocate(perm(1:p)) ; perm = 0
+            allocate(R(p, p)) ; R = zero_rdp
+            allocate(perm(p)) ; perm = 0
             call qr(Q, R, perm, info)
             do i = 1, p
                 call X(i)%add(Q(i))
@@ -3119,14 +3132,7 @@ subroutine krylov_schur_rsp(n, X, H, select_eigs)
         !! Krylov basis.
         real(sp), intent(inout) :: H(:, :)
         !! Upper Hessenberg matrix.
-        interface
-            function selector(lambda) result(out)
-                import sp
-                complex(sp), intent(in) :: lambda(:)
-                logical                       :: out(size(lambda))
-            end function selector
-        end interface
-        procedure(selector) :: select_eigs
+        procedure(eigvals_select_sp) :: select_eigs
         !! Procedure to select the eigenvalues to move in the upper left-block.
 
         !--------------------------------------
@@ -3170,9 +3176,7 @@ end function selector
 
         call X(n+1)%axpby(zero_rsp, X(kdim+1), one_rsp)
 
-        do i = n+2, size(X)
-            call X(i)%zero()
-        enddo
+        call zero_basis(X(n+2:))
 
         ! Update the Hessenberg matrix.
         b = matmul(H(kdim+1, :), Z)
@@ -3192,14 +3196,7 @@ subroutine krylov_schur_rdp(n, X, H, select_eigs)
         !! Krylov basis.
         real(dp), intent(inout) :: H(:, :)
         !! Upper Hessenberg matrix.
-        interface
-            function selector(lambda) result(out)
-                import dp
-                complex(dp), intent(in) :: lambda(:)
-                logical                       :: out(size(lambda))
-            end function selector
-        end interface
-        procedure(selector) :: select_eigs
+        procedure(eigvals_select_dp) :: select_eigs
         !! Procedure to select the eigenvalues to move in the upper left-block.
 
         !--------------------------------------
@@ -3243,9 +3240,7 @@ end function selector
 
         call X(n+1)%axpby(zero_rdp, X(kdim+1), one_rdp)
 
-        do i = n+2, size(X)
-            call X(i)%zero()
-        enddo
+        call zero_basis(X(n+2:))
 
         ! Update the Hessenberg matrix.
         b = matmul(H(kdim+1, :), Z)
@@ -3265,14 +3260,7 @@ subroutine krylov_schur_csp(n, X, H, select_eigs)
         !! Krylov basis.
         complex(sp), intent(inout) :: H(:, :)
         !! Upper Hessenberg matrix.
-        interface
-            function selector(lambda) result(out)
-                import sp
-                complex(sp), intent(in) :: lambda(:)
-                logical                       :: out(size(lambda))
-            end function selector
-        end interface
-        procedure(selector) :: select_eigs
+        procedure(eigvals_select_sp) :: select_eigs
         !! Procedure to select the eigenvalues to move in the upper left-block.
 
         !--------------------------------------
@@ -3316,9 +3304,7 @@ end function selector
 
         call X(n+1)%axpby(zero_csp, X(kdim+1), one_csp)
 
-        do i = n+2, size(X)
-            call X(i)%zero()
-        enddo
+        call zero_basis(X(n+2:))
 
         ! Update the Hessenberg matrix.
         b = matmul(H(kdim+1, :), Z)
@@ -3338,14 +3324,7 @@ subroutine krylov_schur_cdp(n, X, H, select_eigs)
         !! Krylov basis.
         complex(dp), intent(inout) :: H(:, :)
         !! Upper Hessenberg matrix.
-        interface
-            function selector(lambda) result(out)
-                import dp
-                complex(dp), intent(in) :: lambda(:)
-                logical                       :: out(size(lambda))
-            end function selector
-        end interface
-        procedure(selector) :: select_eigs
+        procedure(eigvals_select_dp) :: select_eigs
         !! Procedure to select the eigenvalues to move in the upper left-block.
 
         !--------------------------------------
@@ -3389,9 +3368,7 @@ end function selector
 
         call X(n+1)%axpby(zero_cdp, X(kdim+1), one_cdp)
 
-        do i = n+2, size(X)
-            call X(i)%zero()
-        enddo
+        call zero_basis(X(n+2:))
 
         ! Update the Hessenberg matrix.
         b = matmul(H(kdim+1, :), Z)

diff --git a/src/BaseKrylov.fypp b/src/BaseKrylov.fypp
@@ -10,8 +10,9 @@ module lightkrylov_BaseKrylov
     use LightKrylov_AbstractVectors
     use LightKrylov_AbstractLinops
     implicit none
+    private
 
-    character*128, parameter, private :: this_module = 'LightKrylov_BaseKrylov'
+    character*128, parameter :: this_module = 'LightKrylov_BaseKrylov'
 
     public :: qr
     public :: apply_permutation_matrix
@@ -22,6 +23,10 @@ module lightkrylov_BaseKrylov
     public :: lanczos_bidiagonalization
     public :: lanczos_tridiagonalization
     public :: krylov_schur
+    public :: double_gram_schmidt_step
+    #:for kind in REAL_KINDS
+    public :: eigvals_select_${kind}$
+    #:endfor
 
     interface qr
         #:for kind, type in RC_KINDS_TYPES
@@ -94,6 +99,20 @@ module lightkrylov_BaseKrylov
         #:endfor
     end interface
 
+    !----------------------------------------------------------
+    !-----     ABSTRACT EIGENVALUE SELECTOR INTERFACE     -----
+    !----------------------------------------------------------
+
+    abstract interface
+        #:for kind in REAL_KINDS
+        function eigvals_select_${kind}$(lambda) result(selected)
+            import ${kind}$
+            complex(${kind}$), intent(in) :: lambda(:)
+            logical                       :: selected(size(lambda))
+        end function eigvals_select_${kind}$
+        #:endfor
+    end interface
+
 contains
 
     !-------------------------------------
@@ -112,9 +131,7 @@ contains
         integer :: i, p, info
 
         ! Zero-out X.
-        do i = 1, size(X)
-            call X(i)%zero()
-        enddo
+        call zero_basis(X)
 
         ! Deals with optional args.
         if(present(X0)) then
@@ -128,8 +145,8 @@ contains
             enddo
 
             ! Orthonormalize.
-            allocate(R(1:p, 1:p)) ; R = zero_r${kind}$
-            allocate(perm(1:p)) ; perm = 0
+            allocate(R(p, p)) ; R = zero_r${kind}$
+            allocate(perm(p)) ; perm = 0
             call qr(Q, R, perm, info)
             do i = 1, p
                 call X(i)%add(Q(i))
@@ -889,14 +906,7 @@ contains
         !! Krylov basis.
         ${type}$, intent(inout) :: H(:, :)
         !! Upper Hessenberg matrix.
-        interface
-            function selector(lambda) result(out)
-                import ${kind}$
-                complex(${kind}$), intent(in) :: lambda(:)
-                logical                       :: out(size(lambda))
-            end function selector
-        end interface
-        procedure(selector) :: select_eigs
+        procedure(eigvals_select_${kind}$) :: select_eigs
         !! Procedure to select the eigenvalues to move in the upper left-block.
 
         !--------------------------------------
@@ -940,9 +950,7 @@ contains
 
         call X(n+1)%axpby(zero_${type[0]}$${kind}$, X(kdim+1), one_${type[0]}$${kind}$)
 
-        do i = n+2, size(X)
-            call X(i)%zero()
-        enddo
+        call zero_basis(X(n+2:))
 
         ! Update the Hessenberg matrix.
         b = matmul(H(kdim+1, :), Z)