base implementation

Author: perazz

src/CMakeLists.txt

@@ -24,6 +24,7 @@ set(fppFiles
     stdlib_linalg_outer_product.fypp
     stdlib_linalg_kronecker.fypp
     stdlib_linalg_cross_product.fypp
+    stdlib_linalg_solve.fypp    
     stdlib_linalg_state.fypp 
     stdlib_optval.fypp
     stdlib_selection.fypp

src/stdlib_linalg_solve.fypp

@@ -0,0 +1,122 @@
+#:include "common.fypp"
+#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
+#:set RHS_SUFFIX = ["one","many"]
+#:set RHS_SYMBOL = [ranksuffix(r) for r in [1,2]]
+#:set RHS_EMPTY  = [emptyranksuffix(r) for r in [1,2]]
+#:set ALL_RHS    = list(zip(RHS_SYMBOL,RHS_SUFFIX,RHS_EMPTY))
+module stdlib_linalg_solve
+     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(type,external)
+     private
+
+     !> Solve a linear system
+     public :: solve
+
+     ! NumPy: solve(a, b)
+     ! Scipy: solve(a, b, lower=False, overwrite_a=False, overwrite_b=False, check_finite=True, assume_a='gen', transposed=False)[source]#
+     ! IMSL: lu_solve(a, b, transpose=False)
+
+     interface solve
+        #:for nd,ndsuf,nde in ALL_RHS
+        #:for rk,rt,ri in RC_KINDS_TYPES
+        module procedure stdlib_linalg_${ri}$solve${ndsuf}$
+        #:endfor
+        #:endfor
+     end interface solve
+
+
+     contains
+
+     #:for nd,ndsuf,nde in ALL_RHS
+     #:for rk,rt,ri in RC_KINDS_TYPES
+     ! Compute the solution to a real system of linear equations A * X = B
+     function stdlib_linalg_${ri}$solve${ndsuf}$(a,b,overwrite_a,err) result(x)
+         !> Input matrix a[n,n]
+         ${rt}$, intent(inout), target :: a(:,:)
+         !> Right hand side vector or array, b[n] or b[n,nrhs]
+         ${rt}$, intent(in) :: b${nd}$
+         !> [optional] Can A data be overwritten and destroyed?
+         logical(lk), optional, intent(in) :: overwrite_a
+         !> [optional] state return flag. On error if not requested, the code will stop
+         type(linalg_state_type), optional, intent(out) :: err
+         !> Result array/matrix x[n] or x[n,nrhs]
+         ${rt}$, allocatable, target :: x${nd}$
+         
+         !> Local variables
+         type(linalg_state_type) :: err0
+         integer(ilp) :: lda,n,ldb,nrhs,info
+         integer(ilp), allocatable :: ipiv(:)
+         logical(lk) :: copy_a
+         ${rt}$, pointer :: xmat(:,:),amat(:,:)
+         character(*), parameter :: this = 'solve'
+
+         !> Problem sizes
+         lda  = size(a,1,kind=ilp)
+         n    = size(a,2,kind=ilp)
+         ldb  = size(b,1,kind=ilp)
+         nrhs = size(b  ,kind=ilp)/ldb
+
+         if (lda<1 .or. n<1 .or. ldb<1 .or. lda/=n .or. ldb/=n) then
+            err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid sizes: a=[',lda,',',n,'],',&
+                                                                       'b=[',ldb,',',nrhs,']')
+            allocate(x${nde}$)
+            goto 1
+         end if
+
+         ! Can A be overwritten? By default, do not overwrite
+         if (present(overwrite_a)) then
+            copy_a = .not.overwrite_a
+         else
+            copy_a = .true._lk
+         endif
+
+         ! Pivot indices
+         allocate(ipiv(n))
+
+         ! Initialize a matrix temporary
+         if (copy_a) then
+            allocate(amat(lda,n),source=a)
+         else
+            amat => a
+         endif
+
+         ! Initialize solution with the rhs
+         allocate(x,source=b)
+         xmat(1:n,1:nrhs) => x
+
+         ! Solve system
+         call gesv(n,nrhs,amat,lda,ipiv,xmat,ldb,info)
+
+         ! Process output
+         select case (info)
+            case (0)
+                ! Success
+            case (-1)
+                err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid problem size n=',n)
+            case (-2)
+                err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid rhs size n=',nrhs)
+            case (-4)
+                err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size a=[',lda,',',n,']')
+            case (-7)
+                err0 = linalg_state_type(this,LINALG_ERROR,'invalid matrix size a=[',lda,',',n,']')
+            case (1:)
+                err0 = linalg_state_type(this,LINALG_ERROR,'singular matrix')
+            case default
+                err0 = linalg_state_type(this,LINALG_INTERNAL_ERROR,'catastrophic error')
+         end select
+
+         if (.not.copy_a) deallocate(amat)
+
+         ! Process output and return
+         1 call linalg_error_handling(err0,err)
+
+     end function stdlib_linalg_${ri}$solve${ndsuf}$
+
+
+     #:endfor
+     #:endfor
+
+end module stdlib_linalg_solve