//
//      CoordMatrix.h   define the matrix stored in a coordinate format
//

#ifndef COORDMATRIX_H
#define COORDMATRIX_H

#include<vector>

#include <assert.h>

using namespace std;


class CoordMatrix
	{
	private:
       // val strore the values in matrix;
		vector<double> val;
	   
		// row is row index of value;
		vector<int> row;

		// col is colume index of value ;
	   vector<int> col;
		
	   // nonZeor is the total number of nonzero value in the matrix;
	    int nonZero;

		// numRow is the dimension of the matrix  generall we use square matrix;
		int numRow;
		
	public:
        // default constructor 
		CoordMatrix():nonZero(0),numRow(0) {val.resize(0); row.resize(0); col.resize(0); row.resize(0); col.resize(0); val.resize(0);};

		// constructor from a two D array
		CoordMatrix(double **aVal)
			{  int R=sizeof aVal /sizeof aVal[0];
				int C=column;
			    for (int i=0; i<R; i++)	{
			    for(int j=0; j<C; j++){
	         	if( *(aVal[i]+j) != 0.0 ) 	{row.push_back(i); col.push_back(j); val.push_back( *(aVal[i]+j) ); }
	 					                      }
			                                      }
				nonZero=row.size();

				numRow=R;
			}
}

Dear friends:
I wrote a spare Matrix class and the Matrix was stored in a coornidate format. Only the row index, colume index and the value of sparse matrix was stored. I want to constructe a CoordMatrix from a 2d array, but the code is wrong , i do not how to transfer a 2d of arbitrary size to the constructor. Could you please give me some advice.
regards

> I wrote a spare Matrix class and the Matrix was stored in a coornidate format.

Consider using a data structure that would allow efficient look up of the element at a given row and col position. One way is to use a std::map<> where the (row,col) tuple is the key.


> if( *(aVal[i]+j) != 0.0 ) This check for a floating point value being non-zero may not give you the results that you anticipate.
See: http://www.parashift.com/c++-faq-lite/newbie.html#faq-29.17


> I want to constructe a CoordMatrix from a 2d array

There is no "standard" 2d array in C or C++. There are a variety of ways of simulating it. So the answer is: it depends. In every case, information about the number of rows and cols would have to be passed to a function/constructor.
For more information: http://forums.devx.com/showthread.php?p=519222

If the array has a static storage duration, the number of elements are constants known at compile time and we could do it this way (assuming a std::map<> implementation of the sparse 2d matrix):

#include <iostream>
#include <cmath>
#include <map>
#include <stdexcept>
#include <iomanip>

struct sparse_matrix_2d
{
    struct index ;

    template< std::size_t ROWS, std::size_t COLS >
        explicit inline sparse_matrix_2d( double (&array)[ROWS][COLS] )
            : nrows(ROWS), ncols(COLS)
        {
            for( std::size_t i = 0 ; i < ROWS ; ++i )
                for( std::size_t j = 0 ; j < COLS ; ++j )
                    if( is_non_zero( array[i][j] ) )
                        non_zero_elements[index(i,j)] = array[i][j] ;
        }

    double at( std::size_t row, std::size_t col ) const
    {
        map_t::const_iterator iter = non_zero_elements.find( index(row,col) ) ;
        return iter != non_zero_elements.end() ? iter->second : 0 ;
    }

    void at( std::size_t row, std::size_t col, double value )
    {
        if( (row>=nrows) || (col>=ncols) ) throw std::out_of_range("invalid index") ;
        else if( is_non_zero(value) ) non_zero_elements[ index(row,col) ] = value ;
    }

    inline friend std::ostream& operator<< ( std::ostream& stm, const sparse_matrix_2d& m )
    {
            for( std::size_t i = 0 ; i < m.nrows ; ++i )
            {
                for( std::size_t j = 0 ; j < m.ncols ; ++j )
                      stm << std::setw(10) << m.at(i,j) ;
                stm << '\n' ;
            }
            return stm ;
    }

    struct index
    {
        index( std::size_t r, std::size_t c ) : row(r), col(c) {}
        const std::size_t row ;
        const std::size_t col ;
        bool operator< ( const index& that ) const
        {
          if( row < that.row ) return true ;
          else if( row == that.row ) return col < that.col ;
          else return false ;
        }
    };

    private:
        std::size_t nrows ;
        std::size_t ncols ;
        typedef std::map< index, double > map_t ;
        map_t non_zero_elements ;
        static const double EPSILON ;
        static inline bool is_non_zero( double d ) { return std::abs(d) > EPSILON ; }
};

const double sparse_matrix_2d::EPSILON = 1.e-8 ;

int main()
{
    double d[3][4] = { {1,0,0,2}, {0,0,3,0}, {4,0,5,0} } ;
    sparse_matrix_2d mtx(d) ;
    std::cout << std::fixed << std::setprecision(2) << std::showpoint << mtx << '\n' ;
}

> I wrote a spare Matrix class and the Matrix was stored in a coornidate format.

Consider using a data structure that would allow efficient look up of the element at a given row and col position. One way is to use a std::map<> where the (row,col) tuple is the key.


> if( *(aVal[i]+j) != 0.0 ) This check for a floating point value being non-zero may not give you the results that you anticipate.
See: http://www.parashift.com/c++-faq-lite/newbie.html#faq-29.17


> I want to constructe a CoordMatrix from a 2d array

There is no "standard" 2d array in C or C++. There are a variety of ways of simulating it. So the answer is: it depends. In every case, information about the number of rows and cols would have to be passed to a function/constructor.
For more information: http://forums.devx.com/showthread.php?p=519222

If the array has a static storage duration, the number of elements are constants known at compile time and we could do it this way (assuming a std::map<> implementation of the sparse 2d matrix):

#include <iostream>
#include <cmath>
#include <map>
#include <stdexcept>
#include <iomanip>

struct sparse_matrix_2d
{
    struct index ;

    template< std::size_t ROWS, std::size_t COLS >
        explicit inline sparse_matrix_2d( double (&array)[ROWS][COLS] )
            : nrows(ROWS), ncols(COLS)
        {
            for( std::size_t i = 0 ; i < ROWS ; ++i )
                for( std::size_t j = 0 ; j < COLS ; ++j )
                    if( is_non_zero( array[i][j] ) )
                        non_zero_elements[index(i,j)] = array[i][j] ;
        }

    double at( std::size_t row, std::size_t col ) const
    {
        map_t::const_iterator iter = non_zero_elements.find( index(row,col) ) ;
        return iter != non_zero_elements.end() ? iter->second : 0 ;
    }

    void at( std::size_t row, std::size_t col, double value )
    {
        if( (row>=nrows) || (col>=ncols) ) throw std::out_of_range("invalid index") ;
        else if( is_non_zero(value) ) non_zero_elements[ index(row,col) ] = value ;
    }

    inline friend std::ostream& operator<< ( std::ostream& stm, const sparse_matrix_2d& m )
    {
            for( std::size_t i = 0 ; i < m.nrows ; ++i )
            {
                for( std::size_t j = 0 ; j < m.ncols ; ++j )
                      stm << std::setw(10) << m.at(i,j) ;
                stm << '\n' ;
            }
            return stm ;
    }

    struct index
    {
        index( std::size_t r, std::size_t c ) : row(r), col(c) {}
        const std::size_t row ;
        const std::size_t col ;
        bool operator< ( const index& that ) const
        {
          if( row < that.row ) return true ;
          else if( row == that.row ) return col < that.col ;
          else return false ;
        }
    };

    private:
        std::size_t nrows ;
        std::size_t ncols ;
        typedef std::map< index, double > map_t ;
        map_t non_zero_elements ;
        static const double EPSILON ;
        static inline bool is_non_zero( double d ) { return std::abs(d) > EPSILON ; }
};

const double sparse_matrix_2d::EPSILON = 1.e-8 ;

int main()
{
    double d[3][4] = { {1,0,0,2}, {0,0,3,0}, {4,0,5,0} } ;
    sparse_matrix_2d mtx(d) ;
    std::cout << std::fixed << std::setprecision(2) << std::showpoint << mtx << '\n' ;

}

Thank you very much for your feedback. it is very helpful to me.Could you please help me to solve this confusion. I defined the operator () to get the matrix value at (i, j),one for the get opearation, and one for the set opreation. if the value at(i, j) is zero, it can not be set again, it will give the error information about " can not set a new value do not in the stored row col. "
I create a matrix MA for the test, and the matrix at the place of (2,2) is zero, when i output its value with "cout<<MA(2,2)", it do not use the const version operator(), it use the second version instead, and give me the error information about " can not set a new value do not in the stored row col. "
Since Cout<<MA(2,2) is a right operation, why it do not use the const version?
Regards.

const double CoordMatrix::operator ( )(const int &aRow, const int &aCol) const
		 {  
		    assert(aRow<numRow && aCol<numRow);

			for (unsigned int i=0; i<row.size(); i++)
			  {
			   	  if(row[i]==aRow && col[i]==aCol ) return val[i];					        
				}
			return 0.0;
		  }


	  ////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////
	   // set the value of matrix at the places arow acol
 double& CoordMatrix::operator ( )(const int &aRow, const int &aCol) 
		 {  
		    assert(aRow<numRow && aCol<numRow);

			for (unsigned int i=0; i<row.size(); i++)
			  { 	  
			     if(row[i]==aRow && col[i]==aCol ) return val[i];	
				}
			cerr<<" can not set a new value do not in the stored row col"<<endl;
			exit(1);
		  }

You should consider using a std::map to store your matrix values, rather than the system of std::vector objects that you currently have. This will make the class smaller, clearer and also make the retrieval of values faster. To simplify vijayan121's excellently comprehensive example, so it's easier to see where the advantage is, a sparse matrix class implemented using maps can be made like this:

#include <map>

template< typename element_type >
class sparse_matrix_2d
{
public:
   typedef std::pair< size_t, size_t > index;
   typedef std::map< index, element_type > matrix_element_data;

   const element_type& getElement( const index& idx ) const
   {   return ( m_data.count( idx ) > 0 ? m_data[ idx ] : 0 );   }

   element_type& getEditableElement( const index& idx )
   {   return const_cast< element_type >( getElement( idx ));    }
   
private:
   matrix_element_data m_data;
   size_t m_rows;
   size_t m_cols;
};

Obviously, you need to add lots of other bits to this, like constructors and set methods etc. but you get the idea :o)

Edited 5 Years Ago by ravenous: n/a

> Since Cout<<MA(2,2) is a right operation, why it do not use the const version?

When a non-static member function is overloaded on the cv qualifiers, the function is being overloaded on the type of the this pointer (which is an implicit parameter to the function). CoordMatrix* gives an exact match for the this pointer for the non-const overload; the const overload requires a conversion from CoordMatrix* to const CoordMatrix* . Normal overload resolution applies and the overload resolves to the one with an exact match.

There is an elementary error in one of the functions I posted earlier; you should be able to figure out what is missing on your on.

// ...
    void at( std::size_t row, std::size_t col, double value )
    {
        if( (row>=nrows) || (col>=ncols) ) throw std::out_of_range("invalid index") ;
        else if( is_non_zero(value) ) non_zero_elements[ index(row,col) ] = value ;
        // TODO: what is being missed here?
    }
    // ...

> Since Cout<<MA(2,2) is a right operation, why it do not use the const version?

When a non-static member function is overloaded on the cv qualifiers, the function is being overloaded on the type of the this pointer (which is an implicit parameter to the function). CoordMatrix* gives an exact match for the this pointer for the non-const overload; the const overload requires a conversion from CoordMatrix* to const CoordMatrix* . Normal overload resolution applies and the overload resolves to the one with an exact match.

There is an elementary error in one of the functions I posted earlier; you should be able to figure out what is missing on your on.

// ...
    void at( std::size_t row, std::size_t col, double value )
    {
        if( (row>=nrows) || (col>=ncols) ) throw std::out_of_range("invalid index") ;
        else if( is_non_zero(value) ) non_zero_elements[ index(row,col) ] = value ;
        // TODO: what is being missed here?
    }
    // ...

Thank you very much for your help. I wrote a new class using the map structure yesterday afternoon following your thinking, it is very efficient and compact, i will use it for my CFD compulation in the future.
By the way, "non_zero_elements[ index(row,col) ]", the index(row, col) use the structure name index directly. why it works? I think we should creat a new instance like " index aIndex(row,col), and then "non_zero_elements[ AIndex(row,col) ]".

Regards.
Your Sincerely.

How to identify the column numbers in the i'th row, i tried to use find algorithm, but it can only return one values? Does STL has this kind of algorithm to perfrom this task.

> "non_zero_elements[ index(row,col) ]", the index(row, col) use the structure name index directly.
> why it works? I think we should creat a new instance like
> " index aIndex(row,col), and then "non_zero_elements[ AIndex(row,col) ]".

See: http://people.cs.vt.edu/~kafura/cs2704/anonymous.html


> How to identify the column numbers in the i'th row,
> i tried to use find algorithm, but it can only return one values?
> Does STL has this kind of algorithm to perfrom this task.

There is no standard algorithm which will do this. But it is quite easy to write a function that would do the job.

If the keys in the map are ordered by the predicate bool operator< ( const index& that ) const given earlier:
that is row by row in ascending order, and within each row, col by col in in ascending order, then:

// returns all column numbers in the i'th row
std::vector<std::size_t> sparse_matrix_2d::column_numbers( std::size_t i ) const
{
    typedef std::map< index, double >::const_iterator iterator ;
    iterator begin = non_zero_elements.lower_bound( index( i, 0 ) ) ;
    iterator end = non_zero_elements.upper_bound( index( i, ncols ) ) ;

    std::vector<std::size_t> col_numbers ;
    for( ; begin != end ; ++begin ) col_numbers.push_back( begin->first.col ) ;

    return col_numbers ;
}

> "non_zero_elements[ index(row,col) ]", the index(row, col) use the structure name index directly.
> why it works? I think we should creat a new instance like
> " index aIndex(row,col), and then "non_zero_elements[ AIndex(row,col) ]".

See: http://people.cs.vt.edu/~kafura/cs2704/anonymous.html


> How to identify the column numbers in the i'th row,
> i tried to use find algorithm, but it can only return one values?
> Does STL has this kind of algorithm to perfrom this task.

There is no standard algorithm which will do this. But it is quite easy to write a function that would do the job.

If the keys in the map are ordered by the predicate bool operator< ( const index& that ) const given earlier:
that is row by row in ascending order, and within each row, col by col in in ascending order, then:

// returns all column numbers in the i'th row
std::vector<std::size_t> sparse_matrix_2d::column_numbers( std::size_t i ) const
{
    typedef std::map< index, double >::const_iterator iterator ;
    iterator begin = non_zero_elements.lower_bound( index( i, 0 ) ) ;
    iterator end = non_zero_elements.upper_bound( index( i, ncols ) ) ;

    std::vector<std::size_t> col_numbers ;
    for( ; begin != end ; ++begin ) col_numbers.push_back( begin->first.col ) ;

    return col_numbers ;
}

Dear vijayan121:
Thank you veru much for your help. I use the column numbers for the sparse matrix and vector multiplication. This operation is repeated for many times in the solution process of partial differnential equation using conjugate gradient method.
Since the spase pattern is not changed in the solution process, it is obviously more effciently to store the column numbers of each row in the mesh generation process befor the solution begin in another "class Mesh" by using "verctor< vector<int> >" structure. But how transfer this row-column matrix in the " class Mesh" to the "class CoordMatrix".
Could you please give me some advices for this. I just move from fortran to C++ for several months. The information exchange between different classes in c++ object programming confused me many times. In the fortran , i use global variables which are can be used by all the other functions.But in the C++ object programming, it seems there are no concept about the global variables.
Regards
Your Sincerely

Another way using std::vector( not tested nor compiled )

template<typename T>
class Array2D{
public:
 typedef std::vector<T> ColumnVector;
 typedef std::vector<ColumnVector> Array2DType;
private:
 Array2DType _array;
public:
 Array2D(const int rowSize, const int colSize, const int initVal = 0)
 : _array( rowSize, ColumnVector(colSize, initVal ) );

 const T& operator()(const int row, const int col)const{
   return _array[row][col];
 }
 T& operator()(const int row, const int col){
   return _array[row][col];
 }
};

Another way is to simulate 2D array by using 1D.

> class Mesh" by using "verctor< vector<int> >" structure.
> But how transfer this row-column matrix in the " class Mesh" to the "class CoordMatrix".

This is one way of doing it.

1. Declare mesh to be a friend of the earlier sparse_matrix_2d .

2. In class mesh:
a. define a constructor that can initialize a mesh from a sparse_matrix_2d .
b. define a function that can populate a sparse_matrix_2d with the data in the mesh .

Code only for expositional purpose:

struct sparse_matrix_2d
{
    // other sparse_matrix_2d stuff
    
    friend struct mesh ;
};

struct mesh
{
    struct mesh_element
    {
        mesh_element( std::size_t c, double v ) : column(c), value(v) {}
        std::size_t column ;
        double value ;
    };

    std::size_t ncols ;
    std::vector< std::vector<mesh_element> > mesh_data ;

    explicit mesh( const sparse_matrix_2d& mtx ) : ncols(mtx.ncols), mesh_data(mtx.nrows)
    {
         const sparse_matrix_2d::map_t& non_zero_elements = mtx.non_zero_elements ;
         typedef sparse_matrix_2d::map_t::const_iterator iterator ;
         iterator end = non_zero_elements.end() ;
         for( iterator iter = non_zero_elements.begin() ; iter != end ; ++iter )
            mesh_data[ iter->first.row ].push_back(
                                        mesh_element( iter->first.col, iter->second ) ) ;

    }

    sparse_matrix_2d& populate( sparse_matrix_2d& mtx )
    {
        mtx.non_zero_elements.clear() ;
        mtx.nrows = mesh_data.size() ;
        mtx.ncols = ncols ;

        for( std::size_t r = 0 ; r < mesh_data.size() ; ++r )
            for( std::size_t c = 0 ; c < mesh_data[r].size() ; ++c )
               mtx.non_zero_elements[ sparse_matrix_2d::index(r,c) ] = mesh_data[r][c].value ;

        return mtx ;
    }

    // other mesh operations
};

> class Mesh" by using "verctor< vector<int> >" structure.
> But how transfer this row-column matrix in the " class Mesh" to the "class CoordMatrix".

This is one way of doing it.

1. Declare mesh to be a friend of the earlier sparse_matrix_2d .

2. In class mesh:
a. define a constructor that can initialize a mesh from a sparse_matrix_2d .
b. define a function that can populate a sparse_matrix_2d with the data in the mesh .

Code only for expositional purpose:

struct sparse_matrix_2d
{
    // other sparse_matrix_2d stuff
    
    friend struct mesh ;
};

struct mesh
{
    struct mesh_element
    {
        mesh_element( std::size_t c, double v ) : column(c), value(v) {}
        std::size_t column ;
        double value ;
    };

    std::size_t ncols ;
    std::vector< std::vector<mesh_element> > mesh_data ;

    explicit mesh( const sparse_matrix_2d& mtx ) : ncols(mtx.ncols), mesh_data(mtx.nrows)
    {
         const sparse_matrix_2d::map_t& non_zero_elements = mtx.non_zero_elements ;
         typedef sparse_matrix_2d::map_t::const_iterator iterator ;
         iterator end = non_zero_elements.end() ;
         for( iterator iter = non_zero_elements.begin() ; iter != end ; ++iter )
            mesh_data[ iter->first.row ].push_back(
                                        mesh_element( iter->first.col, iter->second ) ) ;

    }

    sparse_matrix_2d& populate( sparse_matrix_2d& mtx )
    {
        mtx.non_zero_elements.clear() ;
        mtx.nrows = mesh_data.size() ;
        mtx.ncols = ncols ;

        for( std::size_t r = 0 ; r < mesh_data.size() ; ++r )
            for( std::size_t c = 0 ; c < mesh_data[r].size() ; ++c )
               mtx.non_zero_elements[ sparse_matrix_2d::index(r,c) ] = mesh_data[r][c].value ;

        return mtx ;
    }

    // other mesh operations
};

Dear vijayan121:
I came with the idear of visit all the column index and matrix values with an iterator. such as
for(iter=begin(ith row); iter!=end(ith row); iter++)
cout<<*iter;
the overloaded * to give the matrix value. the ++ to move the pointer to every columnes in the ith row.
I write a new class, but it seeme the idear doesnot work. could you please give me some ideas.

//
//      CoordMatrixIterator.h   define the matrix stored in a coordinate format
//

#ifndef COORDMATRIXITERATOR_H
#define COORDMATRIXITERATOR_H
#include "assert.h"

#include<iomanip>
#include<iostream>
#include<map>
#include<algorithm>
#include <valarray>
#include"CoordMatrixIterator.h"
using namespace std;
class CoordMatrixIterator
			{
			private:
   			   CoordMatrixMap A;
			public:
				CoordMatrixIterator(const CoordMatrixMap &aA){A=aA;};
				~CoordMatrixIterator(){};
				
				const  CoordMatrixMap::Coord& begin(const int& aRow) const
				{
				    typedef map< CoordMatrixMap::Coord, double> ::const_iterator iter;
				    iter begin=A.val.lower_bound( CoordMatrixMap::Coord(aRow, 0) );
					return begin->first;
				}
				
				const CoordMatrixMap::Coord& end(const int& aRow) const
				{
				    typedef map< CoordMatrixMap::Coord, double> ::const_iterator iter;
		                    iter end=A.val.upper_bound( CoordMatrixMap::Coord(aRow, 5) );
					return end->first;
				}
                
				CoordMatrixMap::Coord operator ++()const
					{    }

   

			  };

#endif

> I came with the idear of visit all the column index and matrix values with an iterator. such as
> for(iter=begin(ith row); iter!=end(ith row); iter++) cout<<*iter;
> the overloaded * to give the matrix value.

As long as the keys in the map are ordered row by row in ascending order, and within each row, col by col in ascending order, you can base your iterator on the map's iterator. Just use std::map<>::lower_bound() and std::map<>::upper_bound() to get begin and end for a particular row.

struct sparse_matrix_2d
{
    // ...

    struct index
    {
        index( std::size_t r, std::size_t c ) : row(r), col(c) {}
        const std::size_t row ;
        const std::size_t col ;
        bool operator< ( const index& that ) const
        {
          if( row < that.row ) return true ;
          else if( row == that.row ) return col < that.col ;
          else return false ;
        }
    };

    typedef std::map< index, double > map_t ;
    struct row_iterator : std::iterator< std::forward_iterator_tag, double >
    {
        explicit row_iterator( map_t::iterator i ) : iter(i) {}

        double& operator* () { return iter->second ; }

        row_iterator& operator++() { ++iter ; return *this ; }
        row_iterator operator++(int)
        { row_iterator before_increment(*this) ; ++iter ; return before_increment ; }

        bool operator== ( const row_iterator& that ) const
        { return iter == that.iter ; }
        bool operator!= ( const row_iterator& that ) const
        { return iter != that.iter ; }

        private:
            map_t::iterator iter ;
    };

    row_iterator row_begin( std::size_t row_num )
    { return row_iterator( non_zero_elements.lower_bound( index(row_num,0) ) ) ; }

    row_iterator row_end( std::size_t row_num )
    { return row_iterator( non_zero_elements.upper_bound( index(row_num,ncols) ) ) ; }

    private:
        std::size_t nrows ;
        std::size_t ncols ;
        map_t non_zero_elements ;

        // ...
};
This article has been dead for over six months. Start a new discussion instead.