OpenMEEG
fast_sparse_matrix.h
Go to the documentation of this file.
1 /*
2 Project Name : OpenMEEG
3 
4 © INRIA and ENPC (contributors: Geoffray ADDE, Maureen CLERC, Alexandre
5 GRAMFORT, Renaud KERIVEN, Jan KYBIC, Perrine LANDREAU, Théodore PAPADOPOULO,
6 Emmanuel OLIVI
7 Maureen.Clerc.AT.inria.fr, keriven.AT.certis.enpc.fr,
8 kybic.AT.fel.cvut.cz, papadop.AT.inria.fr)
9 
10 The OpenMEEG software is a C++ package for solving the forward/inverse
11 problems of electroencephalography and magnetoencephalography.
12 
13 This software is governed by the CeCILL-B license under French law and
14 abiding by the rules of distribution of free software. You can use,
15 modify and/ or redistribute the software under the terms of the CeCILL-B
16 license as circulated by CEA, CNRS and INRIA at the following URL
17 "http://www.cecill.info".
18 
19 As a counterpart to the access to the source code and rights to copy,
20 modify and redistribute granted by the license, users are provided only
21 with a limited warranty and the software's authors, the holders of the
22 economic rights, and the successive licensors have only limited
23 liability.
24 
25 In this respect, the user's attention is drawn to the risks associated
26 with loading, using, modifying and/or developing or reproducing the
27 software by the user in light of its specific status of free software,
28 that may mean that it is complicated to manipulate, and that also
29 therefore means that it is reserved for developers and experienced
30 professionals having in-depth computer knowledge. Users are therefore
31 encouraged to load and test the software's suitability as regards their
32 requirements in conditions enabling the security of their systems and/or
33 data to be ensured and, more generally, to use and operate it in the
34 same conditions as regards security.
35 
36 The fact that you are presently reading this means that you have had
37 knowledge of the CeCILL-B license and that you accept its terms.
38 */
39 
40 #pragma once
41 
42 #include "OpenMEEGMathsConfig.h"
43 #include "vector.h"
44 #include "sparse_matrix.h"
45 
46 namespace OpenMEEG {
47 
48  class OPENMEEGMATHS_EXPORT FastSparseMatrix
49  {
50  public:
51 
52  inline friend std::ostream& operator<<(std::ostream& f,const FastSparseMatrix &M);
53 
54  protected:
55 
56  double *tank;
57  size_t *js;
58  size_t *rowindex;
59  size_t m_nlin;
60  size_t m_ncol;
61 
62  inline void alloc(size_t nl, size_t nc, size_t nz);
63  inline void destroy();
64 
65  public:
66  inline FastSparseMatrix();
67  inline FastSparseMatrix(size_t n,size_t p, size_t sp);
68  inline FastSparseMatrix( const SparseMatrix &M);
69  inline FastSparseMatrix( const FastSparseMatrix &M);
70  inline ~FastSparseMatrix() {destroy();}
71  inline size_t nlin() const ;
72  inline size_t ncol() const ;
73  inline void write(std::ostream& f) const;
74  inline void read(std::istream& f);
75 
76  inline double operator()(size_t i,size_t j) const;
77  inline double& operator()(size_t i,size_t j);
78  inline Vector operator * (const Vector &v) const;
79  inline void operator =( const FastSparseMatrix &M);
80 
81  inline double& operator[](size_t i) {return tank[i];};
82 
83  inline void info() const;
84 
85  };
86 
87  inline std::ostream& operator<<(std::ostream& f,const FastSparseMatrix &M)
88  {
89  size_t nz = M.rowindex[M.nlin()];
90  f << M.nlin() << " " << M.ncol() << std::endl;
91  f << nz << std::endl;
92  for(size_t i=0;i<M.nlin();i++)
93  {
94  for(size_t j=M.rowindex[i];j<M.rowindex[i+1];j++)
95  {
96  f<<(long unsigned int)i<<"\t"<<(long unsigned int)M.js[j]<<"\t"<<M.tank[j]<<std::endl;
97  }
98  }
99  return f;
100  }
101 
102  inline void FastSparseMatrix::info() const {
103  if ((nlin() == 0) && (ncol() == 0)) {
104  std::cout << "Matrix Empty" << std::endl;
105  return;
106  }
107 
108  std::cout << "Dimensions : " << nlin() << " x " << ncol() << std::endl;
109  std::cout << *this;
110  }
111 
113  {
114  alloc(1,1,1);
115  }
116 
117  inline FastSparseMatrix::FastSparseMatrix(size_t n,size_t p, size_t sp=1)
118  {
119  alloc(n,p,sp);
120  }
121 
123  {
124  destroy();
125  alloc(M.m_nlin,M.m_ncol,M.rowindex[M.m_nlin]);
126  memcpy(tank,M.tank,sizeof(double)*M.rowindex[M.m_nlin]);
127  memcpy(js,M.js,sizeof(size_t)*M.rowindex[M.m_nlin]);
128  memcpy(rowindex,M.rowindex,sizeof(size_t)*(M.m_nlin+1));
129  }
130 
132  {
133  tank=new double[M.size()];
134  js=new size_t[M.size()];
135  rowindex=new size_t[M.nlin()+1];
136  m_nlin=(size_t)M.nlin();
137  m_ncol=(size_t)M.ncol();
138 
139  // we fill a data structure faster for computation
141  int cnt = 0;
142  size_t current_line = (size_t)-1;
143  for( it = M.begin(); it != M.end(); ++it) {
144  size_t i = it->first.first;
145  size_t j = it->first.second;
146  double val = it->second;
147  tank[cnt] = val;
148  js[cnt] = j;
149  if(i != current_line) {
150  for(size_t k = current_line+1; k <= i; ++k) {
151  rowindex[k]=cnt;
152  }
153  current_line = i;
154  }
155  cnt++;
156  }
157 
158  for(size_t k = current_line+1; k <= M.nlin(); ++k) {
159  rowindex[k]=M.size();
160  }
161 
162  }
163 
164  inline void FastSparseMatrix::write(std::ostream& f) const
165  {
166  size_t nz=rowindex[m_nlin];
167  f.write((const char*)&m_nlin,(std::streamsize)sizeof(size_t));
168  f.write((const char*)&m_ncol,(std::streamsize)sizeof(size_t));
169  f.write((const char*)&nz,(std::streamsize)sizeof(size_t));
170  f.write((const char*)tank,(std::streamsize)(sizeof(double)*nz));
171  f.write((const char*)js,(std::streamsize)(sizeof(size_t)*nz));
172  f.write((const char*)rowindex,(std::streamsize)(sizeof(size_t)*m_nlin));
173  }
174 
175  inline void FastSparseMatrix::read(std::istream& f)
176  {
177  destroy();
178  size_t nz;
179  f.read((char*)&m_nlin,(std::streamsize)sizeof(size_t));
180  f.read((char*)&m_ncol,(std::streamsize)sizeof(size_t));
181  f.read((char*)&nz,(std::streamsize)sizeof(size_t));
182  alloc(m_nlin,m_ncol,nz);
183  f.read((char*)tank,(std::streamsize)(sizeof(double)*nz));
184  f.read((char*)js,(std::streamsize)(sizeof(size_t)*nz));
185  f.read((char*)rowindex,(std::streamsize)(sizeof(size_t)*m_nlin));
186  }
187 
188  inline void FastSparseMatrix::alloc(size_t nl, size_t nc, size_t nz)
189  {
190  m_nlin=nl;
191  m_ncol=nc;
192  tank=new double[nz];
193  js=new size_t[nz];
194  rowindex=new size_t[nl+1];
195  rowindex[nl]=nz;
196  }
197 
199  {
200  delete[] tank;
201  delete[] js;
202  delete[] rowindex;
203  }
204 
206  {
207  alloc(m.m_nlin,m.m_ncol,m.rowindex[m.m_nlin]);
208  memcpy(tank,m.tank,sizeof(double)*m.rowindex[m.m_nlin]);
209  memcpy(js,m.js,sizeof(size_t)*m.rowindex[m.m_nlin]);
210  memcpy(rowindex,m.rowindex,sizeof(size_t)*(m.m_nlin+1));
211  }
212 
213  inline size_t FastSparseMatrix::nlin() const {return (size_t)m_nlin;}
214 
215  inline size_t FastSparseMatrix::ncol() const {return (size_t)m_ncol;}
216 
217  inline double FastSparseMatrix::operator()(size_t i,size_t j) const
218  {
219  for(size_t k=rowindex[i];k<rowindex[i+1];k++)
220  {
221  if(js[k]<j) continue;
222  else if(js[k]==j) return tank[k];
223  else break;
224  }
225 
226  return 0;
227  }
228 
229  inline double& FastSparseMatrix::operator()(size_t i,size_t j)
230  {
231  for(size_t k=rowindex[i];k<rowindex[i+1];k++)
232  {
233  if(js[k]<j) continue;
234  else if(js[k]==j) return tank[k];
235  else break;
236  }
237 
238  std::cerr<<"FastSparseMatrix : double& operator()(size_t i,size_t j) can't add element"<<std::endl;
239  exit(1);
240  }
241 
243  {
244  Vector result(m_nlin); result.set(0);
245  double *pt_result=&result(0);
246  Vector *_v=(Vector *)&v;
247  double *pt_vect=&(*_v)(0);
248 
249  for(size_t i=0;i<m_nlin;i++)
250  {
251  double& total=pt_result[i];
252  for(size_t j=rowindex[i];j<rowindex[i+1];j++) {
253  total+=tank[j]*pt_vect[js[j]];
254  }
255  }
256  return result;
257  }
258 }
OpenMEEG::FastSparseMatrix::rowindex
size_t * rowindex
Definition: fast_sparse_matrix.h:58
OpenMEEG::FastSparseMatrix::m_ncol
size_t m_ncol
Definition: fast_sparse_matrix.h:60
sparse_matrix.h
OpenMEEG::FastSparseMatrix::FastSparseMatrix
FastSparseMatrix()
Definition: fast_sparse_matrix.h:112
OpenMEEG::FastSparseMatrix::alloc
void alloc(size_t nl, size_t nc, size_t nz)
Definition: fast_sparse_matrix.h:188
OpenMEEG::SparseMatrix::size
size_t size() const
Definition: sparse_matrix.h:89
OpenMEEG::FastSparseMatrix::destroy
void destroy()
Definition: fast_sparse_matrix.h:198
OpenMEEG::Vector::set
void set(double x)
OpenMEEG::FastSparseMatrix::operator()
double operator()(size_t i, size_t j) const
Definition: fast_sparse_matrix.h:217
OpenMEEG::LinOpInfo::nlin
size_t nlin() const
Definition: linop.h:77
OpenMEEG
Definition: analytics.h:45
OpenMEEG::SparseMatrix
Definition: sparse_matrix.h:62
OpenMEEG::FastSparseMatrix::~FastSparseMatrix
~FastSparseMatrix()
Definition: fast_sparse_matrix.h:70
OpenMEEG::FastSparseMatrix::read
void read(std::istream &f)
Definition: fast_sparse_matrix.h:175
OpenMEEG::FastSparseMatrix::info
void info() const
Definition: fast_sparse_matrix.h:102
OpenMEEG::SparseMatrix::const_iterator
std::map< std::pair< size_t, size_t >, double >::const_iterator const_iterator
Definition: sparse_matrix.h:67
OpenMEEG::FastSparseMatrix::tank
double * tank
Definition: fast_sparse_matrix.h:56
OpenMEEG::FastSparseMatrix::ncol
size_t ncol() const
Definition: fast_sparse_matrix.h:215
OpenMEEG::FastSparseMatrix
Definition: fast_sparse_matrix.h:49
OpenMEEG::Vector
Definition: vector.h:56
OpenMEEG::FastSparseMatrix::write
void write(std::ostream &f) const
Definition: fast_sparse_matrix.h:164
OpenMEEG::LinOpInfo::ncol
virtual size_t ncol() const
Definition: linop.h:80
OpenMEEG::FastSparseMatrix::js
size_t * js
Definition: fast_sparse_matrix.h:57
OpenMEEGMathsConfig.h
OpenMEEG::SparseMatrix::end
const_iterator end() const
Definition: sparse_matrix.h:94
vector.h
OpenMEEG::FastSparseMatrix::operator*
Vector operator*(const Vector &v) const
Definition: fast_sparse_matrix.h:242
OpenMEEG::SparseMatrix::begin
const_iterator begin() const
Definition: sparse_matrix.h:93
OpenMEEG::operator<<
std::ostream & operator<<(std::ostream &os, const Conductivity< REP > &m)
Definition: PropertiesSpecialized.h:63
OpenMEEG::FastSparseMatrix::operator[]
double & operator[](size_t i)
Definition: fast_sparse_matrix.h:81
OpenMEEG::operator*
Vect3 operator*(const double &d, const Vect3 &v)
Definition: vect3.h:151
OpenMEEG::FastSparseMatrix::m_nlin
size_t m_nlin
Definition: fast_sparse_matrix.h:59
OpenMEEG::FastSparseMatrix::operator=
void operator=(const FastSparseMatrix &M)
Definition: fast_sparse_matrix.h:122
OpenMEEG::FastSparseMatrix::nlin
size_t nlin() const
Definition: fast_sparse_matrix.h:213