fmatvec 0.0.0

sparse_matrix.h
1 /* Copyright (C) 2003-2005 Martin Förg
2
3 * This library is free software; you can redistribute it and/or
4 * modify it under the terms of the GNU Lesser General Public
5 * License as published by the Free Software Foundation; either
6 * version 2.1 of the License, or (at your option) any later version.
7 *
8 * This library is distributed in the hope that it will be useful,
9 * but WITHOUT ANY WARRANTY; without even the implied warranty of
10 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
11 * Lesser General Public License for more details.
12 *
13 * You should have received a copy of the GNU Lesser General Public
14 * License along with this library; if not, write to the Free Software
15 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
16 *
17 * Contact:
18 * martin.o.foerg@googlemail.com
19 *
20 */
21
22 #ifndef sparse_matrix_h
23 #define sparse_matrix_h
24
25 #include "matrix.h"
26 #include "square_matrix.h"
27 #include "types.h"
28 #include "_memory.h"
29
30 namespace fmatvec {
31
43 template <class AT>
44 class Matrix<Sparse,Ref,Ref,AT> {
45
46 protected:
47
49
50 Memory<AT> memEle;
51 Memory<int> memI, memJ;
52 AT *ele;
53 int *I{nullptr}, *J{nullptr};
54 int m{0}, n{0}, k{0};
55
56 template <class Type, class Row, class Col> inline Matrix<Sparse,Ref,Ref,AT>& copy(const Matrix<Type,Row,Col,AT> &A);
57 template <class Row> inline Matrix<Sparse,Ref,Ref,AT>& copy(const Matrix<Symmetric,Row,Row,AT> &A);
58 inline Matrix<Sparse,Ref,Ref,AT>& copy(const Matrix<Sparse,Ref,Ref,AT> &A);
59
61
62
63 public:
64 static constexpr bool isVector {false};
65
66 using value_type = AT;
67 using shape_type = Sparse;
68
73 explicit Matrix() : memEle(), ele(0) { }
74
75 explicit Matrix(int m_, int n_, int k_, Noinit) : memEle(k_), memI(m_+1), memJ(k_), ele((AT*)memEle.get()), I((int*)memI.get()), J((int*)memJ.get()), m(m_), n(n_), k(k_) { }
76 explicit Matrix(int m_, int n_, int k_, Init ini=INIT, const AT &a=AT()) : memEle(k_), memI(m_+1), memJ(k_), ele((AT*)memEle.get()), I((int*)memI.get()), J((int*)memJ.get()), m(m_), n(n_), k(k_) { init(a); }
77
83 Matrix(const Matrix<Sparse,Ref,Ref,AT> &A) : memEle(A.k), memI(A.n+1), memJ(A.k), ele((AT*)memEle.get()), I((int*)memI.get()), J((int*)memJ.get()), m(A.m), n(A.n), k(A.k) {
84 copy(A);
85 }
86
87 explicit Matrix(int m_, int n_, int k_, int *I_, int *J_, Init ini=INIT, const AT &a=AT()) : memEle(k_), memI(), memJ(), ele((AT*)memEle.get()), I(I_), J(J_), m(m_), n(n_), k(k_) {
88 init(a);
89 }
90
93 ~Matrix() = default;
94
95 Matrix<Sparse,Ref,Ref,AT>& resize(int m_, int n_, int k_, Noinit) {
96 m = m_; n = n_; k = k_;
97 memEle.resize(k);
98 memI.resize(m+1);
99 memJ.resize(k);
100 ele = (AT*)memEle.get();
101 I = (int*)memI.get();
102 J = (int*)memJ.get();
103 return *this;
104 }
105
106 Matrix<Sparse,Ref,Ref,AT>& resize(int m, int n, int k, Init ini=INIT, const AT &a=AT()) { return resize(m,n,k,Noinit()).init(a); }
107
109 void resize(int m, int n) {
110 resize(m,n,m*n);
111 }
112
119 inline Matrix<Sparse,Ref,Ref,AT>& operator=(const Matrix<Sparse,Ref,Ref,AT> &A) {
120 FMATVEC_ASSERT(m == A.m, AT);
121 FMATVEC_ASSERT(n == A.n, AT);
122 FMATVEC_ASSERT(k == A.k, AT);
123 return copy(A);
124 }
125
132 template<class Type, class Row, class Col>
133 inline Matrix<Sparse,Ref,Ref,AT>& operator=(const Matrix<Type,Row,Col,AT> &A) {
134 FMATVEC_ASSERT(m == A.rows(), AT);
135 FMATVEC_ASSERT(n == A.cols(), AT);
136 FMATVEC_ASSERT(k == A.nonZeroElements(), AT);
137 return copy(A);
138 }
139
146 inline Matrix<Sparse,Ref,Ref,AT>& operator&=(Matrix<Sparse,Ref,Ref,AT> &A) {
147 m=A.m;
148 n=A.n;
149 k=A.k;
150 memEle = A.memEle;
151 memI = A.memI;
152 memJ = A.memJ;
153 ele = A.ele;
154 I = A.I;
155 J = A.J;
156 k = A.k;
157 return *this;
158 }
159
166 template<class Type, class Row, class Col>
167 inline Matrix<Sparse,Ref,Ref,AT>& operator<<=(const Matrix<Type,Row,Col,AT> &A) {
168 int k_ = A.nonZeroElements();
169 if(m!=A.rows() || k!=k_) resize(A.rows(),A.cols(),k_,NONINIT);
170 return copy(A);
171 }
172
177 const AT* operator()() const {
178 return ele;
179 }
180
186 AT* operator()() {
187 return ele;
188 }
189
190 AT& operator()(int i, int j) {
191 FMATVEC_ASSERT(i>=0, AT);
192 FMATVEC_ASSERT(j>=0, AT);
193 FMATVEC_ASSERT(i<m, AT);
194 FMATVEC_ASSERT(j<n, AT);
195
196 return e(i,j);
197 }
198
199 const AT& operator()(int i, int j) const {
200 FMATVEC_ASSERT(i>=0, AT);
201 FMATVEC_ASSERT(j>=0, AT);
202 FMATVEC_ASSERT(i<m, AT);
203 FMATVEC_ASSERT(j<n, AT);
204
205 return e(i,j);
206 }
207
208 AT& e(int i, int j) {
209 return ele[pos(i,j)];
210 }
211
212 const AT& e(int i, int j) const {
213 return ele[pos(i,j)];
214 }
215
216 int pos(int i, int j) const {
217 for(int k=I[i]; k<I[i+1]; k++)
218 if(J[k]==j)
219 return k;
220 return 0;
221 }
222
223 int ldim() {
224 throw std::runtime_error("Matrix<Sparse, Ref, Ref, AT>::ldim() cannot be called.");
225 }
226
231 const int* Ip() const {
232 return I;
233 }
234
239 int* Ip() {
240 return I;
241 }
242
247 const int* Jp() const {
248 return J;
249 }
250
255 int* Jp() {
256 return J;
257 }
258
263 int rows() const {return m;}
264
269 int cols() const {return n;}
270
273 const CBLAS_ORDER blasOrder() const {
274 throw std::runtime_error("Matrix<Sparse, Ref, Ref, AT>::blasOrder() cannot be called.");
275 }
276
284 Matrix<Sparse,Ref,Ref,AT>& init(const AT &val);
285 inline Matrix<Sparse,Ref,Ref,AT>& init(Init, const AT &a=AT()) { return init(a); }
286 inline Matrix<Sparse,Ref,Ref,AT>& init(Noinit, const AT &a=AT()) { return *this; }
287
288 int nonZeroElements() const { return k; }
289 };
290 // ------------------------- Constructors -------------------------------------
291 // ----------------------------------------------------------------------------
292
293 template <class AT> template <class Type, class Row, class Col>
294 inline Matrix<Sparse,Ref,Ref,AT>& Matrix<Sparse,Ref,Ref,AT>::copy(const Matrix<Type,Row,Col,AT> &A) {
295 int k=0;
296 int i;
297 for(i=0; i<A.rows(); i++) {
298 ele[k]=A.e(i,i);
299 J[k]=i;
300 I[i]=k++;
301 for(int j=0; j<i; j++) {
302 // TODO eps
303 if(fabs(A.e(i,j))>1e-16) {
304 ele[k]=A.e(i,j);
305 J[k++]=j;
306 }
307 }
308 for(int j=i+1; j<A.cols(); j++) {
309 // TODO eps
310 if(fabs(A.e(i,j))>1e-16) {
311 ele[k]=A.e(i,j);
312 J[k++]=j;
313 }
314 }
315 }
316 if(n) I[i]=k;
317 return *this;
318 }
319
320 template <class AT> template <class Row>
321 inline Matrix<Sparse,Ref,Ref,AT>& Matrix<Sparse,Ref,Ref,AT>::copy(const Matrix<Symmetric,Row,Row,AT> &A) {
322 int k=0;
323 int i;
324 for(i=0; i<A.size(); i++) {
325 ele[k]=A.ej(i,i);
326 J[k]=i;
327 I[i]=k++;
328 for(int j=0; j<i; j++) {
329 // TODO eps
330 if(fabs(A.ei(i,j))>1e-16) {
331 ele[k]=A.ei(i,j);
332 J[k++]=j;
333 }
334 }
335 for(int j=i+1; j<A.size(); j++) {
336 // TODO eps
337 if(fabs(A.ej(i,j))>1e-16) {
338 ele[k]=A.ej(i,j);
339 J[k++]=j;
340 }
341 }
342 }
343 if(n) I[i]=k;
344 return *this;
345 }
346
347 template <class AT>
348 inline Matrix<Sparse,Ref,Ref,AT>& Matrix<Sparse,Ref,Ref,AT>::copy(const Matrix<Sparse,Ref,Ref,AT> &A) {
349 for(int i=0; i<=m; i++) {
350 I[i] = A.I[i];
351 }
352 for(int i=0; i<k; i++) {
353 ele[i] = A.ele[i];
354 J[i] = A.J[i];
355 }
356 return *this;
357 }
358
359 template <class AT>
360 Matrix<Sparse,Ref,Ref,AT>& Matrix<Sparse,Ref,Ref,AT>::init(const AT& val) {
361 for(int i=0; i<k; i++) {
362 ele[i] = val;
363 }
364 return *this;
365 }
366
367 }
368
369 #endif
fmatvec::Matrix< Sparse, Ref, Ref, AT >::Matrix
Matrix()
Standard constructor.
Definition: sparse_matrix.h:73
fmatvec::Matrix
This is the basic matrix class for arbitrary matrices.
Definition: fmatvec.h:42
fmatvec::Matrix::operator()
AT & operator()(int i, int j)
Standard constructor.
Definition: matrix.h:84
fmatvec::Init
Definition: types.h:89
fmatvec::Matrix< Sparse, Ref, Ref, AT >::Jp
const int * Jp() const
Pointer operator.
Definition: sparse_matrix.h:247
fmatvec::Matrix< Sparse, Ref, Ref, AT >::operator()
AT * operator()()
Pointer operator.
Definition: sparse_matrix.h:186
fmatvec::Matrix< Sparse, Ref, Ref, AT >::rows
int rows() const
Number of rows.
Definition: sparse_matrix.h:263
fmatvec::Matrix< Sparse, Ref, Ref, AT >
This is a matrix class for sparse quadratic matrices.
Definition: sparse_matrix.h:44
fmatvec::Matrix< Sparse, Ref, Ref, AT >::blasOrder
const CBLAS_ORDER blasOrder() const
Storage convention.
Definition: sparse_matrix.h:273
fmatvec::Ref
Definition: types.h:98
fmatvec::Matrix< Sparse, Ref, Ref, AT >::Ip
int * Ip()
Pointer operator.
Definition: sparse_matrix.h:239
fmatvec::Matrix< Sparse, Ref, Ref, AT >::operator&=
Matrix< Sparse, Ref, Ref, AT > & operator&=(Matrix< Sparse, Ref, Ref, AT > &A)
Reference operator.
Definition: sparse_matrix.h:146
fmatvec::Matrix< Sparse, Ref, Ref, AT >::Jp
int * Jp()
Pointer operator.
Definition: sparse_matrix.h:255
fmatvec::Matrix< Sparse, Ref, Ref, AT >::cols
int cols() const
Number of columns.
Definition: sparse_matrix.h:269
fmatvec::Matrix::cols
int cols() const
Number of columns.
fmatvec::Matrix< Sparse, Ref, Ref, AT >::operator=
Matrix< Sparse, Ref, Ref, AT > & operator=(const Matrix< Sparse, Ref, Ref, AT > &A)
Assignment operator.
Definition: sparse_matrix.h:119
fmatvec::Matrix< Sparse, Ref, Ref, AT >::Ip
const int * Ip() const
Pointer operator.
Definition: sparse_matrix.h:231
fmatvec
Namespace fmatvec.
Definition: _memory.cc:28
fmatvec::Sparse
Shape class for sparse matrices.
Definition: types.h:152
fmatvec::Matrix< Sparse, Ref, Ref, AT >::Matrix
Matrix(const Matrix< Sparse, Ref, Ref, AT > &A)
Copy Constructor.
Definition: sparse_matrix.h:83
fmatvec::Noinit
Definition: types.h:88
fmatvec::Matrix< Sparse, Ref, Ref, AT >::operator()
const AT * operator()() const
Pointer operator.
Definition: sparse_matrix.h:177
fmatvec::Matrix< Sparse, Ref, Ref, AT >::resize
void resize(int m, int n)
Resize a sparse matrix.
Definition: sparse_matrix.h:109
fmatvec::Matrix< Sparse, Ref, Ref, AT >::operator=
Matrix< Sparse, Ref, Ref, AT > & operator=(const Matrix< Type, Row, Col, AT > &A)
Assignment operator.
Definition: sparse_matrix.h:133
fmatvec::Matrix::rows
int rows() const
Number of rows.