OpenMEEG
fast_sparse_matrix.h
Go to the documentation of this file.
1 // Project Name: OpenMEEG (http://openmeeg.github.io)
2 // © INRIA and ENPC under the French open source license CeCILL-B.
3 // See full copyright notice in the file LICENSE.txt
4 // If you make a copy of this file, you must either:
5 // - provide also LICENSE.txt and modify this header to refer to it.
6 // - replace this header by the LICENSE.txt content.
7 
8 #pragma once
9 
10 #include "OpenMEEGMathsConfig.h"
11 #include "vector.h"
12 #include "sparse_matrix.h"
13 
14 namespace OpenMEEG {
15 
16  class OPENMEEGMATHS_EXPORT FastSparseMatrix
17  {
18  public:
19 
20  inline friend std::ostream& operator<<(std::ostream& f,const FastSparseMatrix &M);
21 
22  protected:
23 
24  double *tank;
25  size_t *js;
26  size_t *rowindex;
27  size_t m_nlin;
28  size_t m_ncol;
29 
30  inline void alloc(size_t nl, size_t nc, size_t nz);
31  inline void destroy();
32 
33  public:
34  inline FastSparseMatrix();
35  inline FastSparseMatrix(size_t n,size_t p, size_t sp);
36  inline FastSparseMatrix( const SparseMatrix &M);
37  inline FastSparseMatrix( const FastSparseMatrix &M);
38  inline ~FastSparseMatrix() {destroy();}
39  inline size_t nlin() const ;
40  inline size_t ncol() const ;
41  inline void write(std::ostream& f) const;
42  inline void read(std::istream& f);
43 
44  inline double operator()(size_t i,size_t j) const;
45  inline double& operator()(size_t i,size_t j);
46  inline Vector operator * (const Vector &v) const;
47  inline void operator =( const FastSparseMatrix &M);
48 
49  inline double& operator[](size_t i) {return tank[i];};
50 
51  inline void info() const;
52 
53  };
54 
55  inline std::ostream& operator<<(std::ostream& f,const FastSparseMatrix &M)
56  {
57  size_t nz = M.rowindex[M.nlin()];
58  f << M.nlin() << " " << M.ncol() << std::endl;
59  f << nz << std::endl;
60  for(size_t i=0;i<M.nlin();i++)
61  {
62  for(size_t j=M.rowindex[i];j<M.rowindex[i+1];j++)
63  {
64  f<<(long unsigned int)i<<"\t"<<(long unsigned int)M.js[j]<<"\t"<<M.tank[j]<<std::endl;
65  }
66  }
67  return f;
68  }
69 
70  inline void FastSparseMatrix::info() const {
71  if ((nlin() == 0) && (ncol() == 0)) {
72  std::cout << "Matrix Empty" << std::endl;
73  return;
74  }
75 
76  std::cout << "Dimensions : " << nlin() << " x " << ncol() << std::endl;
77  std::cout << *this;
78  }
79 
81  {
82  alloc(1,1,1);
83  }
84 
85  inline FastSparseMatrix::FastSparseMatrix(size_t n,size_t p, size_t sp=1)
86  {
87  alloc(n,p,sp);
88  }
89 
91  {
92  destroy();
93  alloc(M.m_nlin,M.m_ncol,M.rowindex[M.m_nlin]);
94  memcpy(tank,M.tank,sizeof(double)*M.rowindex[M.m_nlin]);
95  memcpy(js,M.js,sizeof(size_t)*M.rowindex[M.m_nlin]);
96  memcpy(rowindex,M.rowindex,sizeof(size_t)*(M.m_nlin+1));
97  }
98 
100  {
101  tank=new double[M.size()];
102  js=new size_t[M.size()];
103  rowindex=new size_t[M.nlin()+1];
104  m_nlin=(size_t)M.nlin();
105  m_ncol=(size_t)M.ncol();
106 
107  // we fill a data structure faster for computation
109  int cnt = 0;
110  size_t current_line = (size_t)-1;
111  for( it = M.begin(); it != M.end(); ++it) {
112  size_t i = it->first.first;
113  size_t j = it->first.second;
114  double val = it->second;
115  tank[cnt] = val;
116  js[cnt] = j;
117  if(i != current_line) {
118  for(size_t k = current_line+1; k <= i; ++k) {
119  rowindex[k]=cnt;
120  }
121  current_line = i;
122  }
123  cnt++;
124  }
125 
126  for(size_t k = current_line+1; k <= M.nlin(); ++k) {
127  rowindex[k]=M.size();
128  }
129 
130  }
131 
132  inline void FastSparseMatrix::write(std::ostream& f) const
133  {
134  size_t nz=rowindex[m_nlin];
135  f.write((const char*)&m_nlin,(std::streamsize)sizeof(size_t));
136  f.write((const char*)&m_ncol,(std::streamsize)sizeof(size_t));
137  f.write((const char*)&nz,(std::streamsize)sizeof(size_t));
138  f.write((const char*)tank,(std::streamsize)(sizeof(double)*nz));
139  f.write((const char*)js,(std::streamsize)(sizeof(size_t)*nz));
140  f.write((const char*)rowindex,(std::streamsize)(sizeof(size_t)*m_nlin));
141  }
142 
143  inline void FastSparseMatrix::read(std::istream& f)
144  {
145  destroy();
146  size_t nz;
147  f.read((char*)&m_nlin,(std::streamsize)sizeof(size_t));
148  f.read((char*)&m_ncol,(std::streamsize)sizeof(size_t));
149  f.read((char*)&nz,(std::streamsize)sizeof(size_t));
150  alloc(m_nlin,m_ncol,nz);
151  f.read((char*)tank,(std::streamsize)(sizeof(double)*nz));
152  f.read((char*)js,(std::streamsize)(sizeof(size_t)*nz));
153  f.read((char*)rowindex,(std::streamsize)(sizeof(size_t)*m_nlin));
154  }
155 
156  inline void FastSparseMatrix::alloc(size_t nl, size_t nc, size_t nz)
157  {
158  m_nlin=nl;
159  m_ncol=nc;
160  tank=new double[nz];
161  js=new size_t[nz];
162  rowindex=new size_t[nl+1];
163  rowindex[nl]=nz;
164  }
165 
167  {
168  delete[] tank;
169  delete[] js;
170  delete[] rowindex;
171  }
172 
174  {
175  alloc(m.m_nlin,m.m_ncol,m.rowindex[m.m_nlin]);
176  memcpy(tank,m.tank,sizeof(double)*m.rowindex[m.m_nlin]);
177  memcpy(js,m.js,sizeof(size_t)*m.rowindex[m.m_nlin]);
178  memcpy(rowindex,m.rowindex,sizeof(size_t)*(m.m_nlin+1));
179  }
180 
181  inline size_t FastSparseMatrix::nlin() const {return (size_t)m_nlin;}
182 
183  inline size_t FastSparseMatrix::ncol() const {return (size_t)m_ncol;}
184 
185  inline double FastSparseMatrix::operator()(size_t i,size_t j) const
186  {
187  for(size_t k=rowindex[i];k<rowindex[i+1];k++)
188  {
189  if(js[k]<j) continue;
190  else if(js[k]==j) return tank[k];
191  else break;
192  }
193 
194  return 0;
195  }
196 
197  inline double& FastSparseMatrix::operator()(size_t i,size_t j)
198  {
199  for(size_t k=rowindex[i];k<rowindex[i+1];k++)
200  {
201  if(js[k]<j) continue;
202  else if(js[k]==j) return tank[k];
203  else break;
204  }
205 
206  throw OpenMEEG::maths::BadSparseOperation("FastSparseMatrix : double& operator()(size_t i,size_t j) can't add element");
207  }
208 
210  {
211  Vector result(m_nlin); result.set(0);
212  double *pt_result=&result(0);
213  Vector *_v=(Vector *)&v;
214  double *pt_vect=&(*_v)(0);
215 
216  for(size_t i=0;i<m_nlin;i++)
217  {
218  double& total=pt_result[i];
219  for(size_t j=rowindex[i];j<rowindex[i+1];j++) {
220  total+=tank[j]*pt_vect[js[j]];
221  }
222  }
223  return result;
224  }
225 }
double & operator[](size_t i)
void write(std::ostream &f) const
Vector operator*(const Vector &v) const
double operator()(size_t i, size_t j) const
void alloc(size_t nl, size_t nc, size_t nz)
void operator=(const FastSparseMatrix &M)
void read(std::istream &f)
Dimension nlin() const
Definition: linop.h:48
virtual Dimension ncol() const
Definition: linop.h:51
const_iterator end() const
Definition: sparse_matrix.h:56
size_t size() const
Definition: sparse_matrix.h:51
const_iterator begin() const
Definition: sparse_matrix.h:55
Tank::const_iterator const_iterator
Definition: sparse_matrix.h:30
void set(const double x)
Vect3 operator*(const double d, const Vect3 &V)
Definition: vect3.h:105
std::ostream & operator<<(std::ostream &os, const Conductivity< REP > &m)