NB. I am curious as to why I can not invert this matrix using floats:
(J602/32 bit/Windows)

diag=. * =...@i.@#

  diag 1 1e8 1e_6
1         0    0
0 100000000    0
0         0 1e_6

  %. diag 1 1e8 1e_6
|domain error
|       %.diag 1 100000000 1e_6

     %. diag 1 1e8 1e_5
1    0      0
0 1e_8      0
0    0 100000

NB. Is there a hardware overflow/underflow that causes the domain error?

NB. My workaround is relying on the nifty fact that x: has an inverse...

  %.&.x: diag 1 1e8 1e_6
1    0       0
0 1e_8       0
0    0 1000000

NB. But this incurs a 2 orders of magnitude penalty:

A=. diag 1e_5,49#1e8

  '%.A' 6!:2~ 10
0.00191876342
  '%.&.x: A' 6!:2~ 10
0.319551787

NB. I generally know when I am about to deal with a diagonal matrix,
NB. and I wouldn't mind just using reciprocal on the diagonal but
NB. I am having problems with some other matrices as well


B=. ".;._2 ] 0 : 0
_1.54198508e_6 _2.11460184e_5  0.00111041321 0.000412974295
0.000309255491 0.000388701673 0.000316571441 0.000276780722
0.000228708457 0.000407042731
_8.53914448e_6 _2.51880232e_5 0.000412974295  0.00100957149
0.000329506981  0.00037242818 0.000295428584 0.000279403683
0.00019269837 0.000399173819
_1.93242844e_5 _1.64241482e_5 0.000309256491 0.000329506981
1e9 0.000329938591 0.000301756474 0.000300975448 0.000128772251
0.000328033069
_7.22436598e_7 _1.89792044e_5 0.000388701673  0.00037242918
0.000329938591            1e9 0.000295252768 0.000278363898
0.000174674958 0.000369625237
2.21070087e_6  _5.1700877e_6 0.000316571441 0.000295428584
0.000301757474 0.000295252768 0.000961677712 0.000267397094
0.000122513186 0.000291930669
8.79414432e_6   5.2053919e_6 0.000276780722 0.000279403683
0.000300975448 0.000278364898 0.000267397094            1e9
0.000112259029 0.000275270547
_8.22583816e_6 _1.38022985e_5 0.000228708457  0.00019269837
0.000128772251 0.000174674958 0.000122514186 0.000112259029
0.000897602311 0.000189206168
_8.46457843e_6 _2.38916908e_5 0.000407042731 0.000399173819
0.000328033069 0.000369625237 0.000291930669 0.000275271547
0.000189206168            1e9
_5.29236623e_6 _2.17598177e_5  0.00036309251 0.000356138358
0.000321507113 0.000340217879 0.000289955034 0.000275345859
0.00016401428 0.000353603807
6.18618386e_5  1.23350153e_5  2.89117576e_5  2.02892574e_5
2.94299546e_5  2.17703052e_5  1.89295016e_5  1.60711228e_5
1.80029803e_5  2.05957395e_5
)

  %.B
|domain error
|       %.B

  %.&.x: B NB. This will take a while
    _4241.7248      _4692.4223  _3.13499965e_9  _2.81091047e_9
_3676.29  _2.45470257e_9     _1622.39661  _2.85952075e_9
21171.6754      17138.5129
    18688.4472      21605.5095    1.3534087e_8   1.27574104e_8
16620.098   1.12221132e_8      6711.78364   1.31282423e_8
_100061.924     _4813.73493
    1361.57405      _107.22163  2.22389472e_11 _5.56934835e_11
_45.6583065  2.26216155e_12     _100.694377 _6.08380988e_11
_1177.24217     _93.3339464
   _22.0803779      1630.88601  4.11708211e_11  8.96697828e_12
60.0458843  4.01543064e_11     _19.0608319 _9.16309002e_12
_1832.78146      63.0928301
_8.44255984e_11 _1.38991708e_10            1e_9 _9.67797959e_23
_1.52048503e_10 _1.14362704e_22 _1.97604825e_11 _8.53398978e_23
_3.11715506e_10  2.67226229e_10
_6.4842665e_11 _6.07890822e_11 _3.22032048e_23            1e_9
_3.59002293e_11 _2.29967399e_23 _2.11895049e_11 _3.67753797e_23
_7.53319406e_10 _6.43113017e_11
   _324.920861     _340.415205 _2.61050606e_10 _2.18675472e_10
1125.70431 _2.31191892e_10     _79.5395589 _2.04882797e_10
454.818729     _66.9829136
_3.38244489e_10 _4.07618444e_10 _2.98114289e_22 _2.50681687e_22
_3.66851989e_10            1e_9 _1.05745794e_10 _2.45434275e_22
1.13105235e_9 _1.11043083e_10
   _49.3419484      12.8850982  5.70451686e_11  1.79199957e_11
66.9694538  4.73585933e_11      1243.02247  1.03082903e_11
_713.269418       102.42037
_3.06205547e_11 _3.39507384e_11  8.71156744e_24 _9.08455674e_24
1.92814882e_11  1.36478111e_23 _1.67490261e_11            1e_9
_9.98473569e_10   4.3043719e_11

  6!:2 '%.&.x: B'
74.0365499
NB. Yikes!!!


NB. Is there any other faster ways to deal with this?
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm

Reply via email to