Gauss Jordan utilities

[salvomic]

hi,
I put here one program of mine (you can use it with credits to me and who collaborated with me and to this forum), useful to make operations with gaussian reduction (elimination), extending the tools that the Prime has already (LU, pivot, ref, RREF, ...).
The program isn't still "perfect", and I'll put there every new amended version as soon as they are ready.
The program has for now these routines:
gaussJordanUtil() : the about with credits
Pivots(m) : a personal algorithm (I'm working still on it) to calculate pivots, using the command pivot() of the Prime, that does a "partial pivoting". It returns a matrix with pivots in diagonal and upper matrix (U)
LU_ext(m) : an extension of the lu() command, to treat also a non square matrix and echelon forms, returning results similar to Pivots(): matrix with pivots, lower matrix (L), upper matrix (U)
LDLt(m) : a routine to perform a transposition of a symmetric and square matrix into L (lower), D (matrix with pivots in diagonal) and U (upper), that's also a transpose of L (so the name L-D-Lt)
LDU(m) : a factorization in L (lower), D (pivots in diagonal) and U (upper) for a square matrix
echelon(m) : simply returns both ref(m) and RREF(m), echelon forms
permMatrix(list) : giving as parameter a list {} or a matrix [] with the number of row to permute, it gives the matrix of permutation.

Enjoy and (or) help to make a better version!

Salvo Micciché (salvomic)

The code:

Code:


calcPivot();
permMatrix();
swapPivot();
zeropivot:=0;
numpivot:=1;
swap:= 0;
gj:=[0,0];
msw:=[0,0];


EXPORT gaussJordanUtil()
// 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;

[compsystems]

Very good application, but the idea is to group your work with directories,
LINEAR ALGEBRA
...
but unfortunately the hpprime still not supported

Code:


DIR LINEAR ALGEBRA 
....
ENDDIR

DIR ....
....
ENDDIR

...

[salvomic]

(06-15-2015 12:32 PM)compsystems Wrote:  Very good application, but the idea is to group your work with directories,
LINEAR ALGEBRA
...
but unfortunately the hpprime still not supported

...

thank you for compliments!
I hope to develop better the application and I put in my agenda also your advice, hoping to could do something soon...

Still I've to develop calculation of a matrix where m is greater than n and few other things...

Salvo