Pivots (Gaussian reduction)

[salvomic]

a different approach is not to use pivot() command to calc Pivots(), but use lu()...
pivot() command calculates partial pivoting element by element, it isn't a complete gaussian elimination, so I had to use the ugly "tricks" to get the job...

LU returns P (permutation, easily transformable into a matrix with permMatrix()), L (lower triangular that contains "multipliers" below the diagonal) and U (upper triangular that should have pivots, reduced, in diagonal).
I need (of Pivots() ) to get a diagonal with pivots *not* reduced (and the product of which is equal to determinant, if the matrix is not singular), as then I get LDU (with LDU() ) and LDLt (with LDLt() ), that require a D matrix with pivots in its diagonal...

I need an hint to extract multipliers from L returned by lu() and make with them the gaussian reduction (with or without permutation, given by the first list of lu())...
I thought a cycle...

i.e. if matrix is M2:= [[2,1,1], [4,-6,0], [-2,7,2]], lu(M2) returns [ [1 2 3] [[1,0,0], [-1,1,0], [2,-1,1]] [[2,1,1], [0,8,3], [0,0,1]] ]
so the permutation matrix is P = [[1,0,0], [0,01], [0,1,0]]
L is the second matrix in the list, and in it we have multipliers -1 (L(2,1)), 2 (L(3,1)), -1 (L(3,2)), so we need to subtract -1 times 1st row from 2nd, 2 times 1st row from 3rd, -1 times 2nd row from 3rd
U should have pivots already in the diagonal, but in this case it has 2,8,1, but manually we get 2, -8, 1 (the actual form of my program returns in fact Pivots(M2) -> {2,-8,1} and the product of pivots must be equal to determinant that for M2 is -16 (2*1*(-8)), while lu() returns {2,8,1} because it consider the matrix with the rows swapped as (1,3,2) and in this case the determinant *is* 16...

The cycle should extract multipliers, perform permutation (multiplying P*M2) and subtractions (following the multipliers) and return Gauss reduction and the required pivots.

Any help?
Salvo

***
EDIT 2
...but this routine (and the new cycle) only reproduces LU and returns ...simply U. I'm outside of the way...

Code:


EXPORT Pivots(m)
  // Gauss-Jordan elimination and pivots
BEGIN
  local r, c, j, k, piv, 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 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

  piv:= list2mat(L1, r);
  M1:= gj;  // Return pivot and reduce matrix with pivots in diagonal
  RETURN {piv, M1};
END;