I want to find the minimal norm, least square solution of an equation
Ax = b
using LAPACK's driver sgelsx in Fortran. In its documentation, it says that the least square solution is stored as b after calling this subroutine but but the dimensions of b and x are generally different, so how can it be identified as the solution?
I have tried calculating a provided example so I can compare the result of my code with the solution, and it's clearly different.
My program is
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
What I get as b (the solution) is [1.55172420 0.620689809 -1.36653137 -0.703225672 -0.371572882], which is supposed to be [2, -0.5].
Problem was solved. It turns out that rcond is a parameter used to set a threshold below which the singular values of A is set to zero, so rcond should be like 0.001 or 0.1 depending on the condition number of A. Unfortunately, I found this explanation in another routine's documentation. I wonder why the author didn't use the same, which is more descriptive, explanation for rcond in sgelsx.