map_matrix.h
Go to the documentation of this file.
1// LIC// ====================================================================
2// LIC// This file forms part of oomph-lib, the object-oriented,
3// LIC// multi-physics finite-element library, available
4// LIC// at http://www.oomph-lib.org.
5// LIC//
6// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26#ifndef OOMPH_MAP_MATRIX_HEADER
27#define OOMPH_MAP_MATRIX_HEADER
28
29
30// Config header generated by autoconfig
31#ifdef HAVE_CONFIG_H
32#include <oomph-lib-config.h>
33#endif
34
35
36#ifdef OOMPH_HAS_MPI
37#include "mpi.h"
38#endif
39
40#include <map>
41
42// oomph-lib headers
43#include "Vector.h"
44#include "oomph_utilities.h"
45
46namespace oomph
47{
48 //================================================================
49 /// MapMatrixMixed is a generalised, STL-map-based, sparse(ish) matrix
50 /// class with mixed indices
51 ///
52 /// The matrix is indexed by indices of type KEY_TYPE_ROW and KEY_TYPE_COL
53 /// and has entries of type VALUE_TYPE.
54 ///
55 /// Careful: If a zero entry is referenced then it is created in
56 /// memory. Therefore this isn't really a practical sparse matrix scheme.
57 /// Do not loop over `all' possible indices as even looking
58 /// at them will inflate the matrix until it occupies as much
59 /// space as a full one -- use (modification of) output routine
60 /// to retrieve all nonzero entries.
61 ///
62 /// However, this is not a serious restriction, as the main purpose
63 /// of this class is to allow non-integer indices.
64 ///
65 /// Example of usage:
66 /// \code
67 ///
68 ///
69 /// // Assume we have a Vector of pointers to objects:
70 /// Vector<Rubbish*> object_pt;
71 ///
72 /// [...]
73 ///
74 /// // Number of objects
75 /// int nentry=object_pt.size();
76 ///
77 /// // Use the pointers to the objects and associated integers as indices
78 /// // in a MapMatrixMixed whose entries are of type int
79 /// MapMatrixMixed<Rubbish*,int,int> like_a_matrix;
80 ///
81 /// for (int i=1;i<nentry;i++)
82 /// {
83 /// for (int j=1;j<nentry;j++)
84 /// {
85 /// int number=100*i+j;
86 /// like_a_matrix(object_pt[i],j)=number;
87 /// }
88 /// }
89 ///
90 /// oomph_info << "Matrix has nnz() " << like_a_matrix.nnz() <<
91 /// " and size() " << like_a_matrix.size() << std::endl;
92 ///
93 /// oomph_info << "\n\n\n Here are the nonzero entries: i, j, a(i,j)\n";
94 ///
95 /// like_a_matrix.output(oomph_info);
96 ///
97 /// // Can be used like a normal matrix:
98 /// like_a_matrix(object_pt[1],20)+=13;
99 /// like_a_matrix(object_pt[1],1)+=13;
100 ///
101 /// oomph_info << "\n\n\n Here are the nonzero entries: i, j, a(i,j)\n";
102 /// like_a_matrix.output(oomph_info);
103 ///
104 /// \endcode
105 ///
106 //================================================================
107 template<class KEY_TYPE_ROW, class KEY_TYPE_COL, class VALUE_TYPE>
109 {
110 public:
111 /// Default (empty) constructor
113
114 /// Broken assignment operator
115 void operator=(const MapMatrixMixed&) = delete;
116
117 /// Typedef to keep the code more readable
118 typedef std::map<KEY_TYPE_COL, VALUE_TYPE> InnerMapMixed;
119
120 /// Typedef to keep the code more readable
121 typedef typename InnerMapMixed::iterator InnerMixedIt;
122
123 /// Typedef to keep the code more readable const version
124 typedef typename InnerMapMixed::const_iterator ConstInnerMixedIt;
125
126 /// Typedef to keep the code more readable
127 typedef std::map<KEY_TYPE_ROW, std::map<KEY_TYPE_COL, VALUE_TYPE>*>
129
130 /// Typedef to keep the code more readable
131 typedef typename OuterMapMixed::iterator OuterMixedIt;
132
133 /// Typedef to keep the code more readable const version
134 typedef typename OuterMapMixed::const_iterator ConstOuterMixedIt;
135
136
137 /// Copy constructor
140 {
141 // Step through the row pointers
142 for (ConstOuterMixedIt it = map_mat.Row_pt.begin();
143 it != map_mat.Row_pt.end();
144 it++)
145 {
146 // Is the row pointer nonzero, i.e. are there any entries in this row?
147 if (it->second != 0)
148 {
149 // Identify the map that holds the entries in this row:
150 InnerMapMixed inner_map = *(it->second);
151
152 // Loop over entries in the row
153 for (ConstInnerMixedIt it2 = inner_map.begin();
154 it2 != inner_map.end();
155 it2++)
156 {
157 // If the entry is nonzero: Copy
158 if (it2->second != 0)
159 {
160 // key1, key2, value
161 (*this)(it->first, it2->first) = it2->second;
162 }
163 }
164 }
165 }
166 }
167
168
169 /// Copy a single column into its own map
170 void copy_column(const KEY_TYPE_COL& j,
171 std::map<KEY_TYPE_ROW, VALUE_TYPE>& copied_map)
172 {
173 // Step through the row pointers
174 for (OuterMixedIt it = Row_pt.begin(); it != Row_pt.end(); it++)
175 {
176 // Is the row pointer nonzero, i.e. are there any entries in this row?
177 if (it->second != 0)
178 {
179 // Identify the map that holds the entries in this row:
180 InnerMapMixed inner_map = *(it->second);
181 // If the desired column of the inner map is non-zero
182 if (inner_map[j] != 0)
183 {
184 // Set the value of the copied map to be the desired column of the
185 // inner map
186 copied_map[it->first] = inner_map[j];
187 }
188 }
189 }
190 }
191
192
193 /// Destructor
195 {
196 // Step through the pointers to row maps
197 for (OuterMixedIt it = Row_pt.begin(); it != Row_pt.end(); it++)
198 {
199 // Is the row pointer nonzero, i.e. are there any entries in this row?
200 if (it->second != 0)
201 {
202 // it->second is the stored object
203 delete it->second;
204 }
205 }
206 }
207
208 /// Wipe all entries
209 void clear()
210 {
211 // Step through the pointers to row maps
212 for (OuterMixedIt it = Row_pt.begin(); it != Row_pt.end(); it++)
213 {
214 // Is the row pointer nonzero, i.e. are there any entries in this row?
215 if (it->second != 0)
216 {
217 // it->second is the stored object: a map which can be cleared!
218 it->second->clear();
219 }
220 }
221 }
222
223
224 /// Return (reference to) entry.
225 /// Careful: If the entry does not exist then it is created and
226 /// set to zero
227 VALUE_TYPE& operator()(const KEY_TYPE_ROW& i, const KEY_TYPE_COL& j)
228 {
229 return *entry_pt(i, j);
230 }
231
232 /// Get an element corresponding to the key (i,j)
233 /// Searches the container for an element with a key equivalent to
234 /// (i,j) and returns the element if found, otherwise the default 0 value
235 /// for the value type is returned. The container is not modified.
236 VALUE_TYPE get(const KEY_TYPE_ROW& i, const KEY_TYPE_COL& j) const
237 {
238 if (Row_pt.count(i) > 0)
239 {
240 // Get the pointer to the row and check the j key
241 InnerMapMixed* inner_map_mixed_pt = Row_pt.find(i)->second;
242 if (inner_map_mixed_pt->count(j) > 0)
243 {
244 return inner_map_mixed_pt->find(j)->second;
245 }
246 else
247 {
248 return VALUE_TYPE(0);
249 }
250 }
251 else
252 {
253 // The key does not exist, return VALUE_TYPE(0)
254 return VALUE_TYPE(0);
255 }
256 }
257
258
259 /// Dump all non-`zero' entries to file.
260 /// Output is in the format
261 /// `i', `j', `entry[i][j]'
262 void output(std::ostream& outfile)
263 {
264 // NOTE:
265 //------
266 // map.first = key
267 // map.second = value
268
269 // Step through the row pointers
270 for (OuterMixedIt it = Row_pt.begin(); it != Row_pt.end(); it++)
271 {
272 // Is the row pointer nonzero, i.e. are there any entries in this row?
273 if (it->second != 0)
274 {
275 // Identify the map that holds the entries in this row:
276 InnerMapMixed inner_map = *(it->second);
277
278 // Loop over entries in the row
279 for (InnerMixedIt it2 = inner_map.begin(); it2 != inner_map.end();
280 it2++)
281 {
282 // If the entry is nonzero: Doc
283 if (it2->second != 0)
284 {
285 // Output key1, key2, value
286 outfile << it->first << " " << it2->first << " " << it2->second
287 << std::endl;
288 }
289 }
290 }
291 }
292 }
293
294 /// Work out number of non-`zero' entries
295 unsigned long nnz()
296 {
297 // Initialise counter for # of nonzero entries
298 unsigned long count = 0;
299
300 // Step through the row pointers
301 for (OuterMixedIt it = Row_pt.begin(); it != Row_pt.end(); it++)
302 {
303 // Is the row pointer nonzero, i.e. are there any entries in this row?
304 if (it->second != 0)
305 {
306 // Identify the map that holds the entries in this row:
307 InnerMapMixed inner_map = *(it->second);
308
309 // Loop over entries in the row
310 for (InnerMixedIt it2 = inner_map.begin(); it2 != inner_map.end();
311 it2++)
312 {
313 // If the entry is nonzero: Increment counter
314 if (it2->second != 0)
315 {
316 count++;
317 }
318 }
319 }
320 }
321 return count;
322 }
323
324 /// Work out number of non-`zero' entries, const version
325 unsigned long nnz() const
326 {
327 // Initialise counter for # of nonzero entries
328 unsigned long count = 0;
329
330 // Step through the row pointers
331 for (ConstOuterMixedIt it = Row_pt.begin(); it != Row_pt.end(); it++)
332 {
333 // Is the row pointer nonzero, i.e. are there any entries in this row?
334 if (it->second != 0)
335 {
336 // Identify the map that holds the entries in this row:
337 InnerMapMixed inner_map = *(it->second);
338
339 // Loop over entries in the row
340 for (ConstInnerMixedIt it2 = inner_map.begin();
341 it2 != inner_map.end();
342 it2++)
343 {
344 // If the entry is nonzero: Increment counter
345 if (it2->second != 0)
346 {
347 count++;
348 }
349 }
350 }
351 }
352 return count;
353 }
354
355
356 /// Work out total number of entries
357 unsigned long size()
358 {
359 // Initialise counter for # of nonzero entries
360 unsigned long count = 0;
361
362 // Step through the row pointers
363 for (OuterMixedIt it = Row_pt.begin(); it != Row_pt.end(); it++)
364 {
365 // Is the row pointer nonzero, i.e. are there any entries in this row?
366 if (it->second != 0)
367 {
368 // Identify the map that holds the entries in this row:
369 InnerMapMixed inner_map = *(it->second);
370
371 // Loop over all (!) entries in the row (incl. zero ones!)
372 for (InnerMixedIt it2 = inner_map.begin(); it2 != inner_map.end();
373 it2++)
374 {
375 count++;
376 }
377 }
378 }
379 return count;
380 }
381
382 /// Work out total number of entries const version
383 unsigned long size() const
384 {
385 // Initialise counter for # of nonzero entries
386 unsigned long count = 0;
387
388 // Step through the row pointers
389 for (ConstOuterMixedIt it = Row_pt.begin(); it != Row_pt.end(); it++)
390 {
391 // Is the row pointer nonzero, i.e. are there any entries in this row?
392 if (it->second != 0)
393 {
394 // Identify the map that holds the entries in this row:
395 InnerMapMixed inner_map = *(it->second);
396
397 // Loop over all (!) entries in the row (incl. zero ones!)
398 for (ConstInnerMixedIt it2 = inner_map.begin();
399 it2 != inner_map.end();
400 it2++)
401 {
402 count++;
403 }
404 }
405 }
406 return count;
407 }
408
409
410 protected:
411 /// Return pointer to entry
412 VALUE_TYPE* entry_pt(const KEY_TYPE_ROW& i, const KEY_TYPE_COL& j)
413 {
414 // There's not a single entry in this row: Entry must be zero.
415 if (Row_pt[i] == 0)
416 {
417 // Create row and entry in row and set the value to zero
418 Row_pt[i] = new std::map<KEY_TYPE_COL, VALUE_TYPE>;
419 (*Row_pt[i])[j] = VALUE_TYPE(0);
420 return &(*Row_pt[i])[j];
421 }
422 // Simply return pointer to existing entry
423 else
424 {
425 return &(*Row_pt[i])[j];
426 }
427 }
428
429
430 /// Here's the generalised matrix structure: A map of pointers to
431 /// the maps that hold the entries in each row.
432 std::map<KEY_TYPE_ROW, std::map<KEY_TYPE_COL, VALUE_TYPE>*> Row_pt;
433 };
434
435
436 /// ////////////////////////////////////////////////////////////////////
437 /// ////////////////////////////////////////////////////////////////////
438 /// ////////////////////////////////////////////////////////////////////
439
440
441 //================================================================
442 /// MapMatrix is a generalised, STL-map-based, sparse(-ish) matrix
443 /// class.
444 ///
445 /// The matrix is indexed by indices of type KEY_TYPE
446 /// and has entries of type VALUE_TYPE. It is a specialisation of the
447 /// class MapMatrixMixed. Please implement future functions in that class.
448 ///
449 /// Careful: If a zero entry is referenced then it is created in
450 /// memory. Therefore this isn't really a practical sparse matrix scheme.
451 /// Do not loop over `all' possible indices as even looking
452 /// at them will inflate the matrix until it occupies as much
453 /// space as a full one -- use (modification of) output routine
454 /// to retrieve all nonzero entries.
455 ///
456 /// However, this is not a serious restriction, as the main purpose
457 /// of this class is to allow non-integer indices.
458 ///
459 /// Example of usage:
460 /// \code
461 ///
462 ///
463 /// // Assume we have a Vector of pointers to objects:
464 /// Vector<Rubbish*> object_pt;
465 ///
466 /// [...]
467 ///
468 /// // Number of objects
469 /// int nentry=object_pt.size();
470 ///
471 /// // Use the pointers to the objects as indices
472 /// // in a MapMatrix whose entries are of type int
473 /// MapMatrix<Rubbish*,int> like_a_matrix;
474 ///
475 /// for (int i=1;i<nentry;i++)
476 /// {
477 /// for (int j=1;j<nentry;j++)
478 /// {
479 /// int number=100*i+j;
480 /// like_a_matrix(object_pt[i],object_pt[j])=number;
481 /// }
482 /// }
483 ///
484 /// oomph_info << "Matrix has nnz() " << like_a_matrix.nnz() <<
485 /// " and size() " << like_a_matrix.size() << std::endl;
486 ///
487 /// oomph_info << "\n\n\n Here are the nonzero entries: i, j, a(i,j)\n";
488 ///
489 /// like_a_matrix.output(oomph_info);
490 ///
491 /// // Can be used like a normal matrix:
492 ///
493 /// //Add to existing entry
494 /// like_a_matrix(object_pt[1],object_pt[2])+=13;
495 ///
496 /// // Add to non-existing entry
497 /// Rubbish* temp_pt=new Rubbish(20);
498 /// like_a_matrix(object_pt[1],temp_pt)+=13;
499 ///
500 /// oomph_info << "\n\n\n Here are the nonzero entries: i, j, a(i,j)\n";
501 /// like_a_matrix.output(oomph_info);
502 ///
503 /// \endcode
504 ///
505 //================================================================
506 template<class KEY_TYPE, class VALUE_TYPE>
507 class MapMatrix : public MapMatrixMixed<KEY_TYPE, KEY_TYPE, VALUE_TYPE>
508 {
509 public:
510 /// Default (empty) constructor
512
513 /// Typedef to keep the code more readable
514 typedef std::map<KEY_TYPE, VALUE_TYPE> InnerMap;
515
516 /// Typedef to keep the code more readable
517 typedef typename InnerMap::iterator InnerIt;
518
519 /// Typedef to keep the code more readable
520 typedef typename InnerMap::const_iterator ConstInnerIt;
521
522 /// Typedef to keep the code more readable
523 typedef std::map<KEY_TYPE, std::map<KEY_TYPE, VALUE_TYPE>*> OuterMap;
524
525 /// Typedef to keep the code more readable
526 typedef typename OuterMap::iterator OuterIt;
527
528 /// Typedef to keep the code more readable
529 typedef typename OuterMap::const_iterator ConstOuterIt;
530
531 /// Copy constructor
533 {
534 // Step through the row pointers
535 for (ConstOuterIt it = map_mat.Row_pt.begin(); it != map_mat.Row_pt.end();
536 it++)
537 {
538 // Is the row pointer nonzero, i.e. are there any entries in this row?
539 if (it->second != 0)
540 {
541 // Identify the map that holds the entries in this row:
542 InnerMap inner_map = *(it->second);
543
544 // Loop over entries in the row
545 for (ConstInnerIt it2 = inner_map.begin(); it2 != inner_map.end();
546 it2++)
547 {
548 // If the entry is nonzero: Copy
549 if (it2->second != 0)
550 {
551 (*this)(it->first, it2->first) = it2->second;
552 }
553 }
554 }
555 }
556 }
557
558 /// Broken assignment operator
559 void operator=(const MapMatrix&) = delete;
560 };
561
562} // namespace oomph
563
564#endif
cstr elem_len * i
Definition: cfortran.h:603
MapMatrixMixed is a generalised, STL-map-based, sparse(ish) matrix class with mixed indices.
Definition: map_matrix.h:109
VALUE_TYPE & operator()(const KEY_TYPE_ROW &i, const KEY_TYPE_COL &j)
Return (reference to) entry. Careful: If the entry does not exist then it is created and set to zero.
Definition: map_matrix.h:227
InnerMapMixed::const_iterator ConstInnerMixedIt
Typedef to keep the code more readable const version.
Definition: map_matrix.h:124
void operator=(const MapMatrixMixed &)=delete
Broken assignment operator.
VALUE_TYPE * entry_pt(const KEY_TYPE_ROW &i, const KEY_TYPE_COL &j)
Return pointer to entry.
Definition: map_matrix.h:412
unsigned long nnz()
Work out number of non-‘zero’ entries.
Definition: map_matrix.h:295
void copy_column(const KEY_TYPE_COL &j, std::map< KEY_TYPE_ROW, VALUE_TYPE > &copied_map)
Copy a single column into its own map.
Definition: map_matrix.h:170
MapMatrixMixed()
Default (empty) constructor.
Definition: map_matrix.h:112
unsigned long nnz() const
Work out number of non-‘zero’ entries, const version.
Definition: map_matrix.h:325
std::map< KEY_TYPE_COL, VALUE_TYPE > InnerMapMixed
Typedef to keep the code more readable.
Definition: map_matrix.h:118
virtual ~MapMatrixMixed()
Destructor.
Definition: map_matrix.h:194
InnerMapMixed::iterator InnerMixedIt
Typedef to keep the code more readable.
Definition: map_matrix.h:121
void output(std::ostream &outfile)
Dump all non-‘zero’ entries to file. Output is in the format ‘i’, ‘j’, ‘entry[i][j]’.
Definition: map_matrix.h:262
OuterMapMixed::const_iterator ConstOuterMixedIt
Typedef to keep the code more readable const version.
Definition: map_matrix.h:134
void clear()
Wipe all entries.
Definition: map_matrix.h:209
std::map< KEY_TYPE_ROW, std::map< KEY_TYPE_COL, VALUE_TYPE > * > Row_pt
Here's the generalised matrix structure: A map of pointers to the maps that hold the entries in each ...
Definition: map_matrix.h:432
unsigned long size() const
Work out total number of entries const version.
Definition: map_matrix.h:383
MapMatrixMixed(const MapMatrixMixed< KEY_TYPE_ROW, KEY_TYPE_COL, VALUE_TYPE > &map_mat)
Copy constructor.
Definition: map_matrix.h:138
OuterMapMixed::iterator OuterMixedIt
Typedef to keep the code more readable.
Definition: map_matrix.h:131
VALUE_TYPE get(const KEY_TYPE_ROW &i, const KEY_TYPE_COL &j) const
Get an element corresponding to the key (i,j) Searches the container for an element with a key equiva...
Definition: map_matrix.h:236
unsigned long size()
Work out total number of entries.
Definition: map_matrix.h:357
std::map< KEY_TYPE_ROW, std::map< KEY_TYPE_COL, VALUE_TYPE > * > OuterMapMixed
Typedef to keep the code more readable.
Definition: map_matrix.h:128
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: map_matrix.h:508
std::map< KEY_TYPE, std::map< KEY_TYPE, VALUE_TYPE > * > OuterMap
Typedef to keep the code more readable.
Definition: map_matrix.h:523
InnerMap::iterator InnerIt
Typedef to keep the code more readable.
Definition: map_matrix.h:517
MapMatrix()
Default (empty) constructor.
Definition: map_matrix.h:511
OuterMap::const_iterator ConstOuterIt
Typedef to keep the code more readable.
Definition: map_matrix.h:529
OuterMap::iterator OuterIt
Typedef to keep the code more readable.
Definition: map_matrix.h:526
InnerMap::const_iterator ConstInnerIt
Typedef to keep the code more readable.
Definition: map_matrix.h:520
MapMatrix(const MapMatrix< KEY_TYPE, VALUE_TYPE > &map_mat)
Copy constructor.
Definition: map_matrix.h:532
void operator=(const MapMatrix &)=delete
Broken assignment operator.
std::map< KEY_TYPE, VALUE_TYPE > InnerMap
Typedef to keep the code more readable.
Definition: map_matrix.h:511
//////////////////////////////////////////////////////////////////// ////////////////////////////////...