deal.II version GIT relicensing-3540-g7552a02177 2025-06-20 13:50:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
matrix_tools.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 1998 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
19
23
24#include <deal.II/fe/fe.h>
26
28
33#include <deal.II/lac/vector.h>
34
36
37#ifdef DEAL_II_WITH_PETSC
41#endif
42
43#ifdef DEAL_II_WITH_TRILINOS
48#endif
49
50#include <algorithm>
51#include <cmath>
52
53
55
56
57
58namespace MatrixTools
59{
60 namespace
61 {
62 template <typename Iterator>
63 bool
64 column_less_than(const typename Iterator::value_type p,
65 const unsigned int column)
66 {
67 return (p.column() < column);
68 }
69 } // namespace
70
71 // TODO:[WB] I don't think that the optimized storage of diagonals is needed
72 // (GK)
73 template <typename number>
74 void
76 const std::map<types::global_dof_index, number> &boundary_values,
78 Vector<number> &solution,
79 Vector<number> &right_hand_side,
80 const bool eliminate_columns)
81 {
82 Assert(matrix.n() == right_hand_side.size(),
83 ExcDimensionMismatch(matrix.n(), right_hand_side.size()));
84 Assert(matrix.n() == solution.size(),
85 ExcDimensionMismatch(matrix.n(), solution.size()));
86 Assert(matrix.n() == matrix.m(),
87 ExcDimensionMismatch(matrix.n(), matrix.m()));
88
89 // if no boundary values are to be applied
90 // simply return
91 if (boundary_values.empty())
92 return;
93
94
95 const types::global_dof_index n_dofs = matrix.m();
96
97 // if a diagonal entry is zero
98 // later, then we use another
99 // number instead. take it to be
100 // the first nonzero diagonal
101 // element of the matrix, or 1 if
102 // there is no such thing
103 number first_nonzero_diagonal_entry = 1;
104 for (unsigned int i = 0; i < n_dofs; ++i)
105 if (matrix.diag_element(i) != number())
106 {
107 first_nonzero_diagonal_entry = matrix.diag_element(i);
108 break;
109 }
110
111
112 typename std::map<types::global_dof_index, number>::const_iterator
113 dof = boundary_values.begin(),
114 endd = boundary_values.end();
115 for (; dof != endd; ++dof)
116 {
117 Assert(dof->first < n_dofs, ExcInternalError());
118
119 const types::global_dof_index dof_number = dof->first;
120 // for each boundary dof:
121
122 // set entries of this line to zero except for the diagonal
123 // entry
124 for (typename SparseMatrix<number>::iterator p =
125 matrix.begin(dof_number);
126 p != matrix.end(dof_number);
127 ++p)
128 if (p->column() != dof_number)
129 p->value() = 0.;
130
131 // set right hand side to
132 // wanted value: if main diagonal
133 // entry nonzero, don't touch it
134 // and scale rhs accordingly. If
135 // zero, take the first main
136 // diagonal entry we can find, or
137 // one if no nonzero main diagonal
138 // element exists. Normally, however,
139 // the main diagonal entry should
140 // not be zero.
141 //
142 // store the new rhs entry to make
143 // the gauss step more efficient
144 number new_rhs;
145 if (matrix.diag_element(dof_number) != number())
146 {
147 new_rhs = dof->second * matrix.diag_element(dof_number);
148 right_hand_side(dof_number) = new_rhs;
149 }
150 else
151 {
152 matrix.set(dof_number, dof_number, first_nonzero_diagonal_entry);
153 new_rhs = dof->second * first_nonzero_diagonal_entry;
154 right_hand_side(dof_number) = new_rhs;
155 }
156
157
158 // if the user wants to have
159 // the symmetry of the matrix
160 // preserved, and if the
161 // sparsity pattern is
162 // symmetric, then do a Gauss
163 // elimination step with the
164 // present row
165 if (eliminate_columns)
166 {
167 // store the only nonzero entry
168 // of this line for the Gauss
169 // elimination step
170 const number diagonal_entry = matrix.diag_element(dof_number);
171
172 // we have to loop over all rows of the matrix which have
173 // a nonzero entry in the column which we work in
174 // presently. if the sparsity pattern is symmetric, then
175 // we can get the positions of these rows cheaply by
176 // looking at the nonzero column numbers of the present
177 // row. we need not look at the first entry of each row,
178 // since that is the diagonal element and thus the present
179 // row
180 for (typename SparseMatrix<number>::iterator q =
181 matrix.begin(dof_number) + 1;
182 q != matrix.end(dof_number);
183 ++q)
184 {
185 const types::global_dof_index row = q->column();
186
187 // find the position of
188 // element
189 // (row,dof_number)
190 bool (*comp)(
192 const unsigned int column) =
193 &column_less_than<typename SparseMatrix<number>::iterator>;
194 const typename SparseMatrix<number>::iterator p =
195 Utilities::lower_bound(matrix.begin(row) + 1,
196 matrix.end(row),
197 dof_number,
198 comp);
199
200 // check whether this line has an entry in the
201 // regarding column (check for ==dof_number and !=
202 // next_row, since if row==dof_number-1, *p is a
203 // past-the-end pointer but points to dof_number
204 // anyway...)
205 //
206 // there should be such an entry! we know this because
207 // we have assumed that the sparsity pattern is
208 // symmetric and we only walk over those rows for
209 // which the current row has a column entry
210 Assert((p != matrix.end(row)) && (p->column() == dof_number),
212 "This function is trying to access an element of the "
213 "matrix that doesn't seem to exist. Are you using a "
214 "nonsymmetric sparsity pattern? If so, you are not "
215 "allowed to set the eliminate_column argument of this "
216 "function, see the documentation."));
217
218 // correct right hand side
219 right_hand_side(row) -=
220 static_cast<number>(p->value()) / diagonal_entry * new_rhs;
221
222 // set matrix entry to zero
223 p->value() = 0.;
224 }
225 }
226
227 // preset solution vector
228 solution(dof_number) = dof->second;
229 }
230 }
231
232
233
234 template <typename number>
235 void
237 const std::map<types::global_dof_index, number> &boundary_values,
239 BlockVector<number> &solution,
240 BlockVector<number> &right_hand_side,
241 const bool eliminate_columns)
242 {
243 const unsigned int blocks = matrix.n_block_rows();
244
245 Assert(matrix.n() == right_hand_side.size(),
246 ExcDimensionMismatch(matrix.n(), right_hand_side.size()));
247 Assert(matrix.n() == solution.size(),
248 ExcDimensionMismatch(matrix.n(), solution.size()));
249 Assert(matrix.n_block_rows() == matrix.n_block_cols(), ExcNotQuadratic());
250 Assert(matrix.get_sparsity_pattern().get_row_indices() ==
251 matrix.get_sparsity_pattern().get_column_indices(),
253 Assert(matrix.get_sparsity_pattern().get_column_indices() ==
254 solution.get_block_indices(),
256 Assert(matrix.get_sparsity_pattern().get_row_indices() ==
257 right_hand_side.get_block_indices(),
259
260 // if no boundary values are to be applied
261 // simply return
262 if (boundary_values.empty())
263 return;
264
265
266 const types::global_dof_index n_dofs = matrix.m();
267
268 // if a diagonal entry is zero
269 // later, then we use another
270 // number instead. take it to be
271 // the first nonzero diagonal
272 // element of the matrix, or 1 if
273 // there is no such thing
274 number first_nonzero_diagonal_entry = 0;
275 for (unsigned int diag_block = 0; diag_block < blocks; ++diag_block)
276 {
277 for (unsigned int i = 0; i < matrix.block(diag_block, diag_block).n();
278 ++i)
279 if (matrix.block(diag_block, diag_block).diag_element(i) != number{})
280 {
281 first_nonzero_diagonal_entry =
282 matrix.block(diag_block, diag_block).diag_element(i);
283 break;
284 }
285 // check whether we have found
286 // something in the present
287 // block
288 if (first_nonzero_diagonal_entry != number{})
289 break;
290 }
291 // nothing found on all diagonal
292 // blocks? if so, use 1.0 instead
293 if (first_nonzero_diagonal_entry == number{})
294 first_nonzero_diagonal_entry = 1;
295
296
297 typename std::map<types::global_dof_index, number>::const_iterator
298 dof = boundary_values.begin(),
299 endd = boundary_values.end();
300 const BlockSparsityPattern &sparsity_pattern =
301 matrix.get_sparsity_pattern();
302
303 // pointer to the mapping between
304 // global and block indices. since
305 // the row and column mappings are
306 // equal, store a pointer on only
307 // one of them
308 const BlockIndices &index_mapping = sparsity_pattern.get_column_indices();
309
310 // now loop over all boundary dofs
311 for (; dof != endd; ++dof)
312 {
313 Assert(dof->first < n_dofs, ExcInternalError());
314
315 // get global index and index
316 // in the block in which this
317 // dof is located
318 const types::global_dof_index dof_number = dof->first;
319 const std::pair<unsigned int, types::global_dof_index> block_index =
320 index_mapping.global_to_local(dof_number);
321
322 // for each boundary dof:
323
324 // set entries of this line
325 // to zero except for the diagonal
326 // entry. Note that the diagonal
327 // entry is always the first one
328 // in a row for square matrices
329 for (unsigned int block_col = 0; block_col < blocks; ++block_col)
330 for (typename SparseMatrix<number>::iterator p =
331 (block_col == block_index.first ?
332 matrix.block(block_index.first, block_col)
333 .begin(block_index.second) +
334 1 :
335 matrix.block(block_index.first, block_col)
336 .begin(block_index.second));
337 p != matrix.block(block_index.first, block_col)
338 .end(block_index.second);
339 ++p)
340 p->value() = 0;
341
342 // set right hand side to
343 // wanted value: if main diagonal
344 // entry nonzero, don't touch it
345 // and scale rhs accordingly. If
346 // zero, take the first main
347 // diagonal entry we can find, or
348 // one if no nonzero main diagonal
349 // element exists. Normally, however,
350 // the main diagonal entry should
351 // not be zero.
352 //
353 // store the new rhs entry to make
354 // the gauss step more efficient
355 number new_rhs;
356 if (matrix.block(block_index.first, block_index.first)
357 .diag_element(block_index.second) != number{})
358 new_rhs =
359 dof->second * matrix.block(block_index.first, block_index.first)
360 .diag_element(block_index.second);
361 else
362 {
363 matrix.block(block_index.first, block_index.first)
364 .diag_element(block_index.second) = first_nonzero_diagonal_entry;
365 new_rhs = dof->second * first_nonzero_diagonal_entry;
366 }
367 right_hand_side.block(block_index.first)(block_index.second) = new_rhs;
368
369
370 // if the user wants to have
371 // the symmetry of the matrix
372 // preserved, and if the
373 // sparsity pattern is
374 // symmetric, then do a Gauss
375 // elimination step with the
376 // present row. this is a
377 // little more complicated for
378 // block matrices.
379 if (eliminate_columns)
380 {
381 // store the only nonzero entry
382 // of this line for the Gauss
383 // elimination step
384 const number diagonal_entry =
385 matrix.block(block_index.first, block_index.first)
386 .diag_element(block_index.second);
387
388 // we have to loop over all
389 // rows of the matrix which
390 // have a nonzero entry in
391 // the column which we work
392 // in presently. if the
393 // sparsity pattern is
394 // symmetric, then we can
395 // get the positions of
396 // these rows cheaply by
397 // looking at the nonzero
398 // column numbers of the
399 // present row.
400 //
401 // note that if we check
402 // whether row @p{row} in
403 // block (r,c) is non-zero,
404 // then we have to check
405 // for the existence of
406 // column @p{row} in block
407 // (c,r), i.e. of the
408 // transpose block
409 for (unsigned int block_row = 0; block_row < blocks; ++block_row)
410 {
411 // get pointers to the sparsity patterns of this block and of
412 // the transpose one
413 const SparsityPattern &this_sparsity =
414 sparsity_pattern.block(block_row, block_index.first);
415
416 SparseMatrix<number> &this_matrix =
417 matrix.block(block_row, block_index.first);
418 SparseMatrix<number> &transpose_matrix =
419 matrix.block(block_index.first, block_row);
420
421 // traverse the row of the transpose block to find the
422 // interesting rows in the present block. don't use the
423 // diagonal element of the diagonal block
424 for (typename SparseMatrix<number>::iterator q =
425 (block_index.first == block_row ?
426 transpose_matrix.begin(block_index.second) + 1 :
427 transpose_matrix.begin(block_index.second));
428 q != transpose_matrix.end(block_index.second);
429 ++q)
430 {
431 // get the number of the column in this row in which a
432 // nonzero entry is. this is also the row of the transpose
433 // block which has an entry in the interesting row
434 const types::global_dof_index row = q->column();
435
436 // find the position of element (row,dof_number) in this
437 // block (not in the transpose one). note that we have to
438 // take care of special cases with square sub-matrices
439 bool (*comp)(
441 const unsigned int column) =
442 &column_less_than<
444
446 this_matrix.end();
447
448 if (this_sparsity.n_rows() == this_sparsity.n_cols())
449 {
450 if (this_matrix.begin(row)->column() ==
451 block_index.second)
452 p = this_matrix.begin(row);
453 else
454 p = Utilities::lower_bound(this_matrix.begin(row) + 1,
455 this_matrix.end(row),
456 block_index.second,
457 comp);
458 }
459 else
460 p = Utilities::lower_bound(this_matrix.begin(row),
461 this_matrix.end(row),
462 block_index.second,
463 comp);
464
465 // check whether this line has an entry in the
466 // regarding column (check for ==dof_number and !=
467 // next_row, since if row==dof_number-1, *p is a
468 // past-the-end pointer but points to dof_number
469 // anyway...)
470 //
471 // there should be such an entry! we know this because
472 // we have assumed that the sparsity pattern is
473 // symmetric and we only walk over those rows for
474 // which the current row has a column entry
475 Assert((p->column() == block_index.second) &&
476 (p != this_matrix.end(row)),
478
479 // correct right hand side
480 right_hand_side.block(block_row)(row) -=
481 number(p->value()) / diagonal_entry * new_rhs;
482
483 // set matrix entry to zero
484 p->value() = 0.;
485 }
486 }
487 }
488
489 // preset solution vector
490 solution.block(block_index.first)(block_index.second) = dof->second;
491 }
492 }
493
494
495
496 template <typename number>
497 void
499 const std::map<types::global_dof_index, number> &boundary_values,
500 const std::vector<types::global_dof_index> &local_dof_indices,
501 FullMatrix<number> &local_matrix,
502 Vector<number> &local_rhs,
503 const bool eliminate_columns)
504 {
505 Assert(local_dof_indices.size() == local_matrix.m(),
506 ExcDimensionMismatch(local_dof_indices.size(), local_matrix.m()));
507 Assert(local_dof_indices.size() == local_matrix.n(),
508 ExcDimensionMismatch(local_dof_indices.size(), local_matrix.n()));
509 Assert(local_dof_indices.size() == local_rhs.size(),
510 ExcDimensionMismatch(local_dof_indices.size(), local_rhs.size()));
511
512 // if there is nothing to do, then exit
513 // right away
514 if (boundary_values.empty())
515 return;
516
517 // otherwise traverse all the dofs used in
518 // the local matrices and vectors and see
519 // what's there to do
520
521 // if we need to treat an entry, then we
522 // set the diagonal entry to its absolute
523 // value. if it is zero, we used to set it
524 // to one, which is a really terrible
525 // choice that can lead to hours of
526 // searching for bugs in programs (I
527 // experienced this :-( ) if the matrix
528 // entries are otherwise very large. this
529 // is so since iterative solvers would
530 // simply not correct boundary nodes for
531 // their correct values since the residual
532 // contributions of their rows of the
533 // linear system is almost zero if the
534 // diagonal entry is one. thus, set it to
535 // the average absolute value of the
536 // nonzero diagonal elements.
537 //
538 // we only compute this value lazily the
539 // first time we need it.
540 number average_diagonal = 0;
541 const unsigned int n_local_dofs = local_dof_indices.size();
542 for (unsigned int i = 0; i < n_local_dofs; ++i)
543 {
544 const typename std::map<types::global_dof_index, number>::const_iterator
545 boundary_value = boundary_values.find(local_dof_indices[i]);
546 if (boundary_value != boundary_values.end())
547 {
548 // remove this row, except for the
549 // diagonal element
550 for (unsigned int j = 0; j < n_local_dofs; ++j)
551 if (i != j)
552 local_matrix(i, j) = 0;
553
554 // replace diagonal entry by its
555 // absolute value to make sure that
556 // everything remains positive, or
557 // by the average diagonal value if
558 // zero
559 if (local_matrix(i, i) == number{})
560 {
561 // if average diagonal hasn't
562 // yet been computed, do so now
563 if (average_diagonal == number{})
564 {
565 unsigned int nonzero_diagonals = 0;
566 for (unsigned int k = 0; k < n_local_dofs; ++k)
567 if (local_matrix(k, k) != number{})
568 {
569 average_diagonal += std::abs(local_matrix(k, k));
570 ++nonzero_diagonals;
571 }
572 if (nonzero_diagonals != 0)
573 average_diagonal /= nonzero_diagonals;
574 else
575 average_diagonal = 0;
576 }
577
578 // only if all diagonal entries
579 // are zero, then resort to the
580 // last measure: choose one
581 if (average_diagonal == number{})
582 average_diagonal = 1.;
583
584 local_matrix(i, i) = average_diagonal;
585 }
586 else
587 local_matrix(i, i) = std::abs(local_matrix(i, i));
588
589 // and replace rhs entry by correct
590 // value
591 local_rhs(i) = local_matrix(i, i) * boundary_value->second;
592
593 // finally do the elimination step
594 // if requested
595 if (eliminate_columns == true)
596 {
597 for (unsigned int row = 0; row < n_local_dofs; ++row)
598 if (row != i)
599 {
600 local_rhs(row) -=
601 local_matrix(row, i) * boundary_value->second;
602 local_matrix(row, i) = 0;
603 }
604 }
605 }
606 }
607 }
608} // namespace MatrixTools
609
610
611
612// explicit instantiations
613#include "numerics/matrix_tools.inst"
614
615
std::pair< unsigned int, size_type > global_to_local(const size_type i) const
SparsityPatternType & block(const size_type row, const size_type column)
const BlockIndices & get_column_indices() const
virtual size_type size() const override
BlockType & block(const unsigned int i)
const BlockIndices & get_block_indices() const
size_type n() const
size_type m() const
const Accessor< number, Constness > & value_type
const_iterator end() const
const_iterator begin() const
size_type n_rows() const
size_type n_cols() const
virtual size_type size() const override
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
static ::ExceptionBase & ExcBlocksDontMatch()
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcNotQuadratic()
static ::ExceptionBase & ExcMessage(std::string arg1)
void apply_boundary_values(const std::map< types::global_dof_index, number > &boundary_values, SparseMatrix< number > &matrix, Vector< number > &solution, Vector< number > &right_hand_side, const bool eliminate_columns=true)
void local_apply_boundary_values(const std::map< types::global_dof_index, number > &boundary_values, const std::vector< types::global_dof_index > &local_dof_indices, FullMatrix< number > &local_matrix, Vector< number > &local_rhs, const bool eliminate_columns)
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition utilities.h:1033
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)