Pivots (Gaussian reduction)
06-07-2015, 11:28 AM (This post was last modified: 06-07-2015 12:55 PM by salvomic.)
Post: #7
 salvomic Senior Member Posts: 1,396 Joined: Jan 2015
RE: Pivots (Gaussian reduction)
(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;

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

 Messages In This Thread Pivots (Gaussian reduction) - salvomic - 06-03-2015, 01:20 PM RE: Pivots (Gaussian elimination) - salvomic - 06-05-2015, 05:06 PM RE: Pivots (Gaussian elimination) - DrD - 06-05-2015, 07:03 PM RE: Pivots (Gaussian elimination) - salvomic - 06-05-2015, 07:29 PM RE: Pivots (Gaussian reduction) - salvomic - 06-06-2015, 09:41 PM RE: Pivots (Gaussian reduction) - DrD - 06-07-2015, 10:42 AM RE: Pivots (Gaussian reduction) - salvomic - 06-07-2015 11:28 AM RE: Pivots (Gaussian reduction) - DrD - 06-07-2015, 12:13 PM RE: Pivots (Gaussian reduction) - salvomic - 06-07-2015, 12:41 PM RE: Pivots (Gaussian reduction) - salvomic - 06-08-2015, 04:02 PM RE: Pivots (Gaussian reduction) - salvomic - 06-08-2015, 10:15 PM RE: Pivots (Gaussian reduction) - DrD - 06-09-2015, 10:55 AM RE: Pivots (Gaussian reduction) - salvomic - 06-09-2015, 12:50 PM RE: Pivots (Gaussian reduction) - salvomic - 06-13-2015, 02:56 PM

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