On Sat, 22 Aug 2009, Eugene McDonnell wrote:
> My dim recollection is that I triangularized Matrix M, giving Matrix
> T, and then the inverse is easy to compute. From that I went backwards
> and got the inverse.
>
> I've just finished looking in my files on Matrix problems, but was
> unable to find anything; very vexing.
(copied from an implementation of j, transliterated into newer j explicit
definition)
--------------8<----------8<------------
qr compute the QR decomposition; rinv computes the inverse of a square upper
triangular matrix
pdt =: +/ . *
en =: 1&{@(,&1 1)@$
q2 =: 3 : 0
n =. en y
m =. >.-: n
a0 =. m{."1 y
a1 =. m}."1 y
t0 =. qr a0
q0 =. >@{. t0
r0 =. >@{: t0
c =. (+|:q0) pdt a1
t1 =. qr a1 - q0 pdt c
q1 =. >@{. t1
r1 =. >@{: t1
(q0,.q1);(r0,.c),(-n){."1 r1
)
norm=: (%:@pdt +)@,
qr =: q2`((% ;&,. ,~...@en@[ $ ]) norm) @. (1&>:@en)
r4 =: 3 : 0
n =. #y
m =. >.-: n
ai =. rinv (m,m){.y
di =. rinv (m,m)}.y
b =. (m,m-n){.y
bx =. - ai pdt b pdt di
(ai,.bx),(-n){."1 di
)
rinv =: r4`% @. (1&>:@#)
minv =: (|....@$ ($,) (r...@] pdt +@|:@[)&>/@qr) " 2
mdiv =: (%...@] +/ . * [) " _ 2
--------------8<----------8<------------
NB. test
minv ?. 5 5$100
0.0281292 0.0239569 _0.0254212 _0.011056 _0.0242152
0.010086 _0.00549055 0.012983 _0.0301618 0.00799012
0.0101669 0.0286513 _0.0146031 _0.00447282 _0.028359
_0.0195192 _0.0143284 0.00738001 0.0143682 0.0268751
_0.0255669 _0.0245614 0.0204128 0.0315047 0.0222872
(%. -: minv) ?. 5 5$100
1
(|: ?. 5 5$100) (%. -: mdiv) ?. 5 5$100
1
--
regards,
====================================================
GPG key 1024D/4434BAB3 2008-08-24
gpg --keyserver subkeys.pgp.net --recv-keys 4434BAB3
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
