Pivots (Gaussian reduction)
06-07-2015, 12:41 PM (This post was last modified: 06-08-2015 10:07 AM by salvomic.)
Post: #9
 salvomic Senior Member Posts: 1,396 Joined: Jan 2015
RE: Pivots (Gaussian reduction)
(06-07-2015 12:13 PM)DrD Wrote:  In your routine: calcPivot(c,r)

FOR j FROM 1 TO r DO
M1:=pivot(M1,1,1); // Changed CAS.pivot(M1,1,1) to pivot(M1,1,1), <<clears syntax error>>

ok, changed it!

now I get right results with the matrix in some post above...

Still wrong result with M9 = [[0,2,1,1,0], [1,-1,2,0,3], [1,1,3,1,2], [3,1,8,2,6], [4,2,-6,-8,0]].
For this I put "echelon()" as function to compare with ref() and RREF(). For now Pivots(M9) gives only the first two rows and the last three with zeros...
Surely, after the row swapping pivot() doesn't offer the input that the program require...
EDIT: with the last version I get for M9 the first 3 rows: why?

More modular version, with a bit more clean code and some others //comments to help ...who want help me to debug
Corrected a bug if the first row was 0...

Code:
 calcPivot(); permMatrix(); swapPivot(); zeropivot:=0; numpivot:=1; 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("LDLt factorization", 20, 60); TEXTOUT_P("LDU (LDV) factorization", 20, 80); TEXTOUT_P("Permutation Matrix from a list", 20, 100); TEXTOUT_P("Echelon: reduced of echelon reduced form", 20, 120); TEXTOUT_P("Program created by Salvo Micciché", 20, 160, 3, RGB(0,255,0)); TEXTOUT_P("Credits: Dale (DrD), Claudio L.", 20, 180, 3, RGB(0,255,0)); WAIT; END; EXPORT Pivots(m)   // Gauss-Jordan elimination and pivots BEGIN   local r, c, piv, temp, len, j;   local permat, lumat, msquare, conf, row;   r:=rowDim(m);   c:=colDim(m);   gj:=MAKEMAT(0,r,c);   msw:=MAKEMAT(0,r,c);   M1:= m; msw:=m; // save matrix in M1 and msw for reference   L1:={0,0}; // dim list of pivots   IF M1(1,1)==0 THEN // find a non-zero pivot for first line or swap     conf:= 100; 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/ less 1st item        END; // if     END; // for     temp:= M1(1); M1(1):=M1(row); M1(row):=temp;     PRINT("Swap row " + 1 + " w/ " + row); PRINT(" ");     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   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 but not the last (product of d = det)      IF (det(permat*m)<>0 AND P<>0) THEN L1(len):= (det(permat*m) / P);       ELSE L1(len):=A; END;     gj(r,c):= L1(len);   ELSE  // matrix (r,c) not square (complete by vector b)     A:= L1(len); // last pivot     msquare:= SUB(m, {1,1}, {r,r}); // subset square for a "pseudo-determinant"     lumat:= lu(msquare);     permat:= permMatrix(lumat(1)); // permutation matrix     P:= product( L1(X),X,1,len-1 );      IF (det(permat*msquare)<>0 AND P<>0) THEN L1(len):= (det(permat*msquare) / P);       ELSE L1(len):=A; END;     gj(r,r):= L1(len); // corrected pivot copied into item r,r     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; // if trick   piv:= list2mat(L1, r);   M1:= gj;  // Return pivot and reduce matrix with pivots in diagonal   RETURN {piv, M1}; END; calcPivot() BEGIN   local j, k, r, c;   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     IF (j<r) THEN  // take minor of the matrix but last row      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; 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;   PRINT("-- Swap row " + numpivot + " w/ " + j);   END; // if    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)