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-2023 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 
46 namespace 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
194  virtual ~MapMatrixMixed()
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
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.
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
VALUE_TYPE * entry_pt(const KEY_TYPE_ROW &i, const KEY_TYPE_COL &j)
Return pointer to entry.
Definition: map_matrix.h:412
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
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
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
//////////////////////////////////////////////////////////////////// ////////////////////////////////...