libMesh
laspack_matrix.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2025 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 // Local includes
21 #include "libmesh/libmesh_config.h"
22 
23 #ifdef LIBMESH_HAVE_LASPACK
24 
25 #include "libmesh/laspack_matrix.h"
26 #include "libmesh/dense_matrix.h"
27 #include "libmesh/dof_map.h"
28 #include "libmesh/sparsity_pattern.h"
29 
30 
31 // C++ Includes
32 #include <memory>
33 
34 
35 namespace libMesh
36 {
37 
38 
39 //-----------------------------------------------------------------------
40 // LaspackMatrix members
41 template <typename T>
43 {
44  // clear data, start over
45  this->clear ();
46 
47  // big trouble if this fails!
48  libmesh_assert(this->_dof_map);
49 
50  const numeric_index_type n_rows = sparsity_pattern.size();
51 
52  // Initialize the _row_start data structure,
53  // allocate storage for the _csr array
54  {
55  std::size_t size = 0;
56 
57  for (numeric_index_type row=0; row<n_rows; row++)
58  size += sparsity_pattern[row].size();
59 
60  _csr.resize (size);
61  _row_start.reserve(n_rows + 1);
62  }
63 
64 
65  // Initialize the _csr data structure.
66  {
67  std::vector<numeric_index_type>::iterator position = _csr.begin();
68 
69  _row_start.push_back (position);
70 
71  for (numeric_index_type row=0; row<n_rows; row++)
72  {
73  // insert the row indices
74  for (const auto & col : sparsity_pattern[row])
75  {
76  libmesh_assert (position != _csr.end());
77  *position = col;
78  ++position;
79  }
80 
81  _row_start.push_back (position);
82  }
83  }
84 
85 
86  // Initialize the matrix
87  libmesh_assert (!this->initialized());
88  this->init ();
89  libmesh_assert (this->initialized());
90  //libMesh::out << "n_rows=" << n_rows << std::endl;
91  //libMesh::out << "m()=" << m() << std::endl;
92  libmesh_assert_equal_to (n_rows, this->m());
93 
94  // Tell the matrix about its structure. Initialize it
95  // to zero.
96  for (numeric_index_type i=0; i<n_rows; i++)
97  {
98  auto rs = _row_start[i];
99 
100  const numeric_index_type length = _row_start[i+1] - rs;
101 
102  Q_SetLen (&_QMat, i+1, length);
103 
104  for (numeric_index_type l=0; l<length; l++)
105  {
106  const numeric_index_type j = *(rs+l);
107 
108  // sanity check
109  //libMesh::out << "m()=" << m() << std::endl;
110  //libMesh::out << "(i,j,l) = (" << i
111  // << "," << j
112  // << "," << l
113  // << ")" << std::endl;
114  //libMesh::out << "pos(i,j)=" << pos(i,j)
115  // << std::endl;
116  libmesh_assert_equal_to (this->pos(i,j), l);
117  Q_SetEntry (&_QMat, i+1, l, j+1, 0.);
118  }
119  }
120 
121  // That's it!
122  //libmesh_here();
123 }
124 
125 
126 
127 template <typename T>
129  const numeric_index_type n_in,
130  const numeric_index_type m_l,
131  const numeric_index_type n_l,
132  const numeric_index_type nnz,
133  const numeric_index_type,
134  const numeric_index_type)
135 {
136  // noz ignored... only used for multiple processors!
137  libmesh_assert_equal_to (m_in, m_l);
138  libmesh_assert_equal_to (n_in, n_l);
139  libmesh_assert_equal_to (m_in, n_in);
140  libmesh_assert_greater (nnz, 0);
141 
142  if ((m_in != m_l) ||
143  (n_in != n_l))
144  libmesh_not_implemented_msg("Laspack does not support distributed matrices");
145 
146  if (m_in != n_in)
147  libmesh_not_implemented_msg("Laspack does not support rectangular matrices");
148 
149  if (nnz < n_in)
150  libmesh_warning("Using inefficient LaspackMatrix allocation via this init() method");
151 
152  const dof_id_type n_rows = m_in;
153 
154  // Laspack doesn't let us assign sparse indices on the fly, so
155  // without a sparsity pattern we're stuck allocating a dense matrix
156  {
157  _csr.resize (m_in * m_in);
158  _row_start.reserve(m_in + 1);
159  }
160 
161  // Initialize the _csr data structure.
162  {
163  std::vector<numeric_index_type>::iterator position = _csr.begin();
164 
165  _row_start.push_back (position);
166 
167  for (numeric_index_type row=0; row<n_rows; row++)
168  {
169  // insert the row indices
170  for (const auto & col : make_range(n_rows))
171  {
172  libmesh_assert (position != _csr.end());
173  *position = col;
174  ++position;
175  }
176 
177  _row_start.push_back (position);
178  }
179  }
180 
181  Q_Constr(&_QMat, const_cast<char *>("Mat"), m_in, _LPFalse, Rowws, Normal, _LPTrue);
182 
183  this->_is_initialized = true;
184 
185  // Tell the matrix about its structure. Initialize it
186  // to zero.
187  for (numeric_index_type i=0; i<n_rows; i++)
188  {
189  auto rs = _row_start[i];
190 
191  const numeric_index_type length = _row_start[i+1] - rs;
192 
193  Q_SetLen (&_QMat, i+1, length);
194 
195  for (numeric_index_type l=0; l<length; l++)
196  {
197  const numeric_index_type j = *(rs+l);
198 
199  libmesh_assert_equal_to (this->pos(i,j), l);
200  Q_SetEntry (&_QMat, i+1, l, j+1, 0.);
201  }
202  }
203 }
204 
205 
206 
207 template <typename T>
209 {
210  // Ignore calls on initialized objects
211  if (this->initialized())
212  return;
213 
214  // We need a sparsity pattern for this!
215  libmesh_assert(this->_sp);
216 
217  // Clear initialized matrices
218  if (this->initialized())
219  this->clear();
220 
221  const numeric_index_type n_rows = this->_dof_map->n_dofs();
222 #ifndef NDEBUG
223  // The following variables are only used for assertions,
224  // so avoid declaring them when asserts are inactive.
225  const numeric_index_type n_cols = n_rows;
226  const numeric_index_type n_l = this->_dof_map->n_dofs_on_processor(0);
227  const numeric_index_type m_l = n_l;
228 #endif
229 
230  // Laspack Matrices only work for uniprocessor cases
231  libmesh_assert_equal_to (n_rows, n_cols);
232  libmesh_assert_equal_to (m_l, n_rows);
233  libmesh_assert_equal_to (n_l, n_cols);
234 
235 #ifndef NDEBUG
236  // The following variables are only used for assertions,
237  // so avoid declaring them when asserts are inactive.
238  const std::vector<numeric_index_type> & n_nz = this->_sp->get_n_nz();
239  const std::vector<numeric_index_type> & n_oz = this->_sp->get_n_oz();
240 #endif
241 
242  // Make sure the sparsity pattern isn't empty
243  libmesh_assert_equal_to (n_nz.size(), n_l);
244  libmesh_assert_equal_to (n_oz.size(), n_l);
245 
246  if (n_rows==0)
247  return;
248 
249  Q_Constr(&_QMat, const_cast<char *>("Mat"), n_rows, _LPFalse, Rowws, Normal, _LPTrue);
250 
251  this->_is_initialized = true;
252 
253  libmesh_assert_equal_to (n_rows, this->m());
254 }
255 
256 
257 
258 template <typename T>
260  const std::vector<numeric_index_type> & rows,
261  const std::vector<numeric_index_type> & cols)
262 
263 {
264  libmesh_assert (this->initialized());
265  unsigned int n_rows = cast_int<unsigned int>(rows.size());
266  unsigned int n_cols = cast_int<unsigned int>(cols.size());
267  libmesh_assert_equal_to (dm.m(), n_rows);
268  libmesh_assert_equal_to (dm.n(), n_cols);
269 
270 
271  for (unsigned int i=0; i<n_rows; i++)
272  for (unsigned int j=0; j<n_cols; j++)
273  this->add(rows[i],cols[j],dm(i,j));
274 }
275 
276 
277 
278 template <typename T>
280 {
281  // There does not seem to be a straightforward way to iterate over
282  // the columns of a LaspackMatrix. So we use some extra storage and
283  // keep track of the column sums while going over the row entries...
284  std::vector<Real> abs_col_sums(this->n());
285 
286  const numeric_index_type n_rows = this->m();
287 
288  for (numeric_index_type row : make_range(n_rows))
289  {
290  auto r_start = _row_start[row];
291 
292  const numeric_index_type len = (_row_start[row+1] - _row_start[row]);
293 
294  // Make sure we agree on the row length
295  libmesh_assert_equal_to (len, Q_GetLen(&_QMat, row+1));
296 
297  for (numeric_index_type l=0; l<len; l++)
298  {
299  const numeric_index_type j = *(r_start + l);
300 
301  // Make sure the data structures are working
302  libmesh_assert_equal_to ((j+1), Q_GetPos (&_QMat, row+1, l));
303 
304  abs_col_sums[j] += std::abs(Q_GetVal (&_QMat, row+1, l));
305  }
306  }
307 
308  return *(std::max_element(abs_col_sums.begin(), abs_col_sums.end()));
309 }
310 
311 
312 
313 template <typename T>
315 {
316  const numeric_index_type n_rows = this->m();
317 
318  Real max_row_sum = 0;
319 
320  for (numeric_index_type row : make_range(n_rows))
321  {
322  Real row_sum = 0;
323 
324  const numeric_index_type len = (_row_start[row+1] - _row_start[row]);
325 
326  // Make sure we agree on the row length
327  libmesh_assert_equal_to (len, Q_GetLen(&_QMat, row+1));
328 
329  for (numeric_index_type l=0; l<len; l++)
330  {
331  // Make sure the data structures are working
332  libmesh_assert_equal_to ((*(_row_start[row] + l)+1), Q_GetPos
333  (&_QMat, row+1, l));
334 
335  row_sum += std::abs(Q_GetVal (&_QMat, row+1, l));
336  }
337 
338  max_row_sum = std::max(max_row_sum, row_sum);
339  }
340 
341  return max_row_sum;
342 }
343 
344 
345 
346 template <typename T>
348 {
349  libmesh_not_implemented();
350 }
351 
352 
353 
354 template <typename T>
356 {
357  libmesh_not_implemented();
358 }
359 
360 
361 
362 template <typename T>
364  std::vector<numeric_index_type> & indices,
365  std::vector<T> & values) const
366 {
367  indices.clear();
368  values.clear();
369 
370  const numeric_index_type len = (_row_start[i+1] - _row_start[i]);
371 
372  // Make sure we agree on the row length
373  libmesh_assert_equal_to (len, Q_GetLen(&_QMat, i+1));
374 
375  auto r_start = _row_start[i];
376 
377  for (numeric_index_type l=0; l<len; l++)
378  {
379  const numeric_index_type j = *(r_start + l);
380 
381  // Make sure the data structures are working
382  libmesh_assert_equal_to ((j+1), Q_GetPos (&_QMat, i+1, l));
383 
384  indices.push_back(j);
385 
386  values.push_back(Q_GetVal (&_QMat, i+1, l));
387  }
388 }
389 
390 
391 
392 template <typename T>
394  SparseMatrix<T>(comm),
395  _closed (false)
396 {
397 }
398 
399 
400 
401 template <typename T>
403 {
404  this->clear ();
405 }
406 
407 
408 
409 template <typename T>
411 {
412  if (this->initialized())
413  {
414  Q_Destr(&_QMat);
415  }
416 
417  _csr.clear();
418  _row_start.clear();
419  _closed = false;
420  this->_is_initialized = false;
421 }
422 
423 
424 
425 template <typename T>
427 {
428  const numeric_index_type n_rows = this->m();
429 
430  for (numeric_index_type row=0; row<n_rows; row++)
431  {
432  auto r_start = _row_start[row];
433 
434  const numeric_index_type len = (_row_start[row+1] - _row_start[row]);
435 
436  // Make sure we agree on the row length
437  libmesh_assert_equal_to (len, Q_GetLen(&_QMat, row+1));
438 
439  for (numeric_index_type l=0; l<len; l++)
440  {
441  const numeric_index_type j = *(r_start + l);
442 
443  // Make sure the data structures are working
444  libmesh_assert_equal_to ((j+1), Q_GetPos (&_QMat, row+1, l));
445 
446  Q_SetEntry (&_QMat, row+1, l, j+1, 0.);
447  }
448  }
449 
450  this->close();
451 }
452 
453 
454 
455 template <typename T>
456 std::unique_ptr<SparseMatrix<T>> LaspackMatrix<T>::zero_clone () const
457 {
458  // Make empty copy with matching comm, initialize, zero, and return.
459  auto mat_copy = std::make_unique<LaspackMatrix<T>>(this->comm());
460  if (this->_dof_map)
461  mat_copy->attach_dof_map(*this->_dof_map);
462  else
463  mat_copy->init(this->m(), this->n(), this->local_m(),
464  this->local_n(), this->local_n());
465 
466  mat_copy->zero();
467 
468  return mat_copy;
469 }
470 
471 
472 
473 template <typename T>
474 std::unique_ptr<SparseMatrix<T>> LaspackMatrix<T>::clone () const
475 {
476  // We don't currently have a faster implementation than making a
477  // zero clone and then filling in the values.
478  auto mat_copy = this->zero_clone();
479  mat_copy->add(1., *this);
480 
481  return mat_copy;
482 }
483 
484 template <typename T>
486 {
487  libmesh_assert (this->initialized());
488 
489  return static_cast<numeric_index_type>(Q_GetDim(const_cast<QMatrix*>(&_QMat)));
490 }
491 
492 
493 
494 template <typename T>
496 {
497  libmesh_assert (this->initialized());
498 
499  return static_cast<numeric_index_type>(Q_GetDim(const_cast<QMatrix*>(&_QMat)));
500 }
501 
502 
503 
504 template <typename T>
506 {
507  return 0;
508 }
509 
510 
511 
512 template <typename T>
514 {
515  return this->m();
516 }
517 
518 
519 
520 template <typename T>
522 {
523  return 0;
524 }
525 
526 
527 
528 template <typename T>
530 {
531  return this->n();
532 }
533 
534 
535 
536 template <typename T>
538  const numeric_index_type j,
539  const T value)
540 {
541  libmesh_assert (this->initialized());
542  libmesh_assert_less (i, this->m());
543  libmesh_assert_less (j, this->n());
544 
545  const numeric_index_type position = this->pos(i,j);
546 
547  // Sanity check
548  libmesh_assert_equal_to (*(_row_start[i]+position), j);
549  libmesh_assert_equal_to ((j+1), Q_GetPos (&_QMat, i+1, position));
550 
551  Q_SetEntry (&_QMat, i+1, position, j+1, value);
552 }
553 
554 
555 
556 template <typename T>
558  const numeric_index_type j,
559  const T value)
560 {
561  libmesh_assert (this->initialized());
562  libmesh_assert_less (i, this->m());
563  libmesh_assert_less (j, this->n());
564 
565  const numeric_index_type position = this->pos(i,j);
566 
567  // Sanity check
568  libmesh_assert_equal_to (*(_row_start[i]+position), j);
569 
570  Q_AddVal (&_QMat, i+1, position, value);
571 }
572 
573 
574 
575 template <typename T>
577  const std::vector<numeric_index_type> & dof_indices)
578 {
579  this->add_matrix (dm, dof_indices, dof_indices);
580 }
581 
582 
583 
584 template <typename T>
585 void LaspackMatrix<T>::add (const T a_in, const SparseMatrix<T> & X_in)
586 {
587  libmesh_assert (this->initialized());
588  libmesh_assert_equal_to (this->m(), X_in.m());
589  libmesh_assert_equal_to (this->n(), X_in.n());
590 
591  const LaspackMatrix<T> * X =
592  cast_ptr<const LaspackMatrix<T> *> (&X_in);
593 
594  _LPNumber a = static_cast<_LPNumber> (a_in);
595 
596  libmesh_assert(X);
597 
598  // loops taken from LaspackMatrix<T>::zero ()
599 
600  const numeric_index_type n_rows = this->m();
601 
602  for (numeric_index_type row=0; row<n_rows; row++)
603  {
604  auto r_start = _row_start[row];
605 
606  const numeric_index_type len = (_row_start[row+1] - _row_start[row]);
607 
608  // Make sure we agree on the row length
609  libmesh_assert_equal_to (len, Q_GetLen(&_QMat, row+1));
610  // compare matrix sparsity structures
611  libmesh_assert_equal_to (len, Q_GetLen(&(X->_QMat), row+1));
612 
613 
614  for (numeric_index_type l=0; l<len; l++)
615  {
616  const numeric_index_type j = *(r_start + l);
617 
618  // Make sure the data structures are working
619  libmesh_assert_equal_to ((j+1), Q_GetPos (&_QMat, row+1, l));
620 
621  const _LPNumber value = a * Q_GetEl(const_cast<QMatrix*>(&(X->_QMat)), row+1, j+1);
622  Q_AddVal (&_QMat, row+1, l, value);
623  }
624  }
625 }
626 
627 
628 
629 
630 template <typename T>
632  const numeric_index_type j) const
633 {
634  libmesh_assert (this->initialized());
635  libmesh_assert_less (i, this->m());
636  libmesh_assert_less (j, this->n());
637 
638  return Q_GetEl (const_cast<QMatrix*>(&_QMat), i+1, j+1);
639 }
640 
641 
642 
643 template <typename T>
645  const numeric_index_type j) const
646 {
647  libmesh_assert_less (i, this->m());
648  libmesh_assert_less (j, this->n());
649  libmesh_assert_less (i+1, _row_start.size());
650  libmesh_assert (_row_start.back() == _csr.end());
651 
652  // note this requires the _csr to be sorted
653  auto p = std::equal_range (_row_start[i], _row_start[i+1], j);
654 
655  // Make sure the row contains the element j
656  libmesh_assert (p.first != p.second);
657 
658  // Make sure the values match
659  libmesh_assert (*p.first == j);
660 
661  // Return the position in the compressed row
662  return std::distance (_row_start[i], p.first);
663 }
664 
665 
666 
667 template <typename T>
669 {
670  libmesh_assert(this->initialized());
671 
672  this->_closed = true;
673 
674  // We've probably changed some entries so we need to tell LASPACK
675  // that cached data is now invalid.
676  *_QMat.DiagElAlloc = _LPFalse;
677  *_QMat.ElSorted = _LPFalse;
678  if (*_QMat.ILUExists)
679  {
680  *_QMat.ILUExists = _LPFalse;
681  Q_Destr(_QMat.ILU);
682  }
683 }
684 
685 
686 
687 //------------------------------------------------------------------
688 // Explicit instantiations
689 template class LIBMESH_EXPORT LaspackMatrix<Number>;
690 
691 } // namespace libMesh
692 
693 
694 #endif // #ifdef LIBMESH_HAVE_LASPACK
virtual void init(const numeric_index_type m, const numeric_index_type n, const numeric_index_type m_l, const numeric_index_type n_l, const numeric_index_type nnz=30, const numeric_index_type noz=10, const numeric_index_type blocksize=1) override
Initialize SparseMatrix with the specified sizes.
virtual std::unique_ptr< SparseMatrix< T > > zero_clone() const override
virtual void get_diagonal(NumericVector< T > &dest) const override
Copies the diagonal part of the matrix into dest.
virtual numeric_index_type col_start() const override
Provides a uniform interface to vector storage schemes for different linear algebra libraries...
Definition: vector_fe_ex5.C:44
unsigned int m() const
The libMesh namespace provides an interface to certain functionality in the library.
Real distance(const Point &p)
virtual void add(const numeric_index_type i, const numeric_index_type j, const T value) override
Add value to the element (i,j).
Generic sparse matrix.
Definition: vector_fe_ex5.C:46
virtual std::unique_ptr< SparseMatrix< T > > clone() const override
numeric_index_type pos(const numeric_index_type i, const numeric_index_type j) const
virtual void get_transpose(SparseMatrix< T > &dest) const override
Copies the transpose of the matrix into dest, which may be *this.
virtual void close() override
Calls the SparseMatrix&#39;s internal assembly routines, ensuring that the values are consistent across p...
dof_id_type numeric_index_type
Definition: id_types.h:99
bool _is_initialized
Flag that tells if init() has been called.
Definition: libmesh.C:315
virtual Real l1_norm() const override
virtual numeric_index_type m() const =0
virtual T operator()(const numeric_index_type i, const numeric_index_type j) const override
void init(triangulateio &t)
Initializes the fields of t to nullptr/0 as necessary.
libmesh_assert(ctx)
virtual Real linfty_norm() const override
virtual numeric_index_type row_start() const override
virtual void zero() override
Set all entries to 0.
virtual void clear() override
Restores the SparseMatrix<T> to a pristine state.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual numeric_index_type m() const override
virtual numeric_index_type n() const override
static const bool value
Definition: xdr_io.C:55
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:140
bool initialized()
Checks that library initialization has been done.
Definition: libmesh.C:334
virtual void update_sparsity_pattern(const SparsityPattern::Graph &) override
Updates the matrix sparsity pattern.
virtual numeric_index_type row_stop() const override
The LaspackMatrix class wraps a QMatrix object from the Laspack library.
unsigned int n() const
virtual void set(const numeric_index_type i, const numeric_index_type j, const T value) override
Set the element (i,j) to value.
Defines a dense matrix for use in Finite Element-type computations.
Definition: dof_map.h:75
virtual void get_row(numeric_index_type i, std::vector< numeric_index_type > &indices, std::vector< T > &values) const override
Get a row from the matrix.
virtual numeric_index_type n() const =0
virtual numeric_index_type col_stop() const override
virtual void add_matrix(const DenseMatrix< T > &dm, const std::vector< numeric_index_type > &rows, const std::vector< numeric_index_type > &cols) override
Add the full matrix dm to the SparseMatrix.
LaspackMatrix(const Parallel::Communicator &comm)
Constructor; initializes the matrix to be empty, without any structure, i.e.
uint8_t dof_id_type
Definition: id_types.h:67
ParallelType
Defines an enum for parallel data structure types.