libMesh
laspack_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/libmesh_config.h"
24 
25 #ifdef LIBMESH_HAVE_LASPACK
26 
27 #include "libmesh/laspack_matrix.h"
28 #include "libmesh/dense_matrix.h"
29 #include "libmesh/dof_map.h"
30 #include "libmesh/sparsity_pattern.h"
31 
32 namespace libMesh
33 {
34 
35 
36 //-----------------------------------------------------------------------
37 // LaspackMatrix members
38 template <typename T>
39 void LaspackMatrix<T>::update_sparsity_pattern (const SparsityPattern::Graph & sparsity_pattern)
40 {
41  // clear data, start over
42  this->clear ();
43 
44  // big trouble if this fails!
45  libmesh_assert(this->_dof_map);
46 
47  const numeric_index_type n_rows = sparsity_pattern.size();
48 
49  // Initialize the _row_start data structure,
50  // allocate storage for the _csr array
51  {
52  std::size_t size = 0;
53 
54  for (numeric_index_type row=0; row<n_rows; row++)
55  size += sparsity_pattern[row].size();
56 
57  _csr.resize (size);
58  _row_start.reserve(n_rows + 1);
59  }
60 
61 
62  // Initialize the _csr data structure.
63  {
64  std::vector<numeric_index_type>::iterator position = _csr.begin();
65 
66  _row_start.push_back (position);
67 
68  for (numeric_index_type row=0; row<n_rows; row++)
69  {
70  // insert the row indices
71  for (SparsityPattern::Row::const_iterator col = sparsity_pattern[row].begin();
72  col != sparsity_pattern[row].end(); ++col)
73  {
74  libmesh_assert (position != _csr.end());
75  *position = *col;
76  ++position;
77  }
78 
79  _row_start.push_back (position);
80  }
81  }
82 
83 
84  // Initialize the matrix
85  libmesh_assert (!this->initialized());
86  this->init ();
87  libmesh_assert (this->initialized());
88  //libMesh::out << "n_rows=" << n_rows << std::endl;
89  //libMesh::out << "m()=" << m() << std::endl;
90  libmesh_assert_equal_to (n_rows, this->m());
91 
92  // Tell the matrix about its structure. Initialize it
93  // to zero.
94  for (numeric_index_type i=0; i<n_rows; i++)
95  {
96  const std::vector<numeric_index_type>::const_iterator
97  rs = _row_start[i];
98 
99  const numeric_index_type length = _row_start[i+1] - rs;
100 
101  Q_SetLen (&_QMat, i+1, length);
102 
103  for (numeric_index_type l=0; l<length; l++)
104  {
105  const numeric_index_type j = *(rs+l);
106 
107  // sanity check
108  //libMesh::out << "m()=" << m() << std::endl;
109  //libMesh::out << "(i,j,l) = (" << i
110  // << "," << j
111  // << "," << l
112  // << ")" << std::endl;
113  //libMesh::out << "pos(i,j)=" << pos(i,j)
114  // << std::endl;
115  libmesh_assert_equal_to (this->pos(i,j), l);
116  Q_SetEntry (&_QMat, i+1, l, j+1, 0.);
117  }
118  }
119 
120  // That's it!
121  //libmesh_here();
122 }
123 
124 
125 
126 template <typename T>
132  const numeric_index_type,
133  const numeric_index_type)
134 {
135  // noz ignored... only used for multiple processors!
136  libmesh_assert_equal_to (m_in, m_l);
137  libmesh_assert_equal_to (n_in, n_l);
138  libmesh_assert_equal_to (m_in, n_in);
139  libmesh_assert_greater (nnz, 0);
140 
141  libmesh_error_msg("ERROR: Only the init() member that uses the DofMap is implemented for Laspack matrices!");
142 
143  this->_is_initialized = true;
144 }
145 
146 
147 
148 template <typename T>
150 {
151  // Ignore calls on initialized objects
152  if (this->initialized())
153  return;
154 
155  // We need the DofMap for this!
156  libmesh_assert(this->_dof_map);
157 
158  // Clear initialized matrices
159  if (this->initialized())
160  this->clear();
161 
162  const numeric_index_type n_rows = this->_dof_map->n_dofs();
163 #ifndef NDEBUG
164  // The following variables are only used for assertions,
165  // so avoid declaring them when asserts are inactive.
166  const numeric_index_type n_cols = n_rows;
167  const numeric_index_type n_l = this->_dof_map->n_dofs_on_processor(0);
168  const numeric_index_type m_l = n_l;
169 #endif
170 
171  // Laspack Matrices only work for uniprocessor cases
172  libmesh_assert_equal_to (n_rows, n_cols);
173  libmesh_assert_equal_to (m_l, n_rows);
174  libmesh_assert_equal_to (n_l, n_cols);
175 
176 #ifndef NDEBUG
177  // The following variables are only used for assertions,
178  // so avoid declaring them when asserts are inactive.
179  const std::vector<numeric_index_type> & n_nz = this->_dof_map->get_n_nz();
180  const std::vector<numeric_index_type> & n_oz = this->_dof_map->get_n_oz();
181 #endif
182 
183  // Make sure the sparsity pattern isn't empty
184  libmesh_assert_equal_to (n_nz.size(), n_l);
185  libmesh_assert_equal_to (n_oz.size(), n_l);
186 
187  if (n_rows==0)
188  return;
189 
190  Q_Constr(&_QMat, const_cast<char *>("Mat"), n_rows, _LPFalse, Rowws, Normal, _LPTrue);
191 
192  this->_is_initialized = true;
193 
194  libmesh_assert_equal_to (n_rows, this->m());
195 }
196 
197 
198 
199 template <typename T>
200 void LaspackMatrix<T>::add_matrix(const DenseMatrix<T> & dm,
201  const std::vector<numeric_index_type> & rows,
202  const std::vector<numeric_index_type> & cols)
203 
204 {
205  libmesh_assert (this->initialized());
206  unsigned int n_rows = cast_int<unsigned int>(rows.size());
207  unsigned int n_cols = cast_int<unsigned int>(cols.size());
208  libmesh_assert_equal_to (dm.m(), n_rows);
209  libmesh_assert_equal_to (dm.n(), n_cols);
210 
211 
212  for (unsigned int i=0; i<n_rows; i++)
213  for (unsigned int j=0; j<n_cols; j++)
214  this->add(rows[i],cols[j],dm(i,j));
215 }
216 
217 
218 
219 template <typename T>
220 void LaspackMatrix<T>::get_diagonal (NumericVector<T> & /*dest*/) const
221 {
222  libmesh_not_implemented();
223 }
224 
225 
226 
227 template <typename T>
228 void LaspackMatrix<T>::get_transpose (SparseMatrix<T> & /*dest*/) const
229 {
230  libmesh_not_implemented();
231 }
232 
233 
234 
235 template <typename T>
236 LaspackMatrix<T>::LaspackMatrix (const Parallel::Communicator & comm) :
237  SparseMatrix<T>(comm),
238  _closed (false)
239 {
240 }
241 
242 
243 
244 template <typename T>
245 LaspackMatrix<T>::~LaspackMatrix ()
246 {
247  this->clear ();
248 }
249 
250 
251 
252 template <typename T>
253 void LaspackMatrix<T>::clear ()
254 {
255  if (this->initialized())
256  {
257  Q_Destr(&_QMat);
258  }
259 
260  _csr.clear();
261  _row_start.clear();
262  _closed = false;
263  this->_is_initialized = false;
264 }
265 
266 
267 
268 template <typename T>
270 {
271  const numeric_index_type n_rows = this->m();
272 
273  for (numeric_index_type row=0; row<n_rows; row++)
274  {
275  const std::vector<numeric_index_type>::const_iterator
276  r_start = _row_start[row];
277 
278  const numeric_index_type len = (_row_start[row+1] - _row_start[row]);
279 
280  // Make sure we agree on the row length
281  libmesh_assert_equal_to (len, Q_GetLen(&_QMat, row+1));
282 
283  for (numeric_index_type l=0; l<len; l++)
284  {
285  const numeric_index_type j = *(r_start + l);
286 
287  // Make sure the data structures are working
288  libmesh_assert_equal_to ((j+1), Q_GetPos (&_QMat, row+1, l));
289 
290  Q_SetEntry (&_QMat, row+1, l, j+1, 0.);
291  }
292  }
293 }
294 
295 
296 
297 template <typename T>
298 numeric_index_type LaspackMatrix<T>::m () const
299 {
300  libmesh_assert (this->initialized());
301 
302  return static_cast<numeric_index_type>(Q_GetDim(const_cast<QMatrix*>(&_QMat)));
303 }
304 
305 
306 
307 template <typename T>
308 numeric_index_type LaspackMatrix<T>::n () const
309 {
310  libmesh_assert (this->initialized());
311 
312  return static_cast<numeric_index_type>(Q_GetDim(const_cast<QMatrix*>(&_QMat)));
313 }
314 
315 
316 
317 template <typename T>
318 numeric_index_type LaspackMatrix<T>::row_start () const
319 {
320  return 0;
321 }
322 
323 
324 
325 template <typename T>
326 numeric_index_type LaspackMatrix<T>::row_stop () const
327 {
328  return this->m();
329 }
330 
331 
332 
333 template <typename T>
334 void LaspackMatrix<T>::set (const numeric_index_type i,
335  const numeric_index_type j,
336  const T value)
337 {
338  libmesh_assert (this->initialized());
339  libmesh_assert_less (i, this->m());
340  libmesh_assert_less (j, this->n());
341 
342  const numeric_index_type position = this->pos(i,j);
343 
344  // Sanity check
345  libmesh_assert_equal_to (*(_row_start[i]+position), j);
346  libmesh_assert_equal_to ((j+1), Q_GetPos (&_QMat, i+1, position));
347 
348  Q_SetEntry (&_QMat, i+1, position, j+1, value);
349 }
350 
351 
352 
353 template <typename T>
354 void LaspackMatrix<T>::add (const numeric_index_type i,
355  const numeric_index_type j,
356  const T value)
357 {
358  libmesh_assert (this->initialized());
359  libmesh_assert_less (i, this->m());
360  libmesh_assert_less (j, this->n());
361 
362  const numeric_index_type position = this->pos(i,j);
363 
364  // Sanity check
365  libmesh_assert_equal_to (*(_row_start[i]+position), j);
366 
367  Q_AddVal (&_QMat, i+1, position, value);
368 }
369 
370 
371 
372 template <typename T>
373 void LaspackMatrix<T>::add_matrix(const DenseMatrix<T> & dm,
374  const std::vector<numeric_index_type> & dof_indices)
375 {
376  this->add_matrix (dm, dof_indices, dof_indices);
377 }
378 
379 
380 
381 template <typename T>
382 void LaspackMatrix<T>::add (const T a_in, SparseMatrix<T> & X_in)
383 {
384  libmesh_assert (this->initialized());
385  libmesh_assert_equal_to (this->m(), X_in.m());
386  libmesh_assert_equal_to (this->n(), X_in.n());
387 
388  LaspackMatrix<T> * X = cast_ptr<LaspackMatrix<T> *> (&X_in);
389  _LPNumber a = static_cast<_LPNumber> (a_in);
390 
391  libmesh_assert(X);
392 
393  // loops taken from LaspackMatrix<T>::zero ()
394 
395  const numeric_index_type n_rows = this->m();
396 
397  for (numeric_index_type row=0; row<n_rows; row++)
398  {
399  const std::vector<numeric_index_type>::const_iterator
400  r_start = _row_start[row];
401 
402  const numeric_index_type len = (_row_start[row+1] - _row_start[row]);
403 
404  // Make sure we agree on the row length
405  libmesh_assert_equal_to (len, Q_GetLen(&_QMat, row+1));
406  // compare matrix sparsity structures
407  libmesh_assert_equal_to (len, Q_GetLen(&(X->_QMat), row+1));
408 
409 
410  for (numeric_index_type l=0; l<len; l++)
411  {
412  const numeric_index_type j = *(r_start + l);
413 
414  // Make sure the data structures are working
415  libmesh_assert_equal_to ((j+1), Q_GetPos (&_QMat, row+1, l));
416 
417  const _LPNumber value = a * Q_GetEl(const_cast<QMatrix*>(&(X->_QMat)), row+1, j+1);
418  Q_AddVal (&_QMat, row+1, l, value);
419  }
420  }
421 }
422 
423 
424 
425 
426 template <typename T>
427 T LaspackMatrix<T>::operator () (const numeric_index_type i,
428  const numeric_index_type j) const
429 {
430  libmesh_assert (this->initialized());
431  libmesh_assert_less (i, this->m());
432  libmesh_assert_less (j, this->n());
433 
434  return Q_GetEl (const_cast<QMatrix*>(&_QMat), i+1, j+1);
435 }
436 
437 
438 
439 template <typename T>
440 numeric_index_type LaspackMatrix<T>::pos (const numeric_index_type i,
441  const numeric_index_type j) const
442 {
443  libmesh_assert_less (i, this->m());
444  libmesh_assert_less (j, this->n());
445  libmesh_assert_less (i+1, _row_start.size());
446  libmesh_assert (_row_start.back() == _csr.end());
447 
448  // note this requires the _csr to be
449  std::pair<std::vector<numeric_index_type>::const_iterator,
450  std::vector<numeric_index_type>::const_iterator> p =
451  std::equal_range (_row_start[i],
452  _row_start[i+1],
453  j);
454 
455  // Make sure the row contains the element j
456  libmesh_assert (p.first != p.second);
457 
458  // Make sure the values match
459  libmesh_assert (*p.first == j);
460 
461  // Return the position in the compressed row
462  return std::distance (_row_start[i], p.first);
463 }
464 
465 
466 
467 //------------------------------------------------------------------
468 // Explicit instantiations
469 template class LaspackMatrix<Number>;
470 
471 } // namespace libMesh
472 
473 
474 #endif // #ifdef LIBMESH_HAVE_LASPACK
The libMesh namespace provides an interface to certain functionality in the library.
const Number zero
.
Definition: libmesh.h:178
Real distance(const Point &p)
libmesh_assert(j)
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
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
PetscErrorCode Vec Mat libmesh_dbg_var(j)
bool _is_initialized
Flag indicating if the data structures have been initialized.
static const bool value
Definition: xdr_io.C:108
bool initialized()
Checks that library initialization has been done.
Definition: libmesh.C:274
bool initialized() const
virtual void clear() libmesh_override
Release all memory and clear data structures.