Line data Source code
1 : /*
2 : * Copyright (c) 2010-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Martin Rupp
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : #ifndef __H__UG__CPU_ALGEBRA__SPARSEMATRIX__
34 : #define __H__UG__CPU_ALGEBRA__SPARSEMATRIX__
35 :
36 : #include <math.h>
37 : #include "common/common.h"
38 : #include "../algebra_common/sparsematrix_util.h"
39 : #include <iostream>
40 : #include <algorithm>
41 : #include "common/util/ostream_util.h"
42 :
43 : #include "../algebra_common/connection.h"
44 : #include "../algebra_common/matrixrow.h"
45 : #include "../common/operations_mat/operations_mat.h"
46 :
47 : #define PROFILE_SPMATRIX(name) PROFILE_BEGIN_GROUP(name, "SparseMatrix algebra")
48 :
49 : #ifndef NDEBUG
50 : #define CHECK_ROW_ITERATORS
51 : #endif
52 :
53 : namespace ug{
54 :
55 : /// \addtogroup cpu_algebra
56 : /// @{
57 :
58 :
59 : // example for the variable CRS storage structure:
60 : // say we have:
61 : // rowStart = 0 3 8
62 : // rowEnd = 3 6 11
63 : // rowMax = 3 8 11
64 : // cols ( | marking end of row): 2 5 6 | 2 6 7 x x| 8 9 10
65 :
66 : // now insert (0 3): row 0 is full (rowEnd[0]==rowMax[0]), copy it to the end, and insert index
67 : // rowStart = 11 3 8
68 : // rowEnd = 15 6 11
69 : // rowMax = 17 8 11
70 : // cols ( | marking end of row): x x x | 2 6 7 x x| 8 9 10 | 2 3 5 6 x x |
71 :
72 : // now insert (1 3): row 1 not full, we can add it
73 : // rowStart 11 3 8
74 : // rowEnd 15 7 11
75 : // rowMax = 17 8 11
76 : // cols : x x x | 2 3 6 7 x | 8 9 10 | 2 3 5 6 x x |
77 :
78 : // defragment:
79 : // rowStart 0 4 8
80 : // rowEnd 4 8 11
81 : // rowMax = 4 8 11
82 : // cols : 2 3 5 6 | 2 3 6 7 | 8 9 10
83 :
84 :
85 : /** SparseMatrix
86 : * \brief sparse matrix for big, variable sparse matrices.
87 : *
88 : * matrix is stored independent row-wise
89 : * When doing discretisation, use the add set and get methods
90 : * for dealing with submatrices of A.
91 : * For other things you can use the row iterators or
92 : * operator()-methods.
93 : *
94 : * \sa matrixrow, CreateAsMultiplyOf
95 : * \param T blocktype
96 : * \param T blocktype
97 : */
98 : template<typename TValueType> class SparseMatrix
99 : {
100 : public:
101 : typedef TValueType value_type;
102 : enum {rows_sorted=true};
103 :
104 : typedef SparseMatrix<value_type> this_type;
105 :
106 : public:
107 : typedef AlgebraicConnection<TValueType> connection;
108 : typedef MatrixRow<this_type> row_type;
109 : typedef ConstMatrixRow<this_type> const_row_type;
110 :
111 : public:
112 : // construction etc
113 : //----------------------
114 :
115 : /// constructor for empty SparseMatrix
116 : SparseMatrix();
117 : /// destructor
118 13 : virtual ~SparseMatrix () {}
119 :
120 : /**
121 : * @brief Clears the matrix vectors and frees their memory
122 : */
123 : void clear_and_free();
124 :
125 : /**
126 : * \brief resizes the SparseMatrix
127 : * \param newRows new nr of rows
128 : * \param newCols new nr of cols
129 : * \return
130 : */
131 : void resize_and_clear(size_t newRows, size_t newCols);
132 : void resize_and_keep_values(size_t newRows, size_t newCols);
133 : void clear_retain_structure();
134 :
135 : /**
136 : * \brief write in a empty SparseMatrix (this) the transpose SparseMatrix of B.
137 : * \param B the matrix of which to create the transpose of
138 : * \param scale an optional scaling
139 : * \return true on success
140 : */
141 : void set_as_transpose_of(const SparseMatrix<value_type> &B, double scale=1.0);
142 : void set_as_transpose_of2(const SparseMatrix<value_type> &B, double scale=1.0);
143 :
144 : /**
145 : * \brief create/recreate this as a copy of SparseMatrix B
146 : * \param B the matrix of which to create a copy of
147 : * \param scale an optional scaling
148 : * \return true on success
149 : */
150 : void set_as_copy_of(const SparseMatrix<value_type> &B, double scale=1.0);
151 : SparseMatrix<value_type> &operator = (const SparseMatrix<value_type> &B)
152 : {
153 0 : set_as_copy_of(B);
154 0 : return *this;
155 : }
156 :
157 :
158 : public:
159 : //! calculate dest = alpha1*v1 + beta1*A*w1 (A = this matrix)
160 : template<typename vector_t>
161 : void axpy(vector_t &dest,
162 : const number &alpha1, const vector_t &v1,
163 : const number &beta1, const vector_t &w1) const;
164 :
165 : //! calculate dest = alpha1*v1 + beta1*A^T*w1 (A = this matrix)
166 : template<typename vector_t>
167 : void axpy_transposed(vector_t &dest,
168 : const number &alpha1, const vector_t &v1,
169 : const number &beta1, const vector_t &w1) const;
170 :
171 : //! calculated dest = beta1*A*w1 . For empty rows, dest will not be changed
172 : template<typename vector_t>
173 : void apply_ignore_zero_rows(vector_t &dest,
174 : const number &beta1, const vector_t &w1) const;
175 :
176 : //! calculated dest = beta1*A*w1 . For empty cols of A (=empty rows of A^T), dest will not be changed
177 : template<typename vector_t>
178 : void apply_transposed_ignore_zero_rows(vector_t &dest,
179 : const number &beta1, const vector_t &w1) const;
180 :
181 : // DEPRECATED!
182 : //! calculate res = A x
183 : // apply is deprecated because of axpy(res, 0.0, res, 1.0, beta, w1)
184 : template<typename Vector_type>
185 : bool apply(Vector_type &res, const Vector_type &x) const
186 : {
187 0 : axpy(res, 0.0, res, 1.0, x);
188 : return true;
189 : }
190 :
191 : //! calculate res = A.T x
192 : // apply is deprecated because of axpy(res, 0.0, res, 1.0, beta, w1)
193 : template<typename Vector_type>
194 : bool apply_transposed(Vector_type &res, const Vector_type &x) const
195 : {
196 0 : axpy_transposed(res, 0.0, res, 1.0, x);
197 : return true;
198 : }
199 :
200 : // matmult_minus is deprecated because of axpy(res, 1.0, res, -1.0, x);
201 : //! calculate res -= A x
202 : template<typename Vector_type>
203 : bool matmul_minus(Vector_type &res, const Vector_type &x) const
204 : {
205 0 : axpy(res, 1.0, res, -1.0, x);
206 0 : return true;
207 : }
208 :
209 :
210 :
211 : /**
212 : * \brief check for isolated condition of an index
213 : * \param i
214 : * \return true if only A[i,i] != 0.0
215 : */
216 : inline bool is_isolated(size_t i) const;
217 :
218 : void scale(double d);
219 : SparseMatrix<value_type> &operator *= (double d) { scale(d); return *this; }
220 :
221 : // submatrix set/get functions
222 : //-------------------------------
223 :
224 : /** Add a local matrix
225 : *
226 : * The local matrix type must declare the following members:
227 : * - num_rows()
228 : * - num_cols()
229 : * - row_index(size_t i)
230 : * - col_index(size_t j)
231 : * - operator()(size_t i, size_t j)
232 : * so that mat(i,j) will go to SparseMat(mat.row_index(i), mat.col_index(j))
233 : * \param mat the whole local matrix type
234 : */
235 : template<typename M>
236 : void add(const M &mat);
237 : template<typename M>
238 : //! set local matrix \sa add
239 : void set(const M &mat);
240 : //! get local matrix \sa add
241 : template<typename M>
242 : void get(M &mat) const;
243 :
244 : // finalizing functions
245 : //----------------------
246 :
247 :
248 :
249 : inline void check_rc(size_t r, size_t c) const
250 : {
251 : UG_ASSERT(r < num_rows() && c < num_cols(), "tried to access element (" << r << ", " << c << ") of " << num_rows() << " x " << num_cols() << " matrix.");
252 : }
253 :
254 : inline void check_row_modifiable(size_t r) const
255 : {
256 : #ifdef CHECK_ROW_ITERATORS
257 : UG_ASSERT(nrOfRowIterators[r] == 0, "row " << r << " is used by an iterator and should not be modified.")
258 : #endif
259 : }
260 :
261 :
262 : //! set matrix to Id*a
263 : void set(double a);
264 :
265 : /** operator() (size_t r, size_t c) const
266 : * access connection (r, c)
267 : * \param r row
268 : * \param c column
269 : * \note if connection (r, c) is not there, returns 0.0
270 : * \return SparseMat(r, c)
271 : */
272 0 : const value_type &operator () (size_t r, size_t c) const
273 : {
274 : check_rc(r, c);
275 22172 : int j=get_index_const(r, c);
276 22172 : if(j == -1)
277 : {
278 0 : static value_type v(0.0);
279 0 : return v;
280 : }
281 : UG_ASSERT(cols[j]==(int)c && j >= rowStart[r] && j < rowEnd[r], "");
282 21291 : return values[j];
283 : }
284 :
285 : /** operator() (size_t r, size_t c) const
286 : * access or create connection (r, c)
287 : * \param r row
288 : * \param c column
289 : * \note (r,c) is added to sparsity pattern if not already there
290 : * use operator()(r,c,bConnectionFound) to prevent
291 : * \return SparseMat(r, c)=0.0 if connection created, otherwise SparseMat(r, c)
292 : */
293 : value_type &operator() (size_t r, size_t c)
294 : {
295 : check_rc(r, c);
296 32199 : int j=get_index(r, c);
297 : UG_ASSERT(j != -1 && cols[j]==(int)c && j >= rowStart[r] && j < rowEnd[r], "");
298 32199 : return values[j];
299 : }
300 :
301 : public:
302 : // row functions
303 :
304 : /**
305 : * set a row of the matrix. all previous content in this row is destroyed (@sa add_matrix_row).
306 : * \param row index of the row to set
307 : * \param c pointer to a array of sorted connections of size nr
308 : * \param nr number of connections in c
309 : * \remark will sort array c
310 : */
311 : void set_matrix_row(size_t row, connection *c, size_t nr);
312 :
313 : /**
314 : * adds the connections c to the matrixrow row.
315 : * if c has a connection con with con.iIndex=i, and the matrix already has a connection (row, i),
316 : * the function will set A(row,i) += con.dValue. otherwise the connection A(row, i) is created
317 : * and set to con.dValue.
318 : * \param row row to add to
319 : * \param c connections ("row") to be added the row.
320 : * \param nr number of connections in array c.
321 : * \return true on success.
322 : * \note you may use double connections in c.
323 : */
324 : void add_matrix_row(size_t row, connection *c, size_t nr);
325 :
326 :
327 : //! calculates dest += alpha * A[row, .] v;
328 : template<typename vector_t>
329 : inline void mat_mult_add_row(size_t row, typename vector_t::value_type &dest, double alpha, const vector_t &v) const;
330 : public:
331 : // accessor functions
332 : //----------------------
333 :
334 : //! returns number of connections of row row.
335 : inline size_t num_connections(size_t i) const
336 : {
337 0 : if(rowStart[i] == -1) return 0;
338 0 : else return rowEnd[i]-rowStart[i];
339 : }
340 :
341 : //! returns number of rows
342 : size_t num_rows() const { return rowEnd.size(); }
343 :
344 : //! returns the number of cols
345 3 : size_t num_cols() const { return m_numCols; }
346 :
347 : //! returns the total number of connections
348 0 : size_t total_num_connections() const { return nnz; }
349 :
350 : public:
351 :
352 : // Iterators
353 : //---------------------------
354 :
355 : // a row_iterator has to suppport
356 : // operator ++, operator +=, index() const, const value_type &value() const, value_type &value()
357 : // a const_row_iterator has to suppport
358 : // operator ++, operator +=, index() const, const value_type &value() const
359 :
360 :
361 : /**
362 : * row_iterator
363 : * iterator over a row
364 : */
365 : class row_iterator
366 : {
367 : SparseMatrix &A;
368 : #ifdef CHECK_ROW_ITERATORS
369 : size_t _row;
370 : #endif
371 : size_t i;
372 : public:
373 : inline void check() const {
374 : #ifdef CHECK_ROW_ITERATORS
375 : A.check_row(_row, i);
376 : #endif
377 : }
378 0 : row_iterator(SparseMatrix &_A, size_t row, size_t _i) : A(_A)
379 : #ifdef CHECK_ROW_ITERATORS
380 : , _row(row)
381 : #endif
382 0 : , i(_i) { A.add_iterator(row); }
383 : row_iterator(row_iterator &&other) : A(other.A),
384 : #ifdef CHECK_ROW_ITERATORS
385 : _row(other._row),
386 : #endif
387 : i(other.i) {
388 : #ifdef CHECK_ROW_ITERATORS
389 : A.add_iterator(other._row);
390 : #else
391 : A.add_iterator(-1);
392 : #endif
393 : }
394 : row_iterator(const row_iterator &other) : A(other.A),
395 : #ifdef CHECK_ROW_ITERATORS
396 : _row(other._row),
397 : #endif
398 : i(other.i) {
399 : #ifdef CHECK_ROW_ITERATORS
400 : A.add_iterator(other._row);
401 : #else
402 : A.add_iterator(-1);
403 : #endif
404 : }
405 : ~row_iterator() {
406 : #ifdef CHECK_ROW_ITERATORS
407 : A.remove_iterator(_row);
408 : #else
409 0 : A.remove_iterator(-1);
410 : #endif
411 0 : }
412 : row_iterator *operator ->() { return this; }
413 0 : value_type &value() { check(); return A.values[i]; }
414 0 : size_t index() const { check(); return A.cols[i]; }
415 : bool operator != (const row_iterator &o) const { return i != o.i; }
416 0 : void operator ++ () { ++i; }
417 : void operator += (int nr) { i+=nr; }
418 : bool operator == (const row_iterator &other) const { return other.i == i; check(); }
419 : };
420 : class const_row_iterator
421 : {
422 : const SparseMatrix &A;
423 : #ifdef CHECK_ROW_ITERATORS
424 : size_t _row; // int?
425 : #endif
426 : size_t i;
427 : public:
428 : inline void check() const {
429 : #ifdef CHECK_ROW_ITERATORS
430 : A.check_row(_row, i);
431 : #endif
432 : }
433 0 : const_row_iterator(const SparseMatrix &_A, size_t row, size_t _i) : A(_A),
434 : #ifdef CHECK_ROW_ITERATORS
435 : _row(row),
436 : #endif
437 0 : i(_i) {A.add_iterator(row);}
438 0 : const_row_iterator(const const_row_iterator &other) : A(other.A),
439 : #ifdef CHECK_ROW_ITERATORS
440 : _row(other._row),
441 : #endif
442 0 : i(other.i) {
443 : #ifdef CHECK_ROW_ITERATORS
444 : A.add_iterator(_row);
445 : #else
446 : A.add_iterator(-1);
447 : #endif
448 : }
449 : const_row_iterator(const_row_iterator&& other) : A(other.A),
450 : #ifdef CHECK_ROW_ITERATORS
451 : _row(other._row),
452 : #endif
453 : i(other.i) {
454 : #ifdef CHECK_ROW_ITERATORS
455 : A.add_iterator(_row);
456 : #else
457 : A.add_iterator(-1);
458 : #endif
459 : }
460 : const_row_iterator& operator=(const_row_iterator const&) = delete;
461 : const_row_iterator& operator=(const_row_iterator&&) = delete;
462 : ~const_row_iterator() {
463 : #ifdef CHECK_ROW_ITERATORS
464 : A.remove_iterator(_row);
465 : #else
466 0 : A.remove_iterator(-1);
467 : #endif
468 0 : }
469 : const const_row_iterator *operator ->() const { return this; }
470 0 : const value_type &value() const { check(); return A.values[i]; }
471 20031 : size_t index() const { check(); return A.cols[i]; }
472 0 : bool operator != (const const_row_iterator &o) const { return i != o.i; }
473 12165 : void operator ++ () { ++i; }
474 : void operator += (int nr) { i+=nr; }
475 : bool operator == (const const_row_iterator &other) const { return other.i == i; }
476 : public: // BUG
477 : // int row() const{return _row;}
478 : size_t idx() const{return i;}
479 : };
480 :
481 :
482 :
483 :
484 0 : row_iterator begin_row(size_t r) { return row_iterator(*this, r, rowStart[r]); }
485 0 : row_iterator end_row(size_t r) { return row_iterator(*this, r, rowEnd[r]); }
486 1445 : const_row_iterator begin_row(size_t r) const { return const_row_iterator(*this, r, rowStart[r]); }
487 13610 : const_row_iterator end_row(size_t r) const { return const_row_iterator(*this, r, rowEnd[r]); }
488 :
489 : row_type get_row(size_t r) { return row_type(*this, r); }
490 : const_row_type get_row(size_t r) const { return const_row_type(*this, r); }
491 :
492 : public:
493 : // connectivity functions
494 : //-------------------------
495 :
496 : bool has_connection(size_t r, size_t c) const
497 : {
498 : check_rc(r, c);
499 : bool bFound;
500 0 : get_connection(r, c, bFound);
501 0 : return bFound;
502 : }
503 :
504 : /**
505 : * \param r index of the row
506 : * \param c index of the column
507 : * \return a const_row_iterator to the connection A(r,c) if existing, otherwise end_row(row)
508 : */
509 : row_iterator get_iterator_or_next(size_t r, size_t c)
510 : {
511 : check_rc(r, c);
512 : if(rowStart[r] == -1 || rowStart[r] == rowEnd[r])
513 : return end_row(r);
514 : else
515 : {
516 : int j=get_index_internal(r, c);
517 : if(j > maxValues) return end_row(r);
518 : else return row_iterator(*this, r, j);
519 : }
520 : }
521 :
522 : /**
523 : * \param r index of the row
524 : * \param c index of the column
525 : * \return a const_row_iterator to the connection A(r,c) if existing, otherwise end_row(row)
526 : */
527 0 : const_row_iterator get_connection(size_t r, size_t c, bool &bFound) const
528 : {
529 : check_rc(r, c);
530 0 : int j=get_index_const(r, c);
531 0 : if(j != -1)
532 : {
533 0 : bFound = true;
534 0 : return const_row_iterator(*this, r, j);
535 : }
536 : else
537 : {
538 0 : bFound = false;
539 : return end_row(r);
540 : }
541 : }
542 : /**
543 : * \param r index of the row
544 : * \param c index of the column
545 : * \return a row_iterator to the connection A(r,c) if existing, otherwise end_row(row)
546 : */
547 0 : row_iterator get_connection(size_t r, size_t c, bool &bFound)
548 : {
549 : check_rc(r, c);
550 0 : int j=get_index_const(r, c);
551 0 : if(j != -1)
552 : {
553 0 : bFound = true;
554 0 : return row_iterator(*this, r, j);
555 : }
556 : else
557 : {
558 0 : bFound = false;
559 : return end_row(r);
560 : }
561 : }
562 :
563 : /**
564 : * \param r index of the row
565 : * \param c index of the column
566 : * \return a const_row_iterator to the connection A(r,c) if existing, otherwise end_row(row)
567 : */
568 : const_row_iterator get_connection(size_t r, size_t c) const
569 : {
570 : bool b;
571 0 : return get_connection(r, c, b);
572 : }
573 : /**
574 : * \param r index of the row
575 : * \param c index of the column
576 : * \return a row_iterator to the connection A(r,c)
577 : * \remark creates connection if necessary.
578 : */
579 : row_iterator get_connection(size_t r, size_t c)
580 : {
581 : check_rc(r, c);
582 : assert(bNeedsValues);
583 : int j=get_index(r, c);
584 : return row_iterator(*this, r, j);
585 : }
586 :
587 :
588 0 : void defragment()
589 : {
590 0 : if(num_rows() != 0 && num_cols() != 0)
591 0 : copyToNewSize(nnz);
592 0 : }
593 :
594 : void defragment() const
595 : {
596 0 : (const_cast<this_type*>(this))->defragment();
597 0 : }
598 :
599 : /**
600 : * copies the matrix to the standard CRS format
601 : * @param numRows (out) num rows of A
602 : * @param numCols (out) num rows of A
603 : * @param argValues (out) value_type vector with non-zero values
604 : * @param argRowStart (out) row i is from argRowStart[i] to argRowStart[i+1]
605 : * @param argColInd (out) argColInd[i] is colum index of nonzero i
606 : */
607 : void copy_crs(size_t &numRows, size_t &numCols,
608 : std::vector<value_type> &argValues, std::vector<int> &argRowStart,
609 : std::vector<int> &argColInd) const
610 : {
611 : numRows = num_rows();
612 : numCols = num_cols();
613 : defragment();
614 : argValues = values;
615 : argRowStart = rowStart;
616 : argColInd = cols;
617 : }
618 :
619 : /**
620 : * returns pointers to CRS format. note that these are only valid as long
621 : * as the matrix is not modified.
622 : * @param numRows (out) num rows of A
623 : * @param numCols (out) num rows of A
624 : * @param pValues (out) value_type vector with non-zero values
625 : * @param pRowStart (out) row i is from pRowStart[i] to pRowStart[i+1]
626 : * @param pColInd (out) pColInd[i] is colum index of nonzero i
627 : */
628 : void get_crs(size_t &numRows, size_t &numCols,
629 : value_type *&pValues, size_t *pRowStart, size_t *pColInd, size_t &nnz) const
630 : {
631 : defragment();
632 : pValues = &values[0];
633 : // FIXME:
634 : //pRowStart = &rowStart[0]; // assigning int * to size_t *
635 : //pColInd = &cols[0]; // assigning int * to size_t *
636 : UG_THROW("SparseMatrix::get_crs() needs to be fixed.");
637 : numRows = num_rows();
638 : numCols = num_cols();
639 : nnz = total_num_connections();
640 : }
641 :
642 : /**
643 : * assigns a reference to the values vector to argument vector
644 : *
645 : * \param argValues vector to be assigned
646 : */
647 : void get_values(std::vector<value_type> &argValues) const
648 : {
649 : defragment();
650 : argValues = values;
651 : }
652 :
653 :
654 : public:
655 : // output functions
656 : //----------------------
657 :
658 : void print(const char * const name = NULL) const;
659 : void printtype() const;
660 :
661 : void print_to_file(const char *filename) const;
662 : void printrow(size_t row) const;
663 :
664 : friend std::ostream &operator<<(std::ostream &out, const SparseMatrix &m)
665 : {
666 : out << "SparseMatrix " //<< m.name
667 : << " [ " << m.num_rows() << " x " << m.num_cols() << " ]";
668 : return out;
669 : }
670 :
671 :
672 0 : void p() const { print(); } // for use in gdb
673 : void pr(size_t row) const {printrow(row); } // for use in gdb
674 :
675 :
676 :
677 :
678 :
679 : private:
680 : // private functions
681 :
682 : void add_iterator(size_t row) const
683 : {
684 : #ifdef CHECK_ROW_ITERATORS
685 : nrOfRowIterators[row]++;
686 : #endif
687 1445 : iIterators++;
688 : }
689 : void remove_iterator(size_t row) const
690 : {
691 : #ifdef CHECK_ROW_ITERATORS
692 : nrOfRowIterators[row]--;
693 : UG_ASSERT(nrOfRowIterators[row] >= 0, row);
694 : #endif
695 13613 : iIterators--;
696 : UG_ASSERT(iIterators >= 0, row);
697 :
698 : }
699 : inline void check_row(size_t row, int i) const
700 : {
701 : UG_ASSERT(i < rowEnd[row] && i >= rowStart[row], "row iterator row " << row << " pos " << i << " out of bounds [" << rowStart[row] << ", " << rowEnd[row] << "]");
702 : }
703 : void assureValuesSize(size_t s);
704 : size_t get_nnz() const { return nnz; }
705 :
706 : private:
707 : // disallowed operations (not defined):
708 : //---------------------------------------
709 : SparseMatrix(SparseMatrix&); ///< disallow copy operator
710 :
711 :
712 : protected:
713 : int get_index_internal(size_t row, int col) const;
714 : int get_index_const(int r, int c) const;
715 : int get_index(int r, int c);
716 : void copyToNewSize(size_t newSize)
717 : {
718 150 : copyToNewSize(newSize, num_cols());
719 150 : }
720 : void copyToNewSize(size_t newSize, size_t maxCols);
721 : void check_fragmentation() const;
722 : int get_nnz_max_cols(size_t maxCols);
723 :
724 : public: // bug
725 : int col(size_t i) const{
726 : assert(i<cols.size());
727 : return cols[i];
728 : }
729 :
730 : #ifndef NDEBUG
731 : int iters() const{
732 : return iIterators;
733 : }
734 : #endif
735 :
736 : protected:
737 : std::vector<int> rowStart;
738 : std::vector<int> rowEnd;
739 : std::vector<int> rowMax;
740 : std::vector<int> cols;
741 : size_t fragmented;
742 : size_t nnz;
743 : bool bNeedsValues;
744 :
745 : std::vector<value_type> values;
746 : int maxValues;
747 : int m_numCols;
748 : mutable int iIterators;
749 :
750 : #ifdef CHECK_ROW_ITERATORS
751 : public:
752 : mutable std::vector<int> nrOfRowIterators;
753 : #endif
754 : };
755 :
756 :
757 : template<typename T>
758 : struct matrix_algebra_type_traits<SparseMatrix<T> >
759 : {
760 : enum{
761 : type=MATRIX_USE_ROW_FUNCTIONS
762 : };
763 : };
764 :
765 : //! calculates dest = alpha1*v1 + beta1 * A1^T *w1;
766 : template<typename vector_t, typename matrix_t>
767 : inline void MatMultTransposedAdd(vector_t &dest,
768 : const number &alpha1, const vector_t &v1,
769 : const number &beta1, const SparseMatrix<matrix_t> &A1, const vector_t &w1)
770 : {
771 : A1.axpy_transposed(dest, alpha1, v1, beta1, w1);
772 : }
773 :
774 :
775 :
776 :
777 :
778 :
779 :
780 :
781 : // end group cpu_algebra
782 : /// \}
783 :
784 : } // namespace ug
785 :
786 : //#include "matrixrow.h"
787 : #include "sparsematrix_impl.h"
788 : #include "sparsematrix_print.h"
789 :
790 : #endif
|
