C++ eigen::SparseMatrix类代码示例

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

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

void SetConstraints(Eigen::SparseMatrix<float>::InnerIterator& it, int index)
{
	if (it.row() == index || it.col() == index)
	{
		it.valueRef() = it.row() == it.col() ? 1.0f : 0.0f;
	}
}

void ProbitNoise::evalModel(Eigen::SparseMatrix<double> & Ytest, const int n, Eigen::VectorXd & predictions, Eigen::VectorXd & predictions_var, const Eigen::MatrixXd &cols, const Eigen::MatrixXd &rows, double mean_rating) {
  const unsigned N = Ytest.nonZeros();
  Eigen::VectorXd pred(N);
  Eigen::VectorXd test(N);

// #pragma omp parallel for schedule(dynamic,8) reduction(+:se, se_avg) <- dark magic :)
  for (int k = 0; k < Ytest.outerSize(); ++k) {
    int idx = Ytest.outerIndexPtr()[k];
    for (Eigen::SparseMatrix<double>::InnerIterator it(Ytest,k); it; ++it) {
     pred[idx] = nCDF(cols.col(it.col()).dot(rows.col(it.row())));
     test[idx] = it.value();

      // https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Online_algorithm
      double pred_avg;
      if (n == 0) {
        pred_avg = pred[idx];
      } else {
        double delta = pred[idx] - predictions[idx];
        pred_avg = (predictions[idx] + delta / (n + 1));
        predictions_var[idx] += delta * (pred[idx] - pred_avg);
      }
      predictions[idx++] = pred_avg;

   }
  }
  auc_test_onesample = auc(pred,test);
  auc_test = auc(predictions, test);
}

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

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

IGL_INLINE void igl::find(
  const Eigen::SparseMatrix<T>& X,
  Eigen::DenseBase<DerivedI> & I,
  Eigen::DenseBase<DerivedJ> & J,
  Eigen::DenseBase<DerivedV> & V)
{
  // Resize outputs to fit nonzeros
  I.derived().resize(X.nonZeros(),1);
  J.derived().resize(X.nonZeros(),1);
  V.derived().resize(X.nonZeros(),1);

  int i = 0;
  // Iterate over outside
  for(int k=0; k<X.outerSize(); ++k)
  {
    // Iterate over inside
    for(typename Eigen::SparseMatrix<T>::InnerIterator it (X,k); it; ++it)
    {
      V(i) = it.value();
      I(i) = it.row();
      J(i) = it.col();
      i++;
    }
  }
}

void LTransform::Create_spMat_U(Eigen::SparseMatrix<double> &spMat_U){
    //*****************************
    //2015-06-29 TYPE=Notes
    //*****************************
    //求解笛卡尔坐标;2.构建控制点矩阵spMat_U;
    //*****************************
    //spMat_U=spMat_V;

    //*****************************
    qDebug() <<"START:spMat_U.insert"<< endl;
    for(int i=0;i<LMT_point.size();i++){

        spMat_U.insert(LMT_point[i].index,0)=LMT_point[i].X;
        spMat_U.insert(LMT_point[i].index,1)=LMT_point[i].Y;
        spMat_U.insert(LMT_point[i].index,2)=LMT_point[i].Z;

    }


    for (int k=0; k<spMat_V.outerSize(); ++k)
        for (Eigen::SparseMatrix<double>::InnerIterator it(spMat_V,k); it; ++it){

            if(!objMesh.is_limitP(LMT_point,it.row())){
                spMat_U.insert(it.row(),it.col())=it.value();
            }

        }
    qDebug() <<"END:spMat_U.insert"<< endl;

}//控制点坐标矩阵;

IGL_INLINE void igl::min(
  const Eigen::SparseMatrix<AType> & A,
  const int dim,
  Eigen::PlainObjectBase<DerivedB> & B,
  Eigen::PlainObjectBase<DerivedI> & I)
{
  const int n = A.cols();
  const int m = A.rows();
  B.resize(dim==1?n:m);
  B.setConstant(std::numeric_limits<typename DerivedB::Scalar>::max());
  I.resize(dim==1?n:m);
  for_each(A,[&B,&I,&dim](int i, int j,const typename DerivedB::Scalar v)
    {
      if(dim == 2)
      {
        std::swap(i,j);
      }
      // Coded as if dim == 1, assuming swap for dim == 2
      if(v < B(j))
      {
        B(j) = v;
        I(j) = i;
      }
    });
  Eigen::VectorXi Z;
  find_zero(A,dim,Z);
  for(int j = 0;j<I.size();j++)
  {
    if(Z(j) != (dim==1?m:n) && 0 < B(j))
    {
      B(j) = 0;
      I(j) = Z(j);
    }
  }
}

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


}

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::slice_into(
  const Eigen::SparseMatrix<T>& X,
  const Eigen::Matrix<int,Eigen::Dynamic,1> & R,
  const Eigen::Matrix<int,Eigen::Dynamic,1> & C,
  Eigen::SparseMatrix<T>& Y)
{

  int xm = X.rows();
  int xn = X.cols();
  assert(R.size() == xm);
  assert(C.size() == xn);
#ifndef NDEBUG
  int ym = Y.size();
  int yn = Y.size();
  assert(R.minCoeff() >= 0);
  assert(R.maxCoeff() < ym);
  assert(C.minCoeff() >= 0);
  assert(C.maxCoeff() < yn);
#endif

  // create temporary dynamic sparse matrix
  Eigen::DynamicSparseMatrix<T, Eigen::RowMajor>  dyn_Y(Y);
  // Iterate over outside
  for(int k=0; k<X.outerSize(); ++k)
  {
    // Iterate over inside
    for(typename Eigen::SparseMatrix<T>::InnerIterator it (X,k); it; ++it)
    {
      dyn_Y.coeffRef(R(it.row()),C(it.col())) = it.value();
    }
  }
  Y = Eigen::SparseMatrix<T>(dyn_Y);
}

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

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