OpenMEEG
Loading...
Searching...
No Matches
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
19namespace 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.
BlockMatrix(const size_t M, const size_t N)
double operator()(const size_t i, const size_t j) const
void set_blocks(const Ranges &rows, const Ranges &cols)
void add_block(const Range &ir, const Range &jr)
double & operator()(const size_t i, const size_t j)
void info() const override
size_t size() const override
const Blocks & blocks() const
const Matrix & block(const unsigned i, const unsigned j) const
Matrix & block(const unsigned i, const unsigned j)
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)
unsigned Index
Definition linop.h:33