From b0a74b111c6a50c703cfabb3861a799686451ce4 Mon Sep 17 00:00:00 2001
From: Jean-Christophe <loiseau.jc@gmail.com>
Date: Thu, 27 Mar 2025 09:19:24 +0100
Subject: [PATCH 1/7] Working implementation of expm.

---
 src/CMakeLists.txt                      |  5 +-
 src/stdlib_linalg.fypp                  | 14 ++++
 src/stdlib_linalg_matrix_functions.fypp | 68 +++++++++++++++++++
 test/linalg/CMakeLists.txt              |  2 +
 test/linalg/test_linalg_expm.fypp       | 90 +++++++++++++++++++++++++
 5 files changed, 177 insertions(+), 2 deletions(-)
 create mode 100644 src/stdlib_linalg_matrix_functions.fypp
 create mode 100644 test/linalg/test_linalg_expm.fypp

diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 933da34de..ed3b67d8c 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -32,16 +32,17 @@ set(fppFiles
     stdlib_linalg_kronecker.fypp
     stdlib_linalg_cross_product.fypp
     stdlib_linalg_eigenvalues.fypp
-    stdlib_linalg_solve.fypp    
+    stdlib_linalg_solve.fypp
     stdlib_linalg_determinant.fypp
     stdlib_linalg_qr.fypp
     stdlib_linalg_inverse.fypp
     stdlib_linalg_pinv.fypp
     stdlib_linalg_norms.fypp
     stdlib_linalg_state.fypp
-    stdlib_linalg_svd.fypp 
+    stdlib_linalg_svd.fypp
     stdlib_linalg_cholesky.fypp
     stdlib_linalg_schur.fypp
+    stdlib_linalg_matrix_functions.fypp
     stdlib_optval.fypp
     stdlib_selection.fypp
     stdlib_sorting.fypp
diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp
index c56f40bed..52e275c7e 100644
--- a/src/stdlib_linalg.fypp
+++ b/src/stdlib_linalg.fypp
@@ -28,6 +28,7 @@ module stdlib_linalg
   public :: eigh
   public :: eigvals
   public :: eigvalsh
+  public :: expm
   public :: eye
   public :: inv
   public :: invert
@@ -1678,6 +1679,19 @@ module stdlib_linalg
       #:endfor            
   end interface mnorm
 
+!> Matrix exponential: function interface
+  interface expm
+    #:for rk,rt,ri in RC_KINDS_TYPES
+    module function expm_${ri}$(A, order) result(E)
+        ${rt}$, intent(in) :: A(:, :)
+        !! On entry, the original matrix. On exit, its exponential.
+        integer(ilp), optional, intent(in) :: order
+        !! Order of the rational approximation.
+        ${rt}$, allocatable :: E(:, :)
+    end function expm_${ri}$
+    #:endfor
+  end interface expm
+
 contains
 
 
diff --git a/src/stdlib_linalg_matrix_functions.fypp b/src/stdlib_linalg_matrix_functions.fypp
new file mode 100644
index 000000000..ccd93d002
--- /dev/null
+++ b/src/stdlib_linalg_matrix_functions.fypp
@@ -0,0 +1,68 @@
+#:include "common.fypp"
+#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
+submodule (stdlib_linalg) stdlib_linalg_matrix_functions
+     use stdlib_linalg_constants
+     use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, &
+         LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR
+     implicit none
+
+    contains
+
+    #:for rk,rt,ri in RC_KINDS_TYPES 
+    module function expm_${ri}$(A, order) result(E)
+        ${rt}$, intent(in) :: A(:, :)
+        integer(ilp), optional, intent(in) :: order
+        ${rt}$, allocatable :: E(:, :)
+
+        ${rt}$, allocatable :: A2(:, :), Q(:, :), X(:, :)
+        real(${rk}$) :: a_norm, c
+        integer(ilp) :: n, ee, k, s
+        logical(lk) :: p
+        integer(ilp) :: p_order
+
+        ! Deal with optional args.
+        p_order = 10 ; if (present(order)) p_order = order
+
+        n = size(A, 1)
+
+        ! Compute the L-infinity norm.
+        a_norm = mnorm(A, "inf")
+
+        ! Determine scaling factor for the matrix.
+        ee = int(log(a_norm) / log(2.0_${rk}$)) + 1
+        s = max(0, ee+1)
+
+        ! Scale the input matrix & initialize polynomial.
+        A2 = A / 2.0_${rk}$**s
+        X = A2
+
+        ! Initialize P & Q and add first step.
+        c = 0.5_${rk}$
+        E = eye(n, mold=1.0_${rk}$) ; E = E + c*A2
+
+        Q = eye(n, mold=1.0_${rk}$) ; 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
+
+        E = matmul(inv(Q), E)
+        do k = 1, s
+            E = matmul(E, E)
+        enddo
+
+        return
+    end function
+    #:endfor
+
+end submodule stdlib_linalg_matrix_functions
diff --git a/test/linalg/CMakeLists.txt b/test/linalg/CMakeLists.txt
index c496bd2b3..a9e5431a3 100644
--- a/test/linalg/CMakeLists.txt
+++ b/test/linalg/CMakeLists.txt
@@ -16,6 +16,7 @@ set(
   "test_linalg_svd.fypp"
   "test_linalg_matrix_property_checks.fypp"
   "test_linalg_sparse.fypp"
+  "test_linalg_expm.fypp"
 )
 fypp_f90("${fyppFlags}" "${fppFiles}" outFiles)
 
@@ -35,3 +36,4 @@ ADDTEST(linalg_schur)
 ADDTEST(linalg_svd)
 ADDTEST(blas_lapack)
 ADDTEST(linalg_sparse)
+ADDTEST(linalg_expm)
diff --git a/test/linalg/test_linalg_expm.fypp b/test/linalg/test_linalg_expm.fypp
new file mode 100644
index 000000000..51d89a2ca
--- /dev/null
+++ b/test/linalg/test_linalg_expm.fypp
@@ -0,0 +1,90 @@
+#:include "common.fypp"
+#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
+! Test Schur decomposition
+module test_linalg_expm
+    use testdrive, only: error_type, check, new_unittest, unittest_type
+    use stdlib_linalg_constants
+    use stdlib_linalg, only: expm, eye, mnorm
+
+    implicit none (type,external)
+    
+    public :: test_expm_computation
+
+    contains
+
+    !> schur decomposition tests
+    subroutine test_expm_computation(tests)
+        !> Collection of tests
+        type(unittest_type), allocatable, intent(out) :: tests(:)
+        
+        allocate(tests(0))
+        
+        #:for rk,rt,ri in RC_KINDS_TYPES
+        tests = [tests, new_unittest("expm_${ri}$",test_expm_${ri}$)]
+        #:endfor
+
+    end subroutine test_expm_computation
+
+    !> Matrix exponential with analytic expression.
+    #:for rk,rt,ri in RC_KINDS_TYPES
+    subroutine test_expm_${ri}$(error)
+        type(error_type), allocatable, intent(out) :: error
+        ! Problem dimension.
+        integer(ilp), parameter :: n = 5, m = 6
+        ! Test matrix.
+        ${rt}$ :: A(n, n), E(n, n), Eref(n, n)
+        real(${rk}$) :: err
+        integer(ilp) :: i, j
+
+        ! Initialize matrix.
+        A = 0.0_${rk}$
+        do i = 1, n-1
+            A(i, i+1) = m*1.0_${rk}$
+        enddo
+
+        ! Reference with analytical exponential.
+        Eref = eye(n, mold=1.0_${rk}$)
+        do i = 1, n-1
+            do j = 1, n-i
+                Eref(i, i+j) = Eref(i, i+j-1)*m/j
+            enddo
+        enddo
+
+        ! Compute matrix exponential.
+        E = expm(A)
+
+        ! Check result.
+        err = mnorm(Eref - E, "inf")
+        call check(error, err < (n**2)*epsilon(1.0_${rk}$), "Analytical matrix exponential.")
+        if (allocated(error)) return
+        return        
+    end subroutine test_expm_${ri}$    
+    #:endfor    
+
+end module test_linalg_expm
+
+program test_expm
+    use, intrinsic :: iso_fortran_env, only : error_unit
+    use testdrive, only : run_testsuite, new_testsuite, testsuite_type
+    use test_linalg_expm, only : test_expm_computation
+    implicit none
+    integer :: stat, is
+    type(testsuite_type), allocatable :: testsuites(:)
+    character(len=*), parameter :: fmt = '("#", *(1x, a))'
+
+    stat = 0
+
+    testsuites = [ &
+        new_testsuite("linalg_expm", test_expm_computation) &
+        ]
+
+    do is = 1, size(testsuites)
+        write(error_unit, fmt) "Testing:", testsuites(is)%name
+        call run_testsuite(testsuites(is)%collect, error_unit, stat)
+    end do
+
+    if (stat > 0) then
+        write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!"
+        error stop
+    end if
+end program test_expm

From fa36f339254ac551561149597b3cbf741899ff30 Mon Sep 17 00:00:00 2001
From: Jean-Christophe <loiseau.jc@gmail.com>
Date: Thu, 27 Mar 2025 12:14:10 +0100
Subject: [PATCH 2/7] Improved implementation + error handling.

---
 src/stdlib_linalg_matrix_functions.fypp | 114 ++++++++++++++++++------
 1 file changed, 88 insertions(+), 26 deletions(-)

diff --git a/src/stdlib_linalg_matrix_functions.fypp b/src/stdlib_linalg_matrix_functions.fypp
index ccd93d002..7bd8b8fcb 100644
--- a/src/stdlib_linalg_matrix_functions.fypp
+++ b/src/stdlib_linalg_matrix_functions.fypp
@@ -1,68 +1,130 @@
 #:include "common.fypp"
 #:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
 submodule (stdlib_linalg) stdlib_linalg_matrix_functions
-     use stdlib_linalg_constants
-     use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, &
+    use stdlib_linalg_constants
+    use stdlib_linalg_lapack, only: gesv
+    use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, &
          LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR
-     implicit none
+    implicit none
 
-    contains
+    #:for rk, rt, ri in (REAL_KINDS_TYPES)
+    ${rt}$, parameter :: zero_${ri}$ = 0._${rk}$
+    ${rt}$, parameter :: one_${ri}$  = 1._${rk}$
+    #:endfor
+    #:for rk, rt, ri in (CMPLX_KINDS_TYPES)
+    ${rt}$, parameter :: zero_${ri}$ = (0._${rk}$, 0._${rk}$)
+    ${rt}$, parameter :: one_${ri}$  = (1._${rk}$, 0._${rk}$)
+    #:endfor
+
+contains
 
     #:for rk,rt,ri in RC_KINDS_TYPES 
-    module function expm_${ri}$(A, order) result(E)
+    module function stdlib_expm_${ri}$(A, order, err) result(E)
+        !> Input matrix A(n, n).
         ${rt}$, intent(in) :: A(:, :)
+        !> [optional] Order of the Pade approximation.
         integer(ilp), optional, intent(in) :: order
+        !> [optional] State return flag.
+        type(linalg_state_type), optional, intent(out) :: err
+        !> Exponential of the input matrix E = exp(A).
         ${rt}$, allocatable :: E(:, :)
 
+        ! Internal variables.
         ${rt}$, allocatable :: A2(:, :), Q(:, :), X(:, :)
-        real(${rk}$) :: a_norm, c
-        integer(ilp) :: n, ee, k, s
-        logical(lk) :: p
-        integer(ilp) :: p_order
+        real(${rk}$)        :: a_norm, c
+        integer(ilp)        :: m, n, ee, k, s, order_, i, j
+        logical(lk)         :: p
+        character(len=*), parameter :: this = "expm"
+        type(linalg_state_type)     :: err0
 
         ! Deal with optional args.
-        p_order = 10 ; if (present(order)) p_order = order
+        order_ = 10 ; if (present(order)) order_ = order
+
+        ! Problem's dimension.
+        m = size(A, 1) ; n = size(A, 2)
 
-        n = size(A, 1)
+        if (m /= n) then
+            err = linalg_state_type(this,LINALG_VALUE_ERROR,'Invalid matrix size A=',[m, n])
+            call linalg_error_handling(err0, err)
+        else if (order_ < 0) then
+            err = linalg_state_type(this, LINALG_VALUE_ERROR, 'Order of Pade approximation &
+                                    needs to be positive, order=', order_)
+            call linalg_error_handling(err0, err)
+        endif
 
         ! Compute the L-infinity norm.
         a_norm = mnorm(A, "inf")
 
         ! Determine scaling factor for the matrix.
         ee = int(log(a_norm) / log(2.0_${rk}$)) + 1
-        s = max(0, ee+1)
+        s  = max(0, ee+1)
 
         ! Scale the input matrix & initialize polynomial.
-        A2 = A / 2.0_${rk}$**s
-        X = A2
+        A2 = A/2.0_${rk}$**s ; X = A2
 
-        ! Initialize P & Q and add first step.
+        ! First step of the Pade approximation.
         c = 0.5_${rk}$
-        E = eye(n, mold=1.0_${rk}$) ; E = E + c*A2
-
-        Q = eye(n, mold=1.0_${rk}$) ; Q = Q - c*A2
+        allocate (E, source=A2) ; allocate (Q, source=A2)
+        do concurrent(i=1:n, j=1:n)
+            E(i, j) =  c*E(i, j) ; if (i == j) E(i, j) = 1.0_${rk}$ + E(i, j) ! E = I + c*A2
+            Q(i, j) = -c*Q(i, j) ; if (i == j) Q(i, j) = 1.0_${rk}$ + Q(i, j) ! Q = I - c*A2
+        enddo
 
         ! Iteratively compute the Pade approximation.
         p = .true.
-        do k = 2, p_order
-            c = c*(p_order - k + 1) / (k * (2*p_order - k + 1))
+        do k = 2, order_
+            c = c * (order_ - k + 1) / (k * (2*order_ - k + 1))
             X = matmul(A2, X)
-            E = E + c*X
+            do concurrent(i=1:n, j=1:n)
+                E(i, j) = E(i, j) + c*X(i, j)       ! E = E + c*X
+            enddo
             if (p) then
-                Q = Q + c*X
+                do concurrent(i=1:n, j=1:n)
+                    Q(i, j) = Q(i, j) + c*X(i, j)   ! Q = Q + c*X
+                enddo
             else
-                Q = Q - c*X
+                do concurrent(i=1:n, j=1:n)
+                    Q(i, j) = Q(i, j) - c*X(i, j)   ! Q = Q - c*X
+                enddo
             endif
             p = .not. p
         enddo
 
-        E = matmul(inv(Q), E)
+        block
+            integer(ilp) :: ipiv(n), info
+            call gesv(n, n, Q, n, ipiv, E, n, info) ! E = inv(Q) @ E
+            call handle_gesv_info(info, n, n, n, err0)
+            call linalg_error_handling(err0, err)
+        end block
+
+        ! This loop should eventually be replaced by a fast matrix_power function.
         do k = 1, s
             E = matmul(E, E)
         enddo
-
         return
-    end function
+    contains
+        elemental subroutine handle_gesv_info(info,lda,n,nrhs,err)
+            integer(ilp), intent(in) :: info,lda,n,nrhs
+            type(linalg_state_type), intent(out) :: err
+            ! Process output
+            select case (info)
+            case (0)
+            ! Success
+            case (-1)
+                err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid problem size n=',n)
+            case (-2)
+                err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid rhs size n=',nrhs)
+            case (-4)
+                err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size a=',[lda,n])
+            case (-7)
+                err = linalg_state_type(this,LINALG_ERROR,'invalid matrix size a=',[lda,n])
+            case (1:)
+                err = linalg_state_type(this,LINALG_ERROR,'singular matrix')
+            case default
+                err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'catastrophic error')
+            end select
+        end subroutine handle_gesv_info
+    end function stdlib_expm_${ri}$
     #:endfor
 
 end submodule stdlib_linalg_matrix_functions

From d8b1857966a038745b342ada4137917ceeb53508 Mon Sep 17 00:00:00 2001
From: Jean-Christophe <loiseau.jc@gmail.com>
Date: Thu, 27 Mar 2025 12:14:20 +0100
Subject: [PATCH 3/7] Added docstring for the interface.

---
 src/stdlib_linalg.fypp | 43 +++++++++++++++++++++++++++++++++++++-----
 1 file changed, 38 insertions(+), 5 deletions(-)

diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp
index 52e275c7e..1d752b223 100644
--- a/src/stdlib_linalg.fypp
+++ b/src/stdlib_linalg.fypp
@@ -1679,16 +1679,49 @@ module stdlib_linalg
       #:endfor            
   end interface mnorm
 
-!> Matrix exponential: function interface
+  !> Matrix exponential: function interface
   interface expm
+    !! version : experimental
+    !!
+    !! Computes the exponential of a matrix using a rational Pade approximation.
+    !!
+    !! ### Description
+    !!
+    !! This interface provides methods for computing the exponential of a matrix
+    !! represented as a standard Fortran rank-2 array. Supported data types include
+    !! `real` and `complex`.
+    !!
+    !! By default, the order of the Pade approximation is set to 10. It can be changed
+    !! via the `order` argument which must be non-negative.
+    !!
+    !! If the input matrix is non-square or the order of the Pade approximation is
+    !! negative, the function returns an error state.
+    !!
+    !! ### Example
+    !!
+    !! ```fortran
+    !!  real(dp) :: A(3, 3), E(3, 3)
+    !!
+    !!  A = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9], [3, 3])
+    !!
+    !!  ! Default Pade approximation of the matrix exponential.
+    !!  E = expm(A)
+    !!
+    !!  ! Pade approximation with specified order.
+    !!  E = expm(A, order=12)
+    !! ```
+    !!
     #:for rk,rt,ri in RC_KINDS_TYPES
-    module function expm_${ri}$(A, order) result(E)
+    module function stdlib_expm_${ri}$(A, order, err) result(E)
+        !> Input matrix a(n, n).
         ${rt}$, intent(in) :: A(:, :)
-        !! On entry, the original matrix. On exit, its exponential.
+        !> [optional] Order of the Pade approximation (default `order=10`)
         integer(ilp), optional, intent(in) :: order
-        !! Order of the rational approximation.
+        !> [optional] State return flag. On error, if not requested, the code will stop.
+        type(linalg_state_type), optional, intent(out) :: err
+        !> Exponential of the input matrix E = exp(A).
         ${rt}$, allocatable :: E(:, :)
-    end function expm_${ri}$
+    end function stdlib_expm_${ri}$
     #:endfor
   end interface expm
 

From 4089d183affda011e2b23ec8ac276b8d95d032c8 Mon Sep 17 00:00:00 2001
From: Jean-Christophe <loiseau.jc@gmail.com>
Date: Thu, 27 Mar 2025 14:23:01 +0100
Subject: [PATCH 4/7] Specs + example.

---
 doc/specs/stdlib_linalg.md      | 34 +++++++++++++++++++++++++++++++++
 example/linalg/CMakeLists.txt   |  1 +
 example/linalg/example_expm.f90 |  7 +++++++
 src/stdlib_linalg.fypp          |  1 +
 4 files changed, 43 insertions(+)
 create mode 100644 example/linalg/example_expm.f90

diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md
index 33fa43ed0..e7afdd912 100644
--- a/doc/specs/stdlib_linalg.md
+++ b/doc/specs/stdlib_linalg.md
@@ -1880,3 +1880,37 @@ If `err` is not present, exceptions trigger an `error stop`.
 {!example/linalg/example_mnorm.f90!}
 ```
 
+## `expm` - Computes the matrix exponential {#expm}
+
+### Status
+
+Experimental
+
+### Description
+
+Given a matrix \(A\), this function compute its matrix exponential \(E = \exp(A)\) using a Pade approximation.
+
+### Syntax
+
+`E = ` [[stdlib_linalg(module):expm(interface)]] `(a [, order, err])`
+
+### Arguments
+
+`a`: Shall be a rank-2 `real` or `complex` array containing the data. It is an `intent(in)` argument.
+
+`order` (optional): Shall be a non-negative `integer` value specifying the order of the Pade approximation. By default `order=10`. It is an `intent(in)` argument. 
+
+`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. 
+### Return value
+
+The returned array `E` contains the Pade approximation of \(\exp(A)\).
+
+If `A` is non-square or `order` is negative, it raise a `LINALG_VALUE_ERROR`.
+If `err` is not present, exceptions trigger an `error stop`.
+
+### Example
+
+```fortran
+{!example/linalg/example_expm.f90!}
+```
+
diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt
index 10f982a02..693fb0308 100644
--- a/example/linalg/CMakeLists.txt
+++ b/example/linalg/CMakeLists.txt
@@ -52,3 +52,4 @@ ADD_EXAMPLE(qr)
 ADD_EXAMPLE(qr_space)
 ADD_EXAMPLE(cholesky)
 ADD_EXAMPLE(chol)
+ADD_EXAMPLE(expm)
diff --git a/example/linalg/example_expm.f90 b/example/linalg/example_expm.f90
new file mode 100644
index 000000000..492b20323
--- /dev/null
+++ b/example/linalg/example_expm.f90
@@ -0,0 +1,7 @@
+program example_expm
+   use stdlib_linalg, only: expm
+   implicit none
+   real :: A(3, 3), E(3, 3)
+   A = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9], [3, 3])
+   E = expm(A)
+end program example_expm
diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp
index 1d752b223..567ae4033 100644
--- a/src/stdlib_linalg.fypp
+++ b/src/stdlib_linalg.fypp
@@ -1684,6 +1684,7 @@ module stdlib_linalg
     !! version : experimental
     !!
     !! Computes the exponential of a matrix using a rational Pade approximation.
+    !! ([Specification](../page/specs/stdlib_linalg.html#expm))
     !!
     !! ### Description
     !!

From c6857bcf5a812954afcf3cabb23c81af0f8e8ce6 Mon Sep 17 00:00:00 2001
From: Jean-Christophe <loiseau.jc@gmail.com>
Date: Sun, 30 Mar 2025 15:00:15 +0200
Subject: [PATCH 5/7] Update doc/specs/stdlib_linalg.md

Co-authored-by: Jeremie Vandenplas <jeremie.vandenplas@gmail.com>
---
 doc/specs/stdlib_linalg.md | 1 +
 1 file changed, 1 insertion(+)

diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md
index e7afdd912..ea67014b8 100644
--- a/doc/specs/stdlib_linalg.md
+++ b/doc/specs/stdlib_linalg.md
@@ -1901,6 +1901,7 @@ Given a matrix \(A\), this function compute its matrix exponential \(E = \exp(A)
 `order` (optional): Shall be a non-negative `integer` value specifying the order of the Pade approximation. By default `order=10`. It is an `intent(in)` argument. 
 
 `err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. 
+
 ### Return value
 
 The returned array `E` contains the Pade approximation of \(\exp(A)\).

From 4310db503b39e6d572921d85ed35a57f993e0bcf Mon Sep 17 00:00:00 2001
From: Jean-Christophe <loiseau.jc@gmail.com>
Date: Mon, 31 Mar 2025 09:06:50 +0200
Subject: [PATCH 6/7] Replace matmul with gemm.

---
 src/stdlib_linalg_matrix_functions.fypp | 9 +++++++--
 1 file changed, 7 insertions(+), 2 deletions(-)

diff --git a/src/stdlib_linalg_matrix_functions.fypp b/src/stdlib_linalg_matrix_functions.fypp
index 7bd8b8fcb..d05cba564 100644
--- a/src/stdlib_linalg_matrix_functions.fypp
+++ b/src/stdlib_linalg_matrix_functions.fypp
@@ -2,6 +2,7 @@
 #:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
 submodule (stdlib_linalg) stdlib_linalg_matrix_functions
     use stdlib_linalg_constants
+    use stdlib_linalg_blas, only: gemm
     use stdlib_linalg_lapack, only: gesv
     use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, &
          LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR
@@ -97,10 +98,14 @@ contains
             call linalg_error_handling(err0, err)
         end block
 
-        ! This loop should eventually be replaced by a fast matrix_power function.
+        ! Matrix squaring.
+        block
+        ${rt}$ :: E_tmp(n, n)
         do k = 1, s
-            E = matmul(E, E)
+            E_tmp = E
+            call gemm("N", "N", n, n, n, one_${ri}$, E_tmp, n, E_tmp, n, zero_${ri}$, E, n)
         enddo
+        end block
         return
     contains
         elemental subroutine handle_gesv_info(info,lda,n,nrhs,err)

From 59ffb202d7626a7b94d834d616955cf63274483a Mon Sep 17 00:00:00 2001
From: Jean-Christophe <loiseau.jc@gmail.com>
Date: Mon, 31 Mar 2025 09:31:18 +0200
Subject: [PATCH 7/7] Error handling tests.

---
 src/stdlib_linalg_matrix_functions.fypp | 17 +++++++-----
 test/linalg/test_linalg_expm.fypp       | 36 +++++++++++++++++++++++++
 2 files changed, 46 insertions(+), 7 deletions(-)

diff --git a/src/stdlib_linalg_matrix_functions.fypp b/src/stdlib_linalg_matrix_functions.fypp
index d05cba564..7962068b9 100644
--- a/src/stdlib_linalg_matrix_functions.fypp
+++ b/src/stdlib_linalg_matrix_functions.fypp
@@ -40,17 +40,20 @@ contains
 
         ! Deal with optional args.
         order_ = 10 ; if (present(order)) order_ = order
+        print *, "inside expm :", order_
 
         ! Problem's dimension.
         m = size(A, 1) ; n = size(A, 2)
 
         if (m /= n) then
-            err = linalg_state_type(this,LINALG_VALUE_ERROR,'Invalid matrix size A=',[m, n])
+            err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'Invalid matrix size A=',[m, n])
             call linalg_error_handling(err0, err)
+            return
         else if (order_ < 0) then
-            err = linalg_state_type(this, LINALG_VALUE_ERROR, 'Order of Pade approximation &
+            err0 = linalg_state_type(this, LINALG_VALUE_ERROR, 'Order of Pade approximation &
                                     needs to be positive, order=', order_)
             call linalg_error_handling(err0, err)
+            return
         endif
 
         ! Compute the L-infinity norm.
@@ -100,11 +103,11 @@ contains
 
         ! Matrix squaring.
         block
-        ${rt}$ :: E_tmp(n, n)
-        do k = 1, s
-            E_tmp = E
-            call gemm("N", "N", n, n, n, one_${ri}$, E_tmp, n, E_tmp, n, zero_${ri}$, E, n)
-        enddo
+            ${rt}$ :: E_tmp(n, n)
+            do k = 1, s
+                E_tmp = E
+                call gemm("N", "N", n, n, n, one_${ri}$, E_tmp, n, E_tmp, n, zero_${ri}$, E, n)
+            enddo
         end block
         return
     contains
diff --git a/test/linalg/test_linalg_expm.fypp b/test/linalg/test_linalg_expm.fypp
index 51d89a2ca..6ff081ff8 100644
--- a/test/linalg/test_linalg_expm.fypp
+++ b/test/linalg/test_linalg_expm.fypp
@@ -5,6 +5,8 @@ module test_linalg_expm
     use testdrive, only: error_type, check, new_unittest, unittest_type
     use stdlib_linalg_constants
     use stdlib_linalg, only: expm, eye, mnorm
+    use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, &
+         LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR
 
     implicit none (type,external)
     
@@ -21,6 +23,7 @@ module test_linalg_expm
         
         #:for rk,rt,ri in RC_KINDS_TYPES
         tests = [tests, new_unittest("expm_${ri}$",test_expm_${ri}$)]
+        tests = [tests, new_unittest("Error-handling expm_${ri}$",test_error_handling_expm_${ri}$)]
         #:endfor
 
     end subroutine test_expm_computation
@@ -61,6 +64,39 @@ module test_linalg_expm
     end subroutine test_expm_${ri}$    
     #:endfor    
 
+    !> Test error handler.
+    #:for rk,rt,ri in RC_KINDS_TYPES
+    subroutine test_error_handling_expm_${ri}$(error)
+        type(error_type), allocatable, intent(out) :: error
+        ! Problem dimension.
+        integer(ilp), parameter :: n = 5, m = 6
+        ! Test matrix.
+        ${rt}$ :: A(n, n), E(n, n)
+        type(linalg_state_type) :: err
+        integer(ilp) :: i
+
+        ! Initialize matrix.
+        A = 0.0_${rk}$
+        do i = 1, n-1
+            A(i, i+1) = m*1.0_${rk}$
+        enddo
+
+        ! Compute matrix exponential.
+        E = expm(A, order=-1, err=err)
+        ! Check result.
+        call check(error, err%error(), "Negative Pade order")
+        if (allocated(error)) return
+
+        ! Compute matrix exponential.
+        E = expm(A(:n, :n-1), err=err)
+        ! Check result.
+        call check(error, err%error(), "Invalid matrix size")
+        if (allocated(error)) return
+
+        return        
+    end subroutine test_error_handling_expm_${ri}$    
+    #:endfor    
+
 end module test_linalg_expm
 
 program test_expm