C++ DMat::numColumns方法代码示例

  void subtractMultipleTo(DMat<RT>& C, 
                          const DMat<RT>& A, 
                          const DMat<RT>& B)
  // C = C - A*B
  {
    typedef typename RT::ElementType ElementType;
    typedef typename DMat<RT>::ConstIterator ConstIterator;
    
    M2_ASSERT(A.numColumns() == B.numRows());
    M2_ASSERT(A.numRows() == C.numRows());
    M2_ASSERT(B.numColumns() == C.numColumns());

    ElementType* result = C.array();

    ElementType tmp;
    A.ring().init(tmp);
    // WARNING: this routine expects the result matrix to be in ROW MAJOR ORDER
    for (size_t i = 0; i<A.numRows(); i++)
      for (size_t j = 0; j<B.numColumns(); j++)
        {
          ConstIterator i1 = A.rowBegin(i);
          ConstIterator iend = A.rowEnd(i);
          ConstIterator j1 = B.columnBegin(j);
          
          while (i1 != iend)
            {
              A.ring().mult(tmp, *i1, *j1);
              A.ring().subtract(*result, *result, tmp);
              ++i1;
              ++j1;
            }
          result++;
        }
    A.ring().clear(tmp);
  }

  void mult(const DMat<RT>& A, 
            const DMat<RT>& B, 
            DMat<RT>& result_product)
  {
    //printf("entering dmat mult\n");
    typedef typename RT::ElementType ElementType;
    typedef typename DMat<RT>::ConstIterator ConstIterator;
    
    M2_ASSERT(A.numColumns() == B.numRows());
    M2_ASSERT(A.numRows() == result_product.numRows());
    M2_ASSERT(B.numColumns() == result_product.numColumns());

    ElementType* result = result_product.array();

    ElementType tmp;
    A.ring().init(tmp);
    // WARNING: this routine expects the result matrix to be in ROW MAJOR ORDER
    for (size_t i = 0; i<A.numRows(); i++)
      for (size_t j = 0; j<B.numColumns(); j++)
        {
          ConstIterator i1 = A.rowBegin(i);
          ConstIterator iend = A.rowEnd(i);
          ConstIterator j1 = B.columnBegin(j);
          
          while (i1 != iend)
            {
              A.ring().mult(tmp, *i1, *j1);
              A.ring().add(*result, *result, tmp);
              ++i1;
              ++j1;
            }
          result++;
        }
    A.ring().clear(tmp);
  }

bool isEqual(const DMat<RT>& A, const DMat<RT>& B)
{
  assert(&A.ring() == &B.ring());
  if (B.numRows() != A.numRows()) return false;
  if (B.numColumns() != A.numColumns()) return false;
  size_t top = A.numRows() * A.numColumns();
  auto elemsA = A.array();
  auto elemsB = B.array();
  for (size_t i = 0; i < top; i++)
    if (!A.ring().is_equal(*elemsA++, *elemsB++)) return false;
  return true;
}

void transpose(const DMat<RT>& A, DMat<RT>& result)
{
  assert(&A != &result);  // these cannot be aliased!
  assert(result.numRows() == A.numColumns());
  assert(result.numColumns() == A.numRows());
  for (size_t c = 0; c < A.numColumns(); ++c)
    {
      auto i = A.columnBegin(c);
      auto j = result.rowBegin(c);
      auto end = A.columnEnd(c);
      for (; i != end; ++i, ++j) A.ring().set(*j, *i);
    }
}

void addInPlace(DMat<RT>& A, const DMat<RT>& B)
// A += B.
{
  assert(&B.ring() == &A.ring());
  assert(B.numRows() == A.numRows());
  assert(B.numColumns() == A.numColumns());

  size_t len = A.numRows() * A.numColumns();
  for (size_t i = 0; i < len; i++)
    {
      A.ring().add(A.array()[i], A.array()[i], B.array()[i]);
    }
}

bool isZero(const DMat<RT>& A)
{
  size_t len = A.numRows() * A.numColumns();
  if (len == 0) return true;
  for (auto t = A.array() + len - 1; t >= A.array(); t--)
    if (!A.ring().is_zero(*t)) return false;
  return true;
}

engine_RawArrayIntPairOrNull rawLQUPFactorizationInPlace(MutableMatrix *A, M2_bool transpose)
{
#ifdef HAVE_FFLAS_FFPACK
  // Suppose A is m x n
  // then we get A = LQUP = LSP, see e.g. http://www.ens-lyon.fr/LIP/Pub/Rapports/RR/RR2006/RR2006-28.pdf
  // P and Q are permutation info using LAPACK's convention:, see
  // http://www.netlib.org/lapack/explore-html/d0/d39/_v_a_r_i_a_n_t_s_2lu_2_r_e_c_2dgetrf_8f.html
  // P is n element permutation on column: size(P)=min(m,n); 
  // for 1 <= i <= min(m,n), col i of the  matrix was interchanged with col P(i).
  // Qt is m element permutation on rows (inverse permutation)
  // for 1 <= i <= min(m,n), col i of the  matrix was interchanged with col P(i).
  A->transpose();

  DMat<M2::ARingZZpFFPACK> *mat = A->coerce< DMat<M2::ARingZZpFFPACK> >();
  if (mat == 0) 
    {
      throw exc::engine_error("LUDivine not defined for this ring");
      //      ERROR("LUDivine not defined for this ring");
      //      return 0;
    }
  size_t nelems = mat->numColumns();
  if (mat->numRows() < mat->numColumns()) nelems = mat->numRows();
 

  std::vector<size_t> P(nelems, -1);
  std::vector<size_t> Qt(nelems, -1);
  
  // ignore return value (rank) of:
  LUdivine(mat->ring().field(),
                       FFLAS::FflasNonUnit,
                       (transpose ? FFLAS::FflasTrans : FFLAS::FflasNoTrans),
                       mat->numRows(),
                       mat->numColumns(),
                       mat->array(),
                       mat->numColumns(),
                       &P[0], 
                       &Qt[0]);

  engine_RawArrayIntPairOrNull result = new engine_RawArrayIntPair_struct;
  result->a = stdvector_to_M2_arrayint(Qt);
  result->b = stdvector_to_M2_arrayint(P);

  return result;
#endif
  return 0;
}

 inline M2_arrayintOrNull rankProfile(const DMat<RT>& A, 
                                      bool row_profile)
 {
   std::vector<size_t> profile;
   if (row_profile)
     {
       // First transpose A
       DMat<RT> B(A.ring(), A.numColumns(), A.numRows());
       MatrixOps::transpose(A,B);
       DMatLinAlg<RT> LUdecomp(B);
       LUdecomp.columnRankProfile(profile);
       return stdvector_to_M2_arrayint(profile);
     }
   else
     {
       DMatLinAlg<RT> LUdecomp(A);
       LUdecomp.columnRankProfile(profile);
       return stdvector_to_M2_arrayint(profile);
     }
 }

 DMat(const DMat<ACoeffRing>& M)
     : mRing(& M.ring())
 {
     fmpq_mat_init(mArray, M.numRows(), M.numColumns());
     fmpq_mat_set(mArray, M.mArray);
 }