std::linalg::matrix_rank_1_update

#include <array>
#include <iostream>
#include <linalg>
#include <mdspan>
#include <vector>

template <class Matrix>
void print_mat(const Matrix& A) {
  for(int i = 0; i < A.extent(0); ++i) {
    for(int j = 0; j < A.extent(1) - 1; ++j) {
      std::cout << A[i, j] << ' ';
    }
    std::cout << A[i, A.extent(1) - 1] << '\n';
  }
}

template <class Vector>
void init_vec(Vector& v) {
  for (int i = 0; i < v.extent(0); ++i) {
    v[i] = i;
  }
}

template <class Matrix>
void init_mat(Matrix& A) {
  for(int i = 0; i < A.extent(0); ++i) {
    for(int j = 0; j < A.extent(1); ++j) {
      A[i,j] = A.extent(1) * i + j;
    }
  }
}

int main()
{
  constexpr size_t N = 4;

  std::vector<double> A_vec(N * N);
  std::vector<double> x_vec(N);
  std::array<double, N> y_vec;

  std::mdspan<
    double,
    std::extents<size_t, N, N>> A(A_vec.data());
  std::mdspan x(x_vec.data(), N);
  std::mdspan y(y_vec.data(), N);

  init_mat(A);
  init_vec(x);
  init_vec(y);

  // (1)
  std::cout << "(1)\n";
  std::linalg::matrix_rank_1_update(
    x,
    y,
    A);
  print_mat(A);

  init_mat(A);

  // (2)
  std::cout << "(2)\n";
  std::linalg::matrix_rank_1_update(
    std::execution::par,
    x,
    y,
    A);
  print_mat(A);

  return 0;
}