subroutine version, non in-place

Author: perazz

doc/specs/stdlib_linalg.md

@@ -948,12 +948,16 @@ The solver is based on LAPACK's `*GETRF` and `*GETRI` backends.
 
 ### Syntax
 
-`call ` [[stdlib_linalg(module):invert(interface)]] `(a, [, pivot] [, err])`
+`call ` [[stdlib_linalg(module):invert(interface)]] `(a, [,inva] [, pivot] [, err])`
 
 ### Arguments
 
-`a`: Shall be a rank-2, square, `real` or `complex` array containing the coefficient matrix. It is an `intent(inout)` argument.
-On output, it is replaced by the inverse of `a`.
+`a`: Shall be a rank-2, square, `real` or `complex` array containing the coefficient matrix. 
+If `inva` is provided, it is an `intent(in)` argument.
+If `inva` is not provided, it is an `intent(inout)` argument: on output, it is replaced by the inverse of `a`. 
+
+`inva` (optional): Shall be a rank-2, square, `real` or `complex` array with the same size, and kind as `a`. 
+On output, it contains the inverse of `a`.
 
 `pivot` (optional): Shall be a rank-1 array of the same kind and matrix dimension as `a`, providing storage for the diagonal pivot indices. It is an `intent(inout)` arguments, and returns the diagonal pivot indices. 
 
@@ -964,6 +968,7 @@ On output, it is replaced by the inverse of `a`.
 Replaces matrix \( A \) with its inverse, \(A^{-1}\).
 
 Raises `LINALG_ERROR` if the matrix is singular or has invalid size.
+Raises `LINALG_VALUE_ERROR` if `inva` and `a` do not have the same size.
 If `err` is not present, exceptions trigger an `error stop`.
 
 ### Example

src/stdlib_linalg.fypp

@@ -604,25 +604,38 @@ module stdlib_linalg
     !!
     !!### Description
     !!     
-    !! This subroutine interface provides a way to compute the inverse of a matrix in-place.    
+    !! This subroutine interface provides a way to compute the inverse of a matrix.    
     !! Supported data types include `real` and `complex`. 
-    !! On completion, matrix `a` is replaced by the inverse matrix. Pre-allocated storage may be provided
-    !! for the array of pivot indices, `pivot`. If all pre-allocated work spaces are provided, no internal 
-    !! memory allocations take place when using this interface.             
+    !! The user may provide a unique matrix argument `a`. In this case, `a` is replaced by the inverse matrix.
+    !! on output. Otherwise, one may provide two separate arguments: an input matrix `a` and an output matrix 
+    !! `inva`. In this case, `a` will not be modified, and the inverse is returned in `inva`.
+    !! Pre-allocated storage may be provided for the array of pivot indices, `pivot`. If all pre-allocated 
+    !! work spaces are provided, no internal memory allocations take place when using this interface.             
     !!
     !!@note The provided subroutines are intended for square matrices.
     !!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``).
     !!      
     #:for rk,rt,ri in RC_KINDS_TYPES
     #:if rk!="xdp" 
-    module subroutine stdlib_linalg_invert_${ri}$(a,pivot,err)
+    module subroutine stdlib_linalg_invert_inplace_${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}$
+    end subroutine stdlib_linalg_invert_inplace_${ri}$
+    ! Compute the square matrix inverse of a
+    module subroutine stdlib_linalg_invert_split_${ri}$(a,inva,pivot,err)
+         !> Input matrix a[n,n].
+         ${rt}$, intent(in) :: a(:,:)
+         !> Inverse matrix a[n,n]. 
+         ${rt}$, intent(out) :: inva(:,:)         
+         !> [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_split_${ri}$
     #:endif
     #:endfor
   end interface invert

src/stdlib_linalg_inverse.fypp

@@ -34,7 +34,7 @@ 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,pivot,err)
+     module subroutine stdlib_linalg_invert_inplace_${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
@@ -92,7 +92,43 @@ submodule (stdlib_linalg) stdlib_linalg_inverse
          if (.not.present(pivot)) deallocate(ipiv)
          call linalg_error_handling(err0,err)
 
-     end subroutine stdlib_linalg_invert_${ri}$
+     end subroutine stdlib_linalg_invert_inplace_${ri}$
+
+     ! Compute the square matrix inverse of a
+     module subroutine stdlib_linalg_invert_split_${ri}$(a,inva,pivot,err)
+         !> Input matrix a[n,n].
+         ${rt}$, intent(in) :: a(:,:)
+         !> Inverse matrix a[n,n]. 
+         ${rt}$, intent(out) :: inva(:,:)         
+         !> [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
+         
+         type(linalg_state_type) :: err0
+         integer(ilp) :: sa(2),sinva(2)
+         
+         sa    = shape(a,kind=ilp)
+         sinva = shape(inva,kind=ilp)
+         
+         if (any(sa/=sinva)) then
+ 
+             err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size: a=',sa,' inva=',sinva)
+            
+         else
+            
+             !> Copy data in 
+             inva = a
+
+             !> Compute matrix inverse
+             call stdlib_linalg_invert_inplace_${ri}$(inva,err=err0)      
+            
+         end if         
+         
+         ! Process output and return
+         call linalg_error_handling(err0,err)   
+         
+     end subroutine stdlib_linalg_invert_split_${ri}$
 
      ! Invert matrix in place
      module function stdlib_linalg_inverse_${ri}$(a,err) result(inva)
@@ -107,7 +143,7 @@ submodule (stdlib_linalg) stdlib_linalg_inverse
          allocate(inva,source=a)
 
          !> Compute matrix inverse
-         call stdlib_linalg_invert_${ri}$(inva,err=err)
+         call stdlib_linalg_invert_inplace_${ri}$(inva,err=err)
 
      end function stdlib_linalg_inverse_${ri}$