Commit c09e8c54 authored by MIJIEUX Thomas's avatar MIJIEUX Thomas

Isolate 'exact breakdown' case in test file

parent d43099b4
# -*- mode: org -*-
# -*- coding: utf-8 -*-
* RESULTS
** Run Examples
#+begin_src bash :session foo :results output :exports both
cd fabulous/build/src/test_core/
./testMatrixMarket 100 500
# most matrix market with all ortho scheme + variable max_space (Krylov search space) parameter:
./testOrthoSchemes
# bcsstk14 with IB CGS_RUHE:
./test_RUHE_MGS_bcsstk_EXACT_BREAKDOWN
#+end_src
** Plot the results
#+begin_src R :results output graphics :file (org-babel-temp-file "figure" ".png") :exports both :width 600 :height 400 :session *R*
args <- commandArgs(TRUE)
if (length(args) < 1) {
stop("need filename as parameter")
}
library(ggplot2)
library(reshape2)
library(tools)
data <- read.table(args[1], header=T)
data$ID <- basename(file_path_sans_ext(args[1]))
if (length(args) >= 2) {
for (i in 2:length(args)) {
d <- read.table(args[i], header=T)
d$ID <- basename(file_path_sans_ext(args[i]))
data <- rbind(data, d)
}
}
id = data$ID;
data$ID1 <- paste(id, "minRes")
data$ID2 <- paste(id, "maxRes")
data$ID3 <- paste(id, "minRealRes")
data$ID4 <- paste(id, "maxRealRes")
p <- ggplot(data, aes(nb_mvp)) +
geom_line(aes(y=minRes, color=ID1)) +
geom_line(aes(y=maxRes, color=ID2)) +
geom_line(aes(y=minRealRes, color=ID3)) +
geom_line(aes(y=maxRealRes, color=ID4))
p + scale_y_log10()
#+end_src
......@@ -143,10 +143,18 @@ private:
Block R{_dim, _nbRHS};
compute_real_residual(A, _solution, B, R); // R <- A*sol - B
auto MinMaxConv = R.check_precision(epsilon, A);
#ifdef FABULOUS_DEBUG_MODE
auto min = std::get<0>(MinMaxConv);
auto max = std::get<1>(MinMaxConv);
printf("real residual=(%g;%g)\n", min, max);
#endif
if (std::get<2>(MinMaxConv)) {
_solution_computed = true;
_cold_restart = false;
return true;
}
_cold_restart = true;
return false;
}
......@@ -272,6 +280,7 @@ private:
}
_block_size_sum.clear();
_block_size_sum.push_back(_nb_direction_kept);
_size_krylov_space = _base.get_nb_vect();
return Lambda0.check_precision(epsilon);
}
......@@ -436,13 +445,16 @@ public:
if (_solution_computed)
X.copy(_solution);
else
else if (_iteration_count)
compute_solution(A, X);
if (_cold_restart)
FABULOUS_WARNING("WOULD HAVE COLD RESTART");
#ifdef FABULOUS_DEBUG_MODE
if (convergence) {
auto mm = eval_current_solution(A, B, X);
printf("real residual=(%g;%g)\n", mm.first, mm.second);
printf("real residual=(%e;%e)\n", mm.first, mm.second);
}
#endif
......
......@@ -46,8 +46,13 @@ public:
auto norm = Q.get_norm(j, A);
R.at(j, j) = norm;
if (norm == 0.0) {
FABULOUS_FATAL_ERROR(
"Rank loss in block candidate for extending the Krylov Space" );
#ifndef FABULOUS_IB_EXACT_BREAKDOWN_CONTINUE
FABULOUS_FATAL_ERROR("Rank loss in block candidate for extending the Krylov Space");
#else
FABULOUS_WARNING("Rank loss in block candidate for extending the Krylov Space");
FABULOUS_WARNING("Continuing since IB can address rank loss problems");
norm = 1.0;
#endif
}
LapackKernI::scal(M, S{1.0} / norm, Q.get_vect(j), 1); // Qj = (1/n)Qj
}
......
......@@ -15,7 +15,10 @@ namespace fabulous {
class OrthogonalizerBlockCOMMON : public Orthogonalizer
{
protected:
using Orthogonalizer::Orthogonalizer; // pull parent constructor
OrthogonalizerBlockCOMMON(const Orthogonalizer &ortho):
Orthogonalizer{ortho}
{
}
/**
* \brief Orthogonalisation process :: Classical Gram Schmidt
......
......@@ -19,7 +19,7 @@ enum class OrthoScheme {
MGS = 0, /*!< Modified Gram Schmidt */
CGS = 1, /*!< Classical Gram Schmidt */
IMGS = 2, /*!< Iterated Modified Gram Schmidt */
ICGS = 3 /*!< Iterated Classical Gram Schmidt */
ICGS = 3, /*!< Iterated Classical Gram Schmidt */
};
/* **************** OPERATOR<< OVERLOAD ********************* */
......
......@@ -60,7 +60,8 @@ private:
S{-1.0}, S{1.0}
);
}
{ // Ortho against P
if (P.get_nb_col() > 0) {
// Ortho against P
A.DotProduct( // C_k = P^{H}*W_k
P.get_nb_col(), 1, P.get_ptr(), ldp, W_k, ldw, C_k, ldc);
LapackKernI::gemm( // W_k = W_k - P * C
......@@ -71,7 +72,8 @@ private:
S{-1.0}, S{1.0}
);
}
{// Ortho against directions in Wj that have already been
if (k > 0) {
// Ortho against directions in Wj that have already been
// ortho-gonalized/normalized : D_k without diagonal =
A.DotProduct( // W_i{i<k} * Wk
k, 1, W.get_ptr(), ldw, W_k, ldw, D_k, ldd);
......@@ -83,14 +85,18 @@ private:
S{-1.0}, S{1.0}
);
}
//Hess[...] = norm(W_k)
auto norm = W.get_norm(k, A);
D_k[k] = norm;
if (norm == 0.0) {
FABULOUS_FATAL_ERROR(
"Rank loss in block candidate for extending the Krylov Space" );
#ifndef FABULOUS_IB_EXACT_BREAKDOWN_CONTINUE
FABULOUS_FATAL_ERROR("Rank loss in block candidate for extending the Krylov Space");
#else
FABULOUS_WARNING("Rank loss in block candidate for extending the Krylov Space");
FABULOUS_WARNING("Continuing since IB can address rank loss problems");
norm = 1.0;
#endif
}
D.at(k, k) = norm;
LapackKernI::scal(dim, S{1.0} / norm, W_k, 1); // W[k] /= norm(W[k])
}
}
......@@ -132,7 +138,8 @@ private:
for (int i = 0 ; i<nb_vect_in_base ; ++i)
h_k[i] += tmp.at(i, 0); // h_k += temp;
}
{// Compute ortho with P
if (P.get_nb_col() != 0 ) {
// Compute ortho with P
Block<S> tmp{P.get_nb_col(), 1};
int ld = tmp.get_leading_dim();
A.DotProduct( // C_k = P^{H}*W_k
......@@ -147,9 +154,10 @@ private:
for (int i = 0; i < P.get_nb_col(); ++i)
C_k[i] += tmp.at(i, 0); // C_k += tmp
}
{//Ortho against directions in Wj that have already been
if (k > 0) {
//Ortho against directions in Wj that have already been
//ortho-gonalized/normalized :
Block<S> tmp{k,1};
Block<S> tmp{k, 1};
int ld = tmp.get_leading_dim();
// D_k without diagonal = W_i{i<k} * Wk
A.DotProduct(k, 1, W.get_ptr(), ldw, W_k, ldw, tmp.get_ptr(), ld);
......@@ -166,11 +174,16 @@ private:
}
//Hess[...] = norm(W_k)
auto norm = W.get_norm(k, A);
D_k[k] = norm;
if (norm == 0.0) {
FABULOUS_FATAL_ERROR(
"Rank loss in block candidate for extending the Krylov Space" );
#ifndef FABULOUS_IB_EXACT_BREAKDOWN_CONTINUE
FABULOUS_FATAL_ERROR("Rank loss in block candidate for extending the Krylov Space");
#else
FABULOUS_WARNING("Rank loss in block candidate for extending the Krylov Space");
FABULOUS_WARNING("Continuing since IB can address rank loss problems");
norm = 1.0;
#endif
}
D.at(k, k) = norm;
LapackKernI::scal(dim, S{1.0} / norm, W_k, 1); // W[k] /= norm(W[k])
}
}
......@@ -221,11 +234,16 @@ private:
}
//Normalization of W_k and storage inside L
auto norm = W.get_norm(k, A);
D_k[k] = norm;
if (norm == 0.0) {
FABULOUS_FATAL_ERROR(
"Rank loss in block candidate for extending the Krylov Space");
#ifndef FABULOUS_IB_EXACT_BREAKDOWN_CONTINUE
FABULOUS_FATAL_ERROR("Rank loss in block candidate for extending the Krylov Space");
#else
FABULOUS_WARNING("Rank loss in block candidate for extending the Krylov Space");
FABULOUS_WARNING("Continuing since IB can address rank loss problems");
norm = 1.0;
#endif
}
D_k[k] = norm;
LapackKernI::scal(dim, S{1.0} / norm, W_k, 1); // W[k] /= norm(W[k])
}
}
......@@ -274,7 +292,6 @@ private:
C_k[n] += tmp;
}
}
//Ortho against already ortho-normalized vectors in W
for (int n = 0; n < k; ++n) {
for (int t = 0; t < _nb_iteration; ++t) {
......@@ -288,17 +305,20 @@ private:
}
// Normalization of Wk and storage inside L
auto norm = W.get_norm(k, A);
D_k[k] = norm;
if (norm == 0.0) {
//We capture the case of current W_k \in span{Vm,P,W_i_{i<k}}
FABULOUS_FATAL_ERROR("Rank loss in block "
"candidate for extending the Krylov Space");
#ifndef FABULOUS_IB_EXACT_BREAKDOWN_CONTINUE
FABULOUS_FATAL_ERROR("Rank loss in block candidate for extending the Krylov Space");
#else
FABULOUS_WARNING("Rank loss in block candidate for extending the Krylov Space");
FABULOUS_WARNING("Continuing since IB can address rank loss problems");
norm = 1.0;
#endif
}
D_k[k] = norm;
LapackKernI::scal(dim, S{1.0} / norm, W_k, 1); // W[k] /= norm(W[k])
}
}
template<class Matrix>
void run(HessIB<S> &L, Base<S> &base, BlockWP<S> &WP, Matrix& A)
{
......
#ifndef FABULOUS_ALGORITHM_H
#define FABULOUS_ALGORITHM_H
#include <cassert>
#include "fabulous/utils/Utils.hpp"
namespace fabulous {
/**
* \brief methods that apply on Blocks (QR factorization)
*/
struct Algorithm {
public:
/**
* \brief compute QR factorization of user-distribution Block
*
* This factorization use the Gram-Schmidt method
*
* \param[in,out] Q the orthogonal part of the QR factorization
* \param[out] R the R part of the QR factorization
* \param[in] A User Callback Object with a DotProduct method implementation
*
* \note R MUST BE already allocated ( Q.get_nb_col() x Q.get_nb_col() )
*/
template< class Matrix, class Block, class S = typename Block::value_type >
static void InPlaceQRFactoMGS_User(Block &Q, Block &R, const Matrix &A)
{
int M = Q.get_nb_row();
int N = Q.get_nb_col();
int LD = Q.get_leading_dim();
assert( R.get_nb_col() == Q.get_nb_col() );
assert( R.get_nb_col() == R.get_nb_row() );
for (int j = 0; j < N; ++j) {
for (int i = 0; i < j; ++i) {
S dot; // dot = DOT(Q_i, Q_j)
A.DotProduct(1, 1, Q.get_vect(i), LD, Q.get_vect(j), LD, &dot, 1);
// Qj = Qj - dot(Qj,Qi)*Qi
LapackKernI::axpy(M, -dot, Q.get_vect(i), 1, Q.get_vect(j), 1);
R.at(i, j) = dot;
}
auto norm = Q.get_norm(j, A);
R.at(j, j) = norm;
if (norm == 0.0) {
FABULOUS_FATAL_ERROR(
"Rank loss in block candidate for extending the Krylov Space" );
}
LapackKernI::scal(M, S{1.0} / norm, Q.get_vect(j), 1); // Qj = (1/n)Qj
}
}
/**
* \brief Compute QR factorization
*
* Compute QR factorization using Lapack kernels.
* The first parameter Q is the input matrix and will be overwritten with the
* Q factor of the factorization.
*
* \param[in,out] Q when the method returns: the Q part of the factorization
* \param[out] R the R part of the factorization
* it is assumed either the user allocated a 'R' Block such as:
* R.get_nb_col() == this->get_nb_col() and
* R.get_nb_row() == this->get_nb_col()
* or the user passed a block of size (0x0) and R will no be computed
*/
template< class Block, class T = typename Block::value_type >
static void InPlaceQRFacto(Block &Q, Block &R)
{
int M = Q.get_nb_row();
int N = Q.get_nb_col();
int LD = Q.get_leading_dim();
if ( R.get_nb_col() != 0 ) {
assert( R.get_nb_col() == Q.get_nb_col() );
assert( R.get_nb_row() == Q.get_nb_col() );
}
std::vector<T> tau;
LapackKernI::geqrf(M, N, Q.get_ptr(), LD, tau);
if (R.get_nb_col() != 0) {
LapackKernI::lacpy( 'U', M, N, Q.get_ptr(), LD,
R.get_ptr(), R.get_leading_dim());
}
LapackKernI::orgqr(M, N, N, Q.get_ptr(), LD, tau.data());
}
};
}; // namespace fabulous
#endif // FABULOUS_ALGORITHM_H
......@@ -23,11 +23,11 @@ inline void error(const std::string &errstr, const std::string &errprefix = "err
inline void warning(const std::string &errstr)
{
std::cerr << Color::warning << "warning: "
<< Color::reset << errstr;
<< Color::reset << errstr << "\n";
}
/** \brief display an debug message */
inline void debug( const std::string & errstr)
inline void debug(const std::string &errstr)
{
#ifdef FABULOUS_DEBUG_MODE
std::cerr << Color::debug << "DEBUG: "
......@@ -37,7 +37,7 @@ inline void debug( const std::string & errstr)
#endif
}
/** \brief display an debug message */
/** \brief display an note message */
inline void note(const std::string &errstr)
{
std::cerr << Color::note << "note: "
......@@ -45,13 +45,13 @@ inline void note(const std::string &errstr)
}
/** \brief display an error message and exit with an error code */
inline void fatal_error(const std::string &errstr)
inline void fatal_error(const std::string &errstr, int error_code = EXIT_FAILURE)
{
error(errstr, "fatal error: ");
#ifdef FABULOUS_DEBUG_MODE
abort();
#endif
exit(EXIT_FAILURE);
exit(error_code);
}
/** \brief strip directory name */
......
......@@ -60,7 +60,7 @@ class Logger {
Timer _timer;
bool _log_real_residual;
int _global_iteration;
int _current_nb_mvp;
int _total_mvp;
std::vector<Entry> _iterations;
std::vector<double> _array;
......@@ -70,17 +70,18 @@ public:
Logger(bool log_real_residual = false):
_log_real_residual(log_real_residual),
_global_iteration(0),
_current_nb_mvp(0)
_total_mvp(0)
{
}
int get_nb_iterations() const { return _iterations.size(); }
int get_nb_mvp() const { return _total_mvp; }
void reset()
{
_iterations.clear();
_global_iteration = 0;
_current_nb_mvp = 0;
_total_mvp = 0;
}
const Entry &get_entry(int i) const { return _iterations.at(i); }
......@@ -97,7 +98,7 @@ public:
/**
* \brief Add an iteration to the list.
*/
void notify_iteration_end(int id, int size_krylov_space, int inSizeBlock,
void notify_iteration_end(int id, int size_krylov_space, int block_size,
std::tuple<P, P, bool> min_max_conv)
{
auto min = std::get<0>(min_max_conv);
......@@ -105,11 +106,11 @@ public:
Entry entry = Entry{
_global_iteration+1, id, size_krylov_space,
_current_nb_mvp+inSizeBlock, inSizeBlock, min, max, 0, 0,
_total_mvp+block_size, block_size, min, max, 0, 0,
_timer.get_time() - _last_timestamp
};
_iterations.push_back(entry);
_current_nb_mvp += inSizeBlock;
_total_mvp += block_size;
++ _global_iteration;
FABULOUS_NOTE("\t\t "<<Color::red<<"MIN="<<Color::reset<<std::scientific<<min
......
......@@ -2,6 +2,7 @@
#define FABULOUS_META_HPP
#include <type_traits>
#include <complex>
namespace fabulous {
namespace meta {
......@@ -50,6 +51,17 @@ Type<T> get_type(const T&)
return Type<T>{};
}
template<class U>
constexpr bool is_complex(const Type<U>&)
{
return false;
}
template<class U>
constexpr bool is_complex(const Type<std::complex<U>>&)
{
return true;
}
}; // namespace meta
......
......@@ -9,27 +9,32 @@ library(ggplot2)
library(reshape2)
library(tools)
data <- read.table(args[1], header=T)
data$ID <- basename(file_path_sans_ext(args[1]))
df <- read.table(args[1], header=T)
df$ID <- basename(file_path_sans_ext(args[1]))
if (length(args) >= 2) {
for (i in 2:length(args)) {
d <- read.table(args[i], header=T)
d$ID <- basename(file_path_sans_ext(args[i]))
data <- rbind(data, d)
df <- rbind(df, d)
}
}
id = data$ID;
data$ID1 <- paste(id, "minRes")
data$ID2 <- paste(id, "maxRes")
data$ID3 <- paste(id, "minRealRes")
data$ID4 <- paste(id, "maxRealRes")
plot_fabulous <- function(df) {
id = df$ID;
df$ID1 <- paste(id, "minRes")
df$ID2 <- paste(id, "maxRes")
df$ID3 <- paste(id, "minRealRes")
df$ID4 <- paste(id, "maxRealRes")
p <- ggplot(data, aes(nb_mvp)) +
geom_line(aes(y=minRes, color=ID1)) +
geom_line(aes(y=maxRes, color=ID2)) +
geom_line(aes(y=minRealRes, color=ID3)) +
geom_line(aes(y=maxRealRes, color=ID4))
p <- ggplot(df, aes(nb_mvp)) +
geom_line(aes(y=minRes, color=ID1)) +
geom_line(aes(y=maxRes, color=ID2)) +
geom_line(aes(y=minRealRes, color=ID3)) +
geom_line(aes(y=maxRealRes, color=ID4))
pl = p + scale_y_log10()
return (pl);
}
p + scale_y_log10()
pl <- plot_fabulous(df)
pl
......@@ -50,6 +50,12 @@ private:
M.at(col, row) = val;
}
bool header_match(const std::vector<std::string> &header, std::string token)
{
return ( std::find(header.begin(), header.end(), token) != header.end() );
}
template<class Matrix>
void LoadMatrix(Matrix &A, std::ifstream &file, const std::string &filename)
{
......@@ -61,11 +67,17 @@ private:
std::string line;
getline(file, line);
auto v = split(line);
if (v[0] != "%%MatrixMarket")
auto header = split(line);
if (header[0] != "%%MatrixMarket")
::fabulous::fatal_error(filename+": unrecognized file format.");
if (std::find(v.begin(), v.end(), "symmetric") != v.end())
if (header_match(header, "symmetric"))
symmetric = true;
using S = typename Matrix::value_type;
if (meta::is_complex(meta::Type<S>{}) && !header_match(header, "complex"))
FABULOUS_FATAL_ERROR("This is a complex parser but file header does have \"complex\"");
else if (!meta::is_complex(meta::Type<S>{}) && !header_match(header, "real"))
FABULOUS_FATAL_ERROR("This is a real parser but file header does have \"real\"");
do {
getline(file,line);
......@@ -74,10 +86,10 @@ private:
} while (commentLine);
std::cout << std::flush;
v = split(line);
int N = std::stoi(v[0]);
int M = std::stoi(v[1]);
int NNZ = std::stoi(v[2]);
auto size = split(line);
int N = std::stoi(size[0]);
int M = std::stoi(size[1]);
int NNZ = std::stoi(size[2]);
//Check square ?
if (N != M) {
......
......@@ -119,8 +119,8 @@ TestRet<P> run_test_BGMRES_filelog(
auto logger = solve(eq, algo, param, ortho, restart);
double elapsed = Timer::get_time() - tic;
int nb = logger.get_nb_iterations();
std::cout<<"Return value / max_ite : " << nb << "/" << maxMVP << std::endl;
int nb = logger.get_nb_mvp();
std::cout<<"Return value / max_mvp : " << nb << "/" << maxMVP << std::endl;
std::stringstream ID, file_ss;
ID << filename <<"_"
......
......@@ -19,6 +19,7 @@ set(FABULOUS_TEST_SRC
testGeneEigenValPbReal.cpp
testBGMResRuhe.cpp
testOrthoSchemes.cpp
test_RUHE_MGS_bcsstk_EXACT_BREAKDOWN.cpp
)
set(FABULOUS_LIBS_FOR_TESTS
......
......@@ -15,40 +15,40 @@ int main(int argc, char *argv[])
if (argc > 1)
nbRHS = std::stoi(args[1]);
using Primary= double;
using Scalar = Primary;
using Matrix = UserInputMatrix<Scalar>;
using Block = Block<Scalar>;
using P= double;
using S = P;
using Matrix = UserInputMatrix<S>;
using Block = Block<S>;
Matrix mat;
MatrixMarketLoader mml;
mml.LoadMatrix(mat, "../data/sherman4.mtx");
int dim = mat.size();
Primary normA = mat.get_inf_norm();
P normA = mat.get_inf_norm();
std::cout<<"Matrix Norm inf ::\t"<<normA<<"\n";
Block XExact{nbRHS, dim};
Block RHS{nbRHS, dim};
Block XExact{dim, nbRHS};
Block RHS{dim, nbRHS};
RandomMatrixLoader random_ld;
random_ld.LoadMatrix(XExact);
mat.MatBlockVect(XExact,RHS);
int maxNbMatVect = 3000, restart = 0, nbEigenVector = 15;
int max_mvp = 3000, max_space = 0, nbEigenVector = 15;
if (args.size() > 2)
restart = std::stoi(args[2]);
max_space = std::stoi(args[2]);
else
restart = dim+1;
max_space = dim+1;
if (args.size() > 3)
nbEigenVector = std::stoi(args[3]);
std::vector<Primary> epsilon{1e-9};
Scalar target{0.0};
std::vector<P> epsilon{1e-9};
S target{0.0};
auto r1 = run_test_BGMRES_filelog(
"sherman4",
mat, RHS, XExact, epsilon,
arnoldi_std(), maxNbMatVect, restart,
arnoldi_std(), max_mvp, max_space,
orthogonalization(OrthoScheme::MGS, OrthoType::BLOCK),
classic_restart()
);
......@@ -56,7 +56,7 @@ int main(int argc, char *argv[])
// auto r2 = run_test_BGMRES_filelog(
// "sherman4_IB",
// mat, RHS, XExact, epsilon,
// arnoldi_ib(), maxNbMatVect, restart,
// arnoldi_ib(), max_mvp, max_space,
// orthogonalization(OrthoScheme::MGS, OrthoType::BLOCK),
// classic_restart()
// );
......@@ -64,7 +64,7 @@ int main(int argc, char *argv[])
// auto r3 = run_test_BGMRES_filelog(
// "sherman4_QR",
// mat, RHS, XExact, epsilon,
// arnoldi_qr(), maxNbMatVect, restart,
// arnoldi_qr(), max_mvp, max_space,
// orthogonalization(OrthoScheme::MGS, OrthoType::BLOCK),
// classic_restart()
// );
......@@ -72,7 +72,7 @@ int main(int argc, char *argv[])
// auto r4 = run_test_BGMRES_filelog(
// "sherman4_DR",
// mat, RHS, XExact, epsilon,
// arnoldi_std(), maxNbMatVect, restart,
// arnoldi_std(), max_mvp, max_space,
// orthogonalization(OrthoScheme::MGS, OrthoType::BLOCK),
// deflated_restart(nbEigenVector, target)
// );
......@@ -82,7 +82,7 @@ int main(int argc, char *argv[])
//Generate | b_i - Ax_i | / |b_i|
Block AxSol{nbRHS,dim};
mat.MatBlockVect(sol,AxSol);
std::vector<Primary> normsAx_iMinusb_i;
std::vector<P> normsAx_iMinusb_i;
for (int i=0; i<nbRHS; ++i) {
//Ax_i - b_i
......
......@@ -15,20 +15,20 @@ int main(int argc, char *argv[])
std::cout << arg << "\n";
int nbRHS = 6;
int inRestart = 200;
int MaxNbMVP = 10000;
int max_space_step = 200;
int max_mvp = 10000;
if (args.size() > 1)
nbRHS = std::stoi(args[1]);
//Matrices
std::vector<std::string> matrices_C{"young1c", "young4c", "bidiagonalmatrix4"};
std::vector<std::string> matrices_C{"young1c", "young4c" /* ,"bidiagonalmatrix4"*/};
std::vector<std::string> matrices_R{"bcsstk14", "sherman4"};
std::string location{"../data/"};
std::string extension{".mtx"};
//Loop over complex matrix
// Loop over complex matrix
for (auto &matrix_name : matrices_C) {