A solution to this problem is to use rationals,
where there is no such thing as "near singular",
only singular or non-singular.
LUcheck x: A
|index error: LU
| p;(=1);p {"1 A[p=.C.(n-1);~.0,(0~:,A)i.1
LUcheck x: |:A
|index error: LU
| p;(=1);p {"1 A[p=.C.(n-1);~.0,(0~:,A)i.1
----- Original Message -----
From: John Randall <[email protected]>
Date: Sunday, March 13, 2011 10:31
Subject: [Jprogramming] LU decomposition essay
To: 'Programming forum' <[email protected]>
> LU may work on singular matrices even when it should not.
> The same is
> true of the getrf family from LAPACK. This is to be
> expected, since a
> small perturbation of a singular matrix is nonsingular. A
> perturbednonsingular matrix remains nonsingular.
>
> In the code below, LUcheck is from the LU decomposition
> essay. The matrix
> A is singular. The matrix B is the transpose of A, and so
> also singular.
> The way both routines do the calculation, A is figured out to be
> nonsingular while B is correctly identified as singular.
>
> require '~addons/math/lapack/lapack.ijs'
> require '~addons/math/lapack/getrf.ijs'
> getrf_z_=:getrf_jlapack_
> ]A=:>:i. 3 3
> 1 2 3
> 4 5 6
> 7 8 9
> B=:|: A
> LUcheck A
> +-----+-----+-----------------+
> |0 1 2|1 0 0|1
> 2 3|
> | |4 1 0|0
> _3 _6|
> | |7 2 1|0 0 _3.55271e_15|
> +-----+-----+-----------------+
> getrf A
> +--------------+-----------------------+-----+
> | 1 0
> 0|7 8 9|3 3 3|
> |0.142857 1 0|0
> 0.857143 1.71429| |
> |0.571429 0.5 1|0 0
> _1.58619e_16| |
> +--------------+-----------------------+-----+
> LUcheck B
> |index error: LU
> | p;(=1);p {"1 A[p=.C.(n-
> 1);~.0,(0~:,A)i.1 getrf B
> info result: 3
> |attention interrupt: error
> | error''
>
>
> Best wishes,
>
> John
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
