Pivots (Gaussian reduction)

[salvomic]

(06-09-2015 10:55 AM)DrD Wrote:  Recognizing that Gauss Jordan elimination can be accomplished in many different ways, ...

Another decomposition for LU(M2) is:

{[[1,0,0],[0.5,1,0],[−0.5,1,1]],[[4,−6,0],[0,4,1],[0,0,1]],[[0,1,0],[1,0,0],[0,0,1]]}, where the matrices are {L,U,P}. Checking: L*U/P ==M2.

There may be another update released before long, and hopefully resolve some of the issues.

yes, I know, in fact I could stop with the LU (and REF) that already Prime has, but I search also for others way, to control some results by books and for study purpose...
My first attempt is no completely sure, as it uses two "tricks", the second attempt give no more nor less U as the LU() command, reproducing it in another way

I made programs more difficult in less time than this one, but this is fascinating me, in some way...

An yes, I control always with Wolfram and with the reverse formula...

At the start I thought that the command pivot() gave the complete reduction, but I can see that applying it more time could get exactly LU or REF or ...RREF.

Anyway, my purpose is to get LDU (LDV) and LDLt, and for now my first attempt is the nearest to them. Maybe.

Salvo

***
The last code with utilities: Pivots, LU_ext (LU also if the matrix is not square), LDLt, LDU, permMatrix, echelon (and a few bugs corrected) is here. It works much better than first and it is conform to some fonts, like the books of Gilbert Strang and a few others, obviously not with *every* matrices...

I promise don't bore you more with this program...

EDIT June 12 2015: every new version will be put in this thread of the Software Library of the Forum

Code:


calcPivot();
permMatrix();
swapPivot();
zeropivot:=0;
numpivot:=1;
swap:= 0;
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("LU extended: for not square mat", 20, 60);
TEXTOUT_P("LDLt factorization", 20, 80);
TEXTOUT_P("LDU (LDV) factorization", 20, 100);
TEXTOUT_P("Permutation Matrix from a list", 20, 120);
TEXTOUT_P("Echelon: reduced of echelon reduced form", 20, 140);
TEXTOUT_P("Program created by Salvo Micciché", 20, 180, 3, RGB(0,255,0));
TEXTOUT_P("Credits: Dale (DrD), Claudio L.", 20, 200, 3, RGB(0,255,0));
WAIT;
END;

EXPORT Pivots(m)
  // Gauss-Jordan elimination and pivots
BEGIN
  local r, c, temp, len, j;
  local permat, lumat, msquare, conf, row;
  r:=rowDim(m);
  c:=colDim(m);
  L1:=MAKELIST(0,X,1,r); // dim list of pivots
  gj:=MAKEMAT(0,r,c);
  msw:=MAKEMAT(0,r,c);
  M1:= m; msw:=m; // save matrix in M1 and msw for reference
  swap:= 0;
  IF M1(1,1)==0 THEN // find a non-zero pivot for first line or swap
    conf:= 1E10; 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/ lesser 1st item 
      END; // if
    END; // for
    temp:= M1(1); M1(1):=M1(row); M1(row):=temp;
    // PRINT; PRINT("Swap row " + 1 + " w/ " + row);
    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
    FOR j FROM 2 TO length(L1) DO // divide by pivots if swapping
   IF L1(j)<>0 THEN gj(j):= gj(j) * L1(j-1)^(-1); END; 
    END; // for
  END; // if
  M1:= gj;
  IF (r==c) THEN L1:= diag(M1); ELSE temp:=SUB(M1,{1,1}, {r,r}); L1:= diag(temp); END;

   len:= length(L1);  // trick 1 to get correct last pivot (product of pivot is determinant)
   P:= product(L1(X),X,1,len-1); // product of pivots but not the last (product of d = det)
   msquare:= SUB(m, {1,1}, {r,r}); // square subset  ("pseudo-determinant" if mat isn't square)
   lumat:= lu(msquare);
   permat:= permMatrix(lumat(1)); // permutation matrix
   IF (det(msquare)<>0 AND P<>0) THEN 
    IF swap==0 THEN L1(len):= det(msquare)/P; ELSE L1(len):= det(permat*msquare)/P; END;
   END; // if
   A:=EVAL(gj(r,r));
   gj(r,r):= L1(len); // corrected pivot copied into item r,r
   IF (r<>c) THEN // mat isn't square
    FOR j FROM r+1 TO c DO
      IF (gj(r,j)<>0) THEN gj(r,j):= gj(r,j)*gj(r,r)/A; END;
    END; // for
   END; // if - end trick

   M1:= gj;  // Return pivots in diagonal and U (upper with pivots)
   IF (r==c) THEN L1:= diag(M1); ELSE temp:=SUB(M1,{1,1}, {r,r}); L1:= diag(temp); END;
   RETURN {diag(L1), M1};
END;

calcPivot()
BEGIN
  local j, k, r, c, temp;
  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 // trick 2 to get echelon item for last but one row
      temp:= IFTE(M1(1,colDim(M1))<>0, M1(1,colDim(M1)), M1(2,colDim(M1))); // correct pseudo-pivot
      gj(r-1,c):= IFTE(temp<>0, temp/gcd(temp), temp);
      temp:= IFTE(M1(1,colDim(M1)-1)<>0, M1(1,colDim(M1)-1), M1(2,colDim(M1)-1));
      gj(r-1,c-1):= IFTE(temp<>0, temp/gcd(temp), temp);
    END; // trick2

  END; // for j
  FOR j FROM 2 TO length(L1) DO    // to divide rows for current pivot 
  IF L1(j-1)<>0 THEN gj(j):= gj(j) * L1(j-1)^(-1); 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;
      swap:= 1;
      // PRINT("-- Swap row " + numpivot + " w/ " + j);
  END; // if
  END; // for 
END;

EXPORT LU_ext(m)
  // Gauss-Jordan elimination and pivots: LU also for rectangular matrices
BEGIN
  local r, c, j, k, temp, plist;
  local permat, Lmat, Umat, Pmat;
  r:=rowDim(m);
  c:=colDim(m);
  L1:=MAKELIST(0,X,1,r); // dim list of pivots
  gj:=MAKEMAT(0,r,c); // building matrix
  gj:= m;
  M1:= SUB(m, {1,1},{r,r}); // here M1 is the square subset of matrix
  temp:= lu(M1);
  plist:= temp(1);
  Lmat:= temp(2);
  Umat:= temp(3);
  Pmat:= permMatrix(plist);
  M1:= Pmat*M1;
  gj:= Pmat*gj;

  FOR j FROM 1 TO r DO  // calc multiplier from Lmat and subtract rows
  FOR k FROM 1 TO j DO
    IF (j<>k) THEN
    IF (Lmat(j,k)<>0) THEN gj(j):= gj(j) - Lmat(j,k)*gj(k); END;
    END; // if
  END; // inner for
  END; // outer for
  FOR j FROM 1 TO r DO // extract pivots
  FOR k FROM 1 TO c DO
    IF (gj(j,k)<>0) THEN L1(j):= gj(j,k); BREAK; END;
  END; // inner for
  END; //outer for
  // adapt last two rows if echelon
  IF (gj(r,c-1)<>0 AND gj(r-1, c-1)<>0 AND gj(r-1, c-2)==0) THEN
    temp:= gj(r,c-1); gj(r):= gj(r) - temp / gj(r-1,c-1)*gj(r-1); 
    Lmat(r,r-1):= temp / gj(r-1,c-1);
    L1(length(L1)):= 0;
  END; 
  PRINT("Permutation matrix is " + plist);
  M1:= gj;  // Return L, U and matrix with pivots in diagonal
  RETURN {diag(L1), Lmat, M1};
END;

EXPORT LDLt(m)
  // Factorizazion LDLt, matrix m must be simmetric: A=At
BEGIN
  local Lmat, LmatT, Dmat, r,c, j;
  r:= rowDim(m);
  c:= colDim(m);
  IF (r<>c) THEN RETURN ("Matrix must be square and symmetric"); END;
  Pivots(m);
  Dmat:= diag(L1);
  Lmat:= M1;
  FOR j FROM 1 TO r DO
  IF L1(j)<>0 THEN Lmat(j):= Lmat(j)/L1(j); END; 
  // 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,c, j;
  local permu, temp;
  r:= rowDim(m);
  c:= colDim(m);
  IF (r<>c) THEN RETURN ("Matrix must be square"); END;
  Pivots(m);
  Dmat:= diag(L1);
  temp:= lu(m);
  Lmat:= temp(2);
  Umat:= temp(3);
  permu:= permMatrix(temp(1));
  FOR j FROM 1 TO r DO
  IF L1(j)<>0 THEN Umat(j):= Umat(j)/L1(j); END; 
  // 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};
END;

EXPORT permMatrix(lst)
  // given a permutation list (like {1,3,2} or [1 3 2]), returns a permutation matrix
BEGIN
  local j, 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;