fmatvec 0.0.0

var_symmetric_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 var_symmetric_matrix_h
23 #define var_symmetric_matrix_h
24
25 #include "types.h"
26 #include "range.h"
27 #include "indices.h"
28
29 namespace fmatvec {
30
39 template <class AT> class Matrix<Symmetric,Var,Var,AT> {
40
41 protected:
42
44
45 int M{0};
46
47 AT *ele;
48
49 template <class Type, class Row, class Col> inline Matrix<Symmetric,Var,Var,AT>& copy(const Matrix<Type,Row,Col,AT> &A);
50 template <class Row> inline Matrix<Symmetric,Var,Var,AT>& copy(const Matrix<Symmetric,Row,Row,AT> &A);
51
53
54 public:
55 static constexpr bool isVector {false};
56 using value_type = AT;
57 using shape_type = Symmetric;
58
63 explicit Matrix() : ele(nullptr) {
64 }
65
66 explicit Matrix(int m, Noinit) : M(m), ele(new AT[M*M]) { }
67 explicit Matrix(int m, Init ini=INIT, const AT &a=AT()) : M(m), ele(new AT[M*M]) { init(a); }
68 explicit Matrix(int m, Eye ini, const AT &a=1) : M(m), ele(new AT[M*M]) { init(ini,a); }
69 explicit Matrix(int m, int n, Noinit) : M(m), ele(new AT[M*M]) { }
70 explicit Matrix(int m, int n, Init ini=INIT, const AT &a=AT()) : M(m), ele(new AT[M*M]) { init(a); }
71 explicit Matrix(int m, int n, Eye ini, const AT &a=1) : M(m), ele(new AT[M*M]) { init(ini,a); }
72
73 // move
74 Matrix(Matrix<Symmetric,Var,Var,AT> &&src) noexcept {
75 M=src.M;
76 src.M=0;
77 ele=src.ele;
78 src.ele=nullptr;
79 }
80 Matrix<Symmetric,Var,Var,AT>& operator=(Matrix<Symmetric,Var,Var,AT> &&src) noexcept {
81 FMATVEC_ASSERT(M == src.M, AT);
82 src.M=0;
83 delete[]ele;
84 ele=src.ele;
85 src.ele=nullptr;
86 return *this;
87 }
88
94 Matrix(const Matrix<Symmetric,Var,Var,AT> &A) : M(A.M), ele(new AT[M*M]) {
95 copy(A);
96 }
97
103 template<class Row>
104 Matrix(const Matrix<Symmetric,Row,Row,AT> &A) : M(A.size()), ele(new AT[M*M]) {
105 copy(A);
106 }
107
113 template<class Type, class Row, class Col>
114 explicit Matrix(const Matrix<Type,Row,Col,AT> &A) : M(A.rows()), ele(new AT[M*M]) {
115 FMATVEC_ASSERT(A.rows() == A.cols(), AT);
116 copy(A);
117 }
118
121 ~Matrix() {
122 delete[] ele;
123 }
124
125 Matrix<Symmetric,Var,Var,AT>& resize(int m, Noinit) {
126 delete[] ele;
127 M = m;
128 ele = new AT[M*M];
129 return *this;
130 }
131
132 Matrix<Symmetric,Var,Var,AT>& resize(int m, Init ini=INIT, const AT &a=AT()) { return resize(m,Noinit()).init(a); }
133
134 Matrix<Symmetric,Var,Var,AT>& resize(int m, Eye ini, const AT &a=1) { return resize(m,Noinit()).init(ini,a); }
135
142 inline Matrix<Symmetric,Var,Var,AT>& operator=(const Matrix<Symmetric,Var,Var,AT> &A) {
143 FMATVEC_ASSERT(M == A.M, AT);
144 return copy(A);
145 }
146
153 template<class Type, class Row, class Col>
154 inline Matrix<Symmetric,Var,Var,AT>& operator=(const Matrix<Type,Row,Col,AT> &A) {
155 FMATVEC_ASSERT(A.rows() == A.cols(), AT);
156 FMATVEC_ASSERT(M == A.rows(), AT);
157 return copy(A);
158 }
159
166 template<class Type, class Row, class Col>
167 inline Matrix<Symmetric,Var,Var,AT>& operator<<=(const Matrix<Type,Row,Col,AT> &A) {
168 FMATVEC_ASSERT(A.rows() == A.cols(), AT);
169 if(M!=A.rows()) resize(A.rows(),NONINIT);
170 return copy(A);
171 }
172 // move
173 inline Matrix<Symmetric,Var,Var,AT>& operator<<=(Matrix<Symmetric,Var,Var,AT> &&src) {
174 FMATVEC_ASSERT(src.rows() == src.cols(), AT);
175 M=src.M;
176 src.M=0;
177 delete[]ele;
178 ele=src.ele;
179 src.ele=nullptr;
180 return *this;
181 }
182
185 void resize(int n, int m) {
186 if(n!=m)
187 throw std::runtime_error("A symmetric matrix cannot have different dimensions for rows and columns.");
188 resize(n);
189 }
190
200 AT& operator()(int i, int j) {
201 FMATVEC_ASSERT(i>=0, AT);
202 FMATVEC_ASSERT(j>=0, AT);
203 FMATVEC_ASSERT(i<M, AT);
204 FMATVEC_ASSERT(j<M, AT);
205
206 return e(i,j);
207 }
208
213 const AT& operator()(int i, int j) const {
214 FMATVEC_ASSERT(i>=0, AT);
215 FMATVEC_ASSERT(j>=0, AT);
216 FMATVEC_ASSERT(i<M, AT);
217 FMATVEC_ASSERT(j<M, AT);
218
219 return e(i,j);
220 }
221
222 AT& ei(int i, int j) {
223 return ele[i*M+j];
224 }
225
226 const AT& ei(int i, int j) const {
227 return ele[i*M+j];
228 }
229
230 AT& ej(int i, int j) {
231 return ele[i+j*M];
232 }
233
234 const AT& ej(int i, int j) const {
235 return ele[i+j*M];
236 }
237
238 AT& e(int i, int j) {
239 return j > i ? ej(i,j) : ei(i,j);
240 }
241
242 const AT& e(int i, int j) const {
243 return j > i ? ej(i,j) : ei(i,j);
244 }
245
251 AT* operator()() {return ele;}
252
257 const AT* operator()() const {return ele;}
258
263 constexpr int size() const {return M;}
264
269 constexpr int rows() const {return M;}
270
275 constexpr int cols() const {return M;}
276
281 int ldim() const {return M;}
282
290 const CBLAS_ORDER blasOrder() const {
291 return CblasRowMajor;
292 }
293
301 const CBLAS_UPLO blasUplo() const {
302 return CblasLower;
303 }
304
312 inline Matrix<Symmetric,Var,Var,AT>& init(const AT &val=AT());
313 inline Matrix<Symmetric,Var,Var,AT>& init(Init, const AT &a=AT()) { return init(a); }
314 inline Matrix<Symmetric,Var,Var,AT>& init(Eye eye, const AT &val=1);
315 inline Matrix<Symmetric,Var,Var,AT>& init(Noinit, const AT &a=AT()) { return *this; }
316
319 inline const Matrix<General,Var,Var,AT> operator()(const Range<Var,Var> &I, const Range<Var,Var> &J) const;
320
323 inline const Matrix<Symmetric,Var,Var,AT> operator()(const Range<Var,Var> &I) const;
324
325 template<class Type, class Row, class Col> inline void set(const Range<Var,Var> &I, const Range<Var,Var> &J, const Matrix<Type,Row,Col,AT> &A);
326 template<class Type, class Row, class Col> inline void add(const Range<Var,Var> &I, const Range<Var,Var> &J, const Matrix<Type,Row,Col,AT> &A);
327
328 template<class Row> inline void set(const Range<Var,Var> &I, const Matrix<Symmetric,Row,Row,AT> &A);
329 template<class Row> inline void add(const Range<Var,Var> &I, const Matrix<Symmetric,Row,Row,AT> &A);
330
331 inline const Matrix<General,Var,Var,AT> operator()(const Indices &I, const Indices &J) const;
332
333 inline const Matrix<Symmetric,Var,Var,AT> operator()(const Indices &I) const;
334
339 explicit inline operator std::vector<std::vector<AT>>() const;
340
346 explicit Matrix(const std::vector<std::vector<AT>> &m);
347
348 inline int nonZeroElements() const;
349 };
350
351 template <class AT>
352 inline Matrix<Symmetric,Var,Var,AT>& Matrix<Symmetric,Var,Var,AT>::init(const AT &val) {
353
354 for(int i=0; i<M; i++)
355 for(int j=i; j<M; j++)
356 ej(i,j) = val;
357
358 return *this;
359 }
360
361 template <class AT>
362 inline Matrix<Symmetric,Var,Var,AT>& Matrix<Symmetric,Var,Var,AT>::init(Eye, const AT &val) {
363
364 for(int i=0; i<size(); i++) {
365 ej(i,i) = val;
366 for(int j=i+1; j<size(); j++) {
367 ej(i,j) = 0;
368 }
369 }
370 return *this;
371 }
372
373 template <class AT>
374 inline const Matrix<General,Var,Var,AT> Matrix<Symmetric,Var,Var,AT>::operator()(const Range<Var,Var> &I, const Range<Var,Var> &J) const {
375 FMATVEC_ASSERT(I.end()<M, AT);
376 FMATVEC_ASSERT(J.end()<M, AT);
377 Matrix<General,Var,Var,AT> A(I.size(),J.size(),NONINIT);
378
379 for(int i=0; i<A.rows(); i++)
380 for(int j=0; j<A.cols(); j++)
381 A.e(i,j) = e(I.start()+i,J.start()+j);
382
383 return A;
384 }
385
386 template <class AT>
387 inline const Matrix<Symmetric,Var,Var,AT> Matrix<Symmetric,Var,Var,AT>::operator()(const Range<Var,Var> &I) const {
388 FMATVEC_ASSERT(I.end()<M, AT);
389
390 Matrix<Symmetric,Var,Var,AT> A(I.size(),NONINIT);
391
392 for(int i=0; i<A.rows(); i++)
393 for(int j=i; j<A.cols(); j++)
394 A.e(i,j) = e(I.start()+i,I.start()+j);
395
396 return A;
397 }
398
399 template <class AT> template<class Type, class Row, class Col>
400 inline void Matrix<Symmetric,Var,Var,AT>::set(const Range<Var,Var> &I, const Range<Var,Var> &J, const Matrix<Type,Row,Col,AT> &A) {
401
402 if(I.start()>=J.start()) FMATVEC_ASSERT(I.start()>=J.end(), AT)
403 else FMATVEC_ASSERT(J.start()>=I.end(), AT);
404 FMATVEC_ASSERT(I.end()<rows(), AT);
405 FMATVEC_ASSERT(J.end()<cols(), AT);
406 FMATVEC_ASSERT(I.size()==A.rows(), AT);
407 FMATVEC_ASSERT(J.size()==A.cols(), AT);
408
409 for(int i=I.start(), ii=0; i<=I.end(); i++, ii++)
410 for(int j=J.start(), jj=0; j<=J.end(); j++, jj++)
411 e(i,j) = A.e(ii,jj);
412 }
413
414 template <class AT> template<class Type, class Row, class Col>
415 inline void Matrix<Symmetric,Var,Var,AT>::add(const Range<Var,Var> &I, const Range<Var,Var> &J, const Matrix<Type,Row,Col,AT> &A) {
416
417 if(I.start()>=J.start()) FMATVEC_ASSERT(I.start()>=J.end(), AT)
418 else FMATVEC_ASSERT(J.start()>=I.end(), AT);
419 FMATVEC_ASSERT(I.end()<rows(), AT);
420 FMATVEC_ASSERT(J.end()<cols(), AT);
421 FMATVEC_ASSERT(I.size()==A.rows(), AT);
422 FMATVEC_ASSERT(J.size()==A.cols(), AT);
423
424 for(int i=I.start(), ii=0; i<=I.end(); i++, ii++)
425 for(int j=J.start(), jj=0; j<=J.end(); j++, jj++)
426 e(i,j) += A.e(ii,jj);
427 }
428
429 template <class AT> template<class Row>
430 inline void Matrix<Symmetric,Var,Var,AT>::set(const Range<Var,Var> &I, const Matrix<Symmetric,Row,Row,AT> &A) {
431
432 FMATVEC_ASSERT(I.end()<size(), AT);
433 FMATVEC_ASSERT(I.size()==A.size(), AT);
434
435 for(int i=I.start(), ii=0; i<=I.end(); i++, ii++)
436 for(int j=i, jj=ii; j<=I.end(); j++, jj++)
437 e(i,j) = A.e(ii,jj);
438 }
439
440 template <class AT> template<class Row>
441 inline void Matrix<Symmetric,Var,Var,AT>::add(const Range<Var,Var> &I, const Matrix<Symmetric,Row,Row,AT> &A) {
442
443 FMATVEC_ASSERT(I.end()<size(), AT);
444 FMATVEC_ASSERT(I.size()==A.size(), AT);
445
446 for(int i=I.start(), ii=0; i<=I.end(); i++, ii++)
447 for(int j=i, jj=ii; j<=I.end(); j++, jj++)
448 e(i,j) += A.e(ii,jj);
449 }
450
451 template <class AT>
452 inline const Matrix<General,Var,Var,AT> Matrix<Symmetric,Var,Var,AT>::operator()(const Indices &I, const Indices &J) const {
453 FMATVEC_ASSERT(I.max()<size(), AT);
454 FMATVEC_ASSERT(J.max()<size(), AT);
455
456 Matrix<General,Var,Var,AT> A(I.size(),J.size(),NONINIT);
457
458 for(int i=0; i<A.rows(); i++)
459 for(int j=0; j<A.cols(); j++)
460 A.e(i,j) = e(I[i],J[j]);
461
462 return A;
463 }
464
465 template <class AT>
466 inline const Matrix<Symmetric,Var,Var,AT> Matrix<Symmetric,Var,Var,AT>::operator()(const Indices &I) const {
467 FMATVEC_ASSERT(I.max()<size(), AT);
468
469 Matrix<Symmetric,Var,Var,AT> A(I.size(),NONINIT);
470
471 for(int i=0; i<A.rows(); i++)
472 for(int j=0; j<A.cols(); j++)
473 A.e(i,j) = e(I[i],I[j]);
474
475 return A;
476 }
477
478 template <class AT>
479 inline Matrix<Symmetric,Var,Var,AT>::operator std::vector<std::vector<AT>>() const {
480 std::vector<std::vector<AT>> ret(rows(),std::vector<AT>(cols()));
481 for(int r=0; r<rows(); r++) {
482 for(int c=r; c<cols(); c++) {
483 ret[r][c]=ej(r,c);
484 if(c>r)
485 ret[c][r]=ej(r,c);
486 }
487 }
488 return ret;
489 }
490
491 template <class AT>
492 inline Matrix<Symmetric,Var,Var,AT>::Matrix(const std::vector<std::vector<AT>> &m) : M(static_cast<int>(m.size())), ele(new AT[M*M]) {
493 for(int r=0; r<rows(); r++) {
494 if(static_cast<int>(m[r].size())!=cols())
495 throw std::runtime_error("The rows of the input have different length.");
496 for(int c=r; c<cols(); c++) {
497 ej(r,c)=m[r][c];
498 if(c>r && abs(m[r][c]-m[c][r])>abs(m[r][c]*1e-13+1e-13))
499 throw std::runtime_error("The input is not symmetric.");
500 }
501 }
502 }
503
505
506 template <class AT> template <class Type, class Row, class Col>
507 inline Matrix<Symmetric,Var,Var,AT>& Matrix<Symmetric,Var,Var,AT>::copy(const Matrix<Type,Row,Col,AT> &A) {
508 for(int i=0; i<M; i++)
509 for(int j=i; j<M; j++)
510 ej(i,j) = A.e(i,j);
511 return *this;
512 }
513
514 template <class AT> template <class Row>
515 inline Matrix<Symmetric,Var,Var,AT>& Matrix<Symmetric,Var,Var,AT>::copy(const Matrix<Symmetric,Row,Row,AT> &A) {
516 for(int i=0; i<M; i++)
517 for(int j=i; j<M; j++)
518 ej(i,j) = A.ej(i,j);
519 return *this;
520 }
521
522 template <class AT>
523 inline int Matrix<Symmetric,Var,Var,AT>::nonZeroElements() const {
524 int k = M;
525 for (int j = 0; j < M; j++) {
526 for (int i = j+1; i < M; i++) {
527 if (fabs(e(i, j)) > 1e-16) {
528 k+=2;
529 }
530 }
531 }
532 return k;
533 }
534
536
537 }
538
539 #endif
fmatvec::Range< Var, Var >::end
int end() const
Last element.
Definition: range.h:91
fmatvec::Matrix< Symmetric, Var, Var, AT >
This is a matrix class for symmetric matrices.
Definition: var_symmetric_matrix.h:39
fmatvec::Matrix
This is the basic matrix class for arbitrary matrices.
Definition: fmatvec.h:42
fmatvec::Matrix< Symmetric, Var, Var, AT >::blasUplo
const CBLAS_UPLO blasUplo() const
Symmetry convention.
Definition: var_symmetric_matrix.h:301
fmatvec::Symmetric
Shape class for symmetric matrices.
Definition: types.h:128
fmatvec::Matrix< General, Var, Var, AT >
This is a matrix class for general matrices.
Definition: var_general_matrix.h:41
fmatvec::Matrix::operator()
AT & operator()(int i, int j)
Standard constructor.
Definition: matrix.h:84
fmatvec::Init
Definition: types.h:89
fmatvec::Var
Definition: types.h:101
fmatvec::Matrix< Symmetric, Var, Var, AT >::rows
constexpr int rows() const
Number of rows.
Definition: var_symmetric_matrix.h:269
fmatvec::Range< Var, Var >::start
int start() const
First element.
Definition: range.h:85
fmatvec::Matrix< Symmetric, Var, Var, AT >::blasOrder
const CBLAS_ORDER blasOrder() const
Storage convention.
Definition: var_symmetric_matrix.h:290
fmatvec::Matrix::cols
int cols() const
Number of columns.
fmatvec::Matrix< Symmetric, Var, Var, AT >::Matrix
Matrix(const Matrix< Symmetric, Var, Var, AT > &A)
Copy Constructor.
Definition: var_symmetric_matrix.h:94
fmatvec::Matrix< Symmetric, Var, Var, AT >::operator()
AT * operator()()
Pointer operator.
Definition: var_symmetric_matrix.h:251
fmatvec::Matrix< Symmetric, Var, Var, AT >::operator=
Matrix< Symmetric, Var, Var, AT > & operator=(const Matrix< Symmetric, Var, Var, AT > &A)
Assignment operator.
Definition: var_symmetric_matrix.h:142
fmatvec::Matrix< Symmetric, Var, Var, AT >::operator()
AT & operator()(int i, int j)
Element operator.
Definition: var_symmetric_matrix.h:200
fmatvec::Matrix< Symmetric, Var, Var, AT >::operator=
Matrix< Symmetric, Var, Var, AT > & operator=(const Matrix< Type, Row, Col, AT > &A)
Assignment operator.
Definition: var_symmetric_matrix.h:154
fmatvec::Matrix< Symmetric, Var, Var, AT >::operator()
const AT * operator()() const
Pointer operator.
Definition: var_symmetric_matrix.h:257
fmatvec
Namespace fmatvec.
Definition: _memory.cc:28
fmatvec::Matrix< Symmetric, Var, Var, AT >::ldim
int ldim() const
Leading dimension.
Definition: var_symmetric_matrix.h:281
fmatvec::Matrix< Symmetric, Var, Var, AT >::~Matrix
~Matrix()
Destructor.
Definition: var_symmetric_matrix.h:121
fmatvec::Matrix< Symmetric, Var, Var, AT >::Matrix
Matrix(const Matrix< Symmetric, Row, Row, AT > &A)
Copy Constructor.
Definition: var_symmetric_matrix.h:104
fmatvec::Matrix< Symmetric, Var, Var, AT >::size
constexpr int size() const
Size.
Definition: var_symmetric_matrix.h:263
fmatvec::Range< Var, Var >::size
int size() const
Size.
Definition: range.h:97
fmatvec::Matrix< Symmetric, Var, Var, AT >::Matrix
Matrix()
Standard constructor.
Definition: var_symmetric_matrix.h:63
fmatvec::Matrix< Symmetric, Var, Var, AT >::Matrix
Matrix(const Matrix< Type, Row, Col, AT > &A)
Copy Constructor.
Definition: var_symmetric_matrix.h:114
fmatvec::Matrix< Symmetric, Var, Var, AT >::cols
constexpr int cols() const
Number of columns.
Definition: var_symmetric_matrix.h:275
fmatvec::Noinit
Definition: types.h:88
fmatvec::Matrix< Symmetric, Var, Var, AT >::resize
void resize(int n, int m)
Definition: var_symmetric_matrix.h:185
fmatvec::Matrix::rows
int rows() const
Number of rows.
fmatvec::Matrix< Symmetric, Var, Var, AT >::operator()
const AT & operator()(int i, int j) const
Element operator.
Definition: var_symmetric_matrix.h:213
fmatvec::Range< Var, Var >
This is an index class for creating submatrices.
Definition: range.h:44
fmatvec::Eye
Definition: types.h:90