computeLU() [5/5]

template<class RingType>

void DMatLUinPlace< RingType >::computeLU ( )

private

Definition at line 230 of file dmat-lu-inplace.hpp.

{
  if (mIsDone) return;
 
  //  std::cout << "computing LU decomposition generic version" << std::endl;
  typename RingType::Element tmp(mLU.ring());
 
  size_t col = 0;  // current column we are working on
  size_t row = 0;  // current row we are working on
  size_t nrows = mLU.numRows();
  size_t ncols = mLU.numColumns();
 
  while (col < ncols && row < nrows)
    {
      // printf("*** in naive row,col = (%ld, %ld) ***\n", row, col);
      // debug_out();
 
      // Step 1: Set the 'upper' values: (row,col)..(nrows-1,col)
      for (size_t r = row; r < nrows; r++)
        {
          for (size_t i = 0; i < row; i++)
            {
              mLU.ring().mult(tmp, mLU.entry(r, i), mLU.entry(i, col));
              mLU.ring().subtract(mLU.entry(r, col), mLU.entry(r, col), tmp);
            }
        }
 
      // printf("after step 1\n");
      // debug_out();
 
      // Step 2: Find a pivot among the elements in step 1.
      //  If one: swap rows if needed
      //  If none, increment 'col', and continue at start of loop
      size_t k = findPivot(row, col);
      if (k == static_cast<size_t>(-1))
        {
          col = col + 1;
          continue;
        }
      // printf("pivot is in row %ld\n", k);
      std::swap(mPerm[row], mPerm[k]);
      if (k != row)
        {
          MatElementaryOps<Mat>::interchange_rows(mLU, k, row);
          mSign = !mSign;
        }
      mPivotColumns.push_back(col);
      const ElementType& pivot = mLU.entry(row, col);
 
      // printf("after step 2:\n");
      // debug_out();
 
      // Step 3A: Set the 'upper' elements in (row,col+1), ..., (row,ncols-1).
      for (size_t c = col + 1; c < ncols; c++)
        {
          for (size_t i = 0; i < row; i++)
            {
              mLU.ring().mult(tmp, mLU.entry(row, i), mLU.entry(i, c));
              mLU.ring().subtract(mLU.entry(row, c), mLU.entry(row, c), tmp);
            }
        }
      // printf("after step 3A:\n");
      // debug_out();
 
      // Step 3B: Set the 'lower' elements in (row+1,row), ..., (nrows-1,row)
      //  from (row+1,col), ..., (nrows-1,col)
      // This just means dividing then by the pivot
      // except, if we have skipped columns for pivots, we must set these
      // elements
      // in column 'row', not 'col'...
      // Step 3C: if row != col, then set these elements to 0:
      //  (row+1,col), ..., (nrows-1,col)
      for (size_t r = row + 1; r < nrows; r++)
        {
          mLU.ring().divide(mLU.entry(r, row), mLU.entry(r, col), pivot);
          if (row < col) ring().set_zero(mLU.entry(r, col));
        }
 
      // printf("after step 3B:\n");
      // debug_out();
 
      row++;
      col++;
    }
 
  mIsDone = true;
}

References findPivot(), mIsDone, mLU, mPerm, mPivotColumns, mSign, ring(), and std::swap().

Referenced by LUinPlace().