On Fri, Apr 30, 2021 at 4:11 AM Imre Patyi <[email protected]> wrote:
> Can anyone tell me which routine should I call in the J addon for LAPACK in
> order to solve a linear least squares problem such as I mentioned?
xGELS i.e. either DGELS (real) or ZGELS (complex).
A simplified J code for your problem (X matrix is not singular and needs no
scaling) can be expressed with mt addon as:
x=:_0.5+(i.1000)%999
b=:^x
X=:(1 o. (x */ >:i.10)) ,. (2 o. (x */ i.11))
load 'math/mt'
NB. number of columns
n=. c_mt_ X
NB. factorize (Q*R = A), Q will be in factored form
QfR=. geqrf_mt_ X
NB. 1) multiply (Q * b) by DORMQR (or ZUNMQR) without forming Q explicitly
NB. 2) solve (R * c2 = head(Q*b)) by xTRSM for c2, by feeding
NB. a matrix which contains R as a system matrix,
NB. here head(noun) means (n {. noun) as xGELS approaches
c2=. QfR ([ trsmlunn_mt_&(n&{.) unmqrlc_mt_) b
>./|(X (+/ . *) c2)-b
7.54952e_15
>./|(((1 o. (y */ (>:i.10))) (+/ . *) ((i.10){c2)) + ((2 o. (y */ (i.
11))) (+/ . *) ((10+i.11){c2)))-(,.^y)
7.99361e_15
Q-factorization (LQ, QL, QR, RQ) algorithms in mt addon are very close to
LAPACK so they are numerically stable.
--
Regards
Igor
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
