Pivots (Gaussian reduction)
06-03-2015, 01:20 PM (This post was last modified: 06-06-2015 01:36 PM by salvomic.)
Post: #1 salvomic Senior Member Posts: 1,394 Joined: Jan 2015
Pivots (Gaussian reduction)
hi everybody,
please, help to debug and improve this program to calculate pivots of a matrix (Gauss-Jordan elimination).
It is important to me to create also another program to calculate LDLt factorization (decomposition) [this is the thread], where for a symmetric matrix L is a lower triangular (Lt its transpose) and D diagonal with pivots...

The program for pivots use the function pivot(), finding recursively the minors of the matrix and applying this function to them, saving pivots and rows in a matrix, and at end giving upper triangular matrix with pivots in diagonal and a list with pivots.
The program can, if a pivot is 0, to swap two rows (this function could be improved) and retry for a new calculation...
But my algorithm is not so good...

I wasn't able to find pivots directly with pivot() function, without using it more times...
Or help to find pivots and reduced form with them in diagonal with ref(), RREF(), pivot(), lu(), diag()...

Salvo

The code:
Code:
 calcPivot(); zeropivot:=0; numpivot:=1; gj:=[0,0]; EXPORT gaussJordan(m)   // Gauss-Jordan elimination and pivots   // Salvo Micciché 2015 BEGIN   local r, c, piv, temp;   r:=rowDim(m);   c:=colDim(m);   gj:=MAKEMAT(0,r,c);   M1:= m;   L1:={0,0};   calcPivot(c,r);   IF (zeropivot==1) THEN   M1:=m;   // M1:= CAS.SWAPROW(M1, numpivot, (numpivot+1));    temp:= M1(numpivot); M1(numpivot):=M1(numpivot+1); M1(numpivot+1):=temp;   calcPivot(c,r);    END; // if   piv:= list2mat(L1, r);   M1:= gj;   // Return pivot, matrix reduced and RREF form   RETURN {piv, M1, ref(m)}; END; calcPivot(c,r) BEGIN   local j, k;   zeropivot:=0;   FOR j FROM 1 TO r DO     M1:=CAS.pivot(M1,1,1);     L1(j):= M1(1,1);   IF (L1(j)==0) THEN zeropivot:=1; numpivot:=j; RETURN {zeropivot, numpivot}; END;   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   M1:= exact(M1./L1(j));   IF L1(j)==0 THEN RETURN "Pivot is 0, division by 0"; END;     END; // for END;

this seems to be good with matrices like [[2,1,1,5],[4,-6,0,-2],[-2,7,2,9]] (or [[2,1,1],[4,-6,0],[-2,7,2]], pivots 2,-8,1), [[2,-1,0],[-1,2,-1],[0,-1,1]] (pivots 2, 3/2 and 4/3)
but has little wrong results with [[1,2,3],[2,4,3],[3,2,-1]]: the pivots should be (1,-4,-3) but I get (1,-4,3/4)
and isn't good for echelon matrices (not RREF reduced), like [[1,3,3,2],[2,6,9,7],[-1,-3,3,4]] that gives [[1,3,3,2],[0,0,0,0],[0,0,0,0]] (pivots 1,0,0?)

EDIT: please, help whit last rows of program: why difference in the last pivot? about what am I wrong? I need a more general solution...

∫aL√0mic (IT9CLU) :: HP Prime 50g 41CX 71b 42s 39s 35s 12C 15C - DM42, DM41X - WP34s Prime Soft. Lib
06-05-2015, 05:06 PM (This post was last modified: 06-05-2015 10:47 PM by salvomic.)
Post: #2 salvomic Senior Member Posts: 1,394 Joined: Jan 2015
RE: Pivots (Gaussian elimination)
A version a bit more "advanced". Now it can handle swapping for rows if there is problem in finding solution (*if* the solution there is, the matrix isn't singular, and so on...).

But I've still not found a real *general* formula: help very appreciated!

When it runs well, the program now returns pivots(), LDU() and LDLt() factorization, and a routine to calculate permutation matrix also.
In this version there isn't still a menu...

The code (sorry for some //comments for debugging):
Code:
 calcPivot(); zeropivot:=0; numpivot:=1; gj:=[0,0]; EXPORT gaussJordan() // Salvo Micciché 2015 // Credits: Dale (DrD), Claudio L. BEGIN RECT_P(); TEXTOUT_P("Gauss Jordan utilities", 20, 10, 5, RGB(0,0,255)); TEXTOUT_P("Pivots and upper matrix", 20, 40); TEXTOUT_P("LDLt factorization", 20, 60); TEXTOUT_P("LDU (LDV) factorization", 20, 80); TEXTOUT_P("Permutation Matrix from a list", 20, 100); TEXTOUT_P("Program created by Salvo Micciché", 20, 150, 3, RGB(0,255,0)); TEXTOUT_P("Credits: Dale (DrD), Claudio L.", 20, 180, 2, RGB(0,255,0)); WAIT; END; EXPORT Pivots(m)   // Gauss-Jordan elimination and pivots   // Salvo Micciché 2015 BEGIN   local r, c, piv, temp;   r:=rowDim(m);   c:=colDim(m);   gj:=MAKEMAT(0,r,c);   M1:= m;   L1:={0,0}; // list of pivots   calcPivot(c,r);   IF (zeropivot==1) THEN   M1:=m;   // M1:= CAS.SWAPROW(M1, numpivot, (numpivot+1));    if (numpivot==r) THEN numpivot:=1; END;  // if pivot is 0 swap rows (if last swap w/ 1)   temp:= M1(numpivot); M1(numpivot):=M1(numpivot+1); M1(numpivot+1):=temp;   calcPivot(c,r);    END; // if   piv:= list2mat(L1, r);   M1:= gj;   // Return pivot, matrix reduced and RREF form   RETURN {piv, M1, ref(m)}; END; calcPivot(c,r) BEGIN   local j, k;   zeropivot:=0;   FOR j FROM 1 TO r DO     M1:=CAS.pivot(M1,1,1);     L1(j):= M1(1,1); // current pivot     IF (L1(j)==0) THEN zeropivot:=1; numpivot:=j; // if pivot is 0 swap rows     PRINT("Swap row " + j + " w/ " + IFTE(j<r, j+1, 1));      RETURN {zeropivot, numpivot}; END;     FOR k FROM 1 TO colDim(M1) DO         gj(j, c-k+1):= M1(1, colDim(M1)-k+1); // make [gj] (upper matrix) with correct rows     END; // inner for     IF (j<r) THEN  // take minor of the matrix     M1:= delrows(M1,1);     M1:= delcols(M1,1);   END; // if   // division or not?   M1:= M1 ./ L1(j); // M2, M3, M4 (also IF (j<r))  //  M1(1):= M1(1) ./ L1(j); // M3, M4, M6   // IF (j>1) THEN M1:= M1 ./ L1(j-1); END; // M3, M4, M7   // without div  M3, M4, M7   IF L1(j)==0 THEN RETURN "Pivot is 0, division by 0"; END;    END; // for END; EXPORT LDLt(m) // Factorizazion LDLt, matrix m must be simmetric: A=At BEGIN Pivots(m); local Lmat, LmatT, Dmat, r, j; r:= rowDim(m); Dmat:= diag(L1); Lmat:= M1; FOR j FROM 1 TO r DO Lmat(j):= Lmat(j)/L1(j); // divide Lt (U) by pivots to get 1 in diagonal END; // for LmatT:= TRN(Lmat); // Return L (lower) matrix, D matrix (pivots in diagonal), Lt (upper) tranpose of L RETURN {LmatT, Dmat, Lmat}; END; EXPORT LDU(m) // LDU factorization BEGIN local Lmat, Umat, Dmat,r, j, temp; local permu; Pivots(m); r:= rowDim(m); Dmat:= diag(L1); temp:= lu(m); Lmat:= temp(2); Umat:= temp(3); permu:= permMatrix(temp(1)); FOR j FROM 1 TO r DO Umat(j):= Umat(j)/L1(j); // divide U by pivots to get 1 in diagonal END; // for PRINT("Permutation matrix is " + temp(1)); // Return L (lower), D matrix (pivots in diagonal), U (upper) and permutation matrix RETURN {Lmat, Dmat, Umat, permu}; END; EXPORT permMatrix(lst) // given a permutation list (like {1,3,2}) and the matrix dimension, returns a permutation matrix BEGIN local j, k, mat, dimen; dimen:= length(lst); mat:= MAKEMAT(0, dimen, dimen); FOR j FROM 1 TO dimen DO mat(j,lst(j)):= 1; END; // for RETURN mat; END;

The examples (see at the end of the program are):
M2: [[2,1,1], [4,-6,0], [-2,7,2]] with pivots: {2,-8,1} -> U [[2,1,1], [0,-8,2], [0,0,1]]
M3: [[1,1,1], [1,2,-2], [2,1,-3]] with pivots: {1,1,8} -> [[1,1,1], [0,1,-3], [0,0,-8]]
M4: [[0,3,-1], [1,1,2], [1,4,1]] with pivots: {1,3,0} -> [[1,1,2], [0,3,-1], [0,0,0]] (this one will swap two times; try also another version: [[1,1,2], [0,3,-1], [1,4,1]]: swap two times, with inverted pivots: without the function for swap in the program it'll give "error: invalid input"...): it's a singular matrix (determinant 0, as the product of pivots)
M6: [[2,-1,0], [-1,2,-1], [0,-1,2]] with pivots: {2, 3/2, 4/3} -> U [[2, -1, 0], [0, 3/3, -1], [0,0, 4/3]]
M7: [[1,2,3], [2,4,3], [3,2,-1]] with pivots: {1, -4, -3} -> [[1,2,3], [0,-4,-10], [0, 0, -3]]
M8: [[1,1,1], [1,1,3], [2,5,8]] with pivots: {1,3,2} -> U [[1,1,1], [0,3,6], [0,0,2]] (with swapping row 2 and 3)

And now the problem of the "non general" formula (please, let me understand where I'm wrong):
Code:
   // M1:= M1./L1(j); // M2, M3, M4 (also IF (j<r))   M1(1):= M1(1)./L1(j); // M3, M4, M6   // IF (j>1) THEN M1:= M1./L1(j-1); END; // M3, M4, M7   // without div  M3, M4, M7
M2 needs division for current pivot in every lines, M6 only in the first, the others (M3, M4, M7, M8)) don't need division (apparently)...

But I'm searching for *one* formula, not more ;-)

Thank you for your patience and effort!

Salvo

EDIT: thanks Dale (DrD) and Claudio L. for some tips!

∫aL√0mic (IT9CLU) :: HP Prime 50g 41CX 71b 42s 39s 35s 12C 15C - DM42, DM41X - WP34s Prime Soft. Lib
06-05-2015, 07:03 PM
Post: #3
 DrD Senior Member Posts: 1,132 Joined: Feb 2014
RE: Pivots (Gaussian elimination)
(06-05-2015 05:06 PM)salvomic Wrote:  M1(1):= M1(1)./L1(j); // M3, M4, M6

The examples (see at the end of the program are):
M2: [[2,1,1], [4,-6,0], [-2,7,2]] with pivots: {2,-8,1} -> U [[2,1,1], [0,-8,2], [0,0,1]]
M3: [[1,1,1], [1,2,-2], [2,1,-3]] with pivots: {1,1,8} -> [[1,1,1], [0,1,-3], [0,0,-8]]
M4: [[0,3,-1], [1,1,2], [1,4,1]] with pivots: {1,3,0} -> [[1,1,2], [0,3,-1], [0,0,0]] (this one will swap two times; try also another version: [[1,1,2], [0,3,-1], [1,4,1]]: swap two times, with inverted pivots: without the function for swap in the program it'll give "error: invalid input"...)
M6: [[2,-1,0], [-1,2,-1], [0,-1,2]] with pivots: {2, 3/2, 4/3} -> U [[2, -1, 0], [0, 3/3, -1], [0,0, 4/3]]
M7: [[1,2,3], [2,4,3], [3,2,-1]] with pivots: {1, -4, -3} -> [[1,2,3], [0,-4,-10], [0, 0, -3]]

And now the problem of the "non general" formula (please, let me understand where I'm wrong):
Code:
   // M1:= M1./L1(j); // M2, M3, M4 (also IF (j<r))   M1(1):= M1(1)./L1(j); // M3, M4, M6   // IF (j>1) THEN M1:= M1./L1(j-1); END; // M3, M4, M7   // without div  M3, M4, M7
M2 needs division for current pivot in every lines, M6 only in the first, the others don't need division (apparently)...

But I'm searching for *one* formula, not more ;-)

Thank you for your patience and effort!
Salvo

In another post Han said to "put spaces around the ./ operator."
06-05-2015, 07:29 PM (This post was last modified: 06-05-2015 10:03 PM by salvomic.)
Post: #4 salvomic Senior Member Posts: 1,394 Joined: Jan 2015
RE: Pivots (Gaussian elimination)
(06-05-2015 07:03 PM)DrD Wrote:  In another post Han said to "put spaces around the ./ operator."

ok, I'm trying it, but doesn't change...
The principal problem is: why M3, M4, M7 seem don't need division, M2 need recursive division, and M6 only for the first row that after will be treated by pivot()?
If I execute manually pivot() -> finding a pivot -> finding minor (reducing first row and column) at the end I need to decide where divide or not also: but how to tell it to the calculator? I'm stuck only in this point...

Salvo

∫aL√0mic (IT9CLU) :: HP Prime 50g 41CX 71b 42s 39s 35s 12C 15C - DM42, DM41X - WP34s Prime Soft. Lib
06-06-2015, 09:41 PM (This post was last modified: 06-07-2015 09:45 AM by salvomic.)
Post: #5 salvomic Senior Member Posts: 1,394 Joined: Jan 2015
RE: Pivots (Gaussian reduction)
ok, this should works much better (with a trick, for now, maybe not the best of all).

Test it, if you like it. Maybe the program has issues with large matrix (col > 3) if c>r (and tell me if you find other bugs).
For example still it fails (after the second row it gives all 0) with a matrix like:
[[0,2,1,1,0], [1,-1,2,0,3], [1,1,3,1,2], [3,1,8,2,6], [4,2,-6,-8,0]] (require row swap, but manually it's solvable); the built in function pivot() supply "too many" 0...
See here for the correct reduction.
For that ref() (function echelon in my program) gives
[[1,-1,2,0,3], [0,1,1/2,1/2,0], [0,0,1,11/17, 12/17], [0,0,0,0,1], [0,0,0,0,0]], that multiplied by {1,2,-17,-1,0} gives
[[1,-1,2,0,3], [0,2,1,1,0], [0,0,-17,-11,-12], [0,0,0,0,1-], [0,0,0,0,0]]
that's the correct result (but I knew the pivots *first*, so I could ...divide)

Salvo

Code:
 calcPivot(); zeropivot:=0; numpivot:=1; gj:=[0,0]; EXPORT gaussJordan() // Salvo Micciché 2015 // Credits: Dale (DrD), Claudio L. BEGIN RECT_P(); TEXTOUT_P("Gauss Jordan utilities", 20, 10, 5, RGB(0,0,255)); TEXTOUT_P("Pivots and upper matrix", 20, 40); TEXTOUT_P("LDLt factorization", 20, 60); TEXTOUT_P("LDU (LDV) factorization", 20, 80); TEXTOUT_P("Permutation Matrix from a list", 20, 100); TEXTOUT_P("Program created by Salvo Micciché", 20, 150, 3, RGB(0,255,0)); TEXTOUT_P("Credits: Dale (DrD), Claudio L.", 20, 180, 2, RGB(0,255,0)); WAIT; END; EXPORT Pivots(m)   // Gauss-Jordan elimination and pivots   // Salvo Micciché 2015 BEGIN   local r, c, piv, temp, len, j;   local permat, lumat, msquare;   r:=rowDim(m);   c:=colDim(m);   gj:=MAKEMAT(0,r,c);   M1:= m;   L1:={0,0}; // list of pivots   calcPivot(c,r);   IF (zeropivot==1) THEN   M1:=m;    if (numpivot==r) THEN numpivot:=1; END;  // if pivot is 0 swap rows (if last swap w/ 1)   // CAS.SWAPROW(M1, numpivot, (numpivot+1));   temp:= M1(numpivot); M1(numpivot):=M1(numpivot+1); M1(numpivot+1):=temp;   calcPivot(c,r);    END; // if   len:= length(L1);  // trick to get correct last pivot (product of pivot is determinant)   IF (r==c) THEN // only if Mat is square     lumat:= lu(m);     permat:= permMatrix(lumat(1)); // permutation matrix     IF (det(permat*m)<>0) THEN L1(len):= (det(permat*m) / (product(L1(X),X,1,len-1))); END;     gj(r,c):= L1(len);   ELSE      A:= L1(len);     B:= gj(r,c);     msquare:= SUB(m, {1,1}, {r,r});     lumat:= lu(msquare);     permat:= permMatrix(lumat(1)); // permutation matrix     IF (det(permat*msquare)<>0) THEN L1(len):= (det(permat*msquare) / (product(L1(X),X,1,len-1))); END;     gj(r,r):= L1(len);      IF (A<>0) THEN gj(r,c):= B*L1(len)/A; ELSE gj(r,c):= B; END;   END;   piv:= list2mat(L1, r);   M1:= gj;   // Return pivot, matrix reduced and RREF form   RETURN {piv, M1}; END; calcPivot(c,r) BEGIN   local j, k;   zeropivot:=0;   FOR j FROM 1 TO r DO     M1:=CAS.pivot(M1,1,1);     L1(j):= M1(1,1); // current pivot     IF (L1(j)==0) THEN zeropivot:=1; numpivot:=j; // if pivot is 0 swap rows     PRINT("Swap row " + j + " w/ " + IFTE(j<r, j+1, 1));      RETURN {zeropivot, numpivot}; END;     FOR k FROM 1 TO colDim(M1) DO     gj(j, c-k+1):= M1(1, colDim(M1)-k+1); // make [gj] (upper matrix) with correct rows     END; // inner for     IF (j<r) THEN  // take minor of the matrix     M1:= delrows(M1,1);     M1:= delcols(M1,1);   END; // if   M1:= M1 ./ L1(j);   IF L1(j)==0 THEN RETURN "Pivot is 0, division by 0"; END;    END; // for END; EXPORT LDLt(m)   // Factorizazion LDLt, matrix m must be simmetric: A=At BEGIN   Pivots(m);   local Lmat, LmatT, Dmat, r, j;   r:= rowDim(m);   Dmat:= diag(L1);   Lmat:= M1;   FOR j FROM 1 TO r DO   Lmat(j):= Lmat(j)/L1(j); // divide Lt (U) by pivots to get 1 in diagonal   END; // for   LmatT:= TRN(Lmat);   // Return L (lower) matrix, D matrix (pivots in diagonal), Lt (upper) transpose of L   RETURN {LmatT, Dmat, Lmat}; END; EXPORT LDU(m) // LDU factorization BEGIN   local Lmat, Umat, Dmat,r, j, temp;   local permu;   Pivots(m);  r:= rowDim(m);  Dmat:= diag(L1);  temp:= lu(m);  Lmat:= temp(2);  Umat:= temp(3);  permu:= permMatrix(temp(1)); FOR j FROM 1 TO r DO   Umat(j):= Umat(j)/L1(j); // divide U by pivots to get 1 in diagonal   END; // for   PRINT("Permutation matrix is " + temp(1));   // Return L (lower), D matrix (pivots in diagonal), U (upper) and permutation matrix   RETURN {Lmat, Dmat, Umat, permu}; END; EXPORT permMatrix(lst)   // given a permutation list (like {1,3,2}) and the matrix dimension, returns a permutation matrix BEGIN   local j, k, mat, dimen;   dimen:= length(lst);   mat:= MAKEMAT(0, dimen, dimen);   FOR j FROM 1 TO dimen DO   mat(j,lst(j)):= 1;   END; // for   RETURN mat; END; EXPORT echelon(m) BEGIN RETURN {ref(m), rref(m)}; END;

∫aL√0mic (IT9CLU) :: HP Prime 50g 41CX 71b 42s 39s 35s 12C 15C - DM42, DM41X - WP34s Prime Soft. Lib
06-07-2015, 10:42 AM
Post: #6
 DrD Senior Member Posts: 1,132 Joined: Feb 2014
RE: Pivots (Gaussian reduction)
Hi Salvo,

The program attached to your last post doesn't compile for me (permMatrix() isn't declared as a subroutine). After fixing that, I get errors with the Pivots sub program. Is this your latest edition?

-Dale-
06-07-2015, 11:28 AM (This post was last modified: 06-07-2015 12:55 PM by salvomic.)
Post: #7 salvomic Senior Member Posts: 1,394 Joined: Jan 2015
RE: Pivots (Gaussian reduction)
(06-07-2015 10:42 AM)DrD Wrote:  Hi Salvo,

The program attached to your last post doesn't compile for me (permMatrix() isn't declared as a subroutine). After fixing that, I get errors with the Pivots sub program. Is this your latest edition?

-Dale-

hi Dale,
try this one (here it compile).
The function are Pivots(), LDLt, LDU, permMatrix, echelon (gaussJordan is the about).
permMatrix() is now exported to be available externally...
I added the same in the program, as here it worked because I had another program with that routine...
permMatrix({1,3,2}) or perMatrix([1,3,2]) give the permutation matrix 3x3 for P32, and so on... i.ex. if you use lu() you get (last firmware) as first element something to [1,3,2], pass it to permMatrix and you'll get [[1,0,0], [0,0,1], [0,1,0]], the permutation matrix. The program use it internally also...

Use Pivots(matrix) to test, for now.
LDLt require a symmetric matrix... (Lt is the traspose of L, D matrix with pivots).
LDU (LDV) use lu() to get L*D*U, where D is the matrix that has diagonal with the pivots and 0 elsewhere...

Code:
 calcPivot(); permMatrix(); zeropivot:=0; numpivot:=1; gj:=[0,0]; EXPORT gaussJordan() // salvomic, Salvo Micciché 2015 // Credits: Dale (DrD), Claudio L. BEGIN RECT_P(); TEXTOUT_P("Gauss Jordan utilities", 20, 10, 5, RGB(0,0,255)); TEXTOUT_P("Pivots and upper matrix", 20, 40); TEXTOUT_P("LDLt factorization", 20, 60); TEXTOUT_P("LDU (LDV) factorization", 20, 80); TEXTOUT_P("Permutation Matrix from a list", 20, 100); TEXTOUT_P("Program created by Salvo Micciché", 20, 150, 3, RGB(0,255,0)); TEXTOUT_P("Credits: Dale (DrD), Claudio L.", 20, 180, 2, RGB(0,255,0)); WAIT; END; EXPORT Pivots(m)   // Gauss-Jordan elimination and pivots   // Salvo Micciché 2015 BEGIN   local r, c, piv, temp, len, j;   local permat, lumat, msquare;   r:=rowDim(m);   c:=colDim(m);   gj:=MAKEMAT(0,r,c);   M1:= m;   L1:={0,0}; // list of pivots   calcPivot(c,r);   IF (zeropivot==1 AND numpivot<>0) THEN     M1:=m;      IF (numpivot==r) THEN numpivot:=1; END;  // if pivot is 0 swap rows (if last swap w/ 1)    FOR j FROM numpivot TO (r-1) DO     IF (numpivot==r) THEN BREAK; END;     IF (M1(j+1,1)==0) THEN      CASE       IF (r<(r-1)) THEN CONTINUE; END;       IF (r>=(r-1)) THEN BREAK; END;       DEFAULT     END; // case      END; // if     // CAS.SWAPROW(M1, numpivot, (numpivot+1));     temp:= M1(numpivot); M1(numpivot):=M1(j+1); M1(j+1):=temp;      calcPivot(c,r);    END; // for    END; // if   len:= length(L1);  // trick to get correct last pivot (product of pivot is determinant)   IF (r==c) THEN // only if Mat is square     A:= L1(len); // last pivot     lumat:= lu(m);     permat:= permMatrix(lumat(1)); // permutation matrix     P:= product(L1(X),X,1,len-1); // product of pivots (= det)     IF (det(permat*m)<>0 AND P<>0) THEN L1(len):= (det(permat*m) / P); END;     gj(r,c):= L1(len);   ELSE      A:= L1(len); // last pivot     // B:= gj(r,c); // last element     P:= product( L1(X),X,1,len-1 );     msquare:= SUB(m, {1,1}, {r,r}); // subset square for a "pseudo-determinant"     lumat:= lu(msquare);     permat:= permMatrix(lumat(1)); // permutation matrix     IF (det(permat*msquare)<>0 AND P<>0) THEN L1(len):= (det(permat*msquare) / P); END;     gj(r,r):= L1(len); // corrected pivot copied into item r,r     // IF (A<>0) THEN gj(r,c):= B*L1(len)/A; ELSE gj(r,c):= B; END;     FOR j FROM (r+1) TO (c) DO       B:= gj(r,j); // current elements after (r,r)       IF (A<>0) THEN gj(r,j):= B*L1(len)/A; ELSE gj(r,j):= B; END;     END; // for   END;   piv:= list2mat(L1, r);   M1:= gj;   // Return pivot and reduce matrix with pivots in diagonal   RETURN {piv, M1}; END; calcPivot(c,r) BEGIN   local j, k;   zeropivot:=0;   FOR j FROM 1 TO r DO     M1:= pivot(M1,1,1);     L1(j):= M1(1,1); // current pivot     IF (L1(j)==0) THEN zeropivot:=1; numpivot:=j; // if pivot is 0 swap rows      PRINT("Swap row " + j + " w/ " + IFTE(j<r, j+1, 1));       RETURN {zeropivot, numpivot};      END;     FOR k FROM 1 TO colDim(M1) DO      gj(j, c-k+1):= M1(1, colDim(M1)-k+1); // make [gj] (upper matrix) with correct rows     END; // inner for     IF (j<r) THEN  // take minor of the matrix      M1:= delrows(M1,1);      M1:= delcols(M1,1);     END; // if   M1:= M1 ./ L1(j); // divide matrix for current pivot   IF L1(j)==0 THEN RETURN "Pivot is 0, division by 0"; END;    END; // for END; EXPORT LDLt(m)   // Factorizazion LDLt, matrix m must be simmetric: A=At BEGIN   Pivots(m);   local Lmat, LmatT, Dmat, r, j;   r:= rowDim(m);   Dmat:= diag(L1);   Lmat:= M1;   FOR j FROM 1 TO r DO   Lmat(j):= Lmat(j)/L1(j); // divide Lt (U) by pivots to get 1 in diagonal   END; // for   LmatT:= TRN(Lmat);   // Return L (lower) matrix, D matrix (pivots in diagonal), Lt (upper) transpose of L   RETURN {LmatT, Dmat, Lmat}; END; EXPORT LDU(m) // LDU factorization BEGIN   local Lmat, Umat, Dmat,r, j, temp;   local permu;   Pivots(m);  r:= rowDim(m);  Dmat:= diag(L1);  temp:= lu(m);  Lmat:= temp(2);  Umat:= temp(3);  permu:= permMatrix(temp(1)); FOR j FROM 1 TO r DO   Umat(j):= Umat(j)/L1(j); // divide U by pivots to get 1 in diagonal   END; // for   PRINT("Permutation matrix is " + temp(1));   // Return L (lower), D matrix (pivots in diagonal), U (upper) and permutation matrix   RETURN {Lmat, Dmat, Umat, permu}; END; EXPORT permMatrix(lst)   // given a permutation list (like {1,3,2} or [1 3 2]), returns a permutation matrix BEGIN   local j, k, mat, dimen;   dimen:= length(lst);   mat:= MAKEMAT(0, dimen, dimen);   FOR j FROM 1 TO dimen DO   mat(j,lst(j)):= 1;   END; // for   RETURN mat; END; EXPORT echelon(m) BEGIN // simply return ref() and RREF() for reference RETURN {ref(m), rref(m)}; END;

∫aL√0mic (IT9CLU) :: HP Prime 50g 41CX 71b 42s 39s 35s 12C 15C - DM42, DM41X - WP34s Prime Soft. Lib
06-07-2015, 12:13 PM
Post: #8
 DrD Senior Member Posts: 1,132 Joined: Feb 2014
RE: Pivots (Gaussian reduction)

FOR j FROM 1 TO r DO
M1:=pivot(M1,1,1); // Changed CAS.pivot(M1,1,1) to pivot(M1,1,1), <<clears syntax error>>
06-07-2015, 12:41 PM (This post was last modified: 06-08-2015 10:07 AM by salvomic.)
Post: #9 salvomic Senior Member Posts: 1,394 Joined: Jan 2015
RE: Pivots (Gaussian reduction)
(06-07-2015 12:13 PM)DrD Wrote:  In your routine: calcPivot(c,r)

FOR j FROM 1 TO r DO
M1:=pivot(M1,1,1); // Changed CAS.pivot(M1,1,1) to pivot(M1,1,1), <<clears syntax error>>

ok, changed it!

now I get right results with the matrix in some post above...

Still wrong result with M9 = [[0,2,1,1,0], [1,-1,2,0,3], [1,1,3,1,2], [3,1,8,2,6], [4,2,-6,-8,0]].
For this I put "echelon()" as function to compare with ref() and RREF(). For now Pivots(M9) gives only the first two rows and the last three with zeros...
Surely, after the row swapping pivot() doesn't offer the input that the program require...
EDIT: with the last version I get for M9 the first 3 rows: why?

More modular version, with a bit more clean code and some others //comments to help ...who want help me to debug Corrected a bug if the first row was 0...

Code:
 calcPivot(); permMatrix(); swapPivot(); zeropivot:=0; numpivot:=1; gj:=[0,0]; msw:=[0,0]; EXPORT gaussJordan() // salvomic, Salvo Micciché 2015 // Credits: Dale (DrD), Claudio L. BEGIN RECT_P(); TEXTOUT_P("Gauss Jordan utilities", 20, 10, 5, RGB(0,0,255)); TEXTOUT_P("Pivots and upper matrix", 20, 40); TEXTOUT_P("LDLt factorization", 20, 60); TEXTOUT_P("LDU (LDV) factorization", 20, 80); TEXTOUT_P("Permutation Matrix from a list", 20, 100); TEXTOUT_P("Echelon: reduced of echelon reduced form", 20, 120); TEXTOUT_P("Program created by Salvo Micciché", 20, 160, 3, RGB(0,255,0)); TEXTOUT_P("Credits: Dale (DrD), Claudio L.", 20, 180, 3, RGB(0,255,0)); WAIT; END; EXPORT Pivots(m)   // Gauss-Jordan elimination and pivots BEGIN   local r, c, piv, temp, len, j;   local permat, lumat, msquare, conf, row;   r:=rowDim(m);   c:=colDim(m);   gj:=MAKEMAT(0,r,c);   msw:=MAKEMAT(0,r,c);   M1:= m; msw:=m; // save matrix in M1 and msw for reference   L1:={0,0}; // dim list of pivots   IF M1(1,1)==0 THEN // find a non-zero pivot for first line or swap     conf:= 100; row:=2;     FOR j FROM 2 TO r DO       IF M1(j,1)<>0 THEN         IF ABS(M1(j,1))<ABS(conf) THEN conf:=M1(j,1); row:=j; END;         // hint for row w/ less 1st item        END; // if     END; // for     temp:= M1(1); M1(1):=M1(row); M1(row):=temp;     PRINT("Swap row " + 1 + " w/ " + row); PRINT(" ");     msw:= M1; // matrix with swapped rows   END; // if find...   calcPivot();   IF (zeropivot==1 AND numpivot<>0) THEN // if calPivot requires swap   M1:= msw;   swapPivot();   calcPivot(); // called again to calc pivots after a swapping of rows   END; // if   len:= length(L1);  // trick to get correct last pivot (product of pivot is determinant)   IF (r==c) THEN // only if Mat is square     A:= L1(len); // last pivot     lumat:= lu(m);     permat:= permMatrix(lumat(1)); // permutation matrix     P:= product(L1(X),X,1,len-1); // product of pivots but not the last (product of d = det)      IF (det(permat*m)<>0 AND P<>0) THEN L1(len):= (det(permat*m) / P);       ELSE L1(len):=A; END;     gj(r,c):= L1(len);   ELSE  // matrix (r,c) not square (complete by vector b)     A:= L1(len); // last pivot     msquare:= SUB(m, {1,1}, {r,r}); // subset square for a "pseudo-determinant"     lumat:= lu(msquare);     permat:= permMatrix(lumat(1)); // permutation matrix     P:= product( L1(X),X,1,len-1 );      IF (det(permat*msquare)<>0 AND P<>0) THEN L1(len):= (det(permat*msquare) / P);       ELSE L1(len):=A; END;     gj(r,r):= L1(len); // corrected pivot copied into item r,r     FOR j FROM (r+1) TO (c) DO          B:= gj(r,j); // current elements after (r,r)             IF (A<>0) THEN gj(r,j):= B*L1(len)/A; ELSE gj(r,j):= B; END;     END; // for   END; // if trick   piv:= list2mat(L1, r);   M1:= gj;  // Return pivot and reduce matrix with pivots in diagonal   RETURN {piv, M1}; END; calcPivot() BEGIN   local j, k, r, c;   r:= rowDim(M1);   c:= colDim(M1);   zeropivot:=0;   FOR j FROM 1 TO r DO     M1:= pivot(M1,1,1);  // partial pivoting first element     L1(j):= M1(1,1); // current pivot     IF (L1(j)==0) THEN   zeropivot:=1; numpivot:=j; // if pivot is 0 swap rows      RETURN {zeropivot, numpivot};      END; // if     FOR k FROM 1 TO colDim(M1) DO      gj(j, c-k+1):= M1(1, colDim(M1)-k+1); // make [gj] (upper matrix) with correct rows     END; // inner for     IF (j<r) THEN  // take minor of the matrix but last row      M1:= delrows(M1,1);      M1:= delcols(M1,1);     END; // if   M1:= M1 ./ L1(j); // divide matrix for current pivot   IF L1(j)==0 THEN RETURN "Pivot is 0, division by 0"; END;    END; // for END; swapPivot() BEGIN   local j, temp, r;   r:= rowDim(M1);     FOR j FROM numpivot+1 TO r DO       IF M1(j,1)<>0 THEN    temp:= M1(numpivot); M1(numpivot):=M1(j); M1(j):=temp;   PRINT("-- Swap row " + numpivot + " w/ " + j);   END; // if    END; // for  END; EXPORT LDLt(m)   // Factorizazion LDLt, matrix m must be simmetric: A=At BEGIN   Pivots(m);   local Lmat, LmatT, Dmat, r, j;   r:= rowDim(m);   Dmat:= diag(L1);   Lmat:= M1;   FOR j FROM 1 TO r DO   Lmat(j):= Lmat(j)/L1(j); // divide Lt (U) by pivots to get 1 in diagonal   END; // for   LmatT:= TRN(Lmat);   // Return L (lower) matrix, D matrix (pivots in diagonal), Lt (upper) transpose of L   RETURN {LmatT, Dmat, Lmat}; END; EXPORT LDU(m) // LDU factorization BEGIN   local Lmat, Umat, Dmat,r, j, temp;   local permu;   Pivots(m);  r:= rowDim(m);  Dmat:= diag(L1);  temp:= lu(m);  Lmat:= temp(2);  Umat:= temp(3);  permu:= permMatrix(temp(1)); FOR j FROM 1 TO r DO   Umat(j):= Umat(j)/L1(j); // divide U by pivots to get 1 in diagonal   END; // for   PRINT("Permutation matrix is " + temp(1));   // Return L (lower), D matrix (pivots in diagonal), U (upper) and permutation matrix   RETURN {Lmat, Dmat, Umat, permu}; END; EXPORT permMatrix(lst)   // given a permutation list (like {1,3,2} or [1 3 2]), returns a permutation matrix BEGIN   local j, k, mat, dimen;   dimen:= length(lst);   mat:= MAKEMAT(0, dimen, dimen);   FOR j FROM 1 TO dimen DO   mat(j,lst(j)):= 1;   END; // for   RETURN mat; END; EXPORT echelon(m) BEGIN // simply return ref() and RREF() for reference RETURN {ref(m), rref(m)}; END;

∫aL√0mic (IT9CLU) :: HP Prime 50g 41CX 71b 42s 39s 35s 12C 15C - DM42, DM41X - WP34s Prime Soft. Lib
06-08-2015, 04:02 PM (This post was last modified: 06-08-2015 07:53 PM by salvomic.)
Post: #10 salvomic Senior Member Posts: 1,394 Joined: Jan 2015
RE: Pivots (Gaussian reduction)
With this new version I get good results for almost 20 matrices (diverse types observed: 2x2, 3x3, 4x4, 5x5, singular, echelon...), among them the M2-M9 from above...

However I used two routines I called "tricks" (see line 55 and 106): I suspect they run as workaround and not as real formulas (if you want, please, control, advise...), but they work almost in the case examined by me.

To have some global variables, for now I used M1, L1, A, B, D, E variables. Expect M1, L1, the other, perhaps, could be changed in local ones...

I'm working for a more clean and better code, without "tricks" maybe EDIT: ok with about 20 types, but if fails (in part) with this "band, 3 diagonal" matrix:
[[2,-1,0,0,0], [-1,2,-1,0,0], [0,-1,2,-1,0], [0,0,-1,2,-1], [0,0,0,-1,2]]
with Pivots() I get {2, 3/2, ⅔, 5/12, 36/5} but the pivots should be {2/1, 3/2, 4/3, 5/4, 6/5}; LDLt gives right L and Lt (transpose) but not D (for the same reason): any help?

Code:
 calcPivot(); permMatrix(); swapPivot(); zeropivot:=0; numpivot:=1; gj:=[0,0]; msw:=[0,0]; EXPORT gaussJordan() // salvomic, Salvo Micciché 2015 // Credits: Dale (DrD), Claudio L. BEGIN RECT_P(); TEXTOUT_P("Gauss Jordan utilities", 20, 10, 5, RGB(0,0,255)); TEXTOUT_P("Pivots and upper matrix", 20, 40); TEXTOUT_P("LDLt factorization", 20, 60); TEXTOUT_P("LDU (LDV) factorization", 20, 80); TEXTOUT_P("Permutation Matrix from a list", 20, 100); TEXTOUT_P("Echelon: reduced of echelon reduced form", 20, 120); TEXTOUT_P("Program created by Salvo Micciché", 20, 160, 3, RGB(0,255,0)); TEXTOUT_P("Credits: Dale (DrD), Claudio L.", 20, 180, 3, RGB(0,255,0)); WAIT; END; EXPORT Pivots(m)   // Gauss-Jordan elimination and pivots BEGIN   local r, c, piv, temp, len, j;   local permat, lumat, msquare, conf, row;   r:=rowDim(m);   c:=colDim(m);   gj:=MAKEMAT(0,r,c);   msw:=MAKEMAT(0,r,c);   M1:= m; msw:=m; // save matrix in M1 and msw for reference   L1:=MAKELIST(0,X,1,r); // dim list of pivots   IF M1(1,1)==0 THEN // find a non-zero pivot for first line or swap     conf:= 100; row:=2;     FOR j FROM 2 TO r DO       IF M1(j,1)<>0 THEN         IF ABS(M1(j,1))<ABS(conf) THEN conf:=M1(j,1); row:=j; END;         // hint for row w/ less 1st item        END; // if     END; // for     temp:= M1(1); M1(1):=M1(row); M1(row):=temp;     PRINT("Swap row " + 1 + " w/ " + row); PRINT(" ");     msw:= M1; // matrix with swapped rows   END; // if find...   calcPivot();   IF (zeropivot==1 AND numpivot<>0) THEN // if calPivot requires swap   M1:= msw;   swapPivot();   calcPivot(); // called again to calc pivots after a swapping of rows   END; // if   len:= length(L1);  // trick to get correct last pivot (product of pivot is determinant)   IF (r==c) THEN // only if Mat is square     A:= L1(len); // last pivot     lumat:= lu(m);     permat:= permMatrix(lumat(1)); // permutation matrix     P:= product(L1(X),X,1,len-1); // product of pivots but not the last (product of d = det)      IF (det(permat*m)<>0 AND P<>0) THEN L1(len):= (det(permat*m) / P);       ELSE L1(len):=A; END;     gj(r,c):= L1(len);     L1(len-1):= IFTE(E==0,D,E); // pseudo-pivot for echelon mat in some cases with 0 for pivot...   ELSE  // matrix (r,c) not square (complete by vector b)     A:= L1(len); // last pivot     msquare:= SUB(m, {1,1}, {r,r}); // subset square for a "pseudo-determinant"     lumat:= lu(msquare);     permat:= permMatrix(lumat(1)); // permutation matrix     P:= product( L1(X),X,1,len-1 );      IF (det(permat*msquare)<>0 AND P<>0) THEN L1(len):= (det(permat*msquare) / P);       ELSE L1(len):=A; END;     gj(r,r):= L1(len); // corrected pivot copied into item r,r     L1(len-1):= IFTE(E==0,D,E); // pseudo-pivot for echelon mat (see calcPivot)     FOR j FROM (r+1) TO (c) DO          B:= gj(r,j); // current elements after (r,r)             IF (A<>0) THEN gj(r,j):= B*L1(len)/A; ELSE gj(r,j):= B; END;     END; // for   END; // if trick   piv:= list2mat(L1, r);   M1:= gj;  // Return pivot and reduce matrix with pivots in diagonal   RETURN {piv, M1}; END; calcPivot() BEGIN   local j, k, r, c;   r:= rowDim(M1);   c:= colDim(M1);   zeropivot:=0;   FOR j FROM 1 TO r DO     M1:= pivot(M1,1,1);  // partial pivoting first element     L1(j):= M1(1,1); // current pivot     IF (L1(j)==0) THEN     zeropivot:=1; numpivot:=j; // if pivot is 0 swap rows       RETURN {zeropivot, numpivot};      END; // if     FOR k FROM 1 TO colDim(M1) DO      gj(j, c-k+1):= M1(1, colDim(M1)-k+1); // make [gj] (upper matrix) with correct rows     END; // inner for (k)     IF (j<r) THEN  // take minor of the matrix but last row      M1:= delrows(M1,1);      M1:= delcols(M1,1);     END; // if     IF (j==(r-2)  AND M1(1,1)==0) THEN // another "trick" to get echelon item for last but one row     D:= IFTE(M1(1,colDim(M1))<>0, M1(1,colDim(M1)), M1(2,colDim(M1))); // correct pseudo-pivot (echelon)      E:= IFTE(M1(1,colDim(M1)-1)<>0, M1(1,colDim(M1)-1), M1(2,colDim(M1)-1));      gj(r-1,c):= D;      gj(r-1,c-1):= E;     ELSE       D:= L1(length(L1)-1); E:= L1(length(L1)-1); // old pivots, if no problem     END; // trick2     M1:= M1 ./ L1(j); // divide matrix for current pivot     IF L1(j)==0 THEN RETURN "Pivot is 0, division by 0"; END;    END; // for j END; swapPivot() BEGIN   local j, temp, r;   r:= rowDim(M1);     FOR j FROM numpivot+1 TO r DO       IF M1(j,1)<>0 THEN    temp:= M1(numpivot); M1(numpivot):=M1(j); M1(j):=temp;   PRINT("-- Swap row " + numpivot + " w/ " + j);   END; // if    END; // for  END; EXPORT LDLt(m)   // Factorizazion LDLt, matrix m must be simmetric: A=At BEGIN   Pivots(m);   local Lmat, LmatT, Dmat, r, j;   r:= rowDim(m);   Dmat:= diag(L1);   Lmat:= M1;   FOR j FROM 1 TO r DO   Lmat(j):= Lmat(j)/L1(j); // divide Lt (U) by pivots to get 1 in diagonal   END; // for   LmatT:= TRN(Lmat);   // Return L (lower) matrix, D matrix (pivots in diagonal), Lt (upper) transpose of L   RETURN {LmatT, Dmat, Lmat}; END; EXPORT LDU(m) // LDU factorization BEGIN   local Lmat, Umat, Dmat,r, j, temp;   local permu;   Pivots(m);  r:= rowDim(m);  Dmat:= diag(L1);  temp:= lu(m);  Lmat:= temp(2);  Umat:= temp(3);  permu:= permMatrix(temp(1)); FOR j FROM 1 TO r DO   Umat(j):= Umat(j)/L1(j); // divide U by pivots to get 1 in diagonal   END; // for   PRINT("Permutation matrix is " + temp(1));   // Return L (lower), D matrix (pivots in diagonal), U (upper) and permutation matrix   RETURN {Lmat, Dmat, Umat, permu}; END; EXPORT permMatrix(lst)   // given a permutation list (like {1,3,2} or [1 3 2]), returns a permutation matrix BEGIN   local j, k, mat, dimen;   dimen:= length(lst);   mat:= MAKEMAT(0, dimen, dimen);   FOR j FROM 1 TO dimen DO   mat(j,lst(j)):= 1;   END; // for   RETURN mat; END; EXPORT echelon(m) BEGIN // simply return ref() and RREF() for reference RETURN {ref(m), rref(m)}; END;

∫aL√0mic (IT9CLU) :: HP Prime 50g 41CX 71b 42s 39s 35s 12C 15C - DM42, DM41X - WP34s Prime Soft. Lib
06-08-2015, 10:15 PM (This post was last modified: 06-09-2015 10:02 AM by salvomic.)
Post: #11 salvomic Senior Member Posts: 1,394 Joined: Jan 2015
RE: Pivots (Gaussian reduction)
a different approach is not to use pivot() command to calc Pivots(), but use lu()...
pivot() command calculates partial pivoting element by element, it isn't a complete gaussian elimination, so I had to use the ugly "tricks" to get the job...

LU returns P (permutation, easily transformable into a matrix with permMatrix()), L (lower triangular that contains "multipliers" below the diagonal) and U (upper triangular that should have pivots, reduced, in diagonal).
I need (of Pivots() ) to get a diagonal with pivots *not* reduced (and the product of which is equal to determinant, if the matrix is not singular), as then I get LDU (with LDU() ) and LDLt (with LDLt() ), that require a D matrix with pivots in its diagonal...

I need an hint to extract multipliers from L returned by lu() and make with them the gaussian reduction (with or without permutation, given by the first list of lu())...
I thought a cycle...

i.e. if matrix is M2:= [[2,1,1], [4,-6,0], [-2,7,2]], lu(M2) returns [ [1 2 3] [[1,0,0], [-1,1,0], [2,-1,1]] [[2,1,1], [0,8,3], [0,0,1]] ]
so the permutation matrix is P = [[1,0,0], [0,01], [0,1,0]]
L is the second matrix in the list, and in it we have multipliers -1 (L(2,1)), 2 (L(3,1)), -1 (L(3,2)), so we need to subtract -1 times 1st row from 2nd, 2 times 1st row from 3rd, -1 times 2nd row from 3rd
U should have pivots already in the diagonal, but in this case it has 2,8,1, but manually we get 2, -8, 1 (the actual form of my program returns in fact Pivots(M2) -> {2,-8,1} and the product of pivots must be equal to determinant that for M2 is -16 (2*1*(-8)), while lu() returns {2,8,1} because it consider the matrix with the rows swapped as (1,3,2) and in this case the determinant *is* 16...

The cycle should extract multipliers, perform permutation (multiplying P*M2) and subtractions (following the multipliers) and return Gauss reduction and the required pivots.

Any help?
Salvo

***
EDIT 2
...but this routine (and the new cycle) only reproduces LU and returns ...simply U. I'm outside of the way...

Code:
 EXPORT Pivots(m)   // Gauss-Jordan elimination and pivots BEGIN   local r, c, j, k, piv, temp, plist;   local permat, Lmat, Umat, Pmat;   r:=rowDim(m);   c:=colDim(m);   L1:=MAKELIST(0,X,1,r); // dim list of pivots   gj:=MAKEMAT(0,r,c); // building matrix   gj:= m;   M1:= SUB(m, {1,1},{r,r}); // here M1 is the square subset of matrix   temp:= lu(M1);   plist:= temp(1);   Lmat:= temp(2);   Umat:= temp(3);   Pmat:= permMatrix(plist);   M1:= Pmat*M1;   gj:= Pmat*gj;   FOR j FROM 1 TO r DO  // calc multiplier and subtract rows     FOR k FROM 1 TO j DO       IF (j<>k) THEN         IF (Lmat(j,k)<>0) THEN gj(j):= gj(j) - Lmat(j,k)*gj(k); END;       END; // if     END; // inner for   END; // outer for   FOR j FROM 1 TO r DO // extract pivots     FOR k FROM 1 TO c DO         IF (gj(j,k)<>0) THEN L1(j):= gj(j,k); BREAK; END;     END; // inner for   END; //outer for   piv:= list2mat(L1, r);   M1:= gj;  // Return pivot and reduce matrix with pivots in diagonal   RETURN {piv, M1}; END;

∫aL√0mic (IT9CLU) :: HP Prime 50g 41CX 71b 42s 39s 35s 12C 15C - DM42, DM41X - WP34s Prime Soft. Lib
06-09-2015, 10:55 AM
Post: #12
 DrD Senior Member Posts: 1,132 Joined: Feb 2014
RE: Pivots (Gaussian reduction)
Recognizing that Gauss Jordan elimination can be accomplished in many different ways, Gauss Jordan, LU decomposition would be dependent on the algorithm used to perform the reduction. With this in mind, check your results with that of Wolfram Alpha once again: LU(M2) from Wolfram Alpha. Specifically, notice that the diagonal of their upper triangular matrix U is [2,-8, 1], yet your decomposition is different. To check the results of the decomposition, L * U / P, should return the original M2 matrix, but in your case it doesn't. The permutation matrix shown, is not correct.

Another decomposition for LU(M2) is:

{[[1,0,0],[0.5,1,0],[−0.5,1,1]],[[4,−6,0],[0,4,1],[0,0,1]],[[0,1,0],[1,0,0],[0,0,1]]}, where the matrices are {L,U,P}. Checking: L*U/P ==M2.

There may be another update released before long, and hopefully resolve some of the issues.
06-09-2015, 12:50 PM (This post was last modified: 06-11-2015 10:45 PM by salvomic.)
Post: #13 salvomic Senior Member Posts: 1,394 Joined: Jan 2015
RE: Pivots (Gaussian reduction)
(06-09-2015 10:55 AM)DrD Wrote:  Recognizing that Gauss Jordan elimination can be accomplished in many different ways, ...

Another decomposition for LU(M2) is:

{[[1,0,0],[0.5,1,0],[−0.5,1,1]],[[4,−6,0],[0,4,1],[0,0,1]],[[0,1,0],[1,0,0],[0,0,1]]}, where the matrices are {L,U,P}. Checking: L*U/P ==M2.

There may be another update released before long, and hopefully resolve some of the issues.

yes, I know, in fact I could stop with the LU (and REF) that already Prime has, but I search also for others way, to control some results by books and for study purpose...
My first attempt is no completely sure, as it uses two "tricks", the second attempt give no more nor less U as the LU() command, reproducing it in another way I made programs more difficult in less time than this one, but this is fascinating me, in some way...

An yes, I control always with Wolfram and with the reverse formula...

At the start I thought that the command pivot() gave the complete reduction, but I can see that applying it more time could get exactly LU or REF or ...RREF.

Anyway, my purpose is to get LDU (LDV) and LDLt, and for now my first attempt is the nearest to them. Maybe.

Salvo

***
The last code with utilities: Pivots, LU_ext (LU also if the matrix is not square), LDLt, LDU, permMatrix, echelon (and a few bugs corrected) is here. It works much better than first and it is conform to some fonts, like the books of Gilbert Strang and a few others, obviously not with *every* matrices...

I promise don't bore you more with this program... EDIT June 12 2015: every new version will be put in this thread of the Software Library of the Forum

Code:
 calcPivot(); permMatrix(); swapPivot(); zeropivot:=0; numpivot:=1; swap:= 0; gj:=[0,0]; msw:=[0,0]; EXPORT gaussJordan() // salvomic, Salvo Micciché 2015 // Credits: Dale (DrD), Claudio L. BEGIN RECT_P(); TEXTOUT_P("Gauss Jordan utilities", 20, 10, 5, RGB(0,0,255)); TEXTOUT_P("Pivots and upper matrix", 20, 40); TEXTOUT_P("LU extended: for not square mat", 20, 60); TEXTOUT_P("LDLt factorization", 20, 80); TEXTOUT_P("LDU (LDV) factorization", 20, 100); TEXTOUT_P("Permutation Matrix from a list", 20, 120); TEXTOUT_P("Echelon: reduced of echelon reduced form", 20, 140); TEXTOUT_P("Program created by Salvo Micciché", 20, 180, 3, RGB(0,255,0)); TEXTOUT_P("Credits: Dale (DrD), Claudio L.", 20, 200, 3, RGB(0,255,0)); WAIT; END; EXPORT Pivots(m)   // Gauss-Jordan elimination and pivots BEGIN   local r, c, temp, len, j;   local permat, lumat, msquare, conf, row;   r:=rowDim(m);   c:=colDim(m);   L1:=MAKELIST(0,X,1,r); // dim list of pivots   gj:=MAKEMAT(0,r,c);   msw:=MAKEMAT(0,r,c);   M1:= m; msw:=m; // save matrix in M1 and msw for reference   swap:= 0;   IF M1(1,1)==0 THEN // find a non-zero pivot for first line or swap     conf:= 1E10; row:=2;     FOR j FROM 2 TO r DO       IF M1(j,1)<>0 THEN         IF ABS(M1(j,1))<ABS(conf) THEN conf:=M1(j,1); row:=j; END;         // hint for row w/ lesser 1st item        END; // if     END; // for     temp:= M1(1); M1(1):=M1(row); M1(row):=temp;     // PRINT; PRINT("Swap row " + 1 + " w/ " + row);     msw:= M1; // matrix with swapped rows   END; // if find...   calcPivot();   IF (zeropivot==1 AND numpivot<>0) THEN // if calPivot requires swap   M1:= msw;   swapPivot();   calcPivot(); // called again to calc pivots after a swapping of rows     FOR j FROM 2 TO length(L1) DO // divide by pivots if swapping    IF L1(j)<>0 THEN gj(j):= gj(j) * L1(j-1)^(-1); END;      END; // for   END; // if   M1:= gj;   IF (r==c) THEN L1:= diag(M1); ELSE temp:=SUB(M1,{1,1}, {r,r}); L1:= diag(temp); END;    len:= length(L1);  // trick 1 to get correct last pivot (product of pivot is determinant)    P:= product(L1(X),X,1,len-1); // product of pivots but not the last (product of d = det)    msquare:= SUB(m, {1,1}, {r,r}); // square subset  ("pseudo-determinant" if mat isn't square)    lumat:= lu(msquare);    permat:= permMatrix(lumat(1)); // permutation matrix    IF (det(msquare)<>0 AND P<>0) THEN      IF swap==0 THEN L1(len):= det(msquare)/P; ELSE L1(len):= det(permat*msquare)/P; END;    END; // if    A:=EVAL(gj(r,r));    gj(r,r):= L1(len); // corrected pivot copied into item r,r    IF (r<>c) THEN // mat isn't square     FOR j FROM r+1 TO c DO       IF (gj(r,j)<>0) THEN gj(r,j):= gj(r,j)*gj(r,r)/A; END;     END; // for    END; // if - end trick    M1:= gj;  // Return pivots in diagonal and U (upper with pivots)    IF (r==c) THEN L1:= diag(M1); ELSE temp:=SUB(M1,{1,1}, {r,r}); L1:= diag(temp); END;    RETURN {diag(L1), M1}; END; calcPivot() BEGIN   local j, k, r, c, temp;   r:= rowDim(M1);   c:= colDim(M1);   zeropivot:=0;   FOR j FROM 1 TO r DO     M1:= pivot(M1,1,1);  // partial pivoting first element     L1(j):= M1(1,1); // current pivot     IF (L1(j)==0) THEN     zeropivot:=1; numpivot:=j; // if pivot is 0 swap rows      RETURN {zeropivot, numpivot};      END; // if     FOR k FROM 1 TO colDim(M1) DO     gj(j, c-k+1):= M1(1, colDim(M1)-k+1); // make [gj] (upper matrix) with correct rows     END; // inner for (k)     IF (j<r) THEN  // take minor of the matrix but last row      M1:= delrows(M1,1);      M1:= delcols(M1,1);     END; // if     IF (j==(r-2) AND M1(1,1)==0) THEN // trick 2 to get echelon item for last but one row       temp:= IFTE(M1(1,colDim(M1))<>0, M1(1,colDim(M1)), M1(2,colDim(M1))); // correct pseudo-pivot       gj(r-1,c):= IFTE(temp<>0, temp/gcd(temp), temp);       temp:= IFTE(M1(1,colDim(M1)-1)<>0, M1(1,colDim(M1)-1), M1(2,colDim(M1)-1));       gj(r-1,c-1):= IFTE(temp<>0, temp/gcd(temp), temp);     END; // trick2   END; // for j   FOR j FROM 2 TO length(L1) DO    // to divide rows for current pivot    IF L1(j-1)<>0 THEN gj(j):= gj(j) * L1(j-1)^(-1); END;    END; // for END; swapPivot() BEGIN   local j, temp, r;   r:= rowDim(M1);   FOR j FROM numpivot+1 TO r DO       IF M1(j,1)<>0 THEN        temp:= M1(numpivot); M1(numpivot):=M1(j); M1(j):=temp;       swap:= 1;       // PRINT("-- Swap row " + numpivot + " w/ " + j);   END; // if   END; // for  END; EXPORT LU_ext(m)   // Gauss-Jordan elimination and pivots: LU also for rectangular matrices BEGIN   local r, c, j, k, temp, plist;   local permat, Lmat, Umat, Pmat;   r:=rowDim(m);   c:=colDim(m);   L1:=MAKELIST(0,X,1,r); // dim list of pivots   gj:=MAKEMAT(0,r,c); // building matrix   gj:= m;   M1:= SUB(m, {1,1},{r,r}); // here M1 is the square subset of matrix   temp:= lu(M1);   plist:= temp(1);   Lmat:= temp(2);   Umat:= temp(3);   Pmat:= permMatrix(plist);   M1:= Pmat*M1;   gj:= Pmat*gj;   FOR j FROM 1 TO r DO  // calc multiplier from Lmat and subtract rows   FOR k FROM 1 TO j DO     IF (j<>k) THEN     IF (Lmat(j,k)<>0) THEN gj(j):= gj(j) - Lmat(j,k)*gj(k); END;     END; // if   END; // inner for   END; // outer for   FOR j FROM 1 TO r DO // extract pivots   FOR k FROM 1 TO c DO     IF (gj(j,k)<>0) THEN L1(j):= gj(j,k); BREAK; END;   END; // inner for   END; //outer for   // adapt last two rows if echelon   IF (gj(r,c-1)<>0 AND gj(r-1, c-1)<>0 AND gj(r-1, c-2)==0) THEN     temp:= gj(r,c-1); gj(r):= gj(r) - temp / gj(r-1,c-1)*gj(r-1);      Lmat(r,r-1):= temp / gj(r-1,c-1);     L1(length(L1)):= 0;   END;    PRINT("Permutation matrix is " + plist);   M1:= gj;  // Return L, U and matrix with pivots in diagonal   RETURN {diag(L1), Lmat, M1}; END; EXPORT LDLt(m)   // Factorizazion LDLt, matrix m must be simmetric: A=At BEGIN   local Lmat, LmatT, Dmat, r,c, j;   r:= rowDim(m);   c:= colDim(m);   IF (r<>c) THEN RETURN ("Matrix must be square and symmetric"); END;   Pivots(m);   Dmat:= diag(L1);   Lmat:= M1;   FOR j FROM 1 TO r DO   IF L1(j)<>0 THEN Lmat(j):= Lmat(j)/L1(j); END;    // divide Lt (U) by pivots to get 1 in diagonal   END; // for   LmatT:= TRN(Lmat);   // Return L (lower) matrix, D matrix (pivots in diagonal), Lt (upper) transpose of L   RETURN {LmatT, Dmat, Lmat}; END; EXPORT LDU(m) // LDU factorization BEGIN   local Lmat, Umat, Dmat,r,c, j;   local permu, temp;   r:= rowDim(m);   c:= colDim(m);   IF (r<>c) THEN RETURN ("Matrix must be square"); END;   Pivots(m);   Dmat:= diag(L1);   temp:= lu(m);   Lmat:= temp(2);   Umat:= temp(3);   permu:= permMatrix(temp(1));   FOR j FROM 1 TO r DO   IF L1(j)<>0 THEN Umat(j):= Umat(j)/L1(j); END;    // divide U by pivots to get 1 in diagonal   END; // for   PRINT("Permutation matrix is " + temp(1));   // Return L (lower), D matrix (pivots in diagonal), U (upper) and permutation matrix   RETURN {Lmat, Dmat, Umat}; END; EXPORT permMatrix(lst)   // given a permutation list (like {1,3,2} or [1 3 2]), returns a permutation matrix BEGIN   local j, mat, dimen;   dimen:= length(lst);   mat:= MAKEMAT(0, dimen, dimen);   FOR j FROM 1 TO dimen DO   mat(j,lst(j)):= 1;   END; // for   RETURN mat; END; EXPORT echelon(m) BEGIN   // simply return ref() and RREF() for reference   RETURN {ref(m), rref(m)}; END;

∫aL√0mic (IT9CLU) :: HP Prime 50g 41CX 71b 42s 39s 35s 12C 15C - DM42, DM41X - WP34s Prime Soft. Lib
06-13-2015, 02:56 PM (This post was last modified: 06-13-2015 03:07 PM by salvomic.)
Post: #14 salvomic Senior Member Posts: 1,394 Joined: Jan 2015
RE: Pivots (Gaussian reduction)
(06-09-2015 12:50 PM)salvomic Wrote:  EDIT June 12 2015: every new version will be put in this thread of the Software Library of the Forum

I need, however two hint (help, advice):
1. the program as is works (sure it could have other bugs, however...), but only if matrix is squared (rxr) or r<c, thank in most range of the cases, but not if r>c (I know that in this case there will be one or more rows dependent, too many equations than variables), but in some cases it is stil important to get elimination.
A first hypothesis should be to add one or more columns with 0 than calculates Pivots() or LU_ext(), eliminate the columns with 0... but this is an arbitrary decision...

2. also elimination (reduction) for columns (instead of for rows) is used, sometimes: for now the program doesn't make it.
I've only a confuse idea how to revert algorithm for this purpose, in simply way...

Thank you
Salvo

∫aL√0mic (IT9CLU) :: HP Prime 50g 41CX 71b 42s 39s 35s 12C 15C - DM42, DM41X - WP34s Prime Soft. Lib
 « Next Oldest | Next Newest »

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