The ranges and accessor example..
Introduction
About the example
The commented program
#include <ginkgo/ginkgo.hpp>
#include <iomanip>
#include <iostream>
LU factorization implementation using Ginkgo ranges For simplicity, we only consider square matrices, and no pivoting.
template <typename Accessor>
note: const means that the range (i.e. the data handler) is constant, not that the underlying data is constant!
{
const auto trail = span{i + 1, A.
length(0)};
note: neither of the lines below need additional memory to store intermediate arrays, all computation is done at the point of assignment
A(trail, i) = A(trail, i) / A(i, i);
caveat: operator * is element-wise multiplication, mmul is matrix multiplication
A(trail, trail) = A(trail, trail) - mmul(A(trail, i), A(i, trail));
}
}
a utility function for printing the factorization on screen
template <typename Accessor>
{
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
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.
use the LU factorization routine defined above to factorize the matrix
print the factorization on screen
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 <ginkgo/ginkgo.hpp>
#include <iomanip>
#include <iostream>
template <typename Accessor>
{
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>
{
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;
ValueType data[] = {
2., 4., 5., -1.0,
4., 11., 12., -1.0,
6., 24., 24., -1.0
};
auto A =
factorize(A);
print_lu(A);
}