Skip to content

Commit 465d230

Browse files
perazzperazz
committed
base implementation
1 parent 892bf1e commit 465d230

File tree

2 files changed

+211
-0
lines changed

2 files changed

+211
-0
lines changed

‎src/CMakeLists.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -21,6 +21,7 @@ set(fppFiles
2121
stdlib_kinds.fypp
2222
stdlib_linalg.fypp
2323
stdlib_linalg_diag.fypp
24+
stdlib_linalg_least_squares.fypp
2425
stdlib_linalg_outer_product.fypp
2526
stdlib_linalg_kronecker.fypp
2627
stdlib_linalg_cross_product.fypp

‎src/stdlib_linalg_least_squares.fypp

Lines changed: 210 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,210 @@
1+
#:include "common.fypp"
2+
#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
3+
#:set RHS_SUFFIX = ["one","many"]
4+
#:set RHS_SYMBOL = [ranksuffix(r) for r in [1,2]]
5+
#:set RHS_EMPTY = [emptyranksuffix(r) for r in [1,2]]
6+
#:set ALL_RHS = list(zip(RHS_SYMBOL,RHS_SUFFIX,RHS_EMPTY))
7+
module stdlib_linalg_least_squares
8+
use stdlib_linalg_constants
9+
use stdlib_linalg_lapack, only: gelsd, stdlib_ilaenv
10+
use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, &
11+
LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR
12+
implicit none(type,external)
13+
private
14+
15+
!> Compute a least squares solution to system Ax=b, i.e. such that the 2-norm abs(b-Ax) is minimized.
16+
public :: lstsq
17+
18+
! NumPy: lstsq(a, b, rcond='warn')
19+
! Scipy: lstsq(a, b, cond=None, overwrite_a=False, overwrite_b=False, check_finite=True, lapack_driver=None)
20+
! IMSL: Result = IMSL_QRSOL(B, [A] [, AUXQR] [, BASIS] [, /DOUBLE] [, QR] [, PIVOT] [, RESIDUAL] [, TOLERANCE])
21+
22+
interface lstsq
23+
#:for nd,ndsuf,nde in ALL_RHS
24+
#:for rk,rt,ri in RC_KINDS_TYPES
25+
module procedure stdlib_linalg_${ri}$_lstsq_${ndsuf}$
26+
#:endfor
27+
#:endfor
28+
end interface lstsq
29+
30+
31+
contains
32+
33+
#:for rk,rt,ri in RC_KINDS_TYPES
34+
! Workspace needed by gesv
35+
subroutine ${ri}$gesv_space(m,n,nrhs,lrwork,liwork,lcwork)
36+
integer(ilp), intent(in) :: m,n,nrhs
37+
integer(ilp), intent(out) :: lrwork,liwork,lcwork
38+
39+
integer(ilp) :: smlsiz,mnmin,nlvl
40+
41+
mnmin = min(m,n)
42+
43+
! Maximum size of the subproblems at the bottom of the computation (~25)
44+
smlsiz = stdlib_ilaenv(9,'${ri}$gelsd',' ',0,0,0,0)
45+
46+
! The exact minimum amount of workspace needed depends on M, N and NRHS. As long as LWORK is at least
47+
nlvl = max(0, ilog2(mnmin/(smlsiz+1))+1)
48+
49+
! Real space
50+
#:if rt.startswith('complex')
51+
lrwork = 10*mnmin+2*mnmin*smlsiz+8*mnmin*nlvl+3*smlsiz*nrhs+max((smlsiz+1)**2,n*(1+nrhs)+2*nrhs)
52+
#:else
53+
lrwork = 12*mnmin+2*mnmin*smlsiz+8*mnmin*nlvl+mnmin*nrhs+(smlsiz+1)**2
54+
#:endif
55+
lrwork = max(1,lrwork)
56+
57+
! Complex space
58+
lcwork = 2*mnmin + nrhs*mnmin
59+
60+
! Integer space
61+
liwork = max(1, 3*mnmin*nlvl+11*mnmin)
62+
63+
! For good performance, the workspace should generally be larger.
64+
lrwork = ceiling(1.25*lrwork,kind=ilp)
65+
lcwork = ceiling(1.25*lcwork,kind=ilp)
66+
liwork = ceiling(1.25*liwork,kind=ilp)
67+
68+
end subroutine ${ri}$gesv_space
69+
70+
#:endfor
71+
72+
#:for nd,ndsuf,nde in ALL_RHS
73+
#:for rk,rt,ri in RC_KINDS_TYPES
74+
75+
! Compute the least-squares solution to a real system of linear equations Ax = B
76+
function stdlib_linalg_${ri}$_lstsq_${ndsuf}$(a,b,cond,overwrite_a,rank,err) result(x)
77+
!> Input matrix a[n,n]
78+
${rt}$, intent(inout), target :: a(:,:)
79+
!> Right hand side vector or array, b[n] or b[n,nrhs]
80+
${rt}$, intent(in) :: b${nd}$
81+
!> [optional] cutoff for rank evaluation: singular values s(i)<=cond*maxval(s) are considered 0.
82+
real(${rk}$), optional, intent(in) :: cond
83+
!> [optional] Can A,b data be overwritten and destroyed?
84+
logical(lk), optional, intent(in) :: overwrite_a
85+
!> [optional] Return rank of A
86+
integer(ilp), optional, intent(out) :: rank
87+
!> [optional] state return flag. On error if not requested, the code will stop
88+
type(linalg_state_type), optional, intent(out) :: err
89+
!> Result array/matrix x[n] or x[n,nrhs]
90+
${rt}$, allocatable, target :: x${nd}$
91+
92+
!> Local variables
93+
type(linalg_state_type) :: err0
94+
integer(ilp) :: m,n,lda,ldb,nrhs,info,mnmin,mnmax,arank,lrwork,liwork,lcwork
95+
integer(ilp), allocatable :: iwork(:)
96+
logical(lk) :: copy_a
97+
real(${rk}$) :: acond,rcond
98+
real(${rk}$), allocatable :: singular(:),rwork(:)
99+
${rt}$, pointer :: xmat(:,:),amat(:,:)
100+
${rt}$, allocatable :: cwork(:)
101+
character(*), parameter :: this = 'lstsq'
102+
103+
!> Problem sizes
104+
m = size(a,1,kind=ilp)
105+
lda = size(a,1,kind=ilp)
106+
n = size(a,2,kind=ilp)
107+
ldb = size(b,1,kind=ilp)
108+
nrhs = size(b ,kind=ilp)/ldb
109+
mnmin = min(m,n)
110+
mnmax = max(m,n)
111+
112+
if (lda<1 .or. n<1 .or. ldb<1 .or. ldb/=m) then
113+
err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid sizes: a=',[lda,n], &
114+
'b=',[ldb,nrhs])
115+
allocate(x${nde}$)
116+
arank = 0
117+
goto 1
118+
end if
119+
120+
! Can A be overwritten? By default, do not overwrite
121+
if (present(overwrite_a)) then
122+
copy_a = .not.overwrite_a
123+
else
124+
copy_a = .true._lk
125+
endif
126+
127+
! Initialize a matrix temporary
128+
if (copy_a) then
129+
allocate(amat(lda,n),source=a)
130+
else
131+
amat => a
132+
endif
133+
134+
! Initialize solution with the rhs
135+
allocate(x,source=b)
136+
xmat(1:n,1:nrhs) => x
137+
138+
! Singular values array (in degreasing order)
139+
allocate(singular(mnmin))
140+
141+
! rcond is used to determine the effective rank of A.
142+
! Singular values S(i) <= RCOND*maxval(S) are treated as zero.
143+
! Use same default value as NumPy
144+
if (present(cond)) then
145+
rcond = cond
146+
else
147+
rcond = epsilon(0.0_${rk}$)*mnmax
148+
endif
149+
if (rcond<0) rcond = epsilon(0.0_${rk}$)*mnmax
150+
151+
! Allocate working space
152+
call ${ri}$gesv_space(m,n,nrhs,lrwork,liwork,lcwork)
153+
#:if rt.startswith('complex')
154+
allocate(rwork(lrwork),cwork(lcwork),iwork(liwork))
155+
#:else
156+
allocate(rwork(lrwork),iwork(liwork))
157+
#:endif
158+
159+
! Solve system using singular value decomposition
160+
#:if rt.startswith('complex')
161+
call gelsd(m,n,nrhs,amat,lda,xmat,ldb,singular,rcond,arank,cwork,lrwork,rwork,iwork,info)
162+
#:else
163+
call gelsd(m,n,nrhs,amat,lda,xmat,ldb,singular,rcond,arank,rwork,lrwork,iwork,info)
164+
#:endif
165+
166+
! The condition number of A in the 2-norm = S(1)/S(min(m,n)).
167+
acond = singular(1)/singular(mnmin)
168+
169+
! Process output
170+
select case (info)
171+
case (0)
172+
! Success
173+
case (:-1)
174+
err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid problem size a=[',lda,',',n,'], b[',ldb,',',nrhs,']')
175+
case (1:)
176+
err0 = linalg_state_type(this,LINALG_ERROR,'SVD did not converge.')
177+
case default
178+
err0 = linalg_state_type(this,LINALG_INTERNAL_ERROR,'catastrophic error')
179+
end select
180+
181+
if (.not.copy_a) deallocate(amat)
182+
183+
! Process output and return
184+
1 call linalg_error_handling(err0,err)
185+
if (present(rank)) rank = arank
186+
187+
end function stdlib_linalg_${ri}$_lstsq_${ndsuf}$
188+
189+
#:endfor
190+
#:endfor
191+
192+
! Simple integer log2 implementation
193+
elemental integer(ilp) function ilog2(x)
194+
integer(ilp), intent(in) :: x
195+
196+
integer(ilp) :: remndr
197+
198+
if (x>0) then
199+
remndr = x
200+
ilog2 = -1_ilp
201+
do while (remndr>0)
202+
ilog2 = ilog2 + 1_ilp
203+
remndr = shiftr(remndr,1)
204+
end do
205+
else
206+
ilog2 = -huge(0_ilp)
207+
endif
208+
end function ilog2
209+
210+
end module stdlib_linalg_least_squares

0 commit comments

Comments
 (0)