fork download
  1. #include <vector>
  2. #include <map>
  3. #include <iostream>
  4.  
  5. using namespace std;
  6.  
  7. template <typename T>
  8. class SparseMatrix
  9. {
  10. private:
  11. SparseMatrix();
  12.  
  13. public:
  14. SparseMatrix(int row, int col);
  15.  
  16. T Get(int row, int col) const;
  17. void Put(int row, int col, T value);
  18.  
  19. int GetRowCount() const;
  20. int GetColCount() const;
  21.  
  22. static void GaussSeidelLinearSolve(const SparseMatrix<T>& A, const SparseMatrix<T>& B, SparseMatrix<T>& X);
  23.  
  24. private:
  25. int dim_row;
  26. int dim_col;
  27.  
  28. vector<map<int, T> > values_by_row;
  29. vector<map<int, T> > values_by_col;
  30. };
  31.  
  32. template <typename T> SparseMatrix<T>::SparseMatrix(int r, int c) : values_by_row(r), values_by_col(c)
  33. {
  34. dim_row = r;
  35. dim_col = c;
  36. }
  37.  
  38. template <typename T> int SparseMatrix<T>::GetRowCount() const
  39. {
  40. return dim_row;
  41. }
  42.  
  43. template <typename T> int SparseMatrix<T>::GetColCount() const
  44. {
  45. return dim_col;
  46. }
  47.  
  48. template <typename T> T SparseMatrix<T>::Get(int r, int c) const
  49. {
  50. typename map<int, T>::const_iterator i = values_by_row[r].find(c);
  51. return (i != values_by_row[r].end()) ? i->second : (T) 0;
  52. }
  53.  
  54. template <typename T> void SparseMatrix<T>::Put(int r, int c, T value)
  55. {
  56. if (value != (T) 0)
  57. {
  58. values_by_row[r][c] = value;
  59. values_by_col[c][r] = value;
  60. }
  61. else
  62. {
  63. values_by_row[r].erase(c);
  64. values_by_col[c].erase(r);
  65. }
  66. }
  67.  
  68. template <typename T> void SparseMatrix<T>::GaussSeidelLinearSolve(const SparseMatrix<T>& A, const SparseMatrix<T>& B, SparseMatrix<T>& X)
  69. {
  70. if (A.GetRowCount() != A.GetColCount() || A.GetColCount() != B.GetRowCount() || A.GetColCount() != X.GetRowCount() ||
  71. B.GetColCount() != 1 || X.GetColCount() != 1) return;
  72.  
  73. for (int z = 0; z < 100; ++z)
  74. {
  75. T aggregate;
  76. for (int r = 0; r < X.GetRowCount(); ++r)
  77. {
  78. aggregate = (T) 0;
  79. for (typename map<int, T>::const_iterator c = A.values_by_row[r].begin(); c != A.values_by_row[r].end(); ++c)
  80. {
  81. if (c->first == r) continue;
  82. aggregate += c->second * X.Get(c->first, 0);
  83. }
  84. X.Put(r, 0, (B.Get(r, 0) - aggregate) / A.Get(r, r));
  85. }
  86. }
  87. }
  88.  
  89. int main()
  90. {
  91. int magnitude = 10000;
  92.  
  93. cout << "Creating A..." << endl;
  94. SparseMatrix<double> A(magnitude, magnitude);
  95. for (int i = 0; i < magnitude; ++i) A.Put(i, i, 4);
  96. cout << "stage 1 complete." << endl;
  97. for (int i = 0; i < magnitude - 1; ++i)
  98. {
  99. A.Put(i, i + 1, -1);
  100. A.Put(i + 1, i, -1);
  101. } // A is now a 10000x10000 matrix similar to the one you use as an example in your question
  102. cout << "state 2 complete." << endl;
  103.  
  104. cout << "Creating B..." << endl;
  105. SparseMatrix<double> B(magnitude, 1);
  106. for (int i = 0; i < magnitude; ++i) B.Put(i, 0, 1.5 * (i - magnitude / 2) + .75);
  107.  
  108. cout << "Creating X..." << endl;
  109. SparseMatrix<double> X(magnitude, 1);
  110.  
  111. cout << "Linear Solving..." << endl;
  112. for (int i = 0; i < 50; ++i) cout << i << ": " << X.Get(i, 0) << endl;
  113. SparseMatrix<double>::GaussSeidelLinearSolve(A, B, X);
  114. for (int i = 0; i < 50; ++i) cout << i << ": " << X.Get(i, 0) << endl;
  115. SparseMatrix<double>::GaussSeidelLinearSolve(A, B, X);
  116. for (int i = 0; i < 50; ++i) cout << i << ": " << X.Get(i, 0) << endl;
  117.  
  118. return 0;
  119. }
Success #stdin #stdout 0.52s 6028KB
stdin
Standard input is empty
stdout
Creating A...
stage 1 complete.
state 2 complete.
Creating B...
Creating X...
Linear Solving...
0: 0
1: 0
2: 0
3: 0
4: 0
5: 0
6: 0
7: 0
8: 0
9: 0
10: 0
11: 0
12: 0
13: 0
14: 0
15: 0
16: 0
17: 0
18: 0
19: 0
20: 0
21: 0
22: 0
23: 0
24: 0
25: 0
26: 0
27: 0
28: 0
29: 0
30: 0
31: 0
32: 0
33: 0
34: 0
35: 0
36: 0
37: 0
38: 0
39: 0
40: 0
41: 0
42: 0
43: 0
44: 0
45: 0
46: 0
47: 0
48: 0
49: 0
0: -2744.72
1: -3479.61
2: -3675.98
3: -3728.04
4: -3741.44
5: -3744.49
6: -3744.75
7: -3744.28
8: -3743.6
9: -3742.87
10: -3742.12
11: -3741.37
12: -3740.62
13: -3739.87
14: -3739.12
15: -3738.37
16: -3737.62
17: -3736.87
18: -3736.12
19: -3735.37
20: -3734.62
21: -3733.87
22: -3733.12
23: -3732.37
24: -3731.62
25: -3730.87
26: -3730.12
27: -3729.37
28: -3728.62
29: -3727.88
30: -3727.12
31: -3726.38
32: -3725.62
33: -3724.88
34: -3724.12
35: -3723.38
36: -3722.62
37: -3721.88
38: -3721.12
39: -3720.38
40: -3719.62
41: -3718.88
42: -3718.12
43: -3717.38
44: -3716.62
45: -3715.88
46: -3715.12
47: -3714.38
48: -3713.62
49: -3712.88
0: -2744.72
1: -3479.61
2: -3675.98
3: -3728.04
4: -3741.44
5: -3744.49
6: -3744.75
7: -3744.28
8: -3743.6
9: -3742.87
10: -3742.12
11: -3741.37
12: -3740.62
13: -3739.87
14: -3739.12
15: -3738.37
16: -3737.62
17: -3736.87
18: -3736.12
19: -3735.37
20: -3734.62
21: -3733.87
22: -3733.12
23: -3732.37
24: -3731.62
25: -3730.87
26: -3730.12
27: -3729.37
28: -3728.62
29: -3727.88
30: -3727.12
31: -3726.38
32: -3725.62
33: -3724.88
34: -3724.12
35: -3723.38
36: -3722.62
37: -3721.88
38: -3721.12
39: -3720.38
40: -3719.62
41: -3718.88
42: -3718.12
43: -3717.38
44: -3716.62
45: -3715.88
46: -3715.12
47: -3714.38
48: -3713.62
49: -3712.88