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

  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 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;
}

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 addTo(DMat<RT>& A, MatrixWindow wA, const DMat<RT>& B, MatrixWindow wB)
{
  assert(wA.sameSize(wB));
  long rA = wA.begin_row;
  long rB = wB.begin_row;
  for (; rA < wA.end_row; ++rA, ++rB)
    {
      long cA = wA.begin_column;
      long cB = wB.begin_column;
      for (; cA < wA.end_column; ++cA, ++cB)
        {
          auto& a = A.entry(rA, cA);
          A.ring().add(a, a, B.entry(rB, cB));
        }
    }
}

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())
 {
   fmpz_mat_init_set(mArray, M.mArray);
 }

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

 void determinantRRR(const DMat<M2::ARingRRR> &M, M2::ARingRRR::elem &result) { M.ring().set_zero(result); }

double ResF4toM2Interface::setDegreeZeroMap(SchreyerFrame& C,
                                       DMat<RingType>& result,
                                       int slanted_degree,
                                       int lev)
// 'result' should be previously initialized, but will be resized.
// return value: -1 means (slanted_degree, lev) is out of range, and the zero matrix was returned.
//   otherwise: the fraction of non-zero elements is returned.
{
  // As above, get the size of the matrix, and 'newcols'
  // Now we loop through the elements of degree 'slanted_degree + lev' at level 'lev'
  const RingType& R = result.ring();
  if (not (lev > 0 and lev <= C.maxLevel()))
    {
      result.resize(0,0);
      return -1;
    }
  assert(lev > 0 and lev <= C.maxLevel());
  int degree = slanted_degree + lev;
  auto& thislevel = C.level(lev);
  int ncols = 0;
  for (auto p=thislevel.begin(); p != thislevel.end(); ++p)
    {
      if (p->mDegree == degree) ncols++;
    }

  auto& prevlevel = C.level(lev-1);
  int* newcomps = new int[prevlevel.size()];
  int nrows = 0;
  for (int i=0; i<prevlevel.size(); i++)
    if (prevlevel[i].mDegree == degree)
      newcomps[i] = nrows++;
    else
      newcomps[i] = -1;

  result.resize(nrows, ncols);

  int col = 0;
  long nnonzeros = 0;  
  for (auto p=thislevel.begin(); p != thislevel.end(); ++p)
    {
      if (p->mDegree != degree) continue;
      auto& f = p->mSyzygy;
      auto end = poly_iter(C.ring(), f, 1);
      auto i = poly_iter(C.ring(), f);
      for ( ; i != end; ++i)
        {
          long comp = C.monoid().get_component(i.monomial());
          if (newcomps[comp] >= 0)
            {
              R.set_from_long(result.entry(newcomps[comp], col), C.gausser().coeff_to_int(i.coefficient()));
              nnonzeros++;
            }
        }
      ++col;
    }
  double frac_nonzero = (nrows*ncols);
  frac_nonzero = static_cast<double>(nnonzeros) / frac_nonzero;

  delete[] newcomps;

  return frac_nonzero;
}

 DMat(const DMat<ACoeffRing>& M)
   : mRing(& M.ring())
 {
   fq_nmod_mat_init_set(mArray, M.mArray, ring().flintContext());
 }