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 : // C++ includes
21 : #include <algorithm> // for std::min
22 : #include <limits>
23 :
24 : // Local Includes
25 : #include "libmesh/dense_subvector.h"
26 : #include "libmesh/dense_vector.h"
27 : #include "libmesh/laspack_vector.h"
28 : #include "libmesh/laspack_matrix.h"
29 : #include "libmesh/int_range.h"
30 :
31 : #ifdef LIBMESH_HAVE_LASPACK
32 :
33 : namespace libMesh
34 : {
35 :
36 : template <typename T>
37 8 : T LaspackVector<T>::sum () const
38 : {
39 4 : libmesh_assert (this->closed());
40 :
41 4 : T _sum = 0;
42 :
43 4 : const numeric_index_type n = this->size();
44 :
45 88 : for (numeric_index_type i=0; i!=n; ++i)
46 80 : _sum += (*this)(i);
47 :
48 8 : return _sum;
49 : }
50 :
51 :
52 :
53 : template <typename T>
54 40 : Real LaspackVector<T>::l1_norm () const
55 : {
56 20 : libmesh_assert (this->closed());
57 :
58 40 : return static_cast<Real>(l1Norm_V(const_cast<QVector*>(&_vec)));
59 : }
60 :
61 :
62 :
63 : template <typename T>
64 16 : Real LaspackVector<T>::l2_norm () const
65 : {
66 8 : libmesh_assert (this->closed());
67 :
68 16 : return static_cast<Real>(l2Norm_V(const_cast<QVector*>(&_vec)));
69 : }
70 :
71 :
72 :
73 : template <typename T>
74 8 : Real LaspackVector<T>::linfty_norm () const
75 : {
76 4 : libmesh_assert (this->closed());
77 :
78 8 : return static_cast<Real>(MaxNorm_V(const_cast<QVector*>(&_vec)));
79 : }
80 :
81 :
82 :
83 : template <typename T>
84 8 : NumericVector<T> & LaspackVector<T>::operator += (const NumericVector<T> & v)
85 : {
86 4 : libmesh_assert (this->closed());
87 :
88 8 : this->add(1., v);
89 :
90 8 : return *this;
91 : }
92 :
93 :
94 :
95 :
96 : template <typename T>
97 16 : NumericVector<T> & LaspackVector<T>::operator -= (const NumericVector<T> & v)
98 : {
99 8 : libmesh_assert (this->closed());
100 :
101 16 : this->add(-1., v);
102 :
103 16 : return *this;
104 : }
105 :
106 :
107 :
108 : template <typename T>
109 8 : NumericVector<T> & LaspackVector<T>::operator *= (const NumericVector<T> & v)
110 : {
111 4 : libmesh_assert_equal_to(size(), v.size());
112 :
113 : #ifndef NDEBUG
114 4 : const bool was_closed = this->_is_closed;
115 : #endif
116 :
117 4 : const numeric_index_type n = this->size();
118 :
119 88 : for (numeric_index_type i=0; i<n; i++)
120 80 : this->set(i, (*this)(i) * v(i));
121 :
122 : // This is an embarrassingly parallel method, but set() isn't in
123 : // general so set() overzealously marked us as non-closed
124 : #ifndef NDEBUG
125 4 : this->_is_closed = was_closed;
126 : #endif
127 :
128 8 : return *this;
129 : }
130 :
131 :
132 :
133 : template <typename T>
134 8 : NumericVector<T> & LaspackVector<T>::operator /= (const NumericVector<T> & v)
135 : {
136 4 : libmesh_assert_equal_to(size(), v.size());
137 :
138 : #ifndef NDEBUG
139 4 : const bool was_closed = this->_is_closed;
140 : #endif
141 :
142 4 : const numeric_index_type n = this->size();
143 :
144 88 : for (numeric_index_type i=0; i<n; i++)
145 80 : this->set(i, (*this)(i) / v(i));
146 :
147 : // This is an embarrassingly parallel method, but set() isn't in
148 : // general so set() overzealously marked us as non-closed
149 : #ifndef NDEBUG
150 4 : this->_is_closed = was_closed;
151 : #endif
152 :
153 8 : return *this;
154 : }
155 :
156 :
157 :
158 : template <typename T>
159 8 : void LaspackVector<T>::reciprocal()
160 : {
161 4 : const numeric_index_type n = this->size();
162 :
163 : #ifndef NDEBUG
164 4 : const bool was_closed = this->_is_closed;
165 : #endif
166 :
167 88 : for (numeric_index_type i=0; i<n; i++)
168 : {
169 40 : T v = (*this)(i);
170 :
171 : // Don't divide by zero!
172 40 : libmesh_assert_not_equal_to (v, T(0));
173 :
174 80 : this->set(i, 1. / v);
175 : }
176 :
177 : // This is an embarrassingly parallel method, but set() isn't in
178 : // general so set() overzealously marked us as non-closed
179 : #ifndef NDEBUG
180 4 : this->_is_closed = was_closed;
181 : #endif
182 8 : }
183 :
184 :
185 :
186 : template <typename T>
187 0 : void LaspackVector<T>::conjugate()
188 : {
189 0 : const numeric_index_type n = this->size();
190 :
191 : #ifndef NDEBUG
192 0 : const bool was_closed = this->_is_closed;
193 : #endif
194 :
195 0 : for (numeric_index_type i=0; i<n; i++)
196 : {
197 0 : T v = (*this)(i);
198 :
199 0 : this->set(i, libmesh_conj(v) );
200 : }
201 :
202 : // This is an embarrassingly parallel method, but set() isn't in
203 : // general so set() overzealously marked us as non-closed
204 : #ifndef NDEBUG
205 0 : this->_is_closed = was_closed;
206 : #endif
207 0 : }
208 :
209 :
210 : template <typename T>
211 8 : void LaspackVector<T>::add (const T v)
212 : {
213 : #ifndef NDEBUG
214 4 : const bool was_closed = this->_is_closed;
215 : #endif
216 :
217 4 : const numeric_index_type n = this->size();
218 :
219 88 : for (numeric_index_type i=0; i<n; i++)
220 80 : this->add (i, v);
221 :
222 : // This is an embarrassingly parallel method, but set() isn't in
223 : // general so set() overzealously marked us as non-closed
224 : #ifndef NDEBUG
225 4 : this->_is_closed = was_closed;
226 : #endif
227 8 : }
228 :
229 :
230 :
231 :
232 : template <typename T>
233 0 : void LaspackVector<T>::add (const NumericVector<T> & v)
234 : {
235 0 : this->add (1., v);
236 0 : }
237 :
238 :
239 :
240 : template <typename T>
241 32 : void LaspackVector<T>::add (const T a, const NumericVector<T> & v_in)
242 : {
243 : // Make sure the vector passed in is really a LaspackVector
244 16 : const LaspackVector * v = cast_ptr<const LaspackVector *>(&v_in);
245 :
246 : #ifndef NDEBUG
247 16 : const bool was_closed = this->_is_closed;
248 : #endif
249 :
250 16 : libmesh_assert(v);
251 16 : libmesh_assert_equal_to (this->size(), v->size());
252 :
253 352 : for (auto i : make_range(v->size()))
254 320 : this->add (i, a*(*v)(i));
255 :
256 : #ifndef NDEBUG
257 16 : this->_is_closed = was_closed;
258 : #endif
259 32 : }
260 :
261 :
262 :
263 : template <typename T>
264 0 : void LaspackVector<T>::add_vector (const NumericVector<T> & vec_in,
265 : const SparseMatrix<T> & mat_in)
266 : {
267 : // Make sure the data passed in are really in Laspack types
268 0 : const LaspackVector<T> * vec = cast_ptr<const LaspackVector<T> *>(&vec_in);
269 0 : const LaspackMatrix<T> * mat = cast_ptr<const LaspackMatrix<T> *>(&mat_in);
270 :
271 0 : libmesh_assert(vec);
272 0 : libmesh_assert(mat);
273 :
274 : // += mat*vec
275 0 : AddAsgn_VV (&_vec, Mul_QV(const_cast<QMatrix*>(&mat->_QMat),
276 0 : const_cast<QVector*>(&vec->_vec)));
277 0 : }
278 :
279 :
280 : template <typename T>
281 0 : void LaspackVector<T>::add_vector_transpose (const NumericVector<T> &,
282 : const SparseMatrix<T> &)
283 : {
284 0 : libmesh_not_implemented();
285 : }
286 :
287 :
288 :
289 : template <typename T>
290 8 : void LaspackVector<T>::scale (const T factor)
291 : {
292 4 : libmesh_assert (this->initialized());
293 :
294 8 : Asgn_VV(&_vec, Mul_SV (factor, &_vec));
295 8 : }
296 :
297 : template <typename T>
298 0 : void LaspackVector<T>::abs()
299 : {
300 0 : libmesh_assert (this->initialized());
301 :
302 0 : const numeric_index_type n = this->size();
303 :
304 0 : for (numeric_index_type i=0; i!=n; ++i)
305 0 : this->set(i,std::abs((*this)(i)));
306 0 : }
307 :
308 : template <typename T>
309 8 : T LaspackVector<T>::dot (const NumericVector<T> & v_in) const
310 : {
311 4 : libmesh_assert (this->initialized());
312 :
313 : // Make sure the NumericVector passed in is really a LaspackVector
314 4 : const LaspackVector<T> * v = cast_ptr<const LaspackVector<T> *>(&v_in);
315 4 : libmesh_assert(v);
316 :
317 12 : return Mul_VV (const_cast<QVector*>(&(this->_vec)),
318 12 : const_cast<QVector*>(&(v->_vec)));
319 : }
320 :
321 :
322 :
323 : template <typename T>
324 : NumericVector<T> &
325 8 : LaspackVector<T>::operator = (const T s)
326 : {
327 4 : libmesh_assert (this->initialized());
328 4 : libmesh_assert (this->closed());
329 :
330 8 : V_SetAllCmp (&_vec, s);
331 :
332 8 : return *this;
333 : }
334 :
335 :
336 :
337 : template <typename T>
338 : NumericVector<T> &
339 20 : LaspackVector<T>::operator = (const NumericVector<T> & v_in)
340 : {
341 : // Make sure the NumericVector passed in is really a LaspackVector
342 10 : const LaspackVector<T> * v =
343 10 : cast_ptr<const LaspackVector<T> *>(&v_in);
344 :
345 10 : libmesh_assert(v);
346 :
347 20 : *this = *v;
348 :
349 20 : return *this;
350 : }
351 :
352 :
353 :
354 : template <typename T>
355 : LaspackVector<T> &
356 24 : LaspackVector<T>::operator = (const LaspackVector<T> & v)
357 : {
358 12 : libmesh_assert (this->initialized());
359 12 : libmesh_assert (v.closed());
360 12 : libmesh_assert_equal_to (this->size(), v.size());
361 :
362 24 : if (v.size() != 0)
363 24 : Asgn_VV (const_cast<QVector*>(&_vec),
364 12 : const_cast<QVector*>(&v._vec)
365 : );
366 :
367 : #ifndef NDEBUG
368 12 : this->_is_closed = true;
369 : #endif
370 :
371 24 : return *this;
372 : }
373 :
374 :
375 :
376 : template <typename T>
377 : NumericVector<T> &
378 0 : LaspackVector<T>::operator = (const std::vector<T> & v)
379 : {
380 : /**
381 : * Case 1: The vector is the same size of
382 : * The global vector. Only add the local components.
383 : */
384 0 : if (this->size() == v.size())
385 0 : for (auto i : index_range(v))
386 0 : this->set (i, v[i]);
387 :
388 : else
389 0 : libmesh_error_msg("this->size() = " << this->size() << " must be equal to v.size() = " << v.size());
390 :
391 0 : return *this;
392 : }
393 :
394 :
395 : template <typename T>
396 0 : void LaspackVector<T>::localize (NumericVector<T> & v_local_in) const
397 : {
398 : // Make sure the NumericVector passed in is really a LaspackVector
399 0 : LaspackVector<T> * v_local =
400 0 : cast_ptr<LaspackVector<T> *>(&v_local_in);
401 :
402 0 : libmesh_assert(v_local);
403 :
404 0 : *v_local = *this;
405 0 : }
406 :
407 :
408 :
409 : template <typename T>
410 0 : void LaspackVector<T>::localize (NumericVector<T> & v_local_in,
411 : const std::vector<numeric_index_type> & libmesh_dbg_var(send_list)) const
412 : {
413 : // Make sure the NumericVector passed in is really a LaspackVector
414 0 : LaspackVector<T> * v_local =
415 0 : cast_ptr<LaspackVector<T> *>(&v_local_in);
416 :
417 0 : libmesh_assert(v_local);
418 0 : libmesh_assert_less_equal (send_list.size(), v_local->size());
419 :
420 0 : *v_local = *this;
421 0 : }
422 :
423 :
424 :
425 : template <typename T>
426 8 : void LaspackVector<T>::localize (std::vector<T> & v_local,
427 : const std::vector<numeric_index_type> & indices) const
428 : {
429 : // LaspackVectors are serial, so we can just copy values
430 12 : v_local.resize(indices.size());
431 :
432 88 : for (auto i : index_range(v_local))
433 120 : v_local[i] = (*this)(indices[i]);
434 8 : }
435 :
436 :
437 :
438 : template <typename T>
439 0 : void LaspackVector<T>::localize (const numeric_index_type libmesh_dbg_var(first_local_idx),
440 : const numeric_index_type libmesh_dbg_var(last_local_idx),
441 : const std::vector<numeric_index_type> & libmesh_dbg_var(send_list))
442 : {
443 0 : libmesh_assert_equal_to (first_local_idx, 0);
444 0 : libmesh_assert_equal_to (last_local_idx+1, this->size());
445 :
446 0 : libmesh_assert_less_equal (send_list.size(), this->size());
447 :
448 : #ifndef NDEBUG
449 0 : this->_is_closed = true;
450 : #endif
451 0 : }
452 :
453 :
454 :
455 : template <typename T>
456 32 : void LaspackVector<T>::localize (std::vector<T> & v_local) const
457 :
458 : {
459 32 : v_local.resize(this->size());
460 :
461 352 : for (auto i : index_range(v_local))
462 320 : v_local[i] = (*this)(i);
463 32 : }
464 :
465 :
466 :
467 : template <typename T>
468 16 : void LaspackVector<T>::localize_to_one (std::vector<T> & v_local,
469 : const processor_id_type libmesh_dbg_var(pid)) const
470 : {
471 8 : libmesh_assert_equal_to (pid, 0);
472 :
473 16 : this->localize (v_local);
474 16 : }
475 :
476 :
477 :
478 : template <typename T>
479 0 : void LaspackVector<T>::pointwise_mult (const NumericVector<T> & /*vec1*/,
480 : const NumericVector<T> & /*vec2*/)
481 : {
482 0 : libmesh_not_implemented();
483 : }
484 :
485 : template <typename T>
486 0 : void LaspackVector<T>::pointwise_divide (const NumericVector<T> & /*vec1*/,
487 : const NumericVector<T> & /*vec2*/)
488 : {
489 0 : libmesh_not_implemented();
490 : }
491 :
492 : template <typename T>
493 0 : Real LaspackVector<T>::max() const
494 : {
495 0 : libmesh_assert (this->initialized());
496 0 : if (!this->size())
497 0 : return -std::numeric_limits<Real>::max();
498 :
499 0 : Real the_max = libmesh_real((*this)(0));
500 :
501 0 : const numeric_index_type n = this->size();
502 :
503 0 : for (numeric_index_type i=1; i<n; i++)
504 0 : the_max = std::max (the_max, libmesh_real((*this)(i)));
505 :
506 0 : return the_max;
507 : }
508 :
509 :
510 :
511 : template <typename T>
512 0 : Real LaspackVector<T>::min () const
513 : {
514 0 : libmesh_assert (this->initialized());
515 0 : if (!this->size())
516 0 : return std::numeric_limits<Real>::max();
517 :
518 0 : Real the_min = libmesh_real((*this)(0));
519 :
520 0 : const numeric_index_type n = this->size();
521 :
522 0 : for (numeric_index_type i=1; i<n; i++)
523 0 : the_min = std::min (the_min, libmesh_real((*this)(i)));
524 :
525 0 : return the_min;
526 : }
527 :
528 :
529 : //------------------------------------------------------------------
530 : // Explicit instantiations
531 : template class LIBMESH_EXPORT LaspackVector<Number>;
532 :
533 : } // namespace libMesh
534 :
535 :
536 : #endif // #ifdef LIBMESH_HAVE_LASPACK
|
