New Github Action for CI

.github/workflows/ci.yml

@@ -0,0 +1,49 @@
+name: ci
+
+on: [push, pull_request]
+
+jobs:
+  test:
+    runs-on: ${{ matrix.os }}
+    strategy:
+      fail-fast: false
+      matrix:
+        include:
+          - os: ubuntu-latest
+            toolchain: {compiler: gcc, version: 12, flags: ['-cpp -O3 -march=native -mtune=native']}
+          - os: ubuntu-latest
+            toolchain: {compiler: gcc, version: 13, flags: ['-cpp -O3 -march=native -mtune=native']}
+          - os: ubuntu-latest
+            toolchain: {compiler: intel-classic, version: '2021.10', flags: ['-fpp -O0 -g']}
+          - os: ubuntu-latest
+            toolchain: {compiler: intel, version: '2024.1', flags: ['-fpp -O0 -g']}
+          # - os: ubuntu-latest
+          #   toolchain: {compiler: nvidia-hpc, version: '24.1', flags: ['-Mpreprocess -Ofast']}
+          - os: windows-latest
+            toolchain: {compiler: gcc, version: 13, flags: ['-cpp -O3 -march=native -mtune=native']}
+          # - os: windows-latest
+          #   toolchain: {compiler: intel-classic, version: '2021.10', flags: ["/fpp /O3"]}
+          # - os: windows-latest
+          #   toolchain: {compiler: intel, version: '2024.1', flags: ["/fpp /O3"]}
+          - os: macos-12
+            toolchain: {compiler: gcc, version: 13, flags: ['-cpp -O3 -march=native -mtune=native']}
+          # - os: macos-12
+          #   toolchain: {compiler: intel-classic, version: '2021.10', flags: ['-fpp -O3 -xhost']}
+
+    steps:
+      - name: Checkout code
+        uses: actions/checkout@v1
+
+      - uses: fortran-lang/setup-fortran@main
+        id: setup-fortran
+        with:
+          compiler: ${{ matrix.toolchain.compiler }}
+          version: ${{ matrix.toolchain.version }}
+
+      - name: Setup Fortran Package Manager
+        uses: fortran-lang/setup-fpm@v5
+        with:
+          github-token: ${{ secrets.GITHUB_TOKEN }}
+
+      - run: |
+          fpm test --flag "${{ join(matrix.toolchain.flags, ' ') }}"

fpm.toml

@@ -11,8 +11,6 @@ auto-executables = true
 auto-tests = true
 auto-examples = true
 module-naming = false
-link = ["blas", "lapack"]
-#link = ["openblas"]
 
 [install]
 library = true

src/AbstractLinops.f90

@@ -8,10 +8,6 @@ module lightkrylov_AbstractLinops
     character*128, parameter :: this_module = 'Lightkrylov_AbstractLinops'
 
     type, abstract, public :: abstract_linop
-    contains
-    private
-        procedure, pass(from), public :: copy
-        generic, public :: assignment(=) => copy
     end type abstract_linop
 
     !------------------------------------------------------------------------------
@@ -390,14 +386,6 @@ end subroutine abstract_matvec_cdp
 
 contains
 
-    subroutine copy(out, from)
-        class(abstract_linop), intent(in) :: from
-        class(abstract_linop), allocatable, intent(out) :: out
-        if (allocated(out)) deallocate(out)
-        allocate(out, source=from)
-        return
-    end subroutine copy
-
     subroutine id_matvec_rsp(self, vec_in, vec_out)
         class(Id_rsp), intent(in) :: self
         class(abstract_vector_rsp), intent(in) :: vec_in

src/AbstractLinops.fypp

@@ -10,10 +10,6 @@ module lightkrylov_AbstractLinops
     character*128, parameter :: this_module = 'Lightkrylov_AbstractLinops'
 
     type, abstract, public :: abstract_linop
-    contains
-    private
-        procedure, pass(from), public :: copy
-        generic, public :: assignment(=) => copy
     end type abstract_linop
 
     #:for kind, type in RC_KINDS_TYPES
@@ -153,14 +149,6 @@ module lightkrylov_AbstractLinops
     #:endfor
 contains
 
-    subroutine copy(out, from)
-        class(abstract_linop), intent(in) :: from
-        class(abstract_linop), allocatable, intent(out) :: out
-        if (allocated(out)) deallocate(out)
-        allocate(out, source=from)
-        return
-    end subroutine copy
-
     #:for kind, type in RC_KINDS_TYPES
     subroutine id_matvec_${type[0]}$${kind}$(self, vec_in, vec_out)
         class(Id_${type[0]}$${kind}$), intent(in) :: self

src/AbstractVectors.f90

@@ -57,10 +57,6 @@ module lightkrylov_AbstractVectors
 
 
     type, abstract, public :: abstract_vector
-    contains
-    private
-        procedure, pass(from), public :: copy => copy_vec
-        generic, public :: assignment(=) => copy
     end type abstract_vector
 
 
@@ -364,14 +360,6 @@ end function abstract_dot_cdp
 
 contains
 
-    subroutine copy_vec(out, from)
-        class(abstract_vector), intent(in)  :: from
-        class(abstract_vector), allocatable, intent(out) :: out
-        if (allocated(out)) deallocate(out)
-        allocate(out, source=from)
-        return
-    end subroutine copy_vec
-
     function norm_rsp(self) result(alpha)
         !! Compute the norm of an `abstract_vector`.
         class(abstract_vector_rsp), intent(in) :: self

src/AbstractVectors.fypp

@@ -48,10 +48,6 @@ module lightkrylov_AbstractVectors
 
 
     type, abstract, public :: abstract_vector
-    contains
-    private
-        procedure, pass(from), public :: copy => copy_vec
-        generic, public :: assignment(=) => copy
     end type abstract_vector
 
 
@@ -135,14 +131,6 @@ module lightkrylov_AbstractVectors
     #:endfor
 contains
 
-    subroutine copy_vec(out, from)
-        class(abstract_vector), intent(in)  :: from
-        class(abstract_vector), allocatable, intent(out) :: out
-        if (allocated(out)) deallocate(out)
-        allocate(out, source=from)
-        return
-    end subroutine copy_vec
-
     #:for kind, type in RC_KINDS_TYPES
     function norm_${type[0]}$${kind}$(self) result(alpha)
         !! Compute the norm of an `abstract_vector`.

src/Constants.f90

@@ -11,7 +11,7 @@ module lightkrylov_constants
 
     integer , parameter, public :: dp = selected_real_kind(15, 307)
     !! Definition of the double precision data type.
-    real(dp), parameter, public :: atol_dp = 10.0_sp ** -precision(1.0_dp)
+    real(dp), parameter, public :: atol_dp = 10.0_dp ** -precision(1.0_dp)
     !! Definition of the absolute tolerance for double precision computations.
     real(dp), parameter, public :: rtol_dp = sqrt(atol_dp)
     !! Definition of the relative tolerance for double precision computations.

src/IterativeSolvers.f90

@@ -5,7 +5,7 @@ module lightkrylov_IterativeSolvers
     use stdlib_sorting, only: sort_index
     use stdlib_optval, only: optval
     use stdlib_io_npy, only: save_npy
-    use stdlib_linalg, only: lstsq
+    use stdlib_linalg, only: lstsq, svd
     use stdlib_stats, only: median
 
     use lightkrylov_constants
@@ -1230,17 +1230,19 @@ subroutine svds_rsp(A, U, S, V, residuals, info, kdim, tolerance, verbosity)
             call lanczos_bidiagonalization(A, Uwrk, Vwrk, B, info, kstart=k, kend=k, verbosity=verbosity, tol=tol)
             call check_info(info, 'lanczos_bidiagonalization', module=this_module, procedure='svds_rsp')
 
-            ! SVD of the k x k bidiagonal matrix.
+            ! SVD of the k x k bidiagonal matrix and residual computation.
             svdvals_wrk = 0.0_sp ; umat = 0.0_sp ; vmat = 0.0_sp
-            call svd(B(:k, :k), umat(:k, :k), svdvals_wrk(:k), vmat(:k, :k))
 
-            ! Compute residuals.
-            beta = B(k+1, k)
-            residuals_wrk(:k) = compute_residual_rsp(beta, vmat(k, :k))
+            if (k > 1) then
+                call svd(B(:k, :k), svdvals_wrk(:k), umat(:k, :k), vmat(:k, :k))
+                vmat(:k, :k) = transpose(vmat(:k, :k))
 
-            ! Check for convergence.
-            conv = count(residuals_wrk(:k) < tol)
-            if (conv >= nsv) exit lanczos
+                residuals_wrk(:k) = compute_residual_rsp(B(k+1, k), vmat(k, :k))
+
+                ! Check for convergence.
+                conv = count(residuals_wrk(:k) < tol)
+                if (conv >= nsv) exit lanczos
+            endif
         enddo lanczos
 
         !--------------------------------
@@ -1325,17 +1327,19 @@ subroutine svds_rdp(A, U, S, V, residuals, info, kdim, tolerance, verbosity)
             call lanczos_bidiagonalization(A, Uwrk, Vwrk, B, info, kstart=k, kend=k, verbosity=verbosity, tol=tol)
             call check_info(info, 'lanczos_bidiagonalization', module=this_module, procedure='svds_rdp')
 
-            ! SVD of the k x k bidiagonal matrix.
+            ! SVD of the k x k bidiagonal matrix and residual computation.
             svdvals_wrk = 0.0_dp ; umat = 0.0_dp ; vmat = 0.0_dp
-            call svd(B(:k, :k), umat(:k, :k), svdvals_wrk(:k), vmat(:k, :k))
 
-            ! Compute residuals.
-            beta = B(k+1, k)
-            residuals_wrk(:k) = compute_residual_rdp(beta, vmat(k, :k))
+            if (k > 1) then
+                call svd(B(:k, :k), svdvals_wrk(:k), umat(:k, :k), vmat(:k, :k))
+                vmat(:k, :k) = transpose(vmat(:k, :k))
 
-            ! Check for convergence.
-            conv = count(residuals_wrk(:k) < tol)
-            if (conv >= nsv) exit lanczos
+                residuals_wrk(:k) = compute_residual_rdp(B(k+1, k), vmat(k, :k))
+
+                ! Check for convergence.
+                conv = count(residuals_wrk(:k) < tol)
+                if (conv >= nsv) exit lanczos
+            endif
         enddo lanczos
 
         !--------------------------------
@@ -1420,17 +1424,19 @@ subroutine svds_csp(A, U, S, V, residuals, info, kdim, tolerance, verbosity)
             call lanczos_bidiagonalization(A, Uwrk, Vwrk, B, info, kstart=k, kend=k, verbosity=verbosity, tol=tol)
             call check_info(info, 'lanczos_bidiagonalization', module=this_module, procedure='svds_csp')
 
-            ! SVD of the k x k bidiagonal matrix.
+            ! SVD of the k x k bidiagonal matrix and residual computation.
             svdvals_wrk = 0.0_sp ; umat = 0.0_sp ; vmat = 0.0_sp
-            call svd(B(:k, :k), umat(:k, :k), svdvals_wrk(:k), vmat(:k, :k))
 
-            ! Compute residuals.
-            beta = B(k+1, k)
-            residuals_wrk(:k) = compute_residual_csp(beta, vmat(k, :k))
+            if (k > 1) then
+                call svd(B(:k, :k), svdvals_wrk(:k), umat(:k, :k), vmat(:k, :k))
+                vmat(:k, :k) = conjg(transpose(vmat(:k, :k)))
 
-            ! Check for convergence.
-            conv = count(residuals_wrk(:k) < tol)
-            if (conv >= nsv) exit lanczos
+                residuals_wrk(:k) = compute_residual_csp(B(k+1, k), vmat(k, :k))
+
+                ! Check for convergence.
+                conv = count(residuals_wrk(:k) < tol)
+                if (conv >= nsv) exit lanczos
+            endif
         enddo lanczos
 
         !--------------------------------
@@ -1515,17 +1521,19 @@ subroutine svds_cdp(A, U, S, V, residuals, info, kdim, tolerance, verbosity)
             call lanczos_bidiagonalization(A, Uwrk, Vwrk, B, info, kstart=k, kend=k, verbosity=verbosity, tol=tol)
             call check_info(info, 'lanczos_bidiagonalization', module=this_module, procedure='svds_cdp')
 
-            ! SVD of the k x k bidiagonal matrix.
+            ! SVD of the k x k bidiagonal matrix and residual computation.
             svdvals_wrk = 0.0_dp ; umat = 0.0_dp ; vmat = 0.0_dp
-            call svd(B(:k, :k), umat(:k, :k), svdvals_wrk(:k), vmat(:k, :k))
 
-            ! Compute residuals.
-            beta = B(k+1, k)
-            residuals_wrk(:k) = compute_residual_cdp(beta, vmat(k, :k))
+            if (k > 1) then
+                call svd(B(:k, :k), svdvals_wrk(:k), umat(:k, :k), vmat(:k, :k))
+                vmat(:k, :k) = conjg(transpose(vmat(:k, :k)))
 
-            ! Check for convergence.
-            conv = count(residuals_wrk(:k) < tol)
-            if (conv >= nsv) exit lanczos
+                residuals_wrk(:k) = compute_residual_cdp(B(k+1, k), vmat(k, :k))
+
+                ! Check for convergence.
+                conv = count(residuals_wrk(:k) < tol)
+                if (conv >= nsv) exit lanczos
+            endif
         enddo lanczos
 
         !--------------------------------

src/IterativeSolvers.fypp

@@ -7,7 +7,7 @@ module lightkrylov_IterativeSolvers
     use stdlib_sorting, only: sort_index
     use stdlib_optval, only: optval
     use stdlib_io_npy, only: save_npy
-    use stdlib_linalg, only: lstsq
+    use stdlib_linalg, only: lstsq, svd
     use stdlib_stats, only: median
 
     use lightkrylov_constants
@@ -466,17 +466,23 @@ contains
             call lanczos_bidiagonalization(A, Uwrk, Vwrk, B, info, kstart=k, kend=k, verbosity=verbosity, tol=tol)
             call check_info(info, 'lanczos_bidiagonalization', module=this_module, procedure='svds_${type[0]}$${kind}$')
 
-            ! SVD of the k x k bidiagonal matrix.
+            ! SVD of the k x k bidiagonal matrix and residual computation.
             svdvals_wrk = 0.0_${kind}$ ; umat = 0.0_${kind}$ ; vmat = 0.0_${kind}$
-            call svd(B(:k, :k), umat(:k, :k), svdvals_wrk(:k), vmat(:k, :k))
 
-            ! Compute residuals.
-            beta = B(k+1, k)
-            residuals_wrk(:k) = compute_residual_${type[0]}$${kind}$(beta, vmat(k, :k))
+            if (k > 1) then
+                call svd(B(:k, :k), svdvals_wrk(:k), umat(:k, :k), vmat(:k, :k))
+                #:if type[0] == "r"
+                vmat(:k, :k) = transpose(vmat(:k, :k))
+                #:else
+                vmat(:k, :k) = conjg(transpose(vmat(:k, :k)))
+                #:endif
 
-            ! Check for convergence.
-            conv = count(residuals_wrk(:k) < tol)
-            if (conv >= nsv) exit lanczos
+                residuals_wrk(:k) = compute_residual_${type[0]}$${kind}$(B(k+1, k), vmat(k, :k))
+
+                ! Check for convergence.
+                conv = count(residuals_wrk(:k) < tol)
+                if (conv >= nsv) exit lanczos
+            endif
         enddo lanczos
 
         !--------------------------------

src/LightKrylov.f90

@@ -28,10 +28,6 @@ module LightKrylov
     public :: cg_sp_opts
     public :: gmres_dp_opts
     public :: cg_dp_opts
-    public :: gmres_sp_opts
-    public :: cg_sp_opts
-    public :: gmres_dp_opts
-    public :: cg_dp_opts
 
     ! AbstractVectors exports.
     public :: abstract_vector

src/LightKrylov.fypp

@@ -26,7 +26,7 @@ module LightKrylov
     public :: dp, atol_dp, rtol_dp
 
     ! Utils exports.
-    #:for kind, type in RC_KINDS_TYPES
+    #:for kind in REAL_KINDS
     public :: gmres_${kind}$_opts
     public :: cg_${kind}$_opts
     #:endfor