Commit dae3024b authored by MIJIEUX Thomas's avatar MIJIEUX Thomas

Fine grained logging at Hessenberg Level

parent d12d53c7
......@@ -3,6 +3,7 @@
#+TITLE: Fabulous installation procedure (Plafrim2)
* INSTALLATION on PLAFRIM2 with chameleon
** With SPACK
*** Install SPACK
......@@ -79,3 +80,31 @@
chmod +x ./install-fabulous.sh
./install-fabulous.sh
#+end_src
* COMPILATION for development purpose
:PROPERTIES:
:header-args:shell: :var WORKDIR=fabulous_root :session *fabulous-build* :results silent
:END:
#+NAME: fabulous_root
| /home/tmijieux/fabulous/ |
You may want to load mkl into environement:
#+BEGIN_SRC shell
source /opt/intel/mkl/bin/mklvars.sh intel64
#+END_SRC
Setup the build system:
#+BEGIN_SRC shell
mkdir -p ${WORKDIR}/build
cd ${WORKDIR}/build
cmake .. -DCMAKE_BUILD_TYPE=Debug \
-DFABULOUS_DEBUG_MODE=ON \
-DFABULOUS_LAPACKE_NANCHECK=OFF # \
# -DCHAMELEON_DIR=$(spack location -i chameleon)
#+END_SRC
Eventually compile fabulous:
#+BEGIN_SRC shell
cd ${WORKDIR}/build
make -j4
#+END_SRC
......@@ -84,16 +84,20 @@ done
Section 4.5, show that the contrary is possible)
#+begin_src R
library(ggplot2)
df <- read.csv2("./build/src/data/res/r=200.res")
df <- rbind(df, read.csv2("./build/src/data/res/r=400.res"))
df <- rbind(df, read.csv2("./build/src/data/res/r=600.res"))
df <- rbind(df, read.csv2("./build/src/data/res/r=800.res"))
df <- rbind(df, read.csv2("./build/src/data/res/r=1000.res"))
df <- read.csv("./build/src/data/res/r=200.res")
df <- rbind(df, read.csv("./build/src/data/res/r=400.res"))
df <- rbind(df, read.csv("./build/src/data/res/r=600.res"))
df <- rbind(df, read.csv("./build/src/data/res/r=800.res"))
df <- rbind(df, read.csv("./build/src/data/res/r=1000.res"))
ggplot(df, aes(x=nb_mvp, y=maxRes, color=name)) +
geom_line() +
geom_hline(aes(yintercept=1e-4, color="threshold")) +
scale_y_log10() + ggtitle("Memory usage influence (young1c IB IMGS-RUHE)")
#+end_src
#+results:
[[file:/tmp/babel-1939RTd/figure1939yUv.png]]
*** Influence of incremental QR factorization
**** run test case
#+begin_src shell
......@@ -105,11 +109,15 @@ mkdir -p ../data/res
**** plot the graphic
#+begin_src R
library(ggplot2)
df <- read.csv2("./build/src/data/res/Full_GELS.res")
df <- rbind(df, read.csv2("./build/src/data/res/GELS_with_Incremental_QR_factorization.res"))
df <- read.csv("./build/src/data/res/Full_GELS.res")
df <- rbind(df, read.csv("./build/src/data/res/GELS_with_Incremental_QR_factorization.res"))
ggplot(df, aes(x=global_iteration, y=least_square_time, color=name)) +
geom_line() + ggtitle("Least Square duration (young1c)")
#+end_src
#+results:
[[file:/tmp/babel-1939RTd/figure1939xoE.png]]
*** Influence of the Algorithm
**** run test case
#+begin_src shell
......@@ -120,19 +128,24 @@ mkdir -p ../data/res
**** plot the graphic
#+begin_src R
library(ggplot2)
df <- read.csv2("./build/src/data/res/STD.res")
df <- rbind(df, read.csv2("./build/src/data/res/STD+DR.res"))
df <- rbind(df, read.csv2("./build/src/data/res/IB.res"))
library(latex2exp)
df <- read.csv("./build/src/data/res/STD.res")
df <- rbind(df, read.csv("./build/src/data/res/STD+DR.res"))
df <- rbind(df, read.csv("./build/src/data/res/IB.res"))
name <- df$name
df$name_min <- paste(name, " minRes")
df$name_max <- paste(name, " maxRes")
df$name_min <- paste(name, " min")
df$name_max <- paste(name, " max")
ggplot(df, aes(x=nb_mvp)) +
geom_line(aes(y=maxRes, color=name_max)) +
geom_line(aes(y=minRes, color=name_min)) +
geom_hline(aes(yintercept=1e-7, color="threshold")) +
scale_y_log10() +
scale_y_log10() + ylab(TeX("$\\eta_b$")) +
ggtitle("coneSphere problem with CGS RUHE (STD, STD+DR and IB algorithm)")
#+end_src
#+results:
[[file:/tmp/babel-1939RTd/figure1939XNE.png]]
*** young1c (nrhs=6, m=90, k=5)
**** run test case
#+begin_src shell
......@@ -155,10 +168,10 @@ mkdir -p ../data/res
#+begin_src R
library(ggplot2)
library(latex2exp)
df <- read.csv2("./build/src/data/res/BGMRES.res")
df <- rbind(df, read.csv2("./build/src/data/res/IB-BGMRES.res"))
df <- rbind(df, read.csv2("./build/src/data/res/BGMRES-DR.res"))
df <- rbind(df, read.csv2("./build/src/data/res/IB-BGMRES-DR.res"))
df <- read.csv("./build/src/data/res/BGMRES.res")
df <- rbind(df, read.csv("./build/src/data/res/IB-BGMRES.res"))
df <- rbind(df, read.csv("./build/src/data/res/BGMRES-DR.res"))
df <- rbind(df, read.csv("./build/src/data/res/IB-BGMRES-DR.res"))
ggplot(df, aes(x=nb_mvp)) +
geom_line(aes(y=maxRes, color=name)) +
geom_line(aes(y=minRes, color=name)) +
......@@ -166,6 +179,8 @@ ggplot(df, aes(x=nb_mvp)) +
ggtitle("young1c (nrhs=6, m=90 k=5)") + ylab(TeX("$\\eta_b$ (min, max)"))
#+end_src
#+results:
[[file:/tmp/babel-15182Ttv/figure15182NeX.png]]
*** lightInTissue (nrhs=6, m=90, k=5)
**** run test case
......@@ -182,10 +197,10 @@ mkdir -p ../data/res
#+begin_src R
library(ggplot2)
library(latex2exp)
df <- read.csv2("./build/src/data/res/BGMRES.res")
df <- rbind(df, read.csv2("./build/src/data/res/IB-BGMRES.res"))
df <- rbind(df, read.csv2("./build/src/data/res/BGMRES-DR.res"))
df <- rbind(df, read.csv2("./build/src/data/res/IB-BGMRES-DR.res"))
df <- read.csv("./build/src/data/res/BGMRES.res")
df <- rbind(df, read.csv("./build/src/data/res/IB-BGMRES.res"))
df <- rbind(df, read.csv("./build/src/data/res/BGMRES-DR.res"))
df <- rbind(df, read.csv("./build/src/data/res/IB-BGMRES-DR.res"))
ggplot(df, aes(x=nb_mvp)) +
geom_line(aes(y=maxRes, color=name)) +
geom_line(aes(y=minRes, color=name)) +
......@@ -207,10 +222,10 @@ mkdir -p ../data/res
#+begin_src R
library(ggplot2)
library(latex2exp)
df <- read.csv2("./build/src/data/res/BGMRES-DR(05).res")
df <- rbind(df, read.csv2("./build/src/data/res/BGMRES-DR(10).res"))
df <- rbind(df, read.csv2("./build/src/data/res/BGMRES-DR(15).res"))
df <- rbind(df, read.csv2("./build/src/data/res/BGMRES-DR(20).res"))
df <- read.csv("./build/src/data/res/BGMRES-DR(05).res")
df <- rbind(df, read.csv("./build/src/data/res/BGMRES-DR(10).res"))
df <- rbind(df, read.csv("./build/src/data/res/BGMRES-DR(15).res"))
df <- rbind(df, read.csv("./build/src/data/res/BGMRES-DR(20).res"))
ggplot(df, aes(x=nb_mvp)) +
geom_line(aes(y=maxRes, color=name)) +
geom_line(aes(y=minRes, color=name)) +
......@@ -223,24 +238,27 @@ ggplot(df, aes(x=nb_mvp)) +
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 IBDR \
-r DEFLATED -p 5 -e 1e-4 -u -o "IB-BGMRES-DR"
./testMatrixMarketChoice -f ../data/New/lightINtissue.mtx -k CSC \
-t BLOCK -s CGS \
-n 6 -m 90 -A IBDR -r DEFLATED -p 5 -e 1e-4 \
-u -o "IB-BGMRES-DR"
#+end_src
**** plot
#+begin_src R
library(ggplot2)
library(latex2exp)
df <- read.csv2("./build/src/data/res/IB-BGMRES-DR.res")
#df <- rbind(df, read.csv2("./build/src/data/res/IB-BGMRES.res"))
#df <- rbind(df, read.csv2("./build/src/data/res/BGMRES-DR.res"))
#df <- rbind(df, read.csv2("./build/src/data/res/IB-BGMRES-DR.res"))
ggplot(df, aes(x=nb_mvp)) +
df <- read.csv("./build/src/data/res/IB-BGMRES-DR.res")
ggplot(df, aes(x=global_iteration)) +
geom_line(aes(y=time, color="total")) +
geom_line(aes(y=least_square_time, color="gels")) +
geom_line(aes(y=mvp_spent, color="mvp")) +
geom_line(aes(y=ortho_spent, color="ortho")) +
ggtitle("lightInTissue (nrhs=6, m=90 k=5)") + ylab(TeX("$\\eta_b$ (min, max)"))
geom_line(aes(y=mvp_time, color="mvp")) +
geom_line(aes(y=ortho_time, color="ortho")) +
geom_line(aes(y=ib_time, color="ib")) +
geom_line(aes(y=(ib_time+ortho_time+mvp_time+least_square_time), color="SUM")) +
ylab(TeX("Time (s)")) +
ggtitle("lightInTissue (nrhs=6, m=90 k=5)")
#+end_src
*** qr ib dr
**** run test
#+begin_src shell
......@@ -254,8 +272,8 @@ mkdir -p ../data/res
#+begin_src R
library(ggplot2)
library(latex2exp)
df <- read.csv2("./build/src/data/res/IB-BGMRES-DR.res")
df <- rbind(df, read.csv2("./build/src/data/res/QR-IB-BGMRES-DR.res"))
df <- read.csv("./build/src/data/res/IB-BGMRES-DR.res")
df <- rbind(df, read.csv("./build/src/data/res/QR-IB-BGMRES-DR.res"))
ggplot(df, aes(x=nb_mvp)) +
geom_line(aes(y=maxRes, color=name)) +
geom_line(aes(y=minRes, color=name)) +
......@@ -275,13 +293,19 @@ mkdir -p ../data/res
#+begin_src R
library(ggplot2)
library(latex2exp)
df <- read.csv2("./build/src/data/res/IB-BGMRES-DR.res")
df <- rbind(df, read.csv2("./build/src/data/res/QR-IB-BGMRES-DR.res"))
df <- read.csv("./build/src/data/res/IB-BGMRES-DR.res")
df <- rbind(df, read.csv("./build/src/data/res/QR-IB-BGMRES-DR.res"))
ggplot(df, aes(x=nb_mvp)) +
geom_line(aes(y=least_square_time, color=name)) +
geom_line(aes(y=least_square_time, color=name)) +
ggtitle("young1c (nrhs=100, m=500 k=5, max_mvp=10000)") + ylab(TeX("$\\eta_b$ (min, max)"))
geom_line(aes(y=least_square_time, color=paste(name,"gels"))) +
geom_line(aes(y=facto_time, color=paste(name, "facto"))) +
geom_line(aes(y=facto_time+least_square_time, color=paste(name, "sum"))) +
ggtitle("young1c (nrhs=100, m=500 k=5, max_mvp=10000)") +
ylab(TeX("$\\eta_b$ (min, max)"))
#+end_src
#+results:
[[file:/tmp/babel-15182Ttv/figure15182m4Q.png]]
*** with chameleon
**** run test
#+begin_src shell
......@@ -304,7 +328,7 @@ export STARPU_FXT_PREFIX=${WORKDIR}/build/
#+begin_src R
library(ggplot2)
library(latex2exp)
df <- read.csv2("./build/src/data/res/CHAM-TL.res")
df <- read.csv("./build/src/data/res/CHAM-TL.res")
ggplot(df, aes(x=global_iteration)) +
geom_line(aes(y=100, color="total")) +
geom_line(aes(y=(least_square_time/time)*100, color="Gels")) +
......@@ -327,12 +351,17 @@ cd ${WORKDIR}/build/src/test_cham/
#+begin_src R
library(ggplot2)
library(latex2exp)
df <- read.csv2("./build/src/data/res/QR.res")
df <- rbind(df, read.csv2("./build/src/data/res/CHAMQR.res"))
df <- read.csv("./build/src/data/res/QR.res")
df <- rbind(df, read.csv("./build/src/data/res/CHAMQR.res"))
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("young1c (nrhs=6, m=90 k=5)") + ylab(TeX("$\\eta_b$ (min, max)"))
#+end_src
** Cleanup results directory
#+BEGIN_SRC shell
mkdir -p ${WORKDIR}/build/data/res
cd ${WORKDIR}/build/data/res
rm * -rf
#+END_SRC
......@@ -127,7 +127,8 @@ private:
void do_ortho(Matrix &A, OrthoParam &orthoparm, Block<S> &W)
{
Orthogonalizer ortho{orthoparm};
ortho.run(_hess, _base, W, A, _logger); // Ortho and filling of L (the (IB) Hessenberg)
// Ortho and filling of L (the (IB) Hessenberg)
ortho.run(_hess, _base, W, A, _logger);
}
template<class Matrix, class Block>
......@@ -140,7 +141,6 @@ private:
_solution = Block{_dim, _nbRHS};
_solution.copy(X);
compute_solution(A, _solution);
Block R{_dim, _nbRHS};
compute_real_residual(A, _solution, B, R); // R <- A*sol - B
......@@ -187,21 +187,12 @@ private:
return std::get<2>(MinMaxConv);
}
std::tuple<P,P,bool>
check_least_square_residual(const std::vector<P> &epsilon)
{
_logger.notify_least_square_check_begin();
auto minmaxconv = _hess.check_least_square_residual(epsilon);
_logger.notify_least_square_check_end();
return minmaxconv;
}
public:
template<class RestartParam>
ArnoldiDR(Logger<P> &logger, Restarter<RestartParam, S> &restarter,
int dim, int nbRHS, int max_krylov_space_size):
_logger(logger),
_hess{max_krylov_space_size, nbRHS},
_hess{max_krylov_space_size, nbRHS, _logger},
_solution_computed(false),
_base(restarter.get_base()),
_dim(dim),
......@@ -213,7 +204,8 @@ public:
{
static_assert(
hessenbergXrestarter<HESSENBERG, RestartParam>::value,
"The Hessenberg used is not compatible with the kind of restarting selected"
"The Hessenberg used is not compatible with the kind "
"of restarting selected"
);
}
......@@ -257,7 +249,7 @@ public:
_base.increase(_nbRHS); // add W to base
//_hess.display_hess_bitmap("Hess Before Least Square");
auto MinMaxConv = check_least_square_residual(epsilon);
auto MinMaxConv = _hess.check_least_square_residual(epsilon);
bool proj_convergence = std::get<2>(MinMaxConv);
convergence = check_real_residual_if(proj_convergence, A, B, X, epsilon);
//check_arnoldi_formula(A, _base, _hess, _nbRHS);
......@@ -278,7 +270,9 @@ public:
_base.reset();
} else if (restarter && !convergence) {
Block<S> PR{_hess.get_nb_vect() + _nbRHS, _nbRHS};
_hess.compute_proj_residual(PR);
Block<S> Y = _hess.alloc_least_square_sol();
_hess.solve_least_square(Y);
_hess.compute_proj_residual(PR, Y);
restarter.save_LS(PR);
Block<S> H = _hess.get_hess_block();
restarter.save_hess(H);
......
......@@ -167,15 +167,6 @@ private:
return false;
}
std::tuple<P,P,bool>
check_least_square_residual(const std::vector<P> &epsilon)
{
_logger.notify_least_square_check_begin();
auto minmaxconv = _F.check_least_square_residual(epsilon);
_logger.notify_least_square_check_end();
return minmaxconv;
}
template<class Matrix>
void do_ortho(Matrix &A, OrthoParam &orthoparm)
{
......@@ -336,7 +327,7 @@ private:
Block<S> Y = _F.alloc_least_square_sol();
Block<S> PR{_F.get_nb_hess_line(), _nbRHS};
_F.solve_least_square(Y);
_F.compute_proj_residual(PR);
_F.compute_proj_residual(PR, Y);
Block<S> U;
_nb_direction_kept = PR.decomposition_SVD(U, epsilon, inv_epsilon);
......@@ -366,7 +357,7 @@ public:
ArnoldiIB(Logger<P> &logger, Restarter<RestartParam, S> &restarter,
int dim, int nbRHS, int max_krylov_space_size):
_logger(logger),
_F{max_krylov_space_size, nbRHS},
_F{max_krylov_space_size, nbRHS, _logger},
_solution_computed{false},
_base{restarter.get_base()},
_dim(dim),
......@@ -432,10 +423,10 @@ public:
// _WP.check_ortho_WP();// _WP.get_W().check_ortho();
//_F.display_hess_extended_bitmap();
auto MinMaxConv = check_least_square_residual(epsilon);
auto MinMaxConv = _F.check_least_square_residual(epsilon);
bool proj_convergence = std::get<2>(MinMaxConv);
convergence = check_real_residual_if(proj_convergence, A, X, B, epsilon);
convergence = check_real_residual_if(proj_convergence,
A, X, B, epsilon);
if (!convergence) {
do_R_criterion(proj_convergence, epsilon, inv_epsilon);
}
......
......@@ -145,13 +145,12 @@ private:
Matrix &A, Block &X, const Block &B,
const std::vector<primary_type> &epsilon)
{
if (!projected_residual_converged)
return false;
_solution = Block{_dim, _nbRHS};
_solution.copy(X);
compute_solution(A, _solution);
_solution_computed = true;
if (!projected_residual_converged)
return false;
Block R{_dim, _nbRHS};
compute_real_residual(A, _solution, B, R); // R <- A*sol - B
......@@ -163,6 +162,7 @@ private:
#endif
if (std::get<2>(MinMaxConv)) {
_solution_computed = true;
_cold_restart = false;
return true;
}
......@@ -170,20 +170,12 @@ private:
return false;
}
std::tuple<P,P,bool>
check_least_square_residual(const std::vector<P> &epsilon)
{
_logger.notify_least_square_check_begin();
auto minmaxconv = _F.check_least_square_residual(epsilon);
_logger.notify_least_square_check_end();
return minmaxconv;
}
template<class Matrix>
void do_ortho(Matrix &A, OrthoParam &orthoparm)
{
Orthogonalizer ortho{orthoparm};
ortho.run(_F, _base, _WP, A, _logger); // Ortho and filling of L (the (IB) Hessenberg)
// Ortho and filling of L (the (IB) Hessenberg)
ortho.run(_F, _base, _WP, A, _logger);
}
template<class Matrix, class Block, class Epsilon, class Restarter>
......@@ -306,8 +298,7 @@ private:
return;
}
_ib_update_scheduled = false;
_logger.notify_ib_begin();
Block<S> &U = _U;
assert( U.get_nb_col() == _nb_direction_kept );
......@@ -322,7 +313,8 @@ private:
qr::OutOfPlaceQRFacto(U1_2, W1_W2, R); // Qr facto of bottom block of U
Block<S> W1, W2;
W1 = W1_W2.sub_block(0, 0, _nbRHS, _nb_direction_kept);
W2 = W1_W2.sub_block(0, _nb_direction_kept, _nbRHS, _nb_direction_discarded);
W2 = W1_W2.sub_block(
0, _nb_direction_kept, _nbRHS, _nb_direction_discarded);
// PART 1
Block<S> V = _base.get_W(_nb_direction_kept);
......@@ -337,7 +329,11 @@ private:
_F.update_phi(W1_W2, _nb_direction_kept);
// L_{j+1,} = | \W_1^{H} | * H_j
// Gj = | \W_2^{H} | * H_j
_F.update_bottom_line(W1_W2, _nb_direction_kept);
_F.update_bottom_line(W1_W2);
_logger.notify_ib_end();
// hook that may perform incremental factorization
_F.after_bottom_line_update(_nb_direction_kept);
}
template<class P>
......@@ -346,8 +342,11 @@ private:
const std::vector<P> &inv_epsilon)
{
Block<S> PR{_F.get_nb_vect()+_nbRHS, _nbRHS};
_F.compute_proj_residual(PR);
Block<S> Y = _F.alloc_least_square_sol();
_F.solve_least_square(Y);
_logger.notify_ib_begin();
_F.compute_proj_residual(PR, Y);
Block<S> U;
_nb_direction_kept = PR.decomposition_SVD(U, epsilon, inv_epsilon);
_nb_direction_discarded = _nbRHS - _nb_direction_kept;
......@@ -371,6 +370,7 @@ private:
}
}
_ib_update_scheduled = true;
_logger.notify_ib_end();
}
/*********************************************************/
......@@ -380,7 +380,7 @@ public:
ArnoldiIBDR(Logger<P> &logger, Restarter<RestartParam, S> &restarter,
int dim, int nbRHS, int max_krylov_space_size):
_logger(logger),
_F{max_krylov_space_size, nbRHS},
_F{max_krylov_space_size, nbRHS, _logger},
_solution_computed{false},
_base{restarter.get_base()},
_dim(dim),
......@@ -444,7 +444,7 @@ public:
W_AMV(A, Vj, _WP); // W <- A*M*Vj: perform matrix multiplication
do_ortho(A, orthoparm);
auto MinMaxConv = check_least_square_residual(epsilon);
auto MinMaxConv = _F.check_least_square_residual(epsilon);
bool proj_convergence = std::get<2>(MinMaxConv);
convergence = check_real_residual_if(proj_convergence, A, X, B, epsilon);
if (!convergence) {
......@@ -468,7 +468,9 @@ public:
_base.reset();
} else if (!convergence && restarter) {
Block<S> PR{_F.get_nb_vect()+_nbRHS, _nbRHS};
_F.compute_proj_residual(PR);
Block<S> Y = _F.alloc_least_square_sol();
_F.solve_least_square(Y);
_F.compute_proj_residual(PR, Y);
restarter.save_LS(PR);
Block<S> H = _F.get_hess_block();
restarter.save_hess(H);
......
......@@ -21,6 +21,7 @@ template<class S> class HessChamQR;
#include "fabulous/utils/Arithmetic.hpp"
#include "fabulous/utils/Error.hpp"
#include "fabulous/utils/Chameleon.hpp"
#include "fabulous/utils/Logger.hpp"
#include "fabulous/kernel/chameleon.hpp"
#include "fabulous/data/Block.hpp"
......@@ -35,7 +36,7 @@ inline std::string make_quad_key(int i, int j, int m, int n)
return ss.str();
}
/**
/*!
* \brief Hessenberg for Incremental QR version with chamelon back-end
*
* This implementation use the 'low level' chameleon API (*_Tile_Async)
......@@ -56,44 +57,46 @@ inline std::string make_quad_key(int i, int j, int m, int n)
template<class S>
class HessChamQR
{
public:
using value_type = typename Arithmetik<S>::value_type;
using primary_type = typename Arithmetik<S>::primary_type;
using P = primary_type;
private:
int _max_nb_iter;
int _mb; /**< maximum number of (big) block line (_nbRHSx_nbRHS) in hessenberg */
int _nb; /**< maximum number of (big) block column (_nbRHSx_nbRHS) in hessenberg */
int _nbRHS; /**< number of right-hand-sides; AND size of one (big) block */
int _nb_block_col; /**< number of block columns in current hessenberg */
int _nb_vect; /**< number of columns in current hessenberg */
std::vector<DescPtr> _tau; /**< tau factors for each QR factorization*/
HessenbergTileMapper<S> _hess; /**< the hessenbergs tiles */
ClassicTileMapper<S> _rhs; /**< the (least square) right-hand-side tiles */
SeqPtr _seq; /**< MORSE sequence object */
Logger<P> &_logger;
int _max_nb_iter; /*!< maximum number of iterations */
int _mb; /*!< max number of (big) block line (_nbRHSx_nbRHS) in hessenberg */
int _nb; /*!< max number of (big) block column (_nbRHSx_nbRHS) in hessenberg */
int _nbRHS; /*!< number of right-hand-sides; AND size of one (big) block */
int _nb_block_col; /*!< number of block columns in current hessenberg */
int _nb_vect; /*!< number of columns in current hessenberg */
std::vector<DescPtr> _tau; /*!< tau factors for each QR factorization*/
HessenbergTileMapper<S> _hess; /*!< the hessenbergs tiles */
ClassicTileMapper<S> _rhs; /*!< the (least square) right-hand-side tiles */
SeqPtr _seq; /*!< MORSE sequence object */
#ifdef FABULOUS_DEBUG_MODE
Block<S> _original; /**< original (not factorized) hessenberg
Block<S> _original; /*!< original (not factorized) hessenberg
stored as Lapack matrix for debug purpose
(more specifically, in the
check_arnoldi_formula() function) */
check_arnoldi_formula() function) */
#endif
typedef std::unordered_map<std::string, MorseDesc<S>> SubMap;
SubMap _sub_hess; /**< keep a reference on
SubMap _sub_hess; /*!< keep a reference on
(HESSENBERG) MorseDesc object to keep them alive */
SubMap _sub_rhs;/**< keep a reference on
SubMap _sub_rhs;/*!< keep a reference on
(least square RHS) MorseDesc object to keep them alive */
bool _last_column_factorized; /**< whether last column was factorized */
bool _solution_computed; /**< whether solution was computed(for last column) */
Block<S> _Y; /**< least square solution */
Block<S> _HR; /**< last column of hessenberg to be filled by orthogonalization */
bool _last_column_factorized; /*!< whether last column was factorized */
bool _solution_computed; /*!< whether solution was computed(for last column) */
Block<S> _Y; /*!< least square solution */
Block<S> _HR; /*!< last column of hessenberg to be filled by orthogonalization */
FABULOUS_DELETE_COPY_MOVE(HessChamQR);
public:
using value_type = typename Arithmetik<S>::value_type;
using primary_type = typename Arithmetik<S>::primary_type;
using P = primary_type;
template< class Block >
template<class Block>
void debug_dump_original(Block &Hm, int m, int n)
{
#ifdef FABULOUS_DEBUG_MODE
......@@ -110,7 +113,8 @@ public:
}
private:
MORSE_desc_t *get_sub(TileMapper<S> &parent, SubMap &map, int i, int j, int m, int n)
MORSE_desc_t *get_sub(TileMapper<S> &parent, SubMap &map,
int i, int j, int m, int n)
{
//FABULOUS_DEBUG("i="<<i<<" j="<<j<<" m="<<m<<" n="<<n);
auto idx = make_quad_key(i, j, m, n);
......@@ -151,14 +155,16 @@ private:
return get_sub_rhs(0, 0, k, 1);
}
/**
* \brief compute qr factorization of the last column of the block hessenberg matrix
/*!
* \brief compute qr factorization of the last column of
* the block hessenberg matrix
*/
void factorize_last_column()
{
if (_last_column_factorized)
return;
_logger.notify_facto_begin();
MORSE_enum trans = Arithmetik<S>::mtrans;
/* STEP 1: Loop over each block in last column to
......@@ -195,10 +201,12 @@ private:
FABULOUS_FATAL_ERROR("return value of ormqrf (RHS) is "<<err);
_last_column_factorized = true;
_logger.notify_facto_end();
}
public:
HessChamQR(int max_krylov_space_size, int nbRHS):
HessChamQR(int max_krylov_space_size, int nbRHS, Logger<P> &logger):
_logger{logger},
_max_nb_iter{ (max_krylov_space_size/nbRHS) +1 },
_mb{ _max_nb_iter+1 },
_nb{ _max_nb_iter },
......@@ -221,7 +229,7 @@ public:
_seq.reset(sequence, SeqDeleter{});
}
/**
/*!
* \brief set the \f$\Lambda_1\f$ block (which is nbRHS x nbRHS)
* (RHS of least square problem)
*/
......@@ -266,12 +274,12 @@ public:
if (_solution_computed)
return;
factorize_last_column();
_logger.notify_least_square_begin();
MORSE_desc_t *R = get_full_tr_matrix();
MORSE_desc_t *Lambda = get_full_rhs_for_trsm();
int k = _nb_block_col;
ClassicTileMapper<S> YY{k, 1, _nbRHS}; // tile pool for RHS copy
MorseDesc<S> yy{YY, 0, 0, k, 1}; // actual descriptor instantation
MORSE_desc_t *YYd = yy.get();
......@@ -284,11 +292,11 @@ public:
MORSE_Sequence_Wait(_seq.get());
MORSE_Tile_to_Lapack(YYd, _Y.get_ptr(), _Y.get_leading_dim());
_solution_computed = true;
_logger.notify_least_square_end();
}
void compute_proj_residual(const Block<S> &)
void compute_proj_residual(const Block<S>&, const Block<S>&)
{
FABULOUS_FATAL_ERROR(
"NOT IMPLEMENTED: HessChamQR not compatible with DeflatedRestarting");
......@@ -301,17 +309,16 @@ public:
int k = _nb_block_col;
assert( _nb_vect == k * _nbRHS );
MORSE_desc_t *Lambda = get_sub_rhs(k, 0, 1, 1);
Block<S> L{_nbRHS, _nbRHS};
MORSE_Sequence_Wait(_seq.get());
MORSE_Tile_to_Lapack(Lambda, L.get_ptr(), L.get_leading_dim());
auto MinMaxConv = L.check_precision(epsilon);
return L.check_precision(epsilon);
return MinMaxConv;
}
void notify_orthogonalization_end()
void notify_ortho_end()
{
int k = _nb_block_col;
assert( _nb_vect == k * _nbRHS );
......@@ -321,7 +328,8 @@ public:
#ifdef FABULOUS_DEBUG_MODE
// _HR.display("ortho col _HR: ");
Block<S> HR_original;
HR_original = _original.sub_block(0, _nb_vect-_nbRHS, _nb_vect+_nbRHS, _nbRHS);
HR_original = _original.sub_block(0, _nb_vect-_nbRHS,
_nb_vect+_nbRHS, _nbRHS);
HR_original.copy(_HR);
#endif
MORSE_desc_t *sub = get_sub_hess(0, k-1, k+1, 1);
......
......@@ -9,11 +9,13 @@ template< class > class HessChamTLDR;
#include "fabulous/kernel/chameleon-toplevel.hpp"
#include "fabulous/data/Block.hpp"
#include "fabulous/utils/Logger.hpp"
namespace fabulous {
/**
* \brief The Hessenberg in the standard BGMres case (with Chameleon TOPLEVEL interface)
/*!
* \brief The Hessenberg in the standard BGMres case
* (with Chameleon TOPLEVEL interface)
*
* Format: <br/>
\verbatim
......@@ -23,7 +25,8 @@ namespace fabulous {
| X |
\endverbatim
*
* @warning the whole matrix is allocated, its size is driven by the restart parameter.
* \warning the whole matrix is allocated, its size is driven by
* the restart parameter.
*/
template< class S >
class HessChamTLDR : public Block<S>
......@@ -32,6 +35,7 @@ public:
FABULOUS_INHERITS_BLOCK(S);
private:
Logger<P> &_logger;
int _nbRHS; /*!< number of right hand sides */
int _nb_vect; /*!< number of vector */
int _nb_eigen_pair; /*!< number of eigen pair in deflated_restarting */
......@@ -64,9 +68,9 @@ public:
}
public:
HessChamTLDR() = default;
HessChamTLDR(int max_krylov_space_size, int nbRHS):
HessChamTLDR(int max_krylov_space_size, int nbRHS, Logger<P> &logger):
super{ max_krylov_space_size+nbRHS, max_krylov_space_size},
_logger{logger},
_nbRHS{nbRHS},
_nb_vect{0},
_nb_eigen_pair{0},
......@@ -110,6 +114,7 @@ public:
if (_solution_computed)