Alien  1.3.0
Developer documentation
Loading...
Searching...
No Matches
SchurBlock.h
1
2//-----------------------------------------------------------------------------
3// Copyright 2000-2023 CEA (www.cea.fr) IFPEN (www.ifpenergiesnouvelles.com)
4// See the top-level COPYRIGHT file for details.
5// SPDX-License-Identifier: Apache-2.0
6//-----------------------------------------------------------------------------
7#pragma once
8
9#ifdef ALIEN_USE_EIGEN3
10#include <Eigen/Core>
11#include <Eigen/LU>
12#endif
13
14/*---------------------------------------------------------------------------*/
15/*---------------------------------------------------------------------------*/
16
17namespace Alien
18{
19
20/*---------------------------------------------------------------------------*/
21/*---------------------------------------------------------------------------*/
22
23class ALIEN_EXPORT SchurBlock1D
24{
25 public:
26#ifdef ALIEN_USE_EIGEN3
27 typedef Eigen::Matrix<Real, Eigen::Dynamic, 1> vector;
28 typedef Eigen::Map<vector> eigen_vector;
29#endif
30 SchurBlock1D(ArrayView<Real> block, Integer size)
31 : m_block(block)
32 , m_size_1(std::min(size, m_block.size()))
33 , m_size_2(m_block.size() - m_size_1)
34#ifdef ALIEN_USE_EIGEN3
35 , m_block_1(m_block.data(), m_size_1)
36 , m_block_2(m_block.data() + m_size_1, m_size_2)
37#endif
38 {
39 }
40
41 virtual ~SchurBlock1D() {}
42
43 ConstArrayView<Real> block() const { return m_block; }
44
45#ifdef ALIEN_USE_EIGEN3
46 eigen_vector& block_1() { return m_block_1; }
47 const eigen_vector& block_1() const { return m_block_1; }
48
49 eigen_vector& block_2() { return m_block_2; }
50 const eigen_vector& block_2() const { return m_block_2; }
51#endif
52 private:
53 ArrayView<Real> m_block;
54
55 Integer m_size_1;
56 Integer m_size_2;
57
58#ifdef ALIEN_USE_EIGEN3
59 eigen_vector m_block_1;
60 eigen_vector m_block_2;
61#endif
62};
63
64/*---------------------------------------------------------------------------*/
65/*---------------------------------------------------------------------------*/
66
67class ALIEN_EXPORT SchurBlock2D
68{
69 public:
70#ifdef ALIEN_USE_EIGEN3
71 typedef Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> matrix;
72 typedef Eigen::OuterStride<> outer_stride;
73 typedef Eigen::Map<matrix, 0, outer_stride> eigen_matrix;
74#endif
75 SchurBlock2D(Array2View<Real> block, Integer size)
76 : m_block(block)
77 , m_sizex_1(std::min(size, m_block.dim1Size()))
78 , m_sizex_2(m_block.dim1Size() - m_sizex_1)
79 , m_sizey_1(std::min(size, m_block.dim2Size()))
80 , m_sizey_2(m_block.dim2Size() - m_sizey_1)
81#ifdef ALIEN_USE_EIGEN3
82 , m_block_11(m_block.unguardedBasePointer(),
83 m_sizex_1, m_sizey_1,
84 outer_stride(m_block.dim2Size()))
85 , m_block_12(m_block.unguardedBasePointer() + m_sizex_1,
86 m_sizex_1, m_sizey_2,
87 outer_stride(m_block.dim2Size()))
88 , m_block_21(m_block.unguardedBasePointer() + (m_sizex_1 * m_block.dim2Size()),
89 m_sizex_2, m_sizey_1,
90 outer_stride(m_block.dim2Size()))
91 , m_block_22(m_block.unguardedBasePointer() + (m_sizex_1 * (m_block.dim2Size() + 1)),
92 m_sizex_2, m_sizey_2,
93 outer_stride(m_block.dim2Size()))
94#endif
95 {
96 }
97
98 virtual ~SchurBlock2D() {}
99
100 ConstArray2View<Real> block() const { return m_block; }
101
102#ifdef ALIEN_USE_EIGEN3
103 eigen_matrix& block_11() { return m_block_11; }
104 const eigen_matrix& block_11() const { return m_block_11; }
105
106 eigen_matrix& block_12() { return m_block_12; }
107 const eigen_matrix& block_12() const { return m_block_12; }
108
109 eigen_matrix& block_21() { return m_block_21; }
110 const eigen_matrix& block_21() const { return m_block_21; }
111
112 eigen_matrix& block_22() { return m_block_22; }
113 const eigen_matrix& block_22() const { return m_block_22; }
114#endif
115 private:
116 Array2View<Real> m_block;
117
118 Integer m_sizex_1;
119 Integer m_sizex_2;
120 Integer m_sizey_1;
121 Integer m_sizey_2;
122
123#ifdef ALIEN_USE_EIGEN3
124 eigen_matrix m_block_11;
125 eigen_matrix m_block_12;
126 eigen_matrix m_block_21;
127 eigen_matrix m_block_22;
128#endif
129};
130
131/*---------------------------------------------------------------------------*/
132/*---------------------------------------------------------------------------*/
133
134class ALIEN_EXPORT SchurAlgo
135{
136 public:
137 SchurAlgo() {}
138
139 static void compute(SchurBlock2D& A, SchurBlock1D& b)
140 {
141 ALIEN_ASSERT((A.block().dim1Size() == A.block().dim2Size()), ("Schur complement must be applied on square block"));
142 ALIEN_ASSERT((A.block().dim1Size() == b.block().size()), ("Blocks size are not equals"));
143 ;
144
145#ifdef ALIEN_USE_EIGEN3
146 typedef SchurBlock2D::eigen_matrix matrix;
147 typedef SchurBlock1D::eigen_vector vector;
148
149 matrix& matrix_11 = A.block_11();
150 matrix& matrix_12 = A.block_12();
151 matrix& matrix_21 = A.block_21();
152 matrix& matrix_22 = A.block_22();
153
154 vector& vector_1 = b.block_1();
155 vector& vector_2 = b.block_2();
156
157 /*
158 A22 => A22^-1
159 A21 => A22^-1 * A 21
160 A11 => A11 - A12 * A22^-1 A21
161
162 b2 => A22^-1 * b2
163 b1 => b1 - A12 * A22^-1 * b2
164 */
165
166 matrix_22 = matrix_22.inverse();
167 matrix_21 = matrix_22 * matrix_21;
168 matrix_11 -= matrix_12 * matrix_21;
169
170 vector_2 = matrix_22 * vector_2;
171 vector_1 -= matrix_12 * vector_2;
172#endif
173 }
174};
175
176/*---------------------------------------------------------------------------*/
177/*---------------------------------------------------------------------------*/
178
179} // namespace Alien
180
181/*---------------------------------------------------------------------------*/
182/*---------------------------------------------------------------------------*/
Alien
-- tab-width: 2; indent-tabs-mode: nil; coding: utf-8-with-signature --
Definition BackEnd.h:17