base implementation

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