Sparse Matrix Implementation

[turbozedd]

The following is my implementation of a sparse matrix class. Feel free to comment on my bad coding style as well as use the code. The input file is the bcsstk14 matrix from http://math.nist.gov/MatrixMarket/

Space savings is huge for this matrix. It has about 40,000 nonzeros, the number of elements in the matrix is about 2000^2=4,000,000. So the sparse matrix takes up about 1% (plus pointers and stuff) of the space required for the dense matrix (array[][]).

There is a multiplication timer in the test file, it comes out with .005 seconds whereas dense multiplication takes about 7 seconds (Linux P4 2.3GHz) for this matrix. Which is a speedup factor of about 1400.

Matrix class header

#ifndef _MATRIX_H
#define	_MATRIX_H

#include <cstdlib>
#include <map>
#include <vector>

template <class T>
class Matrix
{
public:
    typedef std::map<size_t, std::map<size_t , T> > mat_t;
    typedef typename mat_t::iterator row_iter;
    typedef std::map<size_t, T> col_t;
    typedef typename col_t::iterator col_iter;

    Matrix(size_t i){ m=i; n=i; }
    Matrix(size_t i, size_t j){ m=i; n=j; }

    inline
    T& operator()(size_t i, size_t j)
    {
        if(i>=m || j>=n) throw;
        return mat[i][j];
    }
    inline
    T operator()(size_t i, size_t j) const
    {
        if(i>=m || j>=n) throw;
        return mat[i][j];
    }

    std::vector<T> operator*(const std::vector<T>& x)
    {  //Computes y=A*x
        if(this->m != x.size()) throw;

        std::vector<T> y(this->m);
        T sum;

        row_iter ii;
        col_iter jj;

        for(ii=this->mat.begin(); ii!=this->mat.end(); ii++){
            sum=0;
            for(jj=(*ii).second.begin(); jj!=(*ii).second.end(); jj++){
                sum += (*jj).second * x[(*jj).first];
            }
            y[(*ii).first]=sum;
        }

        return y;
    }

    void printMat()
    {
        row_iter ii;
        col_iter jj;
        for(ii=this->mat.begin(); ii!=this->mat.end(); ii++){
            for( jj=(*ii).second.begin(); jj!=(*ii).second.end(); jj++){
                std::cout << (*ii).first << ' ';
                std::cout << (*jj).first << ' ';
                std::cout << (*jj).second << std::endl;
            }
        } std::cout << std::endl;
    }
    
protected:
    Matrix(){}

private:
    mat_t mat;
    size_t m;
    size_t n;
};

#endif	/* _MATRIX_H */

test program

#include <cstdlib>
#include <iostream>
#include <fstream>
#include <cmath>
#include <time.h>
#include "Matrix.h"

using namespace std;

int init(Matrix<long double> &A)
{
    ifstream fin("matin.txt");
    int n,length,i,j;
    if(fin.is_open())
    {
        fin >> n;
	fin >> length;
        int k=0;
        while(k<length){ fin >> i >> j; fin >> A(i-1,j-1); k++; }
        fin.close();
    }
    else return 0;
    return n;
}

int main()
{
    Matrix<long double> A(1806);
    vector<long double> x(1806,1), y;
    long double sum;
    int n,k=0;
    clock_t t0,tf;
    cout.precision(17);

    n=init(A);
    
    //timed multiplication
    t0=clock();
    for(int i=0;i<1000;i++){ y=A*x; }
    tf=clock()-t0;
    
    //output last element
    cout << y[1805]
    //output time to multiply
    cout << "\n\n" << (double)tf/((double)(CLOCKS_PER_SEC)*(1000.0));

    //do naive multiply for last row to check last element of y
    long double tmp=0;
    for(int i=0;i<1806;i++) tmp+=A(1805,i);
    cout << "\n" << tmp <<"\n";

    return 0;
}

[jsmith]

Your #include guard does nothing. The #endif needs to be at the end of the file.

You should throw something; a generic throw can only be caught through ellipsis (which catches everything).

operator* should take its parameter by const reference.

You should write an ostream operator for output instead of printMat().

Why not write accessors on Matrix that allow the user to figure out the dimensions of the matrix?

[turbozedd]

Thanks for the comments jsmith.

I didn't really know what the #ifndef was for, my IDE generated it.

I need to read more about throwing...

I had operator* with const parameter, but I guess it disappeared during debug.

I had planned on doing ostream operator, but I haven't gotten to it yet.

I forgot about getdim().

I'll add changes to my original post.

[helios]

I'd like to suggest that the matrix be a std::map<std::pair<size_t,size_t>,T>

[turbozedd]

Then I have to store a column index up to n times for each element in a row. My implementation is basically Compressed Row Storage see http://www.netlib.org/linalg/html_templates/node90.html
This saves more space by adding a few more pointers.

This is almost the same as what MatLab uses.

[helios]

Are you sure the extra overhead of more trees doesn't offset whatever saving you get from that?

[turbozedd]

No... but how would you execute a fast MatVec multiply with std::map<std::pair<size_t,size_t>,T>? I need to have pointers to each row for mine to work.

And how can I determine the amount of memory being used, during execution that is. I tried doing some sizof() things but it only gives the size of the pointer. Is there a recursive sizeof() ? I do not really want to do a memory analysis by hand.

[helios]

Hmm... Well, if you need to rapidly access the elements of a single column/row, then your approach is indeed faster.

You can use a system monitor to check how much memory the program is using. For Windows, the best one I know of is Process Explorer.

[turbozedd]

I was hoping for something a bit more specific (yes indeed that was a pun). But I suppose if I have a big enough matrix it will dominate the memory slot.