My Project
Mat.h
Go to the documentation of this file.
1 #ifndef _MAT_H
2 #define _MAT_H
3 
4 #include <vector>
5 #include <stdexcept>
6 #include "Vec.h"
7 
9 class Mat {
10 public:
11  Mat(size_t nu) : nu{nu} {}
12  virtual double& operator()(size_t i, size_t j) = 0;
13  virtual double operator()(size_t i, size_t j) const = 0;
14  virtual double dot(const std::vector<double> &x, size_t r) const = 0;
15  virtual size_t memsize() = 0;
16  const size_t nu;
17 };
18 
19 
20 /* dense matrix, only good for small matrixes!*/
21 class DenseMat : public Mat {
22 public:
23 
24  DenseMat(size_t nu) : Mat{nu} {
25  // memory for coefficients
26  a = new double*[nu];
27  for (size_t i=0;i<nu;i++)
28  a[i] = new double[nu];
29 
30  // initialize all to zero
31  for (size_t i=0;i<nu;i++)
32  for (size_t j=0;j<nu;j++)
33  a[i][j] = 0;
34  }
35 
37  for (size_t i=0;i<nu;i++)
38  delete[] a[i];
39  delete[] a;
40  }
41 
42  // returns reference to the A[i,j] coefficient in the full matrix
43  double& operator()(size_t i, size_t j) {return a[i][j];}
44  double operator()(size_t i, size_t j) const {return a[i][j];}
45 
46  // returns dot product of matrix row r with vector x
47  double dot(const std::vector<double> &x, size_t r) const {
48  double sum = 0;
49  for (size_t c=0;c<nu;c++)
50  sum+=a[r][c]*x[c];
51  return sum;
52  }
53 
54  // return memory size in bytes
55  size_t memsize() {return nu*nu*sizeof(double);}
56 
57 protected:
58  double **a; //coefficients, [nu][nu]
59 };
60 
61 
62 /* sparse matrix with up to 7 non-zero values per row*/
63 class SparseMat : public Mat {
64 public:
65 
66  SparseMat(size_t nu) : Mat{nu} {
67 
68  // memory for coefficients
69  a = new double*[nu];
70  for (size_t i=0;i<nu;i++)
71  a[i] = new double[max_vals];
72 
73  // memory for column indexes
74  c = new int*[nu];
75  for (size_t i=0;i<nu;i++)
76  c[i] = new int[max_vals];
77 
78  // clear data, column set to -1 are not set
79  for (size_t i=0;i<nu;i++)
80  for (size_t j=0;j<max_vals;j++) {
81  c[i][j] = -1;
82  a[i][j] = 0;
83  }
84  }
85 
87  for (size_t i=0;i<nu;i++)
88  delete[] a[i];
89  delete[] a;
90 
91  for (size_t i=0;i<nu;i++)
92  delete[] c[i];
93  delete[] c;
94 
95  }
96 
97  // returns reference to the A[i,j] coefficient in the full matrix
98  double& operator()(size_t i, size_t j) {
99  //search for the sparse column corresponding to full matrix column j
100  for (int v=0;v<max_vals;v++) {
101  //did we reach an empty slot? If so, make it correspond to column j
102  if (c[i][v]<0) c[i][v] = j;
103 
104  // does this sparse column map to j? If so, return the coeff
105  if (c[i][v]==j) return a[i][v];
106  }
107 
108  //getting here implies that all max_val slots are already occupied by columns other than j
109  std::runtime_error("Sparse matrix too small!");
110  }
111 
112  // identical to the function above but for read-only access
113  double operator()(size_t i, size_t j) const {
114  //search for the sparse column corresponding to full matrix column j
115  for (int v=0;v<max_vals;v++) {
116  //did we reach an empty slot? If so, make it correspond to column j
117  if (c[i][v]<0) c[i][v] = j;
118 
119  // does this sparse column map to j? If so, return the coeff
120  if (c[i][v]==j) return a[i][v];
121  }
122 
123  //getting here implies that all max_val slots are already occupied by columns other than j
124  std::runtime_error("Sparse matrix too small!");
125  }
126 
127  // returns dot product of matrix row r with vector x
128  double dot(const std::vector<double> &x, size_t r) const {
129  double sum = 0;
130  // loop up to max_vals time and until c[r][v] becomes negative
131  for (size_t v=0;v<max_vals && c[r][v]>=0;v++)
132  sum+=a[r][v]*x[c[r][v]]; // c[v] is effective the "j" in full matrix
133  return sum;
134  }
135 
136  // return memory size in bytes
137  size_t memsize() {return nu*max_vals*sizeof(double);}
138 
139 protected:
140  static constexpr size_t max_vals = 7;
141  double **a;
142  int **c;
143 };
144 
145 
146 
147 #endif
virtual double dot(const std::vector< double > &x, size_t r) const =0
virtual double & operator()(size_t i, size_t j)=0
int ** c
columns in full matrix, [nu][7], -1 if not set
Definition: Mat.h:142
double dot(const std::vector< double > &x, size_t r) const
Definition: Mat.h:128
double & operator()(size_t i, size_t j)
Definition: Mat.h:43
~DenseMat()
Definition: Mat.h:36
DenseMat(size_t nu)
Definition: Mat.h:24
Mat(size_t nu)
Definition: Mat.h:11
size_t memsize()
matrix size in bytes
Definition: Mat.h:137
size_t memsize()
matrix size in bytes
Definition: Mat.h:55
double operator()(size_t i, size_t j) const
Definition: Mat.h:44
double ** a
Definition: Mat.h:58
double dot(const std::vector< double > &x, size_t r) const
Definition: Mat.h:47
double operator()(size_t i, size_t j) const
Definition: Mat.h:113
Definition: Mat.h:63
double & operator()(size_t i, size_t j)
Definition: Mat.h:98
~SparseMat()
Definition: Mat.h:86
matrix base class, specifies the interface to be implemented
Definition: Mat.h:9
Definition: Mat.h:21
SparseMat(size_t nu)
Definition: Mat.h:66
const size_t nu
number of rows (unknowns)
Definition: Mat.h:16
double ** a
coefficients, [nu][7]
Definition: Mat.h:141
virtual size_t memsize()=0
matrix size in bytes