LUincremental() [1/2]

template<typename RT>

M2_arrayintOrNull MatrixOps::LUincremental	(	std::vector< size_t > &	P,
		DMat< RT > &	LU,
		const DMat< RT > &	v,
		int	m )

Definition at line 567 of file mat-linalg.hpp.

{
  size_t n = LU.numRows();
 
  // copy permuted v to m-th column of LU
  // TODO: change to iter
  for (size_t j = 0; j < n; j++)
    LU.ring().set(LU.entry(j, m), v.entry(P[j], 0));
 
  // reduce the m-th column of LU and forward solve
  DMat<RT> x{LU.ring(), n, 1};
  triangularSolve(LU, x, m, 1);
  // place solution of forward solve in U
  // TODO: change to iter
  for (size_t i = 0; i < m; i++)
    LU.ring().set(LU.entry(i, m), x.entry(i, 0));
 
  // look for a pivot in L
  int pivotPosition = -1;
  // TODO: change to iter
  for (size_t j = m; j < n; j++)
    if (!LU.ring().is_zero(LU.entry(j, m)))
      {
        pivotPosition = j;
        break;
      }
  // if no pivot found, return
  if (pivotPosition == -1) return stdvector_to_M2_arrayint(P);
  // otherwise swap rows and update P
  MatElementaryOps<DMat<RT>>::interchange_rows(LU, pivotPosition, m);
  std::swap(P[pivotPosition], P[m]);
 
  // scale column of L
  // TODO: change to iter
  for (int j = m + 1; j < n; j++)
    LU.ring().divide(LU.entry(j, m), LU.entry(j, m), LU.entry(m, m));
 
  return stdvector_to_M2_arrayint(P);
}

References DMat< ACoeffRing >::entry(), LU(), stdvector_to_M2_arrayint(), std::swap(), triangularSolve(), and x.