Gauss Jordan utilities
06-11-2015, 10:42 PM (This post was last modified: 04-25-2016 03:13 PM by salvomic.)
Post: #1
 salvomic Senior Member Posts: 1,366 Joined: Jan 2015
Gauss Jordan utilities
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;

∫aL√0mic (IT9CLU), HP Prime 50g 41CX 71b 42s 12C 15C - DM42 WP34s :: Prime Soft. Lib
06-15-2015, 12:32 PM
Post: #2
 compsystems Senior Member Posts: 1,194 Joined: Dec 2013
RE: Gauss Jordan utilities
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 ...
06-15-2015, 12:39 PM
Post: #3
 salvomic Senior Member Posts: 1,366 Joined: Jan 2015
RE: Gauss Jordan utilities
(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

∫aL√0mic (IT9CLU), HP Prime 50g 41CX 71b 42s 12C 15C - DM42 WP34s :: Prime Soft. Lib
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)