OpenMEEG
Loading...
Searching...
No Matches
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
14namespace 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();
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 {
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)
std::ostream & operator<<(std::ostream &os, const Conductivity< REP > &m)