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 48 : 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 24 : libmesh_assert_equal_to (m_in, m_l);
138 24 : libmesh_assert_equal_to (n_in, n_l);
139 24 : libmesh_assert_equal_to (m_in, n_in);
140 24 : libmesh_assert_greater (nnz, 0);
141 :
142 72 : if ((m_in != m_l) ||
143 48 : (n_in != n_l))
144 0 : libmesh_not_implemented_msg("Laspack does not support distributed matrices");
145 :
146 48 : if (m_in != n_in)
147 0 : libmesh_not_implemented_msg("Laspack does not support rectangular matrices");
148 :
149 24 : if (nnz < n_in)
150 : libmesh_warning("Using inefficient LaspackMatrix allocation via this init() method");
151 :
152 24 : 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 48 : _csr.resize (m_in * m_in);
158 48 : _row_start.reserve(m_in + 1);
159 : }
160 :
161 : // Initialize the _csr data structure.
162 : {
163 24 : std::vector<numeric_index_type>::iterator position = _csr.begin();
164 :
165 48 : _row_start.push_back (position);
166 :
167 240 : for (numeric_index_type row=0; row<n_rows; row++)
168 : {
169 : // insert the row indices
170 960 : for (const auto & col : make_range(n_rows))
171 : {
172 384 : libmesh_assert (position != _csr.end());
173 768 : *position = col;
174 384 : ++position;
175 : }
176 :
177 192 : _row_start.push_back (position);
178 : }
179 : }
180 :
181 48 : Q_Constr(&_QMat, const_cast<char *>("Mat"), m_in, _LPFalse, Rowws, Normal, _LPTrue);
182 :
183 48 : this->_is_initialized = true;
184 :
185 : // Tell the matrix about its structure. Initialize it
186 : // to zero.
187 240 : for (numeric_index_type i=0; i<n_rows; i++)
188 : {
189 192 : auto rs = _row_start[i];
190 :
191 192 : const numeric_index_type length = _row_start[i+1] - rs;
192 :
193 192 : Q_SetLen (&_QMat, i+1, length);
194 :
195 960 : for (numeric_index_type l=0; l<length; l++)
196 : {
197 768 : const numeric_index_type j = *(rs+l);
198 :
199 384 : libmesh_assert_equal_to (this->pos(i,j), l);
200 768 : Q_SetEntry (&_QMat, i+1, l, j+1, 0.);
201 : }
202 : }
203 48 : }
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 : Real LaspackMatrix<T>::linfty_norm () const
315 : {
316 0 : const numeric_index_type n_rows = this->m();
317 :
318 0 : Real max_row_sum = 0;
319 :
320 0 : for (numeric_index_type row : make_range(n_rows))
321 : {
322 0 : Real row_sum = 0;
323 :
324 0 : const numeric_index_type len = (_row_start[row+1] - _row_start[row]);
325 :
326 : // Make sure we agree on the row length
327 0 : libmesh_assert_equal_to (len, Q_GetLen(&_QMat, row+1));
328 :
329 0 : for (numeric_index_type l=0; l<len; l++)
330 : {
331 : // Make sure the data structures are working
332 0 : libmesh_assert_equal_to ((*(_row_start[row] + l)+1), Q_GetPos
333 : (&_QMat, row+1, l));
334 :
335 0 : row_sum += std::abs(Q_GetVal (&_QMat, row+1, l));
336 : }
337 :
338 0 : max_row_sum = std::max(max_row_sum, row_sum);
339 : }
340 :
341 0 : return max_row_sum;
342 : }
343 :
344 :
345 :
346 : template <typename T>
347 0 : void LaspackMatrix<T>::get_diagonal (NumericVector<T> & /*dest*/) const
348 : {
349 0 : libmesh_not_implemented();
350 : }
351 :
352 :
353 :
354 : template <typename T>
355 0 : void LaspackMatrix<T>::get_transpose (SparseMatrix<T> & /*dest*/) const
356 : {
357 0 : libmesh_not_implemented();
358 : }
359 :
360 :
361 :
362 : template <typename T>
363 64 : void LaspackMatrix<T>::get_row(numeric_index_type i,
364 : std::vector<numeric_index_type> & indices,
365 : std::vector<T> & values) const
366 : {
367 32 : indices.clear();
368 32 : values.clear();
369 :
370 128 : const numeric_index_type len = (_row_start[i+1] - _row_start[i]);
371 :
372 : // Make sure we agree on the row length
373 32 : libmesh_assert_equal_to (len, Q_GetLen(&_QMat, i+1));
374 :
375 32 : auto r_start = _row_start[i];
376 :
377 320 : for (numeric_index_type l=0; l<len; l++)
378 : {
379 256 : const numeric_index_type j = *(r_start + l);
380 :
381 : // Make sure the data structures are working
382 128 : libmesh_assert_equal_to ((j+1), Q_GetPos (&_QMat, i+1, l));
383 :
384 256 : indices.push_back(j);
385 :
386 256 : values.push_back(Q_GetVal (&_QMat, i+1, l));
387 : }
388 64 : }
389 :
390 :
391 :
392 : template <typename T>
393 44 : LaspackMatrix<T>::LaspackMatrix (const Parallel::Communicator & comm) :
394 : SparseMatrix<T>(comm),
395 44 : _closed (false)
396 : {
397 44 : }
398 :
399 :
400 :
401 : template <typename T>
402 88 : LaspackMatrix<T>::~LaspackMatrix ()
403 : {
404 44 : this->clear ();
405 110 : }
406 :
407 :
408 :
409 : template <typename T>
410 48 : void LaspackMatrix<T>::clear ()
411 : {
412 48 : if (this->initialized())
413 : {
414 48 : Q_Destr(&_QMat);
415 : }
416 :
417 24 : _csr.clear();
418 24 : _row_start.clear();
419 48 : _closed = false;
420 48 : this->_is_initialized = false;
421 48 : }
422 :
423 :
424 :
425 : template <typename T>
426 28 : void LaspackMatrix<T>::zero ()
427 : {
428 28 : const numeric_index_type n_rows = this->m();
429 :
430 140 : for (numeric_index_type row=0; row<n_rows; row++)
431 : {
432 112 : auto r_start = _row_start[row];
433 :
434 168 : const numeric_index_type len = (_row_start[row+1] - _row_start[row]);
435 :
436 : // Make sure we agree on the row length
437 56 : libmesh_assert_equal_to (len, Q_GetLen(&_QMat, row+1));
438 :
439 560 : for (numeric_index_type l=0; l<len; l++)
440 : {
441 448 : const numeric_index_type j = *(r_start + l);
442 :
443 : // Make sure the data structures are working
444 224 : libmesh_assert_equal_to ((j+1), Q_GetPos (&_QMat, row+1, l));
445 :
446 448 : Q_SetEntry (&_QMat, row+1, l, j+1, 0.);
447 : }
448 : }
449 :
450 28 : this->close();
451 28 : }
452 :
453 :
454 :
455 : template <typename T>
456 16 : std::unique_ptr<SparseMatrix<T>> LaspackMatrix<T>::zero_clone () const
457 : {
458 : // Make empty copy with matching comm, initialize, zero, and return.
459 24 : auto mat_copy = std::make_unique<LaspackMatrix<T>>(this->comm());
460 16 : if (this->_dof_map)
461 0 : mat_copy->attach_dof_map(*this->_dof_map);
462 : else
463 16 : mat_copy->init(this->m(), this->n(), this->local_m(),
464 : this->local_n(), this->local_n());
465 :
466 16 : mat_copy->zero();
467 :
468 24 : return mat_copy;
469 : }
470 :
471 :
472 :
473 : template <typename T>
474 12 : 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 12 : auto mat_copy = this->zero_clone();
479 12 : mat_copy->add(1., *this);
480 :
481 12 : return mat_copy;
482 : }
483 :
484 : template <typename T>
485 976 : numeric_index_type LaspackMatrix<T>::m () const
486 : {
487 882 : libmesh_assert (this->initialized());
488 :
489 976 : return static_cast<numeric_index_type>(Q_GetDim(const_cast<QMatrix*>(&_QMat)));
490 : }
491 :
492 :
493 :
494 : template <typename T>
495 922 : numeric_index_type LaspackMatrix<T>::n () const
496 : {
497 855 : libmesh_assert (this->initialized());
498 :
499 922 : return static_cast<numeric_index_type>(Q_GetDim(const_cast<QMatrix*>(&_QMat)));
500 : }
501 :
502 :
503 :
504 : template <typename T>
505 338 : numeric_index_type LaspackMatrix<T>::row_start () const
506 : {
507 338 : return 0;
508 : }
509 :
510 :
511 :
512 : template <typename T>
513 70 : numeric_index_type LaspackMatrix<T>::row_stop () const
514 : {
515 70 : return this->m();
516 : }
517 :
518 :
519 :
520 : template <typename T>
521 76 : numeric_index_type LaspackMatrix<T>::col_start () const
522 : {
523 76 : return 0;
524 : }
525 :
526 :
527 :
528 : template <typename T>
529 48 : numeric_index_type LaspackMatrix<T>::col_stop () const
530 : {
531 48 : return this->n();
532 : }
533 :
534 :
535 :
536 : template <typename T>
537 128 : void LaspackMatrix<T>::set (const numeric_index_type i,
538 : const numeric_index_type j,
539 : const T value)
540 : {
541 64 : libmesh_assert (this->initialized());
542 64 : libmesh_assert_less (i, this->m());
543 64 : libmesh_assert_less (j, this->n());
544 :
545 128 : const numeric_index_type position = this->pos(i,j);
546 :
547 : // Sanity check
548 64 : libmesh_assert_equal_to (*(_row_start[i]+position), j);
549 64 : libmesh_assert_equal_to ((j+1), Q_GetPos (&_QMat, i+1, position));
550 :
551 128 : Q_SetEntry (&_QMat, i+1, position, j+1, value);
552 128 : }
553 :
554 :
555 :
556 : template <typename T>
557 192 : void LaspackMatrix<T>::add (const numeric_index_type i,
558 : const numeric_index_type j,
559 : const T value)
560 : {
561 96 : libmesh_assert (this->initialized());
562 96 : libmesh_assert_less (i, this->m());
563 96 : libmesh_assert_less (j, this->n());
564 :
565 192 : const numeric_index_type position = this->pos(i,j);
566 :
567 : // Sanity check
568 96 : libmesh_assert_equal_to (*(_row_start[i]+position), j);
569 :
570 192 : Q_AddVal (&_QMat, i+1, position, value);
571 192 : }
572 :
573 :
574 :
575 : template <typename T>
576 0 : void LaspackMatrix<T>::add_matrix(const DenseMatrix<T> & dm,
577 : const std::vector<numeric_index_type> & dof_indices)
578 : {
579 0 : this->add_matrix (dm, dof_indices, dof_indices);
580 0 : }
581 :
582 :
583 :
584 : template <typename T>
585 20 : void LaspackMatrix<T>::add (const T a_in, const SparseMatrix<T> & X_in)
586 : {
587 10 : libmesh_assert (this->initialized());
588 10 : libmesh_assert_equal_to (this->m(), X_in.m());
589 10 : libmesh_assert_equal_to (this->n(), X_in.n());
590 :
591 10 : const LaspackMatrix<T> * X =
592 10 : cast_ptr<const LaspackMatrix<T> *> (&X_in);
593 :
594 10 : _LPNumber a = static_cast<_LPNumber> (a_in);
595 :
596 10 : libmesh_assert(X);
597 :
598 : // loops taken from LaspackMatrix<T>::zero ()
599 :
600 20 : const numeric_index_type n_rows = this->m();
601 :
602 100 : for (numeric_index_type row=0; row<n_rows; row++)
603 : {
604 80 : auto r_start = _row_start[row];
605 :
606 120 : const numeric_index_type len = (_row_start[row+1] - _row_start[row]);
607 :
608 : // Make sure we agree on the row length
609 40 : libmesh_assert_equal_to (len, Q_GetLen(&_QMat, row+1));
610 : // compare matrix sparsity structures
611 40 : libmesh_assert_equal_to (len, Q_GetLen(&(X->_QMat), row+1));
612 :
613 :
614 400 : for (numeric_index_type l=0; l<len; l++)
615 : {
616 320 : const numeric_index_type j = *(r_start + l);
617 :
618 : // Make sure the data structures are working
619 160 : libmesh_assert_equal_to ((j+1), Q_GetPos (&_QMat, row+1, l));
620 :
621 320 : const _LPNumber value = a * Q_GetEl(const_cast<QMatrix*>(&(X->_QMat)), row+1, j+1);
622 320 : Q_AddVal (&_QMat, row+1, l, value);
623 : }
624 : }
625 20 : }
626 :
627 :
628 :
629 :
630 : template <typename T>
631 128 : T LaspackMatrix<T>::operator () (const numeric_index_type i,
632 : const numeric_index_type j) const
633 : {
634 64 : libmesh_assert (this->initialized());
635 64 : libmesh_assert_less (i, this->m());
636 64 : libmesh_assert_less (j, this->n());
637 :
638 128 : return Q_GetEl (const_cast<QMatrix*>(&_QMat), i+1, j+1);
639 : }
640 :
641 :
642 :
643 : template <typename T>
644 704 : numeric_index_type LaspackMatrix<T>::pos (const numeric_index_type i,
645 : const numeric_index_type j) const
646 : {
647 544 : libmesh_assert_less (i, this->m());
648 544 : libmesh_assert_less (j, this->n());
649 544 : libmesh_assert_less (i+1, _row_start.size());
650 544 : libmesh_assert (_row_start.back() == _csr.end());
651 :
652 : // note this requires the _csr to be sorted
653 864 : auto p = std::equal_range (_row_start[i], _row_start[i+1], j);
654 :
655 : // Make sure the row contains the element j
656 544 : libmesh_assert (p.first != p.second);
657 :
658 : // Make sure the values match
659 544 : libmesh_assert (*p.first == j);
660 :
661 : // Return the position in the compressed row
662 864 : return std::distance (_row_start[i], p.first);
663 : }
664 :
665 :
666 :
667 : template <typename T>
668 48 : void LaspackMatrix<T>::close()
669 : {
670 24 : libmesh_assert(this->initialized());
671 :
672 48 : 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 48 : *_QMat.DiagElAlloc = _LPFalse;
677 48 : *_QMat.ElSorted = _LPFalse;
678 48 : if (*_QMat.ILUExists)
679 : {
680 0 : *_QMat.ILUExists = _LPFalse;
681 0 : Q_Destr(_QMat.ILU);
682 : }
683 48 : }
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
|
