linalg: solve

doc/specs/stdlib_linalg.md

@@ -600,6 +600,56 @@ Specifically, upper Hessenberg matrices satisfy `a_ij = 0` when `j < i-1`, and l
 {!example/linalg/example_is_hessenberg.f90!}
 ```
 
+## `solve` - Solves a linear matrix equation or a linear system of scalar equations. 
+
+### Status
+
+Experimental
+
+### Description
+
+This function computes the solution to a linear matrix equation \( A \cdot x = b \), where \( A \) is a square, full-rank, `real` or `complex` matrix.
+
+Result vector `x` returns the exact solution to within numerical precision, provided that the matrix is not ill-conditioned. The solver is based on LAPACK's `*GESV` backends.
+
+### Syntax
+
+`Pure` interface:
+
+`x = ` [[stdlib_linalg(module):solve(interface)]] `(a, b)`
+
+Expert interface:
+
+`x = ` [[stdlib_linalg(module):solve(interface)]] `(a, b, [, overwrite_a], err])`
+
+### Arguments
+
+Two 
+
+`a`: Shall be a rank-2 `real` or `complex` square array containing the coefficient matrix. It is normally an `intent(in)` argument. If `overwrite_a=.true.`, it is an `intent(inout)` argument and is destroyed by the call. 
+
+`b`: Shall be a rank-1 array of the same kind as `a`, containing the right-hand-side vector. It is an `intent(in)` argument.
+
+`overwrite_a` (optional): Shall be an input logical flag. if `.true.`, input matrix a will be used as temporary storage and overwritten. This avoids internal data allocation. This is an `intent(in)` argument.
+
+`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument: the function is not `pure` if this argument is requested.
+
+### Return value
+
+Returns an array value that represents the solution to the linear system of equations.
+
+Raises `LINALG_ERROR` if the matrix is singular to working precision.
+Raises `LINALG_VALUE_ERROR` if the matrix and rhs vectors have invalid/incompatible sizes.
+Exceptions trigger an `error stop`.
+
+### Example
+
+```fortran
+{!example/linalg/example_solve1.f90!}
+
+{!example/linalg/example_solve2.f90!}
+```
+
 ## `det` - Computes the determinant of a square matrix
 
 ### Status

example/linalg/CMakeLists.txt

@@ -18,6 +18,7 @@ ADD_EXAMPLE(state1)
 ADD_EXAMPLE(state2)
 ADD_EXAMPLE(blas_gemv)
 ADD_EXAMPLE(lapack_getrf)
+ADD_EXAMPLE(solve1)
+ADD_EXAMPLE(solve2)
 ADD_EXAMPLE(determinant)
 ADD_EXAMPLE(determinant2)
-

example/linalg/example_solve1.f90

@@ -0,0 +1,26 @@
+program example_solve1
+  use stdlib_linalg_constants, only: sp
+  use stdlib_linalg, only: solve, linalg_state_type 
+  implicit none
+
+  real(sp), allocatable :: A(:,:),b(:),x(:)
+
+  ! Solve a system of 3 linear equations: 
+  !  4x + 3y + 2z =  25
+  ! -2x + 2y + 3z = -10
+  !  3x - 5y + 2z =  -4
+
+  ! Note: Fortran is column-major! -> transpose 
+  A = transpose(reshape([ 4, 3, 2, &
+                         -2, 2, 3, &
+                          3,-5, 2], [3,3])) 
+  b = [25,-10,-4]
+
+  ! Get coefficients of y = coef(1) + x*coef(2) + x^2*coef(3)
+  x = solve(A,b) 
+
+  print *, 'solution: ',x
+  ! 5.0, 3.0, -2.0
+
+end program example_solve1
+

example/linalg/example_solve2.f90

@@ -0,0 +1,26 @@
+program example_solve2
+  use stdlib_linalg_constants, only: sp
+  use stdlib_linalg, only: solve, linalg_state_type 
+  implicit none
+
+  complex(sp), allocatable :: A(:,:),b(:),x(:)
+
+  ! Solve a system of 3 complex linear equations: 
+  !  2x +      iy + 2z = (5-i)
+  ! -ix + (4-3i)y + 6z = i
+  !  4x +      3y +  z = 1
+
+  ! Note: Fortran is column-major! -> transpose 
+  A = transpose(reshape([(2.0, 0.0),(0.0, 1.0),(2.0,0.0), &
+                         (0.0,-1.0),(4.0,-3.0),(6.0,0.0), &
+                         (4.0, 0.0),(3.0, 0.0),(1.0,0.0)] , [3,3])) 
+  b = [(5.0,-1.0),(0.0,1.0),(1.0,0.0)]
+
+  ! Get coefficients of y = coef(1) + x*coef(2) + x^2*coef(3)
+  x = solve(A,b) 
+
+  print *, 'solution: ',x
+  !   (1.0947,0.3674) (-1.519,-0.4539) (1.1784,-0.1078)
+
+end program example_solve2
+

include/common.fypp

@@ -27,11 +27,20 @@
 #:set REAL_KINDS = REAL_KINDS + ["qp"]
 #:endif
 
+#! BLAS/LAPACK initials for each real kind
+#:set REAL_INIT = ["s", "d"]
+#:if WITH_XDP
+#:set REAL_INIT = REAL_INIT + ["x"]
+#:endif
+#:if WITH_QP
+#:set REAL_INIT = REAL_INIT + ["q"]
+#:endif
+
 #! Real types to be considered during templating
 #:set REAL_TYPES = ["real({})".format(k) for k in REAL_KINDS]
 
 #! Collected (kind, type) tuples for real types
-#:set REAL_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES))
+#:set REAL_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_INIT))
 
 #! Complex kinds to be considered during templating
 #:set CMPLX_KINDS = ["sp", "dp"]
@@ -42,11 +51,20 @@
 #:set CMPLX_KINDS = CMPLX_KINDS + ["qp"]
 #:endif
 
+#! BLAS/LAPACK initials for each complex kind
+#:set CMPLX_INIT = ["c", "z"]
+#:if WITH_XDP
+#:set CMPLX_INIT = CMPLX_INIT + ["y"]
+#:endif
+#:if WITH_QP
+#:set CMPLX_INIT = CMPLX_INIT + ["w"]
+#:endif
+
 #! Complex types to be considered during templating
 #:set CMPLX_TYPES = ["complex({})".format(k) for k in CMPLX_KINDS]
 
-#! Collected (kind, type) tuples for complex types
-#:set CMPLX_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES))
+#! Collected (kind, type, initial) tuples for complex types
+#:set CMPLX_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_INIT))
 
 #! Integer kinds to be considered during templating
 #:set INT_KINDS = ["int8", "int16", "int32", "int64"]
@@ -109,6 +127,17 @@
 #{if rank > 0}#(${":" + ",:" * (rank - 1)}$)#{endif}#
 #:enddef
 
+#! Generates an empty array rank suffix.
+#!
+#! Args:
+#!     rank (int): Rank of the variable
+#!
+#! Returns:
+#!     Empty array rank suffix string (e.g. (0,0) if rank = 2)
+#!
+#:def emptyranksuffix(rank)
+#{if rank > 0}#(${"0" + ",0" * (rank - 1)}$)#{endif}#
+#:enddef
 
 #! Joins stripped lines with given character string
 #!

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_determinant.fypp
     stdlib_linalg_state.fypp 
     stdlib_optval.fypp

src/stdlib_linalg.fypp

@@ -1,21 +1,26 @@
 #:include "common.fypp"
 #:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
 #:set RCI_KINDS_TYPES = RC_KINDS_TYPES + INT_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
   !!Provides a support for various linear algebra procedures
   !! ([Specification](../page/specs/stdlib_linalg.html))
-  use stdlib_kinds, only: sp, dp, xdp, qp, lk, &
-    int8, int16, int32, int64
+  use stdlib_kinds, only: xdp, int8, int16, int32, int64
+  use stdlib_linalg_constants, only: sp, dp, qp, lk, ilp
   use stdlib_error, only: error_stop
   use stdlib_optval, only: optval
-  use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling
+  use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling  
   implicit none
   private
 
   public :: det
   public :: operator(.det.)
   public :: diag
   public :: eye
+  public :: solve
   public :: trace
   public :: outer_product
   public :: kronecker_product
@@ -221,6 +226,54 @@ module stdlib_linalg
     #:endfor
   end interface is_hessenberg
 
+  ! Solve linear system system Ax=b.
+  interface solve
+     !! version: experimental 
+     !!
+     !! Solves the linear system \( A \cdot x = b \) for the unknown vector \( x \) from a square matrix \( A \). 
+     !! ([Specification](../page/specs/stdlib_linalg.html#solve-solves-a-linear-matrix-equation-or-a-linear-system-of-scalar-equations))
+     !!
+     !!### Summary 
+     !! Interface for solving a linear system arising from a general matrix.
+     !!
+     !!### Description
+     !! 
+     !! This interface provides methods for computing the solution of a linear matrix system.
+     !! Supported data types include `real` and `complex`. No assumption is made on the matrix 
+     !! structure. 
+     !! The function can solve simultaneously either one (from a 1-d right-hand-side vector `b(:)`) 
+     !! or several (from a 2-d right-hand-side vector `b(:,:)`) systems.
+     !! 
+     !!@note The solution is based on LAPACK's generic LU decomposition based solvers `*GESV`.
+     !!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``).
+     !!    
+     #:for nd,ndsuf,nde in ALL_RHS
+     #:for rk,rt,ri in RC_KINDS_TYPES
+     #:if rk!="xdp"
+     module 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), intent(out) :: err
+         !> Result array/matrix x[n] or x[n,nrhs]
+         ${rt}$, allocatable, target :: x${nd}$
+     end function stdlib_linalg_${ri}$_solve_${ndsuf}$
+     pure module function stdlib_linalg_${ri}$_pure_solve_${ndsuf}$(a,b) result(x)
+         !> Input matrix a[n,n]
+         ${rt}$, intent(in), target :: a(:,:)
+         !> Right hand side vector or array, b[n] or b[n,nrhs]
+         ${rt}$, intent(in) :: b${nd}$
+         !> Result array/matrix x[n] or x[n,nrhs]
+         ${rt}$, allocatable, target :: x${nd}$
+     end function stdlib_linalg_${ri}$_pure_solve_${ndsuf}$
+     #:endif
+     #:endfor
+     #:endfor    
+  end interface solve
 
   interface det
     !! version: experimental