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_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 :
46 : namespace ug
47 : {
48 :
49 : /// \addtogroup lib_algebra
50 : /// @{
51 :
52 : ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
53 : // CreateAsMultiplyOf:
54 : //-------------------------
55 : /**
56 : * \brief Calculates M = A*B*C.
57 : * \param M (out) Matrix M, M = A*B*C$
58 : * \param A (in) Matrix A
59 : * \param B (in) Matrix B
60 : * \param C (in) Matrix C
61 : *
62 : * Complete formula for calculating M=A*B*C:
63 : * \f[
64 : * M_{ij} = \sum_{kl} A_{ik} * B_{kl} * C_{lj}
65 : * \f]
66 : * Calculation is done on row-basis without
67 : * a temporary BC or AB matrix. This has shown to be much faster than
68 : * implementations with temporary matrices due to cache effects.
69 : * We also added an improved way of storing the results of the calculation:
70 : * when we go through connections of B and C and want to add the connection
71 : * (i, j) to M, we need to know if this connection already exists. For this we have
72 : * an array posInConnections, needs n=A.num_rows() memory.
73 : * posInConnections[i]: index in the connections for current row (if not in row: -1)
74 : * tried this also with std::map, but took 1511.53 ms instead of 393.972 ms
75 : * searching in the connections array is also slower
76 : */
77 : template<typename ABC_type, typename A_type, typename B_type, typename C_type>
78 : void 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 :
151 : ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
152 : // AddMultiplyOf:
153 : //-------------------------
154 : /**
155 : * \brief Calculates M += A*B*C.
156 : * \param M (in/out) Matrix M, M += A*B*C$
157 : * \param A (in) Matrix A
158 : * \param B (in) Matrix B
159 : * \param C (in) Matrix C
160 : *
161 : * Complete formula for calculating M=A*B*C:
162 : * \f[
163 : * M_{ij} += \sum_{kl} A_{ik} * B_{kl} * C_{lj}
164 : * \f]
165 : */
166 : template<typename ABC_type, typename A_type, typename B_type, typename C_type>
167 0 : void 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 0 : if(M.num_rows() != A.num_rows())
175 0 : UG_THROW("AddMultiplyOf: row sizes mismatch: M.num_rows = "<<
176 : M.num_rows()<<", A.num_rows = "<<A.num_rows());
177 0 : if(M.num_cols() != C.num_cols())
178 0 : 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 0 : 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 0 : for(size_t i=0; i < A.num_rows(); i++)
199 : {
200 0 : row.clear();
201 : cAiterator itAikEnd = A.end_row(i);
202 0 : for(cAiterator itAik = A.begin_row(i); itAik != itAikEnd; ++itAik)
203 : {
204 0 : if(itAik.value() == 0.0) continue;
205 :
206 : size_t k = itAik.index();
207 : cBiterator itBklEnd = B.end_row(k);
208 0 : for(cBiterator itBkl = B.begin_row(k); itBkl != itBklEnd; ++itBkl)
209 : {
210 0 : if(itBkl.value() == 0.0) continue;
211 : size_t l = itBkl.index();
212 : // ab = A_{ik} * B_{kl}
213 0 : AssignMult(ab, itAik.value(), itBkl.value());
214 :
215 : cCiterator itCljEnd = C.end_row(l);
216 0 : for(cCiterator itClj = C.begin_row(l); itClj != itCljEnd; ++itClj)
217 : {
218 0 : if(itClj.value() == 0.0) continue;
219 0 : AddMult( row (itClj.index() ), ab, itClj.value());
220 : }
221 : }
222 : }
223 :
224 0 : if(epsilonTruncation != 0.0)
225 : {
226 : double m=0;
227 0 : for(RowIterator it = row.begin(); it != row.end(); ++it)
228 : {
229 0 : double d = BlockNorm(it->value());
230 0 : if(d > m) m = d;
231 : }
232 0 : m *= epsilonTruncation;
233 : con2.clear();
234 0 : for(RowIterator it = row.begin(); it != row.end(); ++it)
235 0 : if( BlockNorm(it->value()) > m )
236 0 : con2.push_back(*it);
237 0 : M.add_matrix_row(i, &con2[0], con2.size());
238 : }
239 : else
240 0 : M.add_matrix_row(i, row.unsorted_raw_ptr(), row.num_connections());
241 : }
242 :
243 0 : }
244 :
245 :
246 : ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
247 : // CreateAsMultiplyOf:
248 : //-------------------------
249 : /**
250 : * \brief Calculates M = A*B.
251 : * \param M (out) Matrix M, M = A*B$
252 : * \param A (in) Matrix A
253 : * \param B (in) Matrix B
254 : * \f$ M_{ij} = \sum_k A_{ik} * B_{kj} \f$
255 : * For implementation details, see also CreateAsMultiplyOf(M, A, B, C).
256 : */
257 : template<typename AB_type, typename A_type, typename B_type>
258 : void 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 :
298 : ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
299 : // MatAdd:
300 : //-------------------------
301 : /**
302 : * \brief Calculates M = A + B
303 : * \param M (out) Matrix M, M = A + B
304 : * \param A (in) Matrix A
305 : * \param B (in) Matrix B
306 : * note: A and/or B may be equal to M.
307 : */
308 : template<typename matrix_type>
309 0 : void 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 0 : if(&M != &A)
317 0 : M.resize_and_clear(A.num_rows(), A.num_cols());
318 :
319 : // types
320 0 : std::vector<typename matrix_type::connection > con; con.reserve(10);
321 :
322 : typename matrix_type::connection c;
323 0 : 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 0 : while(itA != endA && itB != endB)
330 : {
331 0 : if(itA.index() == itB.index())
332 : {
333 0 : c.dValue = alpha1 * itA.value() + alpha2 * itB.value();
334 0 : c.iIndex = itA.index();
335 : ++itA; ++itB;
336 : }
337 0 : else if (itA.index() < itB.index())
338 : {
339 0 : c.dValue = itA.value();
340 0 : c.dValue *= alpha1;
341 0 : c.iIndex = itA.index();
342 : ++itA;
343 : }
344 : else
345 : {
346 0 : c.dValue = itB.value();
347 0 : c.dValue *= alpha2;
348 0 : c.iIndex = itB.index();
349 : ++itB;
350 : }
351 0 : con.push_back(c);
352 : }
353 0 : while(itA != endA)
354 : {
355 0 : c.dValue = itA.value();
356 0 : c.dValue *= alpha1;
357 0 : c.iIndex = itA.index();
358 : ++itA;
359 0 : con.push_back(c);
360 : }
361 0 : while(itB != endB)
362 : {
363 0 : c.dValue = itB.value();
364 0 : c.dValue *= alpha2;
365 0 : c.iIndex = itB.index();
366 : ++itB;
367 0 : con.push_back(c);
368 : }
369 :
370 0 : M.set_matrix_row(i, &con[0], con.size());
371 : }
372 0 : M.defragment();
373 0 : }
374 :
375 : /**
376 : * \brief Calculates M = A + B
377 : * \param M (out) Matrix M, M = A + B
378 : * \param A (in) Matrix A
379 : * \param B (in) Matrix B
380 : * note: A and/or B may be equal to M.
381 : */
382 : template<typename matrix_type>
383 0 : void 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 0 : if(&M != &A)
391 0 : M.resize_and_clear(A.num_rows(), A.num_cols());
392 :
393 : // types
394 0 : std::vector<typename matrix_type::connection > con; con.reserve(10);
395 :
396 : typename matrix_type::connection c;
397 0 : 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 0 : 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 0 : for(size_t alpha = 0; alpha < (size_t) GetRows(Aii); ++alpha)
413 : {
414 0 : if (IsDirichletRow(A, i, alpha)) BlockRef(Pii, alpha, alpha) = 0.0;
415 : }
416 :
417 :
418 : // proceed as usual
419 0 : while(itA != endA && itB != endB)
420 : {
421 : // add this value: Bij:=Pii*Bij
422 0 : value_type Bij(Pii*itB.value());
423 : // UG_LOG("Bij =" << Bij << "," << "Pii=" << Pii);
424 :
425 0 : if(itA.index() == itB.index())
426 : {
427 0 : c.dValue = alpha1 * itA.value() + alpha2 * Bij;
428 0 : c.iIndex = itA.index();
429 : ++itA; ++itB;
430 : }
431 0 : else if (itA.index() < itB.index())
432 : {
433 0 : c.dValue = itA.value();
434 0 : c.dValue *= alpha1;
435 0 : c.iIndex = itA.index();
436 : ++itA;
437 : }
438 : else
439 : {
440 : c.dValue = Bij;
441 0 : c.dValue *= alpha2;
442 0 : c.iIndex = itB.index();
443 : ++itB;
444 : }
445 0 : con.push_back(c);
446 : }
447 0 : while(itA != endA)
448 : {
449 0 : c.dValue = itA.value();
450 0 : c.dValue *= alpha1;
451 0 : c.iIndex = itA.index();
452 : ++itA;
453 0 : con.push_back(c);
454 : }
455 0 : while(itB != endB)
456 : {
457 0 : value_type Bij(Pii*itB.value());
458 :
459 : c.dValue = Bij;
460 0 : c.dValue *= alpha2;
461 0 : c.iIndex = itB.index();
462 : ++itB;
463 0 : con.push_back(c);
464 : }
465 : }
466 0 : M.set_matrix_row(i, &con[0], con.size());
467 : }
468 0 : M.defragment();
469 0 : }
470 :
471 : template<typename TSparseMatrix>
472 0 : void GetNeighborhood_worker(const TSparseMatrix &A, size_t node, size_t depth, std::vector<size_t> &indices, std::vector<bool> &bVisited)
473 : {
474 0 : if(depth==0) return;
475 : size_t iSizeBefore = indices.size();
476 : typename TSparseMatrix::const_row_iterator itEnd = A.end_row(node);
477 0 : for(typename TSparseMatrix::const_row_iterator it = A.begin_row(node); it != itEnd; ++it)
478 : {
479 0 : if(it.value() == 0) continue;
480 0 : if(bVisited[it.index()] == false)
481 : {
482 :
483 : bVisited[it.index()] = true;
484 0 : indices.push_back(it.index());
485 : }
486 : }
487 :
488 0 : if(depth==1) return;
489 : size_t iSizeAfter = indices.size();
490 0 : for(size_t i=iSizeBefore; i<iSizeAfter; i++)
491 0 : GetNeighborhood_worker(A, indices[i], depth-1, indices, bVisited);
492 : }
493 :
494 : template<typename TSparseMatrix>
495 0 : void 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 0 : if(bVisited[node] == false)
501 : {
502 : bVisited[node] = true;
503 0 : indices.push_back(node);
504 : }
505 0 : GetNeighborhood_worker(A, node, depth, indices, bVisited);
506 :
507 0 : if(bResetVisitedFlags)
508 0 : for(size_t i=0; i<indices.size(); i++)
509 0 : bVisited[indices[i]] = false;
510 0 : }
511 :
512 : template<typename TSparseMatrix>
513 : void 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 :
521 : template<typename TSparseMatrix>
522 : void 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 :
533 : template<typename TSparseMatrix>
534 : 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)
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 :
554 : template<typename TSparseMatrix>
555 : void 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 :
595 : template<typename TSparseMatrix>
596 : void 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 :
603 : ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
604 : // GetNeighborhood:
605 : //-------------------------
606 : /**
607 : * \brief gets the neighborhood of a node in the connectivity graph of a SparseMatrix.
608 : * \param A (in) Matrix A
609 : * \param node (in) the node where to start
610 : * \param depth (in) the depth of neighborhood. 0 = empty.
611 : * \param indices (out) the indices of the neighbors
612 : * \param posInConnections array to speed up computation. Has to be posInConnections[i] = 0 for all i=0..A.num_rows(). Can be NULL.
613 : * \note the node itself is only included if there is a connection from node to node.
614 : */
615 : #if 0
616 : template<typename T>
617 : void 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 :
678 : ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
679 : // IsCloseToBoundary:
680 : //-------------------------
681 : /**
682 : * \brief determines if a node is close to a unconnected node in the connectivity graph of a SparseMatrix.
683 : * \param A (in) Matrix A
684 : * \param node (in) the node where to start
685 : * \param distance (in) up to which distance "close" is.
686 : * \return if there is a distance long path in graph(A) to an unconnected node, true. otherwise false.
687 : */
688 : template<typename TSparseMatrix>
689 : bool 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 :
701 : ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
702 : /**
703 : * set value for row for entry (i,alpha).
704 : * \param A (in) Matrix A
705 : * \param i (in) row to set
706 : * \param alpha the alpha index
707 : * \param val the value to be set
708 : */
709 : template <typename TSparseMatrix>
710 0 : void 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 0 : 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 0 : for(size_t beta = 0; beta < (size_t) GetCols(block); ++beta)
720 : {
721 0 : BlockRef(block, alpha, beta) = val;
722 : }
723 : }
724 0 : }
725 :
726 : ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
727 : /**
728 : * set value for col for entry (i,alpha).
729 : * \param A (in) Matrix A
730 : * \param i (in) col to set
731 : * \param alpha the alpha index
732 : * \param val the value to be set
733 : */
734 : template <typename TSparseMatrix>
735 0 : void SetCol(TSparseMatrix &A, size_t i, size_t alpha, number val = 0.0)
736 : {
737 : typedef typename TSparseMatrix::value_type value_type;
738 0 : 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 0 : for(size_t beta = 0; beta < (size_t) GetRows(block); ++beta)
743 : {
744 0 : BlockRef(block, beta, alpha) = val;
745 : }
746 : }
747 0 : }
748 :
749 : ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
750 : // SetRow:
751 : //-------------------------
752 : /**
753 : * set value for (block-)row i.
754 : * \param A (in) Matrix A
755 : * \param i (in) row to scales
756 : * \param val (in) value to be set
757 : */
758 : template <typename TSparseMatrix>
759 0 : void SetRow(TSparseMatrix& A, size_t i, number val = 0.0)
760 : {
761 : typedef typename TSparseMatrix::row_iterator iterator;
762 0 : iterator itEnd = A.end_row(i);
763 0 : for(iterator conn = A.begin_row(i); conn != itEnd; ++conn)
764 0 : conn.value() = val;
765 0 : }
766 :
767 : ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
768 : // ScaleRow:
769 : //-------------------------
770 : /**
771 : * scales (block-)row i.
772 : * \param A (in) Matrix A
773 : * \param i (in) row to scales
774 : * \param fac (in) Scaling factor
775 : */
776 : template <typename TSparseMatrix>
777 : void 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 :
786 : ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
787 : // SetDirichletRow:
788 : //-------------------------
789 : /**
790 : * set Dirichlet row for entry (i,alpha).
791 : * \param A (in) Matrix A
792 : * \param i (in) row to set dirichlet, that is A(i,i)(alpha, alpha) = 1.0, A(i,k)(alpha, beta) = 0.0 for all (k, beta) != (i, alpha)$.
793 : * \param alpha the alpha index
794 : */
795 : template <typename TSparseMatrix>
796 0 : void 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 0 : BlockRef(A(i,i), alpha, alpha) = 1.0;
801 :
802 : iterator itEnd = A.end_row(i);
803 0 : 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 0 : for(size_t beta = 0; beta < (size_t) GetCols(block); ++beta)
808 : {
809 0 : if(conn.index() != i) BlockRef(block, alpha, beta) = 0.0;
810 0 : else if(beta != alpha) BlockRef(block, alpha, beta) = 0.0;
811 : }
812 : }
813 0 : }
814 :
815 : //! Evaluates 'true', iff corresponding row is Dirichlet
816 : template <typename TSparseMatrix>
817 0 : bool 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 0 : 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 0 : 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 0 : for(size_t beta = 0; beta < (size_t) GetCols(block); ++beta)
833 : {
834 0 : sum += BlockRef(block, alpha, beta) * BlockRef(block, alpha, beta);
835 : }
836 : }
837 0 : return (sum==1.0);
838 : }
839 :
840 :
841 : ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
842 : // SetDirichletRow:
843 : //-------------------------
844 : /**
845 : * set Dirichlet row for block i.
846 : * \param A (in) Matrix A
847 : * \param i (in) row to set dirichlet, that is A(i,i) = 1.0, A(i,k) = 0.0 for all k != i.
848 : */
849 : template <typename TSparseMatrix>
850 : void 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 :
863 : ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
864 : // SetDirichletRow:
865 : //-------------------------
866 : /**
867 : * set Dirichlet row for block i.
868 : * \param[in/out] Matrix A
869 : * \param[in] vIndex vector of row indices to set dirichlet, that is A(i,i) = 1.0, A(i,k) = 0.0 for all k != i.
870 : */
871 : template <typename TSparseMatrix>
872 : void 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 :
894 : template<typename TSparseMatrix, class TOStream>
895 : void 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 :
918 : template <typename TSparseMatrix, class TIStream>
919 : void 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 :
945 : template<typename TSparseMatrix>
946 : void 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 :
958 : template<typename TSparseMatrix, typename vector_t>
959 : bool 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 :
1006 : template<typename TSparseMatrix, typename vector_t>
1007 : bool 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 :
1040 : /// returns the number of non-zeroes (!= number of connections)
1041 : template<typename TSparseMatrix>
1042 : size_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 :
1055 : /// returns max number of non-zero connections in rows
1056 : template<typename TSparseMatrix>
1057 : size_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 :
1076 : template<typename TSparseMatrix>
1077 : bool 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 :
1093 : template<typename TSparseMatrix>
1094 : bool 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 :
1115 : template<typename TVector>
1116 : bool 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
1139 : template<typename TSparseMatrix>
1140 : bool CheckDiagonalInvertible(const ParallelMatrix<TSparseMatrix> &m)
1141 : {
1142 : return AllProcsTrue(CheckDiagonalInvertible((TSparseMatrix&)m), m.layouts()->proc_comm());
1143 : }
1144 :
1145 : template<typename TVector>
1146 : bool CheckVectorInvertible(const ParallelVector<TVector> &v)
1147 : {
1148 : return AllProcsTrue(CheckVectorInvertible((TVector&)v), v.layouts()->proc_comm());
1149 : }
1150 : #endif
1151 :
1152 :
1153 : template<typename TSparseMatrix>
1154 : struct DenseMatrixFromSparseMatrix
1155 : {
1156 : typedef DenseMatrix<VariableArray2<typename TSparseMatrix::value_type> > type;
1157 : };
1158 :
1159 : template<typename TSparseMatrix>
1160 : typename DenseMatrixFromSparseMatrix<TSparseMatrix>::type &
1161 : GetDenseFromSparse(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 :
1183 : template<typename TSparseMatrix>
1184 : size_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 0 : return S.num_rows() * nrOfRows;
1189 : }
1190 :
1191 : template<typename TDoubleType, typename TSparseMatrix>
1192 0 : void GetDoubleFromSparseBlock(TDoubleType &A, const TSparseMatrix &S)
1193 : {
1194 : const size_t nrOfRows = block_traits<typename TSparseMatrix::value_type>::static_num_rows;
1195 0 : for(size_t r=0; r<S.num_rows(); r++)
1196 0 : for(typename TSparseMatrix::const_row_iterator it = S.begin_row(r); it != S.end_row(r); ++it)
1197 : {
1198 0 : size_t rr = r*nrOfRows;
1199 0 : size_t cc = it.index()*nrOfRows;
1200 0 : for(size_t r2=0; r2<nrOfRows; r2++)
1201 0 : for(size_t c2=0; c2<nrOfRows; c2++)
1202 0 : A(rr + r2, cc + c2) = BlockRef(it.value(), r2, c2);
1203 : }
1204 0 : }
1205 :
1206 : template<typename TDenseType, typename TSparseMatrix>
1207 0 : size_t GetDenseDoubleFromSparse(TDenseType &A, const TSparseMatrix &S)
1208 : {
1209 : size_t N = GetDoubleSize(S);
1210 0 : A.resize(0,0);
1211 0 : A.resize(N, N);
1212 0 : GetDoubleFromSparseBlock(A, S);
1213 0 : return N;
1214 : }
1215 :
1216 : template<typename TDoubleSparse, typename TSparseMatrix>
1217 0 : size_t GetDoubleSparseFromBlockSparse(TDoubleSparse &A, const TSparseMatrix &S)
1218 : {
1219 : size_t N = GetDoubleSize(S);
1220 :
1221 0 : A.resize_and_clear(N, N);
1222 0 : GetDoubleFromSparseBlock(A, S);
1223 0 : A.defragment();
1224 0 : return N;
1225 : }
1226 :
1227 :
1228 :
1229 :
1230 : /// @}
1231 : } // end namespace ug
1232 :
1233 :
1234 : #endif
|
