laspack_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/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"
31 
32 namespace libMesh
33 {
34 
35 
36 //-----------------------------------------------------------------------
37 // LaspackMatrix members
38 template <typename T>
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  // Initize 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>
127 void LaspackMatrix<T>::init (const numeric_index_type libmesh_dbg_var(m_in),
128  const numeric_index_type libmesh_dbg_var(n_in),
129  const numeric_index_type libmesh_dbg_var(m_l),
130  const numeric_index_type libmesh_dbg_var(n_l),
131  const numeric_index_type libmesh_dbg_var(nnz),
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 
142  libMesh::err << "ERROR: Only the init() member that uses the" << std::endl
143  << "DofMap is implemented for Laspack matrices!" << std::endl;
144  libmesh_error();
145 
146  this->_is_initialized = true;
147 }
148 
149 
150 
151 template <typename T>
153 {
154  // Ignore calls on initialized objects
155  if (this->initialized())
156  return;
157 
158  // We need the DofMap for this!
159  libmesh_assert(this->_dof_map);
160 
161  // Clear intialized matrices
162  if (this->initialized())
163  this->clear();
164 
165  const numeric_index_type n_rows = this->_dof_map->n_dofs();
166 #ifndef NDEBUG
167  // The following variables are only used for assertions,
168  // so avoid declaring them when asserts are inactive.
169  const numeric_index_type n_cols = n_rows;
170  const numeric_index_type n_l = this->_dof_map->n_dofs_on_processor(0);
171  const numeric_index_type m_l = n_l;
172 #endif
173 
174  // Laspack Matrices only work for uniprocessor cases
175  libmesh_assert_equal_to (n_rows, n_cols);
176  libmesh_assert_equal_to (m_l, n_rows);
177  libmesh_assert_equal_to (n_l, n_cols);
178 
179 #ifndef NDEBUG
180  // The following variables are only used for assertions,
181  // so avoid declaring them when asserts are inactive.
182  const std::vector<numeric_index_type>& n_nz = this->_dof_map->get_n_nz();
183  const std::vector<numeric_index_type>& n_oz = this->_dof_map->get_n_oz();
184 #endif
185 
186  // Make sure the sparsity pattern isn't empty
187  libmesh_assert_equal_to (n_nz.size(), n_l);
188  libmesh_assert_equal_to (n_oz.size(), n_l);
189 
190  if (n_rows==0)
191  return;
192 
193  Q_Constr(&_QMat, const_cast<char*>("Mat"), n_rows, _LPFalse, Rowws, Normal, _LPTrue);
194 
195  this->_is_initialized = true;
196 
197  libmesh_assert_equal_to (n_rows, this->m());
198 }
199 
200 
201 
202 template <typename T>
204  const std::vector<numeric_index_type>& rows,
205  const std::vector<numeric_index_type>& cols)
206 
207 {
208  libmesh_assert (this->initialized());
209  unsigned int n_rows = libmesh_cast_int<unsigned int>(rows.size());
210  unsigned int n_cols = libmesh_cast_int<unsigned int>(cols.size());
211  libmesh_assert_equal_to (dm.m(), n_rows);
212  libmesh_assert_equal_to (dm.n(), n_cols);
213 
214 
215  for (unsigned int i=0; i<n_rows; i++)
216  for (unsigned int j=0; j<n_cols; j++)
217  this->add(rows[i],cols[j],dm(i,j));
218 }
219 
220 
221 
222 template <typename T>
224 {
225  libmesh_not_implemented();
226 }
227 
228 
229 
230 template <typename T>
232 {
233  libmesh_not_implemented();
234 }
235 
236 
237 
238 template <typename T>
240  SparseMatrix<T>(comm),
241  _closed (false)
242 {
243 }
244 
245 
246 
247 template <typename T>
249 {
250  this->clear ();
251 }
252 
253 
254 
255 template <typename T>
257 {
258  if (this->initialized())
259  {
260  Q_Destr(&_QMat);
261  }
262 
263  _csr.clear();
264  _row_start.clear();
265  _closed = false;
266  this->_is_initialized = false;
267 }
268 
269 
270 
271 template <typename T>
273 {
274  const numeric_index_type n_rows = this->m();
275 
276  for (numeric_index_type row=0; row<n_rows; row++)
277  {
278  const std::vector<numeric_index_type>::const_iterator
279  r_start = _row_start[row];
280 
281  const numeric_index_type len = (_row_start[row+1] - _row_start[row]);
282 
283  // Make sure we agree on the row length
284  libmesh_assert_equal_to (len, Q_GetLen(&_QMat, row+1));
285 
286  for (numeric_index_type l=0; l<len; l++)
287  {
288  const numeric_index_type j = *(r_start + l);
289 
290  // Make sure the data structures are working
291  libmesh_assert_equal_to ((j+1), Q_GetPos (&_QMat, row+1, l));
292 
293  Q_SetEntry (&_QMat, row+1, l, j+1, 0.);
294  }
295  }
296 }
297 
298 
299 
300 template <typename T>
302 {
303  libmesh_assert (this->initialized());
304 
305  return static_cast<numeric_index_type>(Q_GetDim(const_cast<QMatrix*>(&_QMat)));
306 }
307 
308 
309 
310 template <typename T>
312 {
313  libmesh_assert (this->initialized());
314 
315  return static_cast<numeric_index_type>(Q_GetDim(const_cast<QMatrix*>(&_QMat)));
316 }
317 
318 
319 
320 template <typename T>
322 {
323  return 0;
324 }
325 
326 
327 
328 template <typename T>
330 {
331  return this->m();
332 }
333 
334 
335 
336 template <typename T>
338  const numeric_index_type j,
339  const T value)
340 {
341  libmesh_assert (this->initialized());
342  libmesh_assert_less (i, this->m());
343  libmesh_assert_less (j, this->n());
344 
345  const numeric_index_type position = this->pos(i,j);
346 
347  // Sanity check
348  libmesh_assert_equal_to (*(_row_start[i]+position), j);
349  libmesh_assert_equal_to ((j+1), Q_GetPos (&_QMat, i+1, position));
350 
351  Q_SetEntry (&_QMat, i+1, position, j+1, value);
352 }
353 
354 
355 
356 template <typename T>
358  const numeric_index_type j,
359  const T value)
360 {
361  libmesh_assert (this->initialized());
362  libmesh_assert_less (i, this->m());
363  libmesh_assert_less (j, this->n());
364 
365  const numeric_index_type position = this->pos(i,j);
366 
367  // Sanity check
368  libmesh_assert_equal_to (*(_row_start[i]+position), j);
369 
370  Q_AddVal (&_QMat, i+1, position, value);
371 }
372 
373 
374 
375 template <typename T>
377  const std::vector<numeric_index_type>& dof_indices)
378 {
379  this->add_matrix (dm, dof_indices, dof_indices);
380 }
381 
382 
383 
384 template <typename T>
385 void LaspackMatrix<T>::add (const T a_in, SparseMatrix<T> &X_in)
386 {
387  libmesh_assert (this->initialized());
388  libmesh_assert_equal_to (this->m(), X_in.m());
389  libmesh_assert_equal_to (this->n(), X_in.n());
390 
391  LaspackMatrix<T>* X = libmesh_cast_ptr<LaspackMatrix<T>*> (&X_in);
392  _LPNumber a = static_cast<_LPNumber> (a_in);
393 
394  libmesh_assert(X);
395 
396  // loops taken from LaspackMatrix<T>::zero ()
397 
398  const numeric_index_type n_rows = this->m();
399 
400  for (numeric_index_type row=0; row<n_rows; row++)
401  {
402  const std::vector<numeric_index_type>::const_iterator
403  r_start = _row_start[row];
404 
405  const numeric_index_type len = (_row_start[row+1] - _row_start[row]);
406 
407  // Make sure we agree on the row length
408  libmesh_assert_equal_to (len, Q_GetLen(&_QMat, row+1));
409  // compare matrix sparsity structures
410  libmesh_assert_equal_to (len, Q_GetLen(&(X->_QMat), row+1));
411 
412 
413  for (numeric_index_type l=0; l<len; l++)
414  {
415  const numeric_index_type j = *(r_start + l);
416 
417  // Make sure the data structures are working
418  libmesh_assert_equal_to ((j+1), Q_GetPos (&_QMat, row+1, l));
419 
420  const _LPNumber value = a * Q_GetEl(const_cast<QMatrix*>(&(X->_QMat)), row+1, j+1);
421  Q_AddVal (&_QMat, row+1, l, value);
422  }
423  }
424 }
425 
426 
427 
428 
429 template <typename T>
431  const numeric_index_type j) const
432 {
433  libmesh_assert (this->initialized());
434  libmesh_assert_less (i, this->m());
435  libmesh_assert_less (j, this->n());
436 
437  return Q_GetEl (const_cast<QMatrix*>(&_QMat), i+1, j+1);
438 }
439 
440 
441 
442 template <typename T>
444  const numeric_index_type j) const
445 {
446  libmesh_assert_less (i, this->m());
447  libmesh_assert_less (j, this->n());
448  libmesh_assert_less (i+1, _row_start.size());
449  libmesh_assert (_row_start.back() == _csr.end());
450 
451  // note this requires the _csr to be
452  std::pair<std::vector<numeric_index_type>::const_iterator,
453  std::vector<numeric_index_type>::const_iterator> p =
454  std::equal_range (_row_start[i],
455  _row_start[i+1],
456  j);
457 
458  // Make sure the row contains the element j
459  libmesh_assert (p.first != p.second);
460 
461  // Make sure the values match
462  libmesh_assert (*p.first == j);
463 
464  // Return the position in the compressed row
465  return std::distance (_row_start[i], p.first);
466 }
467 
468 
469 
470 //------------------------------------------------------------------
471 // Explicit instantiations
472 template class LaspackMatrix<Number>;
473 
474 } // namespace libMesh
475 
476 
477 #endif // #ifdef LIBMESH_HAVE_LASPACK

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

Hosted By:
SourceForge.net Logo