C++ SparseMatrix类代码示例

std::vector<float> sparseVecMult(const SparseMatrix &matrix, std::vector<float> vec, CL::TinyCL &tiny){
	if (vec.size() != matrix.dim){
		std::cout << "Vector size doesn't match matrix dimensions" << std::endl;
		return std::vector<float>();
	}

	cl::Program program = tiny.loadProgram("../res/sparseMatVec.cl");
	cl::Kernel kernel = tiny.loadKernel(program, "sparseMatVec");

	int nVals = matrix.elements.size();
	int *rows = new int[nVals];
	int *cols = new int[nVals];
	float *vals = new float[nVals];
	matrix.getRaw(rows, cols, vals);

	cl::Buffer rowBuf = tiny.buffer(CL::MEM::READ_ONLY, nVals * sizeof(float), rows);
	cl::Buffer colBuf = tiny.buffer(CL::MEM::READ_ONLY, nVals * sizeof(float), cols);
	cl::Buffer valBuf = tiny.buffer(CL::MEM::READ_ONLY, nVals * sizeof(float), vals);
	cl::Buffer vecBuf = tiny.buffer(CL::MEM::READ_ONLY, vec.size() * sizeof(float), &vec[0]);
	cl::Buffer resBuf = tiny.buffer(CL::MEM::WRITE_ONLY, vec.size() * sizeof(float));

	kernel.setArg(0, sizeof(int), &nVals);
	kernel.setArg(1, rowBuf);
	kernel.setArg(2, colBuf);
	kernel.setArg(3, valBuf);
	kernel.setArg(4, vecBuf);
	kernel.setArg(5, resBuf);

	tiny.runKernel(kernel, cl::NullRange, matrix.dim);

	std::vector<float> result;
	result.resize(matrix.dim);
	tiny.readData(resBuf, matrix.dim * sizeof(float), &result[0], NULL, true);

	delete[] rows;
	delete[] cols;
	delete[] vals;

	return result;
}

void Fill( SparseMatrix<T>& A, T alpha )
{
    EL_DEBUG_CSE
    const Int m = A.Height();
    const Int n = A.Width();
    A.Resize( m, n );
    Zero( A );
    if( alpha != T(0) )
    {
        A.Reserve( m*n ); 
        for( Int i=0; i<m; ++i )
            for( Int j=0; j<n; ++j )
                A.QueueUpdate( i, j, alpha );
        A.ProcessQueues();
    }
}

void MatrixMarket( const SparseMatrix<T>& A, string basename="matrix" )
{
    EL_DEBUG_CSE
    
    string filename = basename + "." + FileExtension(MATRIX_MARKET);
    ofstream file( filename.c_str(), std::ios::binary );
    if( !file.is_open() )
        RuntimeError("Could not open ",filename);

    // Write the header
    // ================
    {
        ostringstream os;
        os << "%%MatrixMarket matrix coordinate ";
        if( IsComplex<T>::value )
            os << "complex "; 
        else
            os << "real ";
        os << "general\n";
        file << os.str();
    }
    
    // Write the size line
    // ===================
    const Int m = A.Height();
    const Int n = A.Width();
    const Int numNonzeros = A.NumEntries();
    file << BuildString(m," ",n," ",numNonzeros,"\n");
    
    // Write the entries
    // =================
    for( Int e=0; e<numNonzeros; ++e )
    {
        const Int i = A.Row(e);
        const Int j = A.Col(e);
        const T value = A.Value(e);
        if( IsComplex<T>::value )
        {
            file << 
              BuildString(i," ",j," ",RealPart(value)," ",ImagPart(value),"\n");
        }
        else
        {
            file << BuildString(i," ",j," ",RealPart(value),"\n");
        }
    }
}

double functionOfW(vector<double>& w, vector<int>& Y, SparseMatrix<double,Eigen::RowMajor>& X)
{
    size_t i;
    double sum=0.0;
    double sum2 = 0.0;
    for (i=0; i<w.size(); i++) {
        sum+=w.at(i)*w.at(i);
    }
    sum = sum*0.5;
    //cout<<"sum = "<<sum<<endl;
    int j;
    for (j=0; j<X.outerSize(); j++) {
        double dot = 0;
         for (SparseMatrix<double,Eigen::RowMajor>::InnerIterator it(X,j); it; ++it) {
            dot+=w.at(it.col())*it.value();
        }
        double temp = exp(-(Y.at(j)*dot));
        sum2+= log(1+temp);
	//cout<<"sum2 = "<<sum2<<endl;
    }
    return sum+sum2;
}

void MatrixUtils::bloc_to_sparse_copyDownTopLeftRight(const SparseBlocMatrix<ContiguousBloc<Element, Index> >& A, SparseMatrix<Element>& B)
{
	uint32 curr_row_B = B.rowdim() - 1;
	
	for(uint32 i=0; i<A.rowBlocDim (); ++i)			//for each row of blocs in A, starting at the last row in B
	{
		if(A[i].empty ())
		{
			curr_row_B -= A.bloc_height();
			continue;
		}

		if(A.isFilledWithEmptyBlocs())
				check_equal_or_raise_exception(A.FirstBlocsColumIndexes[i], 0);

		for(uint32 j=0;j<A[i].size (); ++j)		//for each bloc in the row
		{
			if(A[i][j].empty ())
			{
				continue;
			}

			uint32 bloc_idx = A.FirstBlocsColumIndexes[i] + (A.bloc_width() * j);

			for(uint16 k=0; k<A.bloc_height () && k <= curr_row_B; ++k)	//for each row in the bloc
			{
				for (uint32 l = 0; l < A[i][j].get_row_size(k); ++l)
				{
					B[curr_row_B-k].push_back(
							typename SparseMatrix<Element>::Row::value_type(
									bloc_idx + A[i][j].pos_in_row(k, l),
									A[i][j].value_in_row(k, l)));
				}
			}
		}

		curr_row_B -= A.bloc_height();
	}
}

  void makeBoundaryMatrix (SparseMatrix<double>& boundaryMatrix,
                           const int M,
                           const int N,
                           const double h,
                           const double * viscosityData) {
    #ifdef DEBUG
      cout << "<Creating " << 3 * M * N - M - N << "x" << 2 * M + 2 * N << " BoundaryMatrix>" << endl;
    #endif

    vector<Triplet<double> > tripletList;

    makeBCLaplacianXBlock (tripletList, 0,                 0,     M, N, h, viscosityData);
    makeBCLaplacianYBlock (tripletList, M * (N - 1),       2 * M, M, N, h, viscosityData);
    makeBCDivXBlock       (tripletList, 2 * M * N - M - N, 0,     M, N, h);
    makeBCDivYBlock       (tripletList, 2 * M * N - M - N, 2 * M, M, N, h);

    boundaryMatrix.setFromTriplets (tripletList.begin(), tripletList.end());

    #ifdef DEBUG
      cout << endl;
    #endif
  }

void COctaveInterface::get_real_sparsematrix(SGSparseVector<float64_t>*& matrix, int32_t& num_feat, int32_t& num_vec)
{
	const octave_value mat_feat=get_arg_increment();
	if (!mat_feat.is_sparse_type() || !(mat_feat.is_double_type()))
		SG_ERROR("Expected Sparse Double Matrix as argument %d\n", m_rhs_counter);

	SparseMatrix sm = mat_feat.sparse_matrix_value ();
	num_vec=sm.cols();
	num_feat=sm.rows();
	int64_t nnz=sm.nelem();

	matrix=new SGSparseVector<float64_t>[num_vec];

	int64_t offset=0;
	for (int32_t i=0; i<num_vec; i++)
	{
		int32_t len=sm.cidx(i+1)-sm.cidx(i);
		matrix[i].vec_index=i;
		matrix[i].num_feat_entries=len;

		if (len>0)
		{
			matrix[i].features=new SGSparseVectorEntry<float64_t>[len];

			for (int32_t j=0; j<len; j++)
			{
				matrix[i].features[j].entry=sm.data(offset);
				matrix[i].features[j].feat_index=sm.ridx(offset);
				offset++;
			}
		}
		else
			matrix[i].features=NULL;
	}
	ASSERT(offset=nnz);
}

  void makeStokesMatrix (SparseMatrix<double>& stokesMatrix,
                         const int M,
                         const int N,
                         const double h,
                         const double * viscosityData) {
    #ifdef DEBUG
      cout << "<Creating " << 3 * M * N - M - N << "x" << 3 * M * N - M - N << " stokesMatrix>" << endl;
    #endif

    vector<Triplet<double> > tripletList;

    makeLaplacianXBlock (tripletList, 0,                 0,                 M, N, h, viscosityData);
    makeLaplacianYBlock (tripletList, M * (N - 1),       M * (N - 1),       M, N, h, viscosityData);
    makeGradXBlock      (tripletList, 0,                 2 * M * N - M - N, M, N, h);
    makeGradYBlock      (tripletList, M * (N - 1),       2 * M * N - M - N, M, N, h);
    makeDivXBlock       (tripletList, 2 * M * N - M - N, 0,                 M, N, h);
    makeDivYBlock       (tripletList, 2 * M * N - M - N, M * (N - 1),       M, N, h);

    stokesMatrix.setFromTriplets (tripletList.begin(), tripletList.end());
    #ifdef DEBUG
      cout << endl;
    #endif
  }

void
Assembly::addCachedJacobian(SparseMatrix<Number> & jacobian)
{
  mooseAssert(_cached_jacobian_rows.size() == _cached_jacobian_cols.size(),
              "Error: Cached data sizes MUST be the same!");

  for(unsigned int i=0; i<_cached_jacobian_rows.size(); i++)
    jacobian.add(_cached_jacobian_rows[i], _cached_jacobian_cols[i], _cached_jacobian_values[i]);

  if (_max_cached_jacobians < _cached_jacobian_values.size())
    _max_cached_jacobians = _cached_jacobian_values.size();

  // Try to be more efficient from now on
  // The 2 is just a fudge factor to keep us from having to grow the vector during assembly
  _cached_jacobian_values.clear();
  _cached_jacobian_values.reserve(_max_cached_jacobians*2);

  _cached_jacobian_rows.clear();
  _cached_jacobian_rows.reserve(_max_cached_jacobians*2);

  _cached_jacobian_cols.clear();
  _cached_jacobian_cols.reserve(_max_cached_jacobians*2);
}

std::pair<unsigned int, Real>
EigenSparseLinearSolver<T>::adjoint_solve (SparseMatrix<T> &matrix_in,
                                           NumericVector<T> &solution_in,
                                           NumericVector<T> &rhs_in,
                                           const double tol,
                                           const unsigned int m_its)
{

  START_LOG("adjoint_solve()", "EigenSparseLinearSolver");

  libmesh_experimental();
  EigenSparseMatrix<T> mat_trans(this->comm());
  matrix_in.get_transpose(mat_trans);

  std::pair<unsigned int, Real> retval = this->solve (mat_trans,
                                                      solution_in,
                                                      rhs_in,
                                                      tol,
                                                      m_its);

  STOP_LOG("adjoint_solve()", "EigenSparseLinearSolver");

  return retval;
}

void SymmetricRuizEquil
( SparseMatrix<F>& A, 
  Matrix<Base<F>>& d,
  Int maxIter, bool progress )
{
    DEBUG_CSE
    typedef Base<F> Real;
    const Int n = A.Height();
    Ones( d, n, 1 );

    Matrix<Real> scales;
    const Int indent = PushIndent();
    for( Int iter=0; iter<maxIter; ++iter )
    {
        // Rescale the columns (and rows)
        // ------------------------------
        ColumnMaxNorms( A, scales );
        EntrywiseMap( scales, function<Real(Real)>(DampScaling<Real>) );
        EntrywiseMap( scales, function<Real(Real)>(SquareRootScaling<Real>) );
        DiagonalScale( LEFT, NORMAL, scales, d );
        SymmetricDiagonalSolve( scales, A );
    }
    SetIndent( indent );
}

void JordanCholesky( SparseMatrix<T>& A, Int n )
{
    DEBUG_ONLY(CSE cse("JordanCholesky"))
    Zeros( A, n, n );
    A.Reserve( 3*n );
    
    for( Int e=0; e<n; ++e )
    {
        if( e == 0 )
            A.QueueUpdate( e, e, T(1) );
        else
            A.QueueUpdate( e, e, T(5) );
        if( e > 0 )
            A.QueueUpdate( e, e-1, T(2) );
        if( e < n-1 )
            A.QueueUpdate( e, e+1, T(2) );
    }
    A.ProcessQueues();
}

SparseMatrix<double,0,int> Qinv2(SparseMatrix<double,0,int>& L)
{
  int n = L.rows();
  SparseMatrix<double,0,int> Sigma;
  Sigma.resize(n,n);
  Sigma.reserve(L.nonZeros()*50);
  int j;
  double Lii, val;
  for (int i=n-1; i>=0; i--)
  {
    for(SparseMatrix<double>::ReverseInnerIterator ij(L,i);ij;--ij)
    {

      j=ij.row();
      val = 0;
      SparseMatrix<double>::InnerIterator iL(L,i);
      SparseMatrix<double>::InnerIterator iS(Sigma,j);
      for(; iL; ++iL)
      {
        if(iL.row() == iL.col()){
          Lii = iL.value();
        } else {
          for(;iS.row() < iL.row(); ++iS){}
          if(iS.row() == iL.row())
            val += iL.value()*iS.value();
        }
      }
      if(i==j){
        Sigma.insert(i,j) = 1/(Lii*Lii) - val/Lii;
      } else {
        Sigma.insert(i,j) = -val/Lii;
        Sigma.insert(j,i) = -val/Lii;
      }
    }
  }
  return Sigma;
}

// Solve the symmetric system: At·A·x = At·b
bool nv::LeastSquaresSolver(const SparseMatrix & A, const FullVector & b, FullVector & x, float epsilon/*1e-5f*/)
{
    nvDebugCheck(A.width() == x.dimension());
    nvDebugCheck(A.height() == b.dimension());
    nvDebugCheck(A.height() >= A.width()); // @@ If height == width we could solve it directly...

    const uint D = A.width();

    SparseMatrix At(A.height(), A.width());
    transpose(A, At);

    FullVector Atb(D);
    //mult(Transposed, A, b, Atb);
    mult(At, b, Atb);

    SparseMatrix AtA(D);
    //mult(Transposed, A, NoTransposed, A, AtA);
    mult(At, A, AtA);

    return SymmetricSolver(AtA, Atb, x, epsilon);
}

void GetMappedDiagonal
( const SparseMatrix<T>& A,
        Matrix<S>& d,
        function<S(const T&)> func,
        Int offset )
{
    EL_DEBUG_CSE
    const Int m = A.Height();
    const Int n = A.Width();
    const T* valBuf = A.LockedValueBuffer();
    const Int* colBuf = A.LockedTargetBuffer();

    const Int iStart = Max(-offset,0);
    const Int jStart = Max( offset,0);

    const Int diagLength = El::DiagonalLength(m,n,offset);
    d.Resize( diagLength, 1 );
    Zero( d );
    S* dBuf = d.Buffer();

    for( Int k=0; k<diagLength; ++k )
    {
        const Int i = iStart + k;
        const Int j = jStart + k;
        const Int thisOff = A.RowOffset(i);
        const Int nextOff = A.RowOffset(i+1);
        auto it = std::lower_bound( colBuf+thisOff, colBuf+nextOff, j );
        if( *it == j )
        {
            const Int e = it-colBuf;
            dBuf[Min(i,j)] = func(valBuf[e]);
        }
        else
            dBuf[Min(i,j)] = func(0);
    }
}