Pivots (Gaussian reduction) - salvomic - 06-03-2015 01:20 PM
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'd like to receive your approach, tips, advice...
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...
RE: Pivots (Gaussian elimination) - salvomic - 06-05-2015 05:06 PM
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!
RE: Pivots (Gaussian elimination) - DrD - 06-05-2015 07:03 PM
(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."
RE: Pivots (Gaussian elimination) - salvomic - 06-05-2015 07:29 PM
(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
RE: Pivots (Gaussian reduction) - salvomic - 06-06-2015 09:41 PM
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;
RE: Pivots (Gaussian reduction) - DrD - 06-07-2015 10:42 AM
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-
RE: Pivots (Gaussian reduction) - salvomic - 06-07-2015 11:28 AM
(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;
RE: Pivots (Gaussian reduction) - DrD - 06-07-2015 12:13 PM
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>>
RE: Pivots (Gaussian reduction) - salvomic - 06-07-2015 12:41 PM
(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;
RE: Pivots (Gaussian reduction) - salvomic - 06-08-2015 04:02 PM
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;
RE: Pivots (Gaussian reduction) - salvomic - 06-08-2015 10:15 PM
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;
RE: Pivots (Gaussian reduction) - DrD - 06-09-2015 10:55 AM
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.
RE: Pivots (Gaussian reduction) - salvomic - 06-09-2015 12:50 PM
(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;
RE: Pivots (Gaussian reduction) - salvomic - 06-13-2015 02:56 PM
(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
|