linear algebra – GSL gsl_linalg_SV_decomp returning U and V reversed

I have the following C code:

void calculate_svd_example() {
  int m = 4;
  int n = 5;
  double M(4)(5) = {{1.0, 0.0, 0.0, 0.0, 2.0},
                     {0.0, 0.0, 3.0, 0.0, 0.0},
                     {0.0, 0.0, 0.0, 0.0, 0.0},
                     {0.0, 2.0, 0.0, 0.0, 0.0}};

  gsl_matrix *mat = gsl_matrix_alloc(m, n);
  for (int i = 0; i < m; i++) {
      for (int j = 0; j < n; j++) {
          double x = M(i)(j);
          gsl_matrix_set(mat, i, j, x);
      }
  }
  printf("M = n");
  pretty_print(mat);
  run_svd(mat);
}

#include <stdio.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <time.h>

#include "../include/run_svd.h"

/*
  gsl_matrix_printf prints a matrix as a column vector.  This function
  prints a matrix in block form.
*/
void pretty_print(const gsl_matrix * M)
{
  // Get the dimension of the matrix.
  int rows = M->size1;
  int cols = M->size2;
  // Now print out the data in a square format.
  for(int i = 0; i < rows; i++){
    for(int j = 0; j < cols; j++){
      printf("%f ", gsl_matrix_get(M, i, j));
    }
    printf("n");
  }
  printf("n");
}

void pretty_print_vector(const gsl_vector * M)
{
  int cols = M->size;
    for(int j = 0; j < cols; j++){
      printf("%f ", gsl_vector_get(M, j));
    }
  printf("n");
}

int run_svd(const gsl_matrix * a) {
  gsl_matrix *aa;
  int m = a->size1;
  int n = a->size2;
  if (m >= n) {
    aa = gsl_matrix_alloc(m, n);
    gsl_matrix_memcpy(aa, a);
  } else {
    aa = gsl_matrix_alloc(n, m);
    // Need to transpose the input
    gsl_matrix_transpose_memcpy(aa, a);
  }

  m = aa->size2;
  gsl_matrix * V = gsl_matrix_alloc(m, m);
  gsl_vector * S = gsl_vector_alloc(m);
  gsl_vector * work = gsl_vector_alloc(m);

  /**
   * On output the matrix A is replaced by U. The diagonal elements of 
   * the singular value matrix S are stored in the vector S. The 
   * singular values are non-negative and form a non-increasing sequence 
   * from S_1 to S_N. The matrix V contains the elements of V in 
   * untransposed form. To form the product U S V^T it is necessary to 
   * take the transpose of V. A workspace of length N is required in 
   * work.
   */
  gsl_linalg_SV_decomp(aa, V, S, work);
  printf("U:n");
  pretty_print(aa);
  printf("S:n");
  pretty_print_vector(S);
  printf("V:n");
  pretty_print(V);

  gsl_matrix_free(V);
  gsl_vector_free(S);
  gsl_vector_free(work);

  return 0;
}

It is giving me the following results:

U:
-0.000000 -0.447214  0.000000 -0.000000 
 0.000000 -0.000000 -1.000000 -0.000000 
-1.000000 -0.000000  0.000000  0.000000 
-0.000000 -0.000000  0.000000  1.000000 
-0.000000 -0.894427  0.000000  0.000000 

S:
3.000000 2.236068 2.000000 0.000000 

V:
-0.000000 -1.000000 -0.000000 0.000000 
-1.000000 -0.000000 -0.000000 0.000000 
-0.000000 -0.000000 -0.000000 1.000000 
-0.000000 -0.000000 -1.000000 0.000000 

Aren’t U and V reversed here? Is this a problem with the GSL code, the GSL documentation, or something I’m doing wrong?

c++ – Implementing GSL synchronized_value. Does not pass tests

Core Guidelines mention a type synchronized_value<T>, which supposedly pairs std::mutex with the internal value. I couldn’t find any implementation both in GSL-Lite and MS GSL libraries.
I implemented my class and want to have it reviewed. There are some problems, for instance, make_synch_value does not compile.

I decided to also delete move constructors of xlock, because it stores a reference and I don’t want to deal with invalid state. Moreover I see no reason for such functionality.
Move constructor of syncrhonized_value is also deleted for the same reason (otherwise it would produce invalid xlock).

I used this class to concurrently call driver’s procedures, but have noticed, that with MinGW device stops responding, however with MSVC everything is alright. So I figured out the implementation for this class is to blame (maybe MinGW optimizes out something idk).

I wrote 3 tests: 2 control with std::lock_guard and std::unique_lock that explicitly block a separate mutex. And 1 with synchronized_value.
I receive strange results. Control tests always pass, but synchronized_value fails 10-100 times out of 10000 concurrent invocations in Debug mode. In Release all seems to be fine. I need to find out the reason (in case of optimization Release is expected to fail, not other way around).

Here is a link to my repository with the library and tests (CTest).

#pragma once

#include <mutex>
#include <type_traits>

template <typename ...>
using void_t = void;

template <typename T, typename Tuple, typename = void_t<>>
struct is_aggregate_constructable_ : std::false_type {};

template <typename T, typename ... Args>
struct is_aggregate_constructable_<T,
        std::tuple<Args...>,
        void_t<decltype(T{Args()...})>> : std::true_type {};

template <typename T, typename ... Args>
struct is_aggregate_constructable :
        is_aggregate_constructable_<T, std::tuple<Args...>>{};

// TODO: add timer?
template<typename T, typename Mutex = std::mutex>
class xlock
{
public:
    using value_type = std::conditional_t<std::is_const<T>::value, const T, T>;
    using mutex_type = Mutex;

    xlock(value_type &refValue, Mutex &mutex)
            : refValue_(refValue),
              uniqueLock_(mutex)
    {}

    xlock(const xlock &other) = delete;
    xlock(xlock &&other) = delete;

    xlock &operator=(const xlock &other) = delete;

    xlock &operator=(xlock &&other) = delete;

    value_type *operator->()
    { return &refValue_; }

    const value_type *operator->() const
    { return &refValue_; }

    value_type &operator*()
    { return refValue_; }

    const value_type &operator*() const
    { return refValue_; }

    operator value_type &()
    { return refValue_; }

    operator const value_type &() const
    { return refValue_; }

private:
    value_type &refValue_;
    std::unique_lock<mutex_type> uniqueLock_;
};

template<typename T, typename Mutex = std::mutex>
class synchronized_value
{
public:
    using value_type = T;
    using mutex_type = Mutex;
    using xlock_t = xlock<value_type, mutex_type>;
    using cxlock_t = xlock<const value_type, mutex_type>;

    template<typename =
    std::enable_if_t<std::is_default_constructible<T>::value>>
    synchronized_value()
    {};

    explicit synchronized_value(T value)
            : value_(std::forward<T>(value))
    {}

    synchronized_value(const synchronized_value&) = delete;
    synchronized_value& operator=(const synchronized_value&) = delete;
    synchronized_value& operator=(synchronized_value&&) = delete;

    xlock_t operator*()
    {
        return {value_, mutex_};
    }

    cxlock_t operator*() const
    {
        return {value_, mutex_};
    }

    xlock_t operator->()
    {
        return {value_, mutex_};
    }

    cxlock_t operator->() const
    {
        return {value_, mutex_};
    }

private:

    mutable mutex_type mutex_;
    value_type value_;

};



// this does not work
template<typename T, typename ... R, typename Mutex = std::mutex,
        typename =
std::enable_if_t<std::is_constructible<T, R...>::value ||
                 is_aggregate_constructable<T, R...>()>>
synchronized_value<T, Mutex> make_synch_value(R &&... args)
{
    return {T(std::forward<R>(args)...)};
}

Testing functions:

// run_check.h

#pragma once

#include <mutex>

struct Dummy
{
    int64_t id;
    bool b;
};

// false if check is successful
bool run_check(Dummy& d, uint64_t assign, size_t intervals)
{
    d.id = assign;

    for (size_t i = 0; i != intervals; ++i)
        ++d.id;

    return bool(d.id - assign - intervals);
}

template <typename Guard>
bool check_mutex(Dummy& d, std::mutex& mutex, uint64_t assign, size_t intervals)
{
    Guard guard{mutex};
    d.id = assign;

    for (size_t i = 0; i != intervals; ++i)
        ++d.id;

    return bool(d.id - assign - intervals);
}

Testing unit:


#include "synchronized_value/synchronized_value.h"

#include "run_check.h"

#include <future>
#include <random>

#include <string>
#include <iostream>

template <typename T1, typename T2>
constexpr bool is_same_v = std::is_same<T1, T2>::value;

int main(int argc, const char **argv)
{
    size_t intervals = 10000;

    if (argc == 2)
    {
        try
        {
            intervals = std::stoll(argv(1));
        }
        catch (const std::invalid_argument &)
        {
            std::cerr << "Invalid argument for intervals: using " << intervals <<
                      " intervals fot test" << std::endl;
        }
    }

    static_assert(is_aggregate_constructable<Dummy, uint64_t, bool>(), "Invalid result");
    static_assert(is_aggregate_constructable<Dummy, int, bool>(), "Invalid result");
    static_assert(is_aggregate_constructable<Dummy, int, int>(), "Invalid result");
    static_assert(!is_aggregate_constructable<Dummy, double, int>(), "Invalid result");

    static_assert(is_same_v<typename xlock<Dummy>::value_type,
            Dummy>, "Invalid result");
    static_assert(is_same_v<typename xlock<const Dummy>::value_type,
            const Dummy>, "Invalid result");

    std::mutex mutex;
    Dummy d{1};
    {
        xlock<Dummy> xlock1{d, mutex};
        xlock1->id;
        Dummy dummy = *xlock1;
        Dummy& refDummy = *xlock1;
        const Dummy& crefDummy = *xlock1;
    }
    {
        xlock<const Dummy> xlock_const{d, mutex};
        xlock_const->id;

        Dummy dummy = *xlock_const;
        // Dummy& refDummy = *xlock_const;
        const Dummy& crefDummy = *xlock_const;
    }
    {
        const xlock<const Dummy> xlock_const{d, mutex};
        xlock_const->id;
    }
    {
        synchronized_value<Dummy> synchValue{d};
        synchValue->id;
        Dummy dummy = *synchValue;
        Dummy& refDummy = *synchValue;
        const Dummy& crefDummy = *synchValue;
    }
    {
        /* this does not compile
        const synchronized_value<Dummy> synchValueConst =
                make_synch_value<Dummy>(4, true);
                */
        const synchronized_value<Dummy> synchValueConst2{d};
        synchValueConst2->id;
        Dummy dummy = *synchValueConst2;
        // Dummy& refDummy = *synchValueConst2;
        const Dummy& crefDummy = *synchValueConst2;
    }

    synchronized_value<Dummy> dummy{{16}};

    bool first_check = run_check(*dummy, 42, 1000);

    std::random_device rd;
    std::mt19937 rng(rd());
    std::uniform_int_distribution<int64_t> dist(-69800, 385697);

    std::vector<std::future<bool>> results;

    for (size_t i = 0; i != 1000; ++i)
    {
        results.push_back(std::async(std::launch::async,
                &run_check,
                std::ref((Dummy&)*dummy),
                dist(rng),
                intervals));
    }

    int result = 0;
    for (size_t i = 0; i != 1000; ++i)
        result += (int)results(i).get();

    return result;
}

Control tests you can find at my repo.

  1. How it can be improved?
  2. I may have missed something. If so, what?
  3. Why does the test fails? (basically, why does it not work)

bash – support in compiling and running gcc; gsl over the terminal

Installed gsl 2.5 on Ubuntu 18.04.

Try to compile and execute a sample_matrix.c Script with

$ gcc -Wall -I / usr / local / include -c sample_matrix.c

runs successfully; continue to display the output,

$ ./a.out

gives an error:

bash: ./a.out: No such file or directory

I'm currently in the folder named gsl in which the sample_matrix.c file is stored, and in that folder I also have the gsl pkg installed

Screenshot of the lists in the folder gsl