Pivots (Gaussian reduction)

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