compressed_row_matrix.hpp

图论必用

					vi[cur_elem] = vi[k];
	                
					wa[c] = cur_elem;
	                
					// now we've used this element...
					cur_elem++;
				}
				else
				{
					// this is a duplicate element
					//vi[wa[c]] += vi[k];
					vi[wa[c]] = b(vi[wa[c]], vi[k]);
				}
			}
	        
			// reset just the entries we used (which are stored in the ci
			// array)
			for (typename traits::nz_index_type k=*ri; k < cur_elem; k++)
			{
				wa[ci[k]] = unused;
			}

		}

		*ri = cur_elem;
	}

    /* ========================================================
     *  Routines to pack storage
     * ===================================================== */


	template <class RowIter, class ColIter, class ValIter >
	void pack_storage(compressed_row_matrix<RowIter, ColIter, ValIter>& m)
	{
		pack_storage(m, std::plus<typename std::iterator_traits<ValIter>::value_type > ());
	}

    /* ========================================================
     *  Routines to sort storage
     * ===================================================== */

	/**
	 * The version with a custom sorting operation.  
	 */
	template <class RowIter, class ColIter, class ValIter, class StrictWeakOrdering>
	void sort_storage(compressed_row_matrix<RowIter, ColIter, ValIter>& m,
		StrictWeakOrdering comp)
	{
        //BOOST_STATIC_ASSERT(0 == 1);
	}

	namespace impl
	{
		template <class ColIter, class ValIter>
		struct crm_col_val_iter_helper_type
		{
			typedef boost::tuple<
					typename std::iterator_traits<ColIter>::value_type, 
					typename std::iterator_traits<ValIter>::value_type >
				value_type;

			typedef boost::tuple<
					typename std::iterator_traits<ColIter>::value_type&, 
					typename std::iterator_traits<ValIter>::value_type& >
				ref_type;
		};

	

		template <class ColIter, class ValIter>
		class crm_col_val_iter
			: public boost::iterator_facade<
				crm_col_val_iter<ColIter, ValIter>,
				typename crm_col_val_iter_helper_type<ColIter, ValIter>::value_type,
				std::random_access_iterator_tag,
				typename crm_col_val_iter_helper_type<ColIter, ValIter>::ref_type,
				typename std::iterator_traits<ColIter>::difference_type>
		{
		public:
			crm_col_val_iter()
			{}

			crm_col_val_iter(ColIter ci, ValIter vi)
				: _ci(ci), _vi(vi)
			{
			}

			ColIter _ci;
			ValIter _vi;


		private:
			friend class boost::iterator_core_access;
            typedef typename std::iterator_traits<ColIter>::difference_type my_difference_type;


			void increment()
			{
				++_ci; ++_vi;
			}

			void decrement()
			{
				--_ci; --_vi;
			}

			bool equal(crm_col_val_iter const& other) const
			{
				return (_ci == other._ci);
			}

			typename crm_col_val_iter_helper_type<ColIter, ValIter>::ref_type dereference() const
			{
				//_tmp = boost::make_tuple(*_ci, *_vi);
				//return (_tmp);
				return (typename crm_col_val_iter_helper_type<ColIter, ValIter>::ref_type(*_ci, *_vi));
			}

			void advance(my_difference_type n)
			{
				_ci += n;
				_vi += n;
			}

			my_difference_type distance_to(crm_col_val_iter const& other) const
			{
				return ( other._ci - _ci);
			}


		};

		template <class ColIter, class ValIter>
		struct crm_col_val_iter_tuple_compare
			: public std::binary_function<
					typename crm_col_val_iter_helper_type<ColIter, ValIter>::value_type,
					typename crm_col_val_iter_helper_type<ColIter, ValIter>::value_type,
					bool>
		{
			typedef typename crm_col_val_iter_helper_type<ColIter, ValIter>::value_type T;
			bool operator()(const  T& t1, const T& t2)
			{
				return (boost::get<0>(t1) < boost::get<0>(t2));
			}
		};

		template <class ColIter, class ValIter>
		crm_col_val_iter<ColIter, ValIter>
		crm_make_col_val_iter(ColIter ci, ValIter vi)
		{
			return crm_col_val_iter<ColIter, ValIter>(ci, vi);
		};
	}


	/**
	 * Sort the storage of a matrix.  This function sorts the entries
	 * in each row by column.
	 */
	template <class RowIter, class ColIter, class ValIter>
	void sort_storage(compressed_row_matrix<RowIter, ColIter, ValIter>& m)
	{
		typedef compressed_row_matrix<RowIter, ColIter, ValIter> Matrix;
		typedef smatrix_traits<Matrix> traits;

		typedef typename traits::index_type itype;
		typedef typename traits::value_type vtype;

		RowIter ri, riend;
		ri = m._rstart;
		riend = m._rend-1;

		ColIter ci = m._cstart;
		ValIter vi = m._vstart;

		for(; ri != riend; ++ri)
		{
			typename traits::nz_index_type pos_start = *ri;
			typename traits::nz_index_type pos_end = *(ri+1);

			std::sort(
				impl::crm_make_col_val_iter(ci + pos_start, vi + pos_start),
				impl::crm_make_col_val_iter(ci + pos_end, vi + pos_end),
				impl::crm_col_val_iter_tuple_compare<ColIter,ValIter>());
		}
	}

    /* ========================================================
     *  Routines to multiply
     * ===================================================== */
	
	/**
	 * Specialize the multiply operator for this matrix type.
	 */
	template <class RowIter, class ColIter, class ValIter, class Iter1, class Iter2>
	void mult(const compressed_row_matrix<RowIter, ColIter, ValIter>& m, Iter1 x, Iter2 y)
	{
		RowIter ri = m._rstart;
		ColIter ci = m._cstart;
		ValIter vi = m._vstart;

		typedef compressed_row_matrix<RowIter, ColIter, ValIter> matrix;
		typedef typename smatrix_traits<matrix>::index_type itype;
		typedef typename smatrix_traits<matrix>::size_type stype;

		stype nr = nrows(m);

		for (itype r=0; r < nr; ++r)
		{
			typedef typename smatrix_traits<matrix>::nz_index_type nzitype;
			typedef typename smatrix_traits<matrix>::value_type vtype;

			vtype ip = 0.0;

			for (nzitype cp = ri[r]; cp < ri[r+1]; ++cp)
			{
				ip += vi[cp]*x[ci[cp]];
			}

			y[r] = ip;
		}
	}

	/**
	 * Specialize the transpose multiply operator for this matrix type.
	 */
	template <class RowIter, class ColIter, class ValIter, class Iter1, class Iter2>
	void trans_mult(const compressed_row_matrix<RowIter, ColIter, ValIter>& m, Iter1 x, Iter2 y)
	{
		RowIter ri = m._rstart;
		ColIter ci = m._cstart;
		ValIter vi = m._vstart;

		typedef compressed_row_matrix<RowIter, ColIter, ValIter> matrix;
		typedef typename smatrix_traits<matrix>::index_type itype;
		typedef typename smatrix_traits<matrix>::size_type stype;

		stype nc = ncols(m);

		// first zero the vector
		for (itype c=0; c < nc; ++c)
		{
			y[c] = 0.0;
		}

		for (itype c=0; c < nc; ++c)
		{
			typedef typename smatrix_traits<matrix>::nz_index_type nzitype;
			typedef typename smatrix_traits<matrix>::value_type vtype;

			vtype cv = x[c];

			for (nzitype cp = ri[c]; cp < ri[c+1]; ++cp)
			{
				y[ci[cp]] += vi[cp]*cv;
			}
		}
	}

    /* ========================================================
     *  Routines to find a value
     * ===================================================== */

    /**
	 * Specialize the value operator for this matrix type.
	 */
	template <class RowIter, class ColIter, class ValIter>
	typename smatrix_traits< compressed_row_matrix<RowIter, ColIter, ValIter> >::value_type value(
        typename smatrix_traits< compressed_row_matrix<RowIter, ColIter, ValIter> >::index_type row,
        typename smatrix_traits< compressed_row_matrix<RowIter, ColIter, ValIter> >::index_type col,
        const compressed_row_matrix<RowIter, ColIter, ValIter>& m)
	{
        typedef compressed_row_matrix<RowIter, ColIter, ValIter> matrix;

		RowIter ri = m._rstart;
		ColIter ci = m._cstart;
		ValIter vi = m._vstart;

        typedef typename smatrix_traits<matrix>::nz_index_type nzi_type;
        
        for (nzi_type cp = ri[row]; cp < ri[row+1]; ++cp)
        {
            if (col == ci[cp])
            {
                return (vi[cp]);
            }
        }

        return (0);
	}

    /**
     * csr_matrix is a more user manageable compressed sparse row matrix type.
     * It is based on the same ideas as the compressed_row_matrix, but designed 
     * to be less general and more specific.  
     *
     * csr_matrix uses compressed_row_matrix for some of its operations.
     *
     * The idea is that an application will use the csr_matrix structure to 
     * manage a sparse matrix and design algorithms that are NOT more 
     * generally applicable. 
     */
    /*template <class IndexType, class ValueType, class SizeType=IndexType>
    struct csr_matrix
    {
        SizeType nrows;
        SizeType ncols;
        SizeType nnz;

        IndexType* ai;

        IndexType* aj;
        ValueType* a;
    };

    template <class IndexType, class ValueType, class SizeType>
    smatrix_traits<csr_matrix<IndexType, ValueType, SizeType> >
    {
        typedef compressed_row_matrix<IndexType*, IndexType*, ValueType*>
            compressed_row_matrix_type;

        typedef IndexType index_type;
        typedef ValueType value_type;
        typedef typename compressed_row_matrix_type::nonzero_iterator nonzero_iterator;
        typedef typename compressed_row_matrix_type::nonzero_descriptor nonzero_descriptor;

	    typedef typename compressed_row_matrix_type::row_iterator row_iterator;
		
	    typedef typename compressed_row_matrix_type::row_nonzero_iterator row_nonzero_iterator;
	    typedef typename compressed_row_matrix_type::row_nonzero_descriptor row_nonzero_descriptor;

	    typedef typename SizeType nz_index_type;
		
	    typedef void column_iterator;
	    //typedef typename Mat::column_nonzero_iterator column_nonzero_iterator;

	    typedef unsymmetric_tag symmetry_category;
	    };
    }

    template <class IndexType, class ValueType, class SizeType>
    SizeType nrows(const csr_matrix<IndexType, ValueType, SizeType>& m)
    { return m.nrows; }

    template <class IndexType, class ValueType, class SizeType>
    SizeType ncols(const csr_matrix<IndexType, ValueType, SizeType>& m)
    { return m.ncols; }

    template <class IndexType, class ValueType, class SizeType>
    SizeType nnz(const csr_matrix<IndexType, ValueType, SizeType>& m)
    { return m.nnz; }*/



}

#if _MSC_VER >= 1400
	// disable the warning for deprecated c++ commands
	#pragma warning( pop )
#endif // _MSC_VER >= 1400


#endif // YASMIC_COMPRESSED_ROW_MATRIX