suanPan
BandMat.hpp
Go to the documentation of this file.
1/*******************************************************************************
2 * Copyright (C) 2017-2024 Theodore Chang
3 *
4 * This program is free software: you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation, either version 3 of the License, or
7 * (at your option) any later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with this program. If not, see <http://www.gnu.org/licenses/>.
16 ******************************************************************************/
29// ReSharper disable CppCStyleCast
30#ifndef BANDMAT_HPP
31#define BANDMAT_HPP
32
33#include "DenseMat.hpp"
34
35template<sp_d T> class BandMat : public DenseMat<T> {
36 static constexpr char TRAN = 'N';
37
38 static T bin;
39
40 const uword s_band;
41
42 int solve_trs(Mat<T>&, Mat<T>&&);
43
44protected:
45 const uword m_rows; // memory block layout
46
47 const uword l_band;
48 const uword u_band;
49
51
52 int direct_solve(Mat<T>&, Mat<T>&&) override;
53
54public:
55 BandMat(const uword in_size, const uword in_l, const uword in_u)
56 : DenseMat<T>(in_size, in_size, (2 * in_l + in_u + 1) * in_size)
57 , s_band(in_l + in_u)
58 , m_rows(2 * in_l + in_u + 1)
59 , l_band(in_l)
60 , u_band(in_u) {
61 if(m_rows >= in_size)
62 suanpan_warning("The storage requirement for the banded matrix is larger than that of a full matrix, consider using a full/sparse matrix instead.\n");
63 }
64
65 unique_ptr<MetaMat<T>> make_copy() override { return std::make_unique<BandMat>(*this); }
66
67 void nullify(const uword K) override {
68 this->factored = false;
69 suanpan::for_each(std::max(K, u_band) - u_band, std::min(this->n_rows, K + l_band + 1), [&](const uword I) { this->memory[I + s_band + K * (m_rows - 1)] = T(0); });
70 suanpan::for_each(std::max(K, l_band) - l_band, std::min(this->n_cols, K + u_band + 1), [&](const uword I) { this->memory[K + s_band + I * (m_rows - 1)] = T(0); });
71 }
72
73 T operator()(const uword in_row, const uword in_col) const override {
74 if(in_row > in_col + l_band || in_row + u_band < in_col) [[unlikely]] return bin = T(0);
75 return this->memory[in_row + s_band + in_col * (m_rows - 1)];
76 }
77
78 T& unsafe_at(const uword in_row, const uword in_col) override {
79 this->factored = false;
80 return this->memory[in_row + s_band + in_col * (m_rows - 1)];
81 }
82
83 T& at(const uword in_row, const uword in_col) override {
84 if(in_row > in_col + l_band || in_row + u_band < in_col) [[unlikely]] return bin = T(0);
85 return this->unsafe_at(in_row, in_col);
86 }
87
88 Mat<T> operator*(const Mat<T>&) const override;
89};
90
91template<sp_d T> T BandMat<T>::bin = T(0);
92
93template<sp_d T> Mat<T> BandMat<T>::operator*(const Mat<T>& X) const {
94 Mat<T> Y(arma::size(X));
95
96 const auto M = static_cast<int>(this->n_rows);
97 const auto N = static_cast<int>(this->n_cols);
98 const auto KL = static_cast<int>(l_band);
99 const auto KU = static_cast<int>(u_band);
100 const auto LDA = static_cast<int>(m_rows);
101 constexpr auto INC = 1;
102 T ALPHA = T(1);
103 T BETA = T(0);
104
105 if constexpr(std::is_same_v<T, float>) {
106 using E = float;
107 suanpan::for_each(X.n_cols, [&](const uword I) { arma_fortran(arma_sgbmv)(&TRAN, &M, &N, &KL, &KU, (E*)&ALPHA, (E*)(this->memptr() + l_band), &LDA, (E*)X.colptr(I), &INC, (E*)&BETA, (E*)Y.colptr(I), &INC); });
108 }
109 else {
110 using E = double;
111 suanpan::for_each(X.n_cols, [&](const uword I) { arma_fortran(arma_dgbmv)(&TRAN, &M, &N, &KL, &KU, (E*)&ALPHA, (E*)(this->memptr() + l_band), &LDA, (E*)X.colptr(I), &INC, (E*)&BETA, (E*)Y.colptr(I), &INC); });
112 }
113
114 return Y;
115}
116
117template<sp_d T> int BandMat<T>::direct_solve(Mat<T>& X, Mat<T>&& B) {
118 if(this->factored) return this->solve_trs(X, std::forward<Mat<T>>(B));
119
120 suanpan_assert([&] { if(this->n_rows != this->n_cols) throw invalid_argument("requires a square matrix"); });
121
122 auto INFO = 0;
123
124 auto N = static_cast<int>(this->n_rows);
125 const auto KL = static_cast<int>(l_band);
126 const auto KU = static_cast<int>(u_band);
127 const auto NRHS = static_cast<int>(B.n_cols);
128 const auto LDAB = static_cast<int>(m_rows);
129 const auto LDB = static_cast<int>(B.n_rows);
130 this->pivot.zeros(N);
131 this->factored = true;
132
133 if constexpr(std::is_same_v<T, float>) {
134 using E = float;
135 arma_fortran(arma_sgbsv)(&N, &KL, &KU, &NRHS, (E*)this->memptr(), &LDAB, this->pivot.memptr(), (E*)B.memptr(), &LDB, &INFO);
136 X = std::move(B);
137 }
138 else if(Precision::FULL == this->setting.precision) {
139 using E = double;
140 arma_fortran(arma_dgbsv)(&N, &KL, &KU, &NRHS, (E*)this->memptr(), &LDAB, this->pivot.memptr(), (E*)B.memptr(), &LDB, &INFO);
141 X = std::move(B);
142 }
143 else {
144 this->s_memory = this->to_float();
145 arma_fortran(arma_sgbtrf)(&N, &N, &KL, &KU, this->s_memory.memptr(), &LDAB, this->pivot.memptr(), &INFO);
146 if(0 == INFO) INFO = this->solve_trs(X, std::forward<Mat<T>>(B));
147 }
148
149 if(0 != INFO)
150 suanpan_error("Error code {} received, the matrix is probably singular.\n", INFO);
151
152 return INFO;
153}
154
155template<sp_d T> int BandMat<T>::solve_trs(Mat<T>& X, Mat<T>&& B) {
156 auto INFO = 0;
157
158 const auto N = static_cast<int>(this->n_rows);
159 const auto KL = static_cast<int>(l_band);
160 const auto KU = static_cast<int>(u_band);
161 const auto NRHS = static_cast<int>(B.n_cols);
162 const auto LDAB = static_cast<int>(m_rows);
163 const auto LDB = static_cast<int>(B.n_rows);
164
165 if constexpr(std::is_same_v<T, float>) {
166 using E = float;
167 arma_fortran(arma_sgbtrs)(&TRAN, &N, &KL, &KU, &NRHS, (E*)this->memptr(), &LDAB, this->pivot.memptr(), (E*)B.memptr(), &LDB, &INFO);
168 X = std::move(B);
169 }
170 else if(Precision::FULL == this->setting.precision) {
171 using E = double;
172 arma_fortran(arma_dgbtrs)(&TRAN, &N, &KL, &KU, &NRHS, (E*)this->memptr(), &LDAB, this->pivot.memptr(), (E*)B.memptr(), &LDB, &INFO);
173 X = std::move(B);
174 }
175 else
176 this->mixed_trs(X, std::forward<Mat<T>>(B), [&](fmat& residual) {
177 arma_fortran(arma_sgbtrs)(&TRAN, &N, &KL, &KU, &NRHS, this->s_memory.memptr(), &LDAB, this->pivot.memptr(), residual.memptr(), &LDB, &INFO);
178 return INFO;
179 });
180
181 if(0 != INFO)
182 suanpan_error("Error code {} received, the matrix is probably singular.\n", INFO);
183
184 return INFO;
185}
186
187#endif
188
A BandMat class that holds matrices.
Definition: BandMat.hpp:35
T & unsafe_at(const uword in_row, const uword in_col) override
Access element without bound check.
Definition: BandMat.hpp:78
BandMat(const uword in_size, const uword in_l, const uword in_u)
Definition: BandMat.hpp:55
T & at(const uword in_row, const uword in_col) override
Access element with bound check.
Definition: BandMat.hpp:83
const uword u_band
Definition: BandMat.hpp:48
unique_ptr< MetaMat< T > > make_copy() override
Definition: BandMat.hpp:65
T operator()(const uword in_row, const uword in_col) const override
Access element (read-only), returns zero if out-of-bound.
Definition: BandMat.hpp:73
const uword l_band
Definition: BandMat.hpp:47
const uword m_rows
Definition: BandMat.hpp:45
void nullify(const uword K) override
Definition: BandMat.hpp:67
A DenseMat class that holds matrices.
Definition: DenseMat.hpp:39
std::unique_ptr< T[]> memory
Definition: DenseMat.hpp:48
const uword n_cols
Definition: MetaMat.hpp:119
const uword n_rows
Definition: MetaMat.hpp:118
bool factored
Definition: MetaMat.hpp:74
Mat< T > operator*(const Mat< T > &) const override
Definition: BandMat.hpp:93
int direct_solve(Mat< T > &, Mat< T > &&) override
Definition: BandMat.hpp:117
void for_each(const IT start, const IT end, F &&FN)
Definition: utility.h:28
#define suanpan_warning(...)
Definition: suanPan.h:308
void suanpan_assert(const std::function< void()> &F)
Definition: suanPan.h:296
#define suanpan_error(...)
Definition: suanPan.h:309