Pivots (Gaussian reduction)

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