triangularSolve() [1/2]

template<typename RT>

void MatrixOps::triangularSolve	(	DMat< RT > &	Lv,
		DMat< RT > &	x,
		int	m,
		int	strategy )

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

{
  // TODO: check rings match
  // TODO: no divide by zero
  // TODO: size of matrices
  // TODO: add tests
  // TODO: for size 0 also
  switch (strategy)
    {
      case 0:
        // TODO: change to iter
        for (size_t i = 0; i < m; i++)
          {
            auto& a = x.entry(i, 0);
            x.ring().divide(a, Lv.entry(i, m), Lv.entry(i, i));
            x.ring().negate(a, a);
            MatElementaryOps<DMat<RT>>::column_op(Lv, m, a, i);
            x.ring().negate(a, a);
          }
        break;
      case 1:
        // TODO: change to iter
        for (size_t i = 0; i < m; i++)
          {
            auto& a = x.entry(i, 0);
            x.ring().negate(a, Lv.entry(i, m));
            MatElementaryOps<DMat<RT>>::column_op(Lv, m, a, i);
            x.ring().negate(a, a);
          }
        break;
      case 2:
        // TODO: change to iter
        for (size_t i = 1; i < m + 1; i++)
          {
            auto& a = x.entry(m - i, 0);
            x.ring().divide(a, Lv.entry(m - i, m), Lv.entry(m - i, m - i));
            x.ring().negate(a, a);
            MatElementaryOps<DMat<RT>>::column_op(Lv, m, a, m - i);
            x.ring().negate(a, a);
          }
        break;
      case 3:
        // TODO: change to iter
        for (size_t i = 1; i < m + 1; i++)
          {
            auto& a = x.entry(m - i, 0);
            x.ring().negate(a, Lv.entry(m - i, m));
            MatElementaryOps<DMat<RT>>::column_op(Lv, m, a, m - i);
            x.ring().negate(a, a);
          }
        break;
    }
}

References DMat< ACoeffRing >::entry(), and x.