fmatvec 0.0.0

symmetric_sparse_matrix.h
1 /* Copyright (C) 2003-2022 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 symmetric_sparse_matrix_h
23 #define symmetric_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<SymmetricSparse,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 n{0}, k{0};
55
56 template <class Type, class Row, class Col> inline Matrix<SymmetricSparse,Ref,Ref,AT>& copy(const Matrix<Type,Row,Col,AT> &A);
57 template <class Row> inline Matrix<SymmetricSparse,Ref,Ref,AT>& copy(const Matrix<Symmetric,Row,Row,AT> &A);
58 inline Matrix<SymmetricSparse,Ref,Ref,AT>& copy(const Matrix<SymmetricSparse,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 = SymmetricSparse;
68
73 explicit Matrix() : memEle(), ele(0) { }
74
75 explicit Matrix(int n_, int k_, Noinit) : memEle(k_), memI(n_+1), memJ(k_), ele((AT*)memEle.get()), I((int*)memI.get()), J((int*)memJ.get()), n(n_), k(k_) { }
76 explicit Matrix(int n_, int k_, Init ini=INIT, const AT &a=AT()) : memEle(k_), memI(n_+1), memJ(k_), ele((AT*)memEle.get()), I((int*)memI.get()), J((int*)memJ.get()), n(n_), k(k_) { init(a); }
77
83 Matrix(const Matrix<SymmetricSparse,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()), n(A.n), k(A.k) {
84 copy(A);
85 }
86
87 explicit Matrix(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_), n(n_), k(k_) {
88 init(a);
89 }
90
93 ~Matrix() = default;
94
95 Matrix<SymmetricSparse,Ref,Ref,AT>& resize(int n_, int k_, Noinit) {
96 n = n_; k = k_;
97 memEle.resize(k);
98 memI.resize(n+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<SymmetricSparse,Ref,Ref,AT>& resize(int n, int k, Init ini=INIT, const AT &a=AT()) { return resize(n,k,Noinit()).init(a); }
107
109 void resize(int n, int m) {
110 if(n!=m)
111 throw std::runtime_error("A symmetric sparse matrix cannot have different dimensions for rows and columns.");
112 resize(n,n,n*n);
113 }
114
121 inline Matrix<SymmetricSparse,Ref,Ref,AT>& operator=(const Matrix<SymmetricSparse,Ref,Ref,AT> &A) {
122 FMATVEC_ASSERT(n == A.n, AT);
123 FMATVEC_ASSERT(k == A.k, AT);
124 return copy(A);
125 }
126
133 template<class Type, class Row, class Col>
134 inline Matrix<SymmetricSparse,Ref,Ref,AT>& operator=(const Matrix<Type,Row,Col,AT> &A) {
135 FMATVEC_ASSERT(A.rows() == A.cols(), AT);
136 FMATVEC_ASSERT(n == A.cols(), AT);
137 FMATVEC_ASSERT(k == A.nonZeroElements(), AT);
138 return copy(A);
139 }
140
147 inline Matrix<SymmetricSparse,Ref,Ref,AT>& operator&=(Matrix<SymmetricSparse,Ref,Ref,AT> &A) {
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<SymmetricSparse,Ref,Ref,AT>& operator<<=(const Matrix<Type,Row,Col,AT> &A) {
168 FMATVEC_ASSERT(A.rows() == A.cols(), AT);
169 int k_ = A.nonZeroElements();
170 if(n!=A.rows() || k!=k_) resize(A.rows(),k_,NONINIT);
171 return copy(A);
172 }
173
178 const AT* operator()() const {
179 return ele;
180 }
181
187 AT* operator()() {
188 return ele;
189 }
190
191 AT& operator()(int i, int j) {
192 FMATVEC_ASSERT(i>=0, AT);
193 FMATVEC_ASSERT(j>=0, AT);
194 FMATVEC_ASSERT(i<n, AT);
195 FMATVEC_ASSERT(j<n, AT);
196
197 return e(i,j);
198 }
199
200 const AT& operator()(int i, int j) const {
201 FMATVEC_ASSERT(i>=0, AT);
202 FMATVEC_ASSERT(j>=0, AT);
203 FMATVEC_ASSERT(i<n, AT);
204 FMATVEC_ASSERT(j<n, AT);
205
206 return e(i,j);
207 }
208
209 AT& e(int i, int j) {
210 return ele[pos(i,j)];
211 }
212
213 const AT& e(int i, int j) const {
214 return ele[pos(i,j)];
215 }
216
217 int pos(int i, int j) const {
218 for(int k=I[i]; k<I[i+1]; k++)
219 if(J[k]==j)
220 return k;
221 return 0;
222 }
223
224 int ldim() {
225 throw std::runtime_error("Matrix<SymmetricSparse, Ref, Ref, AT>::ldim() cannot be called.");
226 }
227
232 const int* Ip() const {
233 return I;
234 }
235
240 int* Ip() {
241 return I;
242 }
243
248 const int* Jp() const {
249 return J;
250 }
251
256 int* Jp() {
257 return J;
258 }
259
264 int size() const {return n;}
265
270 int rows() const {return n;}
271
276 int cols() const {return n;}
277
280 const CBLAS_ORDER blasOrder() const {
281 throw std::runtime_error("Matrix<SymmetricSparse, Ref, Ref, AT>::blasOrder() cannot be called.");
282 }
283
291 Matrix<SymmetricSparse,Ref,Ref,AT>& init(const AT &val);
292 inline Matrix<SymmetricSparse,Ref,Ref,AT>& init(Init, const AT &a=AT()) { return init(a); }
293 inline Matrix<SymmetricSparse,Ref,Ref,AT>& init(Noinit, const AT &a=AT()) { return *this; }
294
295 int nonZeroElements() const { return k; }
296 };
297 // ------------------------- Constructors -------------------------------------
298 // ----------------------------------------------------------------------------
299
300 template <class AT> template <class Type, class Row, class Col>
301 inline Matrix<SymmetricSparse,Ref,Ref,AT>& Matrix<SymmetricSparse,Ref,Ref,AT>::copy(const Matrix<Type,Row,Col,AT> &A) {
302 int k=0;
303 int i;
304 for(i=0; i<A.rows(); i++) {
305 ele[k]=A.e(i,i);
306 J[k]=i;
307 I[i]=k++;
308 for(int j=0; j<i; 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 for(int j=i+1; j<A.cols(); j++) {
316 // TODO eps
317 if(fabs(A.e(i,j))>1e-16) {
318 ele[k]=A.e(i,j);
319 J[k++]=j;
320 }
321 }
322 }
323 if(n) I[i]=k;
324 return *this;
325 }
326
327 template <class AT> template <class Row>
328 inline Matrix<SymmetricSparse,Ref,Ref,AT>& Matrix<SymmetricSparse,Ref,Ref,AT>::copy(const Matrix<Symmetric,Row,Row,AT> &A) {
329 int k=0;
330 int i;
331 for(i=0; i<A.size(); i++) {
332 ele[k]=A.ej(i,i);
333 J[k]=i;
334 I[i]=k++;
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<SymmetricSparse,Ref,Ref,AT>& Matrix<SymmetricSparse,Ref,Ref,AT>::copy(const Matrix<SymmetricSparse,Ref,Ref,AT> &A) {
349 for(int i=0; i<=n; 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<SymmetricSparse,Ref,Ref,AT>& Matrix<SymmetricSparse,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< SymmetricSparse, Ref, Ref, AT >::resize
void resize(int n, int m)
Resize a symmetric sparse matrix.
Definition: symmetric_sparse_matrix.h:109
fmatvec::Matrix
This is the basic matrix class for arbitrary matrices.
Definition: fmatvec.h:42
fmatvec::Matrix< SymmetricSparse, Ref, Ref, AT >
This is a matrix class for sparse quadratic matrices.
Definition: symmetric_sparse_matrix.h:44
fmatvec::Matrix::operator()
AT & operator()(int i, int j)
Standard constructor.
Definition: matrix.h:84
fmatvec::Matrix< SymmetricSparse, Ref, Ref, AT >::Matrix
Matrix(const Matrix< SymmetricSparse, Ref, Ref, AT > &A)
Copy Constructor.
Definition: symmetric_sparse_matrix.h:83
fmatvec::Init
Definition: types.h:89
fmatvec::Matrix< SymmetricSparse, Ref, Ref, AT >::rows
int rows() const
Number of rows.
Definition: symmetric_sparse_matrix.h:270
fmatvec::Matrix< SymmetricSparse, Ref, Ref, AT >::blasOrder
const CBLAS_ORDER blasOrder() const
Storage convention.
Definition: symmetric_sparse_matrix.h:280
fmatvec::Matrix< SymmetricSparse, Ref, Ref, AT >::cols
int cols() const
Number of columns.
Definition: symmetric_sparse_matrix.h:276
fmatvec::Matrix< SymmetricSparse, Ref, Ref, AT >::operator&=
Matrix< SymmetricSparse, Ref, Ref, AT > & operator&=(Matrix< SymmetricSparse, Ref, Ref, AT > &A)
Reference operator.
Definition: symmetric_sparse_matrix.h:147
fmatvec::Matrix< SymmetricSparse, Ref, Ref, AT >::Matrix
Matrix()
Standard constructor.
Definition: symmetric_sparse_matrix.h:73
fmatvec::Matrix< SymmetricSparse, Ref, Ref, AT >::size
int size() const
Size.
Definition: symmetric_sparse_matrix.h:264
fmatvec::Matrix< SymmetricSparse, Ref, Ref, AT >::Jp
const int * Jp() const
Pointer operator.
Definition: symmetric_sparse_matrix.h:248
fmatvec::SymmetricSparse
Shape class for sparse matrices.
Definition: types.h:160
fmatvec::Ref
Definition: types.h:98
fmatvec::Matrix< SymmetricSparse, Ref, Ref, AT >::operator=
Matrix< SymmetricSparse, Ref, Ref, AT > & operator=(const Matrix< Type, Row, Col, AT > &A)
Assignment operator.
Definition: symmetric_sparse_matrix.h:134
fmatvec::Matrix< SymmetricSparse, Ref, Ref, AT >::operator()
const AT * operator()() const
Pointer operator.
Definition: symmetric_sparse_matrix.h:178
fmatvec::Matrix< SymmetricSparse, Ref, Ref, AT >::Ip
int * Ip()
Pointer operator.
Definition: symmetric_sparse_matrix.h:240
fmatvec::Matrix::cols
int cols() const
Number of columns.
fmatvec::Matrix< SymmetricSparse, Ref, Ref, AT >::Jp
int * Jp()
Pointer operator.
Definition: symmetric_sparse_matrix.h:256
fmatvec
Namespace fmatvec.
Definition: _memory.cc:28
fmatvec::Matrix< SymmetricSparse, Ref, Ref, AT >::Ip
const int * Ip() const
Pointer operator.
Definition: symmetric_sparse_matrix.h:232
fmatvec::Matrix< SymmetricSparse, Ref, Ref, AT >::operator()
AT * operator()()
Pointer operator.
Definition: symmetric_sparse_matrix.h:187
fmatvec::Noinit
Definition: types.h:88
fmatvec::Matrix< SymmetricSparse, Ref, Ref, AT >::operator=
Matrix< SymmetricSparse, Ref, Ref, AT > & operator=(const Matrix< SymmetricSparse, Ref, Ref, AT > &A)
Assignment operator.
Definition: symmetric_sparse_matrix.h:121
fmatvec::Matrix::rows
int rows() const
Number of rows.