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_IMPL__
34 : #define __H__UG__CPU_ALGEBRA__SPARSEMATRIX_IMPL__
35 :
36 : // creation etc
37 :
38 : #include <fstream>
39 : #include <cstring>
40 :
41 : #include "lib_algebra/common/operations_vec.h"
42 : #include "common/profiler/profiler.h"
43 : #include "sparsematrix.h"
44 : #include <vector>
45 : #include <algorithm>
46 :
47 :
48 :
49 : namespace ug{
50 :
51 : // defined in C99, and sometimes part of <math.h>, <cmath>
52 : template <class __T> inline bool
53 : iszero (__T __val)
54 : {
55 0 : return __val == 0;
56 : }
57 :
58 : template<typename T>
59 13 : SparseMatrix<T>::SparseMatrix()
60 : {
61 : PROFILE_SPMATRIX(SparseMatrix_constructor);
62 13 : bNeedsValues = true;
63 13 : iIterators=0;
64 13 : nnz = 0;
65 13 : m_numCols = 0;
66 13 : maxValues = 0;
67 13 : cols.resize(32);
68 13 : if(bNeedsValues) values.resize(32);
69 13 : }
70 :
71 : template<typename T>
72 : void SparseMatrix<T>::clear_and_free()
73 : {
74 : std::vector<int>().swap(rowStart);
75 : std::vector<int>().swap(rowMax);
76 : std::vector<int>().swap(rowEnd);
77 : m_numCols = 0;
78 : nnz = 0;
79 :
80 : std::vector<int>().swap(cols);
81 : std::vector<value_type>().swap(values);
82 : maxValues = 0;
83 :
84 : #ifdef CHECK_ROW_ITERATORS
85 : std::vector<int>().swap(nrOfRowIterators);
86 : #endif
87 : }
88 :
89 :
90 : template<typename T>
91 6 : void SparseMatrix<T>::resize_and_clear(size_t newRows, size_t newCols)
92 : {
93 : PROFILE_SPMATRIX(SparseMatrix_resize_and_clear);
94 6 : rowStart.clear(); rowStart.resize(newRows+1, -1);
95 6 : rowMax.clear(); rowMax.resize(newRows);
96 6 : rowEnd.clear(); rowEnd.resize(newRows, -1);
97 6 : m_numCols = newCols;
98 6 : nnz = 0;
99 :
100 6 : cols.clear(); cols.resize(newRows);
101 : values.clear();
102 6 : if(bNeedsValues) values.resize(newRows);
103 6 : maxValues = 0;
104 :
105 : #ifdef CHECK_ROW_ITERATORS
106 : nrOfRowIterators.clear();
107 : nrOfRowIterators.resize(newRows, 0);
108 : #endif
109 6 : }
110 :
111 : template<typename T>
112 0 : void SparseMatrix<T>::resize_and_keep_values(size_t newRows, size_t newCols)
113 : {
114 : PROFILE_SPMATRIX(SparseMatrix_resize_and_keep_values);
115 : //UG_LOG("SparseMatrix resize " << newRows << "x" << newCols << "\n");
116 0 : if(newRows == 0 && newCols == 0)
117 0 : return resize_and_clear(0,0);
118 :
119 0 : if(newRows != num_rows())
120 : {
121 : size_t oldrows = num_rows();
122 0 : rowStart.resize(newRows+1, -1);
123 0 : rowMax.resize(newRows);
124 0 : rowEnd.resize(newRows, -1);
125 : #ifdef CHECK_ROW_ITERATORS
126 : nrOfRowIterators.resize(newRows, 0);
127 : #endif
128 0 : for(size_t i=oldrows; i<newRows; i++)
129 : {
130 0 : rowStart[i] = -1;
131 0 : rowEnd[i] = -1;
132 : }
133 : }
134 0 : if((int)newCols < m_numCols)
135 0 : copyToNewSize(get_nnz_max_cols(newCols), newCols);
136 :
137 0 : m_numCols = newCols;
138 : }
139 :
140 :
141 : template<typename T>
142 0 : void SparseMatrix<T>::clear_retain_structure()
143 : {
144 0 : std::fill(values.begin(), values.end(), value_type(0));
145 0 : }
146 :
147 :
148 : template<typename T>
149 0 : void SparseMatrix<T>::set_as_transpose_of(const SparseMatrix<value_type> &B, double scale)
150 : {
151 : PROFILE_SPMATRIX(SparseMatrix_set_as_transpose_of);
152 0 : resize_and_clear(B.num_cols(), B.num_rows());
153 : /*rowStart.resize(B.num_cols(), 0);
154 : rowMax.resize(B.num_cols(), 0);
155 : rowEnd.resize(B.num_cols(), 0);
156 : nnz = B.nnz;
157 : bNeedsValues = B.bNeedsValues;
158 : if(bNeedsValues) values.resize(nnz);
159 : cols.resize(nnz);
160 : maxValues = B.maxValues;
161 : fragmented = 0;
162 : m_numCols = B.num_rows();
163 : size_t r, c;
164 : for(r=0; r<num_rows(); r++)
165 : for(const_row_iterator it = B.begin_row(r); it != B.end_row(r); ++it)
166 : rowMax[it.index()]++;
167 : rowEnd[0] = rowMax[0];
168 : rowEnd[0] = 0;
169 : for(c=1; c<B.num_cols(); c++)
170 : {
171 : rowStart[c] = rowMax[c-1];
172 : rowMax[c] = rowStart[c]+rowMax[c];
173 : rowEnd[c] = rowStart[c];
174 : }*/
175 :
176 0 : for(size_t r=0; r<B.num_rows(); r++)
177 : {
178 : const_row_iterator itEnd = B.end_row(r);
179 0 : for(const_row_iterator it = B.begin_row(r); it != itEnd; ++it)
180 0 : operator()(it.index(), r) = MatrixTranspose(scale*it.value());
181 : }
182 : // todo: sort rows
183 0 : }
184 :
185 : template<typename T>
186 0 : void SparseMatrix<T>::set_as_transpose_of2(const SparseMatrix<value_type> &B, double scale)
187 : {
188 : PROFILE_SPMATRIX(SparseMatrix_set_as_transpose_of2);
189 0 : resize_and_clear(B.num_cols(), B.num_rows());
190 : assert(num_cols() == B.num_rows());
191 : assert(num_rows() == B.num_cols());
192 : size_t r, c;
193 :
194 0 : for(r=0; r<num_rows(); ++r){
195 : assert(rowMax[r] == 0);
196 : }
197 :
198 0 : std::vector<int> rowsize(num_rows());
199 0 : nnz = 0;
200 0 : bNeedsValues = B.bNeedsValues; // allocate value array
201 :
202 0 : maxValues = B.maxValues;
203 0 : fragmented = 0;
204 :
205 0 : for(r=0; r<num_cols(); r++){
206 0 : for(const_row_iterator it = B.begin_row(r); it != B.end_row(r); ++it){
207 0 : if(iszero(it.value())){
208 : }else{
209 0 : ++nnz;
210 0 : ++rowMax[it.index()];
211 : }
212 : }
213 : }
214 :
215 0 : rowStart[0] = 0;
216 0 : rowEnd[0] = 0;
217 0 : for(c=1; c<num_rows(); ++c) {
218 0 : rowStart[c] = rowMax[c-1];
219 0 : rowMax[c] = rowStart[c]+rowMax[c];
220 0 : rowEnd[c] = rowStart[c];
221 : }
222 :
223 0 : cols.resize(nnz);
224 0 : if(bNeedsValues){
225 0 : values.resize(nnz);
226 : }else{
227 : }
228 :
229 0 : for(r=0; r<num_cols(); ++r){
230 0 : for(const_row_iterator it = B.begin_row(r); it != B.end_row(r); ++it){
231 0 : if(iszero(it.value())){
232 : }else{
233 0 : int idx = rowEnd[it.index()]++;
234 0 : if(bNeedsValues){
235 0 : values[idx] = scale * it.value();
236 : }else{
237 : }
238 0 : cols[idx] = r;
239 : }
240 : }
241 : }
242 :
243 : #ifndef NDEBUG
244 : for(r=0; r<num_rows(); ++r){
245 : assert(rowMax[r] == rowEnd[r]);
246 : }
247 :
248 : assert(nnz <= B.nnz);
249 :
250 : for(r=1; r<num_rows(); ++r) {
251 : assert(rowEnd[r] == rowMax[r]);
252 : }
253 : #endif
254 :
255 0 : }
256 :
257 : template<typename T>
258 : template<typename vector_t>
259 0 : inline void SparseMatrix<T>::mat_mult_add_row(size_t row, typename vector_t::value_type &dest, double alpha, const vector_t &v) const
260 : {
261 :
262 0 : size_t itEnd=rowEnd[row];
263 0 : for(size_t rowIt=rowStart[row]; rowIt != itEnd; ++rowIt)
264 0 : MatMultAdd(dest, 1.0, dest, alpha, values[rowIt], v[cols[rowIt]]);
265 :
266 : //for(const_row_iterator conn = begin_row(row); conn != end_row(row); ++conn)
267 : //MatMultAdd(dest, 1.0, dest, alpha, conn.value(), v[conn.index()]);
268 0 : }
269 :
270 :
271 : template<typename T>
272 : template<typename vector_t>
273 0 : void SparseMatrix<T>::apply_ignore_zero_rows(vector_t &dest,
274 : const number &beta1, const vector_t &w1) const
275 : {
276 0 : for(size_t i=0; i < num_rows(); i++)
277 : {
278 0 : size_t rowIt=rowStart[i];
279 0 : size_t itEnd=rowEnd[i];
280 0 : if(rowIt == itEnd)
281 0 : continue;
282 :
283 0 : MatMult(dest[i], beta1, values[rowIt], w1[cols[rowIt]]);
284 0 : for(++rowIt; rowIt != itEnd; ++rowIt)
285 : // res[i] += conn.value() * x[conn.index()];
286 0 : MatMultAdd(dest[i], 1.0, dest[i], beta1, values[rowIt], w1[cols[rowIt]]);
287 : }
288 0 : }
289 :
290 :
291 : // calculate dest = alpha1*v1 + beta1*A*w1 (A = this matrix)
292 : template<typename T>
293 : template<typename vector_t>
294 0 : void SparseMatrix<T>::axpy(vector_t &dest,
295 : const number &alpha1, const vector_t &v1,
296 : const number &beta1, const vector_t &w1) const
297 : {
298 : PROFILE_SPMATRIX(SparseMatrix_axpy);
299 : check_fragmentation();
300 0 : if(alpha1 == 0.0)
301 : {
302 0 : for(size_t i=0; i < num_rows(); i++)
303 : {
304 0 : size_t rowIt=rowStart[i];
305 0 : size_t itEnd=rowEnd[i];
306 0 : if(rowIt == itEnd)
307 : {
308 0 : dest[i] = 0.0;
309 0 : continue;
310 : }
311 0 : MatMult(dest[i], beta1, values[rowIt], w1[cols[rowIt]]);
312 0 : for(++rowIt; rowIt != itEnd; ++rowIt)
313 : // res[i] += conn.value() * x[conn.index()];
314 0 : MatMultAdd(dest[i], 1.0, dest[i], beta1, values[rowIt], w1[cols[rowIt]]);
315 : }
316 : }
317 0 : else if(&dest == &v1)
318 : {
319 0 : if(alpha1 != 1.0) {
320 0 : for(size_t i=0; i < num_rows(); i++)
321 : {
322 0 : dest[i] *= alpha1;
323 0 : mat_mult_add_row(i, dest[i], beta1, w1);
324 : }
325 : }
326 : else
327 0 : for(size_t i=0; i < num_rows(); i++)
328 0 : mat_mult_add_row(i, dest[i], beta1, w1);
329 :
330 : }
331 : else
332 : {
333 0 : for(size_t i=0; i < num_rows(); i++)
334 : {
335 0 : VecScaleAssign(dest[i], alpha1, v1[i]);
336 0 : mat_mult_add_row(i, dest[i], beta1, w1);
337 : }
338 : }
339 0 : }
340 :
341 : // calculate dest = alpha1*v1 + beta1*A^T*w1 (A = this matrix)
342 : template<typename T>
343 : template<typename vector_t>
344 0 : void SparseMatrix<T>::axpy_transposed(vector_t &dest,
345 : const number &alpha1, const vector_t &v1,
346 : const number &beta1, const vector_t &w1) const
347 : {
348 : PROFILE_SPMATRIX(SparseMatrix_axpy_transposed);
349 : check_fragmentation();
350 0 : if(&dest == &v1) {
351 0 : if(alpha1 == 0.0)
352 : dest.set(0.0);
353 0 : else if(alpha1 != 1.0)
354 : dest *= alpha1;
355 : }
356 0 : else if(alpha1 == 0.0)
357 : dest.set(0.0);
358 : else
359 : VecScaleAssign(dest, alpha1, v1);
360 :
361 0 : for(size_t i=0; i<num_rows(); i++)
362 : {
363 :
364 0 : size_t itEnd=rowEnd[i];
365 0 : for(size_t rowIt=rowStart[i]; rowIt != itEnd; ++rowIt)
366 : // dest[conn.index()] += beta1 * conn.value() * w1[i];
367 0 : if(values[rowIt] != 0.0)
368 0 : MatMultTransposedAdd(dest[cols[rowIt]], 1.0, dest[cols[rowIt]], beta1, values[rowIt], w1[i]);
369 : }
370 0 : }
371 :
372 :
373 : template<typename T>
374 : template<typename vector_t>
375 : void SparseMatrix<T>::apply_transposed_ignore_zero_rows(vector_t &dest,
376 : const number &beta1, const vector_t &w1) const
377 : {
378 : for(size_t i=0; i<num_rows(); i++)
379 : {
380 : row_iterator itEnd = end_row(i);
381 : for(const_row_iterator conn = begin_row(i); conn != itEnd; ++conn)
382 : dest[conn.index()] = 0.0;
383 : }
384 :
385 : for(size_t i=0; i<num_rows(); i++)
386 : {
387 : size_t itEnd=rowEnd[i];
388 : for(size_t rowIt=rowStart[i]; rowIt != itEnd; ++rowIt)
389 : // dest[conn.index()] += beta1 * conn.value() * w1[i];
390 : if(values[rowIt] != 0.0)
391 : MatMultTransposedAdd(dest[cols[rowIt]], 1.0, dest[cols[rowIt]], beta1, values[rowIt], w1[i]);
392 : }
393 : }
394 :
395 :
396 : template<typename T>
397 : void SparseMatrix<T>::set(double a)
398 : {
399 : PROFILE_SPMATRIX(SparseMatrix_set);
400 : check_fragmentation();
401 : for(size_t row=0; row<num_rows(); row++)
402 : {
403 : row_iterator itEnd = end_row(row);
404 : for(row_iterator it = begin_row(row); it != itEnd; ++it)
405 : {
406 : if(it.index() == row)
407 : it.value() = a;
408 : else
409 : it.value() = 0.0;
410 : }
411 : }
412 : }
413 :
414 :
415 : template<typename T>
416 0 : inline bool SparseMatrix<T>::is_isolated(size_t i) const
417 : {
418 : UG_ASSERT(i < num_rows() && i >= 0, *this << ": " << i << " out of bounds.");
419 :
420 : const_row_iterator itEnd = end_row(i);
421 0 : for(const_row_iterator it = begin_row(i); it != itEnd; ++it)
422 0 : if(it.index() != i && it.value() != 0.0)
423 : return false;
424 0 : return true;
425 : }
426 :
427 :
428 : template<typename T>
429 0 : void SparseMatrix<T>::set_matrix_row(size_t row, connection *c, size_t nr)
430 : {
431 : /*PROFILE_SPMATRIX(SparseMatrix_set_matrix_row);
432 : bool bSorted=false;
433 : for(size_t i=0; bSorted && i+1 < nr; i++)
434 : bSorted = c[i].iIndex < c[i+1].iIndex;
435 : if(bSorted)
436 : {
437 : int start;
438 : if(rowStart[row] == -1 || rowMax[row] - rowStart[row] < (int)nr)
439 : {
440 : assureValuesSize(maxValues+nr);
441 : start = maxValues;
442 : rowMax[row] = start+nr;
443 : rowStart[row] = start;
444 : maxValues+=nr;
445 : }
446 : else
447 : start = rowStart[row];
448 : rowEnd[row] = start+nr;
449 : for(size_t i=0; i<nr; i++)
450 : {
451 : cols[start+i] = c[i].iIndex;
452 : values[start+i] = c[i].dValue;
453 : }
454 : }
455 : else*/
456 : {
457 0 : for(size_t i=0; i<nr; i++)
458 0 : operator()(row, c[i].iIndex) = c[i].dValue;
459 : }
460 0 : }
461 :
462 : template<typename T>
463 0 : void SparseMatrix<T>::add_matrix_row(size_t row, connection *c, size_t nr)
464 : {
465 : //PROFILE_SPMATRIX(SparseMatrix_add_matrix_row);
466 0 : for(size_t i=0; i<nr; i++)
467 0 : operator()(row, c[i].iIndex) += c[i].dValue;
468 0 : }
469 :
470 :
471 : template<typename T>
472 0 : void SparseMatrix<T>::set_as_copy_of(const SparseMatrix<T> &B, double scale)
473 : {
474 : //PROFILE_SPMATRIX(SparseMatrix_set_as_copy_of);
475 0 : resize_and_clear(B.num_rows(), B.num_cols());
476 0 : for(size_t i=0; i < B.num_rows(); i++)
477 : {
478 : const_row_iterator itEnd = B.end_row(i);
479 0 : for(const_row_iterator it = B.begin_row(i); it != itEnd; ++it)
480 0 : operator()(i, it.index()) = scale*it.value();
481 : }
482 0 : }
483 :
484 :
485 :
486 : template<typename T>
487 0 : void SparseMatrix<T>::scale(double d)
488 : {
489 : //PROFILE_SPMATRIX(SparseMatrix_scale);
490 0 : for(size_t i=0; i < num_rows(); i++)
491 : {
492 : row_iterator itEnd = end_row(i);
493 0 : for(row_iterator it = begin_row(i); it != itEnd; ++it)
494 0 : it.value() *= d;
495 : }
496 0 : }
497 :
498 :
499 :
500 :
501 :
502 : //======================================================================================================
503 : // submatrix set/get
504 :
505 : template<typename T>
506 : template<typename M>
507 : void SparseMatrix<T>::add(const M &mat)
508 : {
509 : for(size_t i=0; i < mat.num_rows(); i++)
510 : {
511 : int r = mat.row_index(i);
512 : for(size_t j=0; j < mat.num_cols(); j++)
513 : {
514 : int c = mat.col_index(j);
515 : (*this)(r,c) += mat(i,j);
516 : }
517 : }
518 : }
519 :
520 :
521 : template<typename T>
522 : template<typename M>
523 : void SparseMatrix<T>::set(const M &mat)
524 : {
525 : for(size_t i=0; i < mat.num_rows(); i++)
526 : {
527 : int r = mat.row_index(i);
528 : for(size_t j=0; j < mat.num_cols(); j++)
529 : {
530 : int c = mat.col_index(j);
531 : (*this)(r,c) = mat(i,j);
532 : }
533 : }
534 : }
535 :
536 : template<typename T>
537 : template<typename M>
538 0 : void SparseMatrix<T>::get(M &mat) const
539 : {
540 0 : for(size_t i=0; i < mat.num_rows(); i++)
541 : {
542 0 : int r = mat.row_index(i);
543 0 : for(size_t j=0; j < mat.num_cols(); j++)
544 : {
545 0 : int c = mat.col_index(j);
546 0 : mat(i,j) = (*this)(r,c);
547 : }
548 : }
549 0 : }
550 :
551 :
552 : template<typename T>
553 51647 : int SparseMatrix<T>::get_index_internal(size_t row, int col) const
554 : {
555 : // PROFILE_SPMATRIX(SP_get_index_internal);
556 : assert(rowStart[row] != -1);
557 51647 : int l = rowStart[row];
558 51647 : int r = rowEnd[row];
559 : int mid=0;
560 161046 : while(l < r)
561 : {
562 126958 : mid = (l+r)/2;
563 126958 : if(cols[mid] < col)
564 77821 : l = mid+1;
565 49137 : else if(cols[mid] > col)
566 31578 : r = mid-1;
567 : else
568 17559 : return mid;
569 : }
570 34088 : mid = (l+r)/2;
571 : int ret;
572 34088 : if(mid < rowStart[row])
573 : ret = rowStart[row];
574 33684 : else if(mid == rowEnd[row] || col <= cols[mid])
575 : ret = mid;
576 40 : else ret = mid+1;
577 : UG_ASSERT(ret <= rowEnd[row] && ret >= rowStart[row], "row iterator row " << row << " pos " << ret << " out of bounds [" << rowStart[row] << ", " << rowEnd[row] << "]");
578 : return ret;
579 : }
580 :
581 :
582 :
583 : template<typename T>
584 22172 : int SparseMatrix<T>::get_index_const(int r, int c) const
585 : {
586 22172 : if(rowStart[r] == -1 || rowStart[r] == rowEnd[r]) return -1;
587 22172 : int index=get_index_internal(r, c);
588 22172 : if(index >= rowStart[r] && index < rowEnd[r] && cols[index] == c)
589 21291 : return index;
590 : else
591 : return -1;
592 : }
593 :
594 :
595 : template<typename T>
596 32211 : int SparseMatrix<T>::get_index(int r, int c)
597 : {
598 : // UG_LOG("get_index " << r << ", " << c << "\n");
599 : // UG_LOG(rowStart[r] << " - " << rowMax[r] << " - " << rowEnd[r] << " - " << cols.size() << " - " << maxValues << "\n");
600 32211 : if(rowStart[r] == -1 || rowStart[r] == rowEnd[r])
601 : {
602 : // UG_LOG("new row\n");
603 : // row did not start, start new row at the end of cols array
604 2884 : assureValuesSize(maxValues+1);
605 2884 : rowStart[r] = maxValues;
606 2884 : rowEnd[r] = maxValues+1;
607 2884 : rowMax[r] = maxValues+1;
608 2884 : if(bNeedsValues) values[maxValues] = 0.0;
609 2884 : cols[maxValues] = c;
610 2884 : maxValues++;
611 2884 : nnz++;
612 2884 : return maxValues-1;
613 : }
614 :
615 : /* for(int i=rowStart[r]; i<rowEnd[r]; i++)
616 : if(cols[i] == c)
617 : return i;*/
618 :
619 : // get the index where (r,c) _should_ be
620 29327 : int index=get_index_internal(r, c);
621 :
622 : if(index < rowEnd[r]
623 29327 : && cols[index] == c)
624 : {
625 : // UG_LOG("found\n");
626 : return index; // we found it
627 : }
628 :
629 : // we did not find it, so we have to add it
630 :
631 : check_row_modifiable(r);
632 :
633 : #ifndef NDEBUG
634 : assert(index == rowEnd[r] || cols[index] > c);
635 : assert(index == rowStart[r] || cols[index-1] < c);
636 : for(int i=rowStart[r]+1; i<rowEnd[r]; i++)
637 : assert(cols[i] > cols[i-1]);
638 : #endif
639 12455 : if(rowEnd[r] == rowMax[r] && rowEnd[r] == maxValues
640 21465 : && maxValues < (int)cols.size())
641 : {
642 : // UG_LOG("at end\n");
643 : // this row is stored at the end, so we can easily make more room
644 20 : rowMax[r]++;
645 20 : maxValues++;
646 : }
647 :
648 21445 : if(rowEnd[r] == rowMax[r])
649 : {
650 : // UG_LOG("renew\n");
651 : // row is full, we need to copy it to another place
652 12435 : int newSize = (rowEnd[r]-rowStart[r])*2;
653 12435 : if(maxValues+newSize > (int)cols.size())
654 : {
655 : // UG_LOG("ass\n");
656 148 : assureValuesSize(maxValues+newSize);
657 148 : index=get_index_internal(r, c);
658 : }
659 : // copy it to the end and insert index
660 12435 : fragmented += rowEnd[r]-rowStart[r];
661 12435 : index = index-rowStart[r]+maxValues;
662 12435 : int j=rowEnd[r]-rowStart[r]+maxValues;
663 12435 : if(rowEnd[r] != 0)
664 95684 : for(int i=rowEnd[r]-1; i>=rowStart[r]; i--, j--)
665 : {
666 95684 : if(j==index) j--;
667 95684 : if(bNeedsValues) values[j] = values[i];
668 95684 : cols[j] = cols[i];
669 95684 : if(i==rowStart[r]) break;
670 : }
671 12435 : rowEnd[r] = maxValues+rowEnd[r]-rowStart[r]+1;
672 12435 : rowStart[r] = maxValues;
673 12435 : rowMax[r] = maxValues+newSize;
674 12435 : maxValues += newSize;
675 : }
676 : else
677 : {
678 : // UG_LOG("enlength\n");
679 : // move part > index so we can insert the index
680 9010 : if(rowEnd[r] != 0)
681 9010 : for(int i=rowEnd[r]-1; i>=index; i--)
682 : {
683 0 : if(bNeedsValues) values[i+1] = values[i];
684 0 : cols[i+1] = cols[i];
685 0 : if(i==index) break;
686 : }
687 9010 : rowEnd[r]++;
688 : }
689 21445 : if(bNeedsValues) values[index] = 0.0;
690 21445 : cols[index] = c;
691 :
692 21445 : nnz++;
693 : #ifndef NDEBUG
694 : assert(index >= rowStart[r] && index < rowEnd[r]);
695 : for(int i=rowStart[r]+1; i<rowEnd[r]; i++)
696 : assert(cols[i] > cols[i-1]);
697 : #endif
698 21445 : return index;
699 :
700 : }
701 :
702 : template<typename T>
703 150 : void SparseMatrix<T>::copyToNewSize(size_t newSize, size_t maxCol)
704 : {
705 : PROFILE_SPMATRIX(SparseMatrix_copyToNewSize);
706 : /*UG_LOG("copyToNewSize: from " << values.size() << " to " << newSize << "\n");
707 : UG_LOG("sizes are " << cols.size() << " and " << values.size() << ", ");
708 : UG_LOG(reset_floats << "capacities are " << cols.capacity() << " and " << values.capacity() << ", NNZ = " << nnz << ", fragmentation = " <<
709 : (1-((double)nnz)/cols.size())*100.0 << "%\n");
710 : if(newSize == nnz) { UG_LOG("Defragmenting to NNZ."); }*/
711 150 : if( (iIterators > 0)
712 150 : || (newSize > values.size() && (100.0*nnz)/newSize < 20 && newSize <= cols.capacity()) )
713 : {
714 : UG_ASSERT(newSize > values.size(), "no nnz-defragmenting while using iterators.");
715 : //UG_LOG("copyToNew not defragmenting because of iterators or low fragmentation.\n");}
716 0 : cols.resize(newSize);
717 0 : cols.resize(cols.capacity());
718 0 : if(bNeedsValues) { values.resize(newSize); values.resize(cols.size()); }
719 0 : return;
720 : }
721 :
722 150 : std::vector<value_type> v(newSize);
723 150 : std::vector<int> c(newSize);
724 : size_t j=0;
725 45174 : for(size_t r=0; r<num_rows(); r++)
726 : {
727 45024 : if(rowStart[r] == -1)
728 2500 : rowStart[r] = rowEnd[r] = rowMax[r] = j;
729 : else
730 : {
731 : size_t start=j;
732 248246 : for(int k=rowStart[r]; k<rowEnd[r]; k++)
733 : {
734 205722 : if(cols[k] < (int)maxCol)
735 : {
736 205722 : if(bNeedsValues) v[j] = values[k];
737 205722 : c[j] = cols[k];
738 205722 : j++;
739 : }
740 : }
741 42524 : rowStart[r] = start;
742 42524 : rowEnd[r] = rowMax[r] = j;
743 : }
744 : }
745 150 : rowStart[num_rows()] = rowEnd[num_rows()-1];
746 150 : fragmented = 0;
747 150 : maxValues = j;
748 150 : if(bNeedsValues) values.swap(v);
749 : cols.swap(c);
750 150 : }
751 :
752 : template<typename T>
753 : void SparseMatrix<T>::check_fragmentation() const
754 : {
755 0 : if((double)nnz/(double)maxValues < 0.9)
756 : defragment();
757 : }
758 :
759 : template<typename T>
760 3032 : void SparseMatrix<T>::assureValuesSize(size_t s)
761 : {
762 3032 : if(s <= cols.size()) return;
763 150 : size_t newSize = nnz*2;
764 150 : if(newSize < s) newSize = s;
765 : copyToNewSize(newSize);
766 :
767 : }
768 :
769 : template<typename T>
770 0 : int SparseMatrix<T>::get_nnz_max_cols(size_t maxCols)
771 : {
772 : int j=0;
773 0 : for(size_t r=0; r<num_rows(); r++)
774 : {
775 0 : if(rowStart[r] == -1) continue;
776 0 : for(int k=rowStart[r]; k<rowEnd[r]; k++)
777 0 : if(cols[k] < (int)maxCols)
778 0 : j++;
779 : }
780 0 : return j;
781 : }
782 :
783 : } // namespace ug
784 :
785 : #endif
786 :
|
