Skip to content

Commit dc707f8

Browse files
authored
linalg: fix typo in least-squares complex space allocation (#830)
2 parents 718ac3a + 89dc992 commit dc707f8

File tree

2 files changed

+45
-3
lines changed

2 files changed

+45
-3
lines changed

src/stdlib_linalg_least_squares.fypp

+2-2
Original file line numberDiff line numberDiff line change
@@ -305,7 +305,7 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares
305305
else
306306
allocate(cwork(lcwork))
307307
endif
308-
ncs = size(iwork,kind=ilp)
308+
ncs = size(cwork,kind=ilp)
309309
#:endif
310310

311311
if (nrs<lrwork .or. nis<liwork#{if rt.startswith('c')}# .or. ncs<lcwork#{endif}# &
@@ -322,7 +322,7 @@ submodule (stdlib_linalg) stdlib_linalg_least_squares
322322
! Solve system using singular value decomposition
323323
call gelsd(m,n,nrhs,amat,lda,xmat,ldb,singular,rcond,arank, &
324324
#:if rt.startswith('complex')
325-
cwork,nrs,rwork,iwork,info)
325+
cwork,ncs,rwork,iwork,info)
326326
#:else
327327
rwork,nrs,iwork,info)
328328
#:endif

test/linalg/test_linalg_lstsq.fypp

+43-1
Original file line numberDiff line numberDiff line change
@@ -4,7 +4,7 @@
44
module test_linalg_least_squares
55
use testdrive, only: error_type, check, new_unittest, unittest_type
66
use stdlib_linalg_constants
7-
use stdlib_linalg, only: lstsq
7+
use stdlib_linalg, only: lstsq,solve_lstsq
88
use stdlib_linalg_state, only: linalg_state_type
99

1010
implicit none (type,external)
@@ -20,6 +20,8 @@ module test_linalg_least_squares
2020
type(unittest_type), allocatable, intent(out) :: tests(:)
2121

2222
allocate(tests(0))
23+
24+
tests = [tests,new_unittest("issue_823",test_issue_823)]
2325

2426
#:for rk,rt,ri in REAL_KINDS_TYPES
2527
#:if rk!="xdp"
@@ -100,6 +102,46 @@ module test_linalg_least_squares
100102

101103
#:endif
102104
#:endfor
105+
106+
! Test issue #823
107+
subroutine test_issue_823(error)
108+
type(error_type), allocatable, intent(out) :: error
109+
110+
! Dimension of the problem.
111+
integer(ilp), parameter :: n = 42
112+
! Data for the least-squares problem.
113+
complex(dp) :: A(n+1, n), b(n+1), x_true(n), x_lstsq(n)
114+
! Internal variables.
115+
real(dp), allocatable :: tmp(:, :, :), tmp_vec(:, :)
116+
! Error handler
117+
type(linalg_state_type) :: state
118+
119+
! Zero-out data.
120+
A = 0.0_dp
121+
b = 0.0_dp
122+
x_lstsq = 0.0_dp
123+
allocate(tmp(n+1, n, 2), tmp_vec(n, 2), source=0.0_dp)
124+
125+
! Generate a random complex least-squares problem of size (n+1, n).
126+
call random_number(tmp)
127+
call random_number(tmp_vec)
128+
129+
A = cmplx(tmp(:, :, 1), tmp(:, :, 2), kind=dp)
130+
x_true = cmplx(tmp_vec(:, 1), tmp_vec(:, 2), kind=dp)
131+
b = matmul(A, x_true)
132+
133+
! Solve the lstsq problem.
134+
call solve_lstsq(A, b, x_lstsq, err=state)
135+
136+
! Check that no segfault occurred
137+
call check(error,state%ok(),'issue 823 returned '//state%print())
138+
if (allocated(error)) return
139+
140+
! Check that least squares are verified
141+
call check(error,all(abs(x_true-x_lstsq)<sqrt(epsilon(0.0_dp))),'issue 823 results')
142+
if (allocated(error)) return
143+
144+
end subroutine test_issue_823
103145

104146
end module test_linalg_least_squares
105147

0 commit comments

Comments
 (0)