option for pre-allocated pivot indices

Author: perazz

src/stdlib_linalg.fypp

@@ -575,9 +575,11 @@ module stdlib_linalg
   interface invert
     #:for rk,rt,ri in RC_KINDS_TYPES
     #:if rk!="xdp" 
-    module subroutine stdlib_linalg_invert_${ri}$(a,err)
+    module subroutine stdlib_linalg_invert_${ri}$(a,pivot,err)
          !> Input matrix a[n,n]
          ${rt}$, intent(inout) :: a(:,:)
+         !> [optional] Storage array for the diagonal pivot indices
+         integer(ilp), optional, intent(inout), target :: pivot(:)         
          !> [optional] state return flag. On error if not requested, the code will stop
          type(linalg_state_type), optional, intent(out) :: err
     end subroutine stdlib_linalg_invert_${ri}$

src/stdlib_linalg_inverse.fypp

@@ -34,31 +34,40 @@ submodule (stdlib_linalg) stdlib_linalg_inverse
      #:for rk,rt,ri in RC_KINDS_TYPES
      #:if rk!="xdp" 
      ! Compute the in-place square matrix inverse of a
-     module subroutine stdlib_linalg_invert_${ri}$(a,err)
-         !> Input matrix a[n,n]
+     module subroutine stdlib_linalg_invert_${ri}$(a,pivot,err)
+         !> Input matrix a[n,n]. On return, A is destroyed and replaced by the inverse
          ${rt}$, intent(inout) :: a(:,:)
+         !> [optional] Storage array for the diagonal pivot indices
+         integer(ilp), optional, intent(inout), target :: pivot(:)
          !> [optional] state return flag. On error if not requested, the code will stop
          type(linalg_state_type), optional, intent(out) :: err
 
          !> Local variables
          type(linalg_state_type) :: err0
-         integer(ilp) :: lda,n,info,nb,lwork
-         integer(ilp), allocatable :: ipiv(:)
+         integer(ilp) :: lda,n,info,nb,lwork,npiv
+         integer(ilp), pointer :: ipiv(:)
          ${rt}$, allocatable :: work(:)
 
          !> Problem sizes
          lda  = size(a,1,kind=ilp)
          n    = size(a,2,kind=ilp)
 
-         if (lda<1 .or. n<1 .or. lda/=n) then
-            err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size: a=[',lda,',',n,']')
+         ! Has a pre-allocated pivots storage array been provided? 
+         if (present(pivot)) then 
+            ipiv => pivot
+         else
+            allocate(ipiv(n))
+         endif
+         npiv = size(ipiv,kind=ilp)
+
+         if (lda<1 .or. n<1 .or. lda/=n .or. npiv<n) then
+            err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size: a=',[lda,n], &
+                                                             ', pivot=',npiv)
+            if (.not.present(pivot)) deallocate(ipiv)
             call linalg_error_handling(err0,err)
             return
          end if
 
-         ! Pivot indices
-         allocate(ipiv(n))
-
          ! Factorize matrix (overwrite result)
          call getrf(lda,n,a,lda,ipiv,info)
 
@@ -80,6 +89,7 @@ submodule (stdlib_linalg) stdlib_linalg_inverse
          call handle_getri_info(info,lda,n,err0)
 
          ! Process output and return
+         if (.not.present(pivot)) deallocate(ipiv)
          call linalg_error_handling(err0,err)
 
      end subroutine stdlib_linalg_invert_${ri}$
@@ -97,7 +107,7 @@ submodule (stdlib_linalg) stdlib_linalg_inverse
          allocate(inva,source=a)
 
          !> Compute matrix inverse
-         call stdlib_linalg_invert_${ri}$(inva,err)
+         call stdlib_linalg_invert_${ri}$(inva,err=err)
 
      end function stdlib_linalg_inverse_${ri}$