sparse_matrix.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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 
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(NULL),
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 = dm.m() / brows.size();
68 
69  libmesh_assert_equal_to (dm.m()%blocksize, 0);
70  libmesh_assert_equal_to (dm.n()%blocksize, 0);
71 
72  std::vector<numeric_index_type> rows, cols;
73 
74  rows.reserve(blocksize*brows.size());
75  cols.reserve(blocksize*bcols.size());
76 
77  for (unsigned int ib=0; ib<brows.size(); ib++)
78  {
79  numeric_index_type i=brows[ib]*blocksize;
80 
81  for (unsigned int v=0; v<blocksize; v++)
82  rows.push_back(i++);
83  }
84 
85  for (unsigned int jb=0; jb<bcols.size(); jb++)
86  {
87  numeric_index_type j=bcols[jb]*blocksize;
88 
89  for (unsigned int v=0; v<blocksize; v++)
90  cols.push_back(j++);
91  }
92 
93  this->add_matrix (dm, rows, cols);
94 }
95 
96 
97 
98 // Full specialization of print method for Complex datatypes
99 template <>
100 void SparseMatrix<Complex>::print(std::ostream& os, const bool sparse) const
101 {
102  // std::complex<>::operator<<() is defined, but use this form
103 
104  if(sparse)
105  {
106  libmesh_not_implemented();
107  }
108 
109  os << "Real part:" << std::endl;
110  for (numeric_index_type i=0; i<this->m(); i++)
111  {
112  for (numeric_index_type j=0; j<this->n(); j++)
113  os << std::setw(8) << (*this)(i,j).real() << " ";
114  os << std::endl;
115  }
116 
117  os << std::endl << "Imaginary part:" << std::endl;
118  for (numeric_index_type i=0; i<this->m(); i++)
119  {
120  for (numeric_index_type j=0; j<this->n(); j++)
121  os << std::setw(8) << (*this)(i,j).imag() << " ";
122  os << std::endl;
123  }
124 }
125 
126 
127 
128 
129 
130 
131 // Full specialization for Real datatypes
132 template <typename T>
135  const SolverPackage solver_package)
136 {
137  // Build the appropriate vector
138  switch (solver_package)
139  {
140 
141 
142 #ifdef LIBMESH_HAVE_LASPACK
143  case LASPACK_SOLVERS:
144  {
146  return ap;
147  }
148 #endif
149 
150 
151 #ifdef LIBMESH_HAVE_PETSC
152  case PETSC_SOLVERS:
153  {
154  AutoPtr<SparseMatrix<T> > ap(new PetscMatrix<T>(comm));
155  return ap;
156  }
157 #endif
158 
159 
160 #ifdef LIBMESH_HAVE_TRILINOS
161  case TRILINOS_SOLVERS:
162  {
164  return ap;
165  }
166 #endif
167 
168 
169 #ifdef LIBMESH_HAVE_EIGEN
170  case EIGEN_SOLVERS:
171  {
173  return ap;
174  }
175 #endif
176 
177 
178  default:
179  libMesh::err << "ERROR: Unrecognized solver package: "
180  << solver_package
181  << std::endl;
182  libmesh_error();
183  }
184 
185  AutoPtr<SparseMatrix<T> > ap(NULL);
186  return ap;
187 }
188 
189 
190 template <typename T>
192  const NumericVector<T>& arg) const
193 {
194  dest.zero();
195  this->vector_mult_add(dest,arg);
196 }
197 
198 
199 
200 template <typename T>
202  const NumericVector<T>& arg) const
203 {
204  /* This functionality is actually implemented in the \p
205  NumericVector class. */
206  dest.add_vector(arg,*this);
207 }
208 
209 
210 
211 template <typename T>
212 void SparseMatrix<T>::zero_rows (std::vector<numeric_index_type> &, T)
213 {
214  /* This functionality isn't implemented or stubbed in every subclass yet */
215  libmesh_not_implemented();
216 }
217 
218 
219 
220 template <typename T>
221 void SparseMatrix<T>::print(std::ostream& os, const bool sparse) const
222 {
223  parallel_object_only();
224 
225  libmesh_assert (this->initialized());
226 
227  if(!this->_dof_map)
228  {
229  os << std::endl << "Error! Trying to print a matrix with no dof_map set!" << std::endl << std::endl;
230  libmesh_error();
231  }
232 
233  // We'll print the matrix from processor 0 to make sure
234  // it's serialized properly
235  if (this->processor_id() == 0)
236  {
237  libmesh_assert_equal_to (this->_dof_map->first_dof(), 0);
238  for (numeric_index_type i=this->_dof_map->first_dof();
239  i!=this->_dof_map->end_dof(); ++i)
240  {
241  if(sparse)
242  {
243  for (numeric_index_type j=0; j<this->n(); j++)
244  {
245  T c = (*this)(i,j);
246  if (c != static_cast<T>(0.0))
247  {
248  os << i << " " << j << " " << c << std::endl;
249  }
250  }
251  }
252  else
253  {
254  for (numeric_index_type j=0; j<this->n(); j++)
255  os << (*this)(i,j) << " ";
256  os << std::endl;
257  }
258  }
259 
260  std::vector<numeric_index_type> ibuf, jbuf;
261  std::vector<T> cbuf;
262  numeric_index_type currenti = this->_dof_map->end_dof();
263  for (processor_id_type p=1; p < this->n_processors(); ++p)
264  {
265  this->comm().receive(p, ibuf);
266  this->comm().receive(p, jbuf);
267  this->comm().receive(p, cbuf);
268  libmesh_assert_equal_to (ibuf.size(), jbuf.size());
269  libmesh_assert_equal_to (ibuf.size(), cbuf.size());
270 
271  if (ibuf.empty())
272  continue;
273  libmesh_assert_greater_equal (ibuf.front(), currenti);
274  libmesh_assert_greater_equal (ibuf.back(), ibuf.front());
275 
276  std::size_t currentb = 0;
277  for (;currenti <= ibuf.back(); ++currenti)
278  {
279  if(sparse)
280  {
281  for (numeric_index_type j=0; j<this->n(); j++)
282  {
283  if (currentb < ibuf.size() &&
284  ibuf[currentb] == currenti &&
285  jbuf[currentb] == j)
286  {
287  os << currenti << " " << j << " " << cbuf[currentb] << std::endl;
288  currentb++;
289  }
290  }
291  }
292  else
293  {
294  for (numeric_index_type j=0; j<this->n(); j++)
295  {
296  if (currentb < ibuf.size() &&
297  ibuf[currentb] == currenti &&
298  jbuf[currentb] == j)
299  {
300  os << cbuf[currentb] << " ";
301  currentb++;
302  }
303  else
304  os << static_cast<T>(0.0) << " ";
305  }
306  os << std::endl;
307  }
308  }
309  }
310  if(!sparse)
311  {
312  for (; currenti != this->m(); ++currenti)
313  {
314  for (numeric_index_type j=0; j<this->n(); j++)
315  os << static_cast<T>(0.0) << " ";
316  os << std::endl;
317  }
318  }
319  }
320  else
321  {
322  std::vector<numeric_index_type> ibuf, jbuf;
323  std::vector<T> cbuf;
324 
325  // We'll assume each processor has access to entire
326  // matrix rows, so (*this)(i,j) is valid if i is a local index.
327  for (numeric_index_type i=this->_dof_map->first_dof();
328  i!=this->_dof_map->end_dof(); ++i)
329  {
330  for (numeric_index_type j=0; j<this->n(); j++)
331  {
332  T c = (*this)(i,j);
333  if (c != static_cast<T>(0.0))
334  {
335  ibuf.push_back(i);
336  jbuf.push_back(j);
337  cbuf.push_back(c);
338  }
339  }
340  }
341  this->comm().send(0,ibuf);
342  this->comm().send(0,jbuf);
343  this->comm().send(0,cbuf);
344  }
345 }
346 
347 
348 
349 //------------------------------------------------------------------
350 // Explicit instantiations
351 template class SparseMatrix<Number>;
352 
353 } // namespace libMesh

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:07 UTC

Hosted By:
SourceForge.net Logo