#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;
}
I2luY2x1ZGUgPHZlY3Rvcj4KI2luY2x1ZGUgPG1hcD4KI2luY2x1ZGUgPGlvc3RyZWFtPgoKdXNpbmcgbmFtZXNwYWNlIHN0ZDsKCnRlbXBsYXRlIDx0eXBlbmFtZSBUPgpjbGFzcyBTcGFyc2VNYXRyaXgKewpwcml2YXRlOgogICAgU3BhcnNlTWF0cml4KCk7CiAgICAKcHVibGljOgogICAgU3BhcnNlTWF0cml4KGludCByb3csIGludCBjb2wpOwogICAgCiAgICBUIEdldChpbnQgcm93LCBpbnQgY29sKSBjb25zdDsKICAgIHZvaWQgUHV0KGludCByb3csIGludCBjb2wsIFQgdmFsdWUpOwogICAgCiAgICBpbnQgR2V0Um93Q291bnQoKSBjb25zdDsKICAgIGludCBHZXRDb2xDb3VudCgpIGNvbnN0OwogICAgCiAgICBzdGF0aWMgdm9pZCBHYXVzc1NlaWRlbExpbmVhclNvbHZlKGNvbnN0IFNwYXJzZU1hdHJpeDxUPiYgQSwgY29uc3QgU3BhcnNlTWF0cml4PFQ+JiBCLCBTcGFyc2VNYXRyaXg8VD4mIFgpOwogICAgCnByaXZhdGU6CiAgICBpbnQgZGltX3JvdzsKICAgIGludCBkaW1fY29sOwogICAgCiAgICB2ZWN0b3I8bWFwPGludCwgVD4gPiB2YWx1ZXNfYnlfcm93OwogICAgdmVjdG9yPG1hcDxpbnQsIFQ+ID4gdmFsdWVzX2J5X2NvbDsKfTsKCnRlbXBsYXRlIDx0eXBlbmFtZSBUPiBTcGFyc2VNYXRyaXg8VD46OlNwYXJzZU1hdHJpeChpbnQgciwgaW50IGMpIDogdmFsdWVzX2J5X3JvdyhyKSwgdmFsdWVzX2J5X2NvbChjKQp7CiAgICBkaW1fcm93ID0gcjsKICAgIGRpbV9jb2wgPSBjOwp9Cgp0ZW1wbGF0ZSA8dHlwZW5hbWUgVD4gaW50IFNwYXJzZU1hdHJpeDxUPjo6R2V0Um93Q291bnQoKSBjb25zdAp7CiAgICByZXR1cm4gZGltX3JvdzsKfQoKdGVtcGxhdGUgPHR5cGVuYW1lIFQ+IGludCBTcGFyc2VNYXRyaXg8VD46OkdldENvbENvdW50KCkgY29uc3QKewogICAgcmV0dXJuIGRpbV9jb2w7Cn0KCnRlbXBsYXRlIDx0eXBlbmFtZSBUPiBUIFNwYXJzZU1hdHJpeDxUPjo6R2V0KGludCByLCBpbnQgYykgY29uc3QKewogICAgdHlwZW5hbWUgbWFwPGludCwgVD46OmNvbnN0X2l0ZXJhdG9yIGkgPSB2YWx1ZXNfYnlfcm93W3JdLmZpbmQoYyk7CiAgICByZXR1cm4gKGkgIT0gdmFsdWVzX2J5X3Jvd1tyXS5lbmQoKSkgPyBpLT5zZWNvbmQgOiAoVCkgMDsKfQoKdGVtcGxhdGUgPHR5cGVuYW1lIFQ+IHZvaWQgU3BhcnNlTWF0cml4PFQ+OjpQdXQoaW50IHIsIGludCBjLCBUIHZhbHVlKQp7CiAgICBpZiAodmFsdWUgIT0gKFQpIDApCiAgICB7CiAgICAgICAgdmFsdWVzX2J5X3Jvd1tyXVtjXSA9IHZhbHVlOwogICAgICAgIHZhbHVlc19ieV9jb2xbY11bcl0gPSB2YWx1ZTsKICAgIH0KICAgIGVsc2UKICAgIHsKICAgICAgICB2YWx1ZXNfYnlfcm93W3JdLmVyYXNlKGMpOwogICAgICAgIHZhbHVlc19ieV9jb2xbY10uZXJhc2Uocik7CiAgICB9Cn0KCnRlbXBsYXRlIDx0eXBlbmFtZSBUPiB2b2lkIFNwYXJzZU1hdHJpeDxUPjo6R2F1c3NTZWlkZWxMaW5lYXJTb2x2ZShjb25zdCBTcGFyc2VNYXRyaXg8VD4mIEEsIGNvbnN0IFNwYXJzZU1hdHJpeDxUPiYgQiwgU3BhcnNlTWF0cml4PFQ+JiBYKQp7CiAgICBpZiAoQS5HZXRSb3dDb3VudCgpICE9IEEuR2V0Q29sQ291bnQoKSB8fCBBLkdldENvbENvdW50KCkgIT0gQi5HZXRSb3dDb3VudCgpIHx8IEEuR2V0Q29sQ291bnQoKSAhPSBYLkdldFJvd0NvdW50KCkgfHwKICAgICAgICBCLkdldENvbENvdW50KCkgIT0gMSB8fCBYLkdldENvbENvdW50KCkgIT0gMSkgcmV0dXJuOwogICAgICAgIAogICAgZm9yIChpbnQgeiA9IDA7IHogPCAxMDA7ICsreikKICAgIHsKICAgICAgICBUIGFnZ3JlZ2F0ZTsKICAgICAgICBmb3IgKGludCByID0gMDsgciA8IFguR2V0Um93Q291bnQoKTsgKytyKQogICAgICAgIHsKICAgICAgICAgICAgYWdncmVnYXRlID0gKFQpIDA7CiAgICAgICAgICAgIGZvciAodHlwZW5hbWUgbWFwPGludCwgVD46OmNvbnN0X2l0ZXJhdG9yIGMgPSBBLnZhbHVlc19ieV9yb3dbcl0uYmVnaW4oKTsgYyAhPSBBLnZhbHVlc19ieV9yb3dbcl0uZW5kKCk7ICsrYykKICAgICAgICAgICAgewogICAgICAgICAgICAgICAgaWYgKGMtPmZpcnN0ID09IHIpIGNvbnRpbnVlOwogICAgICAgICAgICAgICAgYWdncmVnYXRlICs9IGMtPnNlY29uZCAqIFguR2V0KGMtPmZpcnN0LCAwKTsKICAgICAgICAgICAgfQogICAgICAgICAgICBYLlB1dChyLCAwLCAoQi5HZXQociwgMCkgLSBhZ2dyZWdhdGUpIC8gQS5HZXQociwgcikpOwogICAgICAgIH0KICAgIH0KfQoKaW50IG1haW4oKQp7CiAgICBpbnQgbWFnbml0dWRlID0gMTAwMDA7CiAgICAKICAgIGNvdXQgPDwgIkNyZWF0aW5nIEEuLi4iIDw8IGVuZGw7CiAgICBTcGFyc2VNYXRyaXg8ZG91YmxlPiBBKG1hZ25pdHVkZSwgbWFnbml0dWRlKTsKICAgIGZvciAoaW50IGkgPSAwOyBpIDwgbWFnbml0dWRlOyArK2kpIEEuUHV0KGksIGksIDQpOwogICAgY291dCA8PCAic3RhZ2UgMSBjb21wbGV0ZS4iIDw8IGVuZGw7CiAgICBmb3IgKGludCBpID0gMDsgaSA8IG1hZ25pdHVkZSAtIDE7ICsraSkKICAgIHsKICAgICAgICBBLlB1dChpLCBpICsgMSwgLTEpOwogICAgICAgIEEuUHV0KGkgKyAxLCBpLCAtMSk7CiAgICB9IC8vIEEgaXMgbm93IGEgMTAwMDB4MTAwMDAgbWF0cml4IHNpbWlsYXIgdG8gdGhlIG9uZSB5b3UgdXNlIGFzIGFuIGV4YW1wbGUgaW4geW91ciBxdWVzdGlvbgogICAgY291dCA8PCAic3RhdGUgMiBjb21wbGV0ZS4iIDw8IGVuZGw7CiAgICAKICAgIGNvdXQgPDwgIkNyZWF0aW5nIEIuLi4iIDw8IGVuZGw7CiAgICBTcGFyc2VNYXRyaXg8ZG91YmxlPiBCKG1hZ25pdHVkZSwgMSk7CiAgICBmb3IgKGludCBpID0gMDsgaSA8IG1hZ25pdHVkZTsgKytpKSBCLlB1dChpLCAwLCAxLjUgKiAoaSAtIG1hZ25pdHVkZSAvIDIpICsgLjc1KTsKICAgIAogICAgY291dCA8PCAiQ3JlYXRpbmcgWC4uLiIgPDwgZW5kbDsKICAgIFNwYXJzZU1hdHJpeDxkb3VibGU+IFgobWFnbml0dWRlLCAxKTsKICAgIAogICAgY291dCA8PCAiTGluZWFyIFNvbHZpbmcuLi4iIDw8IGVuZGw7CiAgICBmb3IgKGludCBpID0gMDsgaSA8IDUwOyArK2kpIGNvdXQgPDwgaSA8PCAiOiAiIDw8IFguR2V0KGksIDApIDw8IGVuZGw7CiAgICBTcGFyc2VNYXRyaXg8ZG91YmxlPjo6R2F1c3NTZWlkZWxMaW5lYXJTb2x2ZShBLCBCLCBYKTsKICAgIGZvciAoaW50IGkgPSAwOyBpIDwgNTA7ICsraSkgY291dCA8PCBpIDw8ICI6ICIgPDwgWC5HZXQoaSwgMCkgPDwgZW5kbDsKICAgIFNwYXJzZU1hdHJpeDxkb3VibGU+OjpHYXVzc1NlaWRlbExpbmVhclNvbHZlKEEsIEIsIFgpOwogICAgZm9yIChpbnQgaSA9IDA7IGkgPCA1MDsgKytpKSBjb3V0IDw8IGkgPDwgIjogIiA8PCBYLkdldChpLCAwKSA8PCBlbmRsOwogICAgCiAgICByZXR1cm4gMDsKfQ==