Pivots (Gaussian reduction)

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