lusolve gives wrong answer

[winstel@…]

If one wants to know the solutions of A*x = b for the matrix A[1][6] = 0,1,0,0,1,0 and the right-hand-side matrix b[1][1] = 0 and one first computes list L = ludecomp(A) and then list Q = lusolve(L[1],L[2],L[3],b), Singular tells that the space of the homogeneous linear space is spanned by the vectors v1 = 0,0,0,0,0,-1 v2 = 0,1,0,0,-1,0 v3 = 0,0,0,-1,0,0 and v4 = 0,0,-1,0,0,0. This result is wrong, for instance the vector -1,0,0,0,0,0 is in the solution space but not in the space spanned by v1,v2,v3 and v4.

Edit by gorzel:

> ring r=0,x,dp;
> matrix A[1][6] = 0,1,0,0,1,0;
> matrix b[1][1] = 0;
> list L = ludecomp(A);
> list Q = lusolve(L[1],L[2],L[3],b);
> print(Q[2]);
0,
0,
0,
0,
0,
0 
> print(Q[3]);
0, 0, 0, 0, 
0, 1, 0, 0, 
0, 0, 0, -1,
0, 0, -1,0, 
0, -1,0, 0, 
-1,0, 0, 0  

Maybe permutation for the coloumns is missing?

[gorzel]

Simplification and assignment of the above bug report:

> ring r=0,x,dp;
> matrix D[1][2] = 1,0;   // this is OK...
> list J = ludecomp(D);
> lusolve(J[1],J[2],J[3],matrix(1));
[1]:
   1
[2]:
   _[1,1]=1
   _[2,1]=0
[3]:
   _[1,1]=0
   _[2,1]=-1
> lusolve(J[1],J[2],J[3],matrix(0));
[1]:
   1
[2]:
   _[1,1]=0
   _[2,1]=0
[3]:
   _[1,1]=0
   _[2,1]=-1

> matrix A[1][2] = 0,1;   // But this is wrong
> list L = ludecomp(A);
> lusolve(L[1],L[2],L[3],matrix(0));
[1]:
   1
[2]:
   _[1,1]=0
   _[2,1]=0
[3]:
   _[1,1]=0     // wrong

> lusolve(L[1],L[2],L[3],matrix(1));
[1]:
   1
[2]:
   _[1,1]=0
   _[2,1]=1
[3]:
   _[1,1]=0   // wrong

The kernel (nullspace) is not properly computed. See also

> matrix C = 0;
> C;
C[1,1]=0
> lusolve(K[1],K[2],K[3],matrix(0));   // wrong
[1]:
   1
[2]:
   _[1,1]=0
[3]:
   _[1,1]=0
> lusolve(K[1],K[2],K[3],matrix(1));   // OK
[1]:
   0

I find it cumbersome to call lusolve in the form lusolve(L[1],L[2],L[3],b). Why not to allow to call lusolve directly in the form

lusolve(A,b);

See luinvers

​http://www.singular.uni-kl.de/Manual/latest/sing_262.htm#SEC302

which offers this additional possibility.

[gorzel]

The documentation needs some improvement:

The fourth argument, b, is expected to be an (m x 1) matrix, i.e., a vector.

list L=lusolve(P,L,U,b); fills the list L with either one entry = 0 
(signaling that A*x=b has no solution),  
or with the three entries 1, x, H, where x is any (n x 1) solution of 
the given linear system,

1.) Rephrase the first sentence to :

The fourth argument, b, is expected to be a vector i.e., an (m x 1) matrix.

or avoid completely the word vector since the data-type vector is not accepted.

2.) Change list L=lusolve(P,L,U,b); , as in the example, to list Q=lusolve(P,L,U,b);

3.) Use another font for the solution vector x or rename it to z, since in this phrase

A*x=b has no solution), or with the three entries 1, x, H, where x is any (n x 1) solution

the last x stands for times but looks the same in the browser.

[steidel]

Liebster Frank,

wärst Du so nett, Dir das Ticket doch noch einmal anzuschauen und zu reparieren.

Vielen Dank, Dein Dich liebendes SINGULAR-Team.

[kroeker@…]

lusolve gives still wrong results and crashes for the zero matrix:

ring r= 0,x,dp;

matrix A;

list DL = ludecomp(A);
def P = DL[1];
def L = DL[2];
def U = DL[3];

(P*A)==(L*U);

list RL=lusolve(P,L,U,matrix(0) );

I suggest to disable lusolve until it is fixed.