#include <iostream>
#include <vector>

template <typename T>
std::ostream& operator<<(std::ostream&, const std::vector<std::vector<T>>&);

template <typename T>
std::ostream& operator<<(std::ostream&, const std::vector<T>&);

int main()
{
   
   // create a 2 dimensional int vector with known dimensions
   std::vector<std::vector<int>> vec { { 1, 2, 3, 4 },
                                       { 5, 6, 7, 8 },
                                       { 9, 10, 11, 12 }};

   std::cout << "Let's verify the sizes of the 2D vector....\n";
   std::cout << "The vector has " << vec.size() << " rows.\n";
   std::cout << "The vector has " << vec[0].size() << " columns.\n\n";

   // display the 2D vector
   std::cout << vec << '\n';

   // row element replacement using operator[] indexing
   for (size_t i { }; i < vec[0].size(); ++i)
   {
      size_t fst_row { };
      size_t lst_row { vec.size() - 1 };

      vec[fst_row][i] = 0;
      vec[lst_row][i] = 0;
   }

   std::cout << vec << '\n';

   // I'm sure there is a less messy way to do this.......

   // row element replacement using iterators
   for (auto itr { vec[0].begin() }; itr != vec[0].end(); ++itr)
   {
      *itr = 20;
   }

   for (auto itr { vec[vec.size() - 1].begin() }; itr != vec[vec.size() - 1].end(); ++itr)
   {
      *itr = 20;
   }

   std::cout << vec;
}

template <typename T>
std::ostream& operator<<(std::ostream& os, const std::vector<std::vector<T>>& v)
{
   for (auto const& x : v)
   {
      os << x << '\n';
   }

   return os;
}


template <typename T>
std::ostream& operator<<(std::ostream& os, const std::vector<T>& v)
{
   for (auto const& x : v)
   {
      os << x << ' ';
   }

   return os;
}

#include <memory_resource>
#include <stdexcept>
#include <cstdint>

#include <stdlib.h>
#include <unistd.h>

// Just for example:
void* fftw_malloc(std::size_t n) 
{ 
  static long const page_size = sysconf(_SC_PAGESIZE);
  if (page_size <= -1) 
    return nullptr;
    
  void* result;
  if (::posix_memalign(&result, page_size, n) != 0)
    return nullptr;
    
  return result;
}

void fftw_free(void* p) { return ::free(p); } 

struct bad_fftw_malloc: std::bad_alloc 
{ char const* what() { return "bad fftw_malloc"; } };  

struct fftw_resource: std::pmr::memory_resource
{
  void* do_allocate(std::size_t bytes, std::size_t const a) override
  {
    void* const ptr_orig = fftw_malloc(bytes);
    if (ptr_orig == nullptr) 
      throw bad_fftw_malloc{};
    
    // Ensure the alignment requirement is already satisfied by ptr
    void* ptr_copy = ptr_orig;
    if (ptr_orig != std::align(a, bytes, ptr_copy, bytes)) 
      throw bad_fftw_malloc{}; 
    
    return ptr_orig;
  }
  
  void do_deallocate(void* p, std::size_t, std::size_t) override 
  { 
    if (p) fftw_free(p); 
  } 
  
  bool do_is_equal(std::pmr::memory_resource const&r ) const noexcept override
  { 
    return this == &r; 
  }
};

std::pmr::memory_resource* fftw_malloc_resource() 
{ static fftw_resource res; return &res; } 

#include <iostream>
#include <vector>

int main()
{  
  // Now `my_vector` is suitable for use with fftw3
  std::pmr::vector<double> my_vector(fftw_malloc_resource());  
  
  int constexpr n = 10;
  for (int i = 0; i < n; ++i) my_vector.push_back(i + 0.5);
  for (int i = 0; i < n; ++i) std::cout << my_vector[i] << '\n'; 
}

#include <valarray>
#include <numeric>
#include <iostream>
#include <iomanip>

// Print `a` as a matrix of dimensions `width` by `height` to standard output
// Assuming row-major storage order.
void print_as_matrix(std::valarray<double> a, int width, int height)
{
  for (int y = 0; y < height; ++y)
  {
     std::cout << "[";
     for (int x = 0; x < width; ++x)
     {
       std::cout << std::fixed << a[y * width + x] << (x < width - 1? ", ": ""); 
     }
     std::cout << "]\n";
  }
}

std::slice column_slice(int i, int width, int height) { return{i, height, width}; }
std::slice row_slice(int i, int width, int) { return{width * i, width, 1}; }

int main()
{
  int constexpr width = 8;
  int constexpr height = 4;
  
  // allocate a width-by-height matrix in a single valarray
  // arbitrarily fill it with data
  std::valarray<double> a(width * height);      
  std::iota(std::begin(a), std::end(a), 1.0); 
  std::sin(a /= (width / std::atan(1) * 4.0));
  
  print_as_matrix(a, width, height);
  std::cout << '\n';
  
  // Change column 3 to all zeros
  a[column_slice(3, width, height)] = 0.0;
  // Change the last row to all ones
  a[row_slice(height - 1, width, height)] = 1.0;
  
  print_as_matrix(a, width, height);
  std::cout << '\n';
}