Post Reply 
LDLt decomposition?
01-26-2015, 12:30 PM (This post was last modified: 02-06-2015 10:42 PM by salvomic.)
Post: #1
LDLt decomposition?
hi all,
in HP Prime we have for a matrix LU decomposition (it give P, L, U so PA=LU) and Cholesky decomposition (for that A=LLt, Lt is TRAN(L))...
I can't find a factorization LDLt where D is a diagonal of pivot.
Is there a way to get it? a program or something?

Thank you!

Salvo
Visit this user's website Find all posts by this user
Quote this message in a reply
06-01-2015, 05:54 PM
Post: #2
RE: LDLt decomposition?
I wonder still for this question and another one, connected (in my case):
With pivot(Mat, n,m) we get pivot for n,m: what's a simple way to get every pivots (Gauss' reduction); I forgot it Smile

thanks!
Salvo

∫aL√0mic (IT9CLU), HP Prime 50g 41CX 71b 42s 12C 15C - DM42 WP34s :: Prime Soft. Lib
Visit this user's website Find all posts by this user
Quote this message in a reply
06-01-2015, 10:51 PM
Post: #3
RE: LDLt decomposition?
...to get Pivots I start with this code:

Code:

EXPORT gaussJordan(m)
// Gauss-Jordan elimination and pivots
// Salvo Micciché 2015
BEGIN
local temp, temp2, gj, r, c, j, piv;
r:=rowDim(m);
c:=colDim(m);
temp:=MAKEMAT(0,r,c);
gj:=MAKEMAT(0,r,c);
piv:=MAKELIST(0,X,1,r);
gj(1):= m(1);
piv(1):=gj(1,1);
    temp:= pivot(m,1,1);
FOR j FROM 2 TO r DO
    temp:=delrows(temp,1);
    temp:=delcols(temp,1);
    temp:=temp/piv(j-1);
    temp:= pivot(temp,1,1);
    piv(j):=temp(1,1);
    gj(j):= temp(1);
END; // for
RETURN {gj, piv};
END;

but I get wrong results (see image): I should get 0 below pivots at left, not ...at right; I'm trying SWAPROW() but it doesn't run...
Any help?


Attached File(s) Thumbnail(s)
   

∫aL√0mic (IT9CLU), HP Prime 50g 41CX 71b 42s 12C 15C - DM42 WP34s :: Prime Soft. Lib
Visit this user's website Find all posts by this user
Quote this message in a reply
06-02-2015, 11:26 AM
Post: #4
RE: LDLt decomposition?
Is diag() useful for you?

Code:

EXPORT gjp()
BEGIN 
  M5:=[[2,1,1,5],[-8,-2,-12,0],[1,2,0,0]];  // Your matrix
  M6:=diag(M5);  //  Diagonal Matrix
  return M6;  
END;


Attached File(s) Thumbnail(s)
   
Find all posts by this user
Quote this message in a reply
06-02-2015, 12:10 PM (This post was last modified: 06-02-2015 02:23 PM by salvomic.)
Post: #5
RE: LDLt decomposition?
(06-02-2015 11:26 AM)DrD Wrote:  Is diag() useful for you?

Code:

EXPORT gjp()
BEGIN 
  M5:=[[2,1,1,5],[-8,-2,-12,0],[1,2,0,0]];  // Your matrix
  M6:=diag(M5);  //  Diagonal Matrix
  return M6;  
END;

hi Dale,
diag() is also ok to get pivots, but after to have "pivot-ized" the original matrix: [[2,1,1,5],[4,-6,0,-2],[-2,7,2,9]]

My goal is to return a transformed matrix and a list with pivots (and also diag() could be ok).

EDIT:
Something like that, Dale:
Code:

EXPORT gaussJordan(m)
// Gauss-Jordan elimination and pivots
// Salvo Micciché 2015
BEGIN
local temp, k, gj, r, c, j, piv;
r:=rowDim(m);
c:=colDim(m);
temp:=MAKEMAT(0,r,c);
gj:=MAKEMAT(0,r,c);
piv:=MAKELIST(0,X,1,r);
temp:= m;
FOR j FROM 1 TO r DO
      temp:=pivot(temp,1,1);
      // gj(j):= temp(1);
    FOR k FROM 1 TO colDim(temp) DO
        gj(j, c-k+1):= temp(1, colDim(temp)-k+1);
    END; // inner for
      piv(j):=temp(1,1);
    IF (j<r) THEN
    temp:= delrows(temp,1);
    temp:= delcols(temp,1);
    END;
      temp:= temp/piv(j);  
END; // for
piv:= list2mat(piv, r);
RETURN {gj, piv};
END;

this runs, but with "singular" matrices it gives still division by 0 and in some cases fails (we need also swapping of some rows...), however now we have a base to calculate LDLt decomposition also...


Attached File(s) Thumbnail(s)
   

∫aL√0mic (IT9CLU), HP Prime 50g 41CX 71b 42s 12C 15C - DM42 WP34s :: Prime Soft. Lib
Visit this user's website Find all posts by this user
Quote this message in a reply
06-02-2015, 05:54 PM (This post was last modified: 06-02-2015 07:20 PM by salvomic.)
Post: #6
RE: LDLt decomposition?
SWAPROW runs well in command line but not in a program? or am I missing something?

This code doesn't swap row 2 with 3...

Code:
EXPORT gaussJordan(m)
// Gauss-Jordan elimination and pivots
// Salvo Micciché 2015
BEGIN
local k, gj, r, c, j, piv;
r:=rowDim(m);
c:=colDim(m);
gj:=MAKEMAT(0,r,c);
piv:=MAKELIST(0,X,1,r);
M1:= m;
FOR j FROM 1 TO r DO
      M1:=pivot(M1,1,1);
      piv(j):=M1(1,1);
    IF (piv(j)==0) THEN 
PRINT("M1 first " + M1);
    M1:=m;
PRINT("M1 then " + M1);
    SWAPROW(M1, j, (j+1)); 
PRINT("M1 and after " + M1);
    CONTINUE; 
    END; // if
    FOR k FROM 1 TO colDim(M1) DO
        gj(j, c-k+1):= M1(1, colDim(M1)-k+1);
    END; // inner for
    IF (j<r) THEN
    M1:= delrows(M1,1);
    M1:= delcols(M1,1);
    END;
    IF piv(j)==0 THEN RETURN "Pivot is 0, division by 0"; END;
      M1:= M1/piv(j);  
END; // for
piv:= list2mat(piv, r);
RETURN {gj, piv};
END;

See the prints in Terminal (debug, hi) with this matrix
[[1,2,3],[2,4,3],[3,2,1]]
Normally it would give division by zero, so we must swap the row with its next...
The result should be [[1,2,3],[0,-4,-10],[0,0,-3]] with pivots d={1,-4,-3}

I'm using an "old C style" tip, and it works:
Code:

 // M1:= SWAPROW(M1, numpivot, (numpivot+1)); 
temp:= M1(numpivot); M1(numpivot):=M1(numpivot+1); M1(numpivot+1):=temp;

∫aL√0mic (IT9CLU), HP Prime 50g 41CX 71b 42s 12C 15C - DM42 WP34s :: Prime Soft. Lib
Visit this user's website Find all posts by this user
Quote this message in a reply
06-02-2015, 10:05 PM
Post: #7
RE: LDLt decomposition?
I tried both SWAPROW() and SWAPCOL() on a 3x3 matrix and neither worked from the command line or HOME program, for me.

-Dale-
Find all posts by this user
Quote this message in a reply
06-02-2015, 10:15 PM
Post: #8
RE: LDLt decomposition?
(06-02-2015 10:05 PM)DrD Wrote:  I tried both SWAPROW() and SWAPCOL() on a 3x3 matrix and neither worked from the command line or HOME program, for me.

-Dale-

Please, could I suspect a bug? Smile

Strange, however...

Salvo

∫aL√0mic (IT9CLU), HP Prime 50g 41CX 71b 42s 12C 15C - DM42 WP34s :: Prime Soft. Lib
Visit this user's website Find all posts by this user
Quote this message in a reply
06-02-2015, 11:05 PM
Post: #9
RE: LDLt decomposition?
How about CAS.SWAPCOL()? Works here.
Find all posts by this user
Quote this message in a reply
06-03-2015, 11:59 AM
Post: #10
RE: LDLt decomposition?
ALL similar matrix commands have the same requirement: ADDROW(), DELROW(), SWAPROW(), etc., all have the same CAS requirement in order to work on the HC and virtual calc. They seem to work ok on the Android app. I sent a bug report in, but, maybe this is intended behavior? If the command names were in lower case, it would serve as a hint that they need CAS handling.

-Dale-
Find all posts by this user
Quote this message in a reply
06-03-2015, 12:19 PM (This post was last modified: 06-03-2015 02:51 PM by salvomic.)
Post: #11
RE: LDLt decomposition?
(06-03-2015 11:59 AM)DrD Wrote:  ALL similar matrix commands have the same requirement: ADDROW(), DELROW(), SWAPROW(), etc., all have the same CAS requirement in order to work on the HC and virtual calc. They seem to work ok on the Android app. I sent a bug report in, but, maybe this is intended behavior? If the command names were in lower case, it would serve as a hint that they need CAS handling.

-Dale-

I agree.
Here I've issues with them in CAS, however, so I think maybe there could be some bugs with the last FW upgrade...

Salvo

EDIT:
here SWAPROW([[1,2,3],[0,0,-3],[0,-4,-10]],2,3) doesn't work nor in CAS neither in Home.

∫aL√0mic (IT9CLU), HP Prime 50g 41CX 71b 42s 12C 15C - DM42 WP34s :: Prime Soft. Lib
Visit this user's website Find all posts by this user
Quote this message in a reply
06-03-2015, 01:22 PM (This post was last modified: 06-05-2015 06:40 PM by salvomic.)
Post: #12
RE: LDLt decomposition?
however, please, see also this thread (for pivots calculation), that could help us to make an actual program to find Gauss elimination procedure, then pivots to use to find eventually LDLt decomposition.
Please, help Smile

EDIT:
working in progress in this thread in the Forum!

∫aL√0mic (IT9CLU), HP Prime 50g 41CX 71b 42s 12C 15C - DM42 WP34s :: Prime Soft. Lib
Visit this user's website Find all posts by this user
Quote this message in a reply
Post Reply 




User(s) browsing this thread: 1 Guest(s)