18#ifndef TRIPLET_FORM_HPP
19#define TRIPLET_FORM_HPP
24template<sp_d data_t, sp_i index_t>
class csc_form;
25template<sp_d data_t, sp_i index_t>
class csr_form;
33 const index_t*
const row_idx;
34 const index_t*
const col_idx;
37 csr_comparator(
const index_t*
const in_row_idx,
const index_t*
const in_col_idx)
39 , col_idx(in_col_idx) {}
41 bool operator()(
const index_t idx_a,
const index_t idx_b)
const {
42 if(row_idx[idx_a] == row_idx[idx_b])
return col_idx[idx_a] < col_idx[idx_b];
43 return row_idx[idx_a] < row_idx[idx_b];
48 const index_t*
const row_idx;
49 const index_t*
const col_idx;
52 csc_comparator(
const index_t*
const in_row_idx,
const index_t*
const in_col_idx)
54 , col_idx(in_col_idx) {}
56 bool operator()(
const index_t idx_a,
const index_t idx_b)
const {
57 if(col_idx[idx_a] == col_idx[idx_b])
return row_idx[idx_a] < row_idx[idx_b];
58 return col_idx[idx_a] < col_idx[idx_b];
63 const data_t bin = data_t(0);
65 using index_ptr = std::unique_ptr<index_t[]>;
66 using data_ptr = std::unique_ptr<data_t[]>;
68 index_ptr row_idx =
nullptr;
69 index_ptr col_idx =
nullptr;
70 data_ptr val_idx =
nullptr;
72 bool csc_sorted =
false;
73 bool csr_sorted =
false;
75 template<sp_d in_dt, sp_i in_it>
void copy_to(
const std::unique_ptr<in_it[]>& new_row_idx,
const std::unique_ptr<in_it[]>& new_col_idx,
const std::unique_ptr<in_dt[]>& new_val_idx,
const index_t
begin,
const index_t row_offset,
const index_t col_offset,
const data_t scalar)
const { copy_to(new_row_idx.get(), new_col_idx.get(), new_val_idx.get(),
begin, row_offset, col_offset, scalar); }
77 template<sp_d in_dt, sp_i in_it>
void copy_to(in_it*
const new_row_idx, in_it*
const new_col_idx, in_dt*
const new_val_idx,
const index_t
begin,
const index_t row_offset,
const index_t col_offset,
const data_t scalar)
const {
79 new_row_idx[I +
begin] = in_it(row_idx[I] + row_offset);
80 new_col_idx[I +
begin] = in_it(col_idx[I] + col_offset);
81 new_val_idx[I +
begin] = in_dt(scalar * val_idx[I]);
88 void reserve(
const index_t in_elem) {
91 access::rw(
n_alloc) = index_t(std::pow(2., std::ceil(std::log2(in_elem)) + 1));
93 index_ptr new_row_idx(
new index_t[
n_alloc]);
94 index_ptr new_col_idx(
new index_t[
n_alloc]);
95 data_ptr new_val_idx(
new data_t[
n_alloc]);
97 copy_to(new_row_idx, new_col_idx, new_val_idx, 0, 0, 0, 1);
99 row_idx = std::move(new_row_idx);
100 col_idx = std::move(new_col_idx);
101 val_idx = std::move(new_val_idx);
104 void invalidate_sorting_flag() { csc_sorted = csr_sorted =
false; }
106 void condense(
bool =
false);
108 void populate_diagonal() {
114 val_idx[
n_elem + I] = data_t(0);
116 access::rw(
n_elem) += t_elem;
117 invalidate_sorting_flag();
124 template<sp_d in_dt, sp_i in_it>
friend class csc_form;
125 template<sp_d in_dt, sp_i in_it>
friend class csr_form;
140 triplet_form(const index_t in_rows, const index_t in_cols, const index_t in_elem = index_t(0))
147 [[nodiscard]]
const index_t*
row_mem()
const {
return row_idx.get(); }
149 [[nodiscard]]
const index_t*
col_mem()
const {
return col_idx.get(); }
151 [[nodiscard]]
const data_t*
val_mem()
const {
return val_idx.get(); }
153 [[nodiscard]] index_t*
row_mem() {
return row_idx.get(); }
155 [[nodiscard]] index_t*
col_mem() {
return col_idx.get(); }
157 [[nodiscard]] data_t*
val_mem() {
return val_idx.get(); }
159 [[nodiscard]] index_t
row(
const index_t I)
const {
return row_idx[I]; }
161 [[nodiscard]] index_t
col(
const index_t I)
const {
return col_idx[I]; }
163 [[nodiscard]] data_t
val(
const index_t I)
const {
return val_idx[I]; }
171 [[nodiscard]] data_t
max()
const {
178 invalidate_sorting_flag();
181 void init(
const index_t in_elem) {
186 void init(
const index_t in_rows,
const index_t in_cols,
const index_t in_elem) {
187 access::rw(
n_rows) = in_rows;
188 access::rw(
n_cols) = in_cols;
193 for(index_t I = 0; I <
n_elem; ++I)
if(
row == row_idx[I] &&
col == col_idx[I])
return val_idx[I];
194 return access::rw(bin) = 0.;
197 data_t&
at(index_t, index_t);
226 void assemble(
const Mat<data_t>&,
const Col<uword>&);
229 template<sp_d in_dt, sp_i in_it>
void assemble(
const triplet_form<in_dt, in_it>& in_mat,
const std::vector<index_t>& row_shift,
const std::vector<index_t>& col_shift,
const std::vector<data_t>& scalar) {
230 suanpan_assert([&] {
if(scalar.size() != row_shift.size() || scalar.size() != col_shift.size())
throw invalid_argument(
"size mismatch detected"); });
232 reserve(
n_elem + index_t(scalar.size()) * index_t(in_mat.
n_elem));
234 for(
size_t I = 0; I < scalar.size(); ++I)
assemble(in_mat, row_shift[I], col_shift[I], scalar[I]);
237 Mat<data_t>
operator*(
const Col<data_t>& in_mat)
const {
238 Mat<data_t> out_mat(in_mat.n_rows, in_mat.n_cols, fill::zeros);
240 for(index_t I = 0; I <
n_elem; ++I) out_mat(row_idx[I]) += val_idx[I] * in_mat(col_idx[I]);
245 Mat<data_t>
operator*(
const Mat<data_t>& in_mat)
const {
246 Mat<data_t> out_mat(in_mat.n_rows, in_mat.n_cols, fill::zeros);
248 for(index_t I = 0; I <
n_elem; ++I) out_mat.row(row_idx[I]) += val_idx[I] * in_mat.row(col_idx[I]);
260 return copy += in_mat;
265 return copy -= in_mat;
271 [[nodiscard]] Col<data_t>
diag()
const;
282 auto last_row = row_idx[0], last_col = col_idx[0];
284 sp_i auto current_pos = index_t(0);
285 sp_d auto last_sum = data_t(0);
287 auto populate = [&] {
289 row_idx[current_pos] = last_row;
290 col_idx[current_pos] = last_col;
291 val_idx[current_pos] = last_sum;
293 last_sum = data_t(0);
296 for(index_t I = 0; I <
n_elem; ++I) {
297 if(last_row != row_idx[I] || last_col != col_idx[I]) {
299 last_row = row_idx[I];
300 last_col = col_idx[I];
302 last_sum += val_idx[I];
307 access::rw(
n_elem) = current_pos;
311 : csc_sorted{in_mat.csc_sorted}
312 , csr_sorted{in_mat.csr_sorted}
316 in_mat.copy_to(row_idx, col_idx, val_idx, 0, 0, 0, 1);
321 : row_idx{std::move(in_mat.row_idx)}
322 , col_idx{std::move(in_mat.col_idx)}
323 , val_idx{std::move(in_mat.val_idx)}
324 , csc_sorted{in_mat.csc_sorted}
325 , csr_sorted{in_mat.csr_sorted}
332 if(
this == &in_mat)
return *
this;
333 csc_sorted = in_mat.csc_sorted;
334 csr_sorted = in_mat.csr_sorted;
338 in_mat.copy_to(row_idx, col_idx, val_idx, 0, 0, 0, 1);
344 if(
this == &in_mat)
return *
this;
345 csc_sorted = in_mat.csc_sorted;
346 csr_sorted = in_mat.csr_sorted;
347 access::rw(
n_rows) = in_mat.n_rows;
348 access::rw(
n_cols) = in_mat.n_cols;
349 access::rw(
n_elem) = in_mat.n_elem;
350 access::rw(
n_alloc) = in_mat.n_alloc;
351 row_idx = std::move(in_mat.row_idx);
352 col_idx = std::move(in_mat.col_idx);
353 val_idx = std::move(in_mat.val_idx);
360 init(index_t(in_mat.n_nonzero));
362 for(
auto I = in_mat.begin(); I != in_mat.end(); ++I)
at(I.row(), I.col()) = *I;
366 : csc_sorted(in_mat.csc_sorted)
367 , csr_sorted(in_mat.csr_sorted)
377 const sp_i auto shift = index_t(base);
379 in_mat.copy_to(row_idx, col_idx, val_idx, 0, shift, shift, 1);
387 invalidate_sorting_flag();
391 return val_idx[access::rw(
n_elem)++] = data_t(0);
395 suanpan_info(
"A sparse matrix in triplet form with size of {} by {}, the density of {:.3f}%.\n",
n_rows,
n_cols,
static_cast<double>(
n_elem) /
static_cast<double>(
n_rows) /
static_cast<double>(
n_cols) * 1E2);
396 if(
n_elem > index_t(1000)) {
400 for(index_t I = 0; I <
n_elem; ++I)
401 suanpan_info(
"({}, {}) ===> {:+.8E}\n", row_idx[I], col_idx[I], val_idx[I]);
405 if(csr_sorted)
return;
407 std::vector<index_t> index(
n_elem);
408 std::iota(index.begin(), index.end(), index_t(0));
412 index_ptr new_row_idx(
new index_t[
n_alloc]);
413 index_ptr new_col_idx(
new index_t[
n_alloc]);
414 data_ptr new_val_idx(
new data_t[
n_alloc]);
417 new_row_idx[I] = row_idx[index[I]];
418 new_col_idx[I] = col_idx[index[I]];
419 new_val_idx[I] = val_idx[index[I]];
422 row_idx = std::move(new_row_idx);
423 col_idx = std::move(new_col_idx);
424 val_idx = std::move(new_val_idx);
431 if(csc_sorted)
return;
433 std::vector<index_t> index(
n_elem);
434 std::iota(index.begin(), index.end(), index_t(0));
438 index_ptr new_row_idx(
new index_t[
n_alloc]);
439 index_ptr new_col_idx(
new index_t[
n_alloc]);
440 data_ptr new_val_idx(
new data_t[
n_alloc]);
443 new_row_idx[I] = row_idx[index[I]];
444 new_col_idx[I] = col_idx[index[I]];
445 new_val_idx[I] = val_idx[index[I]];
448 row_idx = std::move(new_row_idx);
449 col_idx = std::move(new_col_idx);
450 val_idx = std::move(new_val_idx);
457 if(in_mat.empty())
return;
459 invalidate_sorting_flag();
461 const auto t_elem =
n_elem + index_t(in_mat.n_elem);
466 row_idx[n_elem + I] = index_t(in_dof(I % in_dof.n_elem));
467 col_idx[n_elem + I] = index_t(in_dof(I / in_dof.n_elem));
468 val_idx[n_elem + I] = in_mat(I);
471 access::rw(
n_elem) = t_elem;
477 invalidate_sorting_flag();
483 in_mat.copy_to(row_idx, col_idx, val_idx,
n_elem, row_shift, col_shift, scalar);
485 access::rw(
n_elem) = t_elem;
496 suanpan_for_each(copy.val_idx.get(), copy.val_idx.get() + copy.n_elem, [](data_t& I) { I = -I; });
498 suanpan_for_each(copy.val_idx.get(), copy.val_idx.get() + copy.n_elem, [=](data_t& I) { I *= data_t(scalar); });
509 suanpan_for_each(copy.val_idx.get(), copy.val_idx.get() + copy.n_elem, [](data_t& I) { I = -I; });
511 suanpan_for_each(copy.val_idx.get(), copy.val_idx.get() + copy.n_elem, [=](data_t& I) { I /= data_t(scalar); });
542 invalidate_sorting_flag();
548 in_mat.copy_to(row_idx, col_idx, val_idx,
n_elem, 0, 0, 1);
550 access::rw(
n_elem) = t_elem;
558 invalidate_sorting_flag();
564 in_mat.copy_to(row_idx, col_idx, val_idx,
n_elem, 0, 0, -1);
566 access::rw(
n_elem) = t_elem;
572 Col<data_t> diag_vec(std::min(
n_rows,
n_cols), fill::zeros);
578 auto out_mat = *
this;
580 suanpan::for_each(out_mat.n_elem, [&](
const index_t I) { out_mat.val_idx[I] *= out_mat.row(I) == out_mat.col(I); });
586 auto out_mat = *
this;
588 suanpan::for_each(out_mat.n_elem, [&](
const index_t I) { out_mat.val_idx[I] *= out_mat.row(I) < out_mat.col(I); });
594 auto out_mat = *
this;
596 suanpan::for_each(out_mat.n_elem, [&](
const index_t I) { out_mat.val_idx[I] *= out_mat.col(I) < out_mat.row(I); });
602 auto out_mat = *
this;
604 suanpan::for_each(out_mat.n_elem, [&](
const index_t I) { out_mat.val_idx[I] *= out_mat.row(I) <= out_mat.col(I); });
610 auto out_mat = *
this;
612 suanpan::for_each(out_mat.n_elem, [&](
const index_t I) { out_mat.val_idx[I] *= out_mat.col(I) <= out_mat.row(I); });
625 return std::forward<triplet_form<data_t, index_t>>(mat_a);
630 return std::forward<triplet_form<data_t, index_t>>(mat_b);
635 return std::forward<triplet_form<data_t, index_t>>(mat_a);
Storage< T >::iterator begin(Storage< T > &S)
Definition: Storage.hpp:200
Definition: triplet_form.hpp:47
bool operator()(const index_t idx_a, const index_t idx_b) const
Definition: triplet_form.hpp:56
csc_comparator(const index_t *const in_row_idx, const index_t *const in_col_idx)
Definition: triplet_form.hpp:52
Definition: triplet_form.hpp:32
bool operator()(const index_t idx_a, const index_t idx_b) const
Definition: triplet_form.hpp:41
csr_comparator(const index_t *const in_row_idx, const index_t *const in_col_idx)
Definition: triplet_form.hpp:37
Definition: suanPan.h:330
Definition: suanPan.h:331
std::enable_if_t<!std::numeric_limits< T >::is_integer, bool > approx_equal(T x, T y, int ulp=2)
Definition: utility.h:60
constexpr T max_element(T start, T end)
Definition: utility.h:39
void for_each(const IT start, const IT end, F &&FN)
Definition: utility.h:28
#define suanpan_info
Definition: suanPan.h:305
#define suanpan_for_each
Definition: suanPan.h:187
void suanpan_assert(const std::function< void()> &F)
Definition: suanPan.h:296
#define suanpan_sort
Definition: suanPan.h:186