Commit f41a940d authored by MIJIEUX Thomas's avatar MIJIEUX Thomas

Add description of used MatrixMarket test files

Minor fix on logger mvp count
parent caab0740
......@@ -2,51 +2,89 @@
# -*- coding: utf-8 -*-
* RESULTS
** Test files
[[file:./data/bcsstk14.mtx][bcsstk14.mtx]]
[[file:./data/sherman4.mtx][sherman4.mtx]]
[[file:./data/sherman5.mtx][sherman5.mtx]]
[[file:./data/young1c.mtx][young1c.mtx]]
[[file:./data/young4c.mtx][young4c.mtx]]
[[file:./data/bidiagonalmatrix1.mtx][bidiagonalmatrix1.mtx]]
[[file:./data/bidiagonalmatrix2.mtx][bidiagonalmatrix2.mtx]]
[[file:./data/bidiagonalmatrix3.mtx][bidiagonalmatrix3.mtx]]
[[file:./data/bidiagonalmatrix4.mtx][bidiagonalmatrix3.mtx]]
*** Downloads
#+begin_src shell :session *downloads* :results silent
files=(
"ftp://math.nist.gov/pub/MatrixMarket2/Harwell-Boeing/bcsstruc2/bcsstk14.mtx.gz"
"ftp://math.nist.gov/pub/MatrixMarket2/Harwell-Boeing/sherman/sherman4.mtx.gz"
"ftp://math.nist.gov/pub/MatrixMarket2/Harwell-Boeing/sherman/sherman5.mtx.gz"
"ftp://math.nist.gov/pub/MatrixMarket2/Harwell-Boeing/acoust/young1c.mtx.gz"
"ftp://math.nist.gov/pub/MatrixMarket2/Harwell-Boeing/acoust/young4c.mtx.gz"
)
mkdir -p ./data/
for i in "${files[@]}"; do
filename=$(basename "$i")
filename="${filename%.*}" # strip gz extension
curl $i | zcat > ./data/$filename
done
#+end_src
** Run Examples
#+begin_src bash :session foo :results output :exports both
cd fabulous/build/src/test_core/
./testMatrixMarket 100 500
cd fabulous/build/src/test_core/
./testMatrixMarket 100 500
# most matrix market with all ortho scheme + variable max_space (Krylov search space) parameter:
./testOrthoSchemes
# 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
# 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()
#+begin_src R :session *R* :results output graphics :file (org-babel-temp-file "figure" ".png") :exports both
#!/usr/bin/env Rscript
args <- commandArgs(TRUE)
if (length(args) < 1) {
stop("need filename as parameter")
}
library(ggplot2)
library(reshape2)
library(tools)
df <- read.table(args[1], header=T)
name <- basename(file_path_sans_ext(args[1]))
df$ID <- name
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]))
df <- rbind(df, d)
}
}
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(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);
}
pl <- plot_fabulous(df)
pl
#+end_src
......@@ -216,7 +216,7 @@ public:
bool convergence = std::get<2>(MinMaxConv);
_size_krylov_space = _base.get_nb_vect();
_logger.notify_iteration_end(
_iteration_count, _size_krylov_space, _nbRHS, MinMaxConv);
_iteration_count, _size_krylov_space, _nb_mvp, MinMaxConv);
_logger.log_real_residual(A, B, X);
while (!convergence && _size_krylov_space < max_krylov_space_size) {
......@@ -237,7 +237,7 @@ public:
bool proj_convergence = std::get<2>(MinMaxConv);
convergence = check_real_residual_if(proj_convergence, A, B, X, epsilon);
_logger.notify_iteration_end(
_iteration_count, _size_krylov_space, _nbRHS, MinMaxConv);
_iteration_count, _size_krylov_space, _nb_mvp, MinMaxConv);
_logger.log_real_residual(A, B, X, _base, _hess);
}
......
......@@ -402,7 +402,7 @@ public:
auto MinMaxConv = inexact_breakdown_on_R0(A, R0, epsilon, inv_epsilon, L, WP);
_logger.notify_iteration_end(
_iteration_count, _size_krylov_space, _nbRHS, MinMaxConv);
_iteration_count, _size_krylov_space, _nb_mvp, MinMaxConv);
_logger.log_real_residual(A, B, X);
if (_nb_direction_kept == 0)
......@@ -439,7 +439,7 @@ public:
}
_size_krylov_space = _base.get_nb_vect();
_logger.notify_iteration_end(
_iteration_count, _size_krylov_space, WP.get_size_W(), MinMaxConv);
_iteration_count, _size_krylov_space, _nb_mvp, MinMaxConv);
_logger.log_real_residual(A, B, X, _base, L);
}
......
......@@ -4,6 +4,7 @@
#include <cassert>
#include "fabulous/utils/Utils.hpp"
#include "fabulous/kernel/LapackInterface.hpp"
namespace fabulous {
......@@ -24,8 +25,8 @@ public:
*
* \note R MUST BE already allocated ( Q.get_nb_col() x Q.get_nb_col() )
*/
template< class Matrix, class S >
static void InPlaceQRFactoMGS_User(Block<S> &Q, Block<S> &R, const Matrix &A)
template< class Matrix, template <class> class Block, class S >
static void InPlaceQRFactoMGS_User(Block<S> &Q, Block<S>&R, const Matrix &A)
{
int M = Q.get_nb_row();
int N = Q.get_nb_col();
......
......@@ -61,6 +61,7 @@ class Logger {
bool _log_real_residual;
int _global_iteration;
int _total_mvp;
int _local_mvp;
std::vector<Entry> _iterations;
std::vector<double> _array;
......@@ -68,9 +69,10 @@ class Logger {
public:
Logger(bool log_real_residual = false):
_log_real_residual(log_real_residual),
_global_iteration(0),
_total_mvp(0)
_log_real_residual{log_real_residual},
_global_iteration{0},
_total_mvp{0},
_local_mvp{0}
{
}
......@@ -80,8 +82,10 @@ public:
void reset()
{
_iterations.clear();
_array.clear();
_global_iteration = 0;
_total_mvp = 0;
_local_mvp = 0;
}
const Entry &get_entry(int i) const { return _iterations.at(i); }
......@@ -93,24 +97,28 @@ public:
std::cout<<"\nStart of Iteration ["
<<Color::note<<iteration_count<<" :: "
<<size_krylov_space<<Color::reset<<"]\n";
if (iteration_count == 0)
_local_mvp = 0;
}
/**
* \brief Add an iteration to the list.
*/
void notify_iteration_end(int id, int size_krylov_space, int block_size,
void notify_iteration_end(int id, int size_krylov_space, int local_mvp,
std::tuple<P, P, bool> min_max_conv)
{
auto min = std::get<0>(min_max_conv);
auto max = std::get<1>(min_max_conv);
int mvp = local_mvp - _local_mvp;
_local_mvp = local_mvp;
_total_mvp += mvp;
Entry entry = Entry{
_global_iteration+1, id, size_krylov_space,
_total_mvp+block_size, block_size, min, max, 0, 0,
_total_mvp, mvp, min, max, 0, 0,
_timer.get_time() - _last_timestamp
};
_iterations.push_back(entry);
_total_mvp += block_size;
++ _global_iteration;
FABULOUS_NOTE("\t\t "<<Color::red<<"MIN="<<Color::reset<<std::scientific<<min
......
......@@ -10,7 +10,8 @@ library(reshape2)
library(tools)
df <- read.table(args[1], header=T)
df$ID <- basename(file_path_sans_ext(args[1]))
name <- basename(file_path_sans_ext(args[1]))
df$ID <- name
if (length(args) >= 2) {
for (i in 2:length(args)) {
......
......@@ -129,8 +129,14 @@ TestRet<P> run_test_BGMRES_filelog(
file_ss << "../data/res/"<< ID.str() << ".res";
std::ofstream file{file_ss.str(), std::ios::out};
if (file.is_open())
if (file.is_open()) {
logger.print(ID.str(), file);
file.close();
} else {
std::stringstream errorsstr;
errorsstr << file_ss.str() <<": " <<strerror(errno);
warning(errorsstr.str());
}
auto minmax = compare_blocks(X, XExact);
return TestRet<P>{elapsed, nb, minmax, maxMVP};
}
......
......@@ -19,6 +19,7 @@ set(FABULOUS_TEST_SRC
testGeneEigenValPbReal.cpp
testBGMResRuhe.cpp
testOrthoSchemes.cpp
testOrthoSchemes2.cpp
test_RUHE_MGS_bcsstk_EXACT_BREAKDOWN.cpp
)
......
......@@ -26,7 +26,7 @@ int main(int argc, char *argv[])
mml.LoadMatrix(mat, filename);
int dim = mat.size();
int nbRHS = 10;
int nbRHS = 100;
int max_space = 0;
int nbEigenVector = 20;
int max_mvp = 10000;
......@@ -62,9 +62,10 @@ int main(int argc, char *argv[])
run_test_BGMRES_filelog(
"young1c",
mat, RHS, XExact, epsilon,
arnoldi_qrdr(), max_mvp, max_space,
orthogonalization() + OrthoScheme::IMGS + OrthoType::BLOCK,
deflated_restart(nbEigenVector, target)
arnoldi_std(), max_mvp, max_space,
orthogonalization() + OrthoScheme::ICGS + OrthoType::BLOCK,
classic_restart()
//deflated_restart(nbEigenVector, target)
);
return 0;
}
......@@ -22,7 +22,7 @@ int main(int argc, char *argv[])
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/"};
......
#include <iostream>
#include <sstream>
#include "../test_common/Test.hpp"
#include "../test_common/UserInputMatrix.hpp"
#include "../test_common/MatrixMarketLoader.hpp"
#include "../test_common/RandomMatrixLoader.hpp"
using namespace fabulous;
namespace {
template<class Algo>
void file_test(Algo algo, int nbRHS)
{
using P = double;
using S = std::complex<P>;
using Matrix = UserInputMatrix<S>;
using Block = Block<S>;
std::vector<P> epsilon = {1e-6};
MatrixMarketLoader mml;
Matrix A;
std::string filename{"../data/young1c.mtx"};
mml.LoadMatrix(A, filename);
int max_space_step = 200;
int max_mvp = 10000;
int dim = A.size();
Block XExact{dim, nbRHS};
Block RHS{dim, nbRHS};
// for (int i=0; i<nbRHS; ++i)
// XExact.at(i,i) = S{1};
RandomMatrixLoader rnld;
rnld.LoadMatrix(XExact);
A.MatBlockVect(XExact, RHS);
for (int j = 0, max_space = max_space_step; j<5; ++j, max_space += max_space_step) {
//Loop over ortho choice
for (int Ort = 0; Ort < 4; ++Ort) {
OrthoScheme scheme = static_cast<OrthoScheme>(Ort);
for (int i = 0; i < 2; ++i) {
OrthoType type = static_cast<OrthoType>(i);
std::stringstream ss;
ss << "young1c_r="<<std::to_string(max_space);
std::cout<<"Computing "<<ss.str()<<"_"<<scheme<<"_"<<type<<"...\n";
run_test_BGMRES_filelog(
ss.str(),
A, RHS, XExact, epsilon,
algo, max_mvp, max_space,
orthogonalization(scheme, type),
classic_restart()
);
std::cout<<"\n------------------------------------------------------------\n\n";
}
}
}
}
};
int main(int argc, char *argv[])
{
std::vector<std::string > args{argv, argv+argc};
for (auto arg : args)
std::cout << arg << "\n";
int nbRHS = 6;
if (args.size() > 1)
nbRHS = std::stoi(args[1]);
file_test(arnoldi_ib(), nbRHS);
file_test(arnoldi_std(), nbRHS);
}
......@@ -13,10 +13,6 @@ int main(int argc, char *argv[])
for (auto arg : args)
std::cout << arg << "\n";
int nbRHS = 16;
if (argc > 1)
nbRHS = std::stoi(args[1]);
using P = double;
using S = P;
using Matrix = UserInputMatrix<S>;
......@@ -30,6 +26,10 @@ int main(int argc, char *argv[])
P normA = mat.get_inf_norm();
std::cout<<"Matrix Norm inf ::\t"<<normA<<"\n";
int nbRHS = 6;
if (argc > 1)
nbRHS = std::stoi(args[1]);
Block XExact{dim, nbRHS};
Block RHS{dim, nbRHS};
......@@ -43,12 +43,11 @@ int main(int argc, char *argv[])
if (args.size() > 2)
max_space = std::stoi(args[2]);
else
max_space = dim+1;
std::vector<P> epsilon{1e-6};
S target{0.0};
max_space = 200;
std::vector<P> epsilon{1e-4};
auto r1 = run_test_BGMRES_filelog(
"sherman4",
auto r1 = run_test_BGMRES(
/* "sherman4", */
mat, RHS, XExact, epsilon,
arnoldi_ib(), max_mvp, max_space,
orthogonalization() + OrthoScheme::CGS + OrthoType::RUHE,
......@@ -56,10 +55,11 @@ int main(int argc, char *argv[])
);
//Display timers
std::cout<<"Elapsed time for each method::\n";
print_elapsed_time(r1, "Standard version");
std::cout <<"Forward error normalized by norm inf of A : Min\t"
<<r1.minmax.first/normA<<" Max\t"<<r1.minmax.second/normA<<"\n";
// std::cout<<"Elapsed time for each method::\n";
// print_elapsed_time(r1, "Standard version");
// std::cout <<"Forward error normalized by norm inf of A : Min\t"
// <<r1.minmax.first/normA<<" Max\t"<<r1.minmax.second/normA<<"\n";
print_logs(r1, "bcctsk_CGS_RUHE_CR_200");
return 0;
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment