From 892bf1efbb53c9aceda4547895fbb11084eb420e Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Sun, 21 Apr 2024 15:34:51 +0200 Subject: [PATCH 01/24] 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 465d230c60a7ce489602fb76d91b934b9d21904d Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Sun, 21 Apr 2024 15:43:21 +0200 Subject: [PATCH 02/24] base implementation --- src/CMakeLists.txt | 1 + src/stdlib_linalg_least_squares.fypp | 210 +++++++++++++++++++++++++++ 2 files changed, 211 insertions(+) create mode 100644 src/stdlib_linalg_least_squares.fypp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 165259db3..7c3294943 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -21,6 +21,7 @@ set(fppFiles stdlib_kinds.fypp stdlib_linalg.fypp stdlib_linalg_diag.fypp + stdlib_linalg_least_squares.fypp stdlib_linalg_outer_product.fypp stdlib_linalg_kronecker.fypp stdlib_linalg_cross_product.fypp diff --git a/src/stdlib_linalg_least_squares.fypp b/src/stdlib_linalg_least_squares.fypp new file mode 100644 index 000000000..94c6f8bd1 --- /dev/null +++ b/src/stdlib_linalg_least_squares.fypp @@ -0,0 +1,210 @@ +#: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_least_squares + use stdlib_linalg_constants + use stdlib_linalg_lapack, only: gelsd, stdlib_ilaenv + use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, & + LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR + implicit none(type,external) + private + + !> Compute a least squares solution to system Ax=b, i.e. such that the 2-norm abs(b-Ax) is minimized. + public :: lstsq + + ! NumPy: lstsq(a, b, rcond='warn') + ! Scipy: lstsq(a, b, cond=None, overwrite_a=False, overwrite_b=False, check_finite=True, lapack_driver=None) + ! IMSL: Result = IMSL_QRSOL(B, [A] [, AUXQR] [, BASIS] [, /DOUBLE] [, QR] [, PIVOT] [, RESIDUAL] [, TOLERANCE]) + + interface lstsq + #:for nd,ndsuf,nde in ALL_RHS + #:for rk,rt,ri in RC_KINDS_TYPES + module procedure stdlib_linalg_${ri}$_lstsq_${ndsuf}$ + #:endfor + #:endfor + end interface lstsq + + + contains + + #:for rk,rt,ri in RC_KINDS_TYPES + ! Workspace needed by gesv + subroutine ${ri}$gesv_space(m,n,nrhs,lrwork,liwork,lcwork) + integer(ilp), intent(in) :: m,n,nrhs + integer(ilp), intent(out) :: lrwork,liwork,lcwork + + integer(ilp) :: smlsiz,mnmin,nlvl + + mnmin = min(m,n) + + ! Maximum size of the subproblems at the bottom of the computation (~25) + smlsiz = stdlib_ilaenv(9,'${ri}$gelsd',' ',0,0,0,0) + + ! The exact minimum amount of workspace needed depends on M, N and NRHS. As long as LWORK is at least + nlvl = max(0, ilog2(mnmin/(smlsiz+1))+1) + + ! Real space + #:if rt.startswith('complex') + lrwork = 10*mnmin+2*mnmin*smlsiz+8*mnmin*nlvl+3*smlsiz*nrhs+max((smlsiz+1)**2,n*(1+nrhs)+2*nrhs) + #:else + lrwork = 12*mnmin+2*mnmin*smlsiz+8*mnmin*nlvl+mnmin*nrhs+(smlsiz+1)**2 + #:endif + lrwork = max(1,lrwork) + + ! Complex space + lcwork = 2*mnmin + nrhs*mnmin + + ! Integer space + liwork = max(1, 3*mnmin*nlvl+11*mnmin) + + ! For good performance, the workspace should generally be larger. + lrwork = ceiling(1.25*lrwork,kind=ilp) + lcwork = ceiling(1.25*lcwork,kind=ilp) + liwork = ceiling(1.25*liwork,kind=ilp) + + end subroutine ${ri}$gesv_space + + #:endfor + + #:for nd,ndsuf,nde in ALL_RHS + #:for rk,rt,ri in RC_KINDS_TYPES + + ! Compute the least-squares solution to a real system of linear equations Ax = B + function stdlib_linalg_${ri}$_lstsq_${ndsuf}$(a,b,cond,overwrite_a,rank,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] cutoff for rank evaluation: singular values s(i)<=cond*maxval(s) are considered 0. + real(${rk}$), optional, intent(in) :: cond + !> [optional] Can A,b data be overwritten and destroyed? + logical(lk), optional, intent(in) :: overwrite_a + !> [optional] Return rank of A + integer(ilp), optional, intent(out) :: rank + !> [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) :: m,n,lda,ldb,nrhs,info,mnmin,mnmax,arank,lrwork,liwork,lcwork + integer(ilp), allocatable :: iwork(:) + logical(lk) :: copy_a + real(${rk}$) :: acond,rcond + real(${rk}$), allocatable :: singular(:),rwork(:) + ${rt}$, pointer :: xmat(:,:),amat(:,:) + ${rt}$, allocatable :: cwork(:) + character(*), parameter :: this = 'lstsq' + + !> Problem sizes + m = size(a,1,kind=ilp) + lda = size(a,1,kind=ilp) + n = size(a,2,kind=ilp) + ldb = size(b,1,kind=ilp) + nrhs = size(b ,kind=ilp)/ldb + mnmin = min(m,n) + mnmax = max(m,n) + + if (lda<1 .or. n<1 .or. ldb<1 .or. ldb/=m) then + err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid sizes: a=',[lda,n], & + 'b=',[ldb,nrhs]) + allocate(x${nde}$) + arank = 0 + 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 + + ! 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 + + ! Singular values array (in degreasing order) + allocate(singular(mnmin)) + + ! rcond is used to determine the effective rank of A. + ! Singular values S(i) <= RCOND*maxval(S) are treated as zero. + ! Use same default value as NumPy + if (present(cond)) then + rcond = cond + else + rcond = epsilon(0.0_${rk}$)*mnmax + endif + if (rcond<0) rcond = epsilon(0.0_${rk}$)*mnmax + + ! Allocate working space + call ${ri}$gesv_space(m,n,nrhs,lrwork,liwork,lcwork) + #:if rt.startswith('complex') + allocate(rwork(lrwork),cwork(lcwork),iwork(liwork)) + #:else + allocate(rwork(lrwork),iwork(liwork)) + #:endif + + ! Solve system using singular value decomposition + #:if rt.startswith('complex') + call gelsd(m,n,nrhs,amat,lda,xmat,ldb,singular,rcond,arank,cwork,lrwork,rwork,iwork,info) + #:else + call gelsd(m,n,nrhs,amat,lda,xmat,ldb,singular,rcond,arank,rwork,lrwork,iwork,info) + #:endif + + ! The condition number of A in the 2-norm = S(1)/S(min(m,n)). + acond = singular(1)/singular(mnmin) + + ! Process output + select case (info) + case (0) + ! Success + case (:-1) + err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid problem size a=[',lda,',',n,'], b[',ldb,',',nrhs,']') + case (1:) + err0 = linalg_state_type(this,LINALG_ERROR,'SVD did not converge.') + 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) + if (present(rank)) rank = arank + + end function stdlib_linalg_${ri}$_lstsq_${ndsuf}$ + + #:endfor + #:endfor + + ! Simple integer log2 implementation + elemental integer(ilp) function ilog2(x) + integer(ilp), intent(in) :: x + + integer(ilp) :: remndr + + if (x>0) then + remndr = x + ilog2 = -1_ilp + do while (remndr>0) + ilog2 = ilog2 + 1_ilp + remndr = shiftr(remndr,1) + end do + else + ilog2 = -huge(0_ilp) + endif + end function ilog2 + +end module stdlib_linalg_least_squares From 74756802efe03093735abe8e34f340e46788fc31 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Sun, 21 Apr 2024 16:15:41 +0200 Subject: [PATCH 03/24] submodule --- src/stdlib_linalg.fypp | 37 +++++++++++++++++++++++--- src/stdlib_linalg_least_squares.fypp | 39 +++++++++------------------- 2 files changed, 46 insertions(+), 30 deletions(-) diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index 3ed905c56..ce214e27f 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 :: lstsq public :: trace public :: outer_product public :: kronecker_product @@ -214,6 +221,30 @@ module stdlib_linalg #:endfor end interface is_hessenberg + ! Least squares solution to system Ax=b, i.e. such that the 2-norm abs(b-Ax) is minimized. + interface lstsq + #:for nd,ndsuf,nde in ALL_RHS + #:for rk,rt,ri in RC_KINDS_TYPES + module function stdlib_linalg_${ri}$_lstsq_${ndsuf}$(a,b,cond,overwrite_a,rank,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] cutoff for rank evaluation: singular values s(i)<=cond*maxval(s) are considered 0. + real(${rk}$), optional, intent(in) :: cond + !> [optional] Can A,b data be overwritten and destroyed? + logical(lk), optional, intent(in) :: overwrite_a + !> [optional] Return rank of A + integer(ilp), optional, intent(out) :: rank + !> [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}$ + end function stdlib_linalg_${ri}$_lstsq_${ndsuf}$ + #:endfor + #:endfor + end interface lstsq + contains diff --git a/src/stdlib_linalg_least_squares.fypp b/src/stdlib_linalg_least_squares.fypp index 94c6f8bd1..b99cc8097 100644 --- a/src/stdlib_linalg_least_squares.fypp +++ b/src/stdlib_linalg_least_squares.fypp @@ -4,35 +4,21 @@ #: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_least_squares +submodule (stdlib_linalg) stdlib_linalg_least_squares +!! Least-squares solution to Ax=b use stdlib_linalg_constants use stdlib_linalg_lapack, only: gelsd, stdlib_ilaenv use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, & LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR implicit none(type,external) - private - - !> Compute a least squares solution to system Ax=b, i.e. such that the 2-norm abs(b-Ax) is minimized. - public :: lstsq - - ! NumPy: lstsq(a, b, rcond='warn') - ! Scipy: lstsq(a, b, cond=None, overwrite_a=False, overwrite_b=False, check_finite=True, lapack_driver=None) - ! IMSL: Result = IMSL_QRSOL(B, [A] [, AUXQR] [, BASIS] [, /DOUBLE] [, QR] [, PIVOT] [, RESIDUAL] [, TOLERANCE]) - - interface lstsq - #:for nd,ndsuf,nde in ALL_RHS - #:for rk,rt,ri in RC_KINDS_TYPES - module procedure stdlib_linalg_${ri}$_lstsq_${ndsuf}$ - #:endfor - #:endfor - end interface lstsq - + + character(*), parameter :: this = 'lstsq' contains #:for rk,rt,ri in RC_KINDS_TYPES ! Workspace needed by gesv - subroutine ${ri}$gesv_space(m,n,nrhs,lrwork,liwork,lcwork) + elemental subroutine ${ri}$gesv_space(m,n,nrhs,lrwork,liwork,lcwork) integer(ilp), intent(in) :: m,n,nrhs integer(ilp), intent(out) :: lrwork,liwork,lcwork @@ -73,11 +59,11 @@ module stdlib_linalg_least_squares #:for rk,rt,ri in RC_KINDS_TYPES ! Compute the least-squares solution to a real system of linear equations Ax = B - function stdlib_linalg_${ri}$_lstsq_${ndsuf}$(a,b,cond,overwrite_a,rank,err) result(x) + module function stdlib_linalg_${ri}$_lstsq_${ndsuf}$(a,b,cond,overwrite_a,rank,err) result(x) !> Input matrix a[n,n] - ${rt}$, intent(inout), target :: a(:,:) + ${rt}$, intent(inout), target :: a(:,:) !> Right hand side vector or array, b[n] or b[n,nrhs] - ${rt}$, intent(in) :: b${nd}$ + ${rt}$, intent(in) :: b${nd}$ !> [optional] cutoff for rank evaluation: singular values s(i)<=cond*maxval(s) are considered 0. real(${rk}$), optional, intent(in) :: cond !> [optional] Can A,b data be overwritten and destroyed? @@ -88,8 +74,8 @@ module stdlib_linalg_least_squares type(linalg_state_type), optional, intent(out) :: err !> 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) :: m,n,lda,ldb,nrhs,info,mnmin,mnmax,arank,lrwork,liwork,lcwork integer(ilp), allocatable :: iwork(:) @@ -98,9 +84,8 @@ module stdlib_linalg_least_squares real(${rk}$), allocatable :: singular(:),rwork(:) ${rt}$, pointer :: xmat(:,:),amat(:,:) ${rt}$, allocatable :: cwork(:) - character(*), parameter :: this = 'lstsq' - !> Problem sizes + ! Problem sizes m = size(a,1,kind=ilp) lda = size(a,1,kind=ilp) n = size(a,2,kind=ilp) @@ -207,4 +192,4 @@ module stdlib_linalg_least_squares endif end function ilog2 -end module stdlib_linalg_least_squares +end submodule stdlib_linalg_least_squares From 285aed23c8cbe231f14900840158ba820b3c514b Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Sun, 21 Apr 2024 16:16:38 +0200 Subject: [PATCH 04/24] remove `xdp` --- src/stdlib_linalg.fypp | 2 ++ src/stdlib_linalg_least_squares.fypp | 4 ++++ 2 files changed, 6 insertions(+) diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index ce214e27f..fd8e93b61 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -225,6 +225,7 @@ module stdlib_linalg interface lstsq #:for nd,ndsuf,nde in ALL_RHS #:for rk,rt,ri in RC_KINDS_TYPES + #:if rk!="xdp" module function stdlib_linalg_${ri}$_lstsq_${ndsuf}$(a,b,cond,overwrite_a,rank,err) result(x) !> Input matrix a[n,n] ${rt}$, intent(inout), target :: a(:,:) @@ -241,6 +242,7 @@ module stdlib_linalg !> Result array/matrix x[n] or x[n,nrhs] ${rt}$, allocatable, target :: x${nd}$ end function stdlib_linalg_${ri}$_lstsq_${ndsuf}$ + #:endif #:endfor #:endfor end interface lstsq diff --git a/src/stdlib_linalg_least_squares.fypp b/src/stdlib_linalg_least_squares.fypp index b99cc8097..48ddbbd54 100644 --- a/src/stdlib_linalg_least_squares.fypp +++ b/src/stdlib_linalg_least_squares.fypp @@ -17,6 +17,7 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares contains #:for rk,rt,ri in RC_KINDS_TYPES + #:if rk!="xdp" ! Workspace needed by gesv elemental subroutine ${ri}$gesv_space(m,n,nrhs,lrwork,liwork,lcwork) integer(ilp), intent(in) :: m,n,nrhs @@ -53,10 +54,12 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares end subroutine ${ri}$gesv_space + #:endif #:endfor #:for nd,ndsuf,nde in ALL_RHS #:for rk,rt,ri in RC_KINDS_TYPES + #:if rk!="xdp" ! Compute the least-squares solution to a real system of linear equations Ax = B module function stdlib_linalg_${ri}$_lstsq_${ndsuf}$(a,b,cond,overwrite_a,rank,err) result(x) @@ -171,6 +174,7 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares end function stdlib_linalg_${ri}$_lstsq_${ndsuf}$ + #:endif #:endfor #:endfor From f7251604660e41b083eaaae2e9ada79546d61df9 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 22 Apr 2024 15:03:52 +0200 Subject: [PATCH 05/24] add test --- test/linalg/CMakeLists.txt | 2 + test/linalg/test_linalg_lstsq.fypp | 96 ++++++++++++++++++++++++++++++ 2 files changed, 98 insertions(+) create mode 100644 test/linalg/test_linalg_lstsq.fypp diff --git a/test/linalg/CMakeLists.txt b/test/linalg/CMakeLists.txt index 3d590a9d2..b7bf47f00 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_lstsq.fypp" "test_linalg_matrix_property_checks.fypp" ) fypp_f90("${fyppFlags}" "${fppFiles}" outFiles) ADDTEST(linalg) ADDTEST(linalg_matrix_property_checks) +ADDTEST(linalg_lstsq) ADDTEST(blas_lapack) diff --git a/test/linalg/test_linalg_lstsq.fypp b/test/linalg/test_linalg_lstsq.fypp new file mode 100644 index 000000000..435d71129 --- /dev/null +++ b/test/linalg/test_linalg_lstsq.fypp @@ -0,0 +1,96 @@ +#:include "common.fypp" +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES +! Test least squares solver +module test_linalg_least_squares + use testdrive, only: error_type, check, new_unittest, unittest_type + use stdlib_linalg_constants + use stdlib_linalg, only: lstsq + use stdlib_linalg_state, only: linalg_state_type + + implicit none (type,external) + private + + public :: test_least_squares + + contains + + + !> Solve sample least squares problems + subroutine test_least_squares(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("lease_squares_${ri}$",test_lstsq_one_${ri}$)] + #:endif + #:endfor + + end subroutine test_least_squares + + !> Simple polynomial fit + #:for rk,rt,ri in REAL_KINDS_TYPES + #:if rk!="xdp" + subroutine test_lstsq_one_${ri}$(error) + type(error_type), allocatable, intent(out) :: error + + type(linalg_state_type) :: state + integer(ilp) :: rank + + !> Example scattered data + ${rt}$, parameter :: x(*) = real([1.0, 2.5, 3.5, 4.0, 5.0, 7.0, 8.5], ${rk}$) + ${rt}$, parameter :: y(*) = real([0.3, 1.1, 1.5, 2.0, 3.2, 6.6, 8.6], ${rk}$) + ${rt}$, parameter :: ab(*) = real([0.20925829, 0.12013861], ${rk}$) + + ${rt}$ :: M(size(x),2),p(2) + + ! Coefficient matrix for polynomial y = a + b*x**2 + M(:,1) = x**0 + M(:,2) = x**2 + + ! Find polynomial + p = lstsq(M,y,rank=rank,err=state) + + call check(error,state%ok(),state%print()) + if (allocated(error)) return + + call check(error, all(abs(p-ab)<1.0e-6_${rk}$), 'data converged') + if (allocated(error)) return + + call check(error, rank==2, 'matrix rank == 2') + if (allocated(error)) return + + end subroutine test_lstsq_one_${ri}$ + + #:endif + #:endfor + +end module test_linalg_least_squares + +program test_lstsq + use, intrinsic :: iso_fortran_env, only : error_unit + use testdrive, only : run_testsuite, new_testsuite, testsuite_type + use test_linalg_least_squares, only : test_least_squares + implicit none + integer :: stat, is + type(testsuite_type), allocatable :: testsuites(:) + character(len=*), parameter :: fmt = '("#", *(1x, a))' + + stat = 0 + + testsuites = [ & + new_testsuite("linalg_least_squares", test_least_squares) & + ] + + do is = 1, size(testsuites) + write(error_unit, fmt) "Testing:", testsuites(is)%name + call run_testsuite(testsuites(is)%collect, error_unit, stat) + end do + + if (stat > 0) then + write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!" + error stop + end if +end program test_lstsq From 71a1fe5d1c58b90fcbc2d7255e6c78e0119e87d1 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 22 Apr 2024 16:44:16 +0200 Subject: [PATCH 06/24] add example --- example/linalg/CMakeLists.txt | 1 + example/linalg/example_lstsq.f90 | 27 +++++++++++++++++++++++++++ 2 files changed, 28 insertions(+) create mode 100644 example/linalg/example_lstsq.f90 diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt index 1a5875502..3ace1b980 100644 --- a/example/linalg/CMakeLists.txt +++ b/example/linalg/CMakeLists.txt @@ -18,3 +18,4 @@ ADD_EXAMPLE(state1) ADD_EXAMPLE(state2) ADD_EXAMPLE(blas_gemv) ADD_EXAMPLE(lapack_getrf) +ADD_EXAMPLE(lstsq) diff --git a/example/linalg/example_lstsq.f90 b/example/linalg/example_lstsq.f90 new file mode 100644 index 000000000..2b36b0577 --- /dev/null +++ b/example/linalg/example_lstsq.f90 @@ -0,0 +1,27 @@ +program example_lstsq + use stdlib_linalg_constants, only: dp + use stdlib_linalg, only: lstsq + use stdlib_linalg_state, only: linalg_state_type + implicit none + type(linalg_state_type) :: err + + integer, allocatable :: x(:),y(:) + real(dp), allocatable :: A(:,:),b(:),coef(:) + + ! Data set + x = [1, 2, 2] + y = [5, 13, 25] + + ! Fit three points using a parabola, least squares method + ! A = [1 x x**2] + A = reshape([[1,1,1],x,x**2],[3,3]) + b = y + + ! Get coefficients of y = coef(1) + x*coef(2) + x^2*coef(3) + coef = lstsq(A,b) + + print *, 'parabola: ',coef + ! parabola: -0.42857142857141695 1.1428571428571503 4.2857142857142811 + + +end program example_lstsq From 2b2ade91ebce1486744d98578e4a230cfbb3268b Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 22 Apr 2024 17:04:24 +0200 Subject: [PATCH 07/24] add specs --- doc/specs/stdlib_linalg.md | 44 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 44 insertions(+) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index defe30758..033970efd 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -599,3 +599,47 @@ Specifically, upper Hessenberg matrices satisfy `a_ij = 0` when `j < i-1`, and l ```fortran {!example/linalg/example_is_hessenberg.f90!} ``` + +## `lstsq` - Computes the least squares solution to a linear matrix equation. + +### Status + +Experimental + +### Description + +This function computes the least-squares solution to a linear matrix equation \( A \cdot x = b \). + +Result vector `x` returns the approximate solution that minimizes the 2-norm \( || A \cdot x - b ||_2 \), i.e., it contains the least-squares solution to the problem. Matrix `A` may be full-rank, over-determined, or under-determined. The solver is based on LAPACK's `*GELSD` backends. + +### Syntax + +`x = ` [[stdlib_linalg(module):lstsq(interface)]] `(a, b, [, cond, overwrite_a, rank, err])` + +### Arguments + +`a`: Shall be a rank-2 square array containing the coefficient matrix. It is an `intent(inout)` argument. + +`b`: Shall be a rank-1 array containing the right-hand-side vector. It is an `intent(in)` argument. + +`cond` (optional): Singular value cut-off threshold for rank evaluation: `s_i >= cond*maxval(s), i=1:rank`. Shall be a scalar, `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. + +`rank` (optional): Shall be an `integer` scalar value, that contains the rank of input matrix `A`. This is an `intent(out)` argument. + +`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. + +### Return value + +Returns an array value that represents the solution to the least squares system. + +Raises `LINALG_ERROR` if the underlying SVD process did not converge. +Raises `LINALG_VALUE_ERROR` if the matrix and rhs vectors have invalid/incompatible sizes. +Exceptions trigger an `error stop`. + +### Example + +```fortran +{!example/linalg/example_lstsq.f90!} +``` From 905692c10bd2ef5a94cdf637715b923edb617a5b Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 22 Apr 2024 17:14:04 +0200 Subject: [PATCH 08/24] document interface --- src/stdlib_linalg.fypp | 16 ++++++++++++++++ src/stdlib_linalg_least_squares.fypp | 18 +++++++++++++++++- 2 files changed, 33 insertions(+), 1 deletion(-) diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index fd8e93b61..24f2b2e66 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -223,6 +223,22 @@ module stdlib_linalg ! Least squares solution to system Ax=b, i.e. such that the 2-norm abs(b-Ax) is minimized. interface lstsq + !! version: experimental + !! + !! Computes the squares solution to system \( A \cdot x = b \). + !! ([Specification](../page/specs/stdlib_linalg.html#det-computes-the-determinant-of-a-square-matrix)) + !! + !!### Summary + !! Interface for computing least squares, i.e. the 2-norm \( || (b-A \cdot x ||_2 \) minimizing solution. + !! + !!### Description + !! + !! This interface provides methods for computing the least squares of a linear matrix system. + !! Supported data types include `real` and `complex`. + !! + !!@note The solution is based on LAPACK's singular value decomposition `*GELSD` methods. + !!@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" diff --git a/src/stdlib_linalg_least_squares.fypp b/src/stdlib_linalg_least_squares.fypp index 48ddbbd54..01bccb380 100644 --- a/src/stdlib_linalg_least_squares.fypp +++ b/src/stdlib_linalg_least_squares.fypp @@ -61,8 +61,24 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares #:for rk,rt,ri in RC_KINDS_TYPES #:if rk!="xdp" - ! Compute the least-squares solution to a real system of linear equations Ax = B + ! Compute the least-squares solution to a real system of linear equations Ax = b module function stdlib_linalg_${ri}$_lstsq_${ndsuf}$(a,b,cond,overwrite_a,rank,err) result(x) + !!### Summary + !! Compute least-squares solution to a real system of linear equations \( Ax = b \) + !! + !!### Description + !! + !! This function computes the least-squares solution of a linear matrix problem. + !! + !! param: a Input matrix of size [m,n]. + !! param: b Right-hand-side vector of size [n] or matrix of size [n,nrhs]. + !! param: cond [optional] Real input threshold indicating that singular values `s_i <= cond*maxval(s)` + !! do not contribute to the matrix rank. + !! param: overwrite_a [optional] Flag indicating if the input matrix can be overwritten. + !! param: rank [optional] integer flag returning matrix rank. + !! param: err [optional] State return flag. + !! return: x Solution vector of size [n] or solution matrix of size [n,nrhs]. + !! !> Input matrix a[n,n] ${rt}$, intent(inout), target :: a(:,:) !> Right hand side vector or array, b[n] or b[n,nrhs] From cd9b3e64b34a7788e917523ef475dd9b749c17ae Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 25 Apr 2024 16:20:53 +0200 Subject: [PATCH 09/24] clarify space; remove `goto`; fix allocation leak --- src/stdlib_linalg_least_squares.fypp | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/stdlib_linalg_least_squares.fypp b/src/stdlib_linalg_least_squares.fypp index 01bccb380..481043c86 100644 --- a/src/stdlib_linalg_least_squares.fypp +++ b/src/stdlib_linalg_least_squares.fypp @@ -30,7 +30,7 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares ! Maximum size of the subproblems at the bottom of the computation (~25) smlsiz = stdlib_ilaenv(9,'${ri}$gelsd',' ',0,0,0,0) - ! The exact minimum amount of workspace needed depends on M, N and NRHS. As long as LWORK is at least + ! The exact minimum amount of workspace needed depends on M, N and NRHS. nlvl = max(0, ilog2(mnmin/(smlsiz+1))+1) ! Real space @@ -47,7 +47,8 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares ! Integer space liwork = max(1, 3*mnmin*nlvl+11*mnmin) - ! For good performance, the workspace should generally be larger. + ! For good performance, the workspace should generally be larger. + ! Allocate 25% more space than strictly needed. lrwork = ceiling(1.25*lrwork,kind=ilp) lcwork = ceiling(1.25*lcwork,kind=ilp) liwork = ceiling(1.25*liwork,kind=ilp) @@ -117,8 +118,8 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid sizes: a=',[lda,n], & 'b=',[ldb,nrhs]) allocate(x${nde}$) - arank = 0 - goto 1 + call linalg_error_handling(err0,err) + if (present(rank)) rank = 0 end if ! Can A be overwritten? By default, do not overwrite @@ -182,10 +183,10 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares err0 = linalg_state_type(this,LINALG_INTERNAL_ERROR,'catastrophic error') end select - 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) if (present(rank)) rank = arank end function stdlib_linalg_${ri}$_lstsq_${ndsuf}$ From 56bf7f7a7ba5eb7751d1d60463fb7f5ab9e50cbb Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 25 Apr 2024 16:22:47 +0200 Subject: [PATCH 10/24] typo --- src/stdlib_linalg_least_squares.fypp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stdlib_linalg_least_squares.fypp b/src/stdlib_linalg_least_squares.fypp index 481043c86..4169d7802 100644 --- a/src/stdlib_linalg_least_squares.fypp +++ b/src/stdlib_linalg_least_squares.fypp @@ -140,7 +140,7 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares allocate(x,source=b) xmat(1:n,1:nrhs) => x - ! Singular values array (in degreasing order) + ! Singular values array (in decreasing order) allocate(singular(mnmin)) ! rcond is used to determine the effective rank of A. From c587e1109ff3fc50da781b5579dde1913482c4ec Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 25 Apr 2024 19:48:03 +0200 Subject: [PATCH 11/24] doc typo --- 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 24f2b2e66..1015afca6 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -226,7 +226,7 @@ module stdlib_linalg !! version: experimental !! !! Computes the squares solution to system \( A \cdot x = b \). - !! ([Specification](../page/specs/stdlib_linalg.html#det-computes-the-determinant-of-a-square-matrix)) + !! ([Specification](../page/specs/stdlib_linalg.html#lstsq-computes-the-least-squares-solution-to-a-linear-matrix-equation)) !! !!### Summary !! Interface for computing least squares, i.e. the 2-norm \( || (b-A \cdot x ||_2 \) minimizing solution. From c1d0100f42babb1823d7a2ae6df10fd24becc514 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 26 Apr 2024 08:33:33 +0200 Subject: [PATCH 12/24] 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 033970efd..af36ccb67 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -622,7 +622,7 @@ Result vector `x` returns the approximate solution that minimizes the 2-norm \( `b`: Shall be a rank-1 array containing the right-hand-side vector. It is an `intent(in)` argument. -`cond` (optional): Singular value cut-off threshold for rank evaluation: `s_i >= cond*maxval(s), i=1:rank`. Shall be a scalar, `intent(in)` argument. +`cond` (optional): Shall be a scalar `real` value cut-off threshold for rank evaluation: `s_i >= cond*maxval(s), i=1:rank`. Shall be a scalar, `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 47d4793068ada2d49e69f5935a40ebd4fe80e777 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 26 Apr 2024 08:33:46 +0200 Subject: [PATCH 13/24] 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 af36ccb67..ffbccaa28 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -624,7 +624,7 @@ Result vector `x` returns the approximate solution that minimizes the 2-norm \( `cond` (optional): Shall be a scalar `real` value cut-off threshold for rank evaluation: `s_i >= cond*maxval(s), i=1:rank`. Shall be a scalar, `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. `rank` (optional): Shall be an `integer` scalar value, that contains the rank of input matrix `A`. This is an `intent(out)` argument. From 6d6f49a12c9596efbed820ddf39b1170274e9ed5 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 26 Apr 2024 08:33:57 +0200 Subject: [PATCH 14/24] 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 ffbccaa28..d928464c4 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -634,7 +634,7 @@ Result vector `x` returns the approximate solution that minimizes the 2-norm \( Returns an array value that represents the solution to the least squares system. -Raises `LINALG_ERROR` if the underlying SVD process did not converge. +Raises `LINALG_ERROR` if the underlying Singular Value Decomposition process did not converge. Raises `LINALG_VALUE_ERROR` if the matrix and rhs vectors have invalid/incompatible sizes. Exceptions trigger an `error stop`. From 37857f012e74c857b602d0b35f7aadeccb4ec0a9 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 26 Apr 2024 08:34:04 +0200 Subject: [PATCH 15/24] 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 d928464c4..f43ef89a7 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -635,7 +635,7 @@ Result vector `x` returns the approximate solution that minimizes the 2-norm \( Returns an array value that represents the solution to the least squares system. Raises `LINALG_ERROR` if the underlying Singular Value Decomposition process did not converge. -Raises `LINALG_VALUE_ERROR` if the matrix and rhs vectors have invalid/incompatible sizes. +Raises `LINALG_VALUE_ERROR` if the matrix and right-hand-side vector have invalid/incompatible sizes. Exceptions trigger an `error stop`. ### Example From 09f1a5b53a075938e699ebc54bd7d62222a1a549 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 26 Apr 2024 08:35:07 +0200 Subject: [PATCH 16/24] remove unused state variable --- example/linalg/example_lstsq.f90 | 2 -- 1 file changed, 2 deletions(-) diff --git a/example/linalg/example_lstsq.f90 b/example/linalg/example_lstsq.f90 index 2b36b0577..4c539fa27 100644 --- a/example/linalg/example_lstsq.f90 +++ b/example/linalg/example_lstsq.f90 @@ -1,9 +1,7 @@ program example_lstsq use stdlib_linalg_constants, only: dp use stdlib_linalg, only: lstsq - use stdlib_linalg_state, only: linalg_state_type implicit none - type(linalg_state_type) :: err integer, allocatable :: x(:),y(:) real(dp), allocatable :: A(:,:),b(:),coef(:) From 63d297ecdc961c44f53b451f409828d69b9f8dda Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 26 Apr 2024 08:46:24 +0200 Subject: [PATCH 17/24] add test with random input --- test/linalg/test_linalg_lstsq.fypp | 34 ++++++++++++++++++++++++++---- 1 file changed, 30 insertions(+), 4 deletions(-) diff --git a/test/linalg/test_linalg_lstsq.fypp b/test/linalg/test_linalg_lstsq.fypp index 435d71129..29cd07519 100644 --- a/test/linalg/test_linalg_lstsq.fypp +++ b/test/linalg/test_linalg_lstsq.fypp @@ -13,7 +13,6 @@ module test_linalg_least_squares public :: test_least_squares contains - !> Solve sample least squares problems subroutine test_least_squares(tests) @@ -24,15 +23,16 @@ module test_linalg_least_squares #:for rk,rt,ri in REAL_KINDS_TYPES #:if rk!="xdp" - tests = [tests,new_unittest("lease_squares_${ri}$",test_lstsq_one_${ri}$)] + tests = [tests,new_unittest("least_squares_${ri}$",test_lstsq_one_${ri}$), & + new_unittest("least_squares_randm_${ri}$",test_lstsq_random_${ri}$)] #:endif #:endfor end subroutine test_least_squares - - !> Simple polynomial fit + #:for rk,rt,ri in REAL_KINDS_TYPES #:if rk!="xdp" + !> Simple polynomial fit subroutine test_lstsq_one_${ri}$(error) type(error_type), allocatable, intent(out) :: error @@ -64,6 +64,32 @@ module test_linalg_least_squares end subroutine test_lstsq_one_${ri}$ + !> Fit from random array + subroutine test_lstsq_random_${ri}$(error) + type(error_type), allocatable, intent(out) :: error + + type(linalg_state_type) :: state + integer(ilp), parameter :: n = 12, m = 3 + ${rt}$ :: xsol(m),x(m),y(n),A(n,m) + + ! Random coefficient matrix and solution + call random_number(A) + call random_number(xsol) + + ! Compute rhs + y = matmul(A,xsol) + + ! Find polynomial + x = lstsq(A,y,err=state) + + call check(error,state%ok(),state%print()) + if (allocated(error)) return + + call check(error, all(abs(x-xsol)<1.0e-6_${rk}$), 'data converged') + if (allocated(error)) return + + end subroutine test_lstsq_random_${ri}$ + #:endif #:endfor From 78a8f6e2b0a2a893450ecf678c57095dab9780fe Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 26 Apr 2024 08:50:51 +0200 Subject: [PATCH 18/24] cleanup --- src/stdlib_linalg_least_squares.fypp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/stdlib_linalg_least_squares.fypp b/src/stdlib_linalg_least_squares.fypp index 4169d7802..662610547 100644 --- a/src/stdlib_linalg_least_squares.fypp +++ b/src/stdlib_linalg_least_squares.fypp @@ -176,7 +176,8 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares case (0) ! Success case (:-1) - err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid problem size a=[',lda,',',n,'], b[',ldb,',',nrhs,']') + err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid problem size a=',[lda,n], & + ', b=',[ldb,nrhs]) case (1:) err0 = linalg_state_type(this,LINALG_ERROR,'SVD did not converge.') case default From dafa081dabbdb1143e696da030d03e5c8a0e865e Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 26 Apr 2024 08:56:07 +0200 Subject: [PATCH 19/24] fix comment --- 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 f43ef89a7..bb74b2b36 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -618,9 +618,9 @@ Result vector `x` returns the approximate solution that minimizes the 2-norm \( ### Arguments -`a`: Shall be a rank-2 square array containing the coefficient matrix. It is an `intent(inout)` argument. +`a`: Shall be a rank-2 `real` or `complex` array containing the coefficient matrix. It is an `intent(inout)` argument. -`b`: Shall be a rank-1 array containing the right-hand-side vector. It is an `intent(in)` argument. +`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. `cond` (optional): Shall be a scalar `real` value cut-off threshold for rank evaluation: `s_i >= cond*maxval(s), i=1:rank`. Shall be a scalar, `intent(in)` argument. From 42f2de1bfb5ae3dd4817a624baaa66be7f49dba4 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 29 Apr 2024 10:14:31 +0200 Subject: [PATCH 20/24] subroutine interface --- src/stdlib_linalg.fypp | 52 ++++++++ src/stdlib_linalg_least_squares.fypp | 181 +++++++++++++++++++++------ 2 files changed, 192 insertions(+), 41 deletions(-) diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index 7384d86b6..ce9be7a0e 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -268,6 +268,58 @@ module stdlib_linalg #:endfor end interface lstsq + ! Least squares solution to system Ax=b, i.e. such that the 2-norm abs(b-Ax) is minimized. + interface solve_lstsq + !! version: experimental + !! + !! Computes the squares solution to system \( A \cdot x = b \). + !! ([Specification](../page/specs/stdlib_linalg.html#XXXXXxXxxxxxxx))xxxxx + !! + !!### Summary + !! Subroutine interface for computing least squares, i.e. the 2-norm \( || (b-A \cdot x ||_2 \) minimizing solution. + !! + !!### Description + !! + !! This interface provides methods for computing the least squares of a linear matrix system using + !! a subroutine. Supported data types include `real` and `complex`. If pre-allocated work spaces + !! are provided, no internal memory allocations take place when using this interface. + !! + !!@note The solution is based on LAPACK's singular value decomposition `*GELSD` methods. + !!@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 subroutine stdlib_linalg_${ri}$_solve_lstsq_${ndsuf}$(a,b,x,real_storage,int_storage,& + #{if rt.startswith('c')}#cmpl_storage,#{endif}#cond,overwrite_a,rank,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] real working storage space + real(${rk}$), optional, intent(inout), target :: real_storage(:) + !> [optional] integer working storage space + integer(ilp), optional, intent(inout), target :: int_storage(:) + #:if rt.startswith('complex') + !> [optional] complex working storage space + ${rt}$, optional, intent(inout), target :: cmpl_storage(:) + #:endif + !> [optional] cutoff for rank evaluation: singular values s(i)<=cond*maxval(s) are considered 0. + real(${rk}$), optional, intent(in) :: cond + !> [optional] Can A,b data be overwritten and destroyed? + logical(lk), optional, intent(in) :: overwrite_a + !> [optional] Return rank of A + integer(ilp), optional, intent(out) :: rank + !> [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_lstsq_${ndsuf}$ + #:endif + #:endfor + #:endfor + end interface solve_lstsq + interface det !! version: experimental !! diff --git a/src/stdlib_linalg_least_squares.fypp b/src/stdlib_linalg_least_squares.fypp index 662610547..ffa7958cc 100644 --- a/src/stdlib_linalg_least_squares.fypp +++ b/src/stdlib_linalg_least_squares.fypp @@ -16,10 +16,28 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares contains + elemental subroutine handle_gelsd_info(info,lda,n,ldb,nrhs,err) + integer(ilp), intent(in) :: info,lda,n,ldb,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 a=',[lda,n], & + ', b=',[ldb,nrhs]) + case (1:) + err = linalg_state_type(this,LINALG_ERROR,'SVD did not converge.') + case default + err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'catastrophic error') + + end subroutine handle_gelsd_info + #:for rk,rt,ri in RC_KINDS_TYPES #:if rk!="xdp" - ! Workspace needed by gesv - elemental subroutine ${ri}$gesv_space(m,n,nrhs,lrwork,liwork,lcwork) + ! Workspace needed by gelsd + elemental subroutine ${ri}$gelsd_space(m,n,nrhs,lrwork,liwork,lcwork) integer(ilp), intent(in) :: m,n,nrhs integer(ilp), intent(out) :: lrwork,liwork,lcwork @@ -53,7 +71,7 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares lcwork = ceiling(1.25*lcwork,kind=ilp) liwork = ceiling(1.25*liwork,kind=ilp) - end subroutine ${ri}$gesv_space + end subroutine ${ri}$gelsd_space #:endif #:endfor @@ -93,17 +111,69 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares !> [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}$ + ${rt}$, allocatable, target :: x${nd}$ + + ! Initialize solution with the shape of the rhs + allocate(x,mold=b) + + call stdlib_linalg_${ri}$_solve_lstsq_${ndsuf}$(a,b,x,& + cond=cond,overwrite_a=overwrite_a,rank=rank,err=err) + + end function stdlib_linalg_${ri}$_lstsq_${ndsuf}$ + + ! Compute the least-squares solution to a real system of linear equations Ax = b + module subroutine stdlib_linalg_${ri}$_solve_lstsq_${ndsuf}$(a,b,x, & + real_storage,int_storage#{if rt.startswith('c')}#,cmpl_storage#{endif}#,cond,overwrite_a,rank,err) + + !!### Summary + !! Compute least-squares solution to a real system of linear equations \( Ax = b \) + !! + !!### Description + !! + !! This function computes the least-squares solution of a linear matrix problem. + !! + !! param: a Input matrix of size [m,n]. + !! param: b Right-hand-side vector of size [n] or matrix of size [n,nrhs]. + !! param: cond [optional] Real input threshold indicating that singular values `s_i <= cond*maxval(s)` + !! do not contribute to the matrix rank. + !! param: overwrite_a [optional] Flag indicating if the input matrix can be overwritten. + !! param: rank [optional] integer flag returning matrix rank. + !! param: err [optional] State return flag. + !! return: x Solution vector of size [n] or solution matrix of size [n,nrhs]. + !! + !> 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] real working storage space + real(${rk}$), optional, intent(inout), target :: real_storage(:) + !> [optional] integer working storage space + integer(ilp), optional, intent(inout), target :: int_storage(:) + #:if rt.startswith('c') + !> [optional] complex working storage space + ${rt}$, optional, intent(inout), target :: cmpl_storage(:) + #:endif + !> [optional] cutoff for rank evaluation: singular values s(i)<=cond*maxval(s) are considered 0. + real(${rk}$), optional, intent(in) :: cond + !> [optional] Can A,b data be overwritten and destroyed? + logical(lk), optional, intent(in) :: overwrite_a + !> [optional] Return rank of A + integer(ilp), optional, intent(out) :: rank + !> [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) :: m,n,lda,ldb,nrhs,info,mnmin,mnmax,arank,lrwork,liwork,lcwork - integer(ilp), allocatable :: iwork(:) + integer(ilp) :: m,n,lda,ldb,nrhs,ldx,nrhsx,info,mnmin,mnmax,arank,lrwork,liwork,lcwork + integer(ilp) :: nrs,nis,ncs + integer(ilp), pointer :: iwork(:) logical(lk) :: copy_a real(${rk}$) :: acond,rcond - real(${rk}$), allocatable :: singular(:),rwork(:) - ${rt}$, pointer :: xmat(:,:),amat(:,:) - ${rt}$, allocatable :: cwork(:) + real(${rk}$), allocatable :: singular(:) + real(${rk}$), pointer :: rwork(:) + ${rt}$, pointer :: xmat(:,:),amat(:,:),cwork(:) ! Problem sizes m = size(a,1,kind=ilp) @@ -111,15 +181,17 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares 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 mnmin = min(m,n) mnmax = max(m,n) - if (lda<1 .or. n<1 .or. ldb<1 .or. ldb/=m) then + if (lda<1 .or. n<1 .or. ldb<1 .or. ldb/=m .or. ldx/=m) then err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid sizes: a=',[lda,n], & - 'b=',[ldb,nrhs]) - allocate(x${nde}$) + 'b=',[ldb,nrhs],' x=',[ldx,nrhsx]) call linalg_error_handling(err0,err) if (present(rank)) rank = 0 + return end if ! Can A be overwritten? By default, do not overwrite @@ -137,7 +209,7 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares endif ! Initialize solution with the rhs - allocate(x,source=b) + x = b xmat(1:n,1:nrhs) => x ! Singular values array (in decreasing order) @@ -153,44 +225,71 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares endif if (rcond<0) rcond = epsilon(0.0_${rk}$)*mnmax - ! Allocate working space - call ${ri}$gesv_space(m,n,nrhs,lrwork,liwork,lcwork) - #:if rt.startswith('complex') - allocate(rwork(lrwork),cwork(lcwork),iwork(liwork)) - #:else - allocate(rwork(lrwork),iwork(liwork)) - #:endif + ! Get working space size + call ${ri}$gelsd_space(m,n,nrhs,lrwork,liwork,lcwork) + + ! Real working space + if (present(real_storage)) then + rwork => real_storage + else + allocate(rwork(lrwork)) + endif + nrs = size(rwork,kind=ilp) - ! Solve system using singular value decomposition - #:if rt.startswith('complex') - call gelsd(m,n,nrhs,amat,lda,xmat,ldb,singular,rcond,arank,cwork,lrwork,rwork,iwork,info) - #:else - call gelsd(m,n,nrhs,amat,lda,xmat,ldb,singular,rcond,arank,rwork,lrwork,iwork,info) - #:endif + ! Integer working space + if (present(int_storage)) then + iwork => int_storage + else + allocate(iwork(liwork)) + endif + nis = size(iwork,kind=ilp) + + #:if rt.startswith('complex') + ! Complex working space + if (present(cmpl_storage)) then + cwork => cmpl_storage + else + allocate(cwork(lcwork)) + endif + ncs = size(iwork,kind=ilp) + #:endif + + if (nrs=',lrwork, & + ', int=',nis,' should be >=',liwork & + #{if rt.startswith('complex')}#,', cmplx=',ncs,' should be >=',lcwork#{endif}#) + + else + ! Solve system using singular value decomposition + call gelsd(m,n,nrhs,amat,lda,xmat,ldb,singular,rcond,arank, & + #:if rt.startswith('complex') + cwork,nrs,rwork,iwork,info) + #:else + rwork,nrs,iwork,info) + #:endif + + endif + ! The condition number of A in the 2-norm = S(1)/S(min(m,n)). acond = singular(1)/singular(mnmin) ! Process output - select case (info) - case (0) - ! Success - case (:-1) - err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid problem size a=',[lda,n], & - ', b=',[ldb,nrhs]) - case (1:) - err0 = linalg_state_type(this,LINALG_ERROR,'SVD did not converge.') - case default - err0 = linalg_state_type(this,LINALG_INTERNAL_ERROR,'catastrophic error') - end select - - if (copy_a) deallocate(amat) + call handle_gelsd_info(info,lda,n,ldb,nrhs,err0) ! Process output and return - call linalg_error_handling(err0,err) + 1 if (copy_a) deallocate(amat) if (present(rank)) rank = arank + if (.not.present(real_storage)) deallocate(rwork) + if (.not.present(int_storage)) deallocate(iwork) + #:if rt.startswith('complex') + if (.not.present(cmpl_storage)) deallocate(cwork) + #:endif + call linalg_error_handling(err0,err) - end function stdlib_linalg_${ri}$_lstsq_${ndsuf}$ + end subroutine stdlib_linalg_${ri}$_solve_lstsq_${ndsuf}$ #:endif #:endfor From 2e370f45438a637c90b608609f6edc0ddbdc6ef6 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 29 Apr 2024 10:34:03 +0200 Subject: [PATCH 21/24] export working space interface --- src/stdlib_linalg.fypp | 38 +++++++++++++- src/stdlib_linalg_least_squares.fypp | 77 ++++++++++++++++++++-------- 2 files changed, 93 insertions(+), 22 deletions(-) diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index ce9be7a0e..614fedb1b 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -21,6 +21,8 @@ module stdlib_linalg public :: diag public :: eye public :: lstsq + public :: lstsq_space + public :: solve_lstsq public :: trace public :: outer_product public :: kronecker_product @@ -273,7 +275,7 @@ module stdlib_linalg !! version: experimental !! !! Computes the squares solution to system \( A \cdot x = b \). - !! ([Specification](../page/specs/stdlib_linalg.html#XXXXXxXxxxxxxx))xxxxx + !! ([Specification](../page/specs/stdlib_linalg.html)) !! !!### Summary !! Subroutine interface for computing least squares, i.e. the 2-norm \( || (b-A \cdot x ||_2 \) minimizing solution. @@ -291,7 +293,7 @@ module stdlib_linalg #:for rk,rt,ri in RC_KINDS_TYPES #:if rk!="xdp" module subroutine stdlib_linalg_${ri}$_solve_lstsq_${ndsuf}$(a,b,x,real_storage,int_storage,& - #{if rt.startswith('c')}#cmpl_storage,#{endif}#cond,overwrite_a,rank,err) + #{if rt.startswith('c')}#cmpl_storage,#{endif}#cond,singvals,overwrite_a,rank,err) !> Input matrix a[n,n] ${rt}$, intent(inout), target :: a(:,:) !> Right hand side vector or array, b[n] or b[n,nrhs] @@ -308,6 +310,8 @@ module stdlib_linalg #:endif !> [optional] cutoff for rank evaluation: singular values s(i)<=cond*maxval(s) are considered 0. real(${rk}$), optional, intent(in) :: cond + !> [optional] list of singular values [min(m,n)], in descending magnitude order, returned by the SVD + real(${rk}$), optional, intent(out), target :: singvals(:) !> [optional] Can A,b data be overwritten and destroyed? logical(lk), optional, intent(in) :: overwrite_a !> [optional] Return rank of A @@ -320,6 +324,36 @@ module stdlib_linalg #:endfor end interface solve_lstsq + ! Return the working array space required by the least squares solver + interface lstsq_space + !! version: experimental + !! + !! Computes the integer, real [, complex] working space required by the least-squares solver + !! ([Specification](../page/specs/stdlib_linalg.html)) + !! + !!### Description + !! + !! This interface provides sizes of integer, real [, complex] working spaces required by the + !! least-squares solver. These sizes can be used to pre-allocated working arrays in case several + !! repeated least-squares solutions to a same system are sought. If pre-allocated working arrays + !! are provided, no internal allocations will take place. + !! + #:for nd,ndsuf,nde in ALL_RHS + #:for rk,rt,ri in RC_KINDS_TYPES + #:if rk!="xdp" + pure module subroutine stdlib_linalg_${ri}$_lstsq_space_${ndsuf}$(a,b,lrwork,liwork#{if rt.startswith('c')}#,lcwork#{endif}#) + !> Input matrix a[m,n] + ${rt}$, intent(inout), target :: a(:,:) + !> Right hand side vector or array, b[n] or b[n,nrhs] + ${rt}$, intent(in) :: b${nd}$ + !> Size of the working space arrays + integer(ilp), intent(out) :: lrwork,liwork#{if rt.startswith('c')}#,lcwork#{endif}# + end subroutine stdlib_linalg_${ri}$_lstsq_space_${ndsuf}$ + #:endif + #:endfor + #:endfor + end interface lstsq_space + interface det !! version: experimental !! diff --git a/src/stdlib_linalg_least_squares.fypp b/src/stdlib_linalg_least_squares.fypp index ffa7958cc..fe1d801e8 100644 --- a/src/stdlib_linalg_least_squares.fypp +++ b/src/stdlib_linalg_least_squares.fypp @@ -31,6 +31,7 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares err = linalg_state_type(this,LINALG_ERROR,'SVD did not converge.') case default err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'catastrophic error') + end select end subroutine handle_gelsd_info @@ -80,6 +81,26 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares #:for rk,rt,ri in RC_KINDS_TYPES #:if rk!="xdp" + ! Compute the integer, real, [complex] working space requested byu the least squares procedure + pure module subroutine stdlib_linalg_${ri}$_lstsq_space_${ndsuf}$(a,b,lrwork,liwork#{if rt.startswith('c')}#,lcwork#{endif}#) + !> Input matrix a[m,n] + ${rt}$, intent(inout), target :: a(:,:) + !> Right hand side vector or array, b[n] or b[n,nrhs] + ${rt}$, intent(in) :: b${nd}$ + !> Size of the working space arrays + integer(ilp), intent(out) :: lrwork,liwork + integer(ilp) #{if rt.startswith('c')}#, intent(out)#{endif}# :: lcwork + + integer(ilp) :: m,n,nrhs + + m = size(a,1,kind=ilp) + n = size(a,2,kind=ilp) + nrhs = size(b,kind=ilp)/size(b,1,kind=ilp) + + call ${ri}$gelsd_space(m,n,nrhs,lrwork,liwork,lcwork) + + end subroutine stdlib_linalg_${ri}$_lstsq_space_${ndsuf}$ + ! Compute the least-squares solution to a real system of linear equations Ax = b module function stdlib_linalg_${ri}$_lstsq_${ndsuf}$(a,b,cond,overwrite_a,rank,err) result(x) !!### Summary @@ -98,7 +119,7 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares !! param: err [optional] State return flag. !! return: x Solution vector of size [n] or solution matrix of size [n,nrhs]. !! - !> Input matrix a[n,n] + !> Input matrix a[m,n] ${rt}$, intent(inout), target :: a(:,:) !> Right hand side vector or array, b[n] or b[n,nrhs] ${rt}$, intent(in) :: b${nd}$ @@ -122,8 +143,8 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares end function stdlib_linalg_${ri}$_lstsq_${ndsuf}$ ! Compute the least-squares solution to a real system of linear equations Ax = b - module subroutine stdlib_linalg_${ri}$_solve_lstsq_${ndsuf}$(a,b,x, & - real_storage,int_storage#{if rt.startswith('c')}#,cmpl_storage#{endif}#,cond,overwrite_a,rank,err) + module subroutine stdlib_linalg_${ri}$_solve_lstsq_${ndsuf}$(a,b,x,real_storage,int_storage, & + #{if rt.startswith('c')}#cmpl_storage,#{endif}# cond,singvals,overwrite_a,rank,err) !!### Summary !! Compute least-squares solution to a real system of linear equations \( Ax = b \) @@ -134,12 +155,18 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares !! !! param: a Input matrix of size [m,n]. !! param: b Right-hand-side vector of size [n] or matrix of size [n,nrhs]. + !! param: x Solution vector of size [n] or solution matrix of size [n,nrhs]. + !! param: real_storage [optional] Real working space + !! param: int_storage [optional] Integer working space + #:if rt.startswith('c') + !! param: cmpl_storage [optional] Complex working space + #:endif !! param: cond [optional] Real input threshold indicating that singular values `s_i <= cond*maxval(s)` !! do not contribute to the matrix rank. + !! param: singvals [optional] Real array of size [min(m,n)] returning a list of singular values. !! param: overwrite_a [optional] Flag indicating if the input matrix can be overwritten. !! param: rank [optional] integer flag returning matrix rank. - !! param: err [optional] State return flag. - !! return: x Solution vector of size [n] or solution matrix of size [n,nrhs]. + !! param: err [optional] State return flag. !! !> Input matrix a[n,n] ${rt}$, intent(inout), target :: a(:,:) @@ -157,6 +184,8 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares #:endif !> [optional] cutoff for rank evaluation: singular values s(i)<=cond*maxval(s) are considered 0. real(${rk}$), optional, intent(in) :: cond + !> [optional] list of singular values [min(m,n)], in descending magnitude order, returned by the SVD + real(${rk}$), optional, intent(out), target :: singvals(:) !> [optional] Can A,b data be overwritten and destroyed? logical(lk), optional, intent(in) :: overwrite_a !> [optional] Return rank of A @@ -167,12 +196,11 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares !! Local variables type(linalg_state_type) :: err0 integer(ilp) :: m,n,lda,ldb,nrhs,ldx,nrhsx,info,mnmin,mnmax,arank,lrwork,liwork,lcwork - integer(ilp) :: nrs,nis,ncs + integer(ilp) :: nrs,nis,ncs,nsvd integer(ilp), pointer :: iwork(:) logical(lk) :: copy_a real(${rk}$) :: acond,rcond - real(${rk}$), allocatable :: singular(:) - real(${rk}$), pointer :: rwork(:) + real(${rk}$), pointer :: rwork(:),singular(:) ${rt}$, pointer :: xmat(:,:),amat(:,:),cwork(:) ! Problem sizes @@ -213,8 +241,14 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares xmat(1:n,1:nrhs) => x ! Singular values array (in decreasing order) - allocate(singular(mnmin)) - + if (present(singvals)) then + singular => singvals + nsvd = size(singular,kind=ilp) + else + allocate(singular(mnmin)) + nsvd = mnmin + endif + ! rcond is used to determine the effective rank of A. ! Singular values S(i) <= RCOND*maxval(S) are treated as zero. ! Use same default value as NumPy @@ -254,12 +288,14 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares ncs = size(iwork,kind=ilp) #:endif - if (nrs=',lrwork, & - ', int=',nis,' should be >=',liwork & - #{if rt.startswith('complex')}#,', cmplx=',ncs,' should be >=',lcwork#{endif}#) + ', int=',nis,' should be >=',liwork, & + #{if rt.startswith('complex')}#', cmplx=',ncs,' should be >=',lcwork, &#{endif}# + ', singv=',nsvd,' should be >=',mnmin) else @@ -270,23 +306,24 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares #:else rwork,nrs,iwork,info) #:endif + + ! The condition number of A in the 2-norm = S(1)/S(min(m,n)). + acond = singular(1)/singular(mnmin) + + ! Process output + call handle_gelsd_info(info,lda,n,ldb,nrhs,err0) endif - ! The condition number of A in the 2-norm = S(1)/S(min(m,n)). - acond = singular(1)/singular(mnmin) - - ! Process output - call handle_gelsd_info(info,lda,n,ldb,nrhs,err0) - ! Process output and return - 1 if (copy_a) deallocate(amat) + if (copy_a) deallocate(amat) if (present(rank)) rank = arank if (.not.present(real_storage)) deallocate(rwork) if (.not.present(int_storage)) deallocate(iwork) #:if rt.startswith('complex') if (.not.present(cmpl_storage)) deallocate(cwork) #:endif + if (.not.present(singvals)) deallocate(singular) call linalg_error_handling(err0,err) end subroutine stdlib_linalg_${ri}$_solve_lstsq_${ndsuf}$ From 70da4c8c48e712bf58edb772f718cc8091869d32 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 29 Apr 2024 19:07:54 +0200 Subject: [PATCH 22/24] 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 614fedb1b..e7eb503a3 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 18151e63878b58b83eb8e7c705f92172766168a2 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Tue, 30 Apr 2024 09:18:56 +0200 Subject: [PATCH 23/24] `solve_lstsq` example program --- example/linalg/CMakeLists.txt | 3 +- .../{example_lstsq.f90 => example_lstsq1.f90} | 5 ++- example/linalg/example_lstsq2.f90 | 40 +++++++++++++++++++ 3 files changed, 45 insertions(+), 3 deletions(-) rename example/linalg/{example_lstsq.f90 => example_lstsq1.f90} (85%) create mode 100644 example/linalg/example_lstsq2.f90 diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt index 4ecef6bc0..729f117d3 100644 --- a/example/linalg/CMakeLists.txt +++ b/example/linalg/CMakeLists.txt @@ -18,6 +18,7 @@ ADD_EXAMPLE(state1) ADD_EXAMPLE(state2) ADD_EXAMPLE(blas_gemv) ADD_EXAMPLE(lapack_getrf) -ADD_EXAMPLE(lstsq) +ADD_EXAMPLE(lstsq1) +ADD_EXAMPLE(lstsq2) ADD_EXAMPLE(determinant) ADD_EXAMPLE(determinant2) diff --git a/example/linalg/example_lstsq.f90 b/example/linalg/example_lstsq1.f90 similarity index 85% rename from example/linalg/example_lstsq.f90 rename to example/linalg/example_lstsq1.f90 index 4c539fa27..dfbc031ea 100644 --- a/example/linalg/example_lstsq.f90 +++ b/example/linalg/example_lstsq1.f90 @@ -1,4 +1,5 @@ -program example_lstsq +! Least-squares solver: functional interface +program example_lstsq1 use stdlib_linalg_constants, only: dp use stdlib_linalg, only: lstsq implicit none @@ -22,4 +23,4 @@ program example_lstsq ! parabola: -0.42857142857141695 1.1428571428571503 4.2857142857142811 -end program example_lstsq +end program example_lstsq1 diff --git a/example/linalg/example_lstsq2.f90 b/example/linalg/example_lstsq2.f90 new file mode 100644 index 000000000..389b20d4c --- /dev/null +++ b/example/linalg/example_lstsq2.f90 @@ -0,0 +1,40 @@ +! Demonstrate expert subroutine interface with pre-allocated arrays +program example_lstsq2 + use stdlib_linalg_constants, only: dp,ilp + use stdlib_linalg, only: solve_lstsq, lstsq_space, linalg_state_type + implicit none + + integer, allocatable :: x(:),y(:) + real(dp), allocatable :: A(:,:),b(:),coef(:),real_space(:),singvals(:) + integer(ilp), allocatable :: int_space(:) + integer(ilp) :: lrwork,liwork,arank + + ! Data set + x = [1, 2, 2] + y = [5, 13, 25] + + ! Fit three points using a parabola, least squares method + ! A = [1 x x**2] + A = reshape([[1,1,1],x,x**2],[3,3]) + b = y + + ! Get storage sizes for the arrays and pre-allocate data + call lstsq_space(A,b,lrwork,liwork) + allocate(coef(size(x)),real_space(lrwork),int_space(liwork),singvals(minval(shape(A)))) + + ! Solve coefficients of y = coef(1) + x*coef(2) + x^2*coef(3) + ! with no internal allocations + call solve_lstsq(A,b,x=coef, & + real_storage=real_space, & + int_storage=int_space, & + singvals=singvals, & + overwrite_a=.true., & + rank=arank) + + print *, 'parabola: ',coef + ! parabola: -0.42857142857141695 1.1428571428571503 4.2857142857142811 + print *, 'rank: ',arank + ! rank: 2 + + +end program example_lstsq2 From c55583d55eddeccd5c9e5da61e2a078d60ae4caf Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Tue, 30 Apr 2024 09:36:21 +0200 Subject: [PATCH 24/24] document `solve_lstsq` and `lstsq_space` --- doc/specs/stdlib_linalg.md | 84 +++++++++++++++++++++++++++- src/stdlib_linalg.fypp | 6 +- src/stdlib_linalg_least_squares.fypp | 2 +- 3 files changed, 86 insertions(+), 6 deletions(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index edb12b36f..a081291d0 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -620,7 +620,7 @@ Result vector `x` returns the approximate solution that minimizes the 2-norm \( `a`: Shall be a rank-2 `real` or `complex` array containing the coefficient matrix. It is an `intent(inout)` argument. -`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 one or more right-hand-side vector(s), each in its leading dimension. It is an `intent(in)` argument. `cond` (optional): Shall be a scalar `real` value cut-off threshold for rank evaluation: `s_i >= cond*maxval(s), i=1:rank`. Shall be a scalar, `intent(in)` argument. @@ -632,6 +632,60 @@ Result vector `x` returns the approximate solution that minimizes the 2-norm \( ### Return value +Returns an array value of the same kind and rank as `b`, containing the solution(s) to the least squares system. + +Raises `LINALG_ERROR` if the underlying Singular Value Decomposition process did not converge. +Raises `LINALG_VALUE_ERROR` if the matrix and right-hand-side vector have invalid/incompatible sizes. +Exceptions trigger an `error stop`. + +### Example + +```fortran +{!example/linalg/example_lstsq1.f90!} +``` + +## `solve_lstsq` - Compute the least squares solution to a linear matrix equation (subroutine interface). + +### Status + +Experimental + +### Description + +This subroutine computes the least-squares solution to a linear matrix equation \( A \cdot x = b \). + +Result vector `x` returns the approximate solution that minimizes the 2-norm \( || A \cdot x - b ||_2 \), i.e., it contains the least-squares solution to the problem. Matrix `A` may be full-rank, over-determined, or under-determined. The solver is based on LAPACK's `*GELSD` backends. + +### Syntax + +`call ` [[stdlib_linalg(module):solve_lstsq(interface)]] `(a, b, x, [, real_storage, int_storage, [cmpl_storage, ] cond, singvals, overwrite_a, rank, err])` + +### Arguments + +`a`: Shall be a rank-2 `real` or `complex` array containing the coefficient matrix. It is an `intent(inout)` argument. + +`b`: Shall be a rank-1 or rank-2 array of the same kind as `a`, containing one or more right-hand-side vector(s), each in its leading dimension. It is an `intent(in)` argument. + +`x`: Shall be an array of same kind and rank as `b`, containing the solution(s) to the least squares system. It is an `intent(inout)` argument. + +`real_storage` (optional): Shall be a `real` rank-1 array of the same kind `a`, providing working storage for the solver. It minimum size can be determined with a call to [[stdlib_linalg(module):lstsq_space(interface)]]. It is an `intent(inout)` argument. + +`int_storage` (optional): Shall be an `integer` rank-1 array, providing working storage for the solver. It minimum size can be determined with a call to [[stdlib_linalg(module):lstsq_space(interface)]]. It is an `intent(inout)` argument. + +`cmpl_storage` (optional): For `complex` systems, it shall be a `complex` rank-1 array, providing working storage for the solver. It minimum size can be determined with a call to [[stdlib_linalg(module):lstsq_space(interface)]]. It is an `intent(inout)` argument. + +`cond` (optional): Shall be a scalar `real` value cut-off threshold for rank evaluation: `s_i >= cond*maxval(s), i=1:rank`. Shall be a scalar, `intent(in)` argument. + +`singvals` (optional): Shall be a `real` rank-1 array of the same kind `a` and size at least `minval(shape(a))`, returning the list of singular values `s(i)>=cond*maxval(s)`, in descending order of magnitude. It is an `intent(out)` 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. + +`rank` (optional): Shall be an `integer` scalar value, that contains the rank of input matrix `A`. This is an `intent(out)` argument. + +`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. + +### Return value + Returns an array value that represents the solution to the least squares system. Raises `LINALG_ERROR` if the underlying Singular Value Decomposition process did not converge. @@ -641,9 +695,35 @@ Exceptions trigger an `error stop`. ### Example ```fortran -{!example/linalg/example_lstsq.f90!} +{!example/linalg/example_lstsq2.f90!} ``` +## `lstsq_space` - Compute internal working space requirements for the least squares solver. + +### Status + +Experimental + +### Description + +This subroutine computes the internal working space requirements for the least-squares solver, [[stdlib_linalg(module):solve_lstsq(interface)]] . + +### Syntax + +`call ` [[stdlib_linalg(module):lstsq_space(interface)]] `(a, b, lrwork, liwork [, lcwork])` + +### Arguments + +`a`: Shall be a rank-2 `real` or `complex` array containing the linear system coefficient matrix. It is an `intent(in)` argument. + +`b`: Shall be a rank-1 or rank-2 array of the same kind as `a`, containing the system's right-hand-side vector(s). It is an `intent(in)` argument. + +`lrwork`: Shall be an `integer` scalar, that returns the minimum array size required for the `real` working storage to this system. + +`liwork`: Shall be an `integer` scalar, that returns the minimum array size required for the `integer` working storage to this system. + +`lcwork` (`complex` `a`, `b`): For a `complex` system, shall be an `integer` scalar, that returns the minimum array size required for the `complex` working storage to this system. + ## `det` - Computes the determinant of a square matrix ### Status diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index e7eb503a3..460aa0593 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -275,7 +275,7 @@ module stdlib_linalg !! version: experimental !! !! Computes the squares solution to system \( A \cdot x = b \). - !! ([Specification](../page/specs/stdlib_linalg.html)) + !! ([Specification](../page/specs/stdlib_linalg.html#solve-lstsq-compute-the-least-squares-solution-to-a-linear-matrix-equation-subroutine-interface)) !! !!### Summary !! Subroutine interface for computing least squares, i.e. the 2-norm \( || (b-A \cdot x ||_2 \) minimizing solution. @@ -329,7 +329,7 @@ module stdlib_linalg !! version: experimental !! !! Computes the integer, real [, complex] working space required by the least-squares solver - !! ([Specification](../page/specs/stdlib_linalg.html)) + !! ([Specification](../page/specs/stdlib_linalg.html#lstsq-space-compute-internal-working-space-requirements-for-the-least-squares-solver)) !! !!### Description !! @@ -343,7 +343,7 @@ module stdlib_linalg #:if rk!="xdp" pure module subroutine stdlib_linalg_${ri}$_lstsq_space_${ndsuf}$(a,b,lrwork,liwork#{if rt.startswith('c')}#,lcwork#{endif}#) !> Input matrix a[m,n] - ${rt}$, intent(inout), target :: a(:,:) + ${rt}$, intent(in), target :: a(:,:) !> Right hand side vector or array, b[n] or b[n,nrhs] ${rt}$, intent(in) :: b${nd}$ !> Size of the working space arrays diff --git a/src/stdlib_linalg_least_squares.fypp b/src/stdlib_linalg_least_squares.fypp index fe1d801e8..0b7c2e139 100644 --- a/src/stdlib_linalg_least_squares.fypp +++ b/src/stdlib_linalg_least_squares.fypp @@ -84,7 +84,7 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares ! Compute the integer, real, [complex] working space requested byu the least squares procedure pure module subroutine stdlib_linalg_${ri}$_lstsq_space_${ndsuf}$(a,b,lrwork,liwork#{if rt.startswith('c')}#,lcwork#{endif}#) !> Input matrix a[m,n] - ${rt}$, intent(inout), target :: a(:,:) + ${rt}$, intent(in), target :: a(:,:) !> Right hand side vector or array, b[n] or b[n,nrhs] ${rt}$, intent(in) :: b${nd}$ !> Size of the working space arrays