+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Software Libraries (/forum-10.html)
+--- Forum: HP Prime Software Library (/forum-15.html)
+--- Thread: Gauss Jordan utilities (/thread-4136.html)
Gauss Jordan utilities - salvomic - 06-11-2015 10:42 PM
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;RE: Gauss Jordan utilities - compsystems - 06-15-2015 12:32 PM
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
...RE: Gauss Jordan utilities - salvomic - 06-15-2015 12:39 PM
(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