solve()

template<class RingType>

bool DMatLinAlg< RingType >::solve	(	const Mat &	B,
		Mat &	X )

Input: B, a matrix, the right hand side of AX=B Output: X, a matrix, solution to the above returns false iff inconsistent

printf("b:\n"); debug_out_list(b, LU.numRows());

printf("y:\n"); debug_out_list(y, rk);

printf("past test for consistency\n");

buffer o; printf("after i=%ld\n", i); displayMat(o, X); printf("%s\n", o.str());

printf("x:\n"); debug_out_list(x, LU.numColumns());

buffer o; printf("after col=%ld\n", col); displayMat(o, X); printf("%s\n", o.str());

Definition at line 340 of file dmat-lu.hpp.

{
  // printf("in dmat-lu: solve\n");
  const Mat& LU = mLUObject.LUinPlace();
  // printf("in dmat-lu: after LU solve\n");
 
  // For each column of B, we solve it separately.
 
  // Step 1: Take a column, and permute it via mPerm, into b.
  // Step 2: Solve the lower triangular system Ly=b
  //  If there is no solution, then return false.
  // Step 3: Solve the upper triangular system Ux = y.
 
  // Sizes of objects here:
  //  A is m x n
  //  L is m x r
  //  U is r x n.  Here, r <= m, and A has rank r.
  //  b is a vector 0..m-1
  //  y is a vector 0..r-1
  //  x is a vector 0..n-1
 
  // printf("entering DMatLinAlg::solve\n");
  const std::vector<size_t>& pivotColumns = mLUObject.pivotColumns();
  const std::vector<size_t>& perm = mLUObject.permutation();
  size_t rk = pivotColumns.size();
 
  typename RingType::Element tmp(ring()), tmp2(ring());
 
  typedef typename RingType::ElementArray ElementArray;
  ElementArray b(ring(), LU.numRows());
  ElementArray y(ring(), rk);
  ElementArray x(ring(), LU.numColumns());
 
  X.resize(LU.numColumns(), B.numColumns());
 
  for (size_t col = 0; col < B.numColumns(); col++)
    {
      // Step 0: erase x (b will be set below, y is also set when needed).
 
      for (size_t i = 0; i < LU.numColumns(); i++) ring().set_zero(x[i]);
 
      // Step 1: set b to be the permuted i-th column of B.
      for (size_t r = 0; r < B.numRows(); r++)
        ring().set(b[r], B.entry(perm[r], col));

 
      // Step 2: Solve Ly=b
      for (size_t i = 0; i < rk; i++)
        {
          ring().set(y[i], b[i]);
          for (size_t j = 0; j < i; j++)
            {
              ring().mult(tmp, LU.entry(i, j), y[j]);
              ring().subtract(y[i], y[i], tmp);
            }
        }

 
      // Step 2B: see if the solution is consistent
      for (size_t i = rk; i < LU.numRows(); i++)
        {
          ring().set(tmp, b[i]);
          for (size_t j = 0; j < rk; j++)
            {
              ring().mult(tmp2, LU.entry(i, j), y[j]);
              ring().subtract(tmp, tmp, tmp2);
            }
          if (!ring().is_zero(tmp))
            {
              // printf("returning false\n");
              return false;
            }
        }

 
      // Step 3: Solve Ux=y
      // and place x back into X as col-th column
      for (long i = rk - 1; i >= 0; --i)
        {
          ring().set(x[i], y[i]);
          for (size_t j = i + 1; j <= rk - 1; j++)
            {
              ring().mult(tmp, LU.entry(i, pivotColumns[j]), x[j]);
              ring().subtract(x[i], x[i], tmp);
            }
          ring().divide(x[i], x[i], LU.entry(i, pivotColumns[i]));
          ring().set(X.entry(pivotColumns[i], col), x[i]);

        }

    }
 
  return true;  // The system seems to have been consistent
}

References DMat< ACoeffRing >::entry(), mLUObject, DMat< ACoeffRing >::numColumns(), DMat< ACoeffRing >::numRows(), DMat< ACoeffRing >::resize(), ring(), and x.

Referenced by inverse(), solveInvertible(), and MatrixOps::solveLinear().