Removed one block and one temporary variable allocation.

Author: loiseaujc

src/stdlib_linalg.fypp

@@ -1693,7 +1693,7 @@ module stdlib_linalg
     !! `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.
+    !! via the `order` argument that 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.
@@ -1738,7 +1738,7 @@ module stdlib_linalg
     !! `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.
+    !! via the `order` argument that 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.

src/stdlib_linalg_matrix_functions.fypp

@@ -30,13 +30,13 @@ contains
     module subroutine stdlib_linalg_${ri}$_expm(A, E, order, err)
         !> Input matrix A(n, n).
         ${rt}$, intent(in) :: A(:, :)
-        !> [optional] Order of the Pade approximation.
+        !> Exponential of the input matrix E = exp(A).
+        ${rt}$, intent(out) :: E(:, :)
+         !> [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}$, intent(out) :: E(:, :)
-        
+       
         type(linalg_state_type) :: err0
         integer(ilp) :: lda, n, lde, ne
 
@@ -68,7 +68,7 @@ contains
         type(linalg_state_type), optional, intent(out) :: err
 
         ! Internal variables.
-        ${rt}$, allocatable     :: A2(:, :), Q(:, :), X(:, :)
+        ${rt}$, allocatable     :: A2(:, :), Q(:, :), X(:, :), X_tmp(:, :)
         real(${rk}$)            :: a_norm, c
         integer(ilp)            :: m, n, ee, k, s, order_, i, j
         logical(lk)             :: p
@@ -105,32 +105,29 @@ contains
             enddo
 
             ! Iteratively compute the Pade approximation.
-            block
-                ${rt}$, allocatable :: X_tmp(:, :)
-                p = .true.
-                do k = 2, order_
-                    c = c * (order_ - k + 1) / (k * (2*order_ - k + 1))
-                    X_tmp = X
-                    #:if rt.startswith('complex')
-                    call gemm("N", "N", n, n, n, one_c${rk}$, A2, n, X_tmp, n, zero_c${rk}$, X, n)
-                    #:else
-                    call gemm("N", "N", n, n, n, one_${rk}$, A2, n, X_tmp, n, zero_${rk}$, X, n)
-                    #:endif
+            p = .true.
+            do k = 2, order_
+                c = c * (order_ - k + 1) / (k * (2*order_ - k + 1))
+                X_tmp = X
+                #:if rt.startswith('complex')
+                call gemm("N", "N", n, n, n, one_c${rk}$, A2, n, X_tmp, n, zero_c${rk}$, X, n)
+                #:else
+                call gemm("N", "N", n, n, n, one_${rk}$, A2, n, X_tmp, n, zero_${rk}$, X, n)
+                #:endif
+                do concurrent(i=1:n, j=1:n)
+                    A(i, j) = A(i, j) + c*X(i, j)       ! E = E + c*X
+                enddo
+                if (p) then
                     do concurrent(i=1:n, j=1:n)
-                        A(i, j) = A(i, j) + c*X(i, j)       ! E = E + c*X
+                        Q(i, j) = Q(i, j) + c*X(i, j)   ! Q = Q + c*X
                     enddo
-                    if (p) then
-                        do concurrent(i=1:n, j=1:n)
-                            Q(i, j) = Q(i, j) + c*X(i, j)   ! Q = Q + c*X
-                        enddo
-                    else
-                        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
-            end block
+                else
+                    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
 
             block
                 integer(ilp) :: ipiv(n), info
@@ -139,17 +136,14 @@ contains
             end block
 
             ! Matrix squaring.
-            block
-                ${rt}$, allocatable :: E_tmp(:, :)
-                do k = 1, s
-                    E_tmp = A
-                    #:if rt.startswith('complex')
-                    call gemm("N", "N", n, n, n, one_c${rk}$, E_tmp, n, E_tmp, n, zero_c${rk}$, A, n)
-                    #:else
-                    call gemm("N", "N", n, n, n, one_${rk}$, E_tmp, n, E_tmp, n, zero_${rk}$, A, n)
-                    #:endif
-                enddo
-            end block
+            do k = 1, s
+                X = A ! Re-use X to minimize allocations.
+                #:if rt.startswith('complex')
+                call gemm("N", "N", n, n, n, one_c${rk}$, X, n, X, n, zero_c${rk}$, A, n)
+                #:else
+                call gemm("N", "N", n, n, n, one_${rk}$, X, n, X, n, zero_${rk}$, A, n)
+                #:endif
+            enddo
         endif
 
         call linalg_error_handling(err0, err)