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