ug4
Loading...
Searching...
No Matches
sparsematrix_util.h
Go to the documentation of this file.
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_UTIL__
34#define __H__UG__CPU_ALGEBRA__SPARSEMATRIX_UTIL__
35
36#include "common/profiler/profiler.h"
37#include "unsorted_sparse_vector.h"
38#include "../small_algebra/small_algebra.h"
39
40#ifdef UG_PARALLEL
41#include "pcl/pcl.h"
42#include "lib_algebra/parallelization/parallel_matrix.h"
43#include "lib_algebra/parallelization/parallel_vector.h"
44#endif
45
46namespace ug
47{
48
51
53// CreateAsMultiplyOf:
54//-------------------------
77template<typename ABC_type, typename A_type, typename B_type, typename C_type>
78void CreateAsMultiplyOf(ABC_type &M, const A_type &A, const B_type &B, const C_type &C, double epsilonTruncation=0.0)
79{
80 PROFILE_FUNC_GROUP("algebra");
81 UG_ASSERT(C.num_rows() == B.num_cols() && B.num_rows() == A.num_cols(), "sizes must match");
82
83 // create output matrix M
84 M.resize_and_clear(A.num_rows(), C.num_cols());
85
86
87 typename block_multiply_traits<typename A_type::value_type, typename B_type::value_type>::ReturnType ab;
88
89 typedef UnsortedSparseVector<typename ABC_type::value_type> RowType;
90 typedef typename RowType::iterator RowIterator;
91 RowType row(C.num_cols());;
92
93 std::vector<typename ABC_type::connection> con2;
94 //typename C_type::value_type cvalue;
95
96 typename ABC_type::connection c;
97
98 typedef typename A_type::const_row_iterator cAiterator;
99 typedef typename B_type::const_row_iterator cBiterator;
100 typedef typename C_type::const_row_iterator cCiterator;
101
102 // do
103 // M_{ij} = \sum_kl A_{ik} * B_{kl} * C_{lj}
104 for(size_t i=0; i < A.num_rows(); i++)
105 {
106 row.clear();
107 cAiterator itAikEnd = A.end_row(i);
108 for(cAiterator itAik = A.begin_row(i); itAik != itAikEnd; ++itAik)
109 {
110 if(itAik.value() == 0.0) continue;
111
112 size_t k = itAik.index();
113 cBiterator itBklEnd = B.end_row(k);
114 for(cBiterator itBkl = B.begin_row(k); itBkl != itBklEnd; ++itBkl)
115 {
116 if(itBkl.value() == 0.0) continue;
117 size_t l = itBkl.index();
118 // ab = A_{ik} * B_{kl}
119 AssignMult(ab, itAik.value(), itBkl.value());
120
121 cCiterator itCljEnd = C.end_row(l);
122 for(cCiterator itClj = C.begin_row(l); itClj != itCljEnd; ++itClj)
123 {
124 if(itClj.value() == 0.0) continue;
125 AddMult( row (itClj.index() ), ab, itClj.value());
126 }
127 }
128 }
129
130 if(epsilonTruncation != 0.0)
131 {
132 double m=0;
133 for(RowIterator it = row.begin(); it != row.end(); ++it)
134 {
135 double d = BlockNorm(it->value());
136 if(d > m) m = d;
137 }
138 m *= epsilonTruncation;
139 con2.clear();
140 for(RowIterator it = row.begin(); it != row.end(); ++it)
141 if( BlockNorm(it->value()) > m )
142 con2.push_back(*it);
143 M.set_matrix_row(i, &con2[0], con2.size());
144 }
145 else
146 M.set_matrix_row(i, row.unsorted_raw_ptr(), row.num_connections());
147 }
148
149}
150
152// AddMultiplyOf:
153//-------------------------
166template<typename ABC_type, typename A_type, typename B_type, typename C_type>
167void AddMultiplyOf(ABC_type &M, const A_type &A, const B_type &B, const C_type &C, double epsilonTruncation=0.0)
168{
169 PROFILE_FUNC_GROUP("algebra");
170 UG_ASSERT(C.num_rows() == B.num_cols(), "sizes must match: nRows(C) =" << C.num_rows()<<"!="<< B.num_cols() <<"=Cols(B)");
171 UG_ASSERT(B.num_rows() == A.num_cols(), "sizes must match: nRows(B) =" << B.num_rows()<<"!="<< A.num_cols() <<"=nCols(A)");
172
173 // check
174 if(M.num_rows() != A.num_rows())
175 UG_THROW("AddMultiplyOf: row sizes mismatch: M.num_rows = "<<
176 M.num_rows()<<", A.num_rows = "<<A.num_rows());
177 if(M.num_cols() != C.num_cols())
178 UG_THROW("AddMultiplyOf: column sizes mismatch: M.num_cols = "<<
179 M.num_cols()<<", C.num_cols = "<<C.num_cols());
180
181 typename block_multiply_traits<typename A_type::value_type, typename B_type::value_type>::ReturnType ab;
182
183 typedef UnsortedSparseVector<typename ABC_type::value_type> RowType;
184 typedef typename RowType::iterator RowIterator;
185 RowType row(C.num_cols());;
186
187 std::vector<typename ABC_type::connection> con2;
188 //typename C_type::value_type cvalue;
189
190 typename ABC_type::connection c;
191
192 typedef typename A_type::const_row_iterator cAiterator;
193 typedef typename B_type::const_row_iterator cBiterator;
194 typedef typename C_type::const_row_iterator cCiterator;
195
196 // do
197 // M_{ij} = \sum_kl A_{ik} * B_{kl} * C_{lj}
198 for(size_t i=0; i < A.num_rows(); i++)
199 {
200 row.clear();
201 cAiterator itAikEnd = A.end_row(i);
202 for(cAiterator itAik = A.begin_row(i); itAik != itAikEnd; ++itAik)
203 {
204 if(itAik.value() == 0.0) continue;
205
206 size_t k = itAik.index();
207 cBiterator itBklEnd = B.end_row(k);
208 for(cBiterator itBkl = B.begin_row(k); itBkl != itBklEnd; ++itBkl)
209 {
210 if(itBkl.value() == 0.0) continue;
211 size_t l = itBkl.index();
212 // ab = A_{ik} * B_{kl}
213 AssignMult(ab, itAik.value(), itBkl.value());
214
215 cCiterator itCljEnd = C.end_row(l);
216 for(cCiterator itClj = C.begin_row(l); itClj != itCljEnd; ++itClj)
217 {
218 if(itClj.value() == 0.0) continue;
219 AddMult( row (itClj.index() ), ab, itClj.value());
220 }
221 }
222 }
223
224 if(epsilonTruncation != 0.0)
225 {
226 double m=0;
227 for(RowIterator it = row.begin(); it != row.end(); ++it)
228 {
229 double d = BlockNorm(it->value());
230 if(d > m) m = d;
231 }
232 m *= epsilonTruncation;
233 con2.clear();
234 for(RowIterator it = row.begin(); it != row.end(); ++it)
235 if( BlockNorm(it->value()) > m )
236 con2.push_back(*it);
237 M.add_matrix_row(i, &con2[0], con2.size());
238 }
239 else
240 M.add_matrix_row(i, row.unsorted_raw_ptr(), row.num_connections());
241 }
242
243}
244
245
247// CreateAsMultiplyOf:
248//-------------------------
257template<typename AB_type, typename A_type, typename B_type>
258void CreateAsMultiplyOf(AB_type &M, const A_type &A, const B_type &B)
259{
260 PROFILE_FUNC_GROUP("algebra");
261 UG_ASSERT(B.num_rows() == A.num_cols(), "sizes must match");
262
263 // create output matrix M
264 M.resize_and_clear(A.num_rows(), B.num_cols());
265
266 // types
267 typedef typename A_type::const_row_iterator cAiterator;
268 typedef typename B_type::const_row_iterator cBiterator;
269 typedef UnsortedSparseVector<typename AB_type::value_type> RowType;
270
271 RowType row(B.num_cols());
272
273 // M_{ij} = \sum_k A_{ik} * B_{kj}
274 for(size_t i=0; i < A.num_rows(); i++)
275 {
276 row.clear();
277 cAiterator itAikEnd = A.end_row(i);
278 for(cAiterator itAik = A.begin_row(i); itAik != itAikEnd; ++itAik)
279 {
280 if(itAik.value() == 0.0) continue;
281 size_t k = itAik.index();
282
283 cBiterator itBklEnd = B.end_row(k);
284 for(cBiterator itBkj = B.begin_row(k); itBkj != itBklEnd; ++itBkj)
285 {
286 if(itBkj.value() == 0.0) continue;
287 size_t j = itBkj.index();
288 AddMult( row(j), itAik.value(), itBkj.value());
289 }
290 }
291
292 M.set_matrix_row(i, row.unsorted_raw_ptr(), row.num_connections());
293 }
294
295}
296
297
299// MatAdd:
300//-------------------------
308template<typename matrix_type>
309void MatAdd(matrix_type &M, number alpha1, const matrix_type &A, number alpha2, const matrix_type &B)
310{
311 PROFILE_FUNC_GROUP("algebra");
312 UG_ASSERT(A.num_rows() == B.num_rows() && A.num_cols() == B.num_cols(), "sizes must match");
313 typedef typename matrix_type::const_row_iterator criterator;
314
315 // create output matrix M
316 if(&M != &A)
317 M.resize_and_clear(A.num_rows(), A.num_cols());
318
319 // types
320 std::vector<typename matrix_type::connection > con; con.reserve(10);
321
322 typename matrix_type::connection c;
323 for(size_t i=0; i < A.num_rows(); i++)
324 {
325 con.clear();
326 criterator itA = A.begin_row(i), endA = A.end_row(i);
327 criterator itB = B.begin_row(i), endB = B.end_row(i);
328
329 while(itA != endA && itB != endB)
330 {
331 if(itA.index() == itB.index())
332 {
333 c.dValue = alpha1 * itA.value() + alpha2 * itB.value();
334 c.iIndex = itA.index();
335 ++itA; ++itB;
336 }
337 else if (itA.index() < itB.index())
338 {
339 c.dValue = itA.value();
340 c.dValue *= alpha1;
341 c.iIndex = itA.index();
342 ++itA;
343 }
344 else
345 {
346 c.dValue = itB.value();
347 c.dValue *= alpha2;
348 c.iIndex = itB.index();
349 ++itB;
350 }
351 con.push_back(c);
352 }
353 while(itA != endA)
354 {
355 c.dValue = itA.value();
356 c.dValue *= alpha1;
357 c.iIndex = itA.index();
358 ++itA;
359 con.push_back(c);
360 }
361 while(itB != endB)
362 {
363 c.dValue = itB.value();
364 c.dValue *= alpha2;
365 c.iIndex = itB.index();
366 ++itB;
367 con.push_back(c);
368 }
369
370 M.set_matrix_row(i, &con[0], con.size());
371 }
372 M.defragment();
373}
374
382template<typename matrix_type>
383void MatAddNonDirichlet(matrix_type &M, number alpha1, const matrix_type &A, number alpha2, const matrix_type &B)
384{
385 PROFILE_FUNC_GROUP("algebra");
386 UG_ASSERT(A.num_rows() == B.num_rows() && A.num_cols() == B.num_cols(), "sizes must match");
387 typedef typename matrix_type::const_row_iterator criterator;
388
389 // create output matrix M
390 if(&M != &A)
391 M.resize_and_clear(A.num_rows(), A.num_cols());
392
393 // types
394 std::vector<typename matrix_type::connection > con; con.reserve(10);
395
396 typename matrix_type::connection c;
397 for(size_t i=0; i < A.num_rows(); i++)
398 {
399 con.clear();
400 {
401 criterator itA = A.begin_row(i), endA = A.end_row(i);
402 criterator itB = B.begin_row(i), endB = B.end_row(i);
403
404 // copy only A for dirichlet rows
405 // -> create a pattern Pii
406 typedef typename matrix_type::value_type value_type;
407 const value_type &Aii = A(i,i);
408 value_type Pii = 1.0;
409
410
411 UG_ASSERT (GetRows(Aii)==GetRows(Pii), "Huhh: Numbers of rows does not match!");
412 for(size_t alpha = 0; alpha < (size_t) GetRows(Aii); ++alpha)
413 {
414 if (IsDirichletRow(A, i, alpha)) BlockRef(Pii, alpha, alpha) = 0.0;
415 }
416
417
418 // proceed as usual
419 while(itA != endA && itB != endB)
420 {
421 // add this value: Bij:=Pii*Bij
422 value_type Bij(Pii*itB.value());
423 // UG_LOG("Bij =" << Bij << "," << "Pii=" << Pii);
424
425 if(itA.index() == itB.index())
426 {
427 c.dValue = alpha1 * itA.value() + alpha2 * Bij;
428 c.iIndex = itA.index();
429 ++itA; ++itB;
430 }
431 else if (itA.index() < itB.index())
432 {
433 c.dValue = itA.value();
434 c.dValue *= alpha1;
435 c.iIndex = itA.index();
436 ++itA;
437 }
438 else
439 {
440 c.dValue = Bij;
441 c.dValue *= alpha2;
442 c.iIndex = itB.index();
443 ++itB;
444 }
445 con.push_back(c);
446 }
447 while(itA != endA)
448 {
449 c.dValue = itA.value();
450 c.dValue *= alpha1;
451 c.iIndex = itA.index();
452 ++itA;
453 con.push_back(c);
454 }
455 while(itB != endB)
456 {
457 value_type Bij(Pii*itB.value());
458
459 c.dValue = Bij;
460 c.dValue *= alpha2;
461 c.iIndex = itB.index();
462 ++itB;
463 con.push_back(c);
464 }
465 }
466 M.set_matrix_row(i, &con[0], con.size());
467 }
468 M.defragment();
469}
470
471template<typename TSparseMatrix>
472void GetNeighborhood_worker(const TSparseMatrix &A, size_t node, size_t depth, std::vector<size_t> &indices, std::vector<bool> &bVisited)
473{
474 if(depth==0) return;
475 size_t iSizeBefore = indices.size();
476 typename TSparseMatrix::const_row_iterator itEnd = A.end_row(node);
477 for(typename TSparseMatrix::const_row_iterator it = A.begin_row(node); it != itEnd; ++it)
478 {
479 if(it.value() == 0) continue;
480 if(bVisited[it.index()] == false)
481 {
482
483 bVisited[it.index()] = true;
484 indices.push_back(it.index());
485 }
486 }
487
488 if(depth==1) return;
489 size_t iSizeAfter = indices.size();
490 for(size_t i=iSizeBefore; i<iSizeAfter; i++)
491 GetNeighborhood_worker(A, indices[i], depth-1, indices, bVisited);
492}
493
494template<typename TSparseMatrix>
495void GetNeighborhood(const TSparseMatrix &A, size_t node, size_t depth, std::vector<size_t> &indices, std::vector<bool> &bVisited, bool bResetVisitedFlags=true)
496{
497 PROFILE_FUNC_GROUP("algebra");
498 indices.clear();
499
500 if(bVisited[node] == false)
501 {
502 bVisited[node] = true;
503 indices.push_back(node);
504 }
505 GetNeighborhood_worker(A, node, depth, indices, bVisited);
506
507 if(bResetVisitedFlags)
508 for(size_t i=0; i<indices.size(); i++)
509 bVisited[indices[i]] = false;
510}
511
512template<typename TSparseMatrix>
513void GetNeighborhood(const TSparseMatrix &A, size_t node, size_t depth, std::vector<size_t> &indices)
514{
515 PROFILE_FUNC_GROUP("algebra");
516 std::vector<bool> bVisited(max(A.num_cols(), A.num_rows()), false);
517 GetNeighborhood(A, node, depth, indices, bVisited, false);
518}
519
520
521template<typename TSparseMatrix>
522void MarkNeighbors(const TSparseMatrix &A, size_t node, size_t depth, std::vector<bool> &bVisited)
523{
524 typename TSparseMatrix::const_row_iterator itEnd = A.end_row(node);
525 for(typename TSparseMatrix::const_row_iterator it = A.begin_row(node); it != itEnd; ++it)
526 {
527 if(it.value() == 0) continue;
528 bVisited[it.index()] = true;
529 MarkNeighbors(A, it.index(), depth, bVisited);
530 }
531}
532
533template<typename TSparseMatrix>
534void GetNeighborhoodHierachy_worker(const TSparseMatrix &A, size_t node, size_t depth, size_t maxdepth, std::vector< std::vector<size_t> > &indices, std::vector<bool> &bVisited)
535{
536 size_t iSizeBefore = indices[depth].size();
537 typename TSparseMatrix::const_row_iterator itEnd = A.end_row(node);
538 for(typename TSparseMatrix::const_row_iterator it = A.begin_row(node); it != itEnd; ++it)
539 {
540 if(it.value() == 0) continue;
541 if(bVisited[it.index()] == false)
542 {
543 bVisited[it.index()] = true;
544 indices[depth].push_back(it.index());
545 }
546 }
547
548 if(depth==maxdepth) return;
549 size_t iSizeAfter = indices[depth].size();
550 for(size_t i=iSizeBefore; i<iSizeAfter; i++)
551 GetNeighborhoodHierachy_worker(A, indices[i], depth+1, maxdepth, indices, bVisited);
552}
553
554template<typename TSparseMatrix>
555void GetNeighborhoodHierachy(const TSparseMatrix &A, size_t node, size_t depth, std::vector< std::vector<size_t> > &indices, std::vector<bool> &bVisited,
556 bool bResetVisitedFlags=true)
557{
558 PROFILE_FUNC_GROUP("algebra");
559 if(indices.size() != depth+1)
560 indices.resize(depth+1);
561 for(size_t i=0; i < depth+1; i++)
562 indices[i].clear();
563
564 bVisited[node] = true;
565 indices[0].push_back(node);
566
567 if(depth==0) return;
568
569 for(size_t d = 0; d < depth; d++)
570 {
571 for(size_t i=0; i<indices[d].size(); i++)
572 {
573 size_t k = indices[d][i];
574 typename TSparseMatrix::const_row_iterator itEnd = A.end_row(k);
575 for(typename TSparseMatrix::const_row_iterator it = A.begin_row(k); it != itEnd; ++it)
576 {
577 if(it.value() == 0) continue;
578 if(bVisited[it.index()] == false)
579 {
580 bVisited[it.index()] = true;
581 indices[d+1].push_back(it.index());
582 }
583 }
584
585 }
586 }
587
588 if(bResetVisitedFlags)
589 for(size_t i=0; i < depth+1; i++)
590 for(size_t j=0; i<indices[j].size(); j++)
591 bVisited[j] = false;
592}
593
594
595template<typename TSparseMatrix>
596void GetNeighborhoodHierachy(const TSparseMatrix &A, size_t node, size_t depth, std::vector< std::vector<size_t> > &indices)
597{
598 std::vector<bool> bVisited(std::max(A.num_cols(), A.num_rows()), false);
599 GetNeighborhoodHierachy(A, node, depth, indices, bVisited, false);
600}
601
602
604// GetNeighborhood:
605//-------------------------
615#if 0
616template<typename T>
617void GetNeighborhood(SparseMatrix<T> &A, size_t node, size_t depth, std::vector<size_t> &indices, int *posInConnections=NULL)
618{
619 // perhaps this is better with recursion
620 indices.clear();
621
622
623
624 vector<typename SparseMatrix<T>::const_row_iterator> iterators;
625 iterators.reserve(depth);
626
627 iterators.push_back( A.begin_row(node) );
628
629 while(iterators.size() != 0)
630 {
631 if(iterators.back().isEnd())
632 iterators.pop_back();
633 else
634 {
635 size_t index = iterators.back().index();
636 ++iterators.back();
637 if(iterators.size() < depth)
638 iterators.push_back( A.begin_row(index) );
639 else
640 {
641 size_t pos;
642 if(posInConnections == NULL)
643 {
644 for(pos=0; pos<indices.size(); pos++)
645 if(indices[pos] == index)
646 break;
647 if(pos == indices.size())
648 indices.push_back(index);
649 }
650 else
651 {
652 pos = posInConnections[index];
653 if(pos == -1)
654 {
655 pos = posInConnections[index] = indices.size();
656 indices.push_back(index);
657 }
658 }
659 // else (count etc.)
660 }
661 }
662 }
663
664 // reset posInConnections
665 if(posInConnections)
666 {
667 for(size_t i=0; i<indices.size(); i++)
668 posInConnections[indices[i]] = -1;
669 }
670
671 // sort indices
672 sort(indices.begin(), indices.end());
673}
674#endif
675
676
677
679// IsCloseToBoundary:
680//-------------------------
688template<typename TSparseMatrix>
689bool IsCloseToBoundary(const TSparseMatrix &A, size_t node, size_t distance)
690{
691 typedef typename TSparseMatrix::const_row_iterator iterator;
692 if(distance == 0) return A.is_isolated(node);
693 bool bFound = false;
694 iterator itEnd = A.end_row(node);
695 for(iterator itA = A.begin_row(node); itA != itEnd && !bFound; ++itA)
696 bFound = IsCloseToBoundary(A, itA.index(), distance-1);
697
698 return bFound;
699}
700
702
709template <typename TSparseMatrix>
710void SetRow(TSparseMatrix &A, size_t i, size_t alpha, number val = 0.0)
711{
712 typedef typename TSparseMatrix::row_iterator iterator;
713 typedef typename TSparseMatrix::value_type value_type;
714 iterator itEnd = A.end_row(i);
715 for(iterator conn = A.begin_row(i); conn != itEnd; ++conn)
716 {
717 value_type& block = conn.value();
718 // block : 1x1 (CPU=1), 3x3 (CPU=3)
719 for(size_t beta = 0; beta < (size_t) GetCols(block); ++beta)
720 {
721 BlockRef(block, alpha, beta) = val;
722 }
723 }
724}
725
727
734template <typename TSparseMatrix>
735void SetCol(TSparseMatrix &A, size_t i, size_t alpha, number val = 0.0)
736{
737 typedef typename TSparseMatrix::value_type value_type;
738 for(size_t row = 0; row != A.num_rows(); ++row)
739 {
740 value_type& block = A(row, i);
741 // block : 1x1 (CPU=1), 3x3 (CPU=3)
742 for(size_t beta = 0; beta < (size_t) GetRows(block); ++beta)
743 {
744 BlockRef(block, beta, alpha) = val;
745 }
746 }
747}
748
750// SetRow:
751//-------------------------
758template <typename TSparseMatrix>
759void SetRow(TSparseMatrix& A, size_t i, number val = 0.0)
760{
761 typedef typename TSparseMatrix::row_iterator iterator;
762 iterator itEnd = A.end_row(i);
763 for(iterator conn = A.begin_row(i); conn != itEnd; ++conn)
764 conn.value() = val;
765}
766
768// ScaleRow:
769//-------------------------
776template <typename TSparseMatrix>
777void ScaleRow(TSparseMatrix& A, size_t i, number fac)
778{
779 typedef typename TSparseMatrix::row_iterator iterator;
780 iterator itEnd = A.end_row(i);
781 for(iterator conn = A.begin_row(i); conn != itEnd; ++conn)
782 conn.value() *= fac;
783}
784
785
787// SetDirichletRow:
788//-------------------------
795template <typename TSparseMatrix>
796void SetDirichletRow(TSparseMatrix &A, size_t i, size_t alpha)
797{
798 typedef typename TSparseMatrix::row_iterator iterator;
799 typedef typename TSparseMatrix::value_type value_type;
800 BlockRef(A(i,i), alpha, alpha) = 1.0;
801
802 iterator itEnd = A.end_row(i);
803 for(iterator conn = A.begin_row(i); conn != itEnd; ++conn)
804 {
805 value_type& block = conn.value();
806 // block : 1x1 (CPU=1), 3x3 (CPU=3)
807 for(size_t beta = 0; beta < (size_t) GetCols(block); ++beta)
808 {
809 if(conn.index() != i) BlockRef(block, alpha, beta) = 0.0;
810 else if(beta != alpha) BlockRef(block, alpha, beta) = 0.0;
811 }
812 }
813}
814
816template <typename TSparseMatrix>
817bool IsDirichletRow(const TSparseMatrix &A, size_t i, size_t alpha)
818{
819 typedef typename TSparseMatrix::const_row_iterator iterator;
820 typedef typename TSparseMatrix::value_type value_type;
821
822 // no Dirichlet row,
823 if (BlockRef(A(i,i), alpha, alpha) != 1.0) return false;
824
825 // check, if row sum equals 1
826 number sum=0.0;
827 iterator itEnd = A.end_row(i);
828 for(iterator conn = A.begin_row(i); conn != itEnd; ++conn)
829 {
830 const value_type& block = conn.value();
831 // block : 1x1 (CPU=1), 3x3 (CPU=3)
832 for(size_t beta = 0; beta < (size_t) GetCols(block); ++beta)
833 {
834 sum += BlockRef(block, alpha, beta) * BlockRef(block, alpha, beta);
835 }
836 }
837 return (sum==1.0);
838}
839
840
842// SetDirichletRow:
843//-------------------------
849template <typename TSparseMatrix>
850void SetDirichletRow(TSparseMatrix& A, size_t i)
851{
852 typedef typename TSparseMatrix::row_iterator iterator;
853 typedef typename TSparseMatrix::value_type value_type;
854 iterator itEnd = A.end_row(i);
855 for(iterator conn = A.begin_row(i); conn != itEnd; ++conn)
856 {
857 value_type& block = conn.value();
858 block = 0.0;
859 }
860 A(i,i) = 1.0;
861}
862
864// SetDirichletRow:
865//-------------------------
871template <typename TSparseMatrix>
872void SetDirichletRow(TSparseMatrix& A, const std::vector<size_t> vIndex)
873{
874 typedef typename TSparseMatrix::row_iterator iterator;
875 typedef typename TSparseMatrix::value_type value_type;
876 std::vector<size_t>::const_iterator iter = vIndex.begin();
877 std::vector<size_t>::const_iterator iterEnd = vIndex.end();
878
879 for(; iter < iterEnd; ++iter)
880 {
881 const size_t i = *iter;
882 UG_ASSERT(i < A.num_rows(), "Index to large in index set.");
883
884 iterator itEnd = A.end_row(i);
885 for(iterator conn = A.begin_row(i); conn != itEnd; ++conn)
886 {
887 value_type& block = conn.value();
888 block = 0.0;
889 }
890 A(i,i) = 1.0;
891 }
892}
893
894template<typename TSparseMatrix, class TOStream>
895void SerializeMatrix(TOStream &buf, const TSparseMatrix &A)
896{
897 typedef typename TSparseMatrix::const_row_iterator iterator;
898 Serialize(buf, A.num_rows());
899 Serialize(buf, A.num_cols());
900
901 for(size_t i=0; i < A.num_rows(); i++)
902 {
903 size_t num_connections = A.num_connections(i);
904
905 // serialize number of connections
906 Serialize(buf, num_connections);
907
908 iterator itEnd = A.end_row(i);
909 for(iterator conn = A.begin_row(i); conn != itEnd; ++conn)
910 {
911 // serialize connection
912 Serialize(buf, conn.index());
913 Serialize(buf, conn.value());
914 }
915 }
916}
917
918template <typename TSparseMatrix, class TIStream>
919void DeserializeMatrix(TIStream& buf, TSparseMatrix &A)
920{
921 size_t numRows, numCols, num_connections;
922
923 Deserialize(buf, numRows);
924 Deserialize(buf, numCols);
925 A.resize_and_clear(numRows, numCols);
926
927 std::vector<typename TSparseMatrix::connection> con; con.reserve(16);
928
929 for(size_t i=0; i < A.num_rows; i++)
930 {
931 Deserialize(buf, num_connections);
932
933 con.resize(num_connections);
934
935 for(size_t j=0; j<num_connections; j++)
936 {
937 Deserialize(buf, con[j].iIndex);
938 Deserialize(buf, con[j].dValue);
939 }
940 A.set_matrix_row(i, &con[0], num_connections);
941 }
942 A.defragment();
943}
944
945template<typename TSparseMatrix>
946void ScaleSparseMatrixCommon(TSparseMatrix &A, double d)
947{
948 typedef typename TSparseMatrix::row_iterator iterator;
949 for(size_t r=0; r < A.num_rows(); r++)
950 {
951 iterator itEnd = A.end_row(r);
952 for(iterator it = A.begin_row(r); it != itEnd; ++it)
953 it.value() *= d;
954 }
955}
956
957
958template<typename TSparseMatrix, typename vector_t>
959bool AxpyCommonSparseMatrix(const TSparseMatrix &A, vector_t &dest,
960 const number &alpha1, const vector_t &v1,
961 const number &beta1, const vector_t &w1)
962{
963 typedef typename TSparseMatrix::const_row_iterator iterator;
964// UG_ASSERT(cols == x.size(), "x: " << x << " has wrong length (should be " << cols << "). A: " << *this);
965// UG_ASSERT(rows == res.size(), "res: " << x << " has wrong length (should be " << rows << "). A: " << *this);
966
967 if(alpha1 == 0.0)
968 {
969 for(size_t i=0; i < A.num_rows(); i++)
970 {
971 iterator conn = A.begin_row(i);
972 iterator itEnd = A.end_row(i);
973 if(conn == itEnd) continue;
974 MatMult(dest[i], beta1, conn.value(), w1[conn.index()]);
975 for(++conn; conn != itEnd; ++conn)
976 // res[i] += conn.value() * x[conn.index()];
977 MatMultAdd(dest[i], 1.0, dest[i], beta1, conn.value(), w1[conn.index()]);
978 }
979 }
980 else if(&dest == &v1)
981 {
982 if(alpha1 != 1.0)
983 for(size_t i=0; i < A.num_rows(); i++)
984 {
985 dest[i] *= alpha1;
986 A.mat_mult_add_row(i, dest[i], beta1, w1);
987 }
988 else
989 for(size_t i=0; i < A.num_rows(); i++)
990 A.mat_mult_add_row(i, dest[i], beta1, w1);
991
992 }
993 else
994 {
995 for(size_t i=0; i < A.num_rows(); i++)
996 {
997 VecScaleAssign(dest[i], alpha1, v1[i]);
998 A.mat_mult_add_row(i, dest[i], beta1, w1);
999 }
1000 }
1001
1002 return true;
1003}
1004
1005
1006template<typename TSparseMatrix, typename vector_t>
1007bool Axpy_transposedCommonSparseMatrix(const TSparseMatrix &A, vector_t &dest,
1008 const number &alpha1, const vector_t &v1,
1009 const number &beta1, const vector_t &w1)
1010{
1011 typedef typename TSparseMatrix::const_row_iterator iterator;
1012// UG_ASSERT(rows == x.size(), "x: " << x << " has wrong length (should be " << rows << "). A: " << *this);
1013// UG_ASSERT(cols == res.size(), "res: " << x << " has wrong length (should be " << cols << "). A: " << *this);
1014
1015 if(&dest == &v1) {
1016 if(alpha1 == 0.0)
1017 dest.set(0.0);
1018 else if(alpha1 != 1.0)
1019 dest *= alpha1;
1020 }
1021 else if(alpha1 == 0.0)
1022 dest.set(0.0);
1023 else
1024 VecScaleAssign(dest, alpha1, v1);
1025
1026 for(size_t i=0; i<A.num_rows(); i++)
1027 {
1028 iterator itEnd = A.end_row(i);
1029 for(iterator conn = A.begin_row(i); conn != itEnd; ++conn)
1030 {
1031 if(conn.value() != 0.0)
1032 // dest[conn.index()] += beta1 * conn.value() * w1[i];
1033 MatMultTransposedAdd(dest[conn.index()], 1.0, dest[conn.index()], beta1, conn.value(), w1[i]);
1034 }
1035 }
1036 return true;
1037}
1038
1039
1041template<typename TSparseMatrix>
1042size_t GetNNZs(const TSparseMatrix &A)
1043{
1044 typedef typename TSparseMatrix::const_row_iterator iterator;
1045 size_t m=0;
1046 for(size_t i=0; i<A.num_rows(); i++)
1047 {
1048 iterator itEnd = A.end_row(i);
1049 for(iterator it = A.begin_row(i); it != itEnd; ++it)
1050 if(it.value() != 0.0) m++;
1051 }
1052 return m;
1053}
1054
1056template<typename TSparseMatrix>
1057size_t GetMaxConnections(const TSparseMatrix &A)
1058{
1059 typedef typename TSparseMatrix::const_row_iterator iterator;
1060 size_t m=0;
1061 for(size_t i=0; i<A.num_rows(); i++)
1062 {
1063 //if(m < A.num_connections(i)) m = A.num_connections(i);
1064
1065 size_t n=0;
1066 iterator itEnd = A.end_row(i);
1067 for(iterator it = A.begin_row(i); it != itEnd; ++it)
1068 if(it.value() != 0.0) n++;
1069 if(m < n) m = n;
1070 }
1071
1072 return m;
1073}
1074
1075
1076template<typename TSparseMatrix>
1077bool CheckRowIterators(const TSparseMatrix &A)
1078{
1079#ifndef NDEBUG
1080 bool iIter=0;
1081
1082 for(size_t i=0; i<A.num_rows(); i++)
1083 {
1084 if (A.nrOfRowIterators[i] !=0) UG_LOG ("CheckRowIterators: Failed for row " << i << std::endl);
1085 iIter += A.nrOfRowIterators[i];
1086 }
1087 return iIter==0;
1088#else
1089 return true;
1090#endif
1091}
1092
1093template<typename TSparseMatrix>
1094bool CheckDiagonalInvertible(const TSparseMatrix &A)
1095{
1096#ifndef NDEBUG
1097 typedef typename block_traits<typename TSparseMatrix::value_type>::inverse_type inverse_type;
1098 bool bsucc=true;
1099 inverse_type inv;
1100 for(size_t i=0; i<A.num_rows(); i++)
1101 {
1102 bool b = GetInverse(inv, A(i,i));
1103 if(!b)
1104 {
1105 UG_LOG("WARNING: entry " << i << " = " << A(i,i) << " not invertible\n");
1106 bsucc = false;
1107 }
1108 }
1109 return bsucc;
1110#else
1111 return true;
1112#endif
1113}
1114
1115template<typename TVector>
1116bool CheckVectorInvertible(const TVector &v)
1117{
1118#ifndef NDEBUG
1119 typedef typename block_traits<typename TVector::value_type>::inverse_type inverse_type;
1120 bool bsucc=true;
1121 inverse_type inv;
1122
1123 for(size_t i=0; i<v.size(); i++)
1124 {
1125 bool b = GetInverse(inv, v[i]);
1126 if(!b)
1127 {
1128 UG_LOG("WARNING: entry " << i << " = " << v[i] << " not invertible\n");
1129 bsucc = false;
1130 }
1131 }
1132 return bsucc;
1133#else
1134 return true;
1135#endif
1136}
1137
1138#ifdef UG_PARALLEL
1139template<typename TSparseMatrix>
1140bool CheckDiagonalInvertible(const ParallelMatrix<TSparseMatrix> &m)
1141{
1142 return AllProcsTrue(CheckDiagonalInvertible((TSparseMatrix&)m), m.layouts()->proc_comm());
1143}
1144
1145template<typename TVector>
1146bool CheckVectorInvertible(const ParallelVector<TVector> &v)
1147{
1148 return AllProcsTrue(CheckVectorInvertible((TVector&)v), v.layouts()->proc_comm());
1149}
1150#endif
1151
1152
1153template<typename TSparseMatrix>
1158
1159template<typename TSparseMatrix>
1160typename DenseMatrixFromSparseMatrix<TSparseMatrix>::type &
1161GetDenseFromSparse(typename DenseMatrixFromSparseMatrix<TSparseMatrix>::type &A, const TSparseMatrix &S)
1162{
1163
1164 typedef typename TSparseMatrix::const_row_iterator sparse_row_iterator;
1165 size_t numRows = S.num_rows();
1166 size_t numCols = S.num_cols();
1167// PROGRESS_START(prog, n, "GetDenseFromSparse " << n);
1168 A.resize(numRows, numCols);
1169 for(size_t r=0; r<numRows; r++)
1170 for(size_t c=0; c<numCols; c++)
1171 A(r, c) = 0;
1172 for(size_t r=0; r<numRows; r++)
1173 {
1174// PROGRESS_UPDATE(prog, r);
1175 sparse_row_iterator itEnd = S.end_row(r);
1176 for(sparse_row_iterator it = S.begin_row(r); it != itEnd; ++it)
1177 A(r, it.index()) = it.value();
1178 }
1179// PROGRESS_FINISH(prog);
1180 return A;
1181}
1182
1183template<typename TSparseMatrix>
1184size_t GetDoubleSize(const TSparseMatrix &S)
1185{
1186 const size_t nrOfRows = block_traits<typename TSparseMatrix::value_type>::static_num_rows;
1187 UG_COND_THROW(nrOfRows != block_traits<typename TSparseMatrix::value_type>::static_num_cols, "only square matrices supported");
1188 return S.num_rows() * nrOfRows;
1189}
1190
1191template<typename TDoubleType, typename TSparseMatrix>
1192void GetDoubleFromSparseBlock(TDoubleType &A, const TSparseMatrix &S)
1193{
1194 const size_t nrOfRows = block_traits<typename TSparseMatrix::value_type>::static_num_rows;
1195 for(size_t r=0; r<S.num_rows(); r++)
1196 for(typename TSparseMatrix::const_row_iterator it = S.begin_row(r); it != S.end_row(r); ++it)
1197 {
1198 size_t rr = r*nrOfRows;
1199 size_t cc = it.index()*nrOfRows;
1200 for(size_t r2=0; r2<nrOfRows; r2++)
1201 for(size_t c2=0; c2<nrOfRows; c2++)
1202 A(rr + r2, cc + c2) = BlockRef(it.value(), r2, c2);
1203 }
1204}
1205
1206template<typename TDenseType, typename TSparseMatrix>
1207size_t GetDenseDoubleFromSparse(TDenseType &A, const TSparseMatrix &S)
1208{
1209 size_t N = GetDoubleSize(S);
1210 A.resize(0,0);
1211 A.resize(N, N);
1212 GetDoubleFromSparseBlock(A, S);
1213 return N;
1214}
1215
1216template<typename TDoubleSparse, typename TSparseMatrix>
1217size_t GetDoubleSparseFromBlockSparse(TDoubleSparse &A, const TSparseMatrix &S)
1218{
1219 size_t N = GetDoubleSize(S);
1220
1221 A.resize_and_clear(N, N);
1222 GetDoubleFromSparseBlock(A, S);
1223 A.defragment();
1224 return N;
1225}
1226
1227
1228
1229
1231} // end namespace ug
1232
1233
1234#endif
SparseMatrix
Definition sparsematrix_interface.h:37
SparseMatrix::begin_row
const_row_iterator begin_row(size_t row) const
ug::DenseMatrix
Definition densematrix.h:57
ug::ParallelMatrix
Wrapper for sequential matrices to handle them in parallel.
Definition parallel_matrix.h:65
ug::ParallelMatrix::layouts
ConstSmartPtr< AlgebraLayouts > layouts() const
returns the algebra layouts
Definition parallel_matrix.h:94
ug::ParallelVector
Definition parallel_vector.h:60
ug::ParallelVector::layouts
ConstSmartPtr< AlgebraLayouts > layouts() const
returns the algebra layouts
Definition parallel_vector.h:97
ug::UnsortedSparseVector
Definition unsorted_sparse_vector.h:60
ug::ScaleRow
void ScaleRow(TSparseMatrix &A, size_t i, number fac)
Definition sparsematrix_util.h:777
ug::GetDoubleFromSparseBlock
void GetDoubleFromSparseBlock(TDoubleType &A, const TSparseMatrix &S)
Definition sparsematrix_util.h:1192
ug::GetDoubleSize
size_t GetDoubleSize(const TSparseMatrix &S)
Definition sparsematrix_util.h:1184
ug::AddMultiplyOf
void AddMultiplyOf(ABC_type &M, const A_type &A, const B_type &B, const C_type &C, double epsilonTruncation=0.0)
Calculates M += A*B*C.
Definition sparsematrix_util.h:167
ug::CheckDiagonalInvertible
bool CheckDiagonalInvertible(const TSparseMatrix &A)
Definition sparsematrix_util.h:1094
ug::CheckVectorInvertible
bool CheckVectorInvertible(const TVector &v)
Definition sparsematrix_util.h:1116
ug::CheckRowIterators
bool CheckRowIterators(const TSparseMatrix &A)
Definition sparsematrix_util.h:1077
ug::GetNeighborhood
void GetNeighborhood(const TSparseMatrix &A, size_t node, size_t depth, std::vector< size_t > &indices, std::vector< bool > &bVisited, bool bResetVisitedFlags=true)
Definition sparsematrix_util.h:495
ug::GetDoubleSparseFromBlockSparse
size_t GetDoubleSparseFromBlockSparse(TDoubleSparse &A, const TSparseMatrix &S)
Definition sparsematrix_util.h:1217
ug::GetDenseFromSparse
DenseMatrixFromSparseMatrix< TSparseMatrix >::type & GetDenseFromSparse(typename DenseMatrixFromSparseMatrix< TSparseMatrix >::type &A, const TSparseMatrix &S)
Definition sparsematrix_util.h:1161
ug::AxpyCommonSparseMatrix
bool AxpyCommonSparseMatrix(const TSparseMatrix &A, vector_t &dest, const number &alpha1, const vector_t &v1, const number &beta1, const vector_t &w1)
Definition sparsematrix_util.h:959
ug::GetNeighborhood_worker
void GetNeighborhood_worker(const TSparseMatrix &A, size_t node, size_t depth, std::vector< size_t > &indices, std::vector< bool > &bVisited)
Definition sparsematrix_util.h:472
ug::MarkNeighbors
void MarkNeighbors(const TSparseMatrix &A, size_t node, size_t depth, std::vector< bool > &bVisited)
Definition sparsematrix_util.h:522
ug::GetNNZs
size_t GetNNZs(const TSparseMatrix &A)
returns the number of non-zeroes (!= number of connections)
Definition sparsematrix_util.h:1042
ug::IsDirichletRow
bool IsDirichletRow(const TSparseMatrix &A, size_t i, size_t alpha)
Evaluates 'true', iff corresponding row is Dirichlet.
Definition sparsematrix_util.h:817
ug::CreateAsMultiplyOf
void CreateAsMultiplyOf(ABC_type &M, const A_type &A, const B_type &B, const C_type &C, double epsilonTruncation=0.0)
Calculates M = A*B*C.
Definition sparsematrix_util.h:78
ug::GetDenseDoubleFromSparse
size_t GetDenseDoubleFromSparse(TDenseType &A, const TSparseMatrix &S)
Definition sparsematrix_util.h:1207
ug::Axpy_transposedCommonSparseMatrix
bool Axpy_transposedCommonSparseMatrix(const TSparseMatrix &A, vector_t &dest, const number &alpha1, const vector_t &v1, const number &beta1, const vector_t &w1)
Definition sparsematrix_util.h:1007
ug::SetDirichletRow
void SetDirichletRow(TSparseMatrix &A, size_t i, size_t alpha)
Definition sparsematrix_util.h:796
ug::IsCloseToBoundary
bool IsCloseToBoundary(const TSparseMatrix &A, size_t node, size_t distance)
gets the neighborhood of a node in the connectivity graph of a SparseMatrix.
Definition sparsematrix_util.h:689
ug::GetNeighborhoodHierachy_worker
void GetNeighborhoodHierachy_worker(const TSparseMatrix &A, size_t node, size_t depth, size_t maxdepth, std::vector< std::vector< size_t > > &indices, std::vector< bool > &bVisited)
Definition sparsematrix_util.h:534
ug::SetRow
void SetRow(TSparseMatrix &A, size_t i, size_t alpha, number val=0.0)
Definition sparsematrix_util.h:710
ug::SetCol
void SetCol(TSparseMatrix &A, size_t i, size_t alpha, number val=0.0)
Definition sparsematrix_util.h:735
ug::MatAddNonDirichlet
void MatAddNonDirichlet(matrix_type &M, number alpha1, const matrix_type &A, number alpha2, const matrix_type &B)
Calculates M = A + B.
Definition sparsematrix_util.h:383
ug::GetNeighborhoodHierachy
void GetNeighborhoodHierachy(const TSparseMatrix &A, size_t node, size_t depth, std::vector< std::vector< size_t > > &indices, std::vector< bool > &bVisited, bool bResetVisitedFlags=true)
Definition sparsematrix_util.h:555
ug::DeserializeMatrix
void DeserializeMatrix(TIStream &buf, TSparseMatrix &A)
Definition sparsematrix_util.h:919
ug::GetMaxConnections
size_t GetMaxConnections(const TSparseMatrix &A)
returns max number of non-zero connections in rows
Definition sparsematrix_util.h:1057
ug::ScaleSparseMatrixCommon
void ScaleSparseMatrixCommon(TSparseMatrix &A, double d)
Definition sparsematrix_util.h:946
ug::SerializeMatrix
void SerializeMatrix(TOStream &buf, const TSparseMatrix &A)
Definition sparsematrix_util.h:895
ug::MatAdd
void MatAdd(matrix_t &mOut, const matrix_t &m1, const matrix_t &m2)
adds two matrices and stores the result in a third one
Definition math_matrix_functions_common_impl.hpp:52
UG_ASSERT
#define UG_ASSERT(expr, msg)
Definition assert.h:70
UG_THROW
#define UG_THROW(msg)
Definition error.h:57
UG_LOG
#define UG_LOG(msg)
Definition log.h:367
UG_COND_THROW
#define UG_COND_THROW(cond, msg)
UG_COND_THROW(cond, msg) : performs a UG_THROW(msg) if cond == true.
Definition error.h:61
number
double number
Definition types.h:124
ug
the ug namespace
ug::GetCols
size_t GetCols(const T &t)
ug::GetInverse
bool GetInverse(block_traits< T >::inverse_type &inv, const T &m)
ug::AssignMult
void AssignMult(A &dest, const B &b, const C &vec)
ug::MatMultTransposedAdd
bool MatMultTransposedAdd(vector_t &dest, const number &alpha1, const vector_t &v1, const number &beta1, const matrix_t &A1, const vector_t &w1)
calculates dest = alpha1*v1 + beta1 * A1 *w1;
Definition operations_mat.h:121
ug::MatMultAdd
bool MatMultAdd(vector_t &dest, const number &alpha1, const vector_t &v1, const number &beta1, const matrix_t &A1, const vector_t &w1)
calculates dest = alpha1*v1 + beta1 * A1 *w1;
Definition operations_mat.h:68
ug::BlockRef
double & BlockRef(T &vec, size_t i)
Definition blocks.h:66
ug::Deserialize
void Deserialize(TIStream &buf, ParallelVector< T > &v)
Deerialize for ParallelVector<T>
Definition restart_bridge.cpp:112
ug::Serialize
void Serialize(TOStream &buf, const ParallelVector< T > &v)
Serialize for ParallelVector<T>
Definition restart_bridge.cpp:103
ug::GetRows
size_t GetRows(const T &t)
ug::MatMult
bool MatMult(vector_t &dest, const number &beta1, const matrix_t &A1, const vector_t &w1)
calculates dest = beta1 * A1;
Definition operations_mat.h:59
ug::AddMult
void AddMult(A &dest, const B &b, const C &vec)
ug::BlockNorm
double BlockNorm(const TYPE &v)
Definition blocks.h:57
ug::VecScaleAssign
void VecScaleAssign(double &dest, double alpha1, const double &v1)
calculates dest = alpha1*v1. for doubles
Definition operations_vec.h:49
parallel_vector.h
parallel_matrix.h
PROFILE_FUNC_GROUP
#define PROFILE_FUNC_GROUP(groups)
Definition profiler.h:258
B_type
class SparseMatrix B_type
value_type
T value_type
Definition sparsematrix_interface.h:2
num_connections
size_t num_connections(size_t row) const
ABC_type
class SparseMatrix ABC_type
Calculates M = A*B*C.
A_type
class SparseMatrix A_type
ug::DenseMatrixFromSparseMatrix
Definition sparsematrix_util.h:1155
ug::DenseMatrixFromSparseMatrix::type
DenseMatrix< VariableArray2< typename TSparseMatrix::value_type > > type
Definition sparsematrix_util.h:1156
ug::ab
Definition file_io_grdecl.cpp:60
ug::block_multiply_traits
Definition smallalgebra_interface.h:42
ug::block_traits
Definition communication_policies.h:58
unsorted_sparse_vector.h