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

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

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());


}

void place::erodeSparse(const Eigen::SparseMatrix<double> &src,
                        Eigen::SparseMatrix<double> &dst) {
  dst = Eigen::SparseMatrix<double>(src.rows(), src.cols());
  std::vector<Eigen::Triplet<double>> tripletList;
  Eigen::MatrixXd srcNS = Eigen::MatrixXd(src);

  for (int i = 0; i < srcNS.cols(); ++i) {
    for (int j = 0; j < srcNS.rows(); ++j) {
      double value = 0.0;
      for (int k = -1; k < 1; ++k) {
        for (int l = -1; l < 1; ++l) {
          if (i + k < 0 || i + k >= srcNS.cols() || j + l < 0 ||
              j + l >= srcNS.rows())
            continue;
          else
            value = std::max(value, srcNS(j + l, i + k));
        }
      }

      if (value != 0)
        tripletList.push_back(Eigen::Triplet<double>(j, i, value));
    }
  }
  dst.setFromTriplets(tripletList.begin(), tripletList.end());
}

ZSparseMatrix Assembler2D::getLaplacianMatrix() {
    typedef Eigen::Triplet<double> T;
    std::vector<T> triplets;

    for (size_t i=0; i<m_mesh->getNbrElements(); ++i)
    {
        VectorI idx = m_mesh->getElement(i);
        assert(idx.size() == 3);

        Eigen::MatrixXd& dN = m_DN[i];

        // Small strain-displacement matrix
        //
        Eigen::MatrixXd B(2,3);
        B << dN(0,0), dN(1,0), dN(2,0),
        dN(0,1), dN(1,1), dN(2,1);

        Eigen::MatrixXd k_el = B.transpose() * B * m_mesh->getElementVolume(i);

        for (size_t j=0; j<3; ++j)
            for (size_t k=0; k<3; ++k)
                triplets.push_back(T(idx[j], idx[k], k_el(j,k)));
    }

    Eigen::SparseMatrix<double> L = Eigen::SparseMatrix<double>(m_mesh->getNbrNodes(), m_mesh->getNbrNodes());
    L.setFromTriplets(triplets.begin(), triplets.end());
    return ZSparseMatrix(L);
}

ZSparseMatrix Assembler2D::getDisplacementStrainMatrix() {
    typedef Eigen::Triplet<double> T;
    std::vector<T> triplets;

    for (size_t i=0; i<m_mesh->getNbrElements(); i++) {
        Eigen::MatrixXd dN = m_DN[i];
        VectorI idx = m_mesh->getElement(i);
        assert(idx.size() == 3);
        double V = m_mesh->getElementVolume(i);

        // e_xx
        size_t row = i * 3;
        for (size_t k=0; k<3; k++) {
            triplets.push_back(T(row, idx[k]*2, dN(k,0)));
        }

        // e_yy
        row++;
        for (size_t k=0; k<3; k++) {
            triplets.push_back(T(row, idx[k]*2+1, dN(k,1)));
        }

        // e_xy
        row++;
        for (size_t k=0; k<3; k++) {
            triplets.push_back(T(row, idx[k]*2  , dN(k,1) / 2.0));
            triplets.push_back(T(row, idx[k]*2+1, dN(k,0) / 2.0));
        }
    }

    Eigen::SparseMatrix<double> B = Eigen::SparseMatrix<double>(3*m_mesh->getNbrElements(), 2*m_mesh->getNbrNodes());
    B.setFromTriplets(triplets.begin(), triplets.end());
    return ZSparseMatrix(B);
}

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

  // 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());
  }

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 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 Mesh::buildAreaMatrix(Eigen::SparseMatrix<double>& A) const
{
    std::vector<Eigen::Triplet<double>> ATriplet;
    
    for (VertexCIter v = vertices.begin(); v != vertices.end(); v++) {
        ATriplet.push_back(Eigen::Triplet<double>(v->index, v->index, v->dualArea()));
    }
    
    A.setFromTriplets(ATriplet.begin(), ATriplet.end());
}

void MathUtil::symmetricMatrixFromTriplets(const std::vector<T> &triplets, Eigen::SparseMatrix<double> &M)
{
    vector<T> augtriplets = triplets;
    for(vector<T>::const_iterator it = triplets.begin(); it != triplets.end(); ++it)
    {
        if(it->row() < it->col())
            augtriplets.push_back(T(it->col(), it->row(), it->value()));
    }
    M.setFromTriplets(augtriplets.begin(), augtriplets.end());
}

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

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

inline
void
identity_operator(
    Eigen::SparseMatrix<double>& result,
    int size)
{
  result.resize(size, size);
  result.reserve(size);
  std::vector<double_triplet_t> matrix_coeffs;
  matrix_coeffs.reserve(size); // #diag
  // diag
  for (int i = 0; i < size; ++i) matrix_coeffs.push_back(double_triplet_t(i,i, 1.));
  result.setFromTriplets(matrix_coeffs.begin(), matrix_coeffs.end());
}