/home/docs/checkouts/readthedocs.org/user_builds/ginkgo-test/checkouts/latest/build/doc/doxygen/examples/three-pt-stencil-solver.hpp Source File

Reference API

• build
• doc
• doxygen
• examples

three-pt-stencil-solver.hpp

1

508 *

509 * #include <iostream>

510 * #include <map>

511 * #include <string>

512 * #include <vector>

513 *

514 * #include <ginkgo/ginkgo.hpp>

515 *

516 *

517 * template <typename ValueType, typename IndexType>

518 * void generate_stencil_matrix(IndexType discretization_points,

519 * IndexType* row_ptrs, IndexType* col_idxs,

520 * ValueType* values)

521 * {

522 * IndexType pos = 0;

523 * const ValueType coefs[] = {-1, 2, -1};

524 * row_ptrs[0] = pos;

525 * for (IndexType i = 0; i < discretization_points; ++i) {

526 * for (auto ofs : {-1, 0, 1}) {

527 * if (0 <= i + ofs && i + ofs < discretization_points) {

528 * values[pos] = coefs[ofs + 1];

529 * col_idxs[pos] = i + ofs;

530 * ++pos;

531 * }

532 * }

533 * row_ptrs[i + 1] = pos;

534 * }

535 * }

536 *

537 *

538 * template <typename Closure, typename ValueType, typename IndexType>

539 * void generate_rhs(IndexType discretization_points, Closure f, ValueType u0,

540 * ValueType u1, ValueType* rhs)

541 * {

542 * const ValueType h = 1.0 / (discretization_points + 1);

543 * for (IndexType i = 0; i < discretization_points; ++i) {

544 * const ValueType xi = ValueType(i + 1) * h;

545 * rhs[i] = -f(xi) * h * h;

546 * }

547 * rhs[0] += u0;

548 * rhs[discretization_points - 1] += u1;

549 * }

550 *

551 *

552 * template <typename ValueType, typename IndexType>

553 * void print_solution(IndexType discretization_points, ValueType u0, ValueType u1,

554 * const ValueType* u)

555 * {

556 * std::cout << u0 << '\n';

557 * for (IndexType i = 0; i < discretization_points; ++i) {

558 * std::cout << u[i] << '\n';

559 * }

560 * std::cout << u1 << std::endl;

561 * }

562 *

563 *

564 * template <typename Closure, typename ValueType, typename IndexType>

565 * gko::remove_complex<ValueType> calculate_error(IndexType discretization_points,

566 * const ValueType* u,

567 * Closure correct_u)

568 * {

569 * const ValueType h = 1.0 / (discretization_points + 1);

570 * gko::remove_complex<ValueType> error = 0.0;

571 * for (IndexType i = 0; i < discretization_points; ++i) {

572 * using std::abs;

573 * const ValueType xi = ValueType(i + 1) * h;

574 * error += abs(u[i] - correct_u(xi)) / abs(correct_u(xi));

575 * }

576 * return error;

577 * }

578 *

579 * template <typename ValueType, typename IndexType>

580 * void solve_system(const std::string& executor_string,

581 * IndexType discretization_points, IndexType* row_ptrs,

582 * IndexType* col_idxs, ValueType* values, ValueType* rhs,

583 * ValueType* u, gko::remove_complex<ValueType> reduction_factor)

584 * {

585 * using vec = gko::matrix::Dense<ValueType>;

586 * using mtx = gko::matrix::Csr<ValueType, IndexType>;

587 * using cg = gko::solver::Cg<ValueType>;

588 * using bj = gko::preconditioner::Jacobi<ValueType, IndexType>;

589 * using val_array = gko::array<ValueType>;

590 * using idx_array = gko::array<IndexType>;

591 * const auto& dp = discretization_points;

592 *

593 * std::map<std::string, std::function<std::shared_ptr<gko::Executor>()>>

594 * exec_map{

595 * {"omp", [] { return gko::OmpExecutor::create(); }},

596 * {"cuda",

597 * [] {

598 * return gko::CudaExecutor::create(0,

599 * gko::OmpExecutor::create());

600 * }},

601 * {"hip",

602 * [] {

603 * return gko::HipExecutor::create(0, gko::OmpExecutor::create());

604 * }},

605 * {"dpcpp",

606 * [] {

607 * return gko::DpcppExecutor::create(0,

608 * gko::OmpExecutor::create());

609 * }},

610 * {"reference", [] { return gko::ReferenceExecutor::create(); }}};

611 *

612 * const auto exec = exec_map.at(executor_string)(); // throws if not valid

613 * const auto app_exec = exec->get_master();

614 *

615 *

616 * auto matrix = mtx::create(exec, gko::dim<2>(dp),

617 * val_array::view(app_exec, 3 * dp - 2, values),

618 * idx_array::view(app_exec, 3 * dp - 2, col_idxs),

619 * idx_array::view(app_exec, dp + 1, row_ptrs));

620 *

621 * auto b = vec::create(exec, gko::dim<2>(dp, 1),

622 * val_array::view(app_exec, dp, rhs), 1);

623 *

624 * auto x = vec::create(app_exec, gko::dim<2>(dp, 1),

625 * val_array::view(app_exec, dp, u), 1);

626 *

627 * auto solver_gen =

628 * cg::build()

629 * .with_criteria(gko::stop::Iteration::build().with_max_iters(

630 * gko::size_type(dp)),

631 * gko::stop::ResidualNorm<ValueType>::build()

632 * .with_reduction_factor(reduction_factor))

633 * .with_preconditioner(bj::build())

634 * .on(exec);

635 * auto solver = solver_gen->generate(gko::give(matrix));

636 *

637 * solver->apply(b, x);

638 * }

639 *

640 *

641 * int main(int argc, char* argv[])

642 * {

643 * using ValueType = double;

644 * using IndexType = int;

645 *

646 * std::cout << gko::version_info::get() << std::endl;

647 *

648 * if (argc == 2 && std::string(argv[1]) == "--help") {

649 * std::cerr << "Usage: " << argv[0]

650 * << " [executor] [DISCRETIZATION_POINTS]" << std::endl;

651 * std::exit(-1);

652 * }

653 *

654 * const auto executor_string = argc >= 2 ? argv[1] : "reference";

655 * const IndexType discretization_points =

656 * argc >= 3 ? std::atoi(argv[2]) : 100;

657 *

658 * auto correct_u = [](ValueType x) { return x * x * x; };

659 * auto f = [](ValueType x) { return ValueType(6) * x; };

660 * auto u0 = correct_u(0);

661 * auto u1 = correct_u(1);

662 *

663 * std::vector<IndexType> row_ptrs(discretization_points + 1);

664 * std::vector<IndexType> col_idxs(3 * discretization_points - 2);

665 * std::vector<ValueType> values(3 * discretization_points - 2);

666 * std::vector<ValueType> rhs(discretization_points);

667 * std::vector<ValueType> u(discretization_points, 0.0);

668 * const gko::remove_complex<ValueType> reduction_factor = 1e-7;

669 *

670 * generate_stencil_matrix(discretization_points, row_ptrs.data(),

671 * col_idxs.data(), values.data());

672 * generate_rhs(discretization_points, f, u0, u1, rhs.data());

673 *

674 * solve_system(executor_string, discretization_points, row_ptrs.data(),

675 * col_idxs.data(), values.data(), rhs.data(), u.data(),

676 * reduction_factor);

677 *

678 * std::cout << "The average relative error is "

679 * << calculate_error(discretization_points, u.data(), correct_u) /

680 * discretization_points

681 * << std::endl;

682 * }

683 * @endcode

684*/

gko::CudaExecutor::create

static std::shared_ptr< CudaExecutor > create(int device_id, std::shared_ptr< Executor > master, bool device_reset, allocation_mode alloc_mode=default_cuda_alloc_mode, CUstream_st *stream=nullptr)

gko::DpcppExecutor::create

static std::shared_ptr< DpcppExecutor > create(int device_id, std::shared_ptr< Executor > master, std::string device_type="all", dpcpp_queue_property property=dpcpp_queue_property::in_order)

gko::HipExecutor::create

static std::shared_ptr< HipExecutor > create(int device_id, std::shared_ptr< Executor > master, bool device_reset, allocation_mode alloc_mode=default_hip_alloc_mode, CUstream_st *stream=nullptr)

gko::OmpExecutor::create

static std::shared_ptr< OmpExecutor > create(std::shared_ptr< CpuAllocatorBase > alloc=std::make_shared< CpuAllocator >())

Definition executor.hpp:1396

gko::array

Definition array.hpp:166

gko::matrix::Csr

Definition csr.hpp:123

gko::matrix::Dense< ValueType >

gko::matrix::Dense< ValueType >::create

static std::unique_ptr< Dense > create(std::shared_ptr< const Executor > exec, const dim< 2 > &size={}, size_type stride=0)

gko::preconditioner::Jacobi

Definition jacobi.hpp:189

gko::solver::Cg

Definition cg.hpp:49

gko::stop::ResidualNorm

Definition residual_norm.hpp:113

gko::version_info::get

static const version_info & get()

Definition version.hpp:139

gko::log::profile_event_category::solver

@ solver

gko::solver::initial_guess_mode::rhs

@ rhs

gko::abs

constexpr std::enable_if_t<!is_complex_s< T >::value, T > abs(const T &x)

Definition math.hpp:931

gko::size_type

std::size_t size_type

Definition types.hpp:89

gko::give

std::remove_reference< OwningPointer >::type && give(OwningPointer &&p)

Definition utils_helper.hpp:247

gko::remove_complex

typename detail::remove_complex_s< T >::type remove_complex

Definition math.hpp:260

gko::dim

Definition dim.hpp:26

---

Generated by  1.9.8