1 // --------------------------------------------------------------------- 2 // 3 // Copyright (C) 2014 - 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_packaged_operation_h 17 #define dealii_packaged_operation_h 18 19 #include <deal.II/base/config.h> 20 21 #include <deal.II/base/exceptions.h> 22 23 #include <deal.II/lac/vector_memory.h> 24 25 #include <functional> 26 27 DEAL_II_NAMESPACE_OPEN 28 29 // Forward declarations: 30 #ifndef DOXYGEN 31 template <typename Number> 32 class Vector; 33 template <typename Range, typename Domain, typename Payload> 34 class LinearOperator; 35 template <typename Range = Vector<double>> 36 class PackagedOperation; 37 #endif 38 39 40 /** 41 * A class to store a computation. 42 * 43 * The PackagedOperation class allows lazy evaluation of expressions involving 44 * vectors and linear operators. This is done by storing the computational 45 * expression and only performing the computation when either the object is 46 * implicitly converted to a vector object, or <code>apply</code> (or 47 * <code>apply_add</code>) is invoked by hand. This avoids unnecessary 48 * temporary storage of intermediate results. 49 * 50 * The class essentially consists of <code>std::function</code> objects that 51 * store the knowledge of how to generate the result of a computation and 52 * store it in a vector: 53 * @code 54 * std::function<void(Range &)> apply; 55 * std::function<void(Range &)> apply_add; 56 * @endcode 57 * 58 * Similar to the LinearOperator class it also has knowledge about how to 59 * initialize a vector of the @p Range space: 60 * @code 61 * std::function<void(Range &, bool)> reinit_vector; 62 * @endcode 63 * 64 * As an example consider the addition of multiple vectors 65 * @code 66 * Vector<double> a, b, c, d; 67 * // .. 68 * Vector<double> result = a + b - c + d; 69 * @endcode 70 * or the computation of a residual $b-Ax$: 71 * @code 72 * SparseMatrix<double> A; 73 * Vector<double> b, x; 74 * // .. 75 * const auto op_a = linear_operator(A); 76 * 77 * auto residual = b - op_a * x; 78 * @endcode 79 * The expression <code>residual</code> is of type 80 * <code>PackagedOperation<Vector<double>></code>. It stores 81 * references to <code>A</code>, <code>b</code> and <code>x</code> and defers 82 * the actual computation until <code>apply</code>, or <code>apply_add</code> 83 * are explicitly invoked, 84 * @code 85 * Vector<double> y; 86 * residual.reinit_vector(y); 87 * residual.apply(y); 88 * residual.apply_add(y); 89 * @endcode 90 * or until the @p PackagedOperation object is implicitly converted: 91 * @code 92 * Vector<double> y; 93 * y = residual; 94 * y += residual; 95 * y -= residual; 96 * @endcode 97 * 98 * @note The step-20 tutorial program has a detailed usage example of the 99 * LinearOperator class. 100 * 101 * 102 * @ingroup LAOperators 103 */ 104 template <typename Range> 105 class PackagedOperation 106 { 107 public: 108 /** 109 * Create an empty PackagedOperation object. All <code>std::function</code> 110 * member objects are initialized with default variants that throw an 111 * exception upon invocation. 112 */ PackagedOperation()113 PackagedOperation() 114 { 115 apply = [](Range &) { 116 Assert(false, 117 ExcMessage( 118 "Uninitialized PackagedOperation<Range>::apply called")); 119 }; 120 121 apply_add = [](Range &) { 122 Assert(false, 123 ExcMessage( 124 "Uninitialized PackagedOperation<Range>::apply_add called")); 125 }; 126 127 reinit_vector = [](Range &, bool) { 128 Assert(false, 129 ExcMessage("Uninitialized PackagedOperation<Range>::reinit_vector " 130 "method called")); 131 }; 132 } 133 134 /** 135 * Default copy constructor. 136 */ 137 PackagedOperation(const PackagedOperation<Range> &) = default; 138 139 /** 140 * Constructor that creates a PackagedOperation object from a reference 141 * vector @p u. The PackagedOperation returns @p u. 142 * 143 * The PackagedOperation object that is created stores a reference to @p u. 144 * Thus, the vector must remain a valid reference for the whole lifetime of 145 * the PackagedOperation object. All changes made on @p u after the creation 146 * of the PackagedOperation object are reflected by the operator object. 147 */ PackagedOperation(const Range & u)148 PackagedOperation(const Range &u) 149 { 150 *this = u; 151 } 152 153 /** 154 * Default copy assignment operator. 155 */ 156 PackagedOperation<Range> & 157 operator=(const PackagedOperation<Range> &) = default; 158 159 /** 160 * Copy assignment operator that creates a PackagedOperation object from a 161 * reference vector @p u. The PackagedOperation returns @p u. 162 * 163 * The PackagedOperation object that is created stores a reference to @p u. 164 * Thus, the vector must remain a valid reference for the whole lifetime of 165 * the PackagedOperation object. All changes made on @p u after the creation 166 * of the PackagedOperation object are reflected by the operator object. 167 */ 168 PackagedOperation<Range> & 169 operator=(const Range &u) 170 { 171 apply = [&u](Range &v) { v = u; }; 172 173 apply_add = [&u](Range &v) { v += u; }; 174 175 reinit_vector = [&u](Range &v, bool omit_zeroing_entries) { 176 v.reinit(u, omit_zeroing_entries); 177 }; 178 179 return *this; 180 } 181 182 /** 183 * Convert a PackagedOperation to its result. 184 * 185 * This conversion operator creates a vector of the Range space and calls 186 * <code>apply()</code> on it. 187 */ Range()188 operator Range() const 189 { 190 Range result_vector; 191 192 reinit_vector(result_vector, /*bool omit_zeroing_entries=*/true); 193 apply(result_vector); 194 195 return result_vector; 196 } 197 198 /** 199 * @name In-place vector space operations 200 */ 201 //@{ 202 203 /** 204 * Addition with a PackagedOperation @p second_comp with the same @p Range. 205 */ 206 PackagedOperation<Range> & 207 operator+=(const PackagedOperation<Range> &second_comp) 208 { 209 *this = *this + second_comp; 210 return *this; 211 } 212 213 /** 214 * Subtraction with a PackagedOperation @p second_comp with the same @p 215 * Range. 216 */ 217 PackagedOperation<Range> & 218 operator-=(const PackagedOperation<Range> &second_comp) 219 { 220 *this = *this - second_comp; 221 return *this; 222 } 223 224 /** 225 * Add a constant @p offset (of the @p Range space) to the result of a 226 * PackagedOperation. 227 */ 228 PackagedOperation<Range> & 229 operator+=(const Range &offset) 230 { 231 *this = *this + PackagedOperation<Range>(offset); 232 return *this; 233 } 234 235 /** 236 * Subtract a constant @p offset (of the @p Range space) from the result of 237 * a PackagedOperation. 238 */ 239 PackagedOperation<Range> & 240 operator-=(const Range &offset) 241 { 242 *this = *this - PackagedOperation<Range>(offset); 243 return *this; 244 } 245 246 /** 247 * Scalar multiplication of the PackagedOperation with a @p number. 248 */ 249 PackagedOperation<Range> & 250 operator*=(typename Range::value_type number) 251 { 252 *this = *this * number; 253 return *this; 254 } 255 //@} 256 257 /** 258 * Store the result of the PackagedOperation in a vector v of the @p Range 259 * space. 260 */ 261 std::function<void(Range &v)> apply; 262 263 /** 264 * Add the result of the PackagedOperation to a vector v of the @p Range 265 * space. 266 */ 267 std::function<void(Range &v)> apply_add; 268 269 /** 270 * Initializes a vector v of the Range space to be directly usable as the 271 * destination parameter in an application of apply, or apply_add. Similar 272 * to the reinit functions of the vector classes, the boolean determines 273 * whether a fast initialization is done, i.e., if it is set to false the 274 * content of the vector is set to 0. 275 */ 276 std::function<void(Range &v, bool omit_zeroing_entries)> reinit_vector; 277 }; 278 279 280 /** 281 * @name Vector space operations 282 */ 283 //@{ 284 285 /** 286 * @relatesalso PackagedOperation 287 * 288 * Addition of two PackagedOperation objects @p first_comp and @p second_comp 289 * given by vector space addition of the corresponding results. 290 * 291 * @ingroup LAOperators 292 */ 293 template <typename Range> 294 PackagedOperation<Range> 295 operator+(const PackagedOperation<Range> &first_comp, 296 const PackagedOperation<Range> &second_comp) 297 { 298 PackagedOperation<Range> return_comp; 299 300 return_comp.reinit_vector = first_comp.reinit_vector; 301 302 // ensure to have valid PackagedOperation objects by catching first_comp and 303 // second_comp by value 304 305 return_comp.apply = [first_comp, second_comp](Range &v) { 306 first_comp.apply(v); 307 second_comp.apply_add(v); 308 }; 309 310 return_comp.apply_add = [first_comp, second_comp](Range &v) { 311 first_comp.apply_add(v); 312 second_comp.apply_add(v); 313 }; 314 315 return return_comp; 316 } 317 318 /** 319 * @relatesalso PackagedOperation 320 * 321 * Subtraction of two PackagedOperation objects @p first_comp and @p 322 * second_comp given by vector space addition of the corresponding results. 323 * 324 * @ingroup LAOperators 325 */ 326 template <typename Range> 327 PackagedOperation<Range> 328 operator-(const PackagedOperation<Range> &first_comp, 329 const PackagedOperation<Range> &second_comp) 330 { 331 PackagedOperation<Range> return_comp; 332 333 return_comp.reinit_vector = first_comp.reinit_vector; 334 335 // ensure to have valid PackagedOperation objects by catching first_comp and 336 // second_comp by value 337 338 return_comp.apply = [first_comp, second_comp](Range &v) { 339 second_comp.apply(v); 340 v *= -1.; 341 first_comp.apply_add(v); 342 }; 343 344 return_comp.apply_add = [first_comp, second_comp](Range &v) { 345 first_comp.apply_add(v); 346 v *= -1.; 347 second_comp.apply_add(v); 348 v *= -1.; 349 }; 350 351 return return_comp; 352 } 353 354 /** 355 * @relatesalso PackagedOperation 356 * 357 * Scalar multiplication of a PackagedOperation objects @p comp with a scalar 358 * @p number given by a scaling PackagedOperation result with @p number. 359 * 360 * @ingroup LAOperators 361 */ 362 template <typename Range> 363 PackagedOperation<Range> operator*(const PackagedOperation<Range> &comp, 364 typename Range::value_type number) 365 { 366 PackagedOperation<Range> return_comp; 367 368 return_comp.reinit_vector = comp.reinit_vector; 369 370 // the trivial case: number is zero 371 if (number == 0.) 372 { 373 return_comp.apply = [](Range &v) { v = 0.; }; 374 375 return_comp.apply_add = [](Range &) {}; 376 } 377 else 378 { 379 return_comp.apply = [comp, number](Range &v) { 380 comp.apply(v); 381 v *= number; 382 }; 383 384 return_comp.apply_add = [comp, number](Range &v) { 385 v /= number; 386 comp.apply_add(v); 387 v *= number; 388 }; 389 } 390 391 return return_comp; 392 } 393 394 /** 395 * @relatesalso PackagedOperation 396 * 397 * Scalar multiplication of a PackagedOperation objects @p comp with a scalar 398 * @p number given by a scaling PackagedOperation result with @p number. 399 * 400 * @ingroup LAOperators 401 */ 402 template <typename Range> 403 PackagedOperation<Range> operator*(typename Range::value_type number, 404 const PackagedOperation<Range> &comp) 405 { 406 return comp * number; 407 } 408 409 /** 410 * @relatesalso PackagedOperation 411 * 412 * Add a constant @p offset (of the @p Range space) to the result of a 413 * PackagedOperation. 414 * 415 * @ingroup LAOperators 416 */ 417 template <typename Range> 418 PackagedOperation<Range> 419 operator+(const PackagedOperation<Range> &comp, const Range &offset) 420 { 421 return comp + PackagedOperation<Range>(offset); 422 } 423 424 /** 425 * @relatesalso PackagedOperation 426 * 427 * Add a constant @p offset (of the @p Range space) to the result of a 428 * PackagedOperation. 429 * 430 * @ingroup LAOperators 431 */ 432 template <typename Range> 433 PackagedOperation<Range> 434 operator+(const Range &offset, const PackagedOperation<Range> &comp) 435 { 436 return PackagedOperation<Range>(offset) + comp; 437 } 438 439 /** 440 * @relatesalso PackagedOperation 441 * 442 * Subtract a constant @p offset (of the @p Range space) from the result of a 443 * PackagedOperation. 444 * 445 * @ingroup LAOperators 446 */ 447 template <typename Range> 448 PackagedOperation<Range> 449 operator-(const PackagedOperation<Range> &comp, const Range &offset) 450 { 451 return comp - PackagedOperation<Range>(offset); 452 } 453 454 455 /** 456 * @relatesalso PackagedOperation 457 * 458 * Subtract a computational result from a constant @p offset (of the @p Range 459 * space). The result is a PackagedOperation object that applies this 460 * computation. 461 * 462 * @ingroup LAOperators 463 */ 464 template <typename Range> 465 PackagedOperation<Range> 466 operator-(const Range &offset, const PackagedOperation<Range> &comp) 467 { 468 return PackagedOperation<Range>(offset) - comp; 469 } 470 471 //@} 472 473 474 /** 475 * @name Creation of a PackagedOperation object 476 */ 477 //@{ 478 479 namespace internal 480 { 481 namespace PackagedOperationImplementation 482 { 483 // Poor man's trait class that determines whether type T is a vector: 484 // FIXME: Implement this as a proper type trait - similar to 485 // isBlockVector 486 487 template <typename T> 488 class has_vector_interface 489 { 490 template <typename C> 491 static std::false_type 492 test(...); 493 494 template <typename C> 495 static std::true_type 496 test(decltype(&C::operator+=), 497 decltype(&C::operator-=), 498 decltype(&C::l2_norm)); 499 500 public: 501 // type is std::true_type if Matrix provides vmult_add and Tvmult_add, 502 // otherwise it is std::false_type 503 504 using type = decltype(test<T>(nullptr, nullptr, nullptr)); 505 }; // namespace 506 } // namespace PackagedOperationImplementation 507 } // namespace internal 508 509 510 /** 511 * @relatesalso PackagedOperation 512 * 513 * Create a PackagedOperation object that stores the addition of two vectors. 514 * 515 * The PackagedOperation object that is created stores a reference to @p u and 516 * @p v. Thus, the vectors must remain valid references for the whole lifetime 517 * of the PackagedOperation object. All changes made on @p u or @p v after the 518 * creation of the PackagedOperation object are reflected by the operator 519 * object. 520 * 521 * @ingroup LAOperators 522 */ 523 524 template <typename Range, 525 typename = typename std::enable_if< 526 internal::PackagedOperationImplementation::has_vector_interface< 527 Range>::type::value>::type> 528 PackagedOperation<Range> 529 operator+(const Range &u, const Range &v) 530 { 531 PackagedOperation<Range> return_comp; 532 533 // ensure to have valid PackagedOperation objects by catching op by value 534 // u is caught by reference 535 536 return_comp.reinit_vector = [&u](Range &x, bool omit_zeroing_entries) { 537 x.reinit(u, omit_zeroing_entries); 538 }; 539 540 return_comp.apply = [&u, &v](Range &x) { 541 x = u; 542 x += v; 543 }; 544 545 return_comp.apply_add = [&u, &v](Range &x) { 546 x += u; 547 x += v; 548 }; 549 550 return return_comp; 551 } 552 553 554 /** 555 * @relatesalso PackagedOperation 556 * 557 * Create a PackagedOperation object that stores the subtraction of two 558 * vectors. 559 * 560 * The PackagedOperation object that is created stores a reference to @p u and 561 * @p v. Thus, the vectors must remain valid references for the whole lifetime 562 * of the PackagedOperation object. All changes made on @p u or @p v after the 563 * creation of the PackagedOperation object are reflected by the operator 564 * object. 565 * 566 * @ingroup LAOperators 567 */ 568 569 template <typename Range, 570 typename = typename std::enable_if< 571 internal::PackagedOperationImplementation::has_vector_interface< 572 Range>::type::value>::type> 573 PackagedOperation<Range> 574 operator-(const Range &u, const Range &v) 575 { 576 PackagedOperation<Range> return_comp; 577 578 // ensure to have valid PackagedOperation objects by catching op by value 579 // u is caught by reference 580 581 return_comp.reinit_vector = [&u](Range &x, bool omit_zeroing_entries) { 582 x.reinit(u, omit_zeroing_entries); 583 }; 584 585 return_comp.apply = [&u, &v](Range &x) { 586 x = u; 587 x -= v; 588 }; 589 590 return_comp.apply_add = [&u, &v](Range &x) { 591 x += u; 592 x -= v; 593 }; 594 595 return return_comp; 596 } 597 598 599 /** 600 * @relatesalso PackagedOperation 601 * 602 * Create a PackagedOperation object that stores the scaling of a vector with 603 * a @p number. 604 * 605 * The PackagedOperation object that is created stores a reference to @p u. 606 * Thus, the vectors must remain valid references for the whole lifetime of 607 * the PackagedOperation object. All changes made on @p u or @p v after the 608 * creation of the PackagedOperation object are reflected by the operator 609 * object. 610 * 611 * @ingroup LAOperators 612 */ 613 template <typename Range, 614 typename = typename std::enable_if< 615 internal::PackagedOperationImplementation::has_vector_interface< 616 Range>::type::value>::type> 617 PackagedOperation<Range> operator*(const Range & u, 618 typename Range::value_type number) 619 { 620 return PackagedOperation<Range>(u) * number; 621 } 622 623 624 /** 625 * @relatesalso PackagedOperation 626 * 627 * Create a PackagedOperation object that stores the scaling of a vector with 628 * a @p number. 629 * 630 * The PackagedOperation object that is created stores a reference to @p u. 631 * Thus, the vectors must remain valid references for the whole lifetime of 632 * the PackagedOperation object. All changes made on @p u or @p v after the 633 * creation of the PackagedOperation object are reflected by the operator 634 * object. 635 * 636 * @ingroup LAOperators 637 */ 638 template <typename Range, 639 typename = typename std::enable_if< 640 internal::PackagedOperationImplementation::has_vector_interface< 641 Range>::type::value>::type> 642 PackagedOperation<Range> operator*(typename Range::value_type number, 643 const Range & u) 644 { 645 return number * PackagedOperation<Range>(u); 646 } 647 648 649 /** 650 * @relatesalso PackagedOperation 651 * 652 * Create a PackagedOperation object from a LinearOperator and a reference to 653 * a vector @p u of the Domain space. The object stores the PackagedOperation 654 * $\text{op} \,u$ (in matrix notation). <code>return</code> 655 * (<code>return_add</code>) are implemented with <code>vmult(__1,u)</code> 656 * (<code>vmult_add(__1,u)</code>). 657 * 658 * The PackagedOperation object that is created stores a reference to @p u. 659 * Thus, the vector must remain a valid reference for the whole lifetime of 660 * the PackagedOperation object. All changes made on @p u after the creation 661 * of the PackagedOperation object are reflected by the operator object. 662 * 663 * @ingroup LAOperators 664 */ 665 template <typename Range, typename Domain, typename Payload> 666 PackagedOperation<Range> 667 operator*(const LinearOperator<Range, Domain, Payload> &op, const Domain &u) 668 { 669 PackagedOperation<Range> return_comp; 670 671 return_comp.reinit_vector = op.reinit_range_vector; 672 673 // ensure to have valid PackagedOperation objects by catching op by value 674 // u is caught by reference 675 676 return_comp.apply = [op, &u](Range &v) { op.vmult(v, u); }; 677 678 return_comp.apply_add = [op, &u](Range &v) { op.vmult_add(v, u); }; 679 680 return return_comp; 681 } 682 683 684 /** 685 * @relatesalso PackagedOperation 686 * 687 * Create a PackagedOperation object from a LinearOperator and a reference to 688 * a vector @p u of the Range space. The object stores the PackagedOperation 689 * $\text{op}^T \,u$ (in matrix notation). <code>return</code> 690 * (<code>return_add</code>) are implemented with <code>Tvmult(__1,u)</code> 691 * (<code>Tvmult_add(__1,u)</code>). 692 * 693 * The PackagedOperation object that is created stores a reference to @p u. 694 * Thus, the vector must remain a valid reference for the whole lifetime of 695 * the PackagedOperation object. All changes made on @p u after the creation 696 * of the PackagedOperation object are reflected by the operator object. 697 * 698 * @ingroup LAOperators 699 */ 700 template <typename Range, typename Domain, typename Payload> 701 PackagedOperation<Domain> 702 operator*(const Range &u, const LinearOperator<Range, Domain, Payload> &op) 703 { 704 PackagedOperation<Range> return_comp; 705 706 return_comp.reinit_vector = op.reinit_domain_vector; 707 708 // ensure to have valid PackagedOperation objects by catching op by value 709 // u is caught by reference 710 711 return_comp.apply = [op, &u](Domain &v) { op.Tvmult(v, u); }; 712 713 return_comp.apply_add = [op, &u](Domain &v) { op.Tvmult_add(v, u); }; 714 715 return return_comp; 716 } 717 718 719 /** 720 * @relatesalso PackagedOperation 721 * 722 * Composition of a PackagedOperation object with a LinearOperator. The object 723 * stores the computation $\text{op} \,comp$ (in matrix notation). 724 * 725 * @ingroup LAOperators 726 */ 727 template <typename Range, typename Domain, typename Payload> 728 PackagedOperation<Range> 729 operator*(const LinearOperator<Range, Domain, Payload> &op, 730 const PackagedOperation<Domain> & comp) 731 { 732 PackagedOperation<Range> return_comp; 733 734 return_comp.reinit_vector = op.reinit_range_vector; 735 736 // ensure to have valid PackagedOperation objects by catching op by value 737 // u is caught by reference 738 739 return_comp.apply = [op, comp](Domain &v) { 740 GrowingVectorMemory<Range> vector_memory; 741 742 typename VectorMemory<Range>::Pointer i(vector_memory); 743 op.reinit_domain_vector(*i, /*bool omit_zeroing_entries =*/true); 744 745 comp.apply(*i); 746 op.vmult(v, *i); 747 }; 748 749 return_comp.apply_add = [op, comp](Domain &v) { 750 GrowingVectorMemory<Range> vector_memory; 751 752 typename VectorMemory<Range>::Pointer i(vector_memory); 753 op.reinit_range_vector(*i, /*bool omit_zeroing_entries =*/true); 754 755 comp.apply(*i); 756 op.vmult_add(v, *i); 757 }; 758 759 return return_comp; 760 } 761 762 763 /** 764 * @relatesalso PackagedOperation 765 * 766 * Composition of a PackagedOperation object with a LinearOperator. The object 767 * stores the computation $\text{op}^T \,comp$ (in matrix notation). 768 * 769 * @ingroup LAOperators 770 */ 771 template <typename Range, typename Domain, typename Payload> 772 PackagedOperation<Domain> 773 operator*(const PackagedOperation<Range> & comp, 774 const LinearOperator<Range, Domain, Payload> &op) 775 { 776 PackagedOperation<Range> return_comp; 777 778 return_comp.reinit_vector = op.reinit_domain_vector; 779 780 // ensure to have valid PackagedOperation objects by catching op by value 781 // u is caught by reference 782 783 return_comp.apply = [op, comp](Domain &v) { 784 GrowingVectorMemory<Range> vector_memory; 785 786 typename VectorMemory<Range>::Pointer i(vector_memory); 787 op.reinit_range_vector(*i, /*bool omit_zeroing_entries =*/true); 788 789 comp.apply(*i); 790 op.Tvmult(v, *i); 791 }; 792 793 return_comp.apply_add = [op, comp](Domain &v) { 794 GrowingVectorMemory<Range> vector_memory; 795 796 typename VectorMemory<Range>::Pointer i(vector_memory); 797 op.reinit_range_vector(*i, /*bool omit_zeroing_entries =*/true); 798 799 comp.apply(*i); 800 op.Tvmult_add(v, *i); 801 }; 802 803 return return_comp; 804 } 805 806 //@} 807 808 DEAL_II_NAMESPACE_CLOSE 809 810 #endif 811