The ginkgo-ranges program

The ranges and accessor example..

Table of contents
Introduction About the example The commented program	Results Comments about programming and debugging The plain program

Introduction

About the example

The commented program

a utility function for printing the factorization on screen

template <typename Accessor>
void print_lu(const gko::range<Accessor>& A)
{
    std::cout << std::setprecision(2) << std::fixed;
    std::cout << "L = [";
    for (int i = 0; i < A.length(0); ++i) {
        std::cout << "\n  ";
        for (int j = 0; j < A.length(1); ++j) {
            std::cout << (i > j ? A(i, j) : (i == j) * 1.) << " ";
        }
    }
    std::cout << "\n]\n\nU = [";
    for (int i = 0; i < A.length(0); ++i) {
        std::cout << "\n  ";
        for (int j = 0; j < A.length(1); ++j) {
            std::cout << (i <= j ? A(i, j) : 0.) << " ";
        }
    }
    std::cout << "\n]" << std::endl;
}
 
 
int main(int argc, char* argv[])
{
    using ValueType = double;
    using IndexType = int;

Print version information

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

Create some test data, add some padding just to demonstrate how to use it with ranges. clang-format off

ValueType data[] = {
    2.,  4.,  5., -1.0,
    4., 11., 12., -1.0,
    6., 24., 24., -1.0
};

clang-format on

Create a 3-by-3 range, with a 2D row-major accessor using data as the underlying storage. Set the stride (a.k.a. "LDA") to 4.

auto A =
    gko::range<gko::accessor::row_major<ValueType, 2>>(data, 3u, 3u, 4u);

use the LU factorization routine defined above to factorize the matrix

factorize(A);

print the factorization on screen

    print_lu(A);
}

Results

This is the expected output:

L = [
  1.00 0.00 0.00
  2.00 1.00 0.00
  3.00 4.00 1.00
]
 
U = [
  2.00 4.00 5.00
  0.00 3.00 2.00
  0.00 0.00 1.00
]

Comments about programming and debugging

The plain program

 
#include <iomanip>
#include <iostream>
 
#include <ginkgo/ginkgo.hpp>
 
 
template <typename Accessor>
void factorize(const gko::range<Accessor>& A)
{
    using gko::span;
    assert(A.length(0) == A.length(1));
    for (gko::size_type i = 0; i < A.length(0) - 1; ++i) {
        const auto trail = span{i + 1, A.length(0)};
        A(trail, i) = A(trail, i) / A(i, i);
        A(trail, trail) = A(trail, trail) - mmul(A(trail, i), A(i, trail));
    }
}
 
 
template <typename Accessor>
void print_lu(const gko::range<Accessor>& A)
{
    std::cout << std::setprecision(2) << std::fixed;
    std::cout << "L = [";
    for (int i = 0; i < A.length(0); ++i) {
        std::cout << "\n  ";
        for (int j = 0; j < A.length(1); ++j) {
            std::cout << (i > j ? A(i, j) : (i == j) * 1.) << " ";
        }
    }
    std::cout << "\n]\n\nU = [";
    for (int i = 0; i < A.length(0); ++i) {
        std::cout << "\n  ";
        for (int j = 0; j < A.length(1); ++j) {
            std::cout << (i <= j ? A(i, j) : 0.) << " ";
        }
    }
    std::cout << "\n]" << std::endl;
}
 
 
int main(int argc, char* argv[])
{
    using ValueType = double;
    using IndexType = int;
 
    std::cout << gko::version_info::get() << std::endl;
 
    ValueType data[] = {
        2.,  4.,  5., -1.0,
        4., 11., 12., -1.0,
        6., 24., 24., -1.0
    };
 
    auto A =
        gko::range<gko::accessor::row_major<ValueType, 2>>(data, 3u, 3u, 4u);
 
    factorize(A);
 
    print_lu(A);
}