OpenMEEG
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 
19 namespace OpenMEEG::maths {
20 
23 
24  class BlockMatrix: public LinOp {
25 
26  typedef std::pair<unsigned,unsigned> Index;
27  typedef std::map<Index,Matrix> Blocks;
28 
29  public:
30 
31  BlockMatrix(): LinOp(0,0,BLOCK,2) { }
32  BlockMatrix(const size_t M,const size_t N): LinOp(M,N,BLOCK,2) { }
33 
34  size_t size() const override {
35  unsigned sz = 0;
36  for (const auto& block : all_blocks)
37  sz += block.second.size();
38  return sz;
39  };
40 
41  void info() const override {
42  if ((nlin()==0) && (ncol()==0)) {
43  std::cout << "Empty matrix" << std::endl;
44  return;
45  }
46 
47  std::cout << "Block matrix" << std::endl;
48  std::cout << "Dimensions: " << nlin() << " x " << ncol() << std::endl;
49  std::cout << "Number of blocks: " << all_blocks.size() << std::endl;
50  std::cout << "Number of coefficients: " << size() << std::endl;
51  }
52 
53  Matrix& block(const unsigned i,const unsigned j) { return all_blocks[{i,j}]; }
54  const Matrix& block(const unsigned i,const unsigned j) const { return all_blocks.at({i,j}); }
55 
56  const Blocks& blocks() const { return all_blocks; }
57 
58  void add_block(const Range& ir,const Range& jr) {
59  const unsigned iind = row_ranges.add(ir);
60  const unsigned jind = col_ranges.add(jr);
61  const Index inds = { iind, jind };
62  all_blocks[inds] = Matrix(ir.length(),jr.length());
63  }
64 
65  void set_blocks(const Ranges& rows,const Ranges& cols) {
66  row_ranges = rows;
67  col_ranges = cols;
68  for (const auto& ir : row_ranges)
69  for (const auto& jr : col_ranges)
70  add_block(ir,jr);
71  }
72 
73  double& operator()(const size_t i,const size_t j) {
74  const Index& ind = find_block_indices(i,j);
75  const size_t inblockindex_i = i-row_ranges[ind.first].start();
76  const size_t inblockindex_j = j-col_ranges[ind.second].start();
77  return all_blocks[ind](inblockindex_i,inblockindex_j);
78  }
79 
80  double operator()(const size_t i,const size_t j) const {
81  const Index& ind = find_block_indices(i,j);
82  const size_t inblockindex_i = i-row_ranges[ind.first].start();
83  const size_t inblockindex_j = j-col_ranges[ind.second].start();
84  return all_blocks.at(ind)(inblockindex_i,inblockindex_j);
85  }
86 
87  private:
88 
89  Index find_block_indices(const Range& ir,const Range& jr) const {
90  const unsigned iind = row_ranges.find_index(ir);
91  const unsigned jind = col_ranges.find_index(jr);
92  return {iind,jind};
93  }
94 
95  Index find_block_indices(const unsigned i,const unsigned j) const {
96  const unsigned iind = row_ranges.find_index(i);
97  const unsigned jind = col_ranges.find_index(j);
98  return { iind, jind };
99  }
100 
101  Ranges row_ranges;
102  Ranges col_ranges;
103  Blocks all_blocks;
104  };
105 
106  inline std::ostream& operator<<(std::ostream& os,const BlockMatrix& bm) {
107  for (const auto& block : bm.blocks())
108  os << "Block " << block.first.first << ',' << block.first.second << std::endl;
109  #if 0
110  os << "Block " << block.first.first << ',' << block.first.second << std::endl
111  << block.second << std::endl;
112  #endif
113  return os;
114  }
115 }
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
Block matrix class Block matrix class.
Definition: block_matrix.h:24
double & operator()(const size_t i, const size_t j)
Definition: block_matrix.h:73
BlockMatrix(const size_t M, const size_t N)
Definition: block_matrix.h:32
Matrix & block(const unsigned i, const unsigned j)
Definition: block_matrix.h:53
const Blocks & blocks() const
Definition: block_matrix.h:56
const Matrix & block(const unsigned i, const unsigned j) const
Definition: block_matrix.h:54
double operator()(const size_t i, const size_t j) const
Definition: block_matrix.h:80
void set_blocks(const Ranges &rows, const Ranges &cols)
Definition: block_matrix.h:65
void add_block(const Range &ir, const Range &jr)
Definition: block_matrix.h:58
void info() const override
Definition: block_matrix.h:41
size_t size() const override
Definition: block_matrix.h:34
size_t length() const
Definition: range.h:24
unsigned add(const Range &range)
Definition: ranges.h:28
unsigned find_index(const size_t ind) const
Definition: ranges.h:39
std::ostream & operator<<(std::ostream &os, const BlockMatrix &bm)
Definition: block_matrix.h:106
unsigned Index
Definition: linop.h:33