sparse_matrix.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 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"
27 #include "libmesh/parallel.h"
28 #include "libmesh/petsc_matrix.h"
29 #include "libmesh/sparse_matrix.h"
31 #include "libmesh/numeric_vector.h"
32 #include "libmesh/auto_ptr.h" // libmesh_make_unique
34 
35 namespace libMesh
36 {
37 
38 
39 //------------------------------------------------------------------
40 // SparseMatrix Methods
41 
42 
43 // Constructor
44 template <typename T>
46  ParallelObject(comm_in),
47  _dof_map(nullptr),
48  _is_initialized(false)
49 {}
50 
51 
52 
53 
54 // default implementation is to fall back to non-blocked method
55 template <typename T>
57  const std::vector<numeric_index_type> & brows,
58  const std::vector<numeric_index_type> & bcols)
59 {
60  libmesh_assert_equal_to (dm.m() / brows.size(), dm.n() / bcols.size());
61 
62  const numeric_index_type blocksize = cast_int<numeric_index_type>
63  (dm.m() / brows.size());
64 
65  libmesh_assert_equal_to (dm.m()%blocksize, 0);
66  libmesh_assert_equal_to (dm.n()%blocksize, 0);
67 
68  std::vector<numeric_index_type> rows, cols;
69 
70  rows.reserve(blocksize*brows.size());
71  cols.reserve(blocksize*bcols.size());
72 
73  for (std::size_t ib=0; ib<brows.size(); ib++)
74  {
75  numeric_index_type i=brows[ib]*blocksize;
76 
77  for (unsigned int v=0; v<blocksize; v++)
78  rows.push_back(i++);
79  }
80 
81  for (std::size_t jb=0; jb<bcols.size(); jb++)
82  {
83  numeric_index_type j=bcols[jb]*blocksize;
84 
85  for (unsigned int v=0; v<blocksize; v++)
86  cols.push_back(j++);
87  }
88 
89  this->add_matrix (dm, rows, cols);
90 }
91 
92 
93 
94 // Full specialization of print method for Complex datatypes
95 template <>
96 void SparseMatrix<Complex>::print(std::ostream & os, const bool sparse) const
97 {
98  // std::complex<>::operator<<() is defined, but use this form
99 
100  if (sparse)
101  {
102  libmesh_not_implemented();
103  }
104 
105  os << "Real part:" << std::endl;
106  for (numeric_index_type i=0; i<this->m(); i++)
107  {
108  for (numeric_index_type j=0; j<this->n(); j++)
109  os << std::setw(8) << (*this)(i,j).real() << " ";
110  os << std::endl;
111  }
112 
113  os << std::endl << "Imaginary part:" << std::endl;
114  for (numeric_index_type i=0; i<this->m(); i++)
115  {
116  for (numeric_index_type j=0; j<this->n(); j++)
117  os << std::setw(8) << (*this)(i,j).imag() << " ";
118  os << std::endl;
119  }
120 }
121 
122 
123 
124 
125 
126 
127 // Full specialization for Real datatypes
128 template <typename T>
129 std::unique_ptr<SparseMatrix<T>>
131  const SolverPackage solver_package)
132 {
133  // Avoid unused parameter warnings when no solver packages are enabled.
134  libmesh_ignore(comm);
135 
136  // Build the appropriate vector
137  switch (solver_package)
138  {
139 
140 #ifdef LIBMESH_HAVE_LASPACK
141  case LASPACK_SOLVERS:
142  return libmesh_make_unique<LaspackMatrix<T>>(comm);
143 #endif
144 
145 
146 #ifdef LIBMESH_HAVE_PETSC
147  case PETSC_SOLVERS:
148  return libmesh_make_unique<PetscMatrix<T>>(comm);
149 #endif
150 
151 
152 #ifdef LIBMESH_TRILINOS_HAVE_EPETRA
153  case TRILINOS_SOLVERS:
154  return libmesh_make_unique<EpetraMatrix<T>>(comm);
155 #endif
156 
157 
158 #ifdef LIBMESH_HAVE_EIGEN
159  case EIGEN_SOLVERS:
160  return libmesh_make_unique<EigenSparseMatrix<T>>(comm);
161 #endif
162 
163  default:
164  libmesh_error_msg("ERROR: Unrecognized solver package: " << solver_package);
165  }
166 }
167 
168 
169 template <typename T>
171  const NumericVector<T> & arg) const
172 {
173  dest.zero();
174  this->vector_mult_add(dest,arg);
175 }
176 
177 
178 
179 template <typename T>
181  const NumericVector<T> & arg) const
182 {
183  /* This functionality is actually implemented in the \p
184  NumericVector class. */
185  dest.add_vector(arg,*this);
186 }
187 
188 
189 
190 template <typename T>
191 void SparseMatrix<T>::zero_rows (std::vector<numeric_index_type> &, T)
192 {
193  /* This functionality isn't implemented or stubbed in every subclass yet */
194  libmesh_not_implemented();
195 }
196 
197 
198 
199 template <typename T>
200 void SparseMatrix<T>::print(std::ostream & os, const bool sparse) const
201 {
202  parallel_object_only();
203 
204  libmesh_assert (this->initialized());
205 
206  if (!this->_dof_map)
207  libmesh_error_msg("Error! Trying to print a matrix with no dof_map set!");
208 
209  // We'll print the matrix from processor 0 to make sure
210  // it's serialized properly
211  if (this->processor_id() == 0)
212  {
213  libmesh_assert_equal_to (this->_dof_map->first_dof(), 0);
214  for (numeric_index_type i=this->_dof_map->first_dof();
215  i!=this->_dof_map->end_dof(); ++i)
216  {
217  if (sparse)
218  {
219  for (numeric_index_type j=0; j<this->n(); j++)
220  {
221  T c = (*this)(i,j);
222  if (c != static_cast<T>(0.0))
223  {
224  os << i << " " << j << " " << c << std::endl;
225  }
226  }
227  }
228  else
229  {
230  for (numeric_index_type j=0; j<this->n(); j++)
231  os << (*this)(i,j) << " ";
232  os << std::endl;
233  }
234  }
235 
236  std::vector<numeric_index_type> ibuf, jbuf;
237  std::vector<T> cbuf;
238  numeric_index_type currenti = this->_dof_map->end_dof();
239  for (processor_id_type p=1; p < this->n_processors(); ++p)
240  {
241  this->comm().receive(p, ibuf);
242  this->comm().receive(p, jbuf);
243  this->comm().receive(p, cbuf);
244  libmesh_assert_equal_to (ibuf.size(), jbuf.size());
245  libmesh_assert_equal_to (ibuf.size(), cbuf.size());
246 
247  if (ibuf.empty())
248  continue;
249  libmesh_assert_greater_equal (ibuf.front(), currenti);
250  libmesh_assert_greater_equal (ibuf.back(), ibuf.front());
251 
252  std::size_t currentb = 0;
253  for (;currenti <= ibuf.back(); ++currenti)
254  {
255  if (sparse)
256  {
257  for (numeric_index_type j=0; j<this->n(); j++)
258  {
259  if (currentb < ibuf.size() &&
260  ibuf[currentb] == currenti &&
261  jbuf[currentb] == j)
262  {
263  os << currenti << " " << j << " " << cbuf[currentb] << std::endl;
264  currentb++;
265  }
266  }
267  }
268  else
269  {
270  for (numeric_index_type j=0; j<this->n(); j++)
271  {
272  if (currentb < ibuf.size() &&
273  ibuf[currentb] == currenti &&
274  jbuf[currentb] == j)
275  {
276  os << cbuf[currentb] << " ";
277  currentb++;
278  }
279  else
280  os << static_cast<T>(0.0) << " ";
281  }
282  os << std::endl;
283  }
284  }
285  }
286  if (!sparse)
287  {
288  for (; currenti != this->m(); ++currenti)
289  {
290  for (numeric_index_type j=0; j<this->n(); j++)
291  os << static_cast<T>(0.0) << " ";
292  os << std::endl;
293  }
294  }
295  }
296  else
297  {
298  std::vector<numeric_index_type> ibuf, jbuf;
299  std::vector<T> cbuf;
300 
301  // We'll assume each processor has access to entire
302  // matrix rows, so (*this)(i,j) is valid if i is a local index.
303  for (numeric_index_type i=this->_dof_map->first_dof();
304  i!=this->_dof_map->end_dof(); ++i)
305  {
306  for (numeric_index_type j=0; j<this->n(); j++)
307  {
308  T c = (*this)(i,j);
309  if (c != static_cast<T>(0.0))
310  {
311  ibuf.push_back(i);
312  jbuf.push_back(j);
313  cbuf.push_back(c);
314  }
315  }
316  }
317  this->comm().send(0,ibuf);
318  this->comm().send(0,jbuf);
319  this->comm().send(0,cbuf);
320  }
321 }
322 
323 
324 
325 //------------------------------------------------------------------
326 // Explicit instantiations
327 template class SparseMatrix<Number>;
328 
329 } // namespace libMesh
SparseMatrix(const Parallel::Communicator &comm)
Definition: sparse_matrix.C:45
EIGEN_SOLVERS
Definition: libmesh.C:246
TRILINOS_SOLVERS
Definition: libmesh.C:244
virtual void add_vector(const T *v, const std::vector< numeric_index_type > &dof_indices)
uint8_t processor_id_type
Definition: id_types.h:99
Provides a uniform interface to vector storage schemes for different linear algebra libraries...
Definition: diff_context.h:40
unsigned int m() const
virtual void zero()=0
void libmesh_ignore(const Args &...)
dof_id_type numeric_index_type
Definition: id_types.h:92
LASPACK_SOLVERS
Definition: libmesh.C:248
An object whose state is distributed along a set of processors.
bool initialized()
Definition: libmesh.C:258
unsigned int n() const
A matrix object used for finite element assembly and numerics.
Definition: dense_matrix.h:54