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