1 // --------------------------------------------------------------------- 2 // 3 // Copyright (C) 2015 - 2020 by the deal.II authors 4 // 5 // This file is part of the deal.II library. 6 // 7 // The deal.II library is free software; you can use it, redistribute 8 // it, and/or modify it under the terms of the GNU Lesser General 9 // Public License as published by the Free Software Foundation; either 10 // version 2.1 of the License, or (at your option) any later version. 11 // The full text of the license can be found in the file LICENSE.md at 12 // the top level directory of deal.II. 13 // 14 // --------------------------------------------------------------------- 15 16 #ifndef dealii_la_vector_templates_h 17 #define dealii_la_vector_templates_h 18 19 #include <deal.II/base/config.h> 20 21 #include <deal.II/lac/la_vector.h> 22 #include <deal.II/lac/vector_operation.h> 23 #include <deal.II/lac/vector_operations_internal.h> 24 25 #include <boost/io/ios_state.hpp> 26 27 #include <iomanip> 28 #include <iostream> 29 30 DEAL_II_NAMESPACE_OPEN 31 32 namespace LinearAlgebra 33 { 34 template <typename Number> 35 void reinit(const size_type size,const bool omit_zeroing_entries)36 Vector<Number>::reinit(const size_type size, const bool omit_zeroing_entries) 37 { 38 ReadWriteVector<Number>::reinit(size, omit_zeroing_entries); 39 } 40 41 42 43 template <typename Number> 44 template <typename Number2> 45 void reinit(const ReadWriteVector<Number2> & in_vector,const bool omit_zeroing_entries)46 Vector<Number>::reinit(const ReadWriteVector<Number2> &in_vector, 47 const bool omit_zeroing_entries) 48 { 49 ReadWriteVector<Number>::reinit(in_vector, omit_zeroing_entries); 50 } 51 52 53 54 template <typename Number> 55 void reinit(const IndexSet & locally_stored_indices,const bool omit_zeroing_entries)56 Vector<Number>::reinit(const IndexSet &locally_stored_indices, 57 const bool omit_zeroing_entries) 58 { 59 ReadWriteVector<Number>::reinit(locally_stored_indices, 60 omit_zeroing_entries); 61 } 62 63 64 65 template <typename Number> 66 void reinit(const VectorSpaceVector<Number> & V,const bool omit_zeroing_entries)67 Vector<Number>::reinit(const VectorSpaceVector<Number> &V, 68 const bool omit_zeroing_entries) 69 { 70 // Check that casting will work. 71 Assert(dynamic_cast<const Vector<Number> *>(&V) != nullptr, 72 ExcVectorTypeNotCompatible()); 73 74 // Downcast V. If fails, throws an exception. 75 const Vector<Number> &down_V = dynamic_cast<const Vector<Number> &>(V); 76 Assert(down_V.size() == this->size(), 77 ExcMessage( 78 "Cannot add two vectors with different numbers of elements")); 79 80 ReadWriteVector<Number>::reinit(down_V, omit_zeroing_entries); 81 } 82 83 84 85 template <typename Number> 86 Vector<Number> & 87 Vector<Number>::operator=(const Vector<Number> &in_vector) 88 { 89 if (PointerComparison::equal(this, &in_vector)) 90 return *this; 91 92 this->thread_loop_partitioner = in_vector.thread_loop_partitioner; 93 if (this->size() != in_vector.size()) 94 this->reinit(in_vector, true); 95 96 dealii::internal::VectorOperations::Vector_copy<Number, Number> copier( 97 in_vector.values.get(), this->values.get()); 98 dealii::internal::VectorOperations::parallel_for( 99 copier, 100 static_cast<size_type>(0), 101 this->size(), 102 this->thread_loop_partitioner); 103 104 return *this; 105 } 106 107 108 109 template <typename Number> 110 template <typename Number2> 111 Vector<Number> & 112 Vector<Number>::operator=(const Vector<Number2> &in_vector) 113 { 114 this->thread_loop_partitioner = in_vector.thread_loop_partitioner; 115 if (this->size() != in_vector.size()) 116 ReadWriteVector<Number>::reinit(in_vector, true); 117 118 dealii::internal::VectorOperations::Vector_copy<Number, Number2> copier( 119 in_vector.values.get(), this->values.get()); 120 dealii::internal::VectorOperations::parallel_for( 121 copier, 122 static_cast<size_type>(0), 123 this->size(), 124 this->thread_loop_partitioner); 125 126 return *this; 127 } 128 129 130 131 template <typename Number> 132 Vector<Number> & 133 Vector<Number>::operator=(const Number s) 134 { 135 Assert(s == static_cast<Number>(0), 136 ExcMessage("Only 0 can be assigned to a vector.")); 137 (void)s; 138 139 dealii::internal::VectorOperations::Vector_set<Number> setter( 140 Number(), this->values.get()); 141 dealii::internal::VectorOperations::parallel_for( 142 setter, 0, this->size(), this->thread_loop_partitioner); 143 144 return *this; 145 } 146 147 148 149 template <typename Number> 150 Vector<Number> & 151 Vector<Number>::operator*=(const Number factor) 152 { 153 AssertIsFinite(factor); 154 155 dealii::internal::VectorOperations::Vectorization_multiply_factor<Number> 156 vector_multiply(this->values.get(), factor); 157 dealii::internal::VectorOperations::parallel_for( 158 vector_multiply, 159 static_cast<size_type>(0), 160 this->size(), 161 this->thread_loop_partitioner); 162 163 return *this; 164 } 165 166 167 168 template <typename Number> 169 Vector<Number> & 170 Vector<Number>::operator/=(const Number factor) 171 { 172 AssertIsFinite(factor); 173 Assert(factor != Number(0.), ExcZero()); 174 this->operator*=(Number(1.) / factor); 175 176 return *this; 177 } 178 179 180 181 template <typename Number> 182 Vector<Number> & 183 Vector<Number>::operator+=(const VectorSpaceVector<Number> &V) 184 { 185 // Check that casting will work. 186 Assert(dynamic_cast<const Vector<Number> *>(&V) != nullptr, 187 ExcVectorTypeNotCompatible()); 188 189 // Downcast V. If fails, throws an exception. 190 const Vector<Number> &down_V = dynamic_cast<const Vector<Number> &>(V); 191 Assert(down_V.size() == this->size(), 192 ExcMessage( 193 "Cannot add two vectors with different numbers of elements")); 194 195 dealii::internal::VectorOperations::Vectorization_add_v<Number> vector_add( 196 this->values.get(), down_V.values.get()); 197 dealii::internal::VectorOperations::parallel_for( 198 vector_add, 0, this->size(), this->thread_loop_partitioner); 199 200 return *this; 201 } 202 203 204 205 template <typename Number> 206 Vector<Number> & 207 Vector<Number>::operator-=(const VectorSpaceVector<Number> &V) 208 { 209 // Check that casting will work. 210 Assert(dynamic_cast<const Vector<Number> *>(&V) != nullptr, 211 ExcVectorTypeNotCompatible()); 212 213 // Downcast V. If fails, throws an exception. 214 const Vector<Number> &down_V = dynamic_cast<const Vector<Number> &>(V); 215 Assert(down_V.size() == this->size(), 216 ExcMessage( 217 "Cannot subtract two vectors with different numbers of elements")); 218 dealii::internal::VectorOperations::Vectorization_subtract_v<Number> 219 vector_subtract(this->values.get(), down_V.values.get()); 220 dealii::internal::VectorOperations::parallel_for( 221 vector_subtract, 0, this->size(), this->thread_loop_partitioner); 222 223 return *this; 224 } 225 226 227 228 template <typename Number> 229 Number Vector<Number>::operator*(const VectorSpaceVector<Number> &V) const 230 { 231 // Check that casting will work. 232 Assert(dynamic_cast<const Vector<Number> *>(&V) != nullptr, 233 ExcVectorTypeNotCompatible()); 234 235 // Downcast V. If fails, throws an exception. 236 const Vector<Number> &down_V = dynamic_cast<const Vector<Number> &>(V); 237 Assert(down_V.size() == this->size(), 238 ExcMessage("Cannot compute the scalar product " 239 "of two vectors with different numbers of elements")); 240 Number sum; 241 dealii::internal::VectorOperations::Dot<Number, Number> dot( 242 this->values.get(), down_V.values.get()); 243 dealii::internal::VectorOperations::parallel_reduce( 244 dot, 0, this->size(), sum, this->thread_loop_partitioner); 245 246 return sum; 247 } 248 249 250 251 template <typename Number> 252 void import(const ReadWriteVector<Number> &,VectorOperation::values,std::shared_ptr<const CommunicationPatternBase>)253 Vector<Number>::import(const ReadWriteVector<Number> &, 254 VectorOperation::values, 255 std::shared_ptr<const CommunicationPatternBase>) 256 { 257 AssertThrow(false, ExcMessage("This function is not implemented.")); 258 } 259 260 261 262 template <typename Number> 263 inline void add(const Number a)264 Vector<Number>::add(const Number a) 265 { 266 AssertIsFinite(a); 267 268 dealii::internal::VectorOperations::Vectorization_add_factor<Number> 269 vector_add(this->values.get(), a); 270 dealii::internal::VectorOperations::parallel_for( 271 vector_add, 0, this->size(), this->thread_loop_partitioner); 272 } 273 274 275 276 template <typename Number> 277 void add(const Number a,const VectorSpaceVector<Number> & V)278 Vector<Number>::add(const Number a, const VectorSpaceVector<Number> &V) 279 { 280 // Check that casting will work. 281 Assert(dynamic_cast<const Vector<Number> *>(&V) != nullptr, 282 ExcVectorTypeNotCompatible()); 283 284 // Downcast V. If fails, throws an exception. 285 const Vector<Number> &down_V = dynamic_cast<const Vector<Number> &>(V); 286 AssertIsFinite(a); 287 Assert(down_V.size() == this->size(), 288 ExcMessage( 289 "Cannot add two vectors with different numbers of elements")); 290 291 dealii::internal::VectorOperations::Vectorization_add_av<Number> 292 vector_add_av(this->values.get(), down_V.values.get(), a); 293 dealii::internal::VectorOperations::parallel_for( 294 vector_add_av, 0, this->size(), this->thread_loop_partitioner); 295 } 296 297 298 299 template <typename Number> 300 void add(const Number a,const VectorSpaceVector<Number> & V,const Number b,const VectorSpaceVector<Number> & W)301 Vector<Number>::add(const Number a, 302 const VectorSpaceVector<Number> &V, 303 const Number b, 304 const VectorSpaceVector<Number> &W) 305 { 306 // Check that casting will work. 307 Assert(dynamic_cast<const Vector<Number> *>(&V) != nullptr, 308 ExcVectorTypeNotCompatible()); 309 // Check that casting will work. 310 Assert(dynamic_cast<const Vector<Number> *>(&W) != nullptr, 311 ExcVectorTypeNotCompatible()); 312 313 // Downcast V. If fails, throws an exception. 314 const Vector<Number> &down_V = dynamic_cast<const Vector<Number> &>(V); 315 // Downcast W. If fails, throws an exception. 316 const Vector<Number> &down_W = dynamic_cast<const Vector<Number> &>(W); 317 AssertIsFinite(a); 318 Assert(down_V.size() == this->size(), 319 ExcMessage( 320 "Cannot add two vectors with different numbers of elements")); 321 AssertIsFinite(b); 322 Assert(down_W.size() == this->size(), 323 ExcMessage( 324 "Cannot add two vectors with different numbers of elements")); 325 326 dealii::internal::VectorOperations::Vectorization_add_avpbw<Number> 327 vector_add( 328 this->values.get(), down_V.values.get(), down_W.values.get(), a, b); 329 dealii::internal::VectorOperations::parallel_for( 330 vector_add, 0, this->size(), this->thread_loop_partitioner); 331 } 332 333 334 335 template <typename Number> 336 void sadd(const Number s,const Number a,const VectorSpaceVector<Number> & V)337 Vector<Number>::sadd(const Number s, 338 const Number a, 339 const VectorSpaceVector<Number> &V) 340 { 341 AssertIsFinite(s); 342 AssertIsFinite(a); 343 344 // Check that casting will work. 345 Assert(dynamic_cast<const Vector<Number> *>(&V) != nullptr, 346 ExcVectorTypeNotCompatible()); 347 348 // Downcast V. It fails, throws an exception. 349 const Vector<Number> &down_V = dynamic_cast<const Vector<Number> &>(V); 350 dealii::internal::VectorOperations::Vectorization_sadd_xav<Number> 351 vector_sadd_xav(this->values.get(), down_V.values.get(), a, s); 352 dealii::internal::VectorOperations::parallel_for( 353 vector_sadd_xav, 0, this->size(), this->thread_loop_partitioner); 354 } 355 356 357 358 template <typename Number> 359 void scale(const VectorSpaceVector<Number> & scaling_factors)360 Vector<Number>::scale(const VectorSpaceVector<Number> &scaling_factors) 361 { 362 // Check that casting will work. 363 Assert(dynamic_cast<const Vector<Number> *>(&scaling_factors) != nullptr, 364 ExcVectorTypeNotCompatible()); 365 366 // Downcast scaling_factors. If fails, throws an exception. 367 const Vector<Number> &down_scaling_factors = 368 dynamic_cast<const Vector<Number> &>(scaling_factors); 369 Assert(down_scaling_factors.size() == this->size(), 370 ExcMessage( 371 "Cannot add two vectors with different numbers of elements")); 372 373 dealii::internal::VectorOperations::Vectorization_scale<Number> 374 vector_scale(this->values.get(), down_scaling_factors.values.get()); 375 dealii::internal::VectorOperations::parallel_for( 376 vector_scale, 0, this->size(), this->thread_loop_partitioner); 377 } 378 379 380 381 template <typename Number> 382 void equ(const Number a,const VectorSpaceVector<Number> & V)383 Vector<Number>::equ(const Number a, const VectorSpaceVector<Number> &V) 384 { 385 AssertIsFinite(a); 386 387 // Check that casting will work. 388 Assert(dynamic_cast<const Vector<Number> *>(&V) != nullptr, 389 ExcVectorTypeNotCompatible()); 390 391 // Downcast V. If fails, throws an exception. 392 const Vector<Number> &down_V = dynamic_cast<const Vector<Number> &>(V); 393 dealii::internal::VectorOperations::Vectorization_equ_au<Number> vector_equ( 394 this->values.get(), down_V.values.get(), a); 395 dealii::internal::VectorOperations::parallel_for( 396 vector_equ, 0, this->size(), this->thread_loop_partitioner); 397 } 398 399 400 401 template <typename Number> 402 bool all_zero()403 Vector<Number>::all_zero() const 404 { 405 Assert(this->size(), ExcEmptyObject()); 406 407 const size_type size = this->size(); 408 for (size_type i = 0; i < size; ++i) 409 if (this->values[i] != Number()) 410 return false; 411 412 return true; 413 } 414 415 416 417 template <typename Number> 418 typename Vector<Number>::value_type mean_value()419 Vector<Number>::mean_value() const 420 { 421 Assert(this->size(), ExcEmptyObject()); 422 423 using real_type = typename VectorSpaceVector<Number>::real_type; 424 value_type sum; 425 dealii::internal::VectorOperations::MeanValue<Number> mean_value( 426 this->values.get()); 427 dealii::internal::VectorOperations::parallel_reduce( 428 mean_value, 0, this->size(), sum, this->thread_loop_partitioner); 429 430 return sum / static_cast<real_type>(this->size()); 431 } 432 433 434 435 template <typename Number> 436 typename VectorSpaceVector<Number>::real_type l1_norm()437 Vector<Number>::l1_norm() const 438 { 439 Assert(this->size(), ExcEmptyObject()); 440 441 using real_type = typename VectorSpaceVector<Number>::real_type; 442 real_type sum; 443 dealii::internal::VectorOperations::Norm1<Number, real_type> norm1( 444 this->values.get()); 445 dealii::internal::VectorOperations::parallel_reduce( 446 norm1, 0, this->size(), sum, this->thread_loop_partitioner); 447 448 return sum; 449 } 450 451 452 453 template <typename Number> 454 typename VectorSpaceVector<Number>::real_type l2_norm()455 Vector<Number>::l2_norm() const 456 { 457 Assert(this->size(), ExcEmptyObject()); 458 459 // if l2_norm()^2 is finite and non-zero, the answer is computed as 460 // std::sqrt(norm_sqr()). If norm_sqrt() is infinite or zero, the l2 norm 461 // might still be finite. In that case, recompute ig (this is a rare case, 462 // so working on the vector twice is uncritical and paid off by the extended 463 // precision) using the BLAS approach with a weight, see e.g. dnrm2.f. 464 using real_type = typename VectorSpaceVector<Number>::real_type; 465 real_type norm_square; 466 dealii::internal::VectorOperations::Norm2<Number, real_type> norm2( 467 this->values.get()); 468 dealii::internal::VectorOperations::parallel_reduce( 469 norm2, 0, this->size(), norm_square, this->thread_loop_partitioner); 470 if (numbers::is_finite(norm_square) && 471 norm_square >= std::numeric_limits<real_type>::min()) 472 return std::sqrt(norm_square); 473 else 474 { 475 real_type scale = 0.; 476 real_type sum = 1.; 477 const size_type size = this->size(); 478 for (size_type i = 0; i < size; ++i) 479 { 480 if (this->values[i] != Number()) 481 { 482 const real_type abs_x = 483 numbers::NumberTraits<Number>::abs(this->values[i]); 484 if (scale < abs_x) 485 { 486 sum = 1. + sum * (scale / abs_x) * (scale / abs_x); 487 scale = abs_x; 488 } 489 else 490 sum += (abs_x / scale) * (abs_x / scale); 491 } 492 } 493 AssertIsFinite(scale * std::sqrt(sum)); 494 return scale * std::sqrt(sum); 495 } 496 } 497 498 499 500 template <typename Number> 501 typename VectorSpaceVector<Number>::real_type linfty_norm()502 Vector<Number>::linfty_norm() const 503 { 504 typename ReadWriteVector<Number>::real_type norm = 0.; 505 const size_type size = this->size(); 506 for (size_type i = 0; i < size; ++i) 507 norm = std::max(std::abs(this->values[i]), norm); 508 509 return norm; 510 } 511 512 513 514 template <typename Number> 515 Number add_and_dot(const Number a,const VectorSpaceVector<Number> & V,const VectorSpaceVector<Number> & W)516 Vector<Number>::add_and_dot(const Number a, 517 const VectorSpaceVector<Number> &V, 518 const VectorSpaceVector<Number> &W) 519 { 520 // Check that casting will work. 521 Assert(dynamic_cast<const Vector<Number> *>(&V) != nullptr, 522 ExcVectorTypeNotCompatible()); 523 // Check that casting will work. 524 Assert(dynamic_cast<const Vector<Number> *>(&W) != nullptr, 525 ExcVectorTypeNotCompatible()); 526 527 // Downcast V. If fails, throws an exception. 528 const Vector<Number> &down_V = dynamic_cast<const Vector<Number> &>(V); 529 // Downcast W. If fails, throws an exception. 530 const Vector<Number> &down_W = dynamic_cast<const Vector<Number> &>(W); 531 AssertIsFinite(a); 532 Assert(down_V.size() == this->size(), 533 ExcMessage( 534 "Cannot add two vectors with different numbers of elements")); 535 Assert(down_W.size() == this->size(), 536 ExcMessage( 537 "Cannot add two vectors with different numbers of elements")); 538 539 Number sum; 540 dealii::internal::VectorOperations::AddAndDot<Number> adder( 541 this->values.get(), down_V.values.get(), down_W.values.get(), a); 542 dealii::internal::VectorOperations::parallel_reduce( 543 adder, 0, this->size(), sum, this->thread_loop_partitioner); 544 AssertIsFinite(sum); 545 546 return sum; 547 } 548 549 550 551 template <typename Number> 552 void print_as_numpy_array(std::ostream & out,const unsigned int precision)553 Vector<Number>::print_as_numpy_array(std::ostream & out, 554 const unsigned int precision) const 555 { 556 AssertThrow(out, ExcIO()); 557 boost::io::ios_flags_saver restore_flags(out); 558 559 out.precision(precision); 560 561 const unsigned int n_elements = this->n_elements(); 562 for (unsigned int i = 0; i < n_elements; ++i) 563 out << this->values[i] << ' '; 564 out << '\n' << std::flush; 565 566 AssertThrow(out, ExcIO()); 567 } 568 569 570 571 template <typename Number> 572 void block_write(std::ostream & out)573 Vector<Number>::block_write(std::ostream &out) const 574 { 575 AssertThrow(out, ExcIO()); 576 577 // Other version of the following 578 // out << size() << std::endl << '['; 579 // Reason: operator<< seems to use some resources that lead to problems in 580 // a multithreaded environment. We convert the size index to 581 // unsigned long long int that is at least 64 bits to be able to output it 582 // on all platforms, since std::uint64_t is not in C. 583 const unsigned long long int sz = this->size(); 584 char buf[16]; 585 586 std::sprintf(buf, "%llu", sz); 587 std::strcat(buf, "\n["); 588 589 out.write(buf, std::strlen(buf)); 590 out.write(reinterpret_cast<const char *>(this->begin()), 591 reinterpret_cast<const char *>(this->end()) - 592 reinterpret_cast<const char *>(this->begin())); 593 594 // out << ']'; 595 const char outro = ']'; 596 out.write(&outro, 1); 597 598 AssertThrow(out, ExcIO()); 599 } 600 601 602 603 template <typename Number> 604 void block_read(std::istream & in)605 Vector<Number>::block_read(std::istream &in) 606 { 607 AssertThrow(in, ExcIO()); 608 609 size_type sz; 610 611 char buf[16]; 612 613 in.getline(buf, 16, '\n'); 614 sz = std::atoi(buf); 615 ReadWriteVector<Number>::reinit(sz, true); 616 617 char c; 618 // in >> c; 619 in.read(&c, 1); 620 AssertThrow(c == '[', ExcIO()); 621 622 in.read(reinterpret_cast<char *>(this->begin()), 623 reinterpret_cast<const char *>(this->end()) - 624 reinterpret_cast<const char *>(this->begin())); 625 626 // in >> c; 627 in.read(&c, 1); 628 AssertThrow(c == ']', ExcIO()); 629 } 630 } // namespace LinearAlgebra 631 632 DEAL_II_NAMESPACE_CLOSE 633 634 #endif 635