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