1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2020 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 #include "libmesh/auto_ptr.h" // libmesh_make_unique
33 #include "libmesh/enum_solver_package.h"
34
35 namespace libMesh
36 {
37
38
39 //------------------------------------------------------------------
40 // SparseMatrix Methods
41
42
43 // Constructor
44 template <typename T>
SparseMatrix(const Parallel::Communicator & comm_in)45 SparseMatrix<T>::SparseMatrix (const Parallel::Communicator & comm_in) :
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>
add_block_matrix(const DenseMatrix<T> & dm,const std::vector<numeric_index_type> & brows,const std::vector<numeric_index_type> & bcols)56 void SparseMatrix<T>::add_block_matrix (const DenseMatrix<T> & dm,
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 (auto & row : brows)
74 {
75 numeric_index_type i = row * blocksize;
76
77 for (unsigned int v=0; v<blocksize; v++)
78 rows.push_back(i++);
79 }
80
81 for (auto & col : bcols)
82 {
83 numeric_index_type j = col * 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 <>
print(std::ostream & os,const bool sparse)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 (auto i : make_range(this->m()))
107 {
108 for (auto j : make_range(this->n()))
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 (auto i : make_range(this->m()))
115 {
116 for (auto j : make_range(this->n()))
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>>
build(const Parallel::Communicator & comm,const SolverPackage solver_package)130 SparseMatrix<T>::build(const Parallel::Communicator & comm,
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>
vector_mult(NumericVector<T> & dest,const NumericVector<T> & arg)170 void SparseMatrix<T>::vector_mult (NumericVector<T> & dest,
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>
vector_mult_add(NumericVector<T> & dest,const NumericVector<T> & arg)180 void SparseMatrix<T>::vector_mult_add (NumericVector<T> & dest,
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>
zero_rows(std::vector<numeric_index_type> &,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>
print(std::ostream & os,const bool sparse)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 libmesh_error_msg_if(!this->_dof_map, "Error! Trying to print a matrix with no dof_map set!");
207
208 // We'll print the matrix from processor 0 to make sure
209 // it's serialized properly
210 if (this->processor_id() == 0)
211 {
212 libmesh_assert_equal_to (this->_dof_map->first_dof(), 0);
213 for (numeric_index_type i=this->_dof_map->first_dof();
214 i!=this->_dof_map->end_dof(); ++i)
215 {
216 if (sparse)
217 {
218 for (auto j : make_range(this->n()))
219 {
220 T c = (*this)(i,j);
221 if (c != static_cast<T>(0.0))
222 {
223 os << i << " " << j << " " << c << std::endl;
224 }
225 }
226 }
227 else
228 {
229 for (auto j : make_range(this->n()))
230 os << (*this)(i,j) << " ";
231 os << std::endl;
232 }
233 }
234
235 std::vector<numeric_index_type> ibuf, jbuf;
236 std::vector<T> cbuf;
237 numeric_index_type currenti = this->_dof_map->end_dof();
238 for (auto p : IntRange<processor_id_type>(1, this->n_processors()))
239 {
240 this->comm().receive(p, ibuf);
241 this->comm().receive(p, jbuf);
242 this->comm().receive(p, cbuf);
243 libmesh_assert_equal_to (ibuf.size(), jbuf.size());
244 libmesh_assert_equal_to (ibuf.size(), cbuf.size());
245
246 if (ibuf.empty())
247 continue;
248 libmesh_assert_greater_equal (ibuf.front(), currenti);
249 libmesh_assert_greater_equal (ibuf.back(), ibuf.front());
250
251 std::size_t currentb = 0;
252 for (;currenti <= ibuf.back(); ++currenti)
253 {
254 if (sparse)
255 {
256 for (numeric_index_type j=0; j<this->n(); j++)
257 {
258 if (currentb < ibuf.size() &&
259 ibuf[currentb] == currenti &&
260 jbuf[currentb] == j)
261 {
262 os << currenti << " " << j << " " << cbuf[currentb] << std::endl;
263 currentb++;
264 }
265 }
266 }
267 else
268 {
269 for (auto j : make_range(this->n()))
270 {
271 if (currentb < ibuf.size() &&
272 ibuf[currentb] == currenti &&
273 jbuf[currentb] == j)
274 {
275 os << cbuf[currentb] << " ";
276 currentb++;
277 }
278 else
279 os << static_cast<T>(0.0) << " ";
280 }
281 os << std::endl;
282 }
283 }
284 }
285 if (!sparse)
286 {
287 for (; currenti != this->m(); ++currenti)
288 {
289 for (numeric_index_type j=0; j<this->n(); j++)
290 os << static_cast<T>(0.0) << " ";
291 os << std::endl;
292 }
293 }
294 }
295 else
296 {
297 std::vector<numeric_index_type> ibuf, jbuf;
298 std::vector<T> cbuf;
299
300 // We'll assume each processor has access to entire
301 // matrix rows, so (*this)(i,j) is valid if i is a local index.
302 for (numeric_index_type i=this->_dof_map->first_dof();
303 i!=this->_dof_map->end_dof(); ++i)
304 {
305 for (auto j : make_range(this->n()))
306 {
307 T c = (*this)(i,j);
308 if (c != static_cast<T>(0.0))
309 {
310 ibuf.push_back(i);
311 jbuf.push_back(j);
312 cbuf.push_back(c);
313 }
314 }
315 }
316 this->comm().send(0,ibuf);
317 this->comm().send(0,jbuf);
318 this->comm().send(0,cbuf);
319 }
320 }
321
322
323
324 //------------------------------------------------------------------
325 // Explicit instantiations
326 template class SparseMatrix<Number>;
327
328 } // namespace libMesh
329