MatrixMarket.hpp

Go to the documentation of this file.

 
//
// Copyright (c) 2025, University of Colorado Denver. All rights reserved.
//
// This file is part of <T>LAPACK.
// <T>LAPACK is free software: you can redistribute it and/or modify it under
// the terms of the BSD 3-Clause license. See the accompanying LICENSE file.
 
#ifndef TLAPACK_MATRIXMARKET_HH
#define TLAPACK_MATRIXMARKET_HH
 
#include <random>
#include <tlapack/base/utils.hpp>
#include <type_traits>
 
namespace tlapack {
 
class PCG32 {
   private:
    const uint64_t a = 6364136223846793005ULL;  
    const uint64_t c = 1442695040888963407ULL;  
    uint64_t state;                             
 
   public:
    PCG32(uint64_t s = 1302) noexcept { seed(s); }
 
    void seed(uint64_t s) noexcept
    {
        state = 0;
        operator()();
        state += s;
        operator()();
    }
 
    static constexpr uint32_t min() noexcept { return 0; }
    static constexpr uint32_t max() noexcept { return UINT32_MAX; }
 
    uint32_t operator()() noexcept
    {
        // Constants for a 64-bit state and 32-bit output
        const uint8_t bits = 64;
        const uint8_t opbits = 5;
        constexpr uint8_t mask = (1 << opbits) - 1;
        constexpr uint8_t xshift = (opbits + 32) / 2;
        constexpr uint8_t bottomspare = (bits - 32) - opbits;
 
        // Random rotation
        uint8_t rot = uint8_t(state >> (bits - opbits)) & mask;
 
        // XOR shift high
        uint32_t result =
            static_cast<uint32_t>((state ^ (state >> xshift)) >> bottomspare);
 
        // Rotate right
        result = (result >> rot) | (result << ((-rot) & mask));
 
        // Updates the state
        state = state * a + c;
 
        return result;
    }
};
 
template <class T,
          class Generator,
          class Distribution,
          enable_if_t<is_real<T>, bool> = true>
T rand_helper(Generator& gen, Distribution& d)
{
    return T(d(gen));
}
 
template <class T,
          class Generator,
          class Distribution,
          enable_if_t<is_complex<T>, bool> = true>
T rand_helper(Generator& gen, Distribution& d)
{
    using real_t = real_type<T>;
    real_t r1(d(gen));
    real_t r2(d(gen));
    return complex_type<real_t>(r1, r2);
}
 
template <class T, class Generator>
T rand_helper(Generator& gen)
{
    using real_t = real_type<T>;
    using dist_t =
        typename std::conditional<(is_same_v<real_t, float> ||
                                   is_same_v<real_t, double> ||
                                   is_same_v<real_t, long double>),
                                  std::uniform_real_distribution<real_t>,
                                  std::uniform_real_distribution<float>>::type;
 
    dist_t d;
    return rand_helper<T>(gen, d);
}
 
struct MatrixMarket {
    template <TLAPACK_MATRIX matrix_t>
    void colmajor_read(matrix_t& A, std::istream& is) const
    {
        using idx_t = size_type<matrix_t>;
 
        const idx_t m = nrows(A);
        const idx_t n = ncols(A);
 
        for (idx_t j = 0; j < n; ++j)
            for (idx_t i = 0; i < m; ++i)
                is >> A(i, j);
    }
 
    template <TLAPACK_MATRIX matrix_t>
    void random(matrix_t& A)
    {
        using T = type_t<matrix_t>;
        using idx_t = size_type<matrix_t>;
 
        const idx_t m = nrows(A);
        const idx_t n = ncols(A);
 
        for (idx_t j = 0; j < n; ++j)
            for (idx_t i = 0; i < m; ++i)
                A(i, j) = rand_helper<T>(gen);
    }
 
    template <TLAPACK_UPLO uplo_t, TLAPACK_MATRIX matrix_t>
    void random(uplo_t uplo, matrix_t& A)
    {
        using T = type_t<matrix_t>;
        using idx_t = size_type<matrix_t>;
 
        const idx_t m = nrows(A);
        const idx_t n = ncols(A);
 
        if (uplo == Uplo::Upper) {
            for (idx_t j = 0; j < n; ++j)
                for (idx_t i = 0; i < m; ++i)
                    if (i <= j)
                        A(i, j) = rand_helper<T>(gen);
                    else
                        A(i, j) = T(float(0xCAFEBABE));
        }
        else {
            for (idx_t j = 0; j < n; ++j)
                for (idx_t i = 0; i < m; ++i)
                    if (i >= j)
                        A(i, j) = rand_helper<T>(gen);
                    else
                        A(i, j) = T(float(0xCAFEBABE));
        }
    }
 
    template <TLAPACK_MATRIX matrix_t,
              class Distribution = typename std::conditional<
                  (is_same_v<type_t<matrix_t>, float> ||
                   is_same_v<type_t<matrix_t>, double> ||
                   is_same_v<type_t<matrix_t>, long double>),
                  std::normal_distribution<type_t<matrix_t>>,
                  std::normal_distribution<double>>::type>
    void randn(matrix_t& A)
    {
        using T = type_t<matrix_t>;
        using idx_t = size_type<matrix_t>;
 
        const idx_t m = nrows(A);
        const idx_t n = ncols(A);
 
        Distribution d(0, 1);
        for (idx_t j = 0; j < n; ++j)
            for (idx_t i = 0; i < m; ++i)
                A(i, j) = rand_helper<T>(gen, d);
    }
 
    template <TLAPACK_UPLO uplo_t,
              TLAPACK_MATRIX matrix_t,
              class Distribution = typename std::conditional<
                  (is_same_v<type_t<matrix_t>, float> ||
                   is_same_v<type_t<matrix_t>, double> ||
                   is_same_v<type_t<matrix_t>, long double>),
                  std::normal_distribution<type_t<matrix_t>>,
                  std::normal_distribution<double>>::type>
    void randn(uplo_t uplo, matrix_t& A)
    {
        using T = type_t<matrix_t>;
        using idx_t = size_type<matrix_t>;
 
        const idx_t m = nrows(A);
        const idx_t n = ncols(A);
 
        Distribution d(0, 1);
        if (uplo == Uplo::Upper) {
            for (idx_t j = 0; j < n; ++j)
                for (idx_t i = 0; i < m; ++i)
                    if (i <= j)
                        A(i, j) = rand_helper<T>(gen, d);
                    else
                        A(i, j) = T(float(0xCAFEBABE));
        }
        else {
            for (idx_t j = 0; j < n; ++j)
                for (idx_t i = 0; i < m; ++i)
                    if (i >= j)
                        A(i, j) = rand_helper<T>(gen, d);
                    else
                        A(i, j) = T(float(0xCAFEBABE));
        }
    }
 
    template <TLAPACK_MATRIX matrix_t>
    void hessenberg(matrix_t& A)
    {
        using T = type_t<matrix_t>;
        using idx_t = size_type<matrix_t>;
 
        const idx_t m = nrows(A);
        const idx_t n = ncols(A);
 
        for (idx_t j = 0; j < n; ++j)
            for (idx_t i = 0; i < m; ++i)
                if (i <= j + 1)
                    A(i, j) = rand_helper<T>(gen);
                else
                    A(i, j) = T(float(0xFA57C0DE));
    }
 
    template <TLAPACK_MATRIX matrix_t>
    void hilbert(matrix_t& A) const
    {
        using T = type_t<matrix_t>;
        using idx_t = size_type<matrix_t>;
 
        const idx_t m = nrows(A);
        const idx_t n = ncols(A);
 
        for (idx_t j = 0; j < n; ++j)
            for (idx_t i = 0; i < m; ++i)
                A(i, j) = T(1) / T(i + j + 1);
    }
 
    template <TLAPACK_MATRIX matrix_t>
    void single_value(matrix_t& A, const type_t<matrix_t>& val) const
    {
        using idx_t = size_type<matrix_t>;
 
        const idx_t m = nrows(A);
        const idx_t n = ncols(A);
 
        for (idx_t j = 0; j < n; ++j)
            for (idx_t i = 0; i < m; ++i)
                A(i, j) = val;
    }
 
    template <TLAPACK_UPLO uplo_t, TLAPACK_MATRIX matrix_t>
    void single_value(uplo_t uplo,
                      matrix_t& A,
                      const type_t<matrix_t>& val) const
    {
        using T = type_t<matrix_t>;
        using idx_t = size_type<matrix_t>;
 
        const idx_t m = nrows(A);
        const idx_t n = ncols(A);
 
        if (uplo == Uplo::Upper) {
            for (idx_t j = 0; j < n; ++j)
                for (idx_t i = 0; i < m; ++i)
                    if (i <= j)
                        A(i, j) = val;
                    else
                        A(i, j) = T(float(0xCAFEBABE));
        }
        else {
            for (idx_t j = 0; j < n; ++j)
                for (idx_t i = 0; i < m; ++i)
                    if (i >= j)
                        A(i, j) = val;
                    else
                        A(i, j) = T(float(0xCAFEBABE));
        }
    }
 
    PCG32 gen;
};
 
}  // namespace tlapack
 
#endif  // TLAPACK_MATRIXMARKET_HH