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

void igl::crouzeix_raviart_massmatrix(
    const Eigen::PlainObjectBase<DerivedV> & V, 
    const Eigen::PlainObjectBase<DerivedF> & F, 
    const Eigen::PlainObjectBase<DerivedE> & E,
    const Eigen::PlainObjectBase<DerivedEMAP> & EMAP,
    Eigen::SparseMatrix<MT> & M)
{
  using namespace Eigen;
  using namespace std;
  assert(F.cols() == 3);
  // Mesh should be edge-manifold
  assert(is_edge_manifold(F));
  // number of elements (triangles)
  int m = F.rows();
  // Get triangle areas
  VectorXd TA;
  doublearea(V,F,TA);
  vector<Triplet<MT> > MIJV(3*m);
  for(int f = 0;f<m;f++)
  {
    for(int c = 0;c<3;c++)
    {
      MIJV[f+m*c] = Triplet<MT>(EMAP(f+m*c),EMAP(f+m*c),TA(f)*0.5);
    }
  }
  M.resize(E.rows(),E.rows());
  M.setFromTriplets(MIJV.begin(),MIJV.end());
}

void LaplacianOperator::computeLaplacianOperator( Eigen::SparseMatrix<double>& laplacianOperator )
{
	laplacianOperator.resize(mMeshVertexCount,mMeshVertexCount);
	laplacianOperator.reserve(Eigen::VectorXi::Constant(mMeshVertexCount,10));
	for (int i = 0; i < mMeshVertexCount; i++)
	{
		/* 如果第i个点没有邻接点，即它是一个孤立的点，那么它的laplacian坐标为0 */
		if( adjacentMatrix.innerVector(i).nonZeros() == 0)
		{
			laplacianOperator.insert(i,i) = 0;
			continue;
		}
		laplacianOperator.insert(i,i) = 1;
#ifdef MY_DEBUG
		int adjCount = 0;
#endif		
		for (Eigen::SparseMatrix<double>::InnerIterator it(adjacentMatrix,i); it; it++)
		{
			if(i != it.row())
			{
				laplacianOperator.insert(i,it.row()) = -1/degreeMatrix(i);				
#ifdef MY_DEBUG
			adjCount++;
			if(adjCount >= 10)
				printf("InnerVector size should expand! CurrentMax:%d.\n",adjCount);
#endif
			}
		}
	}
}

IGL_INLINE void igl::PolyVectorFieldFinder<DerivedV, DerivedF>::computeCoefficientLaplacian(int n, Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &D)
{
  std::vector<Eigen::Triplet<std::complex<typename DerivedV::Scalar> >> tripletList;

  // For every non-border edge
  for (unsigned eid=0; eid<numE; ++eid)
  {
    if (!isBorderEdge[eid])
    {
      int fid0 = E2F(eid,0);
      int fid1 = E2F(eid,1);

      tripletList.push_back(Eigen::Triplet<std::complex<typename DerivedV::Scalar> >(fid0,
                                           fid0,
                                           std::complex<typename DerivedV::Scalar>(1.)));
      tripletList.push_back(Eigen::Triplet<std::complex<typename DerivedV::Scalar> >(fid1,
                                           fid1,
                                           std::complex<typename DerivedV::Scalar>(1.)));
      tripletList.push_back(Eigen::Triplet<std::complex<typename DerivedV::Scalar> >(fid0,
                                           fid1,
                                                                                     -1.*std::polar(1.,-1.*n*K[eid])));
      tripletList.push_back(Eigen::Triplet<std::complex<typename DerivedV::Scalar> >(fid1,
                                           fid0,
                                                                                     -1.*std::polar(1.,1.*n*K[eid])));

    }
  }
  D.resize(numF,numF);
  D.setFromTriplets(tripletList.begin(), tripletList.end());


}

  // fmap case
  void create_matrix(const paracel::list_type<paracel::str_type> & linelst,
                     Eigen::SparseMatrix<double, Eigen::RowMajor> & blk_mtx,
                     paracel::dict_type<size_t, paracel::str_type> & rm, 
                     paracel::dict_type<size_t, paracel::str_type> & cm,
                     paracel::dict_type<size_t, int> & dm,
                     paracel::dict_type<size_t, int> & col_dm) {

    paracel::scheduler scheduler(m_comm, pattern, mix); // TODO
    // hash lines into slotslst
    auto result = scheduler.lines_organize(linelst, parserfunc);
    std::cout << "procs " << m_comm.get_rank() << " slotslst generated" << std::endl;
    m_comm.sync();
    // alltoall exchange
    auto stf = scheduler.exchange(result);
    std::cout << "procs " << m_comm.get_rank() << " get desirable lines" << std::endl;
    m_comm.sync();
    // mapping inds to ids, get rmap, cmap, std_new...
    paracel::list_type<std::tuple<size_t, size_t, double> > stf_new;
    scheduler.index_mapping(stf, stf_new, rm, cm, dm, col_dm);
    std::cout << "procs " << m_comm.get_rank() << " index mapping" << std::endl;
    // create block sparse matrix
    paracel::list_type<eigen_triple> nonzero_tpls;
    for(auto & tpl : stf_new) {
      nonzero_tpls.push_back(eigen_triple(std::get<0>(tpl), std::get<1>(tpl), std::get<2>(tpl)));
    }
    blk_mtx.resize(rm.size(), cm.size());
    blk_mtx.setFromTriplets(nonzero_tpls.begin(), nonzero_tpls.end());
  }

IGL_INLINE void igl::sparse(
  const IndexVector & I,
  const IndexVector & J,
  const ValueVector & V,
  const size_t m,
  const size_t n,
  Eigen::SparseMatrix<T>& X)
{
  using namespace std;
  using namespace Eigen;
  assert((int)I.maxCoeff() < (int)m);
  assert((int)I.minCoeff() >= 0);
  assert((int)J.maxCoeff() < (int)n);
  assert((int)J.minCoeff() >= 0);
  assert(I.size() == J.size());
  assert(J.size() == V.size());
  // Really we just need .size() to be the same, but this is safer
  assert(I.rows() == J.rows());
  assert(J.rows() == V.rows());
  assert(I.cols() == J.cols());
  assert(J.cols() == V.cols());

  vector<Triplet<T> > IJV;
  IJV.reserve(I.size());
  for(int x = 0;x<I.size();x++)
  {
    IJV.push_back(Triplet<T >(I(x),J(x),V(x)));
  }
  X.resize(m,n);
  X.setFromTriplets(IJV.begin(),IJV.end());
}

IGL_INLINE void igl::in_element(
  const Eigen::MatrixXd & V,
  const Eigen::MatrixXi & Ele,
  const Eigen::MatrixXd & Q,
  const InElementAABB & aabb,
  Eigen::SparseMatrix<double> & I)
{
  using namespace std;
  using namespace Eigen;
  using namespace igl;
  const int Qr = Q.rows();
  std::vector<Triplet<double> > IJV;
  IJV.reserve(Qr);
#pragma omp parallel for
  for(int e = 0;e<Qr;e++)
  {
    // find all
    const auto R = aabb.find(V,Ele,Q.row(e),false);
    for(const auto r : R)
    {
#pragma omp critical
      IJV.push_back(Triplet<double>(e,r,1));
    }
  }
  I.resize(Qr,Ele.rows());
  I.setFromTriplets(IJV.begin(),IJV.end());
}

void CodeAtlas::FuzzyDependBuilder::buildChildDepend( QMultiHash<QString, ChildPack>& childList , 
													 Eigen::SparseMatrix<double>& vtxEdgeMat,
													 Eigen::VectorXd&             edgeWeightVec,
													 QList<FuzzyDependAttr::DependPair>& dependPair)
{
	QStringList codeList;
	QVector<ChildPack*> childPackPtr;
	for (QMultiHash<QString, ChildPack>::Iterator pChild = childList.begin();
		pChild != childList.end(); ++pChild)
	{
		codeList.push_back(pChild.value().m_code);
		childPackPtr.push_back(&pChild.value());
	}

	std::vector<Triplet> tripletArray;
	QVector<double>		 edgeWeightArray;

	// add dependency edges between child nodes
	int ithSrc = 0;
	for (QMultiHash<QString, ChildPack>::Iterator pChild = childList.begin();
		pChild != childList.end(); ++pChild, ++ithSrc)
	{
		// for each child, find number of occurrences of this child's name
		ChildPack& srcChild = pChild.value();
		const QString& srcName = pChild.key();
		QVector<int> occurence;
		WordExtractor::findOccurrence(srcName, codeList, occurence);

		for (int ithTar = 0; ithTar < childPackPtr.size(); ++ithTar)
		{
			int nOccur = occurence[ithTar];
			if (nOccur == 0 || ithTar == ithSrc)
				continue;

			ChildPack& tarChild = *childPackPtr[ithTar];

			SymbolEdge::Ptr pEdge = SymbolTree::getOrAddEdge(srcChild.m_pNode, tarChild.m_pNode, SymbolEdge::EDGE_FUZZY_DEPEND);
			pEdge->clear();
			pEdge->strength() = nOccur;

			int nEdge = tripletArray.size()/2;
			tripletArray.push_back(Triplet(srcChild.m_index, nEdge ,1.0));
			tripletArray.push_back(Triplet(tarChild.m_index, nEdge ,-1.0));

			edgeWeightArray.push_back(nOccur);
			dependPair.push_back(FuzzyDependAttr::DependPair(srcChild.m_pNode, tarChild.m_pNode, nOccur));
		}
	}

	if (int nEdges = tripletArray.size()/2)
	{
		vtxEdgeMat.resize(childList.size(),nEdges);
		vtxEdgeMat.setFromTriplets(tripletArray.begin(), tripletArray.end());
		edgeWeightVec.resize(nEdges);
		memcpy(edgeWeightVec.data(), edgeWeightArray.data(), sizeof(double)* nEdges);
	}
}

void testEigen(int m, int n, int nnz, std::vector<int>& rows, std::vector<int>& cols,
		std::vector<double>& values, double* matB){

	double start, stop, time_to_solve, time_to_build;
    double tol=1e-9;
    Eigen::SparseMatrix<double> A;

    std::vector< Eigen::Triplet<double> > trips;
    trips.reserve(m * n);

    for (int k = 0; k < nnz; k++){
    	double _val = values[k];
    	int i = rows[k];
    	int j = cols[k];

    	if (fabs(_val) > tol){
    		trips.push_back(Eigen::Triplet<double>(i-1,j-1,_val));
        }
    }



    //NOTE: setFromTriples() accumulates contributions to the same (i,j)!
    A.resize(m, n);
    start = second();
    A.setFromTriplets(trips.begin(), trips.end());
    stop = second();
    time_to_build = stop - start;

	Eigen::SparseLU< Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int> > solverLU;



    Eigen::VectorXd b; b.resize(m);
    for (int i = 0; i < m; i++ ) b(i) = matB[i];

	printf("\nProcessing in Eigen using LU...\n");
	start = second();
	solverLU.compute(A);
	Eigen::VectorXd X = solverLU.solve(b);
	stop = second();
	time_to_solve = stop - start;

    Eigen::VectorXd ax = A * X;
    Eigen::VectorXd bMinusAx = b - ax;

	double h_r[m];
    for (int i=0; i<m; i++) h_r[i]=bMinusAx(i);

    double r_inf = vec_norminf(m, h_r);

    printf("(Eigen) |b - A*x| = %E \n", r_inf);
    printf("(Eigen) Time to build(sec): %f\n", time_to_build);
    printf("(Eigen) Time (sec): %f\n", time_to_solve);
}

inline
void
space_operator(
    Eigen::SparseMatrix<double>& result,
    Eigen::SparseMatrix<double>& laplace,
    const double multiplier,
    Eigen::SparseMatrix<double>& unit_matrix)
{
  result.resize(unit_matrix.rows(), unit_matrix.cols());
  result = laplace*multiplier+unit_matrix;
}

IGL_INLINE void igl::repdiag(
  const Eigen::SparseMatrix<T>& A,
  const int d,
  Eigen::SparseMatrix<T>& B)
{
  using namespace std;
  using namespace Eigen;
  int m = A.rows();
  int n = A.cols();

  vector<Triplet<T> > IJV;
  IJV.reserve(A.nonZeros()*d);
  // Loop outer level
  for (int k=0; k<A.outerSize(); ++k)
  {
    // loop inner level
    for (typename Eigen::SparseMatrix<T>::InnerIterator it(A,k); it; ++it)
    {
      for(int i = 0;i<d;i++)
      {
        IJV.push_back(Triplet<T>(i*m+it.row(),i*n+it.col(),it.value()));
      }
    }
  }
  B.resize(m*d,n*d);
  B.setFromTriplets(IJV.begin(),IJV.end());
  

  // Q: Why is this **Very** slow?

  //int m = A.rows();
  //int n = A.cols();

  //B.resize(m*d,n*d);
  //// Reserve enough space for new non zeros
  //B.reserve(d*A.nonZeros());

  //// loop over reps
  //for(int i=0;i<d;i++)
  //{
  //  // Loop outer level
  //  for (int k=0; k<A.outerSize(); ++k)
  //  {
  //    // loop inner level
  //    for (typename Eigen::SparseMatrix<T>::InnerIterator it(A,k); it; ++it)
  //    {
  //      B.insert(i*m+it.row(),i*n+it.col()) = it.value();
  //    }
  //  }
  //}
  //B.makeCompressed();
}

IGL_INLINE void igl::adjacency_matrix(
  const Eigen::PlainObjectBase<DerivedF> & F, 
  Eigen::SparseMatrix<T>& A)
{
  using namespace std;
  using namespace Eigen;
  typedef typename DerivedF::Scalar Index;

  typedef Triplet<T> IJV;
  vector<IJV > ijv;
  ijv.reserve(F.size()*2);
  // Loop over faces
  for(int i = 0;i<F.rows();i++)
  {
    // Loop over this face
    for(int j = 0;j<F.cols();j++)
    {
      // Get indices of edge: s --> d
      Index s = F(i,j);
      Index d = F(i,(j+1)%F.cols());
      ijv.push_back(IJV(s,d,1));
      ijv.push_back(IJV(d,s,1));
    }
  }

  const Index n = F.maxCoeff()+1;
  A.resize(n,n);
  switch(F.cols())
  {
    case 3:
      A.reserve(6*(F.maxCoeff()+1));
      break;
    case 4:
      A.reserve(26*(F.maxCoeff()+1));
      break;
  }
  A.setFromTriplets(ijv.begin(),ijv.end());

  // Force all non-zeros to be one

  // Iterate over outside
  for(int k=0; k<A.outerSize(); ++k)
  {
    // Iterate over inside
    for(typename Eigen::SparseMatrix<T>::InnerIterator it (A,k); it; ++it)
    {
      assert(it.value() != 0);
      A.coeffRef(it.row(),it.col()) = 1;
    }
  }
}

IGL_INLINE bool igl::GeneralPolyVectorFieldFinder<DerivedV, DerivedF>::
                     solve(const Eigen::VectorXi &isConstrained,
                           const Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic> &cfW,
                           const Eigen::VectorXi &rootsIndex,
                           Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic> &output)
{

  // polynomial is of the form:
  // z^(2n) +
  // -c[0]z^(2n-1) +
  // c[1]z^(2n-2) +
  // -c[2]z^(2n-3) +
  // ... +
  // (-1)^n c[n-1]

  std::vector<Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1>> coeffs(n,Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1>::Zero(numF, 1));

  for (int i =0; i<n; ++i)
  {
    int degree = i+1;

    Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic,1> Ck;
    getGeneralCoeffConstraints(isConstrained,
                               cfW,
                               i,
                               rootsIndex,
                               Ck);

    Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > DD;
    computeCoefficientLaplacian(degree, DD);
    Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > f; f.resize(numF,1);

    if (isConstrained.sum() == numF)
      coeffs[i] = Ck;
    else
      minQuadWithKnownMini(DD, f, isConstrained, Ck, coeffs[i]);
  }

  std::vector<Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 2> > pv;
  setFieldFromGeneralCoefficients(coeffs, pv);

  output.setZero(numF,3*n);
  for (int fi=0; fi<numF; ++fi)
  {
    const Eigen::Matrix<typename DerivedV::Scalar, 1, 3> &b1 = B1.row(fi);
    const Eigen::Matrix<typename DerivedV::Scalar, 1, 3> &b2 = B2.row(fi);
    for (int i=0; i<n; ++i)
      output.block(fi,3*i, 1, 3) = pv[i](fi,0)*b1 + pv[i](fi,1)*b2;
  }
  return true;
}

void Mesh::buildMassMatrix(const VectorXd &q, Eigen::SparseMatrix<double> &M) const
{
    M.resize(numdofs(), numdofs());
    vector<Tr> entries;
    for(OMMesh::VertexIter vi = mesh_->vertices_begin(); vi != mesh_->vertices_end(); ++vi)
    {
        int vidx = vi.handle().idx();
        double area = barycentricDualArea(q, vidx);
        for(int i=0; i<3; i++)
            entries.push_back(Tr(3*vidx+i, 3*vidx+i, area));
    }

    M.setFromTriplets(entries.begin(), entries.end());
}

inline void igl::PlanarizerShapeUp<DerivedV, DerivedF>::assembleSelector(int fi,
                                                                            Eigen::SparseMatrix<typename DerivedV::Scalar > &S)
{
  
  std::vector<Eigen::Triplet<typename DerivedV::Scalar>> tripletList;
  for (int fvi = 0; fvi< ni; fvi++)
  {
    int vi = Fin(fi,fvi);
    tripletList.push_back(Eigen::Triplet<typename DerivedV::Scalar>(3*fvi+0,3*vi+0,1.));
    tripletList.push_back(Eigen::Triplet<typename DerivedV::Scalar>(3*fvi+1,3*vi+1,1.));
    tripletList.push_back(Eigen::Triplet<typename DerivedV::Scalar>(3*fvi+2,3*vi+2,1.));
  }
  
  S.resize(3*ni,3*numV);
  S.setFromTriplets(tripletList.begin(), tripletList.end());
  
}

Eigen::SparseMatrix<double> Condi2Joint(Eigen::SparseMatrix<double> Condi, Eigen::SparseVector<double> Pa)
{	// second dimension of Condi is the parent
	Eigen::SparseMatrix<double> Joint;
	Joint.resize(Condi.rows(), Condi.cols());

	for (int cols = 0; cols < Condi.cols(); cols++)
	{
		Eigen::SparseVector<double> tmp_vec = Condi.block(0, cols, Condi.rows(), 1)*Pa.coeff(cols);
		for (int id_rows = 0; id_rows < tmp_vec.size(); id_rows++)
		{
			Joint.coeffRef(id_rows, cols) = tmp_vec.coeff(id_rows);
		}

	}
	Joint.prune(TOLERANCE);
	return Joint;

}