kernel()

template<class RingType>

size_t DMatLinAlg< RingType >::kernel

(

Mat &

X

)

Output: X, a matrix, columns form a basis of Ax=0 returns dim of nullspace

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

{
  const Mat& LU = mLUObject.LUinPlace();
  const std::vector<size_t>& pivotColumns = mLUObject.pivotColumns();
 
  typename RingType::Element tmp(ring()), tmp2(ring());
 
  // First, let's set X to be the correct size:
  X.resize(LU.numColumns(), LU.numColumns() - pivotColumns.size());
 
  size_t col = 0;
  size_t nextpivotidx = 0;
  size_t nextpivotcol =
      (pivotColumns.size() > 0 ? pivotColumns[0] : LU.numColumns());
  size_t colX = 0;
  while (colX < X.numColumns())
    {
      if (col == nextpivotcol)
        {
          col++;
          nextpivotidx++;
          nextpivotcol =
              (nextpivotidx < pivotColumns.size() ? pivotColumns[nextpivotidx]
                                                  : LU.numColumns());
          continue;
        }
      // At this point, we are ready to create a column of X.
      ring().set_from_long(X.entry(col, colX), -1);
      // Now we loop through and set the elements in the rows of X = pivot
      // columns.
      for (long p = nextpivotidx - 1; p >= 0; p--)
        {
          // set X.entry(pivotColumns[p], colX)
          ring().set(tmp, LU.entry(p, col));
          for (size_t i = nextpivotidx - 1; i >= p + 1; i--)
            {
              ring().mult(tmp2,
                          LU.entry(p, pivotColumns[i]),
                          X.entry(pivotColumns[i], colX));
              ring().subtract(tmp, tmp, tmp2);
            }
          ring().divide(tmp, tmp, LU.entry(p, pivotColumns[p]));
          ring().set(X.entry(pivotColumns[p], colX), tmp);
        }
      colX++;
      col++;
    }
  return X.numColumns();
}

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

Referenced by MatrixOps::nullSpace().