From d3889c89f04385c820e70d12ce7a6c5f6a705248 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 27 May 2024 22:06:55 +0200 Subject: [PATCH 1/4] fix: 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 65c7c0b9f..fe5b211ee 100644 --- a/src/stdlib_linalg_least_squares.fypp +++ b/src/stdlib_linalg_least_squares.fypp @@ -305,7 +305,7 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares else allocate(cwork(lcwork)) endif - ncs = size(iwork,kind=ilp) + ncs = size(cwork,kind=ilp) #:endif if (nrs Date: Mon, 27 May 2024 22:18:46 +0200 Subject: [PATCH 2/4] add test program --- test/linalg/test_linalg_lstsq.fypp | 37 +++++++++++++++++++++++++++++- 1 file changed, 36 insertions(+), 1 deletion(-) diff --git a/test/linalg/test_linalg_lstsq.fypp b/test/linalg/test_linalg_lstsq.fypp index cf57bb178..2d544944b 100644 --- a/test/linalg/test_linalg_lstsq.fypp +++ b/test/linalg/test_linalg_lstsq.fypp @@ -4,7 +4,7 @@ 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, only: lstsq,solve_lstsq use stdlib_linalg_state, only: linalg_state_type implicit none (type,external) @@ -20,6 +20,8 @@ module test_linalg_least_squares type(unittest_type), allocatable, intent(out) :: tests(:) allocate(tests(0)) + + tests = [tests,new_unittest("issue_823",test_issue_823)] #:for rk,rt,ri in REAL_KINDS_TYPES #:if rk!="xdp" @@ -100,6 +102,39 @@ module test_linalg_least_squares #:endif #:endfor + + ! Test issue #823 + subroutine test_issue_823(error) + type(error_type), allocatable, intent(out) :: error + + ! Dimension of the problem. + integer(ilp), parameter :: n = 42 + ! Data for the least-squares problem. + complex(dp) :: A(n+1, n), b(n+1), x_true(n), x_lstsq(n+1) + ! Internal variables. + real(dp), allocatable :: tmp(:, :, :), tmp_vec(:, :) + ! Error handler + type(linalg_state_type) :: state + + ! Zero-out data. + A = 0.0_dp ; b = 0.0_dp ; x_true = 0.0_dp ; x_lstsq = 0.0_dp + allocate(tmp(n+1, n, 2)) ; tmp = 0.0_dp + allocate(tmp_vec(n, 2)) ; tmp_vec = 0.0_dp + + ! Generate a random complex least-squares problem of size (n+1, n). + call random_number(tmp) ; call random_number(tmp_vec) + A = cmplx(tmp(:, :, 1), tmp(:, :, 2), kind=dp) + x_true = cmplx(tmp_vec(:, 1), tmp_vec(:, 2), kind=dp) + b = matmul(A, x_true) + + ! Solve the lstsq problem. + call solve_lstsq(A, b, x_lstsq, err=state) + + ! Check that no segfault occurred + call check(error,state%ok(),'issue 823 returned '//state%print()) + if (allocated(error)) return + + end subroutine test_issue_823 end module test_linalg_least_squares From 81ebb3f5ecae958f5c389a61a84e405451537a42 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Tue, 28 May 2024 07:16:52 +0200 Subject: [PATCH 3/4] ncs in gelsd --- 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 fe5b211ee..ccb885319 100644 --- a/src/stdlib_linalg_least_squares.fypp +++ b/src/stdlib_linalg_least_squares.fypp @@ -322,7 +322,7 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares ! 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) + cwork,ncs,rwork,iwork,info) #:else rwork,nrs,iwork,info) #:endif From 89dc99273cb3a86e901362a72154c39288a6eeda Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Tue, 28 May 2024 08:57:43 +0200 Subject: [PATCH 4/4] cleanup test --- test/linalg/test_linalg_lstsq.fypp | 21 ++++++++++++++------- 1 file changed, 14 insertions(+), 7 deletions(-) diff --git a/test/linalg/test_linalg_lstsq.fypp b/test/linalg/test_linalg_lstsq.fypp index 2d544944b..045ac843d 100644 --- a/test/linalg/test_linalg_lstsq.fypp +++ b/test/linalg/test_linalg_lstsq.fypp @@ -110,22 +110,25 @@ module test_linalg_least_squares ! Dimension of the problem. integer(ilp), parameter :: n = 42 ! Data for the least-squares problem. - complex(dp) :: A(n+1, n), b(n+1), x_true(n), x_lstsq(n+1) + complex(dp) :: A(n+1, n), b(n+1), x_true(n), x_lstsq(n) ! Internal variables. real(dp), allocatable :: tmp(:, :, :), tmp_vec(:, :) ! Error handler type(linalg_state_type) :: state ! Zero-out data. - A = 0.0_dp ; b = 0.0_dp ; x_true = 0.0_dp ; x_lstsq = 0.0_dp - allocate(tmp(n+1, n, 2)) ; tmp = 0.0_dp - allocate(tmp_vec(n, 2)) ; tmp_vec = 0.0_dp + A = 0.0_dp + b = 0.0_dp + x_lstsq = 0.0_dp + allocate(tmp(n+1, n, 2), tmp_vec(n, 2), source=0.0_dp) ! Generate a random complex least-squares problem of size (n+1, n). - call random_number(tmp) ; call random_number(tmp_vec) - A = cmplx(tmp(:, :, 1), tmp(:, :, 2), kind=dp) + call random_number(tmp) + call random_number(tmp_vec) + + A = cmplx(tmp(:, :, 1), tmp(:, :, 2), kind=dp) x_true = cmplx(tmp_vec(:, 1), tmp_vec(:, 2), kind=dp) - b = matmul(A, x_true) + b = matmul(A, x_true) ! Solve the lstsq problem. call solve_lstsq(A, b, x_lstsq, err=state) @@ -134,6 +137,10 @@ module test_linalg_least_squares call check(error,state%ok(),'issue 823 returned '//state%print()) if (allocated(error)) return + ! Check that least squares are verified + call check(error,all(abs(x_true-x_lstsq)