OpenMEEG
Loading...
Searching...
No Matches
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
20namespace OpenMEEG::maths {
21
24
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.
void add_block(const Range &ir, const Range &jr)
const Matrix & block(const unsigned i, const unsigned j) const
Matrix & block(const unsigned i, const unsigned j)
double & operator()(const size_t i, const size_t j)
double operator()(const size_t i, const size_t j) const
unsigned Index
Definition linop.h:33