Reference documentation for deal.II version GIT 6dd81a2340 2022-11-29 05:55:02+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\}}\)
trilinos_sparse_matrix.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2008 - 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
9 // Public License as published by the Free Software Foundation; either
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 
18 
19 #ifdef DEAL_II_WITH_TRILINOS
20 
22 # include <deal.II/base/utilities.h>
23 
30 
32 # include <boost/container/small_vector.hpp>
34 
35 # ifdef DEAL_II_TRILINOS_WITH_EPETRAEXT
36 # include <EpetraExt_MatrixMatrix.h>
37 # endif
38 # include <Epetra_Export.h>
39 # include <Teuchos_RCP.hpp>
40 # include <ml_epetra_utils.h>
41 # include <ml_struct.h>
42 
43 # include <memory>
44 
46 
47 namespace TrilinosWrappers
48 {
49  namespace internal
50  {
51  template <typename VectorType>
52  typename VectorType::value_type *
53  begin(VectorType &V)
54  {
55  return V.begin();
56  }
57 
58  template <typename VectorType>
59  const typename VectorType::value_type *
60  begin(const VectorType &V)
61  {
62  return V.begin();
63  }
64 
65  template <typename VectorType>
66  typename VectorType::value_type *
67  end(VectorType &V)
68  {
69  return V.end();
70  }
71 
72  template <typename VectorType>
73  const typename VectorType::value_type *
74  end(const VectorType &V)
75  {
76  return V.end();
77  }
78 
79  template <>
80  double *
82  {
83  return V.trilinos_vector()[0];
84  }
85 
86  template <>
87  const double *
89  {
90  return V.trilinos_vector()[0];
91  }
92 
93  template <>
94  double *
96  {
97  return V.trilinos_vector()[0] + V.trilinos_vector().MyLength();
98  }
99 
100  template <>
101  const double *
103  {
104  return V.trilinos_vector()[0] + V.trilinos_vector().MyLength();
105  }
106 
107 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
108  template <typename Number>
109  Number *
111  {
112  return V.trilinos_vector().getDataNonConst().get();
113  }
114 
115  template <typename Number>
116  const Number *
118  {
119  return V.trilinos_vector().getData().get();
120  }
121 
122  template <typename Number>
123  Number *
125  {
126  return V.trilinos_vector().getDataNonConst().get() +
127  V.trilinos_vector().getLocalLength();
128  }
129 
130  template <typename Number>
131  const Number *
133  {
134  return V.trilinos_vector().getData().get() +
135  V.trilinos_vector().getLocalLength();
136  }
137 # endif
138  } // namespace internal
139 
140 
141  namespace SparseMatrixIterators
142  {
143  void
145  {
146  // if we are asked to visit the past-the-end line, then simply
147  // release all our caches and go on with life.
148  //
149  // do the same if the row we're supposed to visit is not locally
150  // owned. this is simply going to make non-locally owned rows
151  // look like they're empty
152  if ((this->a_row == matrix->m()) ||
153  (matrix->in_local_range(this->a_row) == false))
154  {
155  colnum_cache.reset();
156  value_cache.reset();
157 
158  return;
159  }
160 
161  // get a representation of the present row
162  int ncols;
164  if (value_cache.get() == nullptr)
165  {
166  value_cache =
167  std::make_shared<std::vector<TrilinosScalar>>(matrix->n());
168  colnum_cache = std::make_shared<std::vector<size_type>>(matrix->n());
169  }
170  else
171  {
172  value_cache->resize(matrix->n());
173  colnum_cache->resize(matrix->n());
174  }
175 
176  int ierr = matrix->trilinos_matrix().ExtractGlobalRowCopy(
177  this->a_row,
178  colnums,
179  ncols,
180  value_cache->data(),
181  reinterpret_cast<TrilinosWrappers::types::int_type *>(
182  colnum_cache->data()));
183  value_cache->resize(ncols);
184  colnum_cache->resize(ncols);
185  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
186 
187  // copy it into our caches if the
188  // line isn't empty. if it is, then
189  // we've done something wrong, since
190  // we shouldn't have initialized an
191  // iterator for an empty line (what
192  // would it point to?)
193  }
194  } // namespace SparseMatrixIterators
195 
196 
197  // The constructor is actually the
198  // only point where we have to check
199  // whether we build a serial or a
200  // parallel Trilinos matrix.
201  // Actually, it does not even matter
202  // how many threads there are, but
203  // only if we use an MPI compiler or
204  // a standard compiler. So, even one
205  // thread on a configuration with
206  // MPI will still get a parallel
207  // interface.
209  : column_space_map(new Epetra_Map(0, 0, Utilities::Trilinos::comm_self()))
210  , matrix(
211  new Epetra_FECrsMatrix(View, *column_space_map, *column_space_map, 0))
212  , last_action(Zero)
213  , compressed(true)
214  {
215  matrix->FillComplete();
216  }
217 
218 
219 
221  const size_type n,
222  const unsigned int n_max_entries_per_row)
223  : column_space_map(
224  new Epetra_Map(static_cast<TrilinosWrappers::types::int_type>(n),
225  0,
226  Utilities::Trilinos::comm_self()))
227  ,
228 
229  // on one processor only, we know how the
230  // columns of the matrix will be
231  // distributed (everything on one
232  // processor), so we can hand in this
233  // information to the constructor. we
234  // can't do so in parallel, where the
235  // information from columns is only
236  // available when entries have been added
237  matrix(new Epetra_FECrsMatrix(
238  Copy,
239  Epetra_Map(static_cast<TrilinosWrappers::types::int_type>(m),
240  0,
241  Utilities::Trilinos::comm_self()),
242  *column_space_map,
243  n_max_entries_per_row,
244  false))
245  , last_action(Zero)
246  , compressed(false)
247  {}
248 
249 
250 
252  const size_type n,
253  const std::vector<unsigned int> &n_entries_per_row)
254  : column_space_map(
255  new Epetra_Map(static_cast<TrilinosWrappers::types::int_type>(n),
256  0,
257  Utilities::Trilinos::comm_self()))
258  , matrix(new Epetra_FECrsMatrix(
259  Copy,
260  Epetra_Map(static_cast<TrilinosWrappers::types::int_type>(m),
261  0,
262  Utilities::Trilinos::comm_self()),
263  *column_space_map,
264  reinterpret_cast<int *>(
265  const_cast<unsigned int *>(n_entries_per_row.data())),
266  false))
267  , last_action(Zero)
268  , compressed(false)
269  {}
270 
271 
272 
273  SparseMatrix::SparseMatrix(const IndexSet & parallel_partitioning,
274  const MPI_Comm & communicator,
275  const unsigned int n_max_entries_per_row)
276  : column_space_map(new Epetra_Map(
277  parallel_partitioning.make_trilinos_map(communicator, false)))
278  , matrix(new Epetra_FECrsMatrix(Copy,
279  *column_space_map,
280  n_max_entries_per_row,
281  false))
282  , last_action(Zero)
283  , compressed(false)
284  {}
285 
286 
287 
288  SparseMatrix::SparseMatrix(const IndexSet &parallel_partitioning,
289  const MPI_Comm &communicator,
290  const std::vector<unsigned int> &n_entries_per_row)
291  : column_space_map(new Epetra_Map(
292  parallel_partitioning.make_trilinos_map(communicator, false)))
293  , matrix(new Epetra_FECrsMatrix(Copy,
294  *column_space_map,
295  reinterpret_cast<int *>(
296  const_cast<unsigned int *>(
297  n_entries_per_row.data())),
298  false))
299  , last_action(Zero)
300  , compressed(false)
301  {}
302 
303 
304 
305  SparseMatrix::SparseMatrix(const IndexSet &row_parallel_partitioning,
306  const IndexSet &col_parallel_partitioning,
307  const MPI_Comm &communicator,
308  const size_type n_max_entries_per_row)
309  : column_space_map(new Epetra_Map(
310  col_parallel_partitioning.make_trilinos_map(communicator, false)))
311  , matrix(new Epetra_FECrsMatrix(
312  Copy,
313  row_parallel_partitioning.make_trilinos_map(communicator, false),
314  n_max_entries_per_row,
315  false))
316  , last_action(Zero)
317  , compressed(false)
318  {}
319 
320 
321 
322  SparseMatrix::SparseMatrix(const IndexSet &row_parallel_partitioning,
323  const IndexSet &col_parallel_partitioning,
324  const MPI_Comm &communicator,
325  const std::vector<unsigned int> &n_entries_per_row)
326  : column_space_map(new Epetra_Map(
327  col_parallel_partitioning.make_trilinos_map(communicator, false)))
328  , matrix(new Epetra_FECrsMatrix(
329  Copy,
330  row_parallel_partitioning.make_trilinos_map(communicator, false),
331  reinterpret_cast<int *>(
332  const_cast<unsigned int *>(n_entries_per_row.data())),
333  false))
334  , last_action(Zero)
335  , compressed(false)
336  {}
337 
338 
339 
341  : column_space_map(new Epetra_Map(sparsity_pattern.domain_partitioner()))
342  , matrix(
343  new Epetra_FECrsMatrix(Copy,
344  sparsity_pattern.trilinos_sparsity_pattern(),
345  false))
346  , last_action(Zero)
347  , compressed(true)
348  {
349  Assert(sparsity_pattern.trilinos_sparsity_pattern().Filled() == true,
350  ExcMessage(
351  "The Trilinos sparsity pattern has not been compressed."));
353  }
354 
355 
356 
358  : column_space_map(std::move(other.column_space_map))
359  , matrix(std::move(other.matrix))
360  , nonlocal_matrix(std::move(other.nonlocal_matrix))
361  , nonlocal_matrix_exporter(std::move(other.nonlocal_matrix_exporter))
362  , last_action(other.last_action)
363  , compressed(other.compressed)
364  {
365  other.last_action = Zero;
366  other.compressed = false;
367  }
368 
369 
370 
371  void
373  {
374  if (this == &rhs)
375  return;
376 
377  nonlocal_matrix.reset();
378  nonlocal_matrix_exporter.reset();
379 
380  // check whether we need to update the whole matrix layout (we have
381  // different maps or if we detect a row where the columns of the two
382  // matrices do not match)
383  bool needs_deep_copy =
384  !matrix->RowMap().SameAs(rhs.matrix->RowMap()) ||
385  !matrix->ColMap().SameAs(rhs.matrix->ColMap()) ||
386  !matrix->DomainMap().SameAs(rhs.matrix->DomainMap()) ||
388  if (!needs_deep_copy)
389  {
390  // Try to copy all the rows of the matrix one by one. In case of error
391  // (i.e., the column indices are different), we need to abort and blow
392  // away the matrix.
393  for (const auto row : locally_owned_range_indices())
394  {
395  const int row_local = matrix->RowMap().LID(
396  static_cast<TrilinosWrappers::types::int_type>(row));
397  Assert((row_local >= 0), ExcAccessToNonlocalRow(row));
398 
399  int n_entries, rhs_n_entries;
400  TrilinosScalar *value_ptr, *rhs_value_ptr;
401  int * index_ptr, *rhs_index_ptr;
402  int ierr = rhs.matrix->ExtractMyRowView(row_local,
403  rhs_n_entries,
404  rhs_value_ptr,
405  rhs_index_ptr);
406  (void)ierr;
407  Assert(ierr == 0, ExcTrilinosError(ierr));
408 
409  ierr = matrix->ExtractMyRowView(row_local,
410  n_entries,
411  value_ptr,
412  index_ptr);
413  Assert(ierr == 0, ExcTrilinosError(ierr));
414 
415  if (n_entries != rhs_n_entries ||
416  std::memcmp(static_cast<void *>(index_ptr),
417  static_cast<void *>(rhs_index_ptr),
418  sizeof(int) * n_entries) != 0)
419  {
420  needs_deep_copy = true;
421  break;
422  }
423 
424  for (int i = 0; i < n_entries; ++i)
425  value_ptr[i] = rhs_value_ptr[i];
426  }
427  }
428 
429  if (needs_deep_copy)
430  {
432  std::make_unique<Epetra_Map>(rhs.trilinos_matrix().DomainMap());
433 
434  // release memory before reallocation
435  matrix = std::make_unique<Epetra_FECrsMatrix>(*rhs.matrix);
436 
437  matrix->FillComplete(*column_space_map, matrix->RowMap());
438  }
439 
440  if (rhs.nonlocal_matrix.get() != nullptr)
442  std::make_unique<Epetra_CrsMatrix>(Copy, rhs.nonlocal_matrix->Graph());
443  }
444 
445 
446 
447  namespace
448  {
450 
451  template <typename SparsityPatternType>
452  void
453  reinit_matrix(const IndexSet & row_parallel_partitioning,
454  const IndexSet & column_parallel_partitioning,
455  const SparsityPatternType & sparsity_pattern,
456  const bool exchange_data,
457  const MPI_Comm & communicator,
458  std::unique_ptr<Epetra_Map> &column_space_map,
459  std::unique_ptr<Epetra_FECrsMatrix> &matrix,
460  std::unique_ptr<Epetra_CrsMatrix> & nonlocal_matrix,
461  std::unique_ptr<Epetra_Export> & nonlocal_matrix_exporter)
462  {
463  // release memory before reallocation
464  matrix.reset();
465  nonlocal_matrix.reset();
466  nonlocal_matrix_exporter.reset();
467 
468  column_space_map = std::make_unique<Epetra_Map>(
469  column_parallel_partitioning.make_trilinos_map(communicator, false));
470 
471  if (column_space_map->Comm().MyPID() == 0)
472  {
473  AssertDimension(sparsity_pattern.n_rows(),
474  row_parallel_partitioning.size());
475  AssertDimension(sparsity_pattern.n_cols(),
476  column_parallel_partitioning.size());
477  }
478 
479  Epetra_Map row_space_map =
480  row_parallel_partitioning.make_trilinos_map(communicator, false);
481 
482  // if we want to exchange data, build a usual Trilinos sparsity pattern
483  // and let that handle the exchange. otherwise, manually create a
484  // CrsGraph, which consumes considerably less memory because it can set
485  // correct number of indices right from the start
486  if (exchange_data)
487  {
488  SparsityPattern trilinos_sparsity;
489  trilinos_sparsity.reinit(row_parallel_partitioning,
490  column_parallel_partitioning,
491  sparsity_pattern,
492  communicator,
493  exchange_data);
494  matrix = std::make_unique<Epetra_FECrsMatrix>(
495  Copy, trilinos_sparsity.trilinos_sparsity_pattern(), false);
496 
497  return;
498  }
499 
500  const size_type first_row = TrilinosWrappers::min_my_gid(row_space_map),
501  last_row =
502  TrilinosWrappers::max_my_gid(row_space_map) + 1;
503  std::vector<int> n_entries_per_row(last_row - first_row);
504 
505  for (size_type row = first_row; row < last_row; ++row)
506  n_entries_per_row[row - first_row] = sparsity_pattern.row_length(row);
507 
508  // The deal.II notation of a Sparsity pattern corresponds to the Epetra
509  // concept of a Graph. Hence, we generate a graph by copying the
510  // sparsity pattern into it, and then build up the matrix from the
511  // graph. This is considerable faster than directly filling elements
512  // into the matrix. Moreover, it consumes less memory, since the
513  // internal reordering is done on ints only, and we can leave the
514  // doubles aside.
515 
516  // for more than one processor, need to specify only row map first and
517  // let the matrix entries decide about the column map (which says which
518  // columns are present in the matrix, not to be confused with the
519  // col_map that tells how the domain dofs of the matrix will be
520  // distributed). for only one processor, we can directly assign the
521  // columns as well. Compare this with bug # 4123 in the Sandia Bugzilla.
522  std::unique_ptr<Epetra_CrsGraph> graph;
523  if (row_space_map.Comm().NumProc() > 1)
524  graph = std::make_unique<Epetra_CrsGraph>(Copy,
525  row_space_map,
526  n_entries_per_row.data(),
527  true);
528  else
529  graph = std::make_unique<Epetra_CrsGraph>(Copy,
530  row_space_map,
531  *column_space_map,
532  n_entries_per_row.data(),
533  true);
534 
535  // This functions assumes that the sparsity pattern sits on all
536  // processors (completely). The parallel version uses an Epetra graph
537  // that is already distributed.
538 
539  // now insert the indices
540  std::vector<TrilinosWrappers::types::int_type> row_indices;
541 
542  for (size_type row = first_row; row < last_row; ++row)
543  {
544  const int row_length = sparsity_pattern.row_length(row);
545  if (row_length == 0)
546  continue;
547 
548  row_indices.resize(row_length, -1);
549  {
550  typename SparsityPatternType::iterator p =
551  sparsity_pattern.begin(row);
552  for (size_type col = 0; p != sparsity_pattern.end(row); ++p, ++col)
553  row_indices[col] = p->column();
554  }
555  graph->Epetra_CrsGraph::InsertGlobalIndices(row,
556  row_length,
557  row_indices.data());
558  }
559 
560  // Eventually, optimize the graph structure (sort indices, make memory
561  // contiguous, etc). note that the documentation of the function indeed
562  // states that we first need to provide the column (domain) map and then
563  // the row (range) map
564  graph->FillComplete(*column_space_map, row_space_map);
565  graph->OptimizeStorage();
566 
567  // check whether we got the number of columns right.
568  AssertDimension(sparsity_pattern.n_cols(),
570  (void)n_global_cols;
571 
572  // And now finally generate the matrix.
573  matrix = std::make_unique<Epetra_FECrsMatrix>(Copy, *graph, false);
574  }
575 
576 
577 
578  // for the non-local graph, we need to circumvent the problem that some
579  // processors will not add into the non-local graph at all: We do not want
580  // to insert dummy elements on >5000 processors because that gets very
581  // slow. Thus, we set a flag in Epetra_CrsGraph that sets the correct
582  // flag. Since it is protected, we need to expose this information by
583  // deriving a class from Epetra_CrsGraph for the purpose of creating the
584  // data structure
585  class Epetra_CrsGraphMod : public Epetra_CrsGraph
586  {
587  public:
588  Epetra_CrsGraphMod(const Epetra_Map &row_map,
589  const int * n_entries_per_row)
590  : Epetra_CrsGraph(Copy, row_map, n_entries_per_row, true)
591  {}
592 
593  void
594  SetIndicesAreGlobal()
595  {
596  this->Epetra_CrsGraph::SetIndicesAreGlobal(true);
597  }
598  };
599 
600 
601 
602  // specialization for DynamicSparsityPattern which can provide us with
603  // more information about the non-locally owned rows
604  template <>
605  void
606  reinit_matrix(const IndexSet & row_parallel_partitioning,
607  const IndexSet & column_parallel_partitioning,
608  const DynamicSparsityPattern &sparsity_pattern,
609  const bool exchange_data,
610  const MPI_Comm & communicator,
611  std::unique_ptr<Epetra_Map> & column_space_map,
612  std::unique_ptr<Epetra_FECrsMatrix> &matrix,
613  std::unique_ptr<Epetra_CrsMatrix> & nonlocal_matrix,
614  std::unique_ptr<Epetra_Export> & nonlocal_matrix_exporter)
615  {
616  matrix.reset();
617  nonlocal_matrix.reset();
618  nonlocal_matrix_exporter.reset();
619 
620  column_space_map = std::make_unique<Epetra_Map>(
621  column_parallel_partitioning.make_trilinos_map(communicator, false));
622 
623  AssertDimension(sparsity_pattern.n_rows(),
624  row_parallel_partitioning.size());
625  AssertDimension(sparsity_pattern.n_cols(),
626  column_parallel_partitioning.size());
627 
628  Epetra_Map row_space_map =
629  row_parallel_partitioning.make_trilinos_map(communicator, false);
630 
631  IndexSet relevant_rows(sparsity_pattern.row_index_set());
632  // serial case
633  if (relevant_rows.size() == 0)
634  {
635  relevant_rows.set_size(
636  TrilinosWrappers::n_global_elements(row_space_map));
637  relevant_rows.add_range(
638  0, TrilinosWrappers::n_global_elements(row_space_map));
639  }
640  relevant_rows.compress();
641  Assert(relevant_rows.n_elements() >=
642  static_cast<unsigned int>(row_space_map.NumMyElements()),
643  ExcMessage(
644  "Locally relevant rows of sparsity pattern must contain "
645  "all locally owned rows"));
646 
647  // check whether the relevant rows correspond to exactly the same map as
648  // the owned rows. In that case, do not create the nonlocal graph and
649  // fill the columns by demand
650  const bool have_ghost_rows = [&]() {
651  std::vector<::types::global_dof_index> indices;
652  relevant_rows.fill_index_vector(indices);
653  Epetra_Map relevant_map(
655  TrilinosWrappers::types::int_type(relevant_rows.n_elements()),
656  (indices.empty() ?
657  nullptr :
658  reinterpret_cast<TrilinosWrappers::types::int_type *>(
659  indices.data())),
660  0,
661  row_space_map.Comm());
662  return !relevant_map.SameAs(row_space_map);
663  }();
664 
665  const unsigned int n_rows = relevant_rows.n_elements();
666  std::vector<TrilinosWrappers::types::int_type> ghost_rows;
667  std::vector<int> n_entries_per_row(row_space_map.NumMyElements());
668  std::vector<int> n_entries_per_ghost_row;
669  for (unsigned int i = 0, own = 0; i < n_rows; ++i)
670  {
671  const TrilinosWrappers::types::int_type global_row =
672  relevant_rows.nth_index_in_set(i);
673  if (row_space_map.MyGID(global_row))
674  n_entries_per_row[own++] = sparsity_pattern.row_length(global_row);
675  else if (sparsity_pattern.row_length(global_row) > 0)
676  {
677  ghost_rows.push_back(global_row);
678  n_entries_per_ghost_row.push_back(
679  sparsity_pattern.row_length(global_row));
680  }
681  }
682 
683  Epetra_Map off_processor_map(-1,
684  ghost_rows.size(),
685  (ghost_rows.size() > 0) ?
686  (ghost_rows.data()) :
687  nullptr,
688  0,
689  row_space_map.Comm());
690 
691  std::unique_ptr<Epetra_CrsGraph> graph;
692  std::unique_ptr<Epetra_CrsGraphMod> nonlocal_graph;
693  if (row_space_map.Comm().NumProc() > 1)
694  {
695  graph =
696  std::make_unique<Epetra_CrsGraph>(Copy,
697  row_space_map,
698  (n_entries_per_row.size() > 0) ?
699  (n_entries_per_row.data()) :
700  nullptr,
701  exchange_data ? false : true);
702  if (have_ghost_rows == true)
703  nonlocal_graph = std::make_unique<Epetra_CrsGraphMod>(
704  off_processor_map, n_entries_per_ghost_row.data());
705  }
706  else
707  graph =
708  std::make_unique<Epetra_CrsGraph>(Copy,
709  row_space_map,
710  *column_space_map,
711  (n_entries_per_row.size() > 0) ?
712  (n_entries_per_row.data()) :
713  nullptr,
714  true);
715 
716  // now insert the indices, select between the right matrix
717  std::vector<TrilinosWrappers::types::int_type> row_indices;
718 
719  for (unsigned int i = 0; i < n_rows; ++i)
720  {
721  const TrilinosWrappers::types::int_type global_row =
722  relevant_rows.nth_index_in_set(i);
723  const int row_length = sparsity_pattern.row_length(global_row);
724  if (row_length == 0)
725  continue;
726 
727  row_indices.resize(row_length, -1);
728  for (int col = 0; col < row_length; ++col)
729  row_indices[col] = sparsity_pattern.column_number(global_row, col);
730 
731  if (row_space_map.MyGID(global_row))
732  graph->InsertGlobalIndices(global_row,
733  row_length,
734  row_indices.data());
735  else
736  {
737  Assert(nonlocal_graph.get() != nullptr, ExcInternalError());
738  nonlocal_graph->InsertGlobalIndices(global_row,
739  row_length,
740  row_indices.data());
741  }
742  }
743 
744  // finalize nonlocal graph and create nonlocal matrix
745  if (nonlocal_graph.get() != nullptr)
746  {
747  // must make sure the IndicesAreGlobal flag is set on all processors
748  // because some processors might not call InsertGlobalIndices (and
749  // we do not want to insert dummy indices on all processors for
750  // large-scale simulations due to the bad impact on performance)
751  nonlocal_graph->SetIndicesAreGlobal();
752  Assert(nonlocal_graph->IndicesAreGlobal() == true,
753  ExcInternalError());
754  nonlocal_graph->FillComplete(*column_space_map, row_space_map);
755  nonlocal_graph->OptimizeStorage();
756 
757  // insert data from nonlocal graph into the final sparsity pattern
758  if (exchange_data)
759  {
760  Epetra_Export exporter(nonlocal_graph->RowMap(), row_space_map);
761  int ierr = graph->Export(*nonlocal_graph, exporter, Add);
762  (void)ierr;
763  Assert(ierr == 0, ExcTrilinosError(ierr));
764  }
765 
766  nonlocal_matrix =
767  std::make_unique<Epetra_CrsMatrix>(Copy, *nonlocal_graph);
768  }
769 
770  graph->FillComplete(*column_space_map, row_space_map);
771  graph->OptimizeStorage();
772 
773  AssertDimension(sparsity_pattern.n_cols(),
775 
776  matrix = std::make_unique<Epetra_FECrsMatrix>(Copy, *graph, false);
777  }
778  } // namespace
779 
780 
781 
782  template <typename SparsityPatternType>
783  void
784  SparseMatrix::reinit(const SparsityPatternType &sparsity_pattern)
785  {
786  reinit_matrix(complete_index_set(sparsity_pattern.n_rows()),
787  complete_index_set(sparsity_pattern.n_cols()),
788  sparsity_pattern,
789  false,
790  MPI_COMM_SELF,
792  matrix,
795  }
796 
797 
798 
799  template <typename SparsityPatternType>
800  inline std::enable_if_t<
801  !std::is_same<SparsityPatternType, ::SparseMatrix<double>>::value>
802  SparseMatrix::reinit(const IndexSet & row_parallel_partitioning,
803  const IndexSet & col_parallel_partitioning,
804  const SparsityPatternType &sparsity_pattern,
805  const MPI_Comm & communicator,
806  const bool exchange_data)
807  {
808  reinit_matrix(row_parallel_partitioning,
809  col_parallel_partitioning,
810  sparsity_pattern,
811  exchange_data,
812  communicator,
814  matrix,
817 
818  // In the end, the matrix needs to be compressed in order to be really
819  // ready.
820  last_action = Zero;
822  }
823 
824 
825 
826  void
827  SparseMatrix::reinit(const SparsityPattern &sparsity_pattern)
828  {
829  matrix.reset();
830  nonlocal_matrix_exporter.reset();
831 
832  // reinit with a (parallel) Trilinos sparsity pattern.
834  std::make_unique<Epetra_Map>(sparsity_pattern.domain_partitioner());
835  matrix = std::make_unique<Epetra_FECrsMatrix>(
836  Copy, sparsity_pattern.trilinos_sparsity_pattern(), false);
837 
838  if (sparsity_pattern.nonlocal_graph.get() != nullptr)
840  std::make_unique<Epetra_CrsMatrix>(Copy,
841  *sparsity_pattern.nonlocal_graph);
842  else
843  nonlocal_matrix.reset();
844 
845  last_action = Zero;
847  }
848 
849 
850 
851  void
852  SparseMatrix::reinit(const SparseMatrix &sparse_matrix)
853  {
854  if (this == &sparse_matrix)
855  return;
856 
858  std::make_unique<Epetra_Map>(sparse_matrix.trilinos_matrix().DomainMap());
859  matrix.reset();
860  nonlocal_matrix_exporter.reset();
861  matrix = std::make_unique<Epetra_FECrsMatrix>(
862  Copy, sparse_matrix.trilinos_sparsity_pattern(), false);
863 
864  if (sparse_matrix.nonlocal_matrix != nullptr)
865  nonlocal_matrix = std::make_unique<Epetra_CrsMatrix>(
866  Copy, sparse_matrix.nonlocal_matrix->Graph());
867  else
868  nonlocal_matrix.reset();
869 
870  last_action = Zero;
872  }
873 
874 
875 
876  template <typename number>
877  inline void
879  const IndexSet & row_parallel_partitioning,
880  const IndexSet & col_parallel_partitioning,
881  const ::SparseMatrix<number> &dealii_sparse_matrix,
882  const MPI_Comm & communicator,
883  const double drop_tolerance,
884  const bool copy_values,
885  const ::SparsityPattern * use_this_sparsity)
886  {
887  if (copy_values == false)
888  {
889  // in case we do not copy values, just
890  // call the other function.
891  if (use_this_sparsity == nullptr)
892  reinit(row_parallel_partitioning,
893  col_parallel_partitioning,
894  dealii_sparse_matrix.get_sparsity_pattern(),
895  communicator,
896  false);
897  else
898  reinit(row_parallel_partitioning,
899  col_parallel_partitioning,
900  *use_this_sparsity,
901  communicator,
902  false);
903  return;
904  }
905 
906  const size_type n_rows = dealii_sparse_matrix.m();
907 
908  AssertDimension(row_parallel_partitioning.size(), n_rows);
909  AssertDimension(col_parallel_partitioning.size(), dealii_sparse_matrix.n());
910 
911  const ::SparsityPattern &sparsity_pattern =
912  (use_this_sparsity != nullptr) ?
913  *use_this_sparsity :
914  dealii_sparse_matrix.get_sparsity_pattern();
915 
916  if (matrix.get() == nullptr || m() != n_rows ||
917  n_nonzero_elements() != sparsity_pattern.n_nonzero_elements())
918  {
919  reinit(row_parallel_partitioning,
920  col_parallel_partitioning,
921  sparsity_pattern,
922  communicator,
923  false);
924  }
925 
926  // fill the values. the same as above: go through all rows of the
927  // matrix, and then all columns. since the sparsity patterns of the
928  // input matrix and the specified sparsity pattern might be different,
929  // need to go through the row for both these sparsity structures
930  // simultaneously in order to really set the correct values.
931  size_type maximum_row_length = matrix->MaxNumEntries();
932  std::vector<size_type> row_indices(maximum_row_length);
933  std::vector<TrilinosScalar> values(maximum_row_length);
934 
935  for (size_type row = 0; row < n_rows; ++row)
936  // see if the row is locally stored on this processor
937  if (row_parallel_partitioning.is_element(row) == true)
938  {
939  ::SparsityPattern::iterator select_index =
940  sparsity_pattern.begin(row);
941  typename ::SparseMatrix<number>::const_iterator it =
942  dealii_sparse_matrix.begin(row);
943  size_type col = 0;
944  if (sparsity_pattern.n_rows() == sparsity_pattern.n_cols())
945  {
946  // optimized diagonal
947  AssertDimension(it->column(), row);
948  if (std::fabs(it->value()) > drop_tolerance)
949  {
950  values[col] = it->value();
951  row_indices[col++] = it->column();
952  }
953  ++select_index;
954  ++it;
955  }
956 
957  while (it != dealii_sparse_matrix.end(row) &&
958  select_index != sparsity_pattern.end(row))
959  {
960  while (select_index->column() < it->column() &&
961  select_index != sparsity_pattern.end(row))
962  ++select_index;
963  while (it->column() < select_index->column() &&
964  it != dealii_sparse_matrix.end(row))
965  ++it;
966 
967  if (it == dealii_sparse_matrix.end(row))
968  break;
969  if (std::fabs(it->value()) > drop_tolerance)
970  {
971  values[col] = it->value();
972  row_indices[col++] = it->column();
973  }
974  ++select_index;
975  ++it;
976  }
977  set(row,
978  col,
979  reinterpret_cast<size_type *>(row_indices.data()),
980  values.data(),
981  false);
982  }
984  }
985 
986 
987 
988  template <typename number>
989  void
991  const ::SparseMatrix<number> &dealii_sparse_matrix,
992  const double drop_tolerance,
993  const bool copy_values,
994  const ::SparsityPattern * use_this_sparsity)
995  {
996  reinit(complete_index_set(dealii_sparse_matrix.m()),
997  complete_index_set(dealii_sparse_matrix.n()),
998  dealii_sparse_matrix,
999  MPI_COMM_SELF,
1000  drop_tolerance,
1001  copy_values,
1002  use_this_sparsity);
1003  }
1004 
1005 
1006 
1007  void
1008  SparseMatrix::reinit(const Epetra_CrsMatrix &input_matrix,
1009  const bool copy_values)
1010  {
1011  Assert(input_matrix.Filled() == true,
1012  ExcMessage("Input CrsMatrix has not called FillComplete()!"));
1013 
1014  column_space_map = std::make_unique<Epetra_Map>(input_matrix.DomainMap());
1015 
1016  const Epetra_CrsGraph *graph = &input_matrix.Graph();
1017 
1018  nonlocal_matrix.reset();
1019  nonlocal_matrix_exporter.reset();
1020  matrix.reset();
1021  matrix = std::make_unique<Epetra_FECrsMatrix>(Copy, *graph, false);
1022 
1023  matrix->FillComplete(*column_space_map, input_matrix.RangeMap(), true);
1024 
1025  if (copy_values == true)
1026  {
1027  // point to the first data entry in the two
1028  // matrices and copy the content
1029  const TrilinosScalar *in_values = input_matrix[0];
1030  TrilinosScalar * values = (*matrix)[0];
1031  const size_type my_nonzeros = input_matrix.NumMyNonzeros();
1032  std::memcpy(values, in_values, my_nonzeros * sizeof(TrilinosScalar));
1033  }
1034 
1035  last_action = Zero;
1037  }
1038 
1039 
1040 
1041  void
1043  {
1044  Epetra_CombineMode mode = last_action;
1045  if (last_action == Zero)
1046  {
1047  if ((operation == ::VectorOperation::add) ||
1048  (operation == ::VectorOperation::unknown))
1049  mode = Add;
1050  else if (operation == ::VectorOperation::insert)
1051  mode = Insert;
1052  else
1053  Assert(
1054  false,
1055  ExcMessage(
1056  "compress() can only be called with VectorOperation add, insert, or unknown"));
1057  }
1058  else
1059  {
1060  Assert(((last_action == Add) &&
1061  (operation != ::VectorOperation::insert)) ||
1062  ((last_action == Insert) &&
1063  (operation != ::VectorOperation::add)),
1064  ExcMessage("Operation and argument to compress() do not match"));
1065  }
1066 
1067  // flush buffers
1068  int ierr;
1069  if (nonlocal_matrix.get() != nullptr && mode == Add)
1070  {
1071  // do only export in case of an add() operation, otherwise the owning
1072  // processor must have set the correct entry
1073  nonlocal_matrix->FillComplete(*column_space_map, matrix->RowMap());
1074  if (nonlocal_matrix_exporter.get() == nullptr)
1076  std::make_unique<Epetra_Export>(nonlocal_matrix->RowMap(),
1077  matrix->RowMap());
1078  ierr =
1080  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1081  ierr = matrix->FillComplete(*column_space_map, matrix->RowMap());
1082  nonlocal_matrix->PutScalar(0);
1083  }
1084  else
1085  ierr =
1086  matrix->GlobalAssemble(*column_space_map, matrix->RowMap(), true, mode);
1087 
1088  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1089 
1090  ierr = matrix->OptimizeStorage();
1091  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1092 
1093  last_action = Zero;
1094 
1095  compressed = true;
1096  }
1097 
1098 
1099 
1100  void
1102  {
1103  // When we clear the matrix, reset
1104  // the pointer and generate an
1105  // empty matrix.
1107  std::make_unique<Epetra_Map>(0, 0, Utilities::Trilinos::comm_self());
1108  matrix = std::make_unique<Epetra_FECrsMatrix>(View, *column_space_map, 0);
1109  nonlocal_matrix.reset();
1110  nonlocal_matrix_exporter.reset();
1111 
1112  matrix->FillComplete();
1113 
1114  compressed = true;
1115  }
1116 
1117 
1118 
1119  void
1121  const TrilinosScalar new_diag_value)
1122  {
1123  Assert(matrix->Filled() == true, ExcMatrixNotCompressed());
1124 
1125  // Only do this on the rows owned
1126  // locally on this processor.
1127  int local_row =
1128  matrix->LRID(static_cast<TrilinosWrappers::types::int_type>(row));
1129  if (local_row >= 0)
1130  {
1132  int * col_indices;
1133  int num_entries;
1134  const int ierr =
1135  matrix->ExtractMyRowView(local_row, num_entries, values, col_indices);
1136  (void)ierr;
1137 
1138  Assert(ierr == 0, ExcTrilinosError(ierr));
1139 
1140  const std::ptrdiff_t diag_index =
1141  std::find(col_indices, col_indices + num_entries, local_row) -
1142  col_indices;
1143 
1144  for (TrilinosWrappers::types::int_type j = 0; j < num_entries; ++j)
1145  if (diag_index != j || new_diag_value == 0)
1146  values[j] = 0.;
1147 
1148  if (diag_index != num_entries)
1149  values[diag_index] = new_diag_value;
1150  }
1151  }
1152 
1153 
1154 
1155  void
1156  SparseMatrix::clear_rows(const std::vector<size_type> &rows,
1157  const TrilinosScalar new_diag_value)
1158  {
1159  for (const auto row : rows)
1160  clear_row(row, new_diag_value);
1161  }
1162 
1163 
1164 
1167  {
1168  // Extract local indices in
1169  // the matrix.
1170  int trilinos_i =
1171  matrix->LRID(static_cast<TrilinosWrappers::types::int_type>(i)),
1172  trilinos_j =
1173  matrix->LCID(static_cast<TrilinosWrappers::types::int_type>(j));
1174  TrilinosScalar value = 0.;
1175 
1176  // If the data is not on the
1177  // present processor, we throw
1178  // an exception. This is one of
1179  // the two tiny differences to
1180  // the el(i,j) call, which does
1181  // not throw any assertions.
1182  if (trilinos_i == -1)
1183  {
1184  Assert(false,
1186  i, j, local_range().first, local_range().second - 1));
1187  }
1188  else
1189  {
1190  // Check whether the matrix has
1191  // already been transformed to local
1192  // indices.
1193  Assert(matrix->Filled(), ExcMatrixNotCompressed());
1194 
1195  // Prepare pointers for extraction
1196  // of a view of the row.
1197  int nnz_present = matrix->NumMyEntries(trilinos_i);
1198  int nnz_extracted;
1199  int * col_indices;
1201 
1202  // Generate the view and make
1203  // sure that we have not generated
1204  // an error.
1205  // TODO Check that col_indices are int and not long long
1206  int ierr = matrix->ExtractMyRowView(trilinos_i,
1207  nnz_extracted,
1208  values,
1209  col_indices);
1210  (void)ierr;
1211  Assert(ierr == 0, ExcTrilinosError(ierr));
1212 
1213  Assert(nnz_present == nnz_extracted,
1214  ExcDimensionMismatch(nnz_present, nnz_extracted));
1215 
1216  // Search the index where we
1217  // look for the value, and then
1218  // finally get it.
1219  const std::ptrdiff_t local_col_index =
1220  std::find(col_indices, col_indices + nnz_present, trilinos_j) -
1221  col_indices;
1222 
1223  // This is actually the only
1224  // difference to the el(i,j)
1225  // function, which means that
1226  // we throw an exception in
1227  // this case instead of just
1228  // returning zero for an
1229  // element that is not present
1230  // in the sparsity pattern.
1231  if (local_col_index == nnz_present)
1232  {
1233  Assert(false, ExcInvalidIndex(i, j));
1234  }
1235  else
1236  value = values[local_col_index];
1237  }
1238 
1239  return value;
1240  }
1241 
1242 
1243 
1245  SparseMatrix::el(const size_type i, const size_type j) const
1246  {
1247  // Extract local indices in
1248  // the matrix.
1249  int trilinos_i =
1250  matrix->LRID(static_cast<TrilinosWrappers::types::int_type>(i)),
1251  trilinos_j =
1252  matrix->LCID(static_cast<TrilinosWrappers::types::int_type>(j));
1253  TrilinosScalar value = 0.;
1254 
1255  // If the data is not on the
1256  // present processor, we can't
1257  // continue. Just print out zero
1258  // as discussed in the
1259  // documentation of this
1260  // function. if you want error
1261  // checking, use operator().
1262  if ((trilinos_i == -1) || (trilinos_j == -1))
1263  return 0.;
1264  else
1265  {
1266  // Check whether the matrix
1267  // already is transformed to
1268  // local indices.
1269  Assert(matrix->Filled(), ExcMatrixNotCompressed());
1270 
1271  // Prepare pointers for extraction
1272  // of a view of the row.
1273  int nnz_present = matrix->NumMyEntries(trilinos_i);
1274  int nnz_extracted;
1275  int * col_indices;
1277 
1278  // Generate the view and make
1279  // sure that we have not generated
1280  // an error.
1281  int ierr = matrix->ExtractMyRowView(trilinos_i,
1282  nnz_extracted,
1283  values,
1284  col_indices);
1285  (void)ierr;
1286  Assert(ierr == 0, ExcTrilinosError(ierr));
1287 
1288  Assert(nnz_present == nnz_extracted,
1289  ExcDimensionMismatch(nnz_present, nnz_extracted));
1290 
1291  // Search the index where we
1292  // look for the value, and then
1293  // finally get it.
1294  const std::ptrdiff_t local_col_index =
1295  std::find(col_indices, col_indices + nnz_present, trilinos_j) -
1296  col_indices;
1297 
1298  // This is actually the only
1299  // difference to the () function
1300  // querying (i,j), where we throw an
1301  // exception instead of just
1302  // returning zero for an element
1303  // that is not present in the
1304  // sparsity pattern.
1305  if (local_col_index == nnz_present)
1306  value = 0;
1307  else
1308  value = values[local_col_index];
1309  }
1310 
1311  return value;
1312  }
1313 
1314 
1315 
1318  {
1319  Assert(m() == n(), ExcNotQuadratic());
1320 
1321 # ifdef DEBUG
1322  // use operator() in debug mode because
1323  // it checks if this is a valid element
1324  // (in parallel)
1325  return operator()(i, i);
1326 # else
1327  // Trilinos doesn't seem to have a
1328  // more efficient way to access the
1329  // diagonal than by just using the
1330  // standard el(i,j) function.
1331  return el(i, i);
1332 # endif
1333  }
1334 
1335 
1336 
1337  unsigned int
1339  {
1340  Assert(row < m(), ExcInternalError());
1341 
1342  // get a representation of the
1343  // present row
1344  int ncols = -1;
1345  int local_row =
1346  matrix->LRID(static_cast<TrilinosWrappers::types::int_type>(row));
1347  Assert((local_row >= 0), ExcAccessToNonlocalRow(row));
1348 
1349  // on the processor who owns this
1350  // row, we'll have a non-negative
1351  // value.
1352  if (local_row >= 0)
1353  {
1354  int ierr = matrix->NumMyRowEntries(local_row, ncols);
1355  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1356  }
1357 
1358  return static_cast<unsigned int>(ncols);
1359  }
1360 
1361 
1362 
1363  void
1364  SparseMatrix::set(const std::vector<size_type> & row_indices,
1365  const std::vector<size_type> & col_indices,
1367  const bool elide_zero_values)
1368  {
1369  Assert(row_indices.size() == values.m(),
1370  ExcDimensionMismatch(row_indices.size(), values.m()));
1371  Assert(col_indices.size() == values.n(),
1372  ExcDimensionMismatch(col_indices.size(), values.n()));
1373 
1374  for (size_type i = 0; i < row_indices.size(); ++i)
1375  set(row_indices[i],
1376  col_indices.size(),
1377  col_indices.data(),
1378  &values(i, 0),
1379  elide_zero_values);
1380  }
1381 
1382 
1383 
1384  void
1386  const std::vector<size_type> & col_indices,
1387  const std::vector<TrilinosScalar> &values,
1388  const bool elide_zero_values)
1389  {
1390  Assert(col_indices.size() == values.size(),
1391  ExcDimensionMismatch(col_indices.size(), values.size()));
1392 
1393  set(row,
1394  col_indices.size(),
1395  col_indices.data(),
1396  values.data(),
1397  elide_zero_values);
1398  }
1399 
1400 
1401 
1402  template <>
1403  void
1404  SparseMatrix::set<TrilinosScalar>(const size_type row,
1405  const size_type n_cols,
1406  const size_type * col_indices,
1407  const TrilinosScalar *values,
1408  const bool elide_zero_values)
1409  {
1410  AssertIndexRange(row, this->m());
1411 
1412  int ierr;
1413  if (last_action == Add)
1414  {
1415  ierr =
1416  matrix->GlobalAssemble(*column_space_map, matrix->RowMap(), true);
1417 
1418  Assert(ierr == 0, ExcTrilinosError(ierr));
1419  }
1420 
1421  last_action = Insert;
1422 
1423  const TrilinosWrappers::types::int_type *col_index_ptr;
1424  const TrilinosScalar * col_value_ptr;
1425  const TrilinosWrappers::types::int_type trilinos_row = row;
1427 
1428  boost::container::small_vector<TrilinosScalar, 200> local_value_array(
1429  n_cols);
1430  boost::container::small_vector<TrilinosWrappers::types::int_type, 200>
1431  local_index_array(n_cols);
1432 
1433  // If we don't elide zeros, the pointers are already available... need to
1434  // cast to non-const pointers as that is the format taken by Trilinos (but
1435  // we will not modify const data)
1436  if (elide_zero_values == false)
1437  {
1438  col_index_ptr =
1439  reinterpret_cast<const TrilinosWrappers::types::int_type *>(
1440  col_indices);
1441  col_value_ptr = values;
1442  n_columns = n_cols;
1443  }
1444  else
1445  {
1446  // Otherwise, extract nonzero values in each row and get the
1447  // respective indices.
1448  col_index_ptr = local_index_array.data();
1449  col_value_ptr = local_value_array.data();
1450 
1451  n_columns = 0;
1452  for (size_type j = 0; j < n_cols; ++j)
1453  {
1454  const double value = values[j];
1455  AssertIsFinite(value);
1456  if (value != 0)
1457  {
1458  local_index_array[n_columns] = col_indices[j];
1459  local_value_array[n_columns] = value;
1460  n_columns++;
1461  }
1462  }
1463 
1464  AssertIndexRange(n_columns, n_cols + 1);
1465  }
1466 
1467 
1468  // If the calling matrix owns the row to which we want to insert values,
1469  // we can directly call the Epetra_CrsMatrix input function, which is much
1470  // faster than the Epetra_FECrsMatrix function. We distinguish between two
1471  // cases: the first one is when the matrix is not filled (i.e., it is
1472  // possible to add new elements to the sparsity pattern), and the second
1473  // one is when the pattern is already fixed. In the former case, we add
1474  // the possibility to insert new values, and in the second we just replace
1475  // data.
1476  if (matrix->RowMap().MyGID(
1477  static_cast<TrilinosWrappers::types::int_type>(row)) == true)
1478  {
1479  if (matrix->Filled() == false)
1480  {
1481  ierr = matrix->Epetra_CrsMatrix::InsertGlobalValues(
1482  row, static_cast<int>(n_columns), col_value_ptr, col_index_ptr);
1483 
1484  // When inserting elements, we do not want to create exceptions in
1485  // the case when inserting non-local data (since that's what we
1486  // want to do right now).
1487  if (ierr > 0)
1488  ierr = 0;
1489  }
1490  else
1491  ierr = matrix->Epetra_CrsMatrix::ReplaceGlobalValues(row,
1492  n_columns,
1493  col_value_ptr,
1494  col_index_ptr);
1495  }
1496  else
1497  {
1498  // When we're at off-processor data, we have to stick with the
1499  // standard Insert/ReplaceGlobalValues function. Nevertheless, the way
1500  // we call it is the fastest one (any other will lead to repeated
1501  // allocation and deallocation of memory in order to call the function
1502  // we already use, which is very inefficient if writing one element at
1503  // a time).
1504  compressed = false;
1505 
1506  if (matrix->Filled() == false)
1507  {
1508  ierr = matrix->InsertGlobalValues(1,
1509  &trilinos_row,
1510  n_columns,
1511  col_index_ptr,
1512  &col_value_ptr,
1513  Epetra_FECrsMatrix::ROW_MAJOR);
1514  if (ierr > 0)
1515  ierr = 0;
1516  }
1517  else
1518  ierr = matrix->ReplaceGlobalValues(1,
1519  &trilinos_row,
1520  n_columns,
1521  col_index_ptr,
1522  &col_value_ptr,
1523  Epetra_FECrsMatrix::ROW_MAJOR);
1524  // use the FECrsMatrix facilities for set even in the case when we
1525  // have explicitly set the off-processor rows because that only works
1526  // properly when adding elements, not when setting them (since we want
1527  // to only touch elements that have been set explicitly, and there is
1528  // no way on the receiving processor to identify them otherwise)
1529  }
1530 
1531  Assert(ierr <= 0, ExcAccessToNonPresentElement(row, col_index_ptr[0]));
1532  AssertThrow(ierr >= 0, ExcTrilinosError(ierr));
1533  }
1534 
1535 
1536 
1537  void
1538  SparseMatrix::add(const std::vector<size_type> & indices,
1540  const bool elide_zero_values)
1541  {
1542  Assert(indices.size() == values.m(),
1543  ExcDimensionMismatch(indices.size(), values.m()));
1544  Assert(values.m() == values.n(), ExcNotQuadratic());
1545 
1546  for (size_type i = 0; i < indices.size(); ++i)
1547  add(indices[i],
1548  indices.size(),
1549  indices.data(),
1550  &values(i, 0),
1551  elide_zero_values);
1552  }
1553 
1554 
1555 
1556  void
1557  SparseMatrix::add(const std::vector<size_type> & row_indices,
1558  const std::vector<size_type> & col_indices,
1560  const bool elide_zero_values)
1561  {
1562  Assert(row_indices.size() == values.m(),
1563  ExcDimensionMismatch(row_indices.size(), values.m()));
1564  Assert(col_indices.size() == values.n(),
1565  ExcDimensionMismatch(col_indices.size(), values.n()));
1566 
1567  for (size_type i = 0; i < row_indices.size(); ++i)
1568  add(row_indices[i],
1569  col_indices.size(),
1570  col_indices.data(),
1571  &values(i, 0),
1572  elide_zero_values);
1573  }
1574 
1575 
1576 
1577  void
1579  const std::vector<size_type> & col_indices,
1580  const std::vector<TrilinosScalar> &values,
1581  const bool elide_zero_values)
1582  {
1583  Assert(col_indices.size() == values.size(),
1584  ExcDimensionMismatch(col_indices.size(), values.size()));
1585 
1586  add(row,
1587  col_indices.size(),
1588  col_indices.data(),
1589  values.data(),
1590  elide_zero_values);
1591  }
1592 
1593 
1594 
1595  void
1597  const size_type n_cols,
1598  const size_type * col_indices,
1599  const TrilinosScalar *values,
1600  const bool elide_zero_values,
1601  const bool /*col_indices_are_sorted*/)
1602  {
1603  AssertIndexRange(row, this->m());
1604  int ierr;
1605  if (last_action == Insert)
1606  {
1607  // TODO: this could lead to a dead lock when only one processor
1608  // calls GlobalAssemble.
1609  ierr =
1610  matrix->GlobalAssemble(*column_space_map, matrix->RowMap(), false);
1611 
1612  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1613  }
1614 
1615  last_action = Add;
1616 
1617  const TrilinosWrappers::types::int_type *col_index_ptr;
1618  const TrilinosScalar * col_value_ptr;
1619  const TrilinosWrappers::types::int_type trilinos_row = row;
1621 
1622  boost::container::small_vector<TrilinosScalar, 100> local_value_array(
1623  n_cols);
1624  boost::container::small_vector<TrilinosWrappers::types::int_type, 100>
1625  local_index_array(n_cols);
1626 
1627  // If we don't elide zeros, the pointers are already available... need to
1628  // cast to non-const pointers as that is the format taken by Trilinos (but
1629  // we will not modify const data)
1630  if (elide_zero_values == false)
1631  {
1632  col_index_ptr =
1633  reinterpret_cast<const TrilinosWrappers::types::int_type *>(
1634  col_indices);
1635  col_value_ptr = values;
1636  n_columns = n_cols;
1637 # ifdef DEBUG
1638  for (size_type j = 0; j < n_cols; ++j)
1639  AssertIsFinite(values[j]);
1640 # endif
1641  }
1642  else
1643  {
1644  // Otherwise, extract nonzero values in each row and the corresponding
1645  // index.
1646  col_index_ptr = local_index_array.data();
1647  col_value_ptr = local_value_array.data();
1648 
1649  n_columns = 0;
1650  for (size_type j = 0; j < n_cols; ++j)
1651  {
1652  const double value = values[j];
1653 
1654  AssertIsFinite(value);
1655  if (value != 0)
1656  {
1657  local_index_array[n_columns] = col_indices[j];
1658  local_value_array[n_columns] = value;
1659  ++n_columns;
1660  }
1661  }
1662 
1663  AssertIndexRange(n_columns, n_cols + 1);
1664  }
1665  // Exit early if there is nothing to do
1666  if (n_columns == 0)
1667  {
1668  return;
1669  }
1670 
1671  // If the calling processor owns the row to which we want to add values, we
1672  // can directly call the Epetra_CrsMatrix input function, which is much
1673  // faster than the Epetra_FECrsMatrix function.
1674  if (matrix->RowMap().MyGID(
1675  static_cast<TrilinosWrappers::types::int_type>(row)) == true)
1676  {
1677  ierr = matrix->Epetra_CrsMatrix::SumIntoGlobalValues(row,
1678  n_columns,
1679  col_value_ptr,
1680  col_index_ptr);
1681  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1682  }
1683  else if (nonlocal_matrix.get() != nullptr)
1684  {
1685  compressed = false;
1686  // this is the case when we have explicitly set the off-processor rows
1687  // and want to create a separate matrix object for them (to retain
1688  // thread-safety)
1689  Assert(nonlocal_matrix->RowMap().LID(
1690  static_cast<TrilinosWrappers::types::int_type>(row)) != -1,
1691  ExcMessage("Attempted to write into off-processor matrix row "
1692  "that has not be specified as being writable upon "
1693  "initialization"));
1694  ierr = nonlocal_matrix->SumIntoGlobalValues(row,
1695  n_columns,
1696  col_value_ptr,
1697  col_index_ptr);
1698  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1699  }
1700  else
1701  {
1702  // When we're at off-processor data, we have to stick with the
1703  // standard SumIntoGlobalValues function. Nevertheless, the way we
1704  // call it is the fastest one (any other will lead to repeated
1705  // allocation and deallocation of memory in order to call the function
1706  // we already use, which is very inefficient if writing one element at
1707  // a time).
1708  compressed = false;
1709 
1710  ierr = matrix->SumIntoGlobalValues(1,
1711  &trilinos_row,
1712  n_columns,
1713  col_index_ptr,
1714  &col_value_ptr,
1715  Epetra_FECrsMatrix::ROW_MAJOR);
1716  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1717  }
1718 
1719 # ifdef DEBUG
1720  if (ierr > 0)
1721  {
1722  std::cout << "------------------------------------------" << std::endl;
1723  std::cout << "Got error " << ierr << " in row " << row << " of proc "
1724  << matrix->RowMap().Comm().MyPID()
1725  << " when trying to add the columns:" << std::endl;
1726  for (TrilinosWrappers::types::int_type i = 0; i < n_columns; ++i)
1727  std::cout << col_index_ptr[i] << " ";
1728  std::cout << std::endl << std::endl;
1729  std::cout << "Matrix row "
1730  << (matrix->RowMap().MyGID(
1731  static_cast<TrilinosWrappers::types::int_type>(row)) ==
1732  false ?
1733  "(nonlocal part)" :
1734  "")
1735  << " has the following indices:" << std::endl;
1736  std::vector<TrilinosWrappers::types::int_type> indices;
1737  const Epetra_CrsGraph * graph =
1738  (nonlocal_matrix.get() != nullptr &&
1739  matrix->RowMap().MyGID(
1740  static_cast<TrilinosWrappers::types::int_type>(row)) == false) ?
1741  &nonlocal_matrix->Graph() :
1742  &matrix->Graph();
1743 
1744  indices.resize(graph->NumGlobalIndices(row));
1745  int n_indices = 0;
1746  graph->ExtractGlobalRowCopy(row,
1747  indices.size(),
1748  n_indices,
1749  indices.data());
1750  AssertDimension(n_indices, indices.size());
1751 
1752  for (TrilinosWrappers::types::int_type i = 0; i < n_indices; ++i)
1753  std::cout << indices[i] << " ";
1754  std::cout << std::endl << std::endl;
1755  Assert(ierr <= 0, ExcAccessToNonPresentElement(row, col_index_ptr[0]));
1756  }
1757 # endif
1758  Assert(ierr >= 0, ExcTrilinosError(ierr));
1759  }
1760 
1761 
1762 
1763  SparseMatrix &
1765  {
1767  compress(
1768  ::VectorOperation::unknown); // TODO: why do we do this? Should we
1769  // not check for is_compressed?
1770 
1771  const int ierr = matrix->PutScalar(d);
1772  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1773  if (nonlocal_matrix.get() != nullptr)
1774  nonlocal_matrix->PutScalar(d);
1775 
1776  return *this;
1777  }
1778 
1779 
1780 
1781  void
1783  {
1784  AssertDimension(rhs.m(), m());
1785  AssertDimension(rhs.n(), n());
1786  AssertDimension(rhs.local_range().first, local_range().first);
1787  AssertDimension(rhs.local_range().second, local_range().second);
1788  Assert(matrix->RowMap().SameAs(rhs.matrix->RowMap()),
1789  ExcMessage("Can only add matrices with same distribution of rows"));
1790  Assert(matrix->Filled() && rhs.matrix->Filled(),
1791  ExcMessage("Addition of matrices only allowed if matrices are "
1792  "filled, i.e., compress() has been called"));
1793 
1794  const bool same_col_map = matrix->ColMap().SameAs(rhs.matrix->ColMap());
1795 
1796  for (const auto row : locally_owned_range_indices())
1797  {
1798  const int row_local = matrix->RowMap().LID(
1799  static_cast<TrilinosWrappers::types::int_type>(row));
1800  Assert((row_local >= 0), ExcAccessToNonlocalRow(row));
1801 
1802  // First get a view to the matrix columns of both matrices. Note that
1803  // the data is in local index spaces so we need to be careful not only
1804  // to compare column indices in case they are derived from the same
1805  // map.
1806  int n_entries, rhs_n_entries;
1807  TrilinosScalar *value_ptr, *rhs_value_ptr;
1808  int * index_ptr, *rhs_index_ptr;
1809  int ierr = rhs.matrix->ExtractMyRowView(row_local,
1810  rhs_n_entries,
1811  rhs_value_ptr,
1812  rhs_index_ptr);
1813  (void)ierr;
1814  Assert(ierr == 0, ExcTrilinosError(ierr));
1815 
1816  ierr =
1817  matrix->ExtractMyRowView(row_local, n_entries, value_ptr, index_ptr);
1818  Assert(ierr == 0, ExcTrilinosError(ierr));
1819  bool expensive_checks = (n_entries != rhs_n_entries || !same_col_map);
1820  if (!expensive_checks)
1821  {
1822  // check if the column indices are the same. If yes, can simply
1823  // copy over the data.
1824  expensive_checks = std::memcmp(static_cast<void *>(index_ptr),
1825  static_cast<void *>(rhs_index_ptr),
1826  sizeof(int) * n_entries) != 0;
1827  if (!expensive_checks)
1828  for (int i = 0; i < n_entries; ++i)
1829  value_ptr[i] += rhs_value_ptr[i] * factor;
1830  }
1831  // Now to the expensive case where we need to check all column indices
1832  // against each other (transformed into global index space) and where
1833  // we need to make sure that all entries we are about to add into the
1834  // lhs matrix actually exist
1835  if (expensive_checks)
1836  {
1837  for (int i = 0; i < rhs_n_entries; ++i)
1838  {
1839  if (rhs_value_ptr[i] == 0.)
1840  continue;
1841  const TrilinosWrappers::types::int_type rhs_global_col =
1842  global_column_index(*rhs.matrix, rhs_index_ptr[i]);
1843  int local_col = matrix->ColMap().LID(rhs_global_col);
1844  int *local_index = Utilities::lower_bound(index_ptr,
1845  index_ptr + n_entries,
1846  local_col);
1847  Assert(local_index != index_ptr + n_entries &&
1848  *local_index == local_col,
1849  ExcMessage(
1850  "Adding the entries from the other matrix "
1851  "failed, because the sparsity pattern "
1852  "of that matrix includes more elements than the "
1853  "calling matrix, which is not allowed."));
1854  value_ptr[local_index - index_ptr] += factor * rhs_value_ptr[i];
1855  }
1856  }
1857  }
1858  }
1859 
1860 
1861 
1862  void
1864  {
1865  // This only flips a flag that tells
1866  // Trilinos that any vmult operation
1867  // should be done with the
1868  // transpose. However, the matrix
1869  // structure is not reset.
1870  int ierr;
1871 
1872  if (!matrix->UseTranspose())
1873  {
1874  ierr = matrix->SetUseTranspose(true);
1875  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1876  }
1877  else
1878  {
1879  ierr = matrix->SetUseTranspose(false);
1880  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1881  }
1882  }
1883 
1884 
1885 
1886  SparseMatrix &
1888  {
1889  const int ierr = matrix->Scale(a);
1890  Assert(ierr == 0, ExcTrilinosError(ierr));
1891  (void)ierr; // removes -Wunused-variable in optimized mode
1892 
1893  return *this;
1894  }
1895 
1896 
1897 
1898  SparseMatrix &
1900  {
1901  Assert(a != 0, ExcDivideByZero());
1902 
1903  const TrilinosScalar factor = 1. / a;
1904 
1905  const int ierr = matrix->Scale(factor);
1906  Assert(ierr == 0, ExcTrilinosError(ierr));
1907  (void)ierr; // removes -Wunused-variable in optimized mode
1908 
1909  return *this;
1910  }
1911 
1912 
1913 
1916  {
1917  Assert(matrix->Filled(), ExcMatrixNotCompressed());
1918  return matrix->NormOne();
1919  }
1920 
1921 
1922 
1925  {
1926  Assert(matrix->Filled(), ExcMatrixNotCompressed());
1927  return matrix->NormInf();
1928  }
1929 
1930 
1931 
1934  {
1935  Assert(matrix->Filled(), ExcMatrixNotCompressed());
1936  return matrix->NormFrobenius();
1937  }
1938 
1939 
1940 
1941  namespace internal
1942  {
1943  namespace SparseMatrixImplementation
1944  {
1945  template <typename VectorType>
1946  inline void
1947  check_vector_map_equality(const Epetra_CrsMatrix &,
1948  const VectorType &,
1949  const VectorType &)
1950  {}
1951 
1952  inline void
1953  check_vector_map_equality(const Epetra_CrsMatrix & m,
1955  const TrilinosWrappers::MPI::Vector &out)
1956  {
1957  Assert(in.trilinos_partitioner().SameAs(m.DomainMap()) == true,
1958  ExcMessage("The column partitioning of a matrix does not match "
1959  "the partitioning of a vector you are trying to "
1960  "multiply it with. Are you multiplying the "
1961  "matrix with a vector that has ghost elements?"));
1962  Assert(out.trilinos_partitioner().SameAs(m.RangeMap()) == true,
1963  ExcMessage("The row partitioning of a matrix does not match "
1964  "the partitioning of a vector you are trying to "
1965  "put the result of a matrix-vector product in. "
1966  "Are you trying to put the product of the "
1967  "matrix with a vector into a vector that has "
1968  "ghost elements?"));
1969  (void)m;
1970  (void)in;
1971  (void)out;
1972  }
1973  } // namespace SparseMatrixImplementation
1974  } // namespace internal
1975 
1976 
1977  template <typename VectorType>
1978  std::enable_if_t<
1979  std::is_same<typename VectorType::value_type, TrilinosScalar>::value>
1980  SparseMatrix::vmult(VectorType &dst, const VectorType &src) const
1981  {
1982  Assert(&src != &dst, ExcSourceEqualsDestination());
1983  Assert(matrix->Filled(), ExcMatrixNotCompressed());
1984  (void)src;
1985  (void)dst;
1986 
1988  src,
1989  dst);
1990  const size_type dst_local_size = internal::end(dst) - internal::begin(dst);
1991  AssertDimension(dst_local_size, matrix->RangeMap().NumMyPoints());
1992  const size_type src_local_size = internal::end(src) - internal::begin(src);
1993  AssertDimension(src_local_size, matrix->DomainMap().NumMyPoints());
1994 
1995  Epetra_MultiVector tril_dst(
1996  View, matrix->RangeMap(), internal::begin(dst), dst_local_size, 1);
1997  Epetra_MultiVector tril_src(View,
1998  matrix->DomainMap(),
1999  const_cast<TrilinosScalar *>(
2000  internal::begin(src)),
2001  src_local_size,
2002  1);
2003 
2004  const int ierr = matrix->Multiply(false, tril_src, tril_dst);
2005  Assert(ierr == 0, ExcTrilinosError(ierr));
2006  (void)ierr; // removes -Wunused-variable in optimized mode
2007  }
2008 
2009 
2010 
2011  template <typename VectorType>
2012  std::enable_if_t<
2013  !std::is_same<typename VectorType::value_type, TrilinosScalar>::value>
2014  SparseMatrix::vmult(VectorType & /*dst*/, const VectorType & /*src*/) const
2015  {
2016  AssertThrow(false, ExcNotImplemented());
2017  }
2018 
2019 
2020 
2021  template <typename VectorType>
2022  std::enable_if_t<
2023  std::is_same<typename VectorType::value_type, TrilinosScalar>::value>
2024  SparseMatrix::Tvmult(VectorType &dst, const VectorType &src) const
2025  {
2026  Assert(&src != &dst, ExcSourceEqualsDestination());
2027  Assert(matrix->Filled(), ExcMatrixNotCompressed());
2028 
2030  dst,
2031  src);
2032  const size_type dst_local_size = internal::end(dst) - internal::begin(dst);
2033  AssertDimension(dst_local_size, matrix->DomainMap().NumMyPoints());
2034  const size_type src_local_size = internal::end(src) - internal::begin(src);
2035  AssertDimension(src_local_size, matrix->RangeMap().NumMyPoints());
2036 
2037  Epetra_MultiVector tril_dst(
2038  View, matrix->DomainMap(), internal::begin(dst), dst_local_size, 1);
2039  Epetra_MultiVector tril_src(View,
2040  matrix->RangeMap(),
2041  const_cast<double *>(internal::begin(src)),
2042  src_local_size,
2043  1);
2044 
2045  const int ierr = matrix->Multiply(true, tril_src, tril_dst);
2046  Assert(ierr == 0, ExcTrilinosError(ierr));
2047  (void)ierr; // removes -Wunused-variable in optimized mode
2048  }
2049 
2050 
2051 
2052  template <typename VectorType>
2053  std::enable_if_t<
2054  !std::is_same<typename VectorType::value_type, TrilinosScalar>::value>
2055  SparseMatrix::Tvmult(VectorType & /*dst*/, const VectorType & /*src*/) const
2056  {
2057  AssertThrow(false, ExcNotImplemented());
2058  }
2059 
2060 
2061 
2062  template <typename VectorType>
2063  void
2064  SparseMatrix::vmult_add(VectorType &dst, const VectorType &src) const
2065  {
2066  Assert(&src != &dst, ExcSourceEqualsDestination());
2067 
2068  // Reinit a temporary vector with fast argument set, which does not
2069  // overwrite the content (to save time).
2070  VectorType tmp_vector;
2071  tmp_vector.reinit(dst, true);
2072  vmult(tmp_vector, src);
2073  dst += tmp_vector;
2074  }
2075 
2076 
2077 
2078  template <typename VectorType>
2079  void
2080  SparseMatrix::Tvmult_add(VectorType &dst, const VectorType &src) const
2081  {
2082  Assert(&src != &dst, ExcSourceEqualsDestination());
2083 
2084  // Reinit a temporary vector with fast argument set, which does not
2085  // overwrite the content (to save time).
2086  VectorType tmp_vector;
2087  tmp_vector.reinit(dst, true);
2088  Tvmult(tmp_vector, src);
2089  dst += tmp_vector;
2090  }
2091 
2092 
2093 
2096  {
2097  Assert(matrix->RowMap().SameAs(matrix->DomainMap()), ExcNotQuadratic());
2098 
2099  MPI::Vector temp_vector;
2100  temp_vector.reinit(v, true);
2101 
2102  vmult(temp_vector, v);
2103  return temp_vector * v;
2104  }
2105 
2106 
2107 
2110  const MPI::Vector &v) const
2111  {
2112  Assert(matrix->RowMap().SameAs(matrix->DomainMap()), ExcNotQuadratic());
2113 
2114  MPI::Vector temp_vector;
2115  temp_vector.reinit(v, true);
2116 
2117  vmult(temp_vector, v);
2118  return u * temp_vector;
2119  }
2120 
2121 
2122 
2125  const MPI::Vector &x,
2126  const MPI::Vector &b) const
2127  {
2128  vmult(dst, x);
2129  dst -= b;
2130  dst *= -1.;
2131 
2132  return dst.l2_norm();
2133  }
2134 
2135 
2136 
2137  namespace internals
2138  {
2140 
2141  void
2142  perform_mmult(const SparseMatrix &inputleft,
2143  const SparseMatrix &inputright,
2144  SparseMatrix & result,
2145  const MPI::Vector & V,
2146  const bool transpose_left)
2147  {
2148  const bool use_vector = (V.size() == inputright.m() ? true : false);
2149  if (transpose_left == false)
2150  {
2151  Assert(inputleft.n() == inputright.m(),
2152  ExcDimensionMismatch(inputleft.n(), inputright.m()));
2153  Assert(inputleft.trilinos_matrix().DomainMap().SameAs(
2154  inputright.trilinos_matrix().RangeMap()),
2155  ExcMessage("Parallel partitioning of A and B does not fit."));
2156  }
2157  else
2158  {
2159  Assert(inputleft.m() == inputright.m(),
2160  ExcDimensionMismatch(inputleft.m(), inputright.m()));
2161  Assert(inputleft.trilinos_matrix().RangeMap().SameAs(
2162  inputright.trilinos_matrix().RangeMap()),
2163  ExcMessage("Parallel partitioning of A and B does not fit."));
2164  }
2165 
2166  result.clear();
2167 
2168  // create a suitable operator B: in case
2169  // we do not use a vector, all we need to
2170  // do is to set the pointer. Otherwise,
2171  // we insert the data from B, but
2172  // multiply each row with the respective
2173  // vector element.
2174  Teuchos::RCP<Epetra_CrsMatrix> mod_B;
2175  if (use_vector == false)
2176  {
2177  mod_B = Teuchos::rcp(const_cast<Epetra_CrsMatrix *>(
2178  &inputright.trilinos_matrix()),
2179  false);
2180  }
2181  else
2182  {
2183  mod_B = Teuchos::rcp(
2184  new Epetra_CrsMatrix(Copy, inputright.trilinos_sparsity_pattern()),
2185  true);
2186  mod_B->FillComplete(inputright.trilinos_matrix().DomainMap(),
2187  inputright.trilinos_matrix().RangeMap());
2188  Assert(inputright.local_range() == V.local_range(),
2189  ExcMessage("Parallel distribution of matrix B and vector V "
2190  "does not match."));
2191 
2192  const int local_N = inputright.local_size();
2193  for (int i = 0; i < local_N; ++i)
2194  {
2195  int N_entries = -1;
2196  double *new_data, *B_data;
2197  mod_B->ExtractMyRowView(i, N_entries, new_data);
2198  inputright.trilinos_matrix().ExtractMyRowView(i,
2199  N_entries,
2200  B_data);
2201  double value = V.trilinos_vector()[0][i];
2202  for (TrilinosWrappers::types::int_type j = 0; j < N_entries; ++j)
2203  new_data[j] = value * B_data[j];
2204  }
2205  }
2206 
2207 
2208  SparseMatrix tmp_result(transpose_left ?
2209  inputleft.locally_owned_domain_indices() :
2210  inputleft.locally_owned_range_indices(),
2211  inputright.locally_owned_domain_indices(),
2212  inputleft.get_mpi_communicator());
2213 
2214 # ifdef DEAL_II_TRILINOS_WITH_EPETRAEXT
2215  EpetraExt::MatrixMatrix::Multiply(inputleft.trilinos_matrix(),
2216  transpose_left,
2217  *mod_B,
2218  false,
2219  const_cast<Epetra_CrsMatrix &>(
2220  tmp_result.trilinos_matrix()));
2221 # else
2222  Assert(false,
2223  ExcMessage("This function requires that the Trilinos "
2224  "installation found while running the deal.II "
2225  "CMake scripts contains the optional Trilinos "
2226  "package 'EpetraExt'. However, this optional "
2227  "part of Trilinos was not found."));
2228 # endif
2229  result.reinit(tmp_result.trilinos_matrix());
2230  }
2231  } // namespace internals
2232 
2233 
2234  void
2236  const SparseMatrix &B,
2237  const MPI::Vector & V) const
2238  {
2239  internals::perform_mmult(*this, B, C, V, false);
2240  }
2241 
2242 
2243 
2244  void
2246  const SparseMatrix &B,
2247  const MPI::Vector & V) const
2248  {
2249  internals::perform_mmult(*this, B, C, V, true);
2250  }
2251 
2252 
2253 
2254  void
2256  {
2257  Assert(false, ExcNotImplemented());
2258  }
2259 
2260 
2261 
2262  // As of now, no particularly neat
2263  // output is generated in case of
2264  // multiple processors.
2265  void
2266  SparseMatrix::print(std::ostream &out,
2267  const bool print_detailed_trilinos_information) const
2268  {
2269  if (print_detailed_trilinos_information == true)
2270  out << *matrix;
2271  else
2272  {
2273  double *values;
2274  int * indices;
2275  int num_entries;
2276 
2277  for (int i = 0; i < matrix->NumMyRows(); ++i)
2278  {
2279  const int ierr =
2280  matrix->ExtractMyRowView(i, num_entries, values, indices);
2281  (void)ierr;
2282  Assert(ierr == 0, ExcTrilinosError(ierr));
2283 
2284  for (TrilinosWrappers::types::int_type j = 0; j < num_entries; ++j)
2285  out << "(" << TrilinosWrappers::global_row_index(*matrix, i)
2286  << ","
2288  << ") " << values[j] << std::endl;
2289  }
2290  }
2291 
2292  AssertThrow(out.fail() == false, ExcIO());
2293  }
2294 
2295 
2296 
2299  {
2300  size_type static_memory =
2301  sizeof(*this) + sizeof(*matrix) + sizeof(*matrix->Graph().DataPtr());
2302  return (
2304  matrix->NumMyNonzeros() +
2305  sizeof(int) * local_size() + static_memory);
2306  }
2307 
2308 
2309 
2310  MPI_Comm
2312  {
2313  const Epetra_MpiComm *mpi_comm =
2314  dynamic_cast<const Epetra_MpiComm *>(&matrix->RangeMap().Comm());
2315  Assert(mpi_comm != nullptr, ExcInternalError());
2316  return mpi_comm->Comm();
2317  }
2318 } // namespace TrilinosWrappers
2319 
2320 
2321 namespace TrilinosWrappers
2322 {
2323  namespace internal
2324  {
2325  namespace LinearOperatorImplementation
2326  {
2328  : use_transpose(false)
2329  , communicator(MPI_COMM_SELF)
2330  , domain_map(IndexSet().make_trilinos_map(communicator.Comm()))
2331  , range_map(IndexSet().make_trilinos_map(communicator.Comm()))
2332  {
2333  vmult = [](Range &, const Domain &) {
2334  Assert(false,
2335  ExcMessage("Uninitialized TrilinosPayload::vmult called "
2336  "(Default constructor)"));
2337  };
2338 
2339  Tvmult = [](Domain &, const Range &) {
2340  Assert(false,
2341  ExcMessage("Uninitialized TrilinosPayload::Tvmult called "
2342  "(Default constructor)"));
2343  };
2344 
2345  inv_vmult = [](Domain &, const Range &) {
2346  Assert(false,
2347  ExcMessage("Uninitialized TrilinosPayload::inv_vmult called "
2348  "(Default constructor)"));
2349  };
2350 
2351  inv_Tvmult = [](Range &, const Domain &) {
2352  Assert(false,
2353  ExcMessage("Uninitialized TrilinosPayload::inv_Tvmult called "
2354  "(Default constructor)"));
2355  };
2356  }
2357 
2358 
2359 
2361  const TrilinosWrappers::SparseMatrix &matrix_exemplar,
2363  : TrilinosPayload(const_cast<Epetra_CrsMatrix &>(
2364  matrix.trilinos_matrix()),
2365  /*op_supports_inverse_operations = */ false,
2366  matrix_exemplar.trilinos_matrix().UseTranspose(),
2367  matrix_exemplar.get_mpi_communicator(),
2368  matrix_exemplar.locally_owned_domain_indices(),
2369  matrix_exemplar.locally_owned_range_indices())
2370  {}
2371 
2372 
2373 
2375  const TrilinosPayload & payload_exemplar,
2377 
2378  : TrilinosPayload(const_cast<Epetra_CrsMatrix &>(
2379  matrix.trilinos_matrix()),
2380  /*op_supports_inverse_operations = */ false,
2381  payload_exemplar.UseTranspose(),
2382  payload_exemplar.get_mpi_communicator(),
2383  payload_exemplar.locally_owned_domain_indices(),
2384  payload_exemplar.locally_owned_range_indices())
2385  {}
2386 
2387 
2388 
2390  const TrilinosWrappers::SparseMatrix & matrix_exemplar,
2391  const TrilinosWrappers::PreconditionBase &preconditioner)
2392  : TrilinosPayload(preconditioner.trilinos_operator(),
2393  /*op_supports_inverse_operations = */ true,
2394  matrix_exemplar.trilinos_matrix().UseTranspose(),
2395  matrix_exemplar.get_mpi_communicator(),
2396  matrix_exemplar.locally_owned_domain_indices(),
2397  matrix_exemplar.locally_owned_range_indices())
2398  {}
2399 
2400 
2401 
2403  const TrilinosWrappers::PreconditionBase &preconditioner_exemplar,
2404  const TrilinosWrappers::PreconditionBase &preconditioner)
2405  : TrilinosPayload(
2406  preconditioner.trilinos_operator(),
2407  /*op_supports_inverse_operations = */ true,
2408  preconditioner_exemplar.trilinos_operator().UseTranspose(),
2409  preconditioner_exemplar.get_mpi_communicator(),
2410  preconditioner_exemplar.locally_owned_domain_indices(),
2411  preconditioner_exemplar.locally_owned_range_indices())
2412  {}
2413 
2414 
2415 
2417  const TrilinosPayload & payload_exemplar,
2418  const TrilinosWrappers::PreconditionBase &preconditioner)
2419  : TrilinosPayload(preconditioner.trilinos_operator(),
2420  /*op_supports_inverse_operations = */ true,
2421  payload_exemplar.UseTranspose(),
2422  payload_exemplar.get_mpi_communicator(),
2423  payload_exemplar.locally_owned_domain_indices(),
2424  payload_exemplar.locally_owned_range_indices())
2425  {}
2426 
2427 
2428 
2430  : vmult(payload.vmult)
2431  , Tvmult(payload.Tvmult)
2432  , inv_vmult(payload.inv_vmult)
2433  , inv_Tvmult(payload.inv_Tvmult)
2434  , use_transpose(payload.use_transpose)
2435  , communicator(payload.communicator)
2436  , domain_map(payload.domain_map)
2437  , range_map(payload.range_map)
2438  {}
2439 
2440 
2441 
2442  // Composite copy constructor
2443  // This is required for PackagedOperations
2445  const TrilinosPayload &second_op)
2446  : use_transpose(false)
2447  , // The combination of operators provides the exact
2448  // definition of the operation
2449  communicator(first_op.communicator)
2450  , domain_map(second_op.domain_map)
2451  , range_map(first_op.range_map)
2452  {}
2453 
2454 
2455 
2458  {
2459  TrilinosPayload return_op(*this);
2460 
2461  return_op.vmult = [](Range &tril_dst, const Range &tril_src) {
2462  tril_dst = tril_src;
2463  };
2464 
2465  return_op.Tvmult = [](Range &tril_dst, const Range &tril_src) {
2466  tril_dst = tril_src;
2467  };
2468 
2469  return_op.inv_vmult = [](Range &tril_dst, const Range &tril_src) {
2470  tril_dst = tril_src;
2471  };
2472 
2473  return_op.inv_Tvmult = [](Range &tril_dst, const Range &tril_src) {
2474  tril_dst = tril_src;
2475  };
2476 
2477  return return_op;
2478  }
2479 
2480 
2481 
2484  {
2485  TrilinosPayload return_op(*this);
2486 
2487  return_op.vmult = [](Range &tril_dst, const Domain &) {
2488  const int ierr = tril_dst.PutScalar(0.0);
2489 
2490  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2491  };
2492 
2493  return_op.Tvmult = [](Domain &tril_dst, const Range &) {
2494  const int ierr = tril_dst.PutScalar(0.0);
2495 
2496  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2497  };
2498 
2499  return_op.inv_vmult = [](Domain &tril_dst, const Range &) {
2500  AssertThrow(false,
2501  ExcMessage("Cannot compute inverse of null operator"));
2502 
2503  const int ierr = tril_dst.PutScalar(0.0);
2504 
2505  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2506  };
2507 
2508  return_op.inv_Tvmult = [](Range &tril_dst, const Domain &) {
2509  AssertThrow(false,
2510  ExcMessage("Cannot compute inverse of null operator"));
2511 
2512  const int ierr = tril_dst.PutScalar(0.0);
2513 
2514  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2515  };
2516 
2517  return return_op;
2518  }
2519 
2520 
2521 
2524  {
2525  TrilinosPayload return_op(*this);
2526  return_op.transpose();
2527  return return_op;
2528  }
2529 
2530 
2531 
2532  IndexSet
2534  {
2535  return IndexSet(domain_map);
2536  }
2537 
2538 
2539 
2540  IndexSet
2542  {
2543  return IndexSet(range_map);
2544  }
2545 
2546 
2547 
2548  MPI_Comm
2550  {
2551  return communicator.Comm();
2552  }
2553 
2554 
2555 
2556  void
2558  {
2560  }
2561 
2562 
2563 
2564  bool
2566  {
2567  return use_transpose;
2568  }
2569 
2570 
2571 
2572  int
2574  {
2575  if (use_transpose != UseTranspose)
2576  {
2581  }
2582  return 0;
2583  }
2584 
2585 
2586 
2587  int
2589  {
2590  // The transposedness of the operations is taken care of
2591  // when we hit the transpose flag.
2592  vmult(Y, X);
2593  return 0;
2594  }
2595 
2596 
2597 
2598  int
2600  {
2601  // The transposedness of the operations is taken care of
2602  // when we hit the transpose flag.
2603  inv_vmult(X, Y);
2604  return 0;
2605  }
2606 
2607 
2608 
2609  const char *
2611  {
2612  return "TrilinosPayload";
2613  }
2614 
2615 
2616 
2617  const Epetra_Comm &
2619  {
2620  return communicator;
2621  }
2622 
2623 
2624 
2625  const Epetra_Map &
2627  {
2628  return domain_map;
2629  }
2630 
2631 
2632 
2633  const Epetra_Map &
2635  {
2636  return range_map;
2637  }
2638 
2639 
2640 
2641  bool
2643  {
2644  return false;
2645  }
2646 
2647 
2648 
2649  double
2651  {
2652  AssertThrow(false, ExcNotImplemented());
2653  return 0.0;
2654  }
2655 
2656 
2657 
2659  operator+(const TrilinosPayload &first_op,
2660  const TrilinosPayload &second_op)
2661  {
2662  using Domain = typename TrilinosPayload::Domain;
2663  using Range = typename TrilinosPayload::Range;
2664  using Intermediate = typename TrilinosPayload::VectorType;
2665  using GVMVectorType = TrilinosWrappers::MPI::Vector;
2666 
2667  Assert(first_op.locally_owned_domain_indices() ==
2668  second_op.locally_owned_domain_indices(),
2669  ExcMessage(
2670  "Operators are set to work on incompatible IndexSets."));
2671  Assert(first_op.locally_owned_range_indices() ==
2672  second_op.locally_owned_range_indices(),
2673  ExcMessage(
2674  "Operators are set to work on incompatible IndexSets."));
2675 
2676  TrilinosPayload return_op(first_op, second_op);
2677 
2678  // Capture by copy so the payloads are always valid
2679  return_op.vmult = [first_op, second_op](Range & tril_dst,
2680  const Domain &tril_src) {
2681  // Duplicated from LinearOperator::operator*
2682  // TODO: Template the constructor on GrowingVectorMemory vector type?
2683  GrowingVectorMemory<GVMVectorType> vector_memory;
2684  VectorMemory<GVMVectorType>::Pointer i(vector_memory);
2685 
2686  // Initialize intermediate vector
2687  const Epetra_Map &first_op_init_map = first_op.OperatorDomainMap();
2688  i->reinit(IndexSet(first_op_init_map),
2689  first_op.get_mpi_communicator(),
2690  /*bool omit_zeroing_entries =*/true);
2691 
2692  // Duplicated from TrilinosWrappers::SparseMatrix::vmult
2693  const size_type i_local_size = i->end() - i->begin();
2694  AssertDimension(i_local_size, first_op_init_map.NumMyPoints());
2695  const Epetra_Map &second_op_init_map = second_op.OperatorDomainMap();
2696  AssertDimension(i_local_size, second_op_init_map.NumMyPoints());
2697  (void)second_op_init_map;
2698  Intermediate tril_int(View,
2699  first_op_init_map,
2700  const_cast<TrilinosScalar *>(i->begin()),
2701  i_local_size,
2702  1);
2703 
2704  // These operators may themselves be transposed or not, so we let them
2705  // decide what the intended outcome is
2706  second_op.Apply(tril_src, tril_int);
2707  first_op.Apply(tril_src, tril_dst);
2708  const int ierr = tril_dst.Update(1.0, tril_int, 1.0);
2709  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2710  };
2711 
2712  return_op.Tvmult = [first_op, second_op](Domain & tril_dst,
2713  const Range &tril_src) {
2714  // Duplicated from LinearOperator::operator*
2715  // TODO: Template the constructor on GrowingVectorMemory vector type?
2716  GrowingVectorMemory<GVMVectorType> vector_memory;
2717  VectorMemory<GVMVectorType>::Pointer i(vector_memory);
2718 
2719  // These operators may themselves be transposed or not, so we let them
2720  // decide what the intended outcome is
2721  // We must first transpose the operators to get the right IndexSets
2722  // for the input, intermediate and result vectors
2723  const_cast<TrilinosPayload &>(first_op).transpose();
2724  const_cast<TrilinosPayload &>(second_op).transpose();
2725 
2726  // Initialize intermediate vector
2727  const Epetra_Map &first_op_init_map = first_op.OperatorRangeMap();
2728  i->reinit(IndexSet(first_op_init_map),
2729  first_op.get_mpi_communicator(),
2730  /*bool omit_zeroing_entries =*/true);
2731 
2732  // Duplicated from TrilinosWrappers::SparseMatrix::vmult
2733  const size_type i_local_size = i->end() - i->begin();
2734  AssertDimension(i_local_size, first_op_init_map.NumMyPoints());
2735  const Epetra_Map &second_op_init_map = second_op.OperatorRangeMap();
2736  AssertDimension(i_local_size, second_op_init_map.NumMyPoints());
2737  (void)second_op_init_map;
2738  Intermediate tril_int(View,
2739  first_op_init_map,
2740  const_cast<TrilinosScalar *>(i->begin()),
2741  i_local_size,
2742  1);
2743 
2744  // These operators may themselves be transposed or not, so we let them
2745  // decide what the intended outcome is
2746  second_op.Apply(tril_src, tril_int);
2747  first_op.Apply(tril_src, tril_dst);
2748  const int ierr = tril_dst.Update(1.0, tril_int, 1.0);
2749  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2750 
2751  // Reset transpose flag
2752  const_cast<TrilinosPayload &>(first_op).transpose();
2753  const_cast<TrilinosPayload &>(second_op).transpose();
2754  };
2755 
2756  return_op.inv_vmult = [first_op, second_op](Domain & tril_dst,
2757  const Range &tril_src) {
2758  // Duplicated from LinearOperator::operator*
2759  // TODO: Template the constructor on GrowingVectorMemory vector type?
2760  GrowingVectorMemory<GVMVectorType> vector_memory;
2761  VectorMemory<GVMVectorType>::Pointer i(vector_memory);
2762 
2763  // Initialize intermediate vector
2764  const Epetra_Map &first_op_init_map = first_op.OperatorRangeMap();
2765  i->reinit(IndexSet(first_op_init_map),
2766  first_op.get_mpi_communicator(),
2767  /*bool omit_zeroing_entries =*/true);
2768 
2769  // Duplicated from TrilinosWrappers::SparseMatrix::vmult
2770  const size_type i_local_size = i->end() - i->begin();
2771  AssertDimension(i_local_size, first_op_init_map.NumMyPoints());
2772  const Epetra_Map &second_op_init_map = second_op.OperatorRangeMap();
2773  AssertDimension(i_local_size, second_op_init_map.NumMyPoints());
2774  (void)second_op_init_map;
2775  Intermediate tril_int(View,
2776  first_op_init_map,
2777  const_cast<TrilinosScalar *>(i->begin()),
2778  i_local_size,
2779  1);
2780 
2781  // These operators may themselves be transposed or not, so we let them
2782  // decide what the intended outcome is
2783  second_op.ApplyInverse(tril_src, tril_int);
2784  first_op.ApplyInverse(tril_src, tril_dst);
2785  const int ierr = tril_dst.Update(1.0, tril_int, 1.0);
2786  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2787  };
2788 
2789  return_op.inv_Tvmult = [first_op, second_op](Range & tril_dst,
2790  const Domain &tril_src) {
2791  // Duplicated from LinearOperator::operator*
2792  // TODO: Template the constructor on GrowingVectorMemory vector type?
2793  GrowingVectorMemory<GVMVectorType> vector_memory;
2794  VectorMemory<GVMVectorType>::Pointer i(vector_memory);
2795 
2796  // These operators may themselves be transposed or not, so we let them
2797  // decide what the intended outcome is
2798  // We must first transpose the operators to get the right IndexSets
2799  // for the input, intermediate and result vectors
2800  const_cast<TrilinosPayload &>(first_op).transpose();
2801  const_cast<TrilinosPayload &>(second_op).transpose();
2802 
2803  // Initialize intermediate vector
2804  const Epetra_Map &first_op_init_map = first_op.OperatorDomainMap();
2805  i->reinit(IndexSet(first_op_init_map),
2806  first_op.get_mpi_communicator(),
2807  /*bool omit_zeroing_entries =*/true);
2808 
2809  // Duplicated from TrilinosWrappers::SparseMatrix::vmult
2810  const size_type i_local_size = i->end() - i->begin();
2811  AssertDimension(i_local_size, first_op_init_map.NumMyPoints());
2812  const Epetra_Map &second_op_init_map = second_op.OperatorDomainMap();
2813  AssertDimension(i_local_size, second_op_init_map.NumMyPoints());
2814  (void)second_op_init_map;
2815  Intermediate tril_int(View,
2816  first_op_init_map,
2817  const_cast<TrilinosScalar *>(i->begin()),
2818  i_local_size,
2819  1);
2820 
2821  // These operators may themselves be transposed or not, so we let them
2822  // decide what the intended outcome is
2823  second_op.ApplyInverse(tril_src, tril_int);
2824  first_op.ApplyInverse(tril_src, tril_dst);
2825  const int ierr = tril_dst.Update(1.0, tril_int, 1.0);
2826  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2827 
2828  // Reset transpose flag
2829  const_cast<TrilinosPayload &>(first_op).transpose();
2830  const_cast<TrilinosPayload &>(second_op).transpose();
2831  };
2832 
2833  return return_op;
2834  }
2835 
2836 
2837 
2838  TrilinosPayload
2839  operator*(const TrilinosPayload &first_op,
2840  const TrilinosPayload &second_op)
2841  {
2842  using Domain = typename TrilinosPayload::Domain;
2843  using Range = typename TrilinosPayload::Range;
2844  using Intermediate = typename TrilinosPayload::VectorType;
2845  using GVMVectorType = TrilinosWrappers::MPI::Vector;
2846 
2848  second_op.locally_owned_range_indices(),
2849  ExcMessage(
2850  "Operators are set to work on incompatible IndexSets."));
2851 
2852  TrilinosPayload return_op(first_op, second_op);
2853 
2854  // Capture by copy so the payloads are always valid
2855  return_op.vmult = [first_op, second_op](Range & tril_dst,
2856  const Domain &tril_src) {
2857  // Duplicated from LinearOperator::operator*
2858  // TODO: Template the constructor on GrowingVectorMemory vector type?
2859  GrowingVectorMemory<GVMVectorType> vector_memory;
2860  VectorMemory<GVMVectorType>::Pointer i(vector_memory);
2861 
2862  // Initialize intermediate vector
2863  const Epetra_Map &first_op_init_map = first_op.OperatorDomainMap();
2864  i->reinit(IndexSet(first_op_init_map),
2865  first_op.get_mpi_communicator(),
2866  /*bool omit_zeroing_entries =*/true);
2867 
2868  // Duplicated from TrilinosWrappers::SparseMatrix::vmult
2869  const size_type i_local_size = i->end() - i->begin();
2870  AssertDimension(i_local_size, first_op_init_map.NumMyPoints());
2871  const Epetra_Map &second_op_init_map = second_op.OperatorRangeMap();
2872  AssertDimension(i_local_size, second_op_init_map.NumMyPoints());
2873  (void)second_op_init_map;
2874  Intermediate tril_int(View,
2875  first_op_init_map,
2876  const_cast<TrilinosScalar *>(i->begin()),
2877  i_local_size,
2878  1);
2879 
2880  // These operators may themselves be transposed or not, so we let them
2881  // decide what the intended outcome is
2882  second_op.Apply(tril_src, tril_int);
2883  first_op.Apply(tril_int, tril_dst);
2884  };
2885 
2886  return_op.Tvmult = [first_op, second_op](Domain & tril_dst,
2887  const Range &tril_src) {
2888  // Duplicated from LinearOperator::operator*
2889  // TODO: Template the constructor on GrowingVectorMemory vector type?
2890  GrowingVectorMemory<GVMVectorType> vector_memory;
2891  VectorMemory<GVMVectorType>::Pointer i(vector_memory);
2892 
2893  // These operators may themselves be transposed or not, so we let them
2894  // decide what the intended outcome is
2895  // We must first transpose the operators to get the right IndexSets
2896  // for the input, intermediate and result vectors
2897  const_cast<TrilinosPayload &>(first_op).transpose();
2898  const_cast<TrilinosPayload &>(second_op).transpose();
2899 
2900  // Initialize intermediate vector
2901  const Epetra_Map &first_op_init_map = first_op.OperatorRangeMap();
2902  i->reinit(IndexSet(first_op_init_map),
2903  first_op.get_mpi_communicator(),
2904  /*bool omit_zeroing_entries =*/true);
2905 
2906  // Duplicated from TrilinosWrappers::SparseMatrix::vmult
2907  const size_type i_local_size = i->end() - i->begin();
2908  AssertDimension(i_local_size, first_op_init_map.NumMyPoints());
2909  const Epetra_Map &second_op_init_map = second_op.OperatorDomainMap();
2910  AssertDimension(i_local_size, second_op_init_map.NumMyPoints());
2911  (void)second_op_init_map;
2912  Intermediate tril_int(View,
2913  first_op_init_map,
2914  const_cast<TrilinosScalar *>(i->begin()),
2915  i_local_size,
2916  1);
2917 
2918  // Apply the operators in the reverse order to vmult
2919  first_op.Apply(tril_src, tril_int);
2920  second_op.Apply(tril_int, tril_dst);
2921 
2922  // Reset transpose flag
2923  const_cast<TrilinosPayload &>(first_op).transpose();
2924  const_cast<TrilinosPayload &>(second_op).transpose();
2925  };
2926 
2927  return_op.inv_vmult = [first_op, second_op](Domain & tril_dst,
2928  const Range &tril_src) {
2929  // Duplicated from LinearOperator::operator*
2930  // TODO: Template the constructor on GrowingVectorMemory vector type?
2931  GrowingVectorMemory<GVMVectorType> vector_memory;
2932  VectorMemory<GVMVectorType>::Pointer i(vector_memory);
2933 
2934  // Initialize intermediate vector
2935  const Epetra_Map &first_op_init_map = first_op.OperatorRangeMap();
2936  i->reinit(IndexSet(first_op_init_map),
2937  first_op.get_mpi_communicator(),
2938  /*bool omit_zeroing_entries =*/true);
2939 
2940  // Duplicated from TrilinosWrappers::SparseMatrix::vmult
2941  const size_type i_local_size = i->end() - i->begin();
2942  AssertDimension(i_local_size, first_op_init_map.NumMyPoints());
2943  const Epetra_Map &second_op_init_map = second_op.OperatorDomainMap();
2944  AssertDimension(i_local_size, second_op_init_map.NumMyPoints());
2945  (void)second_op_init_map;
2946  Intermediate tril_int(View,
2947  first_op_init_map,
2948  const_cast<TrilinosScalar *>(i->begin()),
2949  i_local_size,
2950  1);
2951 
2952  // Apply the operators in the reverse order to vmult
2953  // and the same order as Tvmult
2954  first_op.ApplyInverse(tril_src, tril_int);
2955  second_op.ApplyInverse(tril_int, tril_dst);
2956  };
2957 
2958  return_op.inv_Tvmult = [first_op, second_op](Range & tril_dst,
2959  const Domain &tril_src) {
2960  // Duplicated from LinearOperator::operator*
2961  // TODO: Template the constructor on GrowingVectorMemory vector type?
2962  GrowingVectorMemory<GVMVectorType> vector_memory;
2963  VectorMemory<GVMVectorType>::Pointer i(vector_memory);
2964 
2965  // These operators may themselves be transposed or not, so we let them
2966  // decide what the intended outcome is
2967  // We must first transpose the operators to get the right IndexSets
2968  // for the input, intermediate and result vectors
2969  const_cast<TrilinosPayload &>(first_op).transpose();
2970  const_cast<TrilinosPayload &>(second_op).transpose();
2971 
2972  // Initialize intermediate vector
2973  const Epetra_Map &first_op_init_map = first_op.OperatorDomainMap();
2974  i->reinit(IndexSet(first_op_init_map),
2975  first_op.get_mpi_communicator(),
2976  /*bool omit_zeroing_entries =*/true);
2977 
2978  // Duplicated from TrilinosWrappers::SparseMatrix::vmult
2979  const size_type i_local_size = i->end() - i->begin();
2980  AssertDimension(i_local_size, first_op_init_map.NumMyPoints());
2981  const Epetra_Map &second_op_init_map = second_op.OperatorRangeMap();
2982  AssertDimension(i_local_size, second_op_init_map.NumMyPoints());
2983  (void)second_op_init_map;
2984  Intermediate tril_int(View,
2985  first_op_init_map,
2986  const_cast<TrilinosScalar *>(i->begin()),
2987  i_local_size,
2988  1);
2989 
2990  // These operators may themselves be transposed or not, so we let them
2991  // decide what the intended outcome is
2992  // Apply the operators in the reverse order to Tvmult
2993  // and the same order as vmult
2994  second_op.ApplyInverse(tril_src, tril_int);
2995  first_op.ApplyInverse(tril_int, tril_dst);
2996 
2997  // Reset transpose flag
2998  const_cast<TrilinosPayload &>(first_op).transpose();
2999  const_cast<TrilinosPayload &>(second_op).transpose();
3000  };
3001 
3002  return return_op;
3003  }
3004 
3005  } // namespace LinearOperatorImplementation
3006  } /* namespace internal */
3007 } /* namespace TrilinosWrappers */
3008 
3009 
3010 
3011 // explicit instantiations
3012 # include "trilinos_sparse_matrix.inst"
3013 
3014 # ifndef DOXYGEN
3015 // TODO: put these instantiations into generic file
3016 namespace TrilinosWrappers
3017 {
3018  template void
3019  SparseMatrix::reinit(const ::SparsityPattern &);
3020 
3021  template void
3023 
3024  template void
3026  const IndexSet &,
3027  const ::SparsityPattern &,
3028  const MPI_Comm &,
3029  const bool);
3030 
3031  template void
3033  const IndexSet &,
3034  const DynamicSparsityPattern &,
3035  const MPI_Comm &,
3036  const bool);
3037 
3038  template void
3039  SparseMatrix::vmult(MPI::Vector &, const MPI::Vector &) const;
3040 
3041  template void
3043  const ::Vector<double> &) const;
3044 
3045  template void
3048  const ::LinearAlgebra::distributed::Vector<double> &) const;
3049 
3050 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
3051  template void
3054  const ::LinearAlgebra::TpetraWrappers::Vector<double> &) const;
3055 
3056  template void
3059  const ::LinearAlgebra::TpetraWrappers::Vector<float> &) const;
3060 # endif
3061 
3062  template void
3066 
3067  template void
3068  SparseMatrix::Tvmult(MPI::Vector &, const MPI::Vector &) const;
3069 
3070  template void
3072  const ::Vector<double> &) const;
3073 
3074  template void
3077  const ::LinearAlgebra::distributed::Vector<double> &) const;
3078 
3079 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
3080  template void
3083  const ::LinearAlgebra::TpetraWrappers::Vector<double> &) const;
3084 
3085  template void
3088  const ::LinearAlgebra::TpetraWrappers::Vector<float> &) const;
3089 # endif
3090 
3091  template void
3095 
3096  template void
3098 
3099  template void
3101  const ::Vector<double> &) const;
3102 
3103  template void
3106  const ::LinearAlgebra::distributed::Vector<double> &) const;
3107 
3108 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
3109  template void
3112  const ::LinearAlgebra::TpetraWrappers::Vector<double> &) const;
3113 
3114  template void
3117  const ::LinearAlgebra::TpetraWrappers::Vector<float> &) const;
3118 # endif
3119 
3120  template void
3124 
3125  template void
3127 
3128  template void
3130  const ::Vector<double> &) const;
3131 
3132  template void
3135  const ::LinearAlgebra::distributed::Vector<double> &) const;
3136 
3137 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
3138  template void
3141  const ::LinearAlgebra::TpetraWrappers::Vector<double> &) const;
3142 
3143  template void
3146  const ::LinearAlgebra::TpetraWrappers::Vector<float> &) const;
3147 # endif
3148 
3149  template void
3153 } // namespace TrilinosWrappers
3154 # endif // DOXYGEN
3155 
3157 
3158 #endif // DEAL_II_WITH_TRILINOS
const IndexSet & row_index_set() const
size_type row_length(const size_type row) const
size_type column_number(const size_type row, const size_type index) const
size_type size() const
Definition: index_set.h:1626
bool is_element(const size_type index) const
Definition: index_set.h:1734
Epetra_Map make_trilinos_map(const MPI_Comm &communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition: index_set.cc:777
size_type n_rows() const
size_type n_cols() const
void reinit(const size_type m, const size_type n, const ArrayView< const unsigned int > &row_lengths)
const Epetra_BlockMap & trilinos_partitioner() const
void reinit(const Vector &v, const bool omit_zeroing_entries=false, const bool allow_different_maps=false)
std::shared_ptr< std::vector< TrilinosScalar > > value_cache
std::shared_ptr< std::vector< size_type > > colnum_cache
void set(const size_type i, const size_type j, const TrilinosScalar value)
std::unique_ptr< Epetra_Map > column_space_map
SparseMatrix & operator*=(const TrilinosScalar factor)
void mmult(SparseMatrix &C, const SparseMatrix &B, const MPI::Vector &V=MPI::Vector()) const
std::unique_ptr< Epetra_Export > nonlocal_matrix_exporter
std::unique_ptr< Epetra_FECrsMatrix > matrix
const Epetra_CrsMatrix & trilinos_matrix() const
::types::global_dof_index size_type
TrilinosScalar matrix_norm_square(const MPI::Vector &v) const
IndexSet locally_owned_range_indices() const
void print(std::ostream &out, const bool write_extended_trilinos_info=false) const
void Tmmult(SparseMatrix &C, const SparseMatrix &B, const MPI::Vector &V=MPI::Vector()) const
void clear_row(const size_type row, const TrilinosScalar new_diag_value=0)
std::pair< size_type, size_type > local_range() const
const Epetra_CrsGraph & trilinos_sparsity_pattern() const
void reinit(const SparsityPatternType &sparsity_pattern)
std::enable_if_t< std::is_same< typename VectorType::value_type, TrilinosScalar >::value > Tvmult(VectorType &dst, const VectorType &src) const
void vmult_add(VectorType &dst, const VectorType &src) const
void compress(::VectorOperation::values operation)
SparseMatrix & operator/=(const TrilinosScalar factor)
void Tvmult_add(VectorType &dst, const VectorType &src) const
void clear_rows(const std::vector< size_type > &rows, const TrilinosScalar new_diag_value=0)
TrilinosScalar el(const size_type i, const size_type j) const
IndexSet locally_owned_domain_indices() const
TrilinosScalar residual(MPI::Vector &dst, const MPI::Vector &x, const MPI::Vector &b) const
bool in_local_range(const size_type index) const
void copy_from(const SparseMatrix &source)
unsigned int row_length(const size_type row) const
std::uint64_t n_nonzero_elements() const
TrilinosScalar diag_element(const size_type i) const
void add(const size_type i, const size_type j, const TrilinosScalar value)
unsigned int local_size() const
TrilinosScalar operator()(const size_type i, const size_type j) const
TrilinosScalar matrix_scalar_product(const MPI::Vector &u, const MPI::Vector &v) const
std::enable_if_t< std::is_same< typename VectorType::value_type, TrilinosScalar >::value > vmult(VectorType &dst, const VectorType &src) const
SparseMatrix & operator=(const SparseMatrix &)=delete
std::unique_ptr< Epetra_CrsMatrix > nonlocal_matrix
std::unique_ptr< Epetra_CrsGraph > nonlocal_graph
const Epetra_Map & domain_partitioner() const
const Epetra_FECrsGraph & trilinos_sparsity_pattern() const
std::function< void(VectorType &, const VectorType &)> inv_Tvmult
virtual int ApplyInverse(const VectorType &Y, VectorType &X) const override
virtual int Apply(const VectorType &X, VectorType &Y) const override
std::function< void(VectorType &, const VectorType &)> inv_vmult
Definition: vector.h:109
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:458
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Definition: config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:459
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
Definition: config.h:512
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Point< 2 > second
Definition: grid_out.cc:4606
Point< 2 > first
Definition: grid_out.cc:4605
static ::ExceptionBase & ExcTrilinosError(int arg1)
static ::ExceptionBase & ExcAccessToNonlocalRow(std::size_t arg1)
static ::ExceptionBase & ExcDivideByZero()
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcInvalidIndex(size_type arg1, size_type arg2)
static ::ExceptionBase & ExcSourceEqualsDestination()
static ::ExceptionBase & ExcScalarAssignmentOnlyForZeroValue()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1501
static ::ExceptionBase & ExcNotImplemented()
#define AssertIsFinite(number)
Definition: exceptions.h:1786
static ::ExceptionBase & ExcTrilinosError(int arg1)
static ::ExceptionBase & ExcAccessToNonLocalElement(size_type arg1, size_type arg2, size_type arg3, size_type arg4)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1695
#define AssertIndexRange(index, range)
Definition: exceptions.h:1760
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcMatrixNotCompressed()
static ::ExceptionBase & ExcAccessToNonPresentElement(size_type arg1, size_type arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcNotQuadratic()
#define AssertThrow(cond, exc)
Definition: exceptions.h:1611
IndexSet complete_index_set(const IndexSet::size_type N)
Definition: index_set.h:1054
Expression fabs(const Expression &x)
typename ::internal::FEInterfaceViews::ViewType< dim, spacedim, Extractor >::type View
@ matrix
Contents is actually a matrix.
static const char V
void swap(MemorySpaceData< Number, MemorySpace > &, MemorySpaceData< Number, MemorySpace > &)
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
TrilinosPayload operator+(const TrilinosPayload &first_op, const TrilinosPayload &second_op)
TrilinosPayload operator*(const TrilinosPayload &first_op, const TrilinosPayload &second_op)
void check_vector_map_equality(const Epetra_CrsMatrix &, const VectorType &, const VectorType &)
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(VectorType &V)
::types::global_dof_index size_type
void perform_mmult(const SparseMatrix &inputleft, const SparseMatrix &inputright, SparseMatrix &result, const MPI::Vector &V, const bool transpose_left)
TrilinosWrappers::types::int_type min_my_gid(const Epetra_BlockMap &map)
TrilinosWrappers::types::int64_type n_global_elements(const Epetra_BlockMap &map)
TrilinosWrappers::types::int_type global_column_index(const Epetra_CrsMatrix &matrix, const ::types::global_dof_index i)
TrilinosWrappers::types::int_type max_my_gid(const Epetra_BlockMap &map)
TrilinosWrappers::types::int_type global_row_index(const Epetra_CrsMatrix &matrix, const ::types::global_dof_index i)
TrilinosWrappers::types::int_type n_global_cols(const Epetra_CrsGraph &graph)
const Epetra_Comm & comm_self()
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition: utilities.h:999
Definition: types.h:32
unsigned int global_dof_index
Definition: types.h:81