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

inline
void
spop_resize::apply(SpMat<typename T1::elem_type>& out, const SpOp<T1, spop_resize>& in)
  {
  arma_extra_debug_sigprint();
  
  out = in.m;
  
  out.resize(in.aux_uword_a, in.aux_uword_b);
  }

void Assembler::operKernel(EOExpr<A> oper,const MeshHandler<ORDER>& mesh,
	                     FiniteElement<Integrator, ORDER>& fe, SpMat& OpMat)
{
	std::vector<coeff> triplets;


  	for(auto t=0; t<mesh.num_triangles(); t++)
  	{
		fe.updateElement(mesh.getTriangle(t));

		// Vector of vertices indices (link local to global indexing system)
		std::vector<UInt> identifiers;
		identifiers.resize(ORDER*3);
		for( auto q=0; q<ORDER*3; q++)
			identifiers[q]=mesh.getTriangle(t)[q].id();


		//localM=localMassMatrix(currentelem);
		for(int i = 0; i < 3*ORDER; i++)
		{
			for(int j = 0; j < 3*ORDER; j++)
			{
				Real s=0;

				for(int l = 0;l < Integrator::NNODES; l++)
				{
					s += oper(fe,i,j,l) * fe.getDet() * fe.getAreaReference()* Integrator::WEIGHTS[l];//(*)
				}
			  triplets.push_back(coeff(identifiers[i],identifiers[j],s));
			}
		}

	}

  	UInt nnodes = mesh.num_nodes();
  	OpMat.resize(nnodes, nnodes);
	OpMat.setFromTriplets(triplets.begin(),triplets.end());
	//cout<<"done!"<<endl;;
}

// 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;
//.........这里部分代码省略.........

void ADMMCut::CalculateQ(const VX _D, SpMat &Q)
{
	// Construct Hessian Matrix for D-Qp problem
	Q.resize(6 * Ns_, Nd_);
	vector<Eigen::Triplet<double>> Q_list;

	for (int i = 0; i < Ns_; i++)
	{
		int u = ptr_dualgraph_->v_orig_id(i);
		WF_edge *edge = ptr_frame_->GetNeighborEdge(u);

		while (edge != NULL)
		{
			int e_id = ptr_dualgraph_->e_dual_id(edge->ID());
			if (e_id != -1)
			{
				int v = edge->pvert_->ID();
				int j = ptr_dualgraph_->v_dual_id(v);

				MX eKuu = ptr_stiff_->eKv(edge->ID());
				MX eKeu = ptr_stiff_->eKe(edge->ID());
				VX Fe = ptr_stiff_->Fe(edge->ID());
				VX Di(6);
				VX Dj(6);

				if (i < Ns_ && j < Ns_)
				{
					for (int k = 0; k < 6; k++)
					{
						Di[k] = _D[6 * i + k];
						Dj[k] = _D[6 * j + k];
					}
				}
				else
				{
					if (i < Ns_)
					{
						for (int k = 0; k < 6; k++)
						{
							Di[k] = _D[6 * i + k];
							Dj[k] = 0;
						}
					}

					if (j < Ns_)
					{
						for (int k = 0; k < 6; k++)
						{
							Di[k] = 0;
							Dj[k] = _D[6 * j + k];
						}
					}
				}
				VX Gamma = eKuu * Di + eKeu * Dj - Fe;

				for (int k = 0; k < 6; k++)
				{
					Q_list.push_back(Triplet<double>(6 * i + k, e_id, Gamma[k]));
				}
			}

			edge = edge->pnext_;
		}
	}

	Q.setFromTriplets(Q_list.begin(), Q_list.end());
}

LevelElimination * LAMGLSSolver::coarsenElimination(const Matrix & finerMatrix) {
cout << "coarsenElimination()" << endl;
    //Per salvare i P ed i q
    std::vector<qPStage*> cStages;
    /* Prealloco spazio per al più SETUP_ELIMINATION_MAX_STAGES così evito
     * possibili riallocazioni */
    cStages.reserve(SETUP_ELIMINATION_MAX_STAGES);

    //Copio la matrice così posso modificarla senza modificare quella del finerLevel
    Matrix A(finerMatrix);
    Matrix Acc; //Per lo Schur Complement System
    Matrix Acf; //Per lo Schur Complement System
    // Il numero di stage di eliminazione eseguiti sulla matrice A
    size_t stageNum = 0;

    /* Vector che vengono riusati*/
    DynamicVector<size_t> Cset; //Insieme dei nodi da non eliminare
    DynamicVector<size_t> Fset; //Insieme dei nodi da eliminare
    DynamicVector<size_t> degree; //Grado di ciascun nodo
    DynamicVector<size_t> candidate; //Nodi candidati all'eliminazione
    DynamicVector<NodeEliminationStatus> status; //Stato dei nodi
    /* P è sempre sparsa anche quando A è densa */
    SpMat P;
    DynamicVector<double> q;

    while (stageNum < SETUP_ELIMINATION_MAX_STAGES) {

        size_t A_rows = A.rows();
        if (A_rows <= MAX_DIRECT_SOLVE_SIZE)
            break;
        /* (1) Calcola il vettore degree della matrice. */
        degree.resize(A_rows);

        for (size_t i = 0; i < A_rows; ++i) {
            degree[i] = A.nonZeros(i);
            if (A(i, i) != 0)//Tolgo l'elemento sulla diagonale
                degree[i]--;
            else
                cerr << "La diagionale i-esima c'ha lelemento a zero!!! ERRORE MORTALE!!!" << endl;
        }

        /* (2) [f c ] = lowDegreeNodes(A,degree,MaxDegree)  */
        candidate.resize(A_rows);
        candidate = 0;
        /* Individuo i nodi candidati (degree[i] <= MAX_DEGREE) */
        size_t cnnZ = 0; //Devo usare questa variabili perché i valori "settati" di candidate comprendono il valore "0" cioè il nodo 0
        for (size_t i = 0; i < A_rows; ++i) {
            if (degree[i] <= SETUP_ELIMINATION_MAX_DEGREE) {
                candidate[cnnZ++] = i;
            }
        }
        /* I primi cnnZ elementi sono stati inizializzati */
        candidate.resize(cnnZ, true);

        status.resize(A_rows);
        status = HIGH_DEGREE; //Tutti gli elementi prendono HIGH_DEGREE
        //    status(candidate)  = 0; % Reset all relevant nodes to "not visited"
        for (size_t i = 0; i < candidate.size(); ++i) {
            status[candidate[i]] = NOT_DECIDED;
        }

        for (size_t k = 0; k < candidate.size(); ++k) {
            lowdegreesweep(A, candidate[k], status); //Template call
        }

        /* Adesso devo creare i vettori F e C, inserendovi i nodi che verranno
         * o meno eliminati */
        size_t nf = 0; //|Fset|
        size_t nc = 0; //|Cset|
        Cset.resize(A_rows, false);
        Fset.resize(A_rows, false);
        for (size_t i = 0; i < A_rows; ++i) {
            if (status[i] == LOW_DEGREE)
                Fset[nf++] = i; //Inserisco il nodo i nell'insieme F
            else
                Cset[nc++] = i; //Lo inserisco invece in C
        }

        /* L'insieme C non può mai essere vuoto, dobbiamo lasciargli almeno
         * un elemento
         */
        if (nc == 0) {
            Cset[nc++] = Fset[--nf];
            Fset[nf] = 0;

        }
        Cset.resize(nc, true); //nc non è mai 0
        Fset.resize(nf, true);

        /* FINE di (2) [f c ] = lowDegreeNodes(A,degree,MaxDegree) */


        if ((nf <= SETUP_ELIMINATION_MIN_ELIM_FRACTION * A_rows)) {
            /* Il coarsening non è abbastanza efficace perché andiamo ad eliminare
             * un numero, nf, di nodi che è inferiore alla minima soglia accettabile
             * di eliminazione. Ci fermiamo senza eliminare.*/
            break;
        }


//.........这里部分代码省略.........