C++ SpMat::reserve方法代码示例

inline void eliminationOperators(SpMat& A, DynamicVector<size_t>& Cset, size_t fnode,
        DynamicVector<double>& q, SpMat& P, size_t& P_col, size_t P_row) {
    /* Inizializziamo la riga P_row con A.nonZeros(fnode) - 1 non nulli*/
    double scalingFactor = 1.0;
    /* Riserviamo spazio in ciascuna colonna di P per un adeguato numero
     * di elementi. Il -1 c'è perché (il reciproco del)l'elemento in
     * diagonale lo mettiamo in q
     */
    P.reserve(P_row, A.nonZeros(fnode) - 1);
    /* Per ciascuna f-riga prendo gli elementi in ogni c-colonna */
    DynamicVector<size_t>::Iterator ccol = Cset.begin();
    for (SpMat::Iterator frow = A.begin(fnode); frow != A.end(fnode); ++frow) {
        if (frow->index() == fnode) { //Elemento della diagonale
            q[P_row] = (1.0 / frow->value()); //Q ha elementi pari al numero di righe di R
            scalingFactor = -(q[P_row]);
        } else if (ccol != Cset.end()) {
            break; //Non ha senso andare avanti se abbiamo finito i ccol
        } else if (frow->index() == (*ccol)) {
            /* Elemento fuori della diagonale ed è anche un c-colonna */
            P.append(P_row, P_col, frow->value());
            P_col++;
            ccol++;
        }
    }

    P.finalize(P_row); //Finalizziamo la riga P_row

    /* Non dimentichiamo di scalare gli elementi della riga corrente
     * per  -scalingFactor */
    for (SpMat::Iterator it = P.begin(P_row); it != P.end(P_row); it++)
        it->value() *= -scalingFactor;
}

inline
void
spglue_kron::apply_noalias(SpMat<eT>& out, const SpMat<eT>& A, const SpMat<eT>& B)
  {
  arma_extra_debug_sigprint();
  
  const uword A_n_rows = A.n_rows;
  const uword A_n_cols = A.n_cols;
  
  const uword B_n_rows = B.n_rows;
  const uword B_n_cols = B.n_cols;
  
  const uword out_n_nonzero = A.n_nonzero * B.n_nonzero;
  
  out.reserve(A_n_rows * B_n_rows, A_n_cols * B_n_cols, out_n_nonzero);
  
  if(out_n_nonzero == 0)  { return; }
  
  access::rw(out.col_ptrs[0]) = 0;
  
  uword count = 0;
  
  for(uword A_col=0; A_col < A_n_cols; ++A_col)
  for(uword B_col=0; B_col < B_n_cols; ++B_col)
    {
    for(uword A_i = A.col_ptrs[A_col]; A_i < A.col_ptrs[A_col+1]; ++A_i)
      {
      const uword out_row = A.row_indices[A_i] * B_n_rows;
      
      const eT A_val = A.values[A_i];
      
      for(uword B_i = B.col_ptrs[B_col]; B_i < B.col_ptrs[B_col+1]; ++B_i)
        {
        access::rw(out.values[count])      = A_val * B.values[B_i];
        access::rw(out.row_indices[count]) = out_row + B.row_indices[B_i];
        
        count++;
        }
      }
    
    access::rw(out.col_ptrs[A_col * B_n_cols + B_col + 1]) = count;
    }
  }

// The Real Core Function doing the actual mesh processing.
bool FilterHarmonicPlugin::applyFilter(QAction * action, MeshDocument & md, RichParameterSet & par, vcg::CallBackPos * cb)
{
	switch(ID(action))
	{
	case FP_SCALAR_HARMONIC_FIELD :
	{
		typedef vcg::GridStaticPtr<CMeshO::VertexType, CMeshO::ScalarType> VertexGrid;

		typedef double                           CoeffScalar; // TODO, when moving the code to a class make it a template (CoeffScalar = double)

		typedef CMeshO::ScalarType               ScalarType;
		typedef CMeshO::CoordType                CoordType;
		typedef CMeshO::VertexType               VertexType;
		typedef CMeshO::FaceType                 FaceType;

		typedef Eigen::Triplet<CoeffScalar>      T;
		typedef Eigen::SparseMatrix<CoeffScalar> SpMat; //sparse matrix type of double


		CMeshO & m = md.mm()->cm;

		vcg::tri::Allocator<CMeshO>::CompactFaceVector(m);
		vcg::tri::Allocator<CMeshO>::CompactVertexVector(m);

		md.mm()->updateDataMask(MeshModel::MM_FACEFACETOPO | MeshModel::MM_VERTMARK);
		vcg::tri::UpdateBounding<CMeshO>::Box(m);
		vcg::tri::UpdateTopology<CMeshO>::FaceFace(m);

		int n  = m.VN();
		int fn = m.FN();

		std::vector<T>             coeffs; // coefficients of the system
		std::map<size_t,CoeffScalar> sums; // row sum of the coefficient

		SpMat laplaceMat; // the system to be solved
		laplaceMat.resize(n, n);

		Log("Generating coefficients.`");
		cb(0, "Generating coefficients...");
		vcg::tri::UpdateFlags<CMeshO>::FaceClearV(m);
		// Iterate over the faces
		for (size_t i = 0; i < m.face.size(); ++i)
		{
			CMeshO::FaceType & f = m.face[i];

			if (f.IsD())
			{
				assert(int(i) == fn);
				break; // TODO FIX the indexing of vertices
			}

			assert(!f.IsV());
			f.SetV();

			// Generate coefficients for each edge
			for (int idx = 0; idx < 3; ++idx)
			{
				CoeffScalar weight;
				WeightInfo res = ComputeWeight<FaceType, CoeffScalar>(f, idx, weight);

				switch (res)
				{
				case EdgeAlreadyVisited : continue;
				case Success            : break;
				case BorderEdge         :
					this->errorMessage = "Mesh not closed, cannot compute harmonic field on mesh containing holes or borders";
					return false;
				default: assert(0);
				}

//				if (weight < 0) weight = 0; // TODO check if negative weight may be an issue

				// Add the weight to the coefficients vector for both the vertices of the considered edge
				size_t v0_idx = vcg::tri::Index(m, f.V0(idx));
				size_t v1_idx = vcg::tri::Index(m, f.V1(idx));

				coeffs.push_back(T(v0_idx, v1_idx, -weight));
				coeffs.push_back(T(v1_idx, v0_idx, -weight));

				// Add the weight to the row sum
				sums[v0_idx] += weight;
				sums[v1_idx] += weight;
			}

			f.SetV();
		}

		// Fill the system matrix
		Log("Fill the system matrix");
		cb(10, "Filling the system matrix...");
		laplaceMat.reserve(coeffs.size());
		for (std::map<size_t,CoeffScalar>::const_iterator it = sums.begin(); it != sums.end(); ++it)
		{
			coeffs.push_back(T(it->first, it->first, it->second));
		}
		laplaceMat.setFromTriplets(coeffs.begin(), coeffs.end());

		// Get the two vertices with value set
		VertexGrid vg;
//.........这里部分代码省略.........

arma_hot
inline
void
spglue_minus::apply_noalias(SpMat<eT>& out, const SpProxy<T1>& pa, const SpProxy<T2>& pb)
  {
  arma_extra_debug_sigprint();
  
  arma_debug_assert_same_size(pa.get_n_rows(), pa.get_n_cols(), pb.get_n_rows(), pb.get_n_cols(), "subtraction");
  
  if(pa.get_n_nonzero() == 0)  { out = pb.Q; out *= eT(-1); return; }
  if(pb.get_n_nonzero() == 0)  { out = pa.Q;                return; }
  
  const uword max_n_nonzero = spglue_elem_helper::max_n_nonzero_plus(pa, pb);
  
  // Resize memory to upper bound
  out.reserve(pa.get_n_rows(), pa.get_n_cols(), max_n_nonzero);
  
  // Now iterate across both matrices.
  typename SpProxy<T1>::const_iterator_type x_it  = pa.begin();
  typename SpProxy<T1>::const_iterator_type x_end = pa.end();
  
  typename SpProxy<T2>::const_iterator_type y_it  = pb.begin();
  typename SpProxy<T2>::const_iterator_type y_end = pb.end();
  
  uword count = 0;
  
  while( (x_it != x_end) || (y_it != y_end) )
    {
    eT out_val;
    
    const uword x_it_row = x_it.row();
    const uword x_it_col = x_it.col();
    
    const uword y_it_row = y_it.row();
    const uword y_it_col = y_it.col();
    
    bool use_y_loc = false;
    
    if(x_it == y_it)
      {
      out_val = (*x_it) - (*y_it);
      
      ++x_it;
      ++y_it;
      }
    else
      {
      if((x_it_col < y_it_col) || ((x_it_col == y_it_col) && (x_it_row < y_it_row))) // if y is closer to the end
        {
        out_val = (*x_it);
        
        ++x_it;
        }
      else
        {
        out_val = -(*y_it);  // take the negative
        
        ++y_it;
        
        use_y_loc = true;
        }
      }
    
    if(out_val != eT(0))
      {
      access::rw(out.values[count]) = out_val;
      
      const uword out_row = (use_y_loc == false) ? x_it_row : y_it_row;
      const uword out_col = (use_y_loc == false) ? x_it_col : y_it_col;
      
      access::rw(out.row_indices[count]) = out_row;
      access::rw(out.col_ptrs[out_col + 1])++;
      ++count;
      }
    }
  
  const uword out_n_cols = out.n_cols;
  
  uword* col_ptrs = access::rwp(out.col_ptrs);
  
  // Fix column pointers to be cumulative.
  for(uword c = 1; c <= out_n_cols; ++c)
    {
    col_ptrs[c] += col_ptrs[c - 1];
    }
  
  if(count < max_n_nonzero)
    {
    if(count <= (max_n_nonzero/2))
      {
      out.mem_resize(count);
      }
    else
      {
      // quick resize without reallocating memory and copying data
      access::rw(         out.n_nonzero) = count;
      access::rw(     out.values[count]) = eT(0);
      access::rw(out.row_indices[count]) = uword(0);
      }
    }
//.........这里部分代码省略.........