Pivots (Gaussian reduction)

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