Pivots (Gaussian reduction)

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