Pivots (Gaussian reduction)
06-05-2015, 05:06 PM (This post was last modified: 06-05-2015 10:47 PM by salvomic.)
Post: #2
 salvomic Senior Member Posts: 1,396 Joined: Jan 2015
RE: Pivots (Gaussian elimination)
A version a bit more "advanced". Now it can handle swapping for rows if there is problem in finding solution (*if* the solution there is, the matrix isn't singular, and so on...).

But I've still not found a real *general* formula: help very appreciated!

When it runs well, the program now returns pivots(), LDU() and LDLt() factorization, and a routine to calculate permutation matrix also.
In this version there isn't still a menu...

The code (sorry for some //comments for debugging):
Code:
 calcPivot(); zeropivot:=0; numpivot:=1; gj:=[0,0]; EXPORT gaussJordan() // 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;   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) THEN   M1:=m;   // M1:= CAS.SWAPROW(M1, numpivot, (numpivot+1));    if (numpivot==r) THEN numpivot:=1; END;  // if pivot is 0 swap rows (if last swap w/ 1)   temp:= M1(numpivot); M1(numpivot):=M1(numpivot+1); M1(numpivot+1):=temp;   calcPivot(c,r);    END; // if   piv:= list2mat(L1, r);   M1:= gj;   // Return pivot, matrix reduced and RREF form   RETURN {piv, M1, ref(m)}; END; calcPivot(c,r) BEGIN   local j, k;   zeropivot:=0;   FOR j FROM 1 TO r DO     M1:=CAS.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   // division or not?   M1:= M1 ./ L1(j); // M2, M3, M4 (also IF (j<r))  //  M1(1):= M1(1) ./ L1(j); // M3, M4, M6   // IF (j>1) THEN M1:= M1 ./ L1(j-1); END; // M3, M4, M7   // without div  M3, M4, M7   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) tranpose 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}) and the matrix dimension, 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;

The examples (see at the end of the program are):
M2: [[2,1,1], [4,-6,0], [-2,7,2]] with pivots: {2,-8,1} -> U [[2,1,1], [0,-8,2], [0,0,1]]
M3: [[1,1,1], [1,2,-2], [2,1,-3]] with pivots: {1,1,8} -> [[1,1,1], [0,1,-3], [0,0,-8]]
M4: [[0,3,-1], [1,1,2], [1,4,1]] with pivots: {1,3,0} -> [[1,1,2], [0,3,-1], [0,0,0]] (this one will swap two times; try also another version: [[1,1,2], [0,3,-1], [1,4,1]]: swap two times, with inverted pivots: without the function for swap in the program it'll give "error: invalid input"...): it's a singular matrix (determinant 0, as the product of pivots)
M6: [[2,-1,0], [-1,2,-1], [0,-1,2]] with pivots: {2, 3/2, 4/3} -> U [[2, -1, 0], [0, 3/3, -1], [0,0, 4/3]]
M7: [[1,2,3], [2,4,3], [3,2,-1]] with pivots: {1, -4, -3} -> [[1,2,3], [0,-4,-10], [0, 0, -3]]
M8: [[1,1,1], [1,1,3], [2,5,8]] with pivots: {1,3,2} -> U [[1,1,1], [0,3,6], [0,0,2]] (with swapping row 2 and 3)

And now the problem of the "non general" formula (please, let me understand where I'm wrong):
Code:
   // M1:= M1./L1(j); // M2, M3, M4 (also IF (j<r))   M1(1):= M1(1)./L1(j); // M3, M4, M6   // IF (j>1) THEN M1:= M1./L1(j-1); END; // M3, M4, M7   // without div  M3, M4, M7
M2 needs division for current pivot in every lines, M6 only in the first, the others (M3, M4, M7, M8)) don't need division (apparently)...

But I'm searching for *one* formula, not more ;-)

Thank you for your patience and effort!

Salvo

EDIT: thanks Dale (DrD) and Claudio L. for some tips!

∫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)