sgelsx не дает правильного решения

Question

Я хочу найти минимальную норму, наименьшее квадратное решение уравнения Ax = b с использованием драйвера 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. К сожалению, я нашел это объяснение в другой процедуре документация. Интересно, почему автор не использовал то же самое, что более наглядно, объяснение для rcond в sgelsx.