36 static constexpr char TRAN =
'N';
38 int solve_trs(Mat<T>&, Mat<T>&&);
46 FullMat(
const uword in_rows,
const uword in_cols)
47 :
DenseMat<
T>(in_rows, in_cols, in_rows * in_cols) {}
49 unique_ptr<MetaMat<T>>
make_copy()
override {
return std::make_unique<FullMat>(*
this); }
57 T operator()(
const uword in_row,
const uword in_col)
const override {
return this->
memory[in_row + in_col * this->
n_rows]; }
59 T&
at(
const uword in_row,
const uword in_col)
override {
64 Mat<T>
operator*(
const Mat<T>&)
const override;
68 Mat<T> C(arma::size(B));
70 const auto M =
static_cast<int>(this->n_rows);
71 const auto N =
static_cast<int>(this->n_cols);
73 T ALPHA =
T(1), BETA =
T(0);
76 constexpr auto INC = 1;
78 if constexpr(std::is_same_v<T, float>) {
80 arma_fortran(arma_sgemv)(&TRAN, &
M, &
N, (
E*)&ALPHA, (
E*)this->memptr(), &
M, (
E*)B.memptr(), &INC, (
E*)&BETA, (
E*)C.memptr(), &INC);
84 arma_fortran(arma_dgemv)(&TRAN, &
M, &
N, (
E*)&ALPHA, (
E*)this->memptr(), &
M, (
E*)B.memptr(), &INC, (
E*)&BETA, (
E*)C.memptr(), &INC);
88 const auto K =
static_cast<int>(B.n_cols);
90 if constexpr(std::is_same_v<T, float>) {
92 arma_fortran(arma_sgemm)(&TRAN, &TRAN, &
M, &
K, &
N, (
E*)&ALPHA, (
E*)this->memptr(), &
M, (
E*)B.memptr(), &
N, (
E*)&BETA, (
E*)C.memptr(), &
M);
96 arma_fortran(arma_dgemm)(&TRAN, &TRAN, &
M, &
K, &
N, (
E*)&ALPHA, (
E*)this->memptr(), &
M, (
E*)B.memptr(), &
N, (
E*)&BETA, (
E*)C.memptr(), &
M);
104 if(this->factored)
return this->solve_trs(X, std::forward<Mat<T>>(B));
106 suanpan_assert([&] {
if(this->n_rows != this->n_cols)
throw invalid_argument(
"requires a square matrix"); });
110 auto N =
static_cast<int>(this->n_rows);
111 const auto NRHS =
static_cast<int>(B.n_cols);
112 const auto LDB =
static_cast<int>(B.n_rows);
113 this->pivot.zeros(
N);
114 this->factored =
true;
116 if constexpr(std::is_same_v<T, float>) {
118 arma_fortran(arma_sgesv)(&
N, &NRHS, (
E*)this->memptr(), &
N, this->pivot.memptr(), (
E*)B.memptr(), &LDB, &INFO);
123 arma_fortran(arma_dgesv)(&
N, &NRHS, (
E*)this->memptr(), &
N, this->pivot.memptr(), (
E*)B.memptr(), &LDB, &INFO);
127 this->s_memory = this->to_float();
128 arma_fortran(arma_sgetrf)(&
N, &
N, this->s_memory.memptr(), &
N, this->pivot.memptr(), &INFO);
129 if(0 == INFO) INFO = this->solve_trs(X, std::forward<Mat<T>>(B));
133 suanpan_error(
"Error code {} received, the matrix is probably singular.\n", INFO);
141 const auto N =
static_cast<int>(this->n_rows);
142 const auto NRHS =
static_cast<int>(B.n_cols);
143 const auto LDB =
static_cast<int>(B.n_rows);
145 if constexpr(std::is_same_v<T, float>) {
147 arma_fortran(arma_sgetrs)(&TRAN, &
N, &NRHS, (
E*)this->memptr(), &
N, this->pivot.memptr(), (
E*)B.memptr(), &LDB, &INFO);
152 arma_fortran(arma_dgetrs)(&TRAN, &
N, &NRHS, (
E*)this->memptr(), &
N, this->pivot.memptr(), (
E*)B.memptr(), &LDB, &INFO);
156 this->mixed_trs(X, std::forward<Mat<T>>(B), [&](fmat& residual) {
157 arma_fortran(arma_sgetrs)(&TRAN, &
N, &NRHS, this->s_memory.memptr(), &
N, this->pivot.memptr(), residual.memptr(), &LDB, &INFO);
A DenseMat class that holds matrices.
Definition: DenseMat.hpp:39
std::unique_ptr< T[]> memory
Definition: DenseMat.hpp:48
A FullMat class that holds matrices.
Definition: FullMat.hpp:35
FullMat(const uword in_rows, const uword in_cols)
Definition: FullMat.hpp:46
T operator()(const uword in_row, const uword in_col) const override
Access element (read-only), returns zero if out-of-bound.
Definition: FullMat.hpp:57
T & at(const uword in_row, const uword in_col) override
Access element with bound check.
Definition: FullMat.hpp:59
void nullify(const uword K) override
Definition: FullMat.hpp:51
unique_ptr< MetaMat< T > > make_copy() override
Definition: FullMat.hpp:49
void for_each(const IT start, const IT end, F &&FN)
Definition: utility.h:28
void suanpan_assert(const std::function< void()> &F)
Definition: suanPan.h:296
#define suanpan_error(...)
Definition: suanPan.h:309