Pivots (Gaussian reduction)

[salvomic]

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...

[salvomic]

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!

[DrD]

(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."

[salvomic]

(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

[salvomic]

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;

[DrD]

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-

[salvomic]

(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;

[DrD]

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>>

[salvomic]

(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;

[salvomic]

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;

[salvomic]

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;

[DrD]

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.

[salvomic]

(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;

[salvomic]

(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