Commit 1f23f2c7 authored by MIJIEUX Thomas's avatar MIJIEUX Thomas

Add COO and CSC matrix + test with LightInTissue

parent 0c0f7d6f
......@@ -55,7 +55,7 @@ make
** Influence of restart(m) parameter
*** run test case
#+begin_src sh :session *TEST* :results silent
cd ${WORKDIR}/build/src/test_core/
cd ${WORKDIR}/build/src/test_basic/
mkdir -p ../data/res
sizes=(200 400 600 800 1000)
for siz in "${sizes[@]}"; do
......@@ -83,7 +83,7 @@ ggplot(df, aes(x=nb_mvp, y=maxRes, color=name)) +
** Influence of incremental QR factorization
*** run test case
#+begin_src sh :session *TEST* :results silent
cd ${WORKDIR}/build/src/test_core/
cd ${WORKDIR}/build/src/test_basic/
mkdir -p ../data/res
./testMatrixMarketChoice -t BLOCK -s CGS -m 900 -A STDDR -u -o "Full_GELS"
./testMatrixMarketChoice -t BLOCK -s CGS -m 900 -A QR -u -o "GELS_with_Incremental_QR_factorization"
......@@ -99,7 +99,7 @@ ggplot(df, aes(x=global_iteration, y=least_square_time, color=name)) +
** Influence of the Algorithm
*** run test case
#+begin_src sh :session *TEST* :results silent
cd ${WORKDIR}/build/src/test_core/
cd ${WORKDIR}/build/src/test_basic/
mkdir -p ../data/res
./testMatrixIW_DR
#+end_src
......@@ -122,7 +122,7 @@ ggplot(df, aes(x=nb_mvp)) +
** young1c (nrhs=6, m=90, k=5)
*** run test case
#+begin_src sh :session *TEST* :results silent
cd ${WORKDIR}/build/src/test_core/
cd ${WORKDIR}/build/src/test_basic/
mkdir -p ../data/res
#nbRHS=6; maxSpace=90
./testMatrixMarketChoice -n 6 -t BLOCK -s MGS -m 90 -e 1e-6 -A QR -u -o "BGMRES"
......@@ -152,10 +152,35 @@ ggplot(df, aes(x=nb_mvp)) +
scale_y_log10() + geom_hline(aes(yintercept=1e-6, color="threshold")) +
ggtitle("young1c (nrhs=6, m=90 k=5)") + ylab(TeX("$\\eta_b$ (min, max)"))
#+end_src
** lightInTissue (nrhs=6, m=90, k=5)
*** run test case
#+begin_src sh :session *TEST* :results silent
cd ${WORKDIR}/build/src/test_basic/
mkdir -p ../data/res
#nbRHS=6; maxSpace=90
./testMatrixMarketChoice -f ../data/New/lightINtissue.mtx -k CSC -n 6 -m 90 -A QR -e 1e-6 -u -o "BGMRES"
./testMatrixMarketChoice -f ../data/New/lightINtissue.mtx -k CSC -n 6 -m 90 -A IB -e 1e-6 -u -o "IB-BGMRES"
./testMatrixMarketChoice -f ../data/New/lightINtissue.mtx -k CSC -n 6 -m 90 -A QRDR -r DEFLATED -p 5 -e 1e-6 -u -o "BGMRES-DR"
./testMatrixMarketChoice -f ../data/New/lightINtissue.mtx -k CSC -n 6 -m 90 -A IBDR -r DEFLATED -p 5 -e 1e-6 -u -o "IB-BGMRES-DR"
#+end_src
*** plot the graphic
#+begin_src R :results output graphics :file (org-babel-temp-file "figure" ".png") :exports both :width 600 :height 400 :session *R*
library(ggplot2)
library(latex2exp)
df <- read.table("./build/src/data/res/BGMRES.res", header=T)
df <- rbind(df, read.table("./build/src/data/res/IB-BGMRES.res", header=T))
df <- rbind(df, read.table("./build/src/data/res/BGMRES-DR.res", header=T))
df <- rbind(df, read.table("./build/src/data/res/IB-BGMRES-DR.res", header=T))
ggplot(df, aes(x=nb_mvp)) +
geom_line(aes(y=maxRes, color=name)) +
geom_line(aes(y=minRes, color=name)) +
scale_y_log10() + geom_hline(aes(yintercept=1e-6, color="threshold")) +
ggtitle("lightInTissue (nrhs=6, m=90 k=5)") + ylab(TeX("$\\eta_b$ (min, max)"))
#+end_src
** Influence of deflated restart (k) parameter (nrhs=6, m=90)
*** run test case
#+begin_src sh :session *TEST* :results silent
cd ${WORKDIR}/build/src/test_core/
cd ${WORKDIR}/build/src/test_basic/
mkdir -p ../data/res
#nbRHS=6; maxSpace=90
./testMatrixMarketChoice -n 6 -t RUHE -s MGS -m 90 -e 1e-6 -A QRDR -r DEFLATED -p 5 -u -o "BGMRES-DR(05)"
......
......@@ -207,6 +207,14 @@ public:
++ _nb_restart;
print_iteration_end_info(arnoldi.get_krylov_space_size());
}
if (_mvp == 0) {
warning(
"No matrix multiplication was performed. "
"You may want to check the algorithm parameters"
);
}
X.scale(normB); B.scale(normB);
return _mvp;
}
......
......@@ -14,6 +14,10 @@ ELSE()
)
ENDIF()
IF(BLA_VENDOR MATCHES "Intel*")
ADD_DEFINITIONS( -DFABULOUS_USE_MKL_SPARSE )
ENDIF(BLA_VENDOR MATCHES "Intel*")
ADD_SUBDIRECTORY(api)
ADD_SUBDIRECTORY(test_basic)
ADD_SUBDIRECTORY(test_api)
......
SET(FABULOUS_BASIC_SRC
testBlockClass.cpp
testBlockSvd.cpp
testBlockWP.cpp
testMatrixMarket.cpp
testMatrixMarket_CSR.cpp
testMatrixMarket_CSC.cpp
testMatrixMarketChoice.cpp
testMatrixMarketReal.cpp
testMatrixMarket_DR.cpp
......
......@@ -49,6 +49,8 @@ const char *gengetopt_args_info_help[] = {
" -r, --restart=ENUM restarting type (possible values=\"CLASSIC\",\n \"DEFLATED\" default=`CLASSIC')",
" -p, --nb-eigen-pair=INT number of eigen pair for deflated restarting\n (default=`20')",
" -P, --precond enable a preconditioner working with diagonally\n dominant matrices (default=off)",
" -j, --parser=ENUM set the matrix parser (MatrixMarket;Innovation\n Works;Random) (possible values=\"MM\", \"IW\",\n \"RD\" default=`MM')",
" -k, --storage=ENUM set the matrix storage type (possible\n values=\"DENSE\", \"CSC\", \"COO\"\n default=`DENSE')",
" -o, --output=STRING set the output filename, if not set the input\n filename will be used",
" -u, --suffix automatically append parameter as a suffix to output\n filename (default=on)",
0
......@@ -77,6 +79,8 @@ const char *cmdline_parser_scheme_values[] = {"CGS", "ICGS", "MGS", "IMGS", 0};
const char *cmdline_parser_arithmetic_values[] = {"REAL_FLOAT", "REAL_DOUBLE", "COMPLEX_FLOAT", "COMPLEX_DOUBLE", 0}; /*< Possible values for arithmetic. */
const char *cmdline_parser_algo_values[] = {"STDDR", "CHAMDR", "QR", "CHAMQR", "QRDR", "IB", "IBDR", 0}; /*< Possible values for algo. */
const char *cmdline_parser_restart_values[] = {"CLASSIC", "DEFLATED", 0}; /*< Possible values for restart. */
const char *cmdline_parser_parser_values[] = {"MM", "IW", "RD", 0}; /*< Possible values for parser. */
const char *cmdline_parser_storage_values[] = {"DENSE", "CSC", "COO", 0}; /*< Possible values for storage. */
static char *
gengetopt_strdup (const char *s);
......@@ -99,6 +103,8 @@ void clear_given (struct gengetopt_args_info *args_info)
args_info->restart_given = 0 ;
args_info->nb_eigen_pair_given = 0 ;
args_info->precond_given = 0 ;
args_info->parser_given = 0 ;
args_info->storage_given = 0 ;
args_info->output_given = 0 ;
args_info->suffix_given = 0 ;
}
......@@ -132,6 +138,10 @@ void clear_args (struct gengetopt_args_info *args_info)
args_info->nb_eigen_pair_arg = 20;
args_info->nb_eigen_pair_orig = NULL;
args_info->precond_flag = 0;
args_info->parser_arg = parser_arg_MM;
args_info->parser_orig = NULL;
args_info->storage_arg = storage_arg_DENSE;
args_info->storage_orig = NULL;
args_info->output_arg = NULL;
args_info->output_orig = NULL;
args_info->suffix_flag = 1;
......@@ -158,8 +168,10 @@ void init_args_info(struct gengetopt_args_info *args_info)
args_info->restart_help = gengetopt_args_info_help[12] ;
args_info->nb_eigen_pair_help = gengetopt_args_info_help[13] ;
args_info->precond_help = gengetopt_args_info_help[14] ;
args_info->output_help = gengetopt_args_info_help[15] ;
args_info->suffix_help = gengetopt_args_info_help[16] ;
args_info->parser_help = gengetopt_args_info_help[15] ;
args_info->storage_help = gengetopt_args_info_help[16] ;
args_info->output_help = gengetopt_args_info_help[17] ;
args_info->suffix_help = gengetopt_args_info_help[18] ;
}
......@@ -256,6 +268,8 @@ cmdline_parser_release (struct gengetopt_args_info *args_info)
free_string_field (&(args_info->algo_orig));
free_string_field (&(args_info->restart_orig));
free_string_field (&(args_info->nb_eigen_pair_orig));
free_string_field (&(args_info->parser_orig));
free_string_field (&(args_info->storage_orig));
free_string_field (&(args_info->output_arg));
free_string_field (&(args_info->output_orig));
......@@ -359,6 +373,10 @@ cmdline_parser_dump(FILE *outfile, struct gengetopt_args_info *args_info)
write_into_file(outfile, "nb-eigen-pair", args_info->nb_eigen_pair_orig, 0);
if (args_info->precond_given)
write_into_file(outfile, "precond", 0, 0 );
if (args_info->parser_given)
write_into_file(outfile, "parser", args_info->parser_orig, cmdline_parser_parser_values);
if (args_info->storage_given)
write_into_file(outfile, "storage", args_info->storage_orig, cmdline_parser_storage_values);
if (args_info->output_given)
write_into_file(outfile, "output", args_info->output_orig, 0);
if (args_info->suffix_given)
......@@ -650,12 +668,14 @@ cmdline_parser_internal (
{ "restart", 1, NULL, 'r' },
{ "nb-eigen-pair", 1, NULL, 'p' },
{ "precond", 0, NULL, 'P' },
{ "parser", 1, NULL, 'j' },
{ "storage", 1, NULL, 'k' },
{ "output", 1, NULL, 'o' },
{ "suffix", 0, NULL, 'u' },
{ 0, 0, 0, 0 }
};
c = getopt_long (argc, argv, "hVt:s:i:f:n:m:M:e:a:A:r:p:Po:u", long_options, &option_index);
c = getopt_long (argc, argv, "hVt:s:i:f:n:m:M:e:a:A:r:p:Pj:k:o:u", long_options, &option_index);
if (c == -1) break; /* Exit from `while (1)' loop. */
......@@ -824,6 +844,30 @@ cmdline_parser_internal (
additional_error))
goto failure;
break;
case 'j': /* set the matrix parser (MatrixMarket;Innovation Works;Random). */
if (update_arg( (void *)&(args_info->parser_arg),
&(args_info->parser_orig), &(args_info->parser_given),
&(local_args_info.parser_given), optarg, cmdline_parser_parser_values, "MM", ARG_ENUM,
check_ambiguity, override, 0, 0,
"parser", 'j',
additional_error))
goto failure;
break;
case 'k': /* set the matrix storage type. */
if (update_arg( (void *)&(args_info->storage_arg),
&(args_info->storage_orig), &(args_info->storage_given),
&(local_args_info.storage_given), optarg, cmdline_parser_storage_values, "DENSE", ARG_ENUM,
check_ambiguity, override, 0, 0,
"storage", 'k',
additional_error))
goto failure;
break;
case 'o': /* set the output filename, if not set the input filename will be used. */
......
......@@ -39,6 +39,8 @@ enum enum_scheme { scheme__NULL = -1, scheme_arg_CGS = 0, scheme_arg_ICGS, schem
enum enum_arithmetic { arithmetic__NULL = -1, arithmetic_arg_REAL_FLOAT = 0, arithmetic_arg_REAL_DOUBLE, arithmetic_arg_COMPLEX_FLOAT, arithmetic_arg_COMPLEX_DOUBLE };
enum enum_algo { algo__NULL = -1, algo_arg_STDDR = 0, algo_arg_CHAMDR, algo_arg_QR, algo_arg_CHAMQR, algo_arg_QRDR, algo_arg_IB, algo_arg_IBDR };
enum enum_restart { restart__NULL = -1, restart_arg_CLASSIC = 0, restart_arg_DEFLATED };
enum enum_parser { parser__NULL = -1, parser_arg_MM = 0, parser_arg_IW, parser_arg_RD };
enum enum_storage { storage__NULL = -1, storage_arg_DENSE = 0, storage_arg_CSC, storage_arg_COO };
/** @brief Where the command line options are stored */
struct gengetopt_args_info
......@@ -83,6 +85,12 @@ struct gengetopt_args_info
const char *nb_eigen_pair_help; /**< @brief number of eigen pair for deflated restarting help description. */
int precond_flag; /**< @brief enable a preconditioner working with diagonally dominant matrices (default=off). */
const char *precond_help; /**< @brief enable a preconditioner working with diagonally dominant matrices help description. */
enum enum_parser parser_arg; /**< @brief set the matrix parser (MatrixMarket;Innovation Works;Random) (default='MM'). */
char * parser_orig; /**< @brief set the matrix parser (MatrixMarket;Innovation Works;Random) original value given at command line. */
const char *parser_help; /**< @brief set the matrix parser (MatrixMarket;Innovation Works;Random) help description. */
enum enum_storage storage_arg; /**< @brief set the matrix storage type (default='DENSE'). */
char * storage_orig; /**< @brief set the matrix storage type original value given at command line. */
const char *storage_help; /**< @brief set the matrix storage type help description. */
char * output_arg; /**< @brief set the output filename, if not set the input filename will be used. */
char * output_orig; /**< @brief set the output filename, if not set the input filename will be used original value given at command line. */
const char *output_help; /**< @brief set the output filename, if not set the input filename will be used help description. */
......@@ -104,6 +112,8 @@ struct gengetopt_args_info
unsigned int restart_given ; /**< @brief Whether restart was given. */
unsigned int nb_eigen_pair_given ; /**< @brief Whether nb-eigen-pair was given. */
unsigned int precond_given ; /**< @brief Whether precond was given. */
unsigned int parser_given ; /**< @brief Whether parser was given. */
unsigned int storage_given ; /**< @brief Whether storage was given. */
unsigned int output_given ; /**< @brief Whether output was given. */
unsigned int suffix_given ; /**< @brief Whether suffix was given. */
......@@ -235,6 +245,8 @@ extern const char *cmdline_parser_scheme_values[]; /**< @brief Possible values
extern const char *cmdline_parser_arithmetic_values[]; /**< @brief Possible values for arithmetic. */
extern const char *cmdline_parser_algo_values[]; /**< @brief Possible values for algo. */
extern const char *cmdline_parser_restart_values[]; /**< @brief Possible values for restart. */
extern const char *cmdline_parser_parser_values[]; /**< @brief Possible values for parser. */
extern const char *cmdline_parser_storage_values[]; /**< @brief Possible values for storage. */
#ifdef __cplusplus
......
......@@ -3,8 +3,6 @@ package "fabulous_test"
purpose "solve equation with multiple rhs"
option "type" t "orthogonalization type"
values="RUHE","BLOCK" enum
optional default="RUHE"
......@@ -49,6 +47,14 @@ option "nb-eigen-pair" p "number of eigen pair for deflated restarting"
option "precond" P "enable a preconditioner working with diagonally dominant matrices"
flag off
option "parser" j "set the matrix parser (MatrixMarket;Innovation Works;Random)"
values="MM", "IW", "RD" enum
optional default="MM"
option "storage" k "set the matrix storage type"
values="DENSE", "CSC", "COO" enum
optional default="DENSE"
option "output" o "set the output filename, if not set the input filename will be used"
string optional
......
#include "../test_common/Test.hpp"
#include "../test_common/UserInputMatrix.hpp"
#include "../test_common/CSCMatrix.hpp"
#include "../test_common/COOMatrix.hpp"
#include "../test_common/MatrixMarketLoader.hpp"
#include "../test_common/RandomMatrixLoader.hpp"
#include "../test_common/MatrixIWLoader.hpp"
#include "cmdline.h"
#include "cmdline.c" // BAD
......@@ -38,11 +42,11 @@ OrthoType get_ortho_type(const struct gengetopt_args_info &info)
return OrthoType::RUHE;
}
template<class S, class Algo, class Restart>
void run_test_algo(const struct gengetopt_args_info &info, Algo algo, Restart restart)
template<class S, class Algo, class Restart, class Loader, class Matrix>
void run_test_storage(const struct gengetopt_args_info &info,
Algo algo, Restart restart, Loader matrix_loader, Type<Matrix>)
{
using P = typename Arithmetik<S>::primary_type;
using Matrix = UserInputMatrix<S>;
using Block = Block<S>;
int iter = info.iteration_arg;
if (iter < 2) {
......@@ -50,10 +54,9 @@ void run_test_algo(const struct gengetopt_args_info &info, Algo algo, Restart re
}
Matrix mat{info.precond_flag != 0};
MatrixMarketLoader mml;
std::string filename{info.file_arg};
std::cout<<"File ::\t"<<filename<<"\n";
mml.LoadMatrix(mat, filename);
matrix_loader.LoadMatrix(mat, filename);
int dim = mat.size();
int nbRHS = info.nrhs_arg;
......@@ -89,6 +92,28 @@ void run_test_algo(const struct gengetopt_args_info &info, Algo algo, Restart re
);
}
template<class S, class Algo, class Restart, class Loader>
void run_test_parser(const struct gengetopt_args_info &info, Algo algo, Restart restart, Loader loader)
{
switch (info.storage_arg) {
case storage_arg_DENSE: run_test_storage<S>(info, algo, restart, loader, Type<UserInputMatrix<S>>{}); break;
case storage_arg_CSC: run_test_storage<S>(info, algo, restart, loader, Type<CSCMatrix<S>>{}); break;
case storage_arg_COO: run_test_storage<S>(info, algo, restart, loader, Type<COOMatrix<S>>{}); break;
default: fatal_error("invalid storage type"); break;
}
}
template<class S, class Algo, class Restart>
void run_test_algo(const struct gengetopt_args_info &info, Algo algo, Restart restart)
{
switch (info.parser_arg) {
case parser_arg_MM: run_test_parser<S>(info, algo, restart, MatrixMarketLoader{}); break;
case parser_arg_IW: run_test_parser<S>(info, algo, restart, MatrixIWLoader{}); break;
case parser_arg_RD: run_test_parser<S>(info, algo, restart, RandomMatrixLoader{}); break;
default: fatal_error("invalid parser type"); break;
}
}
template<class S, class Restart>
void run_test_restart(const struct gengetopt_args_info &info, Restart restart)
{
......
......@@ -2,21 +2,19 @@
#include "../test_common/UserInputMatrix.hpp"
#include "../test_common/MatrixMarketLoader.hpp"
#include "../test_common/RandomMatrixLoader.hpp"
#include "../test_common/CSRMatrix.hpp"
#include "../test_common/CSCMatrix.hpp"
#include "../test_common/COOMatrix.hpp"
using namespace fabulous;
int main(int argc, char *argv[])
{
std::vector<std::string> args{argv, argv+argc};
for (auto arg : args)
std::cout << arg << "\n";
atexit(abort);
using P = double;
using S = std::complex<P>;
using Matrix = Z_CSRMatrix;
using Matrix = CSCMatrix<S>;
using Block = Block<S>;
std::string filename{"../data/young1c.mtx"};
......
#ifndef FABULOUS_COO_MATRIX_HPP
#define FABULOUS_COO_MATRIX_HPP
#include <memory>
#include <vector>
#include <complex>
namespace fabulous {
template<class> class UserInputMatrix;
};
#include "fabulous/data/Block.hpp"
#include "fabulous/kernel/LapackInterface.hpp"
#include "fabulous/kernel/QR.hpp"
#include "fabulous/utils/Arithmetic.hpp"
#include "./sparse.hpp"
namespace fabulous {
/**
* Matrix Class: need to implement MatBlockVect()
* useRightPreCond(), preCondBlockVect() and DotProduct()
* member function in order to works with BGMRes algorithm
*
* An abstract class could be used to enforce that.
*/
template<class S>
class COOMatrix
{
public:
using value_type = typename Arithmetik<S>::value_type;
using primary_type = typename Arithmetik<S>::primary_type;
private:
int _m;
std::vector<int> _row;
std::vector<int> _col;
std::vector<S> _val;
public:
S& at(int i, int j)
{
_val.emplace_back(S{0.0});
_row.emplace_back(i+1);
_col.emplace_back(j+1);
return _val.back();
}
const S& at(int, int) const
{
FABULOUS_FATAL_ERROR("badat");
static S dummy;
return dummy;
}
COOMatrix(bool precond=false)
{
if (precond) {
FABULOUS_WARNING("precond parameter ignored with COOMatrix");
}
}
void resize(int m, int n, int nnz)
{
FABULOUS_DEBUG("resize");
assert(m == n);
FABULOUS_DEBUG("m="<<m<<" n="<< n<<" nnz="<<nnz);
FABULOUS_ASSERT( m == n );
_m = m;
_val.clear();
_row.clear();
_col.clear();
_val.reserve(nnz);
_row.reserve(nnz);
_col.reserve(nnz);
}
primary_type get_inf_norm() const
{
FABULOUS_FATAL_ERROR("");
return -1.0;
}
/* ---------------------------------------------------------------- */
/* --------------- CALLBACKS -------------------------------------- */
/* ---------------------------------------------------------------- */
int size() const { return _m; }
int get_nb_row() const { return size(); }
int get_nb_col() const { return size(); }
bool useRightPrecond() const { return false; }
void MatBlockVect(const Block<S> &X, Block<S> &B,
S alpha = S{1.0}, S beta = S{0.0}) const
{
assert( B.get_nb_col() == X.get_nb_col() );
assert( X.get_nb_row() == size() );
assert( B.get_nb_row() == size() );
int M = size();
int N = X.get_nb_col();
int K = M;
int NNZ = _val.size();
SparseKernI::coomm(
M, N, K, alpha,
_val.data(), _row.data(), _col.data(), NNZ,
X.get_ptr(), X.get_leading_dim(),
beta, B.get_ptr(), B.get_leading_dim()
);
}
void QRFacto(Block<S> &Q, Block<S> &R) const
{
QR::InPlaceQRFactoMGS_User(Q, R, *this);
//Algorithm::InPlaceQRFacto(Q, R);
}
/**
* \brief Preconditioner for diagonally dominant matrix
*
* This preconditionner works well for diagonally dominant matrices.
* M^{-1} == diagonal_matrix(inverse_of_diagonal_element(A))
* M^{-1} * A ~= I
*/
void PrecondBlockVect(Block<S> &, Block<S> &) const
{
}
void DotProduct(int M, int N,
const S *A, int lda,
const S *B, int ldb,
S *C, int ldc) const
{
int dim = size();
LapackKernI::Tgemm(M, N, dim, A, lda, B, ldb, C, ldc, S{1.0}, S{0.0});
}
}; // end class UserInputMatrix
}; // end namespace fabulous
#endif // FABULOUS_USER_INPUT_MATRIX_HPP
#ifndef FABULOUS_CSR_MATRIX_HPP
#define FABULOUS_CSR_MATRIX_HPP
#ifndef FABULOUS_CSC_MATRIX_HPP
#define FABULOUS_CSC_MATRIX_HPP
#include <memory>
#include <vector>
#include <complex>
#define MKL_Complex16 std::complex<double>
#include <mkl_spblas.h>
namespace fabulous {
template<class> class UserInputMatrix;
......@@ -18,6 +13,8 @@ template<class> class UserInputMatrix;
#include "fabulous/kernel/QR.hpp"
#include "fabulous/utils/Arithmetic.hpp"
#include "./sparse.hpp"
namespace fabulous {
/**
......@@ -27,10 +24,9 @@ namespace fabulous {
*
* An abstract class could be used to enforce that.
*/
class Z_CSRMatrix
template<class S>
class CSCMatrix
{
using S = std::complex<double>;
public:
using value_type = typename Arithmetik<S>::value_type;
using primary_type = typename Arithmetik<S>::primary_type;
......@@ -40,31 +36,51 @@ private:
std::vector<int> _indices;
std::vector<int> _row_ptr;
int _last_j;
public:
S& at(int i, int j)
{
if (j < _last_j) {
FABULOUS_FATAL_ERROR("matrix was not filled in CSR order");
FABULOUS_FATAL_ERROR("matrix was not filled in CSC order");
}
_values.emplace_back(S{0.0});
_indices.emplace_back(i);
int k = _values.size();
_indices.emplace_back(i+1);
if (j > _last_j) {
_row_ptr.back() = k;
_row_ptr.emplace_back(k+1);
} else {
_row_ptr.back() = k;
if (j > _last_j+1) {
FABULOUS_FATAL_ERROR("empty column in source matrix");
}
// FABULOUS_DEBUG("newline idx="<<_row_ptr.back());
_row_ptr.emplace_back(_row_ptr.back());
}
++ _row_ptr.back();
// FABULOUS_DEBUG("last_ptr="<<_row_ptr.back());
_last_j = j;
// FABULOUS_DEBUG("i="<<i<<" j="<<j);
return _values.back();
}
Z_CSRMatrix():_last_j{0} {}
const S& at(int, int) const
{
FABULOUS_FATAL_ERROR("bad at");
static S dummy;
return dummy;
}
CSCMatrix(bool precond=false):_last_j{0}
{
if (precond) {
FABULOUS_WARNING("precond parameter ignored with COOMatrix");
}
}
void resize(int m, int n, int nnz)
{
FABULOUS_DEBUG("resize");
assert(m == n);
FABULOUS_DEBUG("m="<<m<<" n="<< n<<" nnz="<<nnz);
FABULOUS_ASSERT( m == n );
_values.clear();
_indices.clear();
......@@ -72,8 +88,11 @@ public:
_values.reserve(nnz);
_indices.reserve(nnz);
FABULOUS_DEBUG("m+1="<<m+1);
_row_ptr.reserve(m+1);
_row_ptr.push_back(0);
_row_ptr.emplace_back(1);
_row_ptr.emplace_back(1);
}
primary_type get_inf_norm() const
......@@ -86,6 +105,9 @@ public:
/* --------------- CALLBACKS -------------------------------------- */
/* ---------------------------------------------------------------- */
int size() const { return _row_ptr.size()-1; }
int get_nb_row() const { return size(); }
int get_nb_col() const { return size(); }
bool useRightPrecond() const { return false; }
void MatBlockVect(const Block<S> &X, Block<S> &B,
......@@ -94,21 +116,18 @@ public:
assert( B.get_nb_col() == X.get_nb_col() );
assert( X.get_nb_row() == size() );
assert( B.get_nb_row() == size() );
assert( _row_ptr.capacity() == _row_ptr.size() );
int M = size();
int N = X.get_nb_col();
char desca[6] = {'G', 0, 'N', 'C', };
char trans = 'N';
int ldx = X.get_leading_dim();
int ldb = B.get_leading_dim();
mkl_zcsrmm(
&trans, &M, &N, &M, &alpha, desca,
_values.data(),
_indices.data(),
int K = M;
SparseKernI::cscmm(
M, N, K, alpha,
_values.data(), _indices.data(),
_row_ptr.data(), _row_ptr.data()+1,
X.get_ptr(), &ldx,
&beta, B.get_ptr(), &ldb
X.get_ptr(), X.get_leading_dim(),
beta, B.get_ptr(), B.get_leading_dim()
);
}
......@@ -135,7 +154,7 @@ public:
S *C, int ldc) const
{
int dim = size();
LapackKernI::Tgemm(M, N, dim, A, lda, B, ldb, C, ldc, S{1.0}, S{0.0} );
LapackKernI::Tgemm(M, N, dim, A, lda, B, ldb, C, ldc, S{1.0}, S{0.0});
}
}; // end class UserInputMatrix
......
......@@ -27,15 +27,15 @@ private:
std::is_same<U, std::complex<T>>::value>::type >
void load_line(Matrix &M, const std::vector<std::string> &v, bool symmetric)
{
int col = std::stoi(v[0]) - 1;
int row = std::stoi(v[1]) - 1;
int i = std::stoi(v[0]) - 1;
int j = std::stoi(v[1]) - 1;
T real = T(std::stod(v[2]));
T imag = T(std::stod(v[3]));
M.at(row, col) = U{real, imag};
M.at(i, j) = U{real, imag};
if (symmetric)
M.at(col, row) = U{real, imag};
M.at(j, i) = U{real, imag};
}
template<class Matrix,
......@@ -59,7 +59,6 @@ private:
return ( std::find(header.begin(), header.end(), token) != header.end() );
}
template<class Matrix>
void LoadMatrix(Matrix &A, std::ifstream &file, const std::string &filename)
{
......
......@@ -37,7 +37,7 @@ public:
RandomMatrixLoader(): _gen(0) {}
template<class Matrix>
void LoadMatrix(Matrix &A)
void LoadMatrix(Matrix &A, const std::string& ="")
{
using S = typename Matrix::value_type;
using P = typename Matrix::primary_type;
......
......@@ -146,8 +146,8 @@ TestRet<P> run_test_BGMRES_filelog(
file.close();
} else {
std::stringstream errorsstr;