Commit ad4bd273 authored by MIJIEUX Thomas's avatar MIJIEUX Thomas

Implement IB+DR

parent 88be30e1
......@@ -50,9 +50,9 @@ endif (FABULOUS_LAPACKE_NANCHECK)
if(FABULOUS_DEBUG_MODE)
set( CMAKE_CXX_FLAGS
"${CMAKE_CXX_FLAGS} -DFABULOUS_DEBUG_MODE -ggdb -g3 -O0 -fmax-errors=2" )
"${CMAKE_CXX_FLAGS} -DFABULOUS_DEBUG_MODE -ggdb -g3 -O0 ") #-ferror-limit=2" )
set( CMAKE_C_FLAGS
"${CMAKE_C_FLAGS} -DFABULOUS_DEBUG_MODE -ggdb -g3 -O0 -fmax-errors=2" )
"${CMAKE_C_FLAGS} -DFABULOUS_DEBUG_MODE -ggdb -g3 -O0 ") #-ferror-limit=2" )
else()
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3 -fmax-errors=2")
set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -O3 -fmax-errors=2")
......
......@@ -4,6 +4,7 @@
#+TITLE: FABuLOuS numericals results
* Get Test files
** File list
[[file:./build/src/data/young1c.mtx][young1c.mtx]]
[[file:./build/src/data/matrices_BEM/MatconeSpherePC_MAIN_MAIN_0][MatconeSpherePC_MAIN_MAIN_0]] (binary)
[[file:./build/src/data/matrices_BEM/coneSpherePC_RHS_MAIN.0.res][coneSpherePC_RHS_MAIN.0.res]] (binary)
......@@ -36,8 +37,8 @@ for i in "${files[@]}"; do
done
#+end_src
* Run Examples
Set this variable to your work directory (where you downloaded fabulous):
* Setup workstation
Set this variable to your work directory (where you downloaded fabulous):
#+begin_src sh :session *TEST* :results silent
export WORKDIR=/home/tmijieux/fabulous # for example
#+end_src
......@@ -50,7 +51,7 @@ cd build
cmake .. # -DCHAMELEON_DIR=$(spack location -i chameleon)
make
#+end_src
* Run Examples
** Impact of restart parameter
*** run test case
#+begin_src sh :session *TEST* :results silent
......@@ -80,6 +81,7 @@ ggplot(df, aes(x=nb_mvp, y=maxRes, color=name)) +
*** run test case
#+begin_src sh :session *TEST* :results silent
cd ${WORKDIR}/build/src/test_core/
mkdir -p ../data/res
./testMatrixMarketChoice -t BLOCK -s CGS -m 900 -A STDDR -u -o "Basic_GELS"
./testMatrixMarketChoice -t BLOCK -s CGS -m 900 -A QR -u -o "QR_factorization"
#+end_src
......@@ -120,10 +122,12 @@ ggplot(df, aes(x=nb_mvp)) +
*** run test case
#+begin_src sh :session *TEST* :results silent
cd ${WORKDIR}/build/src/test_core/
mkdir -p ../data/res
#nbRHS=6; maxSpace=90
./testMatrixMarketChoice -n 6 -t BLOCK -s MGS -m 90 -e 1e-6 -u -o "BGMRES"
./testMatrixMarketChoice -n 6 -t BLOCK -s MGS -m 90 -e 1e-6 -A IB -u -o "IB-BGMRES"
./testMatrixMarketChoice -n 6 -t RUHE -s MGS -m 90 -e 1e-6 -r DEFLATED -p 5 -u -o "BGMRES-DR"
./testMatrixMarketChoice -n 6 -t RUHE -s MGS -m 90 -e 1e-6 -A IBDR -r DEFLATED -p 5 -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*
......@@ -132,6 +136,7 @@ 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)) +
......
......@@ -18,12 +18,11 @@
*** DONE Add timer for orthogonalization, least square and matrix vector product
*** DONE Add global timer
*** TODO Improve timer semantics and logs
*** TODO Add flops/s counter
** DONE Improve RESULTS.org
Eventually, anyone must be able to gather all interesting results into RESULTS.org
just by evaluating code blocks from RESULTS.org and/or tangling RESULTS.org
* TODO Implement IB+DR
* DONE Implement IB+DR
* TODO Implement IB+DR+QR
* TODO iterated orthogonalization stop criterion
......
......@@ -5,6 +5,7 @@
#include <fabulous/arnoldi/ArnoldiIB.hpp>
#include <fabulous/arnoldi/ArnoldiDR.hpp>
#include <fabulous/arnoldi/ArnoldiIBDR.hpp>
#include <fabulous/hessenberg/HessDR.hpp>
#include <fabulous/hessenberg/HessQR.hpp>
......@@ -12,6 +13,7 @@
#include <fabulous/hessenberg/HessChamTLDR.hpp>
#include <fabulous/hessenberg/HessChamQR.hpp>
#include <fabulous/hessenberg/HessIB.hpp>
#include <fabulous/hessenberg/HessIBDR.hpp>
namespace fabulous {
......@@ -19,6 +21,7 @@ template<class S> using ARNOLDI_STD = ArnoldiDR<HessDR, S>;
template<class S> using ARNOLDI_QR = ArnoldiDR<HessQR, S>;
template<class S> using ARNOLDI_QRDR = ArnoldiDR<HessQRDR, S>;
template<class S> using ARNOLDI_IB = ArnoldiIB<HessIB, S>;
template<class S> using ARNOLDI_IBDR = ArnoldiIBDR<HessIBDR, S>;
#ifdef FABULOUS_USE_CHAMELEON
template<class S> using ARNOLDI_CHAM_QR = ArnoldiDR<HessChamQR, S>;
......@@ -29,6 +32,7 @@ inline Template<ARNOLDI_STD> arnoldi_std() { return Template<ARNOLDI_STD>{}; }
inline Template<ARNOLDI_QR> arnoldi_qr() { return Template<ARNOLDI_QR>{}; }
inline Template<ARNOLDI_QRDR> arnoldi_qrdr() { return Template<ARNOLDI_QRDR>{}; }
inline Template<ARNOLDI_IB> arnoldi_ib() { return Template<ARNOLDI_IB>{}; }
inline Template<ARNOLDI_IBDR> arnoldi_ibdr() { return Template<ARNOLDI_IBDR>{}; }
#ifdef FABULOUS_USE_CHAMELEON
inline Template<ARNOLDI_CHAM_QR> arnoldi_cham_qr()
......@@ -57,6 +61,11 @@ std::ostream& operator<<(std::ostream &o, const Template<ARNOLDI_IB> &)
return (o << "IB");
}
std::ostream& operator<<(std::ostream &o, const Template<ARNOLDI_IBDR> &)
{
return (o << "IBDR");
}
#ifdef FABULOUS_USE_CHAMELEON
std::ostream& operator<<(std::ostream &o, const Template<ARNOLDI_CHAM_QR> &)
......
......@@ -10,7 +10,7 @@ template< class S > class BGMRes;
#include "fabulous/data/Block.hpp"
#include "fabulous/utils/Utils.hpp"
#include "fabulous/orthogonalization/Arnoldi_Ortho.hpp"
#include "fabulous/orthogonalization/OrthoParam.hpp"
#include "fabulous/restart/DeflatedRestart.hpp"
#include "fabulous/algo/Equation.hpp"
#include "fabulous/algo/Parameters.hpp"
......@@ -67,7 +67,7 @@ private:
template< class Matrix, class RestartParam >
void print_start_info(
int dim, int nbRHS, int maxMVP, int max_krylov_space_size,
const Orthogonalizer &ortho, const RestartParam &restart_param,
const OrthoParam &ortho, const RestartParam &restart_param,
const std::vector<P> &epsilon, Matrix &A, std::ostream &o = std::cout)
{
Arithmetik<S> arith;
......@@ -120,7 +120,7 @@ public:
int dim, int nbRHS, S *B_, int ldb, S *X_, int ldx,
const std::vector<P> &epsilon,
Algo algo, const int max_mvp, const int max_space,
Orthogonalizer ortho, RestartParam restart_param)
OrthoParam ortho, RestartParam restart_param)
{
Block<S> X{dim, nbRHS, ldx, X_};
Block<S> B{dim, nbRHS, ldb, B_};
......@@ -133,7 +133,7 @@ public:
template<class Algo, class Equation, class RestartParam>
int solve(Equation &eq, Algo &algo, Parameters param,
Orthogonalizer ortho, RestartParam restart_param)
OrthoParam ortho, RestartParam restart_param)
{
return solve(eq._A, eq._X, eq._B, algo,
param._max_mvp, param._max_space,
......@@ -154,7 +154,7 @@ public:
* \param max_krylov_space_size maximum size of Krylov space
* \param[in] epsilon Target accuracy
* \param[in] restart_param DeflatedRestart or ClassicRestart instance
* \param ortho Orthogonalizer
* \param ortho OrthoParam
* (through QR factorization) or vector wise arnoldi. <br/>
* (In distributed, only the vector wise will works)
*
......@@ -165,7 +165,7 @@ public:
int solve( Matrix &A, Block<S> &X, Block<S> &B,
Algo&, const int max_mvp, const int max_krylov_space_size,
const std::vector<P> &epsilon,
Orthogonalizer ortho, RestartParam restart_param )
OrthoParam ortho, RestartParam restart_param )
{
check_params(A, X, B);
const int nbRHS = B.get_nb_col();
......
......@@ -5,7 +5,7 @@
#include "fabulous/algo/Parameters.hpp"
#include "fabulous/algo/BGMRes.hpp"
#include "fabulous/orthogonalization/Arnoldi_Ortho.hpp"
#include "fabulous/orthogonalization/OrthoParam.hpp"
#include "fabulous/utils/Logger.hpp"
namespace fabulous {
......@@ -17,7 +17,7 @@ template<class Algo, class Equation, class Restart,
class S = typename Equation::value_type,
class P = typename Equation::primary_type >
inline Logger<P> solve(Equation &eq, const Algo &algo, Parameters param,
Orthogonalizer ortho, Restart restart)
OrthoParam ortho, Restart restart)
{
BGMRes<S> solver;
solver.solve(eq, algo, param, ortho, restart);
......@@ -31,7 +31,7 @@ template<class Algo, class Equation,
class S = typename Equation::value_type,
class P = typename Equation::primary_type >
inline Logger<P> solve(Equation &eq, const Algo &algo, Parameters param,
Orthogonalizer ortho = orthogonalization())
OrthoParam ortho = orthogonalization())
{
BGMRes<S> solver;
solver.solve(eq, algo, param, ortho, classic_restart());
......
......@@ -9,10 +9,10 @@ template<template<class> class HESSENBERG, class S > class ArnoldiDR;
#include "fabulous/data/Block.hpp"
#include "fabulous/utils/Utils.hpp"
#include "fabulous/utils/Logger.hpp"
#include "fabulous/utils/TypeInfo.hpp"
#include "fabulous/utils/Traits.hpp"
#include "fabulous/restart/DeflatedRestart.hpp"
#include "fabulous/orthogonalization/Arnoldi_Ortho.hpp"
#include "fabulous/orthogonalization/Orthogonalizer.hpp"
namespace fabulous {
......@@ -205,28 +205,26 @@ public:
*/
template<class Matrix, class Restarter>
bool run( const Matrix &A, Block<S> &X, const Block<S> &B,
const int max_krylov_space_size,
const std::vector<P> &epsilon,
Restarter &restarter,
Orthogonalizer &ortho )
const int max_krylov_space_size, const std::vector<P> &epsilon,
Restarter &restarter, OrthoParam &orthoparm )
{
print_header(max_krylov_space_size);
Block<S> R1{}; // R1 = RHS of projected problem
std::tuple<P,P,bool> MinMaxConv;
_logger.notify_iteration_begin(_iteration_count, _size_krylov_space);
if (restarter && _base.get_nb_block() > 0) {
restarter.run(R1, _hess);
MinMaxConv = restarter.run(_hess, epsilon);
//check_arnoldi_formula(A, _base, _hess);
} else { //No DR, or first run of arnoldi
_base.reset();
Block<S> R0 = _base.get_W(_nbRHS);
compute_real_residual(A, X, B, R0);
_base.increase(_nbRHS);
R1.realloc(_nbRHS, _nbRHS);
Block<S> R1{_nbRHS, _nbRHS}; // R1 = RHS of projected problem
A.QRFacto(R0, R1); /* auto minmaxR0 = R0.get_min_max_norm(A); */
_hess.set_rhs(R1);
MinMaxConv = R1.check_precision(epsilon);
}
auto MinMaxConv = R1.check_precision(epsilon);
bool convergence = std::get<2>(MinMaxConv);
_size_krylov_space = _base.get_nb_vect();
_logger.notify_iteration_end(
......@@ -241,6 +239,8 @@ public:
Block<S> W = _base.get_W(_nbRHS); // W (candidate)
W_AMV(A, Vj, W); // W := A*M*V_j
Orthogonalizer ortho{orthoparm};
_logger.notify_ortho_begin();
ortho.run(_hess, _base, W, A);
_logger.notify_ortho_end();
......@@ -252,6 +252,9 @@ public:
auto MinMaxConv = 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);
_logger.notify_iteration_end(
_iteration_count, _size_krylov_space, _nb_mvp, MinMaxConv);
_logger.log_real_residual(A, B, X, _base, _hess);
......
......@@ -9,7 +9,11 @@ template<template<class> class HESSENBERG, class S> class ArnoldiIB;
#include "fabulous/data/BlockWP.hpp"
#include "fabulous/utils/Utils.hpp"
#include "fabulous/utils/Logger.hpp"
#include "fabulous/utils/TypeInfo.hpp"
#include "fabulous/utils/Traits.hpp"
#include "fabulous/restart/DeflatedRestart.hpp"
#include "fabulous/orthogonalization/Orthogonalizer.hpp"
#include "fabulous/orthogonalization/OrthoParam.hpp"
namespace fabulous {
......@@ -390,8 +394,7 @@ public:
bool run( Matrix &A, Block<S> &X, Block<S> &B,
const int max_krylov_space_size,
const std::vector<P> &epsilon,
Restarter&,
Orthogonalizer &ortho)
Restarter&, OrthoParam &orthoparm)
{
bool convergence = false;
print_header(max_krylov_space_size);
......@@ -433,6 +436,8 @@ public:
W_AMV(A, Vj, WP); // W <- A*M*Vj: perform matrix multiplication
// Ortho and filling of L (the (IB) Hessenberg)
Orthogonalizer ortho{orthoparm};
_logger.notify_ortho_begin();
ortho.run(_F, _base, WP, A);
_logger.notify_ortho_end();
......
This diff is collapsed.
......@@ -257,6 +257,13 @@ public:
std::copy(o.get_vect(j), o.get_vect(j)+M, get_vect(j));
}
Block copy()
{
Block cpy{_m, _n};
cpy.copy(*this);
return cpy;
}
/* \brief pretty print the block values in a tabular format
*/
void display(std::string name ="", std::ostream &o = std::cout) const
......
......@@ -22,6 +22,12 @@ private:
public:
FABULOUS_INHERITS_BLOCK(S);
BlockWP():
super{},
_cursor(0)
{
}
BlockWP(int dim, int nbRHS):
super{dim, nbRHS},
_cursor{0}
......
......@@ -90,7 +90,7 @@ public:
Block<S> get_H1new()
{
assert( _nb_eigen_pair == _nb_vect );
assert( _nb_vect == _nb_eigen_pair );
return this->sub_block(0, 0, _nb_eigen_pair+_nbRHS, _nb_eigen_pair);
}
......
......@@ -281,7 +281,7 @@ public:
get_ptr(), get_leading_dim(),
_YY.get_ptr(), _YY.get_leading_dim(),
PR.get_ptr(), PR.get_nb_row(),
S{1.0}, S{-1.0}
S{-1.0}, S{1.0}
);
}
......
......@@ -30,12 +30,14 @@ private:
int _nb_vect; /*!< number of vector inside the hessenberg */
int _nb_eigen_pair; /*!< number of eigen pair in deflated_restarting */
bool _restarted;
std::vector<int> _block_size; /*!< size of each block */
std::vector<int> _block_size_sum; /*!< sum of size of each block */
std::shared_ptr<Block<S>> _phi; /*!< needed to compute the rhs for Least square.
Ref: IB BGMRes Dr : Page 10. */
Block<S> _lambda1; /*!< initial value of RHS (projected problem) */
Block<S> _Lambda; /*!< rhs of least square problem */
bool _solution_computed; /*!< whether solution (for current size) was computed */
Block<S> _YY; /*!< least square solution */
......@@ -59,7 +61,10 @@ public:
_nb_block_col{0},
_nb_vect{0},
_nb_eigen_pair{0},
_restarted{false},
_phi(nullptr),
_lambda1{},
_Lambda{},
_solution_computed{false},
_YY{}
{
......@@ -77,9 +82,9 @@ public:
{
//Check
if (lambda.get_nb_col() != _nbRHS)
std::cout<<"Lambda nb_col is wrong while initiating lambda in Hess\n";
if (lambda.get_nb_row() != _nbRHS)
std::cout<<"Lambda nb_row is wrong while initiating lambda in Hess\n";
FABULOUS_FATAL_ERROR("Lambda nb_col is wrong while initiating lambda in Hess\n");
if (lambda.get_nb_row() != _nbRHS + _nb_eigen_pair)
FABULOUS_FATAL_ERROR("Lambda nb_row is wrong while initiating lambda in Hess\n");
_lambda1 = lambda;
}
......@@ -100,6 +105,26 @@ public:
_phi->at(i, i) = S{1.0};
}
/**
* \brief Initialize \f$ \Phi \f$ (refer to Annex, Eq 2.5)
* \param p1 number of direction kept in the inexact breakdown on R0
*/
void init_phi_restart(int p1)
{
assert(p1 == _nb_eigen_pair);
assert( _lambda1.get_nb_row() == p1+_nbRHS );
_restarted = true;
if (_phi)
FABULOUS_FATAL_ERROR("THIS FUNCTION CANNOT BE CALLED TWICE!");
int N = _lambda1.get_nb_col();
_phi.reset(new Block<S>{N + p1, N});
for (int i = 0; i < N; ++i)
_phi->at(p1+i, i) = S{1.0};
}
/**
* \brief Update \f$\Phi\f$ if IB happened on R0
*
......@@ -115,21 +140,21 @@ public:
/* If _phi has been set, it means that an I.B. happened
* on R0, and we have a modified RHS for least square
* that depends on _phi */
int N = _nb_vect + _nbRHS;
assert( _phi->get_nb_row() == N );
int nj = _block_size_sum.back();
int M = _phi->get_nb_row() + p_jplus1;
Block<S> *new_phi = new Block<S>{M, _nbRHS};
Block<S> *new_phi = new Block<S>{N+p_jplus1, _nbRHS};
// Init first n_j rows
Block<S> subphi = new_phi->sub_block(0, 0, nj, _nbRHS);
Block<S> subphi = new_phi->sub_block(0, 0, _nb_vect, _nbRHS);
subphi.copy(*_phi);
// Compute newphi(nj+1:nj+p,:) := W1_W2^{H} * _phi(nj+1:nj+p,:)
LapackKernI::Tgemm(
_nbRHS, _nbRHS, _nbRHS,
W1_W2.get_ptr(), W1_W2.get_leading_dim(),
_phi->get_ptr(nj, 0), _phi->get_leading_dim(),
new_phi->get_ptr(nj, 0), new_phi->get_leading_dim(),
_phi->get_ptr(_nb_vect, 0), _phi->get_leading_dim(),
new_phi->get_ptr(_nb_vect, 0), new_phi->get_leading_dim(),
S{1.0}, S{0.0}
);
_phi.reset(new_phi);
......@@ -147,13 +172,40 @@ public:
return;
}
LapackKernI::gemm( // RHS <- Phi * Lambda1
_phi->get_nb_row(), _nbRHS, _nbRHS,
_phi->get_ptr(), _phi->get_leading_dim(),
_lambda1.get_ptr(), _lambda1.get_leading_dim(),
Lambda.get_ptr(), Lambda.get_leading_dim(),
S{1.0}, S{0.0}
);
if (_restarted) {
assert( Lambda.get_nb_row() == _nb_vect + _nbRHS );
assert( Lambda.get_nb_row() == _phi->get_nb_row() );
int p1 = _nb_eigen_pair;
int p = _nbRHS;
// restart update method :: Lambda1: (k+p, p)
Block<S> L1sub = _lambda1.sub_block(0, 0, p1, p);
Lambda.copy(L1sub);
LapackKernI::gemm(
p1, p, p,
_phi->get_ptr(), _phi->get_leading_dim(),
_lambda1.get_ptr()+p1, _lambda1.get_leading_dim(),
Lambda.get_ptr(), Lambda.get_leading_dim(),
S{1.0}, S{1.0}
);
LapackKernI::gemm(
Lambda.get_nb_row()-p1, p, p,
_phi->get_ptr()+p1, _phi->get_leading_dim(),
_lambda1.get_ptr()+p1, _lambda1.get_leading_dim(),
Lambda.get_ptr()+p1, Lambda.get_leading_dim(),
S{1.0}, S{0.0}
);
} else {
// breakdown on R0 restart method :: Lambda1: (p, p)
LapackKernI::gemm( // RHS <- Phi * Lambda1
Lambda.get_nb_row(), _nbRHS, _nbRHS,
_phi->get_ptr(), _phi->get_leading_dim(),
_lambda1.get_ptr(), _lambda1.get_leading_dim(),
Lambda.get_ptr(), Lambda.get_leading_dim(),
S{1.0}, S{0.0}
);
}
}
/**
......@@ -190,6 +242,13 @@ public:
return (_block_size_sum.back() + block_size <= _max_krylov_space_size);
}
void notify_orthogonalization_end() const {/*nop*/}
Block<S> get_hess_block()
{
return this->sub_block(0, 0, _nb_vect+_nbRHS, _nb_vect);
}
/**
* \brief increase the size of the Hessenberg
*
......@@ -284,7 +343,7 @@ public:
get_ptr(), get_leading_dim(),
_YY.get_ptr(), _YY.get_leading_dim(),
PR.get_ptr(), PR.get_nb_row(),
S{1.0}, S{-1.0}
S{-1.0}, S{1.0}
);
}
......@@ -358,6 +417,12 @@ public:
_nb_vect+nb_discarded_direction, offset, block_size, block_size);
}
Block<S> get_F1new()
{
assert( _nb_vect == _nb_eigen_pair );
return this->sub_block(0, 0, _nb_eigen_pair+_nbRHS, _nb_eigen_pair);
}
void reserve_DR(int nb_eigen_pair)
{
// CHECK this is called just after hessenberg initialization:
......@@ -369,8 +434,6 @@ public:
_nb_eigen_pair = nb_eigen_pair;
increase(nb_eigen_pair);
}
void notify_orthogonalization_end() const {/*nop*/}
};
}; // namespace fabulous
......
......@@ -147,48 +147,43 @@ FABULOUS_ARITHMETIC_LIST(FABULOUS_SPECIALIZE_MQR);
// Real versions, with a tricks to split the complex vector into 2 real's one
template<>
inline int LapackKernI::ggev(int n, float *A, int lda,
inline int LapackKernI::ggev(int n,
float *A, int lda,
float *B, int ldb,
std::complex<float> *alpha,
float *beta,
float *vl, int ldvl,
float *vr, int ldvr)
{
// Create temporary real vectors to hold real part and imag part of alpha
float *alphar = new float[n];
float *alphai = new float[n];
memset(alphar,0,sizeof(float)*n);
memset(alphai,0,sizeof(float)*n);
// temporary real vectors to hold real part and imag part of alpha
std::vector<float> alphar, alphai;
alphar.resize(n, 0.0);
alphai.resize(n, 0.0);
int out = LAPACKE_sggev(LAPACK_COL_MAJOR,'N','V',
n, A, lda,
B, ldb,
alphar,alphai,
alphar.data(), alphai.data(),
beta,
vl,ldvl,
vr,ldvr);
std::complex<float> J{0,1};
for (int i=0; i<n; ++i) // Copy back inside alpha
alpha[i] = alphar[i] + J * alphai[i];
delete [] alphar;
delete [] alphai;
return out;
}
// Real versions, with a tricks to split the complex vector into 2 real's one
template<>
inline int LapackKernI::ggev(int n, double *A, int lda,
inline int LapackKernI::ggev(int n,
double *A, int lda,
double *B, int ldb,
std::complex<double> *alpha,
double *beta,
double *vl, int ldvl ,
double *vr, int ldvr)
{
//Create temporary real vectors to hold real part and imag part of
//alpha
// temporary real vectors to hold real part and imag part of alpha
std::vector<double> alphar, alphai;
alphar.resize(n, 0.0);
alphai.resize(n, 0.0);
......@@ -210,12 +205,13 @@ inline int LapackKernI::ggev(int n, double *A, int lda,
}
template<>
inline int LapackKernI::ggev(int n, std::complex<float> * a,int lda,
std::complex<float> * b, int ldb,
std::complex<float> * alpha,
std::complex<float> * beta,
std::complex<float> * vl,int ldvl ,
std::complex<float> * vr, int ldvr)
inline int LapackKernI::ggev(int n,
std::complex<float> *a, int lda,
std::complex<float> *b, int ldb,
std::complex<float> *alpha,
std::complex<float> *beta,
std::complex<float> *vl, int ldvl,
std::complex<float> *vr, int ldvr)
{
return LAPACKE_cggev(LAPACK_COL_MAJOR,'N','V',
n,a,lda,
......@@ -227,11 +223,12 @@ inline int LapackKernI::ggev(int n, std::complex<float> * a,int lda,
template<>
inline int LapackKernI::ggev(int n, std::complex<double> * a,int lda,
inline int LapackKernI::ggev(int n,
std::complex<double> *a, int lda,
std::complex<double> *b, int ldb,
std::complex<double> *alpha,
std::complex<double> *beta,
std::complex<double> *vl,int ldvl ,
std::complex<double> *vl, int ldvl,
std::complex<double> *vr, int ldvr)
{
return LAPACKE_zggev(LAPACK_COL_MAJOR,'N','V',
......
#ifndef ARNOLDI_ORTHO_BLOCK_H
#define ARNOLDI_ORTHO_BLOCK_H
#ifndef ARNOLDI_ORTHO_BLOCK_HPP
#define ARNOLDI_ORTHO_BLOCK_HPP
#ifndef FABULOUS_ORTHO_IMPL_INCLUDE
#error "Never include this header directly; include 'Arnoldi_Ortho.hpp' instead"
#endif
#include "./Arnoldi_Ortho_BLOCK_common.hpp"
#include "fabulous/utils/Traits.hpp"
#include "fabulous/data/Base.hpp"
#include "fabulous/data/Block.hpp"
#include "fabulous/orthogonalization/BLOCK_common.hpp"
namespace fabulous {
/**
* \brief Orthogonalization methods for BLOCK variants WITHOUT Inexact Breakdown
*/
template<class HESS>
class OrthogonalizerBlock : public OrthogonalizerBlockCOMMON
template< class HESS >
class OrthogonalizerBlockSTD : public OrthogonalizerBlockCOMMON
{
private:
friend class Orthogonalizer;
OrthogonalizerBlock(const Orthogonalizer &param):
OrthogonalizerBlockSTD(const OrthoParam &param):
OrthogonalizerBlockCOMMON{param}
{
}
......@@ -46,4 +45,4 @@ private:
}; // namespace fabulous
#endif // ARNOLDI_ORTHO_BLOCK_H
#endif // ARNOLDI_ORTHO_BLOCK_HPP
#ifndef ARNOLDI_ORTHO_BLOCK_IB_H
#define ARNOLDI_ORTHO_BLOCK_IB_H
#ifndef ARNOLDI_ORTHO_BLOCK_IB_HPP
#define ARNOLDI_ORTHO_BLOCK_IB_HPP
#ifndef FABULOUS_ORTHO_IMPL_INCLUDE
#error "Never include this header directly; include 'Arnoldi_Ortho.hpp' instead"
#endif
#include "./Arnoldi_Ortho_BLOCK_common.hpp"
#include "fabulous/utils/Traits.hpp"
#include "fabulous/data/Base.hpp"
#include "fabulous/data/BlockWP.hpp"
#include "fabulous/orthogonalization/BLOCK_common.hpp"
namespace fabulous {
/**
* \brief Orthogonalization methods for BLOCK variants WITH Inexact Breakdown
*/
template<class S>
class OrthogonalizerBlock<HessIB<S>> : public OrthogonalizerBlockCOMMON
template< class HESS >
class OrthogonalizerBlockIB : public OrthogonalizerBlockCOMMON
{
private:
friend class Orthogonalizer;
OrthogonalizerBlock(const Orthogonalizer &param):
OrthogonalizerBlockIB(const OrthoParam &param):
OrthogonalizerBlockCOMMON{param}
{
}
......@@ -31,7 +30,8 @@ private:
* \param[in,out] WP candidate block in Inexact breakdown case
* \param ortho_scheme the orthogonalization scheme,
*/
void run(HessIB<S> &L, Base<S> &base, BlockWP<S> &WP)
template<class S>
void run(HESS &L, Base<S> &base, BlockWP<S> &WP)
{
L.increase(WP.get_size_W());
......@@ -48,4 +48,4 @@ private: