From 06bfe4be804cb9e8619d4e77cd1d740087113e7c Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 25 Apr 2024 14:20:29 +0200 Subject: [PATCH 01/29] templated BLAS/LAPACK initials --- include/common.fypp | 35 ++++++++++++++++++++++++++++++++--- 1 file changed, 32 insertions(+), 3 deletions(-) diff --git a/include/common.fypp b/include/common.fypp index 1683239e4..ed1cf2b4e 100644 --- a/include/common.fypp +++ b/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 #! From a2afe6b6083d493a32f003c842ec027d2c2882e8 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 25 Apr 2024 14:49:23 +0200 Subject: [PATCH 02/29] base implementation --- src/CMakeLists.txt | 1 + src/stdlib_linalg_solve.fypp | 122 +++++++++++++++++++++++++++++++++++ 2 files changed, 123 insertions(+) create mode 100644 src/stdlib_linalg_solve.fypp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 165259db3..7d8a6a014 100644 --- a/src/CMakeLists.txt +++ b/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 diff --git a/src/stdlib_linalg_solve.fypp b/src/stdlib_linalg_solve.fypp new file mode 100644 index 000000000..4de00e463 --- /dev/null +++ b/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 From d929077ae7e9289e69cf617f42c64b7cab1e5b44 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 25 Apr 2024 14:50:27 +0200 Subject: [PATCH 03/29] exclude `xdp` --- src/stdlib_linalg_solve.fypp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/stdlib_linalg_solve.fypp b/src/stdlib_linalg_solve.fypp index 4de00e463..b088e956d 100644 --- a/src/stdlib_linalg_solve.fypp +++ b/src/stdlib_linalg_solve.fypp @@ -22,7 +22,9 @@ module stdlib_linalg_solve interface solve #:for nd,ndsuf,nde in ALL_RHS #:for rk,rt,ri in RC_KINDS_TYPES + #:if rk!="xdp" module procedure stdlib_linalg_${ri}$solve${ndsuf}$ + #:endif #:endfor #:endfor end interface solve @@ -32,6 +34,7 @@ module stdlib_linalg_solve #:for nd,ndsuf,nde in ALL_RHS #:for rk,rt,ri in RC_KINDS_TYPES + #:if rk!="xdp" ! 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] @@ -115,7 +118,7 @@ module stdlib_linalg_solve end function stdlib_linalg_${ri}$solve${ndsuf}$ - + #:endif #:endfor #:endfor From 5c817c8ee3cd544f2028eb49f28d57e945a621b3 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 25 Apr 2024 14:59:16 +0200 Subject: [PATCH 04/29] `pure` interfaces --- src/stdlib_linalg_solve.fypp | 118 +++++++++++++++++++++++++++-------- 1 file changed, 91 insertions(+), 27 deletions(-) diff --git a/src/stdlib_linalg_solve.fypp b/src/stdlib_linalg_solve.fypp index b088e956d..59fad0f49 100644 --- a/src/stdlib_linalg_solve.fypp +++ b/src/stdlib_linalg_solve.fypp @@ -10,7 +10,6 @@ module stdlib_linalg_solve 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 @@ -24,13 +23,40 @@ module stdlib_linalg_solve #:for rk,rt,ri in RC_KINDS_TYPES #:if rk!="xdp" module procedure stdlib_linalg_${ri}$solve${ndsuf}$ + module procedure stdlib_linalg_${ri}$_pure_solve${ndsuf}$ #:endif #:endfor #:endfor end interface solve - + + + character(*), parameter :: this = 'solve' contains + + elemental subroutine handle_gesv_info(info,lda,n,nrhs,err) + integer(ilp), intent(in) :: info,lda,n,nrhs + type(linalg_state_type), intent(out) :: err + + ! Process output + select case (info) + case (0) + ! Success + case (-1) + err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid problem size n=',n) + case (-2) + err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid rhs size n=',nrhs) + case (-4) + err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size a=',[lda,n]) + case (-7) + err = linalg_state_type(this,LINALG_ERROR,'invalid matrix size a=',[lda,n]) + case (1:) + err = linalg_state_type(this,LINALG_ERROR,'singular matrix') + case default + err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'catastrophic error') + end select + + end subroutine handle_gesv_info #:for nd,ndsuf,nde in ALL_RHS #:for rk,rt,ri in RC_KINDS_TYPES @@ -44,7 +70,7 @@ module stdlib_linalg_solve !> [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 + type(linalg_state_type), intent(out) :: err !> Result array/matrix x[n] or x[n,nrhs] ${rt}$, allocatable, target :: x${nd}$ @@ -53,8 +79,7 @@ module stdlib_linalg_solve integer(ilp) :: lda,n,ldb,nrhs,info integer(ilp), allocatable :: ipiv(:) logical(lk) :: copy_a - ${rt}$, pointer :: xmat(:,:),amat(:,:) - character(*), parameter :: this = 'solve' + ${rt}$, pointer :: xmat(:,:),amat(:,:) !> Problem sizes lda = size(a,1,kind=ilp) @@ -62,11 +87,12 @@ module stdlib_linalg_solve 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,']') + if (any([lda,n,ldb]<1) .or. any([lda,ldb]/=n)) then + err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid sizes: a=',[lda,n], & + ', b=',[ldb,nrhs]) allocate(x${nde}$) - goto 1 + call linalg_error_handling(err0,err) + return end if ! Can A be overwritten? By default, do not overwrite @@ -94,30 +120,68 @@ module stdlib_linalg_solve 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 + call handle_gesv_info(info,lda,n,nrhs,err0) - if (.not.copy_a) deallocate(amat) + if (copy_a) deallocate(amat) ! Process output and return - 1 call linalg_error_handling(err0,err) + call linalg_error_handling(err0,err) end function stdlib_linalg_${ri}$solve${ndsuf}$ + ! Compute the solution to a real system of linear equations A * X = B (pure interface) + pure 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}$ + + !> Local variables + type(linalg_state_type) :: err0 + integer(ilp) :: lda,n,ldb,nrhs,info + integer(ilp), allocatable :: ipiv(:) + ${rt}$, pointer :: xmat(:,:) + ${rt}$, allocatable :: 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 (any([lda,n,ldb]<1) .or. any([lda,ldb]/=n)) then + err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid sizes: a=',[lda,n], & + ', b=',[ldb,nrhs]) + allocate(x${nde}$) + call linalg_error_handling(err0) + return + end if + + ! Pivot indices + allocate(ipiv(n)) + + ! Initialize a matrix temporary + allocate(amat,source=a) + + ! 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 + call handle_gesv_info(info,lda,n,nrhs,err0) + + deallocate(amat) + + ! Process output and return + call linalg_error_handling(err0) + + end function stdlib_linalg_${ri}$_pure_solve${ndsuf}$ #:endif #:endfor #:endfor From 4cef2acc7032a2e283128d04a9b498cba4a8d9cb Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 25 Apr 2024 15:02:07 +0200 Subject: [PATCH 05/29] cleanup names --- src/stdlib_linalg_solve.fypp | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/src/stdlib_linalg_solve.fypp b/src/stdlib_linalg_solve.fypp index 59fad0f49..970d484be 100644 --- a/src/stdlib_linalg_solve.fypp +++ b/src/stdlib_linalg_solve.fypp @@ -22,8 +22,8 @@ module stdlib_linalg_solve #:for nd,ndsuf,nde in ALL_RHS #:for rk,rt,ri in RC_KINDS_TYPES #:if rk!="xdp" - module procedure stdlib_linalg_${ri}$solve${ndsuf}$ - module procedure stdlib_linalg_${ri}$_pure_solve${ndsuf}$ + module procedure stdlib_linalg_${ri}$_solve_${ndsuf}$ + module procedure stdlib_linalg_${ri}$_pure_solve_${ndsuf}$ #:endif #:endfor #:endfor @@ -62,7 +62,7 @@ module stdlib_linalg_solve #:for rk,rt,ri in RC_KINDS_TYPES #:if rk!="xdp" ! 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) + 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] @@ -127,10 +127,10 @@ module stdlib_linalg_solve ! Process output and return call linalg_error_handling(err0,err) - end function stdlib_linalg_${ri}$solve${ndsuf}$ + end function stdlib_linalg_${ri}$_solve_${ndsuf}$ ! Compute the solution to a real system of linear equations A * X = B (pure interface) - pure function stdlib_linalg_${ri}$_pure_solve${ndsuf}$(a,b) result(x) + pure 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] @@ -144,7 +144,6 @@ module stdlib_linalg_solve integer(ilp), allocatable :: ipiv(:) ${rt}$, pointer :: xmat(:,:) ${rt}$, allocatable :: amat(:,:) - character(*), parameter :: this = 'solve' !> Problem sizes lda = size(a,1,kind=ilp) @@ -181,7 +180,7 @@ module stdlib_linalg_solve ! Process output and return call linalg_error_handling(err0) - end function stdlib_linalg_${ri}$_pure_solve${ndsuf}$ + end function stdlib_linalg_${ri}$_pure_solve_${ndsuf}$ #:endif #:endfor #:endfor From 7b7c051c08486cfe64a81bdff99c795aafc0e336 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 25 Apr 2024 15:07:31 +0200 Subject: [PATCH 06/29] submodule --- src/stdlib_linalg.fypp | 43 +++++++++++++++++++++++++++++++++--- src/stdlib_linalg_solve.fypp | 36 ++++++++---------------------- 2 files changed, 49 insertions(+), 30 deletions(-) diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index 3ed905c56..81b57c1df 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -1,17 +1,24 @@ #:include "common.fypp" -#:set RCI_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES + INT_KINDS_TYPES +#: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, & - 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 implicit none private public :: diag public :: eye + public :: solve public :: trace public :: outer_product public :: kronecker_product @@ -214,6 +221,36 @@ module stdlib_linalg #:endfor end interface is_hessenberg + ! Solve linear system system Ax=b. + interface solve + #: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 + contains diff --git a/src/stdlib_linalg_solve.fypp b/src/stdlib_linalg_solve.fypp index 970d484be..b44141ca5 100644 --- a/src/stdlib_linalg_solve.fypp +++ b/src/stdlib_linalg_solve.fypp @@ -4,31 +4,13 @@ #: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 +submodule (stdlib_linalg) stdlib_linalg_solve +!! Solve linear system Ax=b 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) - - !> 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 - #:if rk!="xdp" - module procedure stdlib_linalg_${ri}$_solve_${ndsuf}$ - module procedure stdlib_linalg_${ri}$_pure_solve_${ndsuf}$ - #:endif - #:endfor - #:endfor - end interface solve - character(*), parameter :: this = 'solve' @@ -62,7 +44,7 @@ module stdlib_linalg_solve #:for rk,rt,ri in RC_KINDS_TYPES #:if rk!="xdp" ! 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) + 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] @@ -74,14 +56,14 @@ module stdlib_linalg_solve !> Result array/matrix x[n] or x[n,nrhs] ${rt}$, allocatable, target :: x${nd}$ - !> Local variables + ! 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(:,:) - !> Problem sizes + ! Problem sizes lda = size(a,1,kind=ilp) n = size(a,2,kind=ilp) ldb = size(b,1,kind=ilp) @@ -130,7 +112,7 @@ module stdlib_linalg_solve end function stdlib_linalg_${ri}$_solve_${ndsuf}$ ! Compute the solution to a real system of linear equations A * X = B (pure interface) - pure function stdlib_linalg_${ri}$_pure_solve_${ndsuf}$(a,b) result(x) + 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] @@ -138,14 +120,14 @@ module stdlib_linalg_solve !> Result array/matrix x[n] or x[n,nrhs] ${rt}$, allocatable, target :: x${nd}$ - !> Local variables + ! Local variables type(linalg_state_type) :: err0 integer(ilp) :: lda,n,ldb,nrhs,info integer(ilp), allocatable :: ipiv(:) ${rt}$, pointer :: xmat(:,:) ${rt}$, allocatable :: amat(:,:) - !> Problem sizes + ! Problem sizes lda = size(a,1,kind=ilp) n = size(a,2,kind=ilp) ldb = size(b,1,kind=ilp) @@ -185,4 +167,4 @@ module stdlib_linalg_solve #:endfor #:endfor -end module stdlib_linalg_solve +end submodule stdlib_linalg_solve From 0be918eb873c540213dae5c985659d2bd84afd52 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 25 Apr 2024 15:31:27 +0200 Subject: [PATCH 07/29] `real` and `complex` tests --- test/linalg/CMakeLists.txt | 2 + test/linalg/test_linalg_solve.fypp | 190 +++++++++++++++++++++++++++++ 2 files changed, 192 insertions(+) create mode 100644 test/linalg/test_linalg_solve.fypp diff --git a/test/linalg/CMakeLists.txt b/test/linalg/CMakeLists.txt index 3d590a9d2..ec9bb5624 100644 --- a/test/linalg/CMakeLists.txt +++ b/test/linalg/CMakeLists.txt @@ -2,10 +2,12 @@ set( fppFiles "test_linalg.fypp" "test_blas_lapack.fypp" + "test_linalg_solve.fypp" "test_linalg_matrix_property_checks.fypp" ) fypp_f90("${fyppFlags}" "${fppFiles}" outFiles) ADDTEST(linalg) ADDTEST(linalg_matrix_property_checks) +ADDTEST(linalg_solve) ADDTEST(blas_lapack) diff --git a/test/linalg/test_linalg_solve.fypp b/test/linalg/test_linalg_solve.fypp new file mode 100644 index 000000000..ff115965b --- /dev/null +++ b/test/linalg/test_linalg_solve.fypp @@ -0,0 +1,190 @@ +#:include "common.fypp" +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES +! Test linear system solver +module test_linalg_solve + use stdlib_linalg_constants + use stdlib_linalg_state + use stdlib_linalg, only: solve + use testdrive, only: error_type, check, new_unittest, unittest_type + + implicit none (type,external) + private + + public :: test_linear_systems + + contains + + !> Solve real and complex linear systems + subroutine test_linear_systems(tests) + !> Collection of tests + type(unittest_type), allocatable, intent(out) :: tests(:) + + allocate(tests(0)) + + #:for rk,rt,ri in REAL_KINDS_TYPES + #:if rk!="xdp" + tests = [tests,new_unittest("solve_${ri}$",test_${ri}$_solve), & + new_unittest("solve_${ri}$_multiple",test_${ri}$_solve_multiple)] + #:endif + #:endfor + + #:for ck,ct,ci in CMPLX_KINDS_TYPES + #:if ck!="xdp" + tests = [tests,new_unittest("solve_complex_${ci}$",test_${ci}$_solve), & + new_unittest("solve_2x2_complex_${ci}$",test_2x2_${ci}$_solve)] + #:endif + #:endfor + + end subroutine test_linear_systems + + #:for rk,rt,ri in REAL_KINDS_TYPES + #:if rk!="xdp" + !> Simple linear system + subroutine test_${ri}$_solve(error) + type(error_type), allocatable, intent(out) :: error + + type(linalg_state_type) :: state + + ${rt}$ :: A(3,3) = transpose(reshape([${rt}$ :: 1, 3, 3, & + 1, 3, 4, & + 1, 4, 3], [3,3])) + ${rt}$ :: b (3) = [${rt}$ :: 1, 4, -1] + ${rt}$ :: res(3) = [${rt}$ :: -2, -2, 3] + ${rt}$ :: x(3) + + x = solve(a,b,err=state) + + call check(error,state%ok(),state%print()) + if (allocated(error)) return + + call check(error, all(abs(x-res) Simple linear system with multiple right hand sides + subroutine test_${ri}$_solve_multiple(error) + type(error_type), allocatable, intent(out) :: error + + type(linalg_state_type) :: state + + ${rt}$ :: A(3,3) = transpose(reshape([${rt}$ :: 1,-1, 2, & + 0, 1, 1, & + 1,-1, 3], [3,3])) + ${rt}$ :: b(3,3) = transpose(reshape([${rt}$ :: 0, 1, 2, & + 1,-2,-1, & + 2, 3,-1], [3,3])) + ${rt}$ :: res(3,3) = transpose(reshape([${rt}$ ::-5,-7,10, & + -1,-4, 2, & + 2, 2,-3], [3,3])) + ${rt}$ :: x(3,3) + + x = solve(a,b,err=state) + + call check(error,state%ok(),state%print()) + if (allocated(error)) return + + call check(error, all(abs(x-res) Complex linear system + !> Militaru, Popa, "On the numerical solving of complex linear systems", + !> Int J Pure Appl Math 76(1), 113-122, 2012. + subroutine test_${ri}$_solve(error) + type(error_type), allocatable, intent(out) :: error + + type(linalg_state_type) :: state + + ${rt}$ :: A(5,5), b(5), res(5), x(5) + integer(ilp) :: i + + ! Fill in linear system + A = (0.0_${rk}$,0.0_${rk}$) + + A(1:2,1) = [(19.73_${rk}$,0.0_${rk}$),(0.0_${rk}$,-0.51_${rk}$)] + A(1:3,2) = [(12.11_${rk}$,-1.0_${rk}$),(32.3_${rk}$,7.0_${rk}$),(0.0_${rk}$,-0.51_${rk}$)] + A(1:4,3) = [(0.0_${rk}$,5.0_${rk}$),(23.07_${rk}$,0.0_${rk}$),(70.0_${rk}$,7.3_${rk}$),(1.0_${rk}$,1.1_${rk}$)] + A(2:5,4) = [(0.0_${rk}$,1.0_${rk}$),(3.95_${rk}$,0.0_${rk}$),(50.17_${rk}$,0.0_${rk}$),(0.0_${rk}$,-9.351_${rk}$)] + A(3:5,5) = [(19.0_${rk}$,31.83_${rk}$),(45.51_${rk}$,0.0_${rk}$),(55.0_${rk}$,0.0_${rk}$)] + + b = [(77.38_${rk}$,8.82_${rk}$),(157.48_${rk}$,19.8_${rk}$),(1175.62_${rk}$,20.69_${rk}$),(912.12_${rk}$,-801.75_${rk}$),(550.0_${rk}$,-1060.4_${rk}$)] + + ! Exact result + res = [(3.3_${rk}$,-1.0_${rk}$),(1.0_${rk}$,0.17_${rk}$),(5.5_${rk}$,0.0_${rk}$),(9.0_${rk}$,0.0_${rk}$),(10.0_${rk}$,-17.75_${rk}$)] + + x = solve(a,b,err=state) + + call check(error,state%ok(),state%print()) + if (allocated(error)) return + + call check(error, all(abs(x-res) 2x2 Complex linear system + !> https://math.stackexchange.com/questions/1996540/solving-linear-equation-systems-with-complex-coefficients-and-variables + subroutine test_2x2_${ri}$_solve(error) + type(error_type), allocatable, intent(out) :: error + + type(linalg_state_type) :: state + + ${rt}$ :: A(2,2), b(2), res(2), x(2) + integer(ilp) :: i + + ! Fill in linear system + A(1,:) = [(+1.0_${rk}$,+1.0_${rk}$),(-1.0_${rk}$,0.0_${rk}$)] + A(2,:) = [(+1.0_${rk}$,-1.0_${rk}$),(+1.0_${rk}$,1.0_${rk}$)] + + b = [(0.0_${rk}$,1.0_${rk}$),(1.0_${rk}$,0.0_${rk}$)] + + ! Exact result + res = [(0.5_${rk}$,0.5_${rk}$),(0.0_${rk}$,0.0_${rk}$)] + + x = solve(a,b,err=state) + + call check(error,state%ok(),state%print()) + if (allocated(error)) return + + call check(error, all(abs(x-res) 0) then + write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!" + error stop + end if +end program test_solve + From 9906c9371579a609bc55c69ad2ce6d937f7dfcb3 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 26 Apr 2024 09:20:27 +0200 Subject: [PATCH 08/29] add `real` and `complex` examples --- example/linalg/CMakeLists.txt | 2 ++ example/linalg/example_solve1.f90 | 26 ++++++++++++++++++++++++++ example/linalg/example_solve2.f90 | 26 ++++++++++++++++++++++++++ src/stdlib_linalg.fypp | 3 +++ 4 files changed, 57 insertions(+) create mode 100644 example/linalg/example_solve1.f90 create mode 100644 example/linalg/example_solve2.f90 diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt index 1a5875502..af811450b 100644 --- a/example/linalg/CMakeLists.txt +++ b/example/linalg/CMakeLists.txt @@ -18,3 +18,5 @@ ADD_EXAMPLE(state1) ADD_EXAMPLE(state2) ADD_EXAMPLE(blas_gemv) ADD_EXAMPLE(lapack_getrf) +ADD_EXAMPLE(solve1) +ADD_EXAMPLE(solve2) diff --git a/example/linalg/example_solve1.f90 b/example/linalg/example_solve1.f90 new file mode 100644 index 000000000..0b5c1bbf1 --- /dev/null +++ b/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 + diff --git a/example/linalg/example_solve2.f90 b/example/linalg/example_solve2.f90 new file mode 100644 index 000000000..a761ae7bd --- /dev/null +++ b/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 + diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index 81b57c1df..ffeb72e62 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -31,6 +31,9 @@ module stdlib_linalg public :: is_triangular public :: is_hessenberg + ! Export linalg error handling + public :: linalg_state_type, linalg_error_handling + interface diag !! version: experimental !! From a5d1b8a1c1f814fc190c611b7113c77366aeab99 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 26 Apr 2024 09:38:24 +0200 Subject: [PATCH 09/29] add specs --- doc/specs/stdlib_linalg.md | 50 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 50 insertions(+) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index defe30758..b5e9769cd 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -599,3 +599,53 @@ Specifically, upper Hessenberg matrices satisfy `a_ij = 0` when `j < i-1`, and l ```fortran {!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!} +``` From c1365ff2506adcb1238060a58f1d0beb289c3c61 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 26 Apr 2024 09:52:11 +0200 Subject: [PATCH 10/29] document interface --- src/stdlib_linalg.fypp | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index ffeb72e62..b282137b4 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -226,6 +226,25 @@ module stdlib_linalg ! 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" From 3de68343860d0cda376b3ab89c9f4ae884b824c3 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 26 Apr 2024 10:02:11 +0200 Subject: [PATCH 11/29] fix resolve conflict --- src/stdlib_linalg.fypp | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index 9dd04d81f..d54b228d3 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -36,9 +36,6 @@ module stdlib_linalg ! Export linalg error handling public :: linalg_state_type, linalg_error_handling - ! Export linalg error handling - public :: linalg_state_type, linalg_error_handling - interface diag !! version: experimental !! From 04f126df1a4752e06c94d335cd6cb6b292e4bf51 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Sat, 27 Apr 2024 14:39:28 +0200 Subject: [PATCH 12/29] Update doc/specs/stdlib_linalg.md Co-authored-by: ZUO Zhihua --- doc/specs/stdlib_linalg.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 8de068818..b737471d7 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -628,7 +628,7 @@ 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. +`b`: Shall be a rank-1 or rank-2 array of the same kind as `a`, containing the right-hand-side vector(s). 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. From 7ae0510b1e83dd5dec947f30667fb993f3223ea9 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Sat, 27 Apr 2024 14:40:11 +0200 Subject: [PATCH 13/29] Update doc/specs/stdlib_linalg.md Co-authored-by: ZUO Zhihua --- doc/specs/stdlib_linalg.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index b737471d7..498e186de 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -620,7 +620,7 @@ Result vector `x` returns the exact solution to within numerical precision, prov Expert interface: -`x = ` [[stdlib_linalg(module):solve(interface)]] `(a, b, [, overwrite_a], err])` +`x = ` [[stdlib_linalg(module):solve(interface)]] `(a, b [, overwrite_a], err)` ### Arguments From ed135d9b3194db6ff986ac3568b038e93e1484e1 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Tue, 30 Apr 2024 10:05:43 +0200 Subject: [PATCH 14/29] Update doc/specs/stdlib_linalg.md Co-authored-by: Jeremie Vandenplas --- doc/specs/stdlib_linalg.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 498e186de..0f28ec08f 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -600,7 +600,7 @@ 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. +## `solve` - Solves a linear matrix equation or a linear system of equations. ### Status From e9bf0208b0f09b00278074264405dbb09df3c571 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Tue, 30 Apr 2024 10:05:51 +0200 Subject: [PATCH 15/29] Update doc/specs/stdlib_linalg.md Co-authored-by: Jeremie Vandenplas --- doc/specs/stdlib_linalg.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 0f28ec08f..1cc103e52 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -630,7 +630,7 @@ Two `b`: Shall be a rank-1 or rank-2 array of the same kind as `a`, containing the right-hand-side vector(s). 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. +`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. From 3265f8f0835ba88937114ca26acb78d3a2fdb080 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Tue, 30 Apr 2024 10:05:58 +0200 Subject: [PATCH 16/29] Update src/stdlib_linalg.fypp Co-authored-by: Jeremie Vandenplas --- src/stdlib_linalg.fypp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index d54b228d3..2839e81d5 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -12,7 +12,7 @@ module stdlib_linalg 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 From 77fc5bd1c5e166cae3e450844ab3df75f359e5f2 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Tue, 30 Apr 2024 10:06:04 +0200 Subject: [PATCH 17/29] Update doc/specs/stdlib_linalg.md Co-authored-by: Jeremie Vandenplas --- doc/specs/stdlib_linalg.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 1cc103e52..1888f2edd 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -632,7 +632,7 @@ Two `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. +`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 provided. ### Return value From b04b3d9b2b5f3377f1fb6e23be6ecd25a156185f Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Tue, 30 Apr 2024 10:06:10 +0200 Subject: [PATCH 18/29] Update doc/specs/stdlib_linalg.md Co-authored-by: Jeremie Vandenplas --- doc/specs/stdlib_linalg.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 1888f2edd..f5c7416c0 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -610,7 +610,7 @@ Experimental 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. +Result vector or array `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 From 8d5e682e40bffc7680aa3be8cc38d839b7e7a86d Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Tue, 30 Apr 2024 10:06:19 +0200 Subject: [PATCH 19/29] Update src/stdlib_linalg_solve.fypp Co-authored-by: Jeremie Vandenplas --- src/stdlib_linalg_solve.fypp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stdlib_linalg_solve.fypp b/src/stdlib_linalg_solve.fypp index b44141ca5..4162ca257 100644 --- a/src/stdlib_linalg_solve.fypp +++ b/src/stdlib_linalg_solve.fypp @@ -111,7 +111,7 @@ submodule (stdlib_linalg) stdlib_linalg_solve end function stdlib_linalg_${ri}$_solve_${ndsuf}$ - ! Compute the solution to a real system of linear equations A * X = B (pure interface) + !> Compute the solution to a real system of linear equations A * X = B (pure interface) pure module function stdlib_linalg_${ri}$_pure_solve_${ndsuf}$(a,b) result(x) !> Input matrix a[n,n] ${rt}$, intent(in), target :: a(:,:) From 4458b88248c2c6fb72f91e0843a98f445d852bb2 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Tue, 30 Apr 2024 10:24:02 +0200 Subject: [PATCH 20/29] fix test --- test/linalg/test_linalg_solve.fypp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/test/linalg/test_linalg_solve.fypp b/test/linalg/test_linalg_solve.fypp index ff115965b..25234a40b 100644 --- a/test/linalg/test_linalg_solve.fypp +++ b/test/linalg/test_linalg_solve.fypp @@ -129,30 +129,30 @@ module test_linalg_solve end subroutine test_${ri}$_solve !> 2x2 Complex linear system - !> https://math.stackexchange.com/questions/1996540/solving-linear-equation-systems-with-complex-coefficients-and-variables subroutine test_2x2_${ri}$_solve(error) type(error_type), allocatable, intent(out) :: error type(linalg_state_type) :: state + + ${rt}$, parameter :: i = (0.0_${rk}$,1.0_${rk}$) ${rt}$ :: A(2,2), b(2), res(2), x(2) - integer(ilp) :: i ! Fill in linear system - A(1,:) = [(+1.0_${rk}$,+1.0_${rk}$),(-1.0_${rk}$,0.0_${rk}$)] - A(2,:) = [(+1.0_${rk}$,-1.0_${rk}$),(+1.0_${rk}$,1.0_${rk}$)] + A(1,:) = [ 1+2*i, 2-i] + A(2,:) = [ 2+i , i] - b = [(0.0_${rk}$,1.0_${rk}$),(1.0_${rk}$,0.0_${rk}$)] + b = [1,-1] ! Exact result - res = [(0.5_${rk}$,0.5_${rk}$),(0.0_${rk}$,0.0_${rk}$)] + res = [(-0.28_${rk}$,-0.04_${rk}$),(0.36_${rk}$,0.48_${rk}$)] x = solve(a,b,err=state) - + call check(error,state%ok(),state%print()) if (allocated(error)) return - call check(error, all(abs(x-res) Date: Tue, 30 Apr 2024 19:31:57 +0200 Subject: [PATCH 21/29] implement `subroutine` interface --- src/stdlib_linalg.fypp | 26 ++++++- src/stdlib_linalg_solve.fypp | 142 +++++++++++++++++------------------ 2 files changed, 96 insertions(+), 72 deletions(-) diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index 2839e81d5..4c964287d 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -21,6 +21,7 @@ module stdlib_linalg public :: diag public :: eye public :: solve + public :: solve_lu public :: trace public :: outer_product public :: kronecker_product @@ -264,7 +265,7 @@ module stdlib_linalg 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(:,:) + ${rt}$, intent(in) :: 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] @@ -275,6 +276,29 @@ module stdlib_linalg #:endfor end interface solve + interface solve_lu + #:for nd,ndsuf,nde in ALL_RHS + #:for rk,rt,ri in RC_KINDS_TYPES + #:if rk!="xdp" + pure module subroutine stdlib_linalg_${ri}$_solve_lu_${ndsuf}$(a,b,x,pivot,overwrite_a,err) + !> 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}$ + !> Result array/matrix x[n] or x[n,nrhs] + ${rt}$, intent(inout), contiguous, target :: x${nd}$ + !> [optional] Storage array for the diagonal pivot indices + integer(ilp), optional, intent(inout), target :: pivot(:) + !> [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 + end subroutine stdlib_linalg_${ri}$_solve_lu_${ndsuf}$ + #:endif + #:endfor + #:endfor + end interface solve_lu + interface det !! version: experimental !! diff --git a/src/stdlib_linalg_solve.fypp b/src/stdlib_linalg_solve.fypp index 4162ca257..8612d5052 100644 --- a/src/stdlib_linalg_solve.fypp +++ b/src/stdlib_linalg_solve.fypp @@ -55,27 +55,73 @@ submodule (stdlib_linalg) stdlib_linalg_solve type(linalg_state_type), intent(out) :: err !> Result array/matrix x[n] or x[n,nrhs] ${rt}$, allocatable, target :: x${nd}$ + + ! Initialize solution shape from the rhs array + allocate(x,mold=b) + + call stdlib_linalg_${ri}$_solve_lu_${ndsuf}$(a,b,x,overwrite_a=overwrite_a,err=err) + + end function stdlib_linalg_${ri}$_solve_${ndsuf}$ + + !> Compute the solution to a real system of linear equations A * X = B (pure interface) + pure module function stdlib_linalg_${ri}$_pure_solve_${ndsuf}$(a,b) result(x) + !> Input matrix a[n,n] + ${rt}$, intent(in) :: 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}$ + + ! Local variables + ${rt}$, allocatable :: amat(:,:) + + ! Copy `a` so it can be intent(in) + allocate(amat,source=a) + + ! Initialize solution shape from the rhs array + allocate(x,mold=b) + + call stdlib_linalg_${ri}$_solve_lu_${ndsuf}$(amat,b,x,overwrite_a=.true.) + + end function stdlib_linalg_${ri}$_pure_solve_${ndsuf}$ + + !> Compute the solution to a real system of linear equations A * X = B (pure interface) + pure module subroutine stdlib_linalg_${ri}$_solve_lu_${ndsuf}$(a,b,x,pivot,overwrite_a,err) + !> 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}$ + !> Result array/matrix x[n] or x[n,nrhs] + ${rt}$, intent(inout), contiguous, target :: x${nd}$ + !> [optional] Storage array for the diagonal pivot indices + integer(ilp), optional, intent(inout), target :: pivot(:) + !> [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 ! Local variables type(linalg_state_type) :: err0 - integer(ilp) :: lda,n,ldb,nrhs,info - integer(ilp), allocatable :: ipiv(:) + integer(ilp) :: lda,n,ldb,ldx,nrhsx,nrhs,info,npiv + integer(ilp), pointer :: ipiv(:) logical(lk) :: copy_a ${rt}$, pointer :: xmat(:,:),amat(:,:) ! 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 (any([lda,n,ldb]<1) .or. any([lda,ldb]/=n)) then - err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid sizes: a=',[lda,n], & - ', b=',[ldb,nrhs]) - allocate(x${nde}$) - call linalg_error_handling(err0,err) - return - end if + lda = size(a,1,kind=ilp) + n = size(a,2,kind=ilp) + ldb = size(b,1,kind=ilp) + nrhs = size(b ,kind=ilp)/ldb + ldx = size(x,1,kind=ilp) + nrhsx = size(x ,kind=ilp)/ldx + + ! 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) ! Can A be overwritten? By default, do not overwrite if (present(overwrite_a)) then @@ -84,8 +130,13 @@ submodule (stdlib_linalg) stdlib_linalg_solve copy_a = .true._lk endif - ! Pivot indices - allocate(ipiv(n)) + if (any([lda,n,ldb]<1) .or. any([lda,ldb,ldx]/=n) .or. nrhsx/=nrhs .or. npiv/=n) then + err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid sizes: a=',[lda,n], & + 'b=',[ldb,nrhs],' x=',[ldx,nrhsx], & + 'pivot=',n) + call linalg_error_handling(err0,err) + return + end if ! Initialize a matrix temporary if (copy_a) then @@ -95,7 +146,7 @@ submodule (stdlib_linalg) stdlib_linalg_solve endif ! Initialize solution with the rhs - allocate(x,source=b) + x = b xmat(1:n,1:nrhs) => x ! Solve system @@ -105,64 +156,13 @@ submodule (stdlib_linalg) stdlib_linalg_solve call handle_gesv_info(info,lda,n,nrhs,err0) if (copy_a) deallocate(amat) + if (.not.present(pivot)) deallocate(ipiv) ! Process output and return call linalg_error_handling(err0,err) - end function stdlib_linalg_${ri}$_solve_${ndsuf}$ - - !> Compute the solution to a real system of linear equations A * X = B (pure interface) - 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}$ - - ! Local variables - type(linalg_state_type) :: err0 - integer(ilp) :: lda,n,ldb,nrhs,info - integer(ilp), allocatable :: ipiv(:) - ${rt}$, pointer :: xmat(:,:) - ${rt}$, allocatable :: amat(:,:) - - ! 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 (any([lda,n,ldb]<1) .or. any([lda,ldb]/=n)) then - err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid sizes: a=',[lda,n], & - ', b=',[ldb,nrhs]) - allocate(x${nde}$) - call linalg_error_handling(err0) - return - end if - - ! Pivot indices - allocate(ipiv(n)) - - ! Initialize a matrix temporary - allocate(amat,source=a) - - ! 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 - call handle_gesv_info(info,lda,n,nrhs,err0) - - deallocate(amat) - - ! Process output and return - call linalg_error_handling(err0) - - end function stdlib_linalg_${ri}$_pure_solve_${ndsuf}$ + end subroutine stdlib_linalg_${ri}$_solve_lu_${ndsuf}$ + #:endif #:endfor #:endfor From 5e0620c3297f0d86e6517805c70786cf9061b52d Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 9 May 2024 08:24:56 +0200 Subject: [PATCH 22/29] specify full-rank --- doc/specs/stdlib_linalg.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index f5c7416c0..da98e2f83 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -636,7 +636,7 @@ Two ### Return value -Returns an array value that represents the solution to the linear system of equations. +For a full-rank matrix, 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. From c9f5f0c6641d9fee1acb70dcb893b66ebf1c097d Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 9 May 2024 08:34:35 +0200 Subject: [PATCH 23/29] document `solve_lu` --- doc/specs/stdlib_linalg.md | 56 ++++++++++++++++++++++++++++++++++++-- 1 file changed, 54 insertions(+), 2 deletions(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index da98e2f83..1c5a8cd4e 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -610,7 +610,9 @@ Experimental 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 or array `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. +Result vector or array `x` returns the exact solution to within numerical precision, provided that the matrix is not ill-conditioned. +An error is returned if the matrix is rank-deficient or singular to working precision. +The solver is based on LAPACK's `*GESV` backends. ### Syntax @@ -640,7 +642,57 @@ For a full-rank matrix, returns an array value that represents the solution to t 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`. +If `err` is not present, exceptions trigger an `error stop`. + +### Example + +```fortran +{!example/linalg/example_solve1.f90!} + +{!example/linalg/example_solve2.f90!} +``` + +## `solve_lu` - Solves a linear matrix equation or a linear system of equations (subroutine interface). + +### Status + +Experimental + +### Description + +This subroutine 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 or array `x` returns the exact solution to within numerical precision, provided that the matrix is not ill-conditioned. +An error is returned if the matrix is rank-deficient or singular to working precision. +The solver is based on LAPACK's `*GESV` backends. + +### Syntax + +`Pure` interface: + +`x = ` [[stdlib_linalg(module):solve_lu(interface)]] `(a, b, x, pivot [, 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 or rank-2 array of the same kind as `a`, containing the right-hand-side vector(s). It is an `intent(in)` argument. + +`x`: Shall be a rank-1 or rank-2 array of the same kind and size as `b`, that returns the solution(s) to the system. It is an `intent(inout)` argument, and must have the `contiguous` property. + +`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. + +### Return value + +For a full-rank matrix, 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. +If `err` is not present, exceptions trigger an `error stop`. ### Example From e75bb2fff71bb7a967fb17f62ddf26cd8f2bdc35 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 9 May 2024 08:34:41 +0200 Subject: [PATCH 24/29] add `solve_lu` test --- example/linalg/CMakeLists.txt | 1 + example/linalg/example_solve3.f90 | 30 ++++++++++++++++++++++++++++++ 2 files changed, 31 insertions(+) create mode 100644 example/linalg/example_solve3.f90 diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt index 1c358e147..60ca821e8 100644 --- a/example/linalg/CMakeLists.txt +++ b/example/linalg/CMakeLists.txt @@ -20,5 +20,6 @@ ADD_EXAMPLE(blas_gemv) ADD_EXAMPLE(lapack_getrf) ADD_EXAMPLE(solve1) ADD_EXAMPLE(solve2) +ADD_EXAMPLE(solve3) ADD_EXAMPLE(determinant) ADD_EXAMPLE(determinant2) diff --git a/example/linalg/example_solve3.f90 b/example/linalg/example_solve3.f90 new file mode 100644 index 000000000..8ee15c091 --- /dev/null +++ b/example/linalg/example_solve3.f90 @@ -0,0 +1,30 @@ +program example_solve3 + use stdlib_linalg_constants, only: sp + use stdlib_linalg, only: solve_lu, linalg_state_type + implicit none + + integer(ilp) :: test + 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)] + + ! Pre-allocate x + allocate(x,source=b) + + ! Call system many times avoiding reallocation + do test=1,100 + call solve_lu(A,b,x) + print "(i2,'-th solution: ',*(1x,f12.6))", test,x + end do + +end program example_solve3 + From 449d0a25c208c5b2a99188a61fcd119e9c86fb9e Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 9 May 2024 08:46:10 +0200 Subject: [PATCH 25/29] add pivot --- doc/specs/stdlib_linalg.md | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 1c5a8cd4e..66a193798 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -668,9 +668,13 @@ The solver is based on LAPACK's `*GESV` backends. ### Syntax -`Pure` interface: +Simple (`Pure`) interface: + +`x = ` [[stdlib_linalg(module):solve_lu(interface)]] `(a, b, x)` -`x = ` [[stdlib_linalg(module):solve_lu(interface)]] `(a, b, x, pivot [, overwrite_a, err])` +Expert (`Pure`) interface: + +`x = ` [[stdlib_linalg(module):solve_lu(interface)]] `(a, b, x [, pivot, overwrite_a, err])` ### Arguments @@ -682,6 +686,8 @@ Two `x`: Shall be a rank-1 or rank-2 array of the same kind and size as `b`, that returns the solution(s) to the system. It is an `intent(inout)` argument, and must have the `contiguous` property. +`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. + `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. From b288520c406939724f3b081b547b0e25ee87c528 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 9 May 2024 08:49:32 +0200 Subject: [PATCH 26/29] cleanup subroutine example; add preallocated pivot --- example/linalg/example_solve3.f90 | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/example/linalg/example_solve3.f90 b/example/linalg/example_solve3.f90 index 8ee15c091..ae92ea5da 100644 --- a/example/linalg/example_solve3.f90 +++ b/example/linalg/example_solve3.f90 @@ -1,9 +1,10 @@ program example_solve3 - use stdlib_linalg_constants, only: sp + use stdlib_linalg_constants, only: sp,ilp use stdlib_linalg, only: solve_lu, linalg_state_type implicit none integer(ilp) :: test + integer(ilp), allocatable :: pivot(:) complex(sp), allocatable :: A(:,:),b(:),x(:) ! Solve a system of 3 complex linear equations: @@ -15,15 +16,16 @@ program example_solve3 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)] ! Pre-allocate x - allocate(x,source=b) + allocate(b(size(A,2)),pivot(size(A,2))) + allocate(x,mold=b) ! Call system many times avoiding reallocation do test=1,100 - call solve_lu(A,b,x) - print "(i2,'-th solution: ',*(1x,f12.6))", test,x + b = test*[(5.0,-1.0),(0.0,1.0),(1.0,0.0)] + call solve_lu(A,b,x,pivot) + print "(i3,'-th solution: ',*(1x,f12.6))", test,x end do end program example_solve3 From 7aab8446848c40af22d0d16e1daf18b8e1bf1a60 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 9 May 2024 08:53:16 +0200 Subject: [PATCH 27/29] document `solve_lu` interface --- src/stdlib_linalg.fypp | 24 +++++++++++++++++++++++- 1 file changed, 23 insertions(+), 1 deletion(-) diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index 4c964287d..9f1ee0c85 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -232,7 +232,7 @@ module stdlib_linalg !! 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)) + !! ([Specification](../page/specs/stdlib_linalg.html#solve-solves-a-linear-matrix-equation-or-a-linear-system-of-equations)) !! !!### Summary !! Interface for solving a linear system arising from a general matrix. @@ -276,7 +276,29 @@ module stdlib_linalg #:endfor end interface solve + ! Solve linear system Ax = b using LU decomposition (subroutine interface). interface solve_lu + !! 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-lu-solves-a-linear-matrix-equation-or-a-linear-system-of-equations-subroutine-interface)) + !! + !!### Summary + !! Subroutine interface for solving a linear system using LU decomposition. + !! + !!### Description + !! + !! This interface provides methods for computing the solution of a linear matrix system using + !! a subroutine. Supported data types include `real` and `complex`. No assumption is made on the matrix + !! structure. Preallocated space for the solution vector `x` is user-provided, and it 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 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" From a05809eba2925e949527ff9cb19f737a05b592b2 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 9 May 2024 10:59:21 +0200 Subject: [PATCH 28/29] typo --- doc/specs/stdlib_linalg.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 8d877ad01..c414c071a 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -669,11 +669,11 @@ The solver is based on LAPACK's `*GESV` backends. Simple (`Pure`) interface: -`x = ` [[stdlib_linalg(module):solve_lu(interface)]] `(a, b, x)` +`call ` [[stdlib_linalg(module):solve_lu(interface)]] `(a, b, x)` Expert (`Pure`) interface: -`x = ` [[stdlib_linalg(module):solve_lu(interface)]] `(a, b, x [, pivot, overwrite_a, err])` +`call ` [[stdlib_linalg(module):solve_lu(interface)]] `(a, b, x [, pivot, overwrite_a, err])` ### Arguments From 5832df5f51b010b984299d28264b3af5216cc7b9 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 9 May 2024 11:21:28 +0200 Subject: [PATCH 29/29] avoid 128-bit random numbers --- test/linalg/test_linalg_lstsq.fypp | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/test/linalg/test_linalg_lstsq.fypp b/test/linalg/test_linalg_lstsq.fypp index 29cd07519..9aba94871 100644 --- a/test/linalg/test_linalg_lstsq.fypp +++ b/test/linalg/test_linalg_lstsq.fypp @@ -70,14 +70,17 @@ module test_linalg_least_squares type(linalg_state_type) :: state integer(ilp), parameter :: n = 12, m = 3 + real :: Arnd(n,m),xrnd(m) ${rt}$ :: xsol(m),x(m),y(n),A(n,m) ! Random coefficient matrix and solution - call random_number(A) - call random_number(xsol) + call random_number(Arnd) + call random_number(xrnd) ! Compute rhs - y = matmul(A,xsol) + A = real(Arnd,${rk}$) + xsol = real(xrnd,${rk}$) + y = matmul(A,xsol) ! Find polynomial x = lstsq(A,y,err=state)