LDLt decomposition?

[salvomic]

(06-02-2015 11:26 AM)DrD Wrote:  Is diag() useful for you?

Code:


EXPORT gjp()
BEGIN 
  M5:=[[2,1,1,5],[-8,-2,-12,0],[1,2,0,0]];  // Your matrix
  M6:=diag(M5);  //  Diagonal Matrix
  return M6;  
END;

hi Dale,
diag() is also ok to get pivots, but after to have "pivot-ized" the original matrix: [[2,1,1,5],[4,-6,0,-2],[-2,7,2,9]]

My goal is to return a transformed matrix and a list with pivots (and also diag() could be ok).

EDIT:
Something like that, Dale:

Code:


EXPORT gaussJordan(m)
// Gauss-Jordan elimination and pivots
// Salvo Micciché 2015
BEGIN
local temp, k, gj, r, c, j, piv;
r:=rowDim(m);
c:=colDim(m);
temp:=MAKEMAT(0,r,c);
gj:=MAKEMAT(0,r,c);
piv:=MAKELIST(0,X,1,r);
temp:= m;
FOR j FROM 1 TO r DO
      temp:=pivot(temp,1,1);
      // gj(j):= temp(1);
    FOR k FROM 1 TO colDim(temp) DO
        gj(j, c-k+1):= temp(1, colDim(temp)-k+1);
    END; // inner for
      piv(j):=temp(1,1);
    IF (j<r) THEN
    temp:= delrows(temp,1);
    temp:= delcols(temp,1);
    END;
      temp:= temp/piv(j);  
END; // for
piv:= list2mat(piv, r);
RETURN {gj, piv};
END;

this runs, but with "singular" matrices it gives still division by 0 and in some cases fails (we need also swapping of some rows...), however now we have a base to calculate LDLt decomposition also...