#include <vector>
#include <map>
#include <iostream>

using namespace std;

template <typename T>
class SparseMatrix
{
private:
    SparseMatrix();
    
public:
    SparseMatrix(int row, int col);
    
    T Get(int row, int col) const;
    void Put(int row, int col, T value);
    
    int GetRowCount() const;
    int GetColCount() const;
    
    static void GaussSeidelLinearSolve(const SparseMatrix<T>& A, const SparseMatrix<T>& B, SparseMatrix<T>& X);
    
private:
    int dim_row;
    int dim_col;
    
    vector<map<int, T> > values_by_row;
    vector<map<int, T> > values_by_col;
};

template <typename T> SparseMatrix<T>::SparseMatrix(int r, int c) : values_by_row(r), values_by_col(c)
{
    dim_row = r;
    dim_col = c;
}

template <typename T> int SparseMatrix<T>::GetRowCount() const
{
    return dim_row;
}

template <typename T> int SparseMatrix<T>::GetColCount() const
{
    return dim_col;
}

template <typename T> T SparseMatrix<T>::Get(int r, int c) const
{
    typename map<int, T>::const_iterator i = values_by_row[r].find(c);
    return (i != values_by_row[r].end()) ? i->second : (T) 0;
}

template <typename T> void SparseMatrix<T>::Put(int r, int c, T value)
{
    if (value != (T) 0)
    {
        values_by_row[r][c] = value;
        values_by_col[c][r] = value;
    }
    else
    {
        values_by_row[r].erase(c);
        values_by_col[c].erase(r);
    }
}

template <typename T> void SparseMatrix<T>::GaussSeidelLinearSolve(const SparseMatrix<T>& A, const SparseMatrix<T>& B, SparseMatrix<T>& X)
{
    if (A.GetRowCount() != A.GetColCount() || A.GetColCount() != B.GetRowCount() || A.GetColCount() != X.GetRowCount() ||
        B.GetColCount() != 1 || X.GetColCount() != 1) return;
        
    for (int z = 0; z < 100; ++z)
    {
        T aggregate;
        for (int r = 0; r < X.GetRowCount(); ++r)
        {
            aggregate = (T) 0;
            for (typename map<int, T>::const_iterator c = A.values_by_row[r].begin(); c != A.values_by_row[r].end(); ++c)
            {
                if (c->first == r) continue;
                aggregate += c->second * X.Get(c->first, 0);
            }
            X.Put(r, 0, (B.Get(r, 0) - aggregate) / A.Get(r, r));
        }
    }
}

int main()
{
    int magnitude = 10000;
    
    cout << "Creating A..." << endl;
    SparseMatrix<double> A(magnitude, magnitude);
    for (int i = 0; i < magnitude; ++i) A.Put(i, i, 4);
    cout << "stage 1 complete." << endl;
    for (int i = 0; i < magnitude - 1; ++i)
    {
        A.Put(i, i + 1, -1);
        A.Put(i + 1, i, -1);
    } // A is now a 10000x10000 matrix similar to the one you use as an example in your question
    cout << "state 2 complete." << endl;
    
    cout << "Creating B..." << endl;
    SparseMatrix<double> B(magnitude, 1);
    for (int i = 0; i < magnitude; ++i) B.Put(i, 0, 1.5 * (i - magnitude / 2) + .75);
    
    cout << "Creating X..." << endl;
    SparseMatrix<double> X(magnitude, 1);
    
    cout << "Linear Solving..." << endl;
    for (int i = 0; i < 50; ++i) cout << i << ": " << X.Get(i, 0) << endl;
    SparseMatrix<double>::GaussSeidelLinearSolve(A, B, X);
    for (int i = 0; i < 50; ++i) cout << i << ": " << X.Get(i, 0) << endl;
    SparseMatrix<double>::GaussSeidelLinearSolve(A, B, X);
    for (int i = 0; i < 50; ++i) cout << i << ": " << X.Get(i, 0) << endl;
    
    return 0;
}