sgelsx 没有产生正确的解决方案

Question

我想求一个方程的最小范数，最小二乘解 Ax = b 在 Fortran 中使用 LAPACK 的驱动程序 sgelsx。在它的documentation里面说调用这个子程序后最小二乘解存储为b，但是b和x的维数一般不一样，怎么能它被确定为解决方案？我已尝试计算提供的示例，以便将我的代码结果与解决方案进行比较，结果明显不同。

我的程序是

program test_svd
implicit none

integer:: m,n,nrhs,LDA,LDB,RANK,INFO,i,nwork
integer, allocatable, dimension(:)::JPVT
real, allocatable, dimension(:):: work
real, allocatable, dimension(:,:):: a

real, allocatable, dimension(:,:):: b
real:: RCOND

m = 5
n = 2
nrhs = 1
LDA = 10
LDB = 10
RCOND = 500.0
nwork = maxval( [minval([m,n])+3*n, 2*minval([m,n])+nrhs ])

allocate(b(LDB,nrhs))
allocate(a(LDA,n))
allocate(JPVT(n))
allocate(work(nwork ))

JPVT = (/ (1.0 , i = 1,n) /)


a(1,1) = 1.0
a(2,1) = 1.0
a(3,1) = 1.0
a(4,1) = 1.0
a(5,1) = 1.0

a(1,2) = -2.0
a(2,2) = -1.0
a(3,2) = 0.0
a(4,2) = 2.0
a(5,2) = 3.0

b(1,1) = 4.0
b(2,1) = 2.0
b(3,1) = 1.0
b(4,1) = 1.0
b(5,1) = 1.0

!print *, a, size(a)
!print *, b, size(b)

call sgelsx(m,n,nrhs,a,LDA,b,LDB,JPVT,RCOND,RANK,work,INFO)

print *, b


end

我得到的b（解决方案）是[1.55172420 0.620689809 -1.36653137 -0.703225672 -0.371572882]，它应该是[2, -0.5]。

Answer 1

问题已解决。事实证明 rcond 是一个用于设置阈值的参数，低于该阈值 A 的奇异值将设置为零，因此 rcond 应该像 0.001 或 0.1 取决于 A 的条件数。不幸的是，我在另一个例程的文档中找到了这个解释。我想知道为什么作者没有在 sgelsx.

中对 rcond 使用相同的更具描述性的解释

sgelsx 没有产生正确的解决方案

sgelsx does not produce the correct solution

fortran

lapack

least-squares