Commit 561fe57c authored by MIJIEUX Thomas's avatar MIJIEUX Thomas

begin implementation of GCR

parent 7e4d3390
......@@ -6,6 +6,7 @@
#include <fabulous/arnoldi/ArnoldiIB.hpp>
#include <fabulous/arnoldi/ArnoldiDR.hpp>
#include <fabulous/arnoldi/ArnoldiIBDR.hpp>
#include <fabulous/arnoldi/Arnoldi2.hpp>
#include <fabulous/hessenberg/HessDR.hpp>
#include <fabulous/hessenberg/HessQR.hpp>
......@@ -18,6 +19,7 @@
namespace fabulous {
namespace bgmres {
template<class S> using ARNOLDI_STD = ArnoldiDR<HessDR<S>, S>;
template<class S> using ARNOLDI_QR = ArnoldiDR<HessQR<S>, S>;
template<class S> using ARNOLDI_QRDR = ArnoldiDR<HessQRDR<S>, S>;
......@@ -30,63 +32,75 @@ template<class S> using ARNOLDI_CHAM_QR = ArnoldiDR<HessChamQR<S>, S>;
template<class S> using ARNOLDI_CHAM_TL = ArnoldiDR<HessChamTLDR<S>, S>;
#endif // FABULOUS_USE_CHAMELEON
inline auto arnoldi_std() { return Algorithm<ARNOLDI_STD>{}; }
inline auto arnoldi_qr() { return Algorithm<ARNOLDI_QR>{}; }
inline auto arnoldi_qrdr() { return Algorithm<ARNOLDI_QRDR>{}; }
inline auto arnoldi_ib() { return Algorithm<ARNOLDI_IB>{}; }
inline auto arnoldi_ibdr() { return Algorithm<ARNOLDI_IBDR>{}; }
inline auto arnoldi_qribdr() { return Algorithm<ARNOLDI_QRIBDR>{}; }
inline auto std() { return Algorithm<ARNOLDI_STD>{}; }
inline auto qr() { return Algorithm<ARNOLDI_QR>{}; }
inline auto qrdr() { return Algorithm<ARNOLDI_QRDR>{}; }
inline auto ib() { return Algorithm<ARNOLDI_IB>{}; }
inline auto ibdr() { return Algorithm<ARNOLDI_IBDR>{}; }
inline auto qribdr() { return Algorithm<ARNOLDI_QRIBDR>{}; }
#ifdef FABULOUS_USE_CHAMELEON
inline auto arnoldi_cham_qr() { return Algorithm<ARNOLDI_CHAM_QR>{}; }
inline auto arnoldi_cham_tl() { return Algorithm<ARNOLDI_CHAM_TL>{}; }
inline auto cham_qr() { return Algorithm<ARNOLDI_CHAM_QR>{}; }
inline auto cham_tl() { return Algorithm<ARNOLDI_CHAM_TL>{}; }
#endif // FABULOUS_USE_CHAMELEON
std::ostream& operator<<(std::ostream &o, const Algorithm<ARNOLDI_STD> &)
inline std::ostream& operator<<(std::ostream &o, const Algorithm<ARNOLDI_STD> &)
{
return (o << "");
}
std::ostream& operator<<(std::ostream &o, const Algorithm<ARNOLDI_QR> &)
inline std::ostream& operator<<(std::ostream &o, const Algorithm<ARNOLDI_QR> &)
{
return (o << "QR-");
}
std::ostream& operator<<(std::ostream &o, const Algorithm<ARNOLDI_QRDR> &)
inline std::ostream& operator<<(std::ostream &o, const Algorithm<ARNOLDI_QRDR> &)
{
return (o << "QR-");
}
std::ostream& operator<<(std::ostream &o, const Algorithm<ARNOLDI_IB> &)
inline std::ostream& operator<<(std::ostream &o, const Algorithm<ARNOLDI_IB> &)
{
return (o << "IB-");
}
std::ostream& operator<<(std::ostream &o, const Algorithm<ARNOLDI_IBDR> &)
inline std::ostream& operator<<(std::ostream &o, const Algorithm<ARNOLDI_IBDR> &)
{
return (o << "IB-");
}
std::ostream& operator<<(std::ostream &o, const Algorithm<ARNOLDI_QRIBDR> &)
inline std::ostream& operator<<(std::ostream &o, const Algorithm<ARNOLDI_QRIBDR> &)
{
return (o << "QR-IB-");
}
#ifdef FABULOUS_USE_CHAMELEON
std::ostream& operator<<(std::ostream &o, const Algorithm<ARNOLDI_CHAM_QR> &)
inline std::ostream& operator<<(std::ostream &o, const Algorithm<ARNOLDI_CHAM_QR> &)
{
return (o << "CHAM_QR");
}
std::ostream& operator<<(std::ostream &o, const Algorithm<ARNOLDI_CHAM_TL> &)
inline std::ostream& operator<<(std::ostream &o, const Algorithm<ARNOLDI_CHAM_TL> &)
{
return (o << "CHAM_TL");
}
#endif // FABULOUS_USE_CHAMELEON
}; // end namespace bgmres
namespace bgcr {
inline auto std() { return Algorithm<::fabulous::bgcr::Arnoldi>{}; }
inline std::ostream& operator<<(std::ostream &o, const Algorithm<::fabulous::bgcr::Arnoldi> &)
{
return (o << "");
}
}; // end namespace bgcr
}; // end namespace fabulous
#endif // FABULOUS_ALGO_TYPE_HPP
#ifndef FABULOUS_BGCR_HPP
#define FABULOUS_BGCR_HPP
#include <algorithm>
#include <iostream>
namespace fabulous {
template<class S> class BGCR;
};
#include "fabulous/data/Block.hpp"
#include "fabulous/data/Base.hpp"
#include "fabulous/utils/Utils.hpp"
#include "fabulous/utils/Logger.hpp"
#include "fabulous/orthogonalization/OrthoParam.hpp"
#include "fabulous/restart/Restarter.hpp"
#include "fabulous/algo/Equation.hpp"
#include "fabulous/algo/Parameters.hpp"
namespace fabulous {
namespace bgcr {
/**
* \brief %Block General Minimum Residual (%BGCR) algorithm
*/
template<class S>
class BGCR
{
public:
using value_type = typename Arithmetik<S>::value_type;
using primary_type = typename Arithmetik<S>::primary_type;
private:
using P = primary_type;
Logger<P> _logger;
int _nb_restart;
int _nb_iteration;
int _mvp;
private:
template<class Matrix>
void print_precond(const Matrix &A, std::ostream &o = std::cerr)
{
using fabulous::Color;
o << "Right Pre Conditionner : [";
if (A.useRightPrecond())
o<<Color::green<< "USED";
else
o<<Color::red<< "NOT USED";
o<<Color::reset<<"]\n";
}
void print_iteration_start_info(
int size_to_span, std::ostream &o = std::cerr)
{
o << "################## Restart number: " <<_nb_restart <<" #################\n";
o << "################## SizeToSpan : " <<size_to_span <<" ##############\n";
}
void print_iteration_end_info(
int size_krylov_space, std::ostream &o = std::cerr)
{
o <<"######################################################\n";
o <<" Restart Done : "<< _nb_restart <<"\n";
o <<" Total iterations : "
<< _nb_iteration<<"(+"<<_nb_restart<<")\n";
o <<" Size of Krylov Space : "<< size_krylov_space <<"\n";
o <<" Total Number of Matrix vector Product : "<< _mvp<<"\n";
o <<"######################################################\n";
}
template<class Matrix, class Algo>
void print_start_info(
int dim, int nbRHS, int maxMVP, int max_krylov_space_size,
const OrthoParam &ortho, const std::vector<P> &epsilon,
Matrix &A, Algo &algo,std::ostream &o = std::cerr)
{
Arithmetik<S> arith;
o << "#######################################################\n";
std::stringstream ss;
ss << algo << "BGCR" ;
int l = (55-(ss.str().length()+2))/2;
o << std::string(l, '#') << " " << ss.str() << " " << std::string(l, '#') << "\n";
o << "#######################################################\n";
o << "Dimension of Problem : "<<dim<<"\n";
o << "Number of Right Hand Side : "<<nbRHS<<"\n";
o << "Maximum nb of Mat Vect product: "<<maxMVP<<"\n";
o << "Max Size Krylov space before restart : "<<max_krylov_space_size<<"\n";
o << "Orthogonalization Scheme : "<<ortho<<"\n";
o << "Arithmetic used : "<<arith<<"\n";
if (epsilon.size() == 1)
o << "Targeted Tolerance : "<<epsilon.front()<<"\n";
else
o << "Multi precision is being used\n";
print_precond(A, o);
o << "#######################################################\n";
}
template <class Matrix >
void check_params(Matrix &A, const Block<S> &X, const Block<S> B)
{
FABULOUS_ASSERT( A.size() == X.get_nb_row() );
FABULOUS_ASSERT( A.size() == B.get_nb_row() );
FABULOUS_ASSERT( X.get_nb_col() == B.get_nb_col() );
}
public:
BGCR(bool log_real_residual = false):
_logger{log_real_residual},
_nb_restart{0},
_nb_iteration{0},
_mvp{0}
{
}
void reset()
{
_nb_restart = 0;
_nb_iteration = 0;
_mvp = 0;
_logger.reset();
}
template< class Algo, class Matrix >
int solve( const Matrix &A,
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,
OrthoParam ortho)
{
Block<S> X{dim, nbRHS, ldx, X_};
Block<S> B{dim, nbRHS, ldb, B_};
return solve(
A, X, B, algo, max_mvp, max_space,
epsilon, ortho
);
}
template<class Algo, class Equation>
int solve(Equation &eq, Algo algo, Parameters param, OrthoParam ortho)
{
return solve(eq._A, eq._X, eq._B, algo,
param._max_mvp, param._max_space,
eq._epsilon, ortho);
}
/**
* \brief BGCR with Inexact Breakdowns, restart and preconditionner
*
* Ref: IB-BGCR-DR (Giraud Jing Agullo): p.21
*
* \param[in] A Matrix, need to provide <br/>
* MatBlockVect(), size(), useRightPreCond(), PrecondBlockVect(), DotProduct()
*
* \param[in,out] X Initial guess solution and solution in ouput
* \param[in] B Block containing Right Hand Side
* \param algo instance of Algorithm as returned by one of the
::fabulous::arnoldi_XX() function family
* \param max_mvp Maximum number of Matrix Vector Product.
* \param max_krylov_space_size maximum size of Krylov space
* \param[in] epsilon Target accuracy
* \param ortho OrthoParam
* (through QR factorization) or vector wise arnoldi. <br/>
* (In distributed, only the vector wise will works)
*
* \return Total Number of Matrix vector product done
* (cumulated size of the spanned Krylov Spaces)
*/
template< class Algo, class Matrix >
int solve( Matrix &A, Block<S> &X, Block<S> &B,
Algo algo, const int max_mvp, const int max_krylov_space_size,
const std::vector<P> &epsilon, OrthoParam ortho )
{
check_params(A, X, B);
const int nbRHS = B.get_nb_col();
const int dim = B.get_nb_row();
reset();
print_start_info( dim, nbRHS, max_mvp, max_krylov_space_size,
ortho, epsilon, A, algo );
std::vector<P> normB, inv_normB;
normB = B.compute_norms(A);
inv_normB = array_inverse(normB);
X.scale(inv_normB); B.scale(inv_normB);
Base<S> V{dim, max_krylov_space_size+nbRHS};
Base<S> Z{dim, max_krylov_space_size+nbRHS};
bool convergence = false;
while (!convergence && _mvp < max_mvp) {
//Compute nb of mat vect product to give to Arnoldi procedure
int size_to_span = std::min(max_krylov_space_size, max_mvp-_mvp);
size_to_span = std::max(size_to_span, nbRHS);
print_iteration_start_info(size_to_span);
if (nbRHS > size_to_span) {
break;
}
using Arnoldi = typename Algo::template algo<S>;
Arnoldi arnoldi{_logger, dim, nbRHS, size_to_span, V, Z};
convergence = arnoldi.run( A, X, B, epsilon, ortho );
_nb_iteration += arnoldi.get_nb_iteration();
int mvp = arnoldi.get_nb_mvp();
_mvp += mvp;
if (mvp == 0) {
break;
}
++ _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;
}
Logger<P> &get_logger()
{
return _logger;
}
}; // end class BGCR
}; // end namespace bgcr
}; // end namespace fabulous
#endif // FABULOUS_BGCR_HPP
......@@ -9,13 +9,16 @@ template<class S> class BGMRes;
};
#include "fabulous/data/Block.hpp"
#include "fabulous/data/Base.hpp"
#include "fabulous/utils/Utils.hpp"
#include "fabulous/utils/Logger.hpp"
#include "fabulous/orthogonalization/OrthoParam.hpp"
#include "fabulous/restart/DeflatedRestart.hpp"
#include "fabulous/restart/Restarter.hpp"
#include "fabulous/algo/Equation.hpp"
#include "fabulous/algo/Parameters.hpp"
namespace fabulous {
namespace bgmres {
/**
* \brief %Block General Minimum Residual (%BGMRes) algorithm
......@@ -231,6 +234,7 @@ public:
}; // end class BGMRes
}; // end namespace bgmres
}; // end namespace fabulous
#endif // FABULOUS_BGMRES_HPP
......@@ -2,7 +2,8 @@
#define FABULOUS_EQUATION_HPP
namespace fabulous {
template<class S> class BGMRes;
namespace bgmres { template<class S> class BGMRes; };
namespace bgcr { template<class S> class BGCR; };
template<class Matrix, class S> class Equation;
};
......@@ -21,7 +22,8 @@ public:
using value_type = typename Arithmetik<S>::value_type;
using primary_type = typename Arithmetik<S>::primary_type;
private:
friend class BGMRes<S>;
friend class bgmres::BGMRes<S>;
friend class bgcr::BGCR<S>;
using P = primary_type;
Matrix &_A;
......
......@@ -2,8 +2,8 @@
#define FABULOUS_PARAMETERS_HPP
namespace fabulous {
template<class S> class BGMRes;
namespace bgmres { template<class S> class BGMRes; };
namespace bgcr { template<class S> class BGCR; };
/**
* \brief Hold parameters for BGMRES algorithm
......@@ -11,7 +11,8 @@ template<class S> class BGMRes;
class Parameters
{
private:
template<class S> friend class BGMRes;
template<class S> friend class bgmres::BGMRes;
template<class S> friend class bgcr::BGCR;
int _max_mvp; /*!< maximum number of matrix vector product */
int _max_space; /*!< maximum size of krylov search space */
......
......@@ -4,11 +4,13 @@
#include "fabulous/algo/Equation.hpp"
#include "fabulous/algo/Parameters.hpp"
#include "fabulous/algo/BGMRes.hpp"
#include "fabulous/algo/BGCR.hpp"
#include "fabulous/orthogonalization/OrthoParam.hpp"
#include "fabulous/utils/Logger.hpp"
namespace fabulous {
namespace bgmres {
/**
* \brief main fabulous solver entry point
......@@ -38,6 +40,25 @@ Logger<P> solve(Equation &eq, const Algo &algo, Parameters param,
return solver.get_logger();
}
}; // end namespace bgmres
namespace bgcr {
/**
* \brief main fabulous solver entry point
*/
template<class Algo, class Equation,
class S = typename Equation::value_type,
class P = typename Equation::primary_type >
Logger<P> solve(Equation &eq, const Algo &algo,
Parameters param, OrthoParam ortho)
{
BGCR<S> solver;
solver.solve(eq, algo, param, ortho);
return solver.get_logger();
}
}; // end namespace bcgr
}; // end namespace fabulous
#endif // FABULOUS_SOLVE_HPP
This diff is collapsed.
......@@ -11,10 +11,11 @@ template<class HESSENBERG, class S> class ArnoldiDR;
#include "fabulous/utils/Logger.hpp"
#include "fabulous/utils/Traits.hpp"
#include "fabulous/restart/DeflatedRestart.hpp"
#include "fabulous/restart/Restarter.hpp"
#include "fabulous/orthogonalization/Orthogonalizer.hpp"
namespace fabulous {
namespace bgmres {
/**
* \brief %Arnoldi iterations
......@@ -284,6 +285,7 @@ public:
};
}; // end namespace bgmres
}; // end namespace fabulous
#endif // FABULOUS_ARNOLDI_DR_HPP
......@@ -10,13 +10,14 @@ template<class HESSENBERG, class S> class ArnoldiIB;
#include "fabulous/utils/Utils.hpp"
#include "fabulous/utils/Logger.hpp"
#include "fabulous/utils/Traits.hpp"
#include "fabulous/restart/DeflatedRestart.hpp"
#include "fabulous/restart/Restarter.hpp"
#include "fabulous/kernel/QR.hpp"
#include "fabulous/orthogonalization/Orthogonalizer.hpp"
#include "fabulous/orthogonalization/OrthoParam.hpp"
namespace fabulous {
namespace bgmres {
/**
* \brief %Arnoldi iterations with inexact breakdowns
......@@ -459,6 +460,7 @@ public:
}; // end class ArnoldiIB
}; // end namespace bgmres
}; // end namespace fabulous
#endif // FABULOUS_ARNOLDI_IB_HPP
......@@ -12,10 +12,11 @@ template<class HESSENBERG, class S> class ArnoldiIBDR;
#include "fabulous/utils/Traits.hpp"
#include "fabulous/orthogonalization/OrthoParam.hpp"
#include "fabulous/orthogonalization/Orthogonalizer.hpp"
#include "fabulous/restart/DeflatedRestart.hpp"
#include "fabulous/restart/Restarter.hpp"
#include "fabulous/kernel/QR.hpp"
namespace fabulous {
namespace bgmres {
/**
* \brief %Arnoldi iterations with Inexact Breakdowns and Deflated Restarting
......@@ -491,6 +492,7 @@ public:
}; // end class ArnoldiIBDR
}; // end namespace bgmres
}; // end namespace fabulous
#endif // FABULOUS_ARNOLDI_IB_DR_HPP
......@@ -92,14 +92,14 @@ public:
_last_flops{0},
_m{m}, _n{n}, _ld{ld}, _ptr{ptr}
{
assert( ld >= m );
assert( ld > 0 );
if (!_ptr) {
_sptr.reset(new S[n*ld], std::default_delete<S[]>());
_ld = ld = std::max(_m, 1);
_sptr.reset(new S[_m*_n], std::default_delete<S[]>());
_ptr = _sptr.get();
std::fill(_ptr, _ptr+n*ld, S{0.0});
}
assert( ld >= m );
assert( ld > 0 );
}
Block(int m, int n, int ld, S* ptr, SPtr sptr):
......@@ -108,14 +108,15 @@ public:
_sptr{sptr},
_ptr{ptr}
{
assert( ld >= m );
assert( ld > 0 );
if (!_ptr) {
_sptr.reset(new S[n*ld], std::default_delete<S[]>());
_ld = ld = std::max(_m, 1);
_sptr.reset(new S[_m*_n], std::default_delete<S[]>());
_ptr = _sptr.get();
std::fill(_ptr, _ptr+n*ld, S{0.0});
}
assert( ld >= m );
assert( ld > 0 );
}
Block(int m, int n, int ld, SPtr sptr):
......@@ -124,14 +125,14 @@ public:
_sptr{sptr},
_ptr{_sptr.get()}
{
assert( ld >= m );
assert( ld > 0 );
if (!_ptr) {
_ld = ld = std::max(_m, 1);
_sptr.reset(new S[n*ld], std::default_delete<S[]>());
_ptr = _sptr.get();
std::fill(_ptr, _ptr+n*ld, S{0.0});
}
assert( ld >= m );
assert( ld > 0 );
}
/*!
......@@ -294,8 +295,11 @@ public:
void check_ortho(std::string name = "B")
{
Block A{_n, _n};
lapacke::Tgemm(*this, *this, A, S{1.0}, S{0.0});
lapacke::Tgemm(_n, _n, _m,
get_ptr(), get_leading_dim(),
get_ptr(), get_leading_dim(),
A.get_ptr(), A.get_leading_dim(),
S{1.0}, S{0.0});
for (int k = 0; k < _n; ++k)
A.at(k, k) -= S{1.0};
auto MM = A.get_min_max_norm();
......@@ -306,12 +310,29 @@ public:
}
/*! check that the block is filled with zeros */
void check_null()
bool check_null()
{
for (int j = 0; j < _n; ++j)
for (int i = 0; i < _m; ++i)
assert( at(i, j) == S{0.0} );
for (int j = 0; j < _n; ++j) {
for (int i = 0; i < _m; ++i) {
FABULOUS_ASSERT( at(i, j) == S{0.0} );
if ( at(i, j) == S{0.0} )
return false;
}
}
FABULOUS_DEBUG("check_null OK");
return true;
}
bool check_nan()
{
for (int j = 0; j < _n; ++j) {
for (int i = 0; i < _m; ++i) {
FABULOUS_ASSERT( not fabulous::isnan(at(i, j)) );
if (fabulous::isnan(at(i,j)))
return false;
}
}
return true;
}
/*! \brief compare norm of each vector against values in epsilon
......
......@@ -33,8 +33,11 @@ int64_t InPlaceQRFactoMGS_User(Block<S> &Q, Block<S> &R, const Matrix &A)
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() );
if (R.get_nb_col() != 0) {
assert( R.get_nb_col() == Q.get_nb_col() );
assert( R.get_nb_col() == R.get_nb_row() );
}
int64_t nb_flops = 0;
for (int j = 0; j < N; ++j) {
......@@ -46,14 +49,17 @@ int64_t InPlaceQRFactoMGS_User(Block<S> &Q, Block<S> &R, const Matrix &A)
// Qj = Qj - dot(Qj,Qi)*Qi
lapacke::axpy(M, -dot, Q.get_vect(i), 1, Q.get_vect(j), 1);
nb_flops += lapacke::axpy_flops<S>(M);
R.at(i, j) = dot;
if (R.get_nb_col() != 0)
R.at(i, j) = dot;
}
auto norm = Q.get_norm(j, A);
nb_flops += Q.get_last_flops();
R.at(j, j) = norm;
if (R.get_nb_col() != 0)
R.at(j, j) = norm;
if (norm == 0.0) {
#ifndef FABULOUS_IB_EXACT_BREAKDOWN_CONTINUE
Q.display();
// Q.display();
FABULOUS_NOTE("column [j="<<j<<"]");
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");
......@@ -107,15 +113,25 @@ int64_t InPlaceQRFacto(Block<S> &Q, Block<S> &R)
}
template<class S>
int64_t OutOfPlaceQRFacto(Block<S> &A, Block<S> &Q, Block<S> &R)
int64_t OutOfPlaceQRFacto(Block<S> &B, Block<S> &Q, Block<S> &R)
{
assert( &Q != &A );
assert( Q.get_ptr() != A.get_ptr() );
FABULOUS_ASSERT( &Q != &B );
FABULOUS_ASSERT( Q.get_ptr() != B.get_ptr() );
Q.copy(A);
Q.copy(B);
return InPlaceQRFacto(Q, R);
}
template<class Matrix, class S>
int64_t OutOfPlaceQRFacto(Block<S> &B, Block<S> &Q, Block<S> &R, Matrix &A)
{
FABULOUS_ASSERT( &Q != &B );
FABULOUS_ASSERT( Q.get_ptr() != B.get_ptr() );
Q.copy(B);
return InPlaceQRFactoMGS_User(Q, R, A);
}
}; // end namespace QR
}; // end namespace fabulous
......
......@@ -276,6 +276,41 @@ void dot(int N,
cblas_zdotc_sub(N, X, incX, Y, incY, ret);
}
inline int trsm(char side, int m, int n, float alpha,
float *A, int lda, float *B, int ldb)
{
const CBLAS_SIDE cside = (side == 'R' || side == 'r') ? CblasRight : CblasLeft;
cblas_strsm(CblasColMajor, cside, CblasUpper, CblasNoTrans, CblasNonUnit,
m, n, alpha, A, lda, B, ldb);
return 0;
}
inline int trsm(char side, int m, int n, double alpha,
double *A, int lda, double *B, int ldb)
{
const CBLAS_SIDE cside = (side == 'R' || side == 'r') ? CblasRight : CblasLeft;
cblas_dtrsm(CblasColMajor, cside, CblasUpper, CblasNoTrans, CblasNonUnit,
m, n, alpha, A, lda, B, ldb);
return 0;
}