OpenMEEG
symm_block_matrix.h
Go to the documentation of this file.
1 // Project Name: OpenMEEG (http://openmeeg.github.io)
2 // © INRIA and ENPC under the French open source license CeCILL-B.
3 // See full copyright notice in the file LICENSE.txt
4 // If you make a copy of this file, you must either:
5 // - provide also LICENSE.txt and modify this header to refer to it.
6 // - replace this header by the LICENSE.txt content.
7 
8 #pragma once
9 
10 #include <iostream>
11 #include <map>
12 #include <algorithm>
13 
14 #include <linop.h>
15 #include <range.h>
16 #include <ranges.h>
17 #include <matrix.h>
18 #include <OMMathExceptions.H>
19 
20 namespace OpenMEEG::maths {
21 
24 
25  class SymmetricBlockMatrix: public LinOp {
26 
27  typedef std::pair<unsigned,unsigned> Index;
28  typedef std::map<Index,Matrix> Blocks;
29 
30  public:
31 
33  SymmetricBlockMatrix(const size_t N): LinOp(N,N,BLOCK_SYMMETRIC,2) { }
34 
35  size_t size() const override {
36  unsigned sz = 0;
37  for (const auto& block : blocks)
38  sz += block.second.size();
39  return sz;
40  };
41 
42  void info() const override {
43  if ((nlin()==0) && (ncol()==0)) {
44  std::cout << "Empty matrix" << std::endl;
45  return;
46  }
47 
48  std::cout << "Symmetric block matrix" << std::endl;
49  std::cout << "Dimensions: " << nlin() << " x " << ncol() << std::endl;
50  std::cout << "Number of blocks: " << blocks.size() << std::endl;
51  std::cout << "Number of coefficients: " << size() << std::endl;
52  }
53 
54  Matrix& block(const unsigned i,const unsigned j) {
55  bool transposed;
56  const Index& ind = find_block_indices(i,j,transposed);
57  if (transposed)
58  throw 1;
59  return blocks[ind];
60  }
61 
62  const Matrix& block(const unsigned i,const unsigned j) const {
63  bool transposed;
64  const Index& ind = find_block_indices(i,j,transposed);
65  if (transposed)
66  throw 1;
67  return blocks.at(ind);
68  }
69 
70  void add_block(const Range& ir,const Range& jr) {
71  bool transposed;
72  const Index& ind = create_block_indices(ir,jr,transposed);
73  const Index size = (transposed) ? Index({jr.length(),ir.length()}) : Index({ir.length(),jr.length()});
74  blocks[ind] = Matrix(size.first,size.second);
75  }
76 
77  void set_blocks(const Ranges& r) {
78  blocks.clear();
79  ranges.clear();
80  for (unsigned i=0; i<r.size(); ++i)
81  for (unsigned j=i; j<r.size(); ++j)
82  add_block(r[i],r[j]);
83  }
84 
85  double& operator()(const size_t i,const size_t j) {
86  bool transposed;
87  const Index& ind = find_block_indices(i,j,transposed);
88  const size_t inblockindex_i = ((!transposed) ? i : j)-ranges[ind.first].start();
89  const size_t inblockindex_j = ((!transposed) ? j : i)-ranges[ind.second].start();
90  return blocks[ind](inblockindex_i,inblockindex_j);
91  }
92 
93  double operator()(const size_t i,const size_t j) const {
94  bool transposed;
95  const Index& ind = find_block_indices(i,j,transposed);
96  const size_t inblockindex_i = ((!transposed) ? i : j)-ranges[ind.first].start();
97  const size_t inblockindex_j = ((!transposed) ? j : i)-ranges[ind.second].start();
98  return blocks.at(ind)(inblockindex_i,inblockindex_j);
99  }
100 
101  private:
102 
103  unsigned find_block_index(const Range& r) const { return ranges.find_index(r); }
104  unsigned find_block_index(const unsigned i) const { return ranges.find_index(i); }
105 
106  unsigned create_block_index(const Range& r) try {
107  const unsigned ind = ranges.find_index(r);
108  return ind;
109  } catch(...) {
110  ranges.push_back(r);
111  return ranges.size()-1;
112  }
113 
114  Index create_block_indices(const Range& ir,const Range& jr,bool& transposed) {
115  transposed = ir.start()>jr.start();
116  const unsigned iind = create_block_index(ir);
117  const unsigned jind = create_block_index(jr);
118  return (transposed) ? Index({jind,iind}) : Index({iind,jind});
119  }
120 
121  Index find_block_indices(const Range& ir,const Range& jr,bool& transposed) const {
122  transposed = ir.start()>jr.start();
123  const unsigned iind = find_block_index(ir);
124  const unsigned jind = find_block_index(jr);
125  return (transposed) ? Index({jind,iind}) : Index({iind,jind});
126  }
127 
128  Index find_block_indices(const unsigned i,const unsigned j,bool& transposed) const {
129  transposed = i>j;
130  const unsigned iind = find_block_index(i);
131  const unsigned jind = find_block_index(j);
132  return (transposed) ? Index({jind,iind}) : Index({iind,jind});
133  }
134 
135  Ranges ranges;
136  Blocks blocks;
137  };
138 }
Dimension nlin() const
Definition: linop.h:48
virtual Dimension ncol() const
Definition: linop.h:51
Matrix class Matrix class.
Definition: matrix.h:28
size_t size() const
Get Matrix size.
Definition: matrix.h:61
size_t length() const
Definition: range.h:24
unsigned find_index(const size_t ind) const
Definition: ranges.h:39
Block symmetric matrix class Block symmetric matrix class.
const Matrix & block(const unsigned i, const unsigned j) const
void add_block(const Range &ir, const Range &jr)
Matrix & block(const unsigned i, const unsigned j)
double operator()(const size_t i, const size_t j) const
double & operator()(const size_t i, const size_t j)
unsigned Index
Definition: linop.h:33