Line data Source code
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>
42 0 : void LaspackMatrix<T>::update_sparsity_pattern (const SparsityPattern::Graph & sparsity_pattern)
43 : {
44 : // clear data, start over
45 0 : this->clear ();
46 :
47 : // big trouble if this fails!
48 0 : libmesh_assert(this->_dof_map);
49 :
50 0 : 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 0 : std::size_t size = 0;
56 :
57 0 : for (numeric_index_type row=0; row<n_rows; row++)
58 0 : size += sparsity_pattern[row].size();
59 :
60 0 : _csr.resize (size);
61 0 : _row_start.reserve(n_rows + 1);
62 : }
63 :
64 :
65 : // Initialize the _csr data structure.
66 : {
67 0 : std::vector<numeric_index_type>::iterator position = _csr.begin();
68 :
69 0 : _row_start.push_back (position);
70 :
71 0 : for (numeric_index_type row=0; row<n_rows; row++)
72 : {
73 : // insert the row indices
74 0 : for (const auto & col : sparsity_pattern[row])
75 : {
76 0 : libmesh_assert (position != _csr.end());
77 0 : *position = col;
78 0 : ++position;
79 : }
80 :
81 0 : _row_start.push_back (position);
82 : }
83 : }
84 :
85 :
86 : // Initialize the matrix
87 0 : libmesh_assert (!this->initialized());
88 0 : this->init ();
89 0 : libmesh_assert (this->initialized());
90 : //libMesh::out << "n_rows=" << n_rows << std::endl;
91 : //libMesh::out << "m()=" << m() << std::endl;
92 0 : libmesh_assert_equal_to (n_rows, this->m());
93 :
94 : // Tell the matrix about its structure. Initialize it
95 : // to zero.
96 0 : for (numeric_index_type i=0; i<n_rows; i++)
97 : {
98 0 : auto rs = _row_start[i];
99 :
100 0 : const numeric_index_type length = _row_start[i+1] - rs;
101 :
102 0 : Q_SetLen (&_QMat, i+1, length);
103 :
104 0 : for (numeric_index_type l=0; l<length; l++)
105 : {
106 0 : 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 0 : libmesh_assert_equal_to (this->pos(i,j), l);
117 0 : Q_SetEntry (&_QMat, i+1, l, j+1, 0.);
118 : }
119 : }
120 :
121 : // That's it!
122 : //libmesh_here();
123 0 : }
124 :
125 :
126 :
127 : template <typename T>
128 44 : void LaspackMatrix<T>::init (const numeric_index_type m_in,
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 22 : libmesh_assert_equal_to (m_in, m_l);
138 22 : libmesh_assert_equal_to (n_in, n_l);
139 22 : libmesh_assert_equal_to (m_in, n_in);
140 22 : libmesh_assert_greater (nnz, 0);
141 :
142 66 : if ((m_in != m_l) ||
143 44 : (n_in != n_l))
144 0 : libmesh_not_implemented_msg("Laspack does not support distributed matrices");
145 :
146 44 : if (m_in != n_in)
147 0 : libmesh_not_implemented_msg("Laspack does not support rectangular matrices");
148 :
149 22 : if (nnz < n_in)
150 : libmesh_warning("Using inefficient LaspackMatrix allocation via this init() method");
151 :
152 22 : 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 44 : _csr.resize (m_in * m_in);
158 44 : _row_start.reserve(m_in + 1);
159 : }
160 :
161 : // Initialize the _csr data structure.
162 : {
163 22 : std::vector<numeric_index_type>::iterator position = _csr.begin();
164 :
165 44 : _row_start.push_back (position);
166 :
167 220 : for (numeric_index_type row=0; row<n_rows; row++)
168 : {
169 : // insert the row indices
170 880 : for (const auto & col : make_range(n_rows))
171 : {
172 352 : libmesh_assert (position != _csr.end());
173 704 : *position = col;
174 352 : ++position;
175 : }
176 :
177 176 : _row_start.push_back (position);
178 : }
179 : }
180 :
181 44 : Q_Constr(&_QMat, const_cast<char *>("Mat"), m_in, _LPFalse, Rowws, Normal, _LPTrue);
182 :
183 44 : this->_is_initialized = true;
184 :
185 : // Tell the matrix about its structure. Initialize it
186 : // to zero.
187 220 : for (numeric_index_type i=0; i<n_rows; i++)
188 : {
189 176 : auto rs = _row_start[i];
190 :
191 176 : const numeric_index_type length = _row_start[i+1] - rs;
192 :
193 176 : Q_SetLen (&_QMat, i+1, length);
194 :
195 880 : for (numeric_index_type l=0; l<length; l++)
196 : {
197 704 : const numeric_index_type j = *(rs+l);
198 :
199 352 : libmesh_assert_equal_to (this->pos(i,j), l);
200 704 : Q_SetEntry (&_QMat, i+1, l, j+1, 0.);
201 : }
202 : }
203 44 : }
204 :
205 :
206 :
207 : template <typename T>
208 0 : void LaspackMatrix<T>::init (const ParallelType)
209 : {
210 : // Ignore calls on initialized objects
211 0 : if (this->initialized())
212 0 : return;
213 :
214 : // We need a sparsity pattern for this!
215 0 : libmesh_assert(this->_sp);
216 :
217 : // Clear initialized matrices
218 0 : if (this->initialized())
219 0 : this->clear();
220 :
221 0 : 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 0 : const numeric_index_type n_cols = n_rows;
226 0 : const numeric_index_type n_l = this->_dof_map->n_dofs_on_processor(0);
227 0 : const numeric_index_type m_l = n_l;
228 : #endif
229 :
230 : // Laspack Matrices only work for uniprocessor cases
231 0 : libmesh_assert_equal_to (n_rows, n_cols);
232 0 : libmesh_assert_equal_to (m_l, n_rows);
233 0 : 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 0 : const std::vector<numeric_index_type> & n_nz = this->_sp->get_n_nz();
239 0 : 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 0 : libmesh_assert_equal_to (n_nz.size(), n_l);
244 0 : libmesh_assert_equal_to (n_oz.size(), n_l);
245 :
246 0 : if (n_rows==0)
247 0 : return;
248 :
249 0 : Q_Constr(&_QMat, const_cast<char *>("Mat"), n_rows, _LPFalse, Rowws, Normal, _LPTrue);
250 :
251 0 : this->_is_initialized = true;
252 :
253 0 : libmesh_assert_equal_to (n_rows, this->m());
254 : }
255 :
256 :
257 :
258 : template <typename T>
259 12 : void LaspackMatrix<T>::add_matrix(const DenseMatrix<T> & dm,
260 : const std::vector<numeric_index_type> & rows,
261 : const std::vector<numeric_index_type> & cols)
262 :
263 : {
264 6 : libmesh_assert (this->initialized());
265 12 : unsigned int n_rows = cast_int<unsigned int>(rows.size());
266 12 : unsigned int n_cols = cast_int<unsigned int>(cols.size());
267 6 : libmesh_assert_equal_to (dm.m(), n_rows);
268 6 : libmesh_assert_equal_to (dm.n(), n_cols);
269 :
270 :
271 60 : for (unsigned int i=0; i<n_rows; i++)
272 240 : for (unsigned int j=0; j<n_cols; j++)
273 384 : this->add(rows[i],cols[j],dm(i,j));
274 12 : }
275 :
276 :
277 :
278 : template <typename T>
279 36 : Real LaspackMatrix<T>::l1_norm () const
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 36 : std::vector<Real> abs_col_sums(this->n());
285 :
286 36 : const numeric_index_type n_rows = this->m();
287 :
288 180 : for (numeric_index_type row : make_range(n_rows))
289 : {
290 144 : auto r_start = _row_start[row];
291 :
292 216 : const numeric_index_type len = (_row_start[row+1] - _row_start[row]);
293 :
294 : // Make sure we agree on the row length
295 72 : libmesh_assert_equal_to (len, Q_GetLen(&_QMat, row+1));
296 :
297 720 : for (numeric_index_type l=0; l<len; l++)
298 : {
299 576 : const numeric_index_type j = *(r_start + l);
300 :
301 : // Make sure the data structures are working
302 288 : libmesh_assert_equal_to ((j+1), Q_GetPos (&_QMat, row+1, l));
303 :
304 576 : abs_col_sums[j] += std::abs(Q_GetVal (&_QMat, row+1, l));
305 : }
306 : }
307 :
308 72 : return *(std::max_element(abs_col_sums.begin(), abs_col_sums.end()));
309 : }
310 :
311 :
312 :
313 : template <typename T>
314 0 : void LaspackMatrix<T>::get_diagonal (NumericVector<T> & /*dest*/) const
315 : {
316 0 : libmesh_not_implemented();
317 : }
318 :
319 :
320 :
321 : template <typename T>
322 0 : void LaspackMatrix<T>::get_transpose (SparseMatrix<T> & /*dest*/) const
323 : {
324 0 : libmesh_not_implemented();
325 : }
326 :
327 :
328 :
329 : template <typename T>
330 64 : void LaspackMatrix<T>::get_row(numeric_index_type i,
331 : std::vector<numeric_index_type> & indices,
332 : std::vector<T> & values) const
333 : {
334 32 : indices.clear();
335 32 : values.clear();
336 :
337 128 : const numeric_index_type len = (_row_start[i+1] - _row_start[i]);
338 :
339 : // Make sure we agree on the row length
340 32 : libmesh_assert_equal_to (len, Q_GetLen(&_QMat, i+1));
341 :
342 32 : auto r_start = _row_start[i];
343 :
344 320 : for (numeric_index_type l=0; l<len; l++)
345 : {
346 256 : const numeric_index_type j = *(r_start + l);
347 :
348 : // Make sure the data structures are working
349 128 : libmesh_assert_equal_to ((j+1), Q_GetPos (&_QMat, i+1, l));
350 :
351 256 : indices.push_back(j);
352 :
353 256 : values.push_back(Q_GetVal (&_QMat, i+1, l));
354 : }
355 64 : }
356 :
357 :
358 :
359 : template <typename T>
360 40 : LaspackMatrix<T>::LaspackMatrix (const Parallel::Communicator & comm) :
361 : SparseMatrix<T>(comm),
362 40 : _closed (false)
363 : {
364 40 : }
365 :
366 :
367 :
368 : template <typename T>
369 80 : LaspackMatrix<T>::~LaspackMatrix ()
370 : {
371 40 : this->clear ();
372 100 : }
373 :
374 :
375 :
376 : template <typename T>
377 44 : void LaspackMatrix<T>::clear ()
378 : {
379 44 : if (this->initialized())
380 : {
381 44 : Q_Destr(&_QMat);
382 : }
383 :
384 22 : _csr.clear();
385 22 : _row_start.clear();
386 44 : _closed = false;
387 44 : this->_is_initialized = false;
388 44 : }
389 :
390 :
391 :
392 : template <typename T>
393 28 : void LaspackMatrix<T>::zero ()
394 : {
395 28 : const numeric_index_type n_rows = this->m();
396 :
397 140 : for (numeric_index_type row=0; row<n_rows; row++)
398 : {
399 112 : auto r_start = _row_start[row];
400 :
401 168 : const numeric_index_type len = (_row_start[row+1] - _row_start[row]);
402 :
403 : // Make sure we agree on the row length
404 56 : libmesh_assert_equal_to (len, Q_GetLen(&_QMat, row+1));
405 :
406 560 : for (numeric_index_type l=0; l<len; l++)
407 : {
408 448 : const numeric_index_type j = *(r_start + l);
409 :
410 : // Make sure the data structures are working
411 224 : libmesh_assert_equal_to ((j+1), Q_GetPos (&_QMat, row+1, l));
412 :
413 448 : Q_SetEntry (&_QMat, row+1, l, j+1, 0.);
414 : }
415 : }
416 :
417 28 : this->close();
418 28 : }
419 :
420 :
421 :
422 : template <typename T>
423 16 : std::unique_ptr<SparseMatrix<T>> LaspackMatrix<T>::zero_clone () const
424 : {
425 : // Make empty copy with matching comm, initialize, zero, and return.
426 24 : auto mat_copy = std::make_unique<LaspackMatrix<T>>(this->comm());
427 16 : if (this->_dof_map)
428 0 : mat_copy->attach_dof_map(*this->_dof_map);
429 : else
430 16 : mat_copy->init(this->m(), this->n(), this->local_m(),
431 : this->local_n(), this->local_n());
432 :
433 16 : mat_copy->zero();
434 :
435 24 : return mat_copy;
436 : }
437 :
438 :
439 :
440 : template <typename T>
441 12 : std::unique_ptr<SparseMatrix<T>> LaspackMatrix<T>::clone () const
442 : {
443 : // We don't currently have a faster implementation than making a
444 : // zero clone and then filling in the values.
445 12 : auto mat_copy = this->zero_clone();
446 12 : mat_copy->add(1., *this);
447 :
448 12 : return mat_copy;
449 : }
450 :
451 : template <typename T>
452 944 : numeric_index_type LaspackMatrix<T>::m () const
453 : {
454 850 : libmesh_assert (this->initialized());
455 :
456 944 : return static_cast<numeric_index_type>(Q_GetDim(const_cast<QMatrix*>(&_QMat)));
457 : }
458 :
459 :
460 :
461 : template <typename T>
462 890 : numeric_index_type LaspackMatrix<T>::n () const
463 : {
464 823 : libmesh_assert (this->initialized());
465 :
466 890 : return static_cast<numeric_index_type>(Q_GetDim(const_cast<QMatrix*>(&_QMat)));
467 : }
468 :
469 :
470 :
471 : template <typename T>
472 338 : numeric_index_type LaspackMatrix<T>::row_start () const
473 : {
474 338 : return 0;
475 : }
476 :
477 :
478 :
479 : template <typename T>
480 70 : numeric_index_type LaspackMatrix<T>::row_stop () const
481 : {
482 70 : return this->m();
483 : }
484 :
485 :
486 :
487 : template <typename T>
488 76 : numeric_index_type LaspackMatrix<T>::col_start () const
489 : {
490 76 : return 0;
491 : }
492 :
493 :
494 :
495 : template <typename T>
496 48 : numeric_index_type LaspackMatrix<T>::col_stop () const
497 : {
498 48 : return this->n();
499 : }
500 :
501 :
502 :
503 : template <typename T>
504 128 : void LaspackMatrix<T>::set (const numeric_index_type i,
505 : const numeric_index_type j,
506 : const T value)
507 : {
508 64 : libmesh_assert (this->initialized());
509 64 : libmesh_assert_less (i, this->m());
510 64 : libmesh_assert_less (j, this->n());
511 :
512 128 : const numeric_index_type position = this->pos(i,j);
513 :
514 : // Sanity check
515 64 : libmesh_assert_equal_to (*(_row_start[i]+position), j);
516 64 : libmesh_assert_equal_to ((j+1), Q_GetPos (&_QMat, i+1, position));
517 :
518 128 : Q_SetEntry (&_QMat, i+1, position, j+1, value);
519 128 : }
520 :
521 :
522 :
523 : template <typename T>
524 192 : void LaspackMatrix<T>::add (const numeric_index_type i,
525 : const numeric_index_type j,
526 : const T value)
527 : {
528 96 : libmesh_assert (this->initialized());
529 96 : libmesh_assert_less (i, this->m());
530 96 : libmesh_assert_less (j, this->n());
531 :
532 192 : const numeric_index_type position = this->pos(i,j);
533 :
534 : // Sanity check
535 96 : libmesh_assert_equal_to (*(_row_start[i]+position), j);
536 :
537 192 : Q_AddVal (&_QMat, i+1, position, value);
538 192 : }
539 :
540 :
541 :
542 : template <typename T>
543 0 : void LaspackMatrix<T>::add_matrix(const DenseMatrix<T> & dm,
544 : const std::vector<numeric_index_type> & dof_indices)
545 : {
546 0 : this->add_matrix (dm, dof_indices, dof_indices);
547 0 : }
548 :
549 :
550 :
551 : template <typename T>
552 20 : void LaspackMatrix<T>::add (const T a_in, const SparseMatrix<T> & X_in)
553 : {
554 10 : libmesh_assert (this->initialized());
555 10 : libmesh_assert_equal_to (this->m(), X_in.m());
556 10 : libmesh_assert_equal_to (this->n(), X_in.n());
557 :
558 10 : const LaspackMatrix<T> * X =
559 10 : cast_ptr<const LaspackMatrix<T> *> (&X_in);
560 :
561 10 : _LPNumber a = static_cast<_LPNumber> (a_in);
562 :
563 10 : libmesh_assert(X);
564 :
565 : // loops taken from LaspackMatrix<T>::zero ()
566 :
567 20 : const numeric_index_type n_rows = this->m();
568 :
569 100 : for (numeric_index_type row=0; row<n_rows; row++)
570 : {
571 80 : auto r_start = _row_start[row];
572 :
573 120 : const numeric_index_type len = (_row_start[row+1] - _row_start[row]);
574 :
575 : // Make sure we agree on the row length
576 40 : libmesh_assert_equal_to (len, Q_GetLen(&_QMat, row+1));
577 : // compare matrix sparsity structures
578 40 : libmesh_assert_equal_to (len, Q_GetLen(&(X->_QMat), row+1));
579 :
580 :
581 400 : for (numeric_index_type l=0; l<len; l++)
582 : {
583 320 : const numeric_index_type j = *(r_start + l);
584 :
585 : // Make sure the data structures are working
586 160 : libmesh_assert_equal_to ((j+1), Q_GetPos (&_QMat, row+1, l));
587 :
588 320 : const _LPNumber value = a * Q_GetEl(const_cast<QMatrix*>(&(X->_QMat)), row+1, j+1);
589 320 : Q_AddVal (&_QMat, row+1, l, value);
590 : }
591 : }
592 20 : }
593 :
594 :
595 :
596 :
597 : template <typename T>
598 128 : T LaspackMatrix<T>::operator () (const numeric_index_type i,
599 : const numeric_index_type j) const
600 : {
601 64 : libmesh_assert (this->initialized());
602 64 : libmesh_assert_less (i, this->m());
603 64 : libmesh_assert_less (j, this->n());
604 :
605 128 : return Q_GetEl (const_cast<QMatrix*>(&_QMat), i+1, j+1);
606 : }
607 :
608 :
609 :
610 : template <typename T>
611 672 : numeric_index_type LaspackMatrix<T>::pos (const numeric_index_type i,
612 : const numeric_index_type j) const
613 : {
614 512 : libmesh_assert_less (i, this->m());
615 512 : libmesh_assert_less (j, this->n());
616 512 : libmesh_assert_less (i+1, _row_start.size());
617 512 : libmesh_assert (_row_start.back() == _csr.end());
618 :
619 : // note this requires the _csr to be sorted
620 832 : auto p = std::equal_range (_row_start[i], _row_start[i+1], j);
621 :
622 : // Make sure the row contains the element j
623 512 : libmesh_assert (p.first != p.second);
624 :
625 : // Make sure the values match
626 512 : libmesh_assert (*p.first == j);
627 :
628 : // Return the position in the compressed row
629 832 : return std::distance (_row_start[i], p.first);
630 : }
631 :
632 :
633 :
634 : template <typename T>
635 48 : void LaspackMatrix<T>::close()
636 : {
637 24 : libmesh_assert(this->initialized());
638 :
639 48 : this->_closed = true;
640 :
641 : // We've probably changed some entries so we need to tell LASPACK
642 : // that cached data is now invalid.
643 48 : *_QMat.DiagElAlloc = _LPFalse;
644 48 : *_QMat.ElSorted = _LPFalse;
645 48 : if (*_QMat.ILUExists)
646 : {
647 0 : *_QMat.ILUExists = _LPFalse;
648 0 : Q_Destr(_QMat.ILU);
649 : }
650 48 : }
651 :
652 :
653 :
654 : //------------------------------------------------------------------
655 : // Explicit instantiations
656 : template class LIBMESH_EXPORT LaspackMatrix<Number>;
657 :
658 : } // namespace libMesh
659 :
660 :
661 : #endif // #ifdef LIBMESH_HAVE_LASPACK
|
