libMesh
sparse_matrix.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 // C++ includes
21 
22 // Local Includes
23 #include "libmesh/dof_map.h"
24 #include "libmesh/dense_matrix.h"
25 #include "libmesh/laspack_matrix.h"
26 #include "libmesh/eigen_sparse_matrix.h"
27 #include "libmesh/parallel.h"
28 #include "libmesh/petsc_matrix.h"
29 #include "libmesh/sparse_matrix.h"
30 #include "libmesh/trilinos_epetra_matrix.h"
31 #include "libmesh/numeric_vector.h"
32 
33 namespace libMesh
34 {
35 
36 
37 //------------------------------------------------------------------
38 // SparseMatrix Methods
39 
40 
41 // Constructor
42 template <typename T>
44  ParallelObject(comm_in),
45  _dof_map(libmesh_nullptr),
46  _is_initialized(false)
47 {}
48 
49 
50 
51 // Destructor
52 template <typename T>
54 {}
55 
56 
57 
58 
59 // default implementation is to fall back to non-blocked method
60 template <typename T>
62  const std::vector<numeric_index_type> & brows,
63  const std::vector<numeric_index_type> & bcols)
64 {
65  libmesh_assert_equal_to (dm.m() / brows.size(), dm.n() / bcols.size());
66 
67  const numeric_index_type blocksize = cast_int<numeric_index_type>
68  (dm.m() / brows.size());
69 
70  libmesh_assert_equal_to (dm.m()%blocksize, 0);
71  libmesh_assert_equal_to (dm.n()%blocksize, 0);
72 
73  std::vector<numeric_index_type> rows, cols;
74 
75  rows.reserve(blocksize*brows.size());
76  cols.reserve(blocksize*bcols.size());
77 
78  for (std::size_t ib=0; ib<brows.size(); ib++)
79  {
80  numeric_index_type i=brows[ib]*blocksize;
81 
82  for (unsigned int v=0; v<blocksize; v++)
83  rows.push_back(i++);
84  }
85 
86  for (std::size_t jb=0; jb<bcols.size(); jb++)
87  {
88  numeric_index_type j=bcols[jb]*blocksize;
89 
90  for (unsigned int v=0; v<blocksize; v++)
91  cols.push_back(j++);
92  }
93 
94  this->add_matrix (dm, rows, cols);
95 }
96 
97 
98 
99 // Full specialization of print method for Complex datatypes
100 template <>
101 void SparseMatrix<Complex>::print(std::ostream & os, const bool sparse) const
102 {
103  // std::complex<>::operator<<() is defined, but use this form
104 
105  if (sparse)
106  {
107  libmesh_not_implemented();
108  }
109 
110  os << "Real part:" << std::endl;
111  for (numeric_index_type i=0; i<this->m(); i++)
112  {
113  for (numeric_index_type j=0; j<this->n(); j++)
114  os << std::setw(8) << (*this)(i,j).real() << " ";
115  os << std::endl;
116  }
117 
118  os << std::endl << "Imaginary part:" << std::endl;
119  for (numeric_index_type i=0; i<this->m(); i++)
120  {
121  for (numeric_index_type j=0; j<this->n(); j++)
122  os << std::setw(8) << (*this)(i,j).imag() << " ";
123  os << std::endl;
124  }
125 }
126 
127 
128 
129 
130 
131 
132 // Full specialization for Real datatypes
133 template <typename T>
136  const SolverPackage solver_package)
137 {
138  // Avoid unused parameter warnings when no solver packages are enabled.
139  libmesh_ignore(comm);
140 
141  // Build the appropriate vector
142  switch (solver_package)
143  {
144 
145 #ifdef LIBMESH_HAVE_LASPACK
146  case LASPACK_SOLVERS:
147  return UniquePtr<SparseMatrix<T>>(new LaspackMatrix<T>(comm));
148 #endif
149 
150 
151 #ifdef LIBMESH_HAVE_PETSC
152  case PETSC_SOLVERS:
154 #endif
155 
156 
157 #ifdef LIBMESH_TRILINOS_HAVE_EPETRA
158  case TRILINOS_SOLVERS:
159  return UniquePtr<SparseMatrix<T>>(new EpetraMatrix<T>(comm));
160 #endif
161 
162 
163 #ifdef LIBMESH_HAVE_EIGEN
164  case EIGEN_SOLVERS:
166 #endif
167 
168  default:
169  libmesh_error_msg("ERROR: Unrecognized solver package: " << solver_package);
170  }
171 
172  libmesh_error_msg("We'll never get here!");
173  return UniquePtr<SparseMatrix<T>>();
174 }
175 
176 
177 template <typename T>
179  const NumericVector<T> & arg) const
180 {
181  dest.zero();
182  this->vector_mult_add(dest,arg);
183 }
184 
185 
186 
187 template <typename T>
189  const NumericVector<T> & arg) const
190 {
191  /* This functionality is actually implemented in the \p
192  NumericVector class. */
193  dest.add_vector(arg,*this);
194 }
195 
196 
197 
198 template <typename T>
199 void SparseMatrix<T>::zero_rows (std::vector<numeric_index_type> &, T)
200 {
201  /* This functionality isn't implemented or stubbed in every subclass yet */
202  libmesh_not_implemented();
203 }
204 
205 
206 
207 template <typename T>
208 void SparseMatrix<T>::print(std::ostream & os, const bool sparse) const
209 {
210  parallel_object_only();
211 
212  libmesh_assert (this->initialized());
213 
214  if (!this->_dof_map)
215  libmesh_error_msg("Error! Trying to print a matrix with no dof_map set!");
216 
217  // We'll print the matrix from processor 0 to make sure
218  // it's serialized properly
219  if (this->processor_id() == 0)
220  {
221  libmesh_assert_equal_to (this->_dof_map->first_dof(), 0);
222  for (numeric_index_type i=this->_dof_map->first_dof();
223  i!=this->_dof_map->end_dof(); ++i)
224  {
225  if (sparse)
226  {
227  for (numeric_index_type j=0; j<this->n(); j++)
228  {
229  T c = (*this)(i,j);
230  if (c != static_cast<T>(0.0))
231  {
232  os << i << " " << j << " " << c << std::endl;
233  }
234  }
235  }
236  else
237  {
238  for (numeric_index_type j=0; j<this->n(); j++)
239  os << (*this)(i,j) << " ";
240  os << std::endl;
241  }
242  }
243 
244  std::vector<numeric_index_type> ibuf, jbuf;
245  std::vector<T> cbuf;
246  numeric_index_type currenti = this->_dof_map->end_dof();
247  for (processor_id_type p=1; p < this->n_processors(); ++p)
248  {
249  this->comm().receive(p, ibuf);
250  this->comm().receive(p, jbuf);
251  this->comm().receive(p, cbuf);
252  libmesh_assert_equal_to (ibuf.size(), jbuf.size());
253  libmesh_assert_equal_to (ibuf.size(), cbuf.size());
254 
255  if (ibuf.empty())
256  continue;
257  libmesh_assert_greater_equal (ibuf.front(), currenti);
258  libmesh_assert_greater_equal (ibuf.back(), ibuf.front());
259 
260  std::size_t currentb = 0;
261  for (;currenti <= ibuf.back(); ++currenti)
262  {
263  if (sparse)
264  {
265  for (numeric_index_type j=0; j<this->n(); j++)
266  {
267  if (currentb < ibuf.size() &&
268  ibuf[currentb] == currenti &&
269  jbuf[currentb] == j)
270  {
271  os << currenti << " " << j << " " << cbuf[currentb] << std::endl;
272  currentb++;
273  }
274  }
275  }
276  else
277  {
278  for (numeric_index_type j=0; j<this->n(); j++)
279  {
280  if (currentb < ibuf.size() &&
281  ibuf[currentb] == currenti &&
282  jbuf[currentb] == j)
283  {
284  os << cbuf[currentb] << " ";
285  currentb++;
286  }
287  else
288  os << static_cast<T>(0.0) << " ";
289  }
290  os << std::endl;
291  }
292  }
293  }
294  if (!sparse)
295  {
296  for (; currenti != this->m(); ++currenti)
297  {
298  for (numeric_index_type j=0; j<this->n(); j++)
299  os << static_cast<T>(0.0) << " ";
300  os << std::endl;
301  }
302  }
303  }
304  else
305  {
306  std::vector<numeric_index_type> ibuf, jbuf;
307  std::vector<T> cbuf;
308 
309  // We'll assume each processor has access to entire
310  // matrix rows, so (*this)(i,j) is valid if i is a local index.
311  for (numeric_index_type i=this->_dof_map->first_dof();
312  i!=this->_dof_map->end_dof(); ++i)
313  {
314  for (numeric_index_type j=0; j<this->n(); j++)
315  {
316  T c = (*this)(i,j);
317  if (c != static_cast<T>(0.0))
318  {
319  ibuf.push_back(i);
320  jbuf.push_back(j);
321  cbuf.push_back(c);
322  }
323  }
324  }
325  this->comm().send(0,ibuf);
326  this->comm().send(0,jbuf);
327  this->comm().send(0,cbuf);
328  }
329 }
330 
331 
332 
333 //------------------------------------------------------------------
334 // Explicit instantiations
335 template class SparseMatrix<Number>;
336 
337 } // namespace libMesh
virtual void add_block_matrix(const DenseMatrix< T > &dm, const std::vector< numeric_index_type > &brows, const std::vector< numeric_index_type > &bcols)
Add the full matrix dm to the SparseMatrix.
Definition: sparse_matrix.C:61
virtual numeric_index_type n() const =0
Encapsulates the MPI_Comm object.
Definition: parallel.h:657
unsigned int n() const
virtual void add_vector(const T *v, const std::vector< numeric_index_type > &dof_indices)
Computes , where v is a pointer and each dof_indices[i] specifies where to add value v[i]...
virtual numeric_index_type m() const =0
SolverPackage
Defines an enum for various linear solver packages.
processor_id_type n_processors() const
uint8_t processor_id_type
Definition: id_types.h:99
void vector_mult(NumericVector< T > &dest, const NumericVector< T > &arg) const
Multiplies the matrix by the NumericVector arg and stores the result in NumericVector dest...
const class libmesh_nullptr_t libmesh_nullptr
Numeric vector.
Definition: dof_map.h:66
void vector_mult_add(NumericVector< T > &dest, const NumericVector< T > &arg) const
Multiplies the matrix by the NumericVector arg and adds the result to the NumericVector dest...
static UniquePtr< SparseMatrix< T > > build(const Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package())
Builds a SparseMatrix<T> using the linear solver package specified by solver_package.
The libMesh namespace provides an interface to certain functionality in the library.
virtual void zero()=0
Set all entries to zero.
libmesh_assert(j)
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
virtual bool initialized() const
dof_id_type first_dof(const processor_id_type proc) const
Definition: dof_map.h:535
dof_id_type end_dof(const processor_id_type proc) const
Definition: dof_map.h:577
virtual void add_matrix(const DenseMatrix< T > &dm, const std::vector< numeric_index_type > &rows, const std::vector< numeric_index_type > &cols)=0
Add the full matrix dm to the SparseMatrix.
dof_id_type numeric_index_type
Definition: id_types.h:92
bool _is_initialized
Flag that tells if init() has been called.
Definition: libmesh.C:255
unsigned int m() const
virtual void zero_rows(std::vector< numeric_index_type > &rows, T diag_value=0.0)
Sets all row entries to 0 then puts diag_value in the diagonal entry.
void send(const unsigned int dest_processor_id, const T &buf, const MessageTag &tag=no_tag) const
Blocking-send to one processor with data-defined type.
This class forms the base class for all other classes that are expected to be implemented in parallel...
void libmesh_ignore(const T &)
DofMap const * _dof_map
The DofMap object associated with this object.
SparseMatrix(const Parallel::Communicator &comm LIBMESH_CAN_DEFAULT_TO_COMMWORLD)
Constructor; initializes the matrix to be empty, without any structure, i.e.
Definition: sparse_matrix.C:43
const Parallel::Communicator & comm() const
void print(std::ostream &os=libMesh::out, const bool sparse=false) const
Print the contents of the matrix to the screen in a uniform style, regardless of matrix/solver packag...
Status receive(const unsigned int dest_processor_id, T &buf, const MessageTag &tag=any_tag) const
Blocking-receive from one processor with data-defined type.
Defines a dense matrix for use in Finite Element-type computations.
Definition: dof_map.h:64
virtual ~SparseMatrix()
Destructor.
Definition: sparse_matrix.C:53
processor_id_type processor_id() const