suanPan
triplet_form.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 ******************************************************************************/
17
18#ifndef TRIPLET_FORM_HPP
19#define TRIPLET_FORM_HPP
20
21#include <Toolbox/utility.h>
22#include <numeric>
23
24template<sp_d data_t, sp_i index_t> class csc_form;
25template<sp_d data_t, sp_i index_t> class csr_form;
26
27enum class SparseBase : short unsigned {
28 ZERO,
29 ONE
30};
31
32template<sp_i index_t> class csr_comparator {
33 const index_t* const row_idx;
34 const index_t* const col_idx;
35
36public:
37 csr_comparator(const index_t* const in_row_idx, const index_t* const in_col_idx)
38 : row_idx(in_row_idx)
39 , col_idx(in_col_idx) {}
40
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];
44 }
45};
46
47template<sp_i index_t> class csc_comparator {
48 const index_t* const row_idx;
49 const index_t* const col_idx;
50
51public:
52 csc_comparator(const index_t* const in_row_idx, const index_t* const in_col_idx)
53 : row_idx(in_row_idx)
54 , col_idx(in_col_idx) {}
55
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];
59 }
60};
61
62template<sp_d data_t, sp_i index_t> class triplet_form final {
63 const data_t bin = data_t(0);
64
65 using index_ptr = std::unique_ptr<index_t[]>;
66 using data_ptr = std::unique_ptr<data_t[]>;
67
68 index_ptr row_idx = nullptr; // index storage
69 index_ptr col_idx = nullptr; // index storage
70 data_ptr val_idx = nullptr; // value storage
71
72 bool csc_sorted = false;
73 bool csr_sorted = false;
74
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); }
76
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 {
78 suanpan::for_each(n_elem, [&](const index_t I) {
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]);
82 });
83 }
84
88 void reserve(const index_t in_elem) {
89 if(in_elem <= n_alloc) return;
90
91 access::rw(n_alloc) = index_t(std::pow(2., std::ceil(std::log2(in_elem)) + 1));
92
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]);
96
97 copy_to(new_row_idx, new_col_idx, new_val_idx, 0, 0, 0, 1);
98
99 row_idx = std::move(new_row_idx);
100 col_idx = std::move(new_col_idx);
101 val_idx = std::move(new_val_idx);
102 }
103
104 void invalidate_sorting_flag() { csc_sorted = csr_sorted = false; }
105
106 void condense(bool = false);
107
108 void populate_diagonal() {
109 const auto t_elem = std::min(n_rows, n_cols);
110 reserve(n_elem + t_elem);
111 suanpan::for_each(t_elem, [&](const index_t I) {
112 row_idx[n_elem + I] = I;
113 col_idx[n_elem + I] = I;
114 val_idx[n_elem + I] = data_t(0);
115 });
116 access::rw(n_elem) += t_elem;
117 invalidate_sorting_flag();
118 }
119
120public:
121 using data_type = data_t;
122 using index_type = index_t;
123
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;
126 template<sp_d in_dt, sp_i in_it> friend class triplet_form;
127
128 const index_t n_rows = 0;
129 const index_t n_cols = 0;
130 const index_t n_elem = 0;
131 const index_t n_alloc = 0;
132
133 triplet_form() = default;
136 triplet_form& operator=(const triplet_form&);
137 triplet_form& operator=(triplet_form&&) noexcept;
138 ~triplet_form() = default;
139
140 triplet_form(const index_t in_rows, const index_t in_cols, const index_t in_elem = index_t(0))
141 : n_rows(in_rows)
142 , n_cols(in_cols) { init(in_elem); }
143
144 template<sp_d in_dt> explicit triplet_form(const SpMat<in_dt>&);
145 template<sp_d in_dt, sp_i in_it> explicit triplet_form(triplet_form<in_dt, in_it>&, SparseBase = SparseBase::ZERO, bool = false);
146
147 [[nodiscard]] const index_t* row_mem() const { return row_idx.get(); }
148
149 [[nodiscard]] const index_t* col_mem() const { return col_idx.get(); }
150
151 [[nodiscard]] const data_t* val_mem() const { return val_idx.get(); }
152
153 [[nodiscard]] index_t* row_mem() { return row_idx.get(); }
154
155 [[nodiscard]] index_t* col_mem() { return col_idx.get(); }
156
157 [[nodiscard]] data_t* val_mem() { return val_idx.get(); }
158
159 [[nodiscard]] index_t row(const index_t I) const { return row_idx[I]; }
160
161 [[nodiscard]] index_t col(const index_t I) const { return col_idx[I]; }
162
163 [[nodiscard]] data_t val(const index_t I) const { return val_idx[I]; }
164
165 [[nodiscard]] bool is_csr_sorted() const { return csr_sorted; }
166
167 [[nodiscard]] bool is_csc_sorted() const { return csc_sorted; }
168
169 [[nodiscard]] bool is_empty() const { return 0 == n_elem; }
170
171 [[nodiscard]] data_t max() const {
172 if(is_empty()) return data_t(0);
173 return *suanpan::max_element(val_idx.get(), val_idx.get() + n_elem);
174 }
175
176 void zeros() {
177 access::rw(n_elem) = 0;
178 invalidate_sorting_flag();
179 }
180
181 void init(const index_t in_elem) {
182 zeros();
183 reserve(in_elem);
184 }
185
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;
189 init(in_elem);
190 }
191
192 data_t operator()(const index_t row, const index_t col) const {
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.;
195 }
196
197 data_t& at(index_t, index_t);
198
199 void print() const;
200
201 void csr_sort();
202 void csc_sort();
203
205 csr_sort();
206 condense(false);
207 }
208
210 csc_sort();
211 condense(false);
212 }
213
215 populate_diagonal();
216 csr_sort();
217 condense(true);
218 }
219
221 populate_diagonal();
222 csc_sort();
223 condense(true);
224 }
225
226 void assemble(const Mat<data_t>&, const Col<uword>&);
227 template<sp_d in_dt, sp_i in_it> void assemble(const triplet_form<in_dt, in_it>&, index_t, index_t, data_t);
228
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"); });
231
232 reserve(n_elem + index_t(scalar.size()) * index_t(in_mat.n_elem));
233
234 for(size_t I = 0; I < scalar.size(); ++I) assemble(in_mat, row_shift[I], col_shift[I], scalar[I]);
235 }
236
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);
239
240 for(index_t I = 0; I < n_elem; ++I) out_mat(row_idx[I]) += val_idx[I] * in_mat(col_idx[I]);
241
242 return out_mat;
243 }
244
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);
247
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]);
249
250 return out_mat;
251 }
252
253 template<sp_d T2> triplet_form operator*(T2) const;
254 template<sp_d T2> triplet_form operator/(T2) const;
255 template<sp_d T2> triplet_form& operator*=(T2);
256 template<sp_d T2> triplet_form& operator/=(T2);
257
258 triplet_form operator+(const triplet_form& in_mat) const {
259 triplet_form copy = *this;
260 return copy += in_mat;
261 }
262
263 triplet_form operator-(const triplet_form& in_mat) const {
264 triplet_form copy = *this;
265 return copy -= in_mat;
266 }
267
270
271 [[nodiscard]] Col<data_t> diag() const;
272 [[nodiscard]] triplet_form diagonal() const;
273 [[nodiscard]] triplet_form strictly_upper() const;
274 [[nodiscard]] triplet_form strictly_lower() const;
275 [[nodiscard]] triplet_form upper() const;
276 [[nodiscard]] triplet_form lower() const;
277};
278
279template<sp_d data_t, sp_i index_t> void triplet_form<data_t, index_t>::condense(const bool full) {
280 if(n_elem < 2) return;
281
282 auto last_row = row_idx[0], last_col = col_idx[0];
283
284 sp_i auto current_pos = index_t(0);
285 sp_d auto last_sum = data_t(0);
286
287 auto populate = [&] {
288 if(suanpan::approx_equal(last_sum, data_t(0)) && (!full || last_row != last_col)) return;
289 row_idx[current_pos] = last_row;
290 col_idx[current_pos] = last_col;
291 val_idx[current_pos] = last_sum;
292 ++current_pos;
293 last_sum = data_t(0);
294 };
295
296 for(index_t I = 0; I < n_elem; ++I) {
297 if(last_row != row_idx[I] || last_col != col_idx[I]) {
298 populate();
299 last_row = row_idx[I];
300 last_col = col_idx[I];
301 }
302 last_sum += val_idx[I];
303 }
304
305 populate();
306
307 access::rw(n_elem) = current_pos;
308}
309
310template<sp_d data_t, sp_i index_t> triplet_form<data_t, index_t>::triplet_form(const triplet_form& in_mat)
311 : csc_sorted{in_mat.csc_sorted}
312 , csr_sorted{in_mat.csr_sorted}
313 , n_rows{in_mat.n_rows}
314 , n_cols{in_mat.n_cols} {
315 init(in_mat.n_alloc);
316 in_mat.copy_to(row_idx, col_idx, val_idx, 0, 0, 0, 1);
317 access::rw(n_elem) = in_mat.n_elem;
318}
319
320template<sp_d data_t, sp_i index_t> triplet_form<data_t, index_t>::triplet_form(triplet_form&& in_mat) noexcept
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}
326 , n_rows{in_mat.n_rows}
327 , n_cols{in_mat.n_cols}
328 , n_elem{in_mat.n_elem}
329 , n_alloc{in_mat.n_alloc} {}
330
332 if(this == &in_mat) return *this;
333 csc_sorted = in_mat.csc_sorted;
334 csr_sorted = in_mat.csr_sorted;
335 access::rw(n_rows) = in_mat.n_rows;
336 access::rw(n_cols) = in_mat.n_cols;
337 init(in_mat.n_alloc);
338 in_mat.copy_to(row_idx, col_idx, val_idx, 0, 0, 0, 1);
339 access::rw(n_elem) = in_mat.n_elem;
340 return *this;
341}
342
343template<sp_d data_t, sp_i index_t> triplet_form<data_t, index_t>& triplet_form<data_t, index_t>::operator=(triplet_form&& in_mat) noexcept {
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);
354 return *this;
355}
356
357template<sp_d data_t, sp_i index_t> template<sp_d in_dt> triplet_form<data_t, index_t>::triplet_form(const SpMat<in_dt>& in_mat)
358 : n_rows(index_t(in_mat.n_rows))
359 , n_cols(index_t(in_mat.n_cols)) {
360 init(index_t(in_mat.n_nonzero));
361
362 for(auto I = in_mat.begin(); I != in_mat.end(); ++I) at(I.row(), I.col()) = *I;
363}
364
365template<sp_d data_t, sp_i index_t> template<sp_d in_dt, sp_i in_it> triplet_form<data_t, index_t>::triplet_form(triplet_form<in_dt, in_it>& in_mat, const SparseBase base, const bool full)
366 : csc_sorted(in_mat.csc_sorted)
367 , csr_sorted(in_mat.csr_sorted)
368 , n_rows(index_t(in_mat.n_rows))
369 , n_cols(index_t(in_mat.n_cols)) {
370 if(in_mat.is_empty()) return;
371
372 init(index_t(in_mat.n_alloc));
373
374 if(full) in_mat.full_csc_condense();
375 else in_mat.csc_condense();
376
377 const sp_i auto shift = index_t(base);
378
379 in_mat.copy_to(row_idx, col_idx, val_idx, 0, shift, shift, 1);
380
381 access::rw(n_elem) = index_t(in_mat.n_elem);
382}
383
384template<sp_d data_t, sp_i index_t> data_t& triplet_form<data_t, index_t>::at(const index_t row, const index_t col) {
385 suanpan_assert([&] { if(row >= n_rows || col >= n_cols) throw invalid_argument("inconsistent size"); });
386
387 invalidate_sorting_flag();
388 reserve(n_elem + 1);
389 row_idx[n_elem] = row;
390 col_idx[n_elem] = col;
391 return val_idx[access::rw(n_elem)++] = data_t(0);
392}
393
394template<sp_d data_t, sp_i index_t> void triplet_form<data_t, index_t>::print() const {
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)) {
397 suanpan_info("More than 1000 elements exist.\n");
398 return;
399 }
400 for(index_t I = 0; I < n_elem; ++I)
401 suanpan_info("({}, {}) ===> {:+.8E}\n", row_idx[I], col_idx[I], val_idx[I]);
402}
403
404template<sp_d data_t, sp_i index_t> void triplet_form<data_t, index_t>::csr_sort() {
405 if(csr_sorted) return;
406
407 std::vector<index_t> index(n_elem);
408 std::iota(index.begin(), index.end(), index_t(0));
409
410 suanpan_sort(index.begin(), index.end(), csr_comparator<index_t>(row_idx.get(), col_idx.get()));
411
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]);
415
416 suanpan::for_each(n_elem, [&](const index_t I) {
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]];
420 });
421
422 row_idx = std::move(new_row_idx);
423 col_idx = std::move(new_col_idx);
424 val_idx = std::move(new_val_idx);
425
426 csr_sorted = true;
427 csc_sorted = false;
428}
429
430template<sp_d data_t, sp_i index_t> void triplet_form<data_t, index_t>::csc_sort() {
431 if(csc_sorted) return;
432
433 std::vector<index_t> index(n_elem);
434 std::iota(index.begin(), index.end(), index_t(0));
435
436 suanpan_sort(index.begin(), index.end(), csc_comparator<index_t>(row_idx.get(), col_idx.get()));
437
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]);
441
442 suanpan::for_each(n_elem, [&](const index_t I) {
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]];
446 });
447
448 row_idx = std::move(new_row_idx);
449 col_idx = std::move(new_col_idx);
450 val_idx = std::move(new_val_idx);
451
452 csc_sorted = true;
453 csr_sorted = false;
454}
455
456template<sp_d data_t, sp_i index_t> void triplet_form<data_t, index_t>::assemble(const Mat<data_t>& in_mat, const Col<uword>& in_dof) {
457 if(in_mat.empty()) return;
458
459 invalidate_sorting_flag();
460
461 const auto t_elem = n_elem + index_t(in_mat.n_elem);
462
463 reserve(t_elem);
464
465 suanpan::for_each(in_mat.n_elem, [&](const uword I) {
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);
469 });
470
471 access::rw(n_elem) = t_elem;
472}
473
474template<sp_d data_t, sp_i index_t> template<sp_d in_dt, sp_i in_it> void triplet_form<data_t, index_t>::assemble(const triplet_form<in_dt, in_it>& in_mat, const index_t row_shift, const index_t col_shift, const data_t scalar) {
475 if(in_mat.is_empty()) return;
476
477 invalidate_sorting_flag();
478
479 const auto t_elem = n_elem + in_mat.n_elem;
480
481 reserve(t_elem);
482
483 in_mat.copy_to(row_idx, col_idx, val_idx, n_elem, row_shift, col_shift, scalar);
484
485 access::rw(n_elem) = t_elem;
486}
487
488template<sp_d data_t, sp_i index_t> template<sp_d T2> triplet_form<data_t, index_t> triplet_form<data_t, index_t>::operator*(const T2 scalar) const {
489 if(suanpan::approx_equal(T2(0), scalar)) return triplet_form(n_rows, n_cols);
490
491 auto copy = *this;
492
493 if(suanpan::approx_equal(T2(1), scalar)) return copy;
494
495 if(suanpan::approx_equal(T2(-1), scalar))
496 suanpan_for_each(copy.val_idx.get(), copy.val_idx.get() + copy.n_elem, [](data_t& I) { I = -I; });
497 else
498 suanpan_for_each(copy.val_idx.get(), copy.val_idx.get() + copy.n_elem, [=](data_t& I) { I *= data_t(scalar); });
499
500 return copy;
501}
502
503template<sp_d data_t, sp_i index_t> template<sp_d T2> triplet_form<data_t, index_t> triplet_form<data_t, index_t>::operator/(const T2 scalar) const {
504 auto copy = *this;
505
506 if(suanpan::approx_equal(T2(1), scalar)) return copy;
507
508 if(suanpan::approx_equal(T2(-1), scalar))
509 suanpan_for_each(copy.val_idx.get(), copy.val_idx.get() + copy.n_elem, [](data_t& I) { I = -I; });
510 else
511 suanpan_for_each(copy.val_idx.get(), copy.val_idx.get() + copy.n_elem, [=](data_t& I) { I /= data_t(scalar); });
512
513 return copy;
514}
515
516template<sp_d data_t, sp_i index_t> template<sp_d T2> triplet_form<data_t, index_t>& triplet_form<data_t, index_t>::operator*=(const T2 scalar) {
517 if(suanpan::approx_equal(T2(1), scalar)) return *this;
518
519 if(suanpan::approx_equal(T2(0), scalar)) zeros();
520 else if(suanpan::approx_equal(T2(-1), scalar))
521 suanpan_for_each(val_idx.get(), val_idx.get() + n_elem, [](data_t& I) { I = -I; });
522 else
523 suanpan_for_each(val_idx.get(), val_idx.get() + n_elem, [=](data_t& I) { I *= data_t(scalar); });
524
525 return *this;
526}
527
528template<sp_d data_t, sp_i index_t> template<sp_d T2> triplet_form<data_t, index_t>& triplet_form<data_t, index_t>::operator/=(const T2 scalar) {
529 if(suanpan::approx_equal(T2(1), scalar)) return *this;
530
531 if(suanpan::approx_equal(T2(-1), scalar))
532 suanpan_for_each(val_idx.get(), val_idx.get() + n_elem, [](data_t& I) { I = -I; });
533 else
534 suanpan_for_each(val_idx.get(), val_idx.get() + n_elem, [=](data_t& I) { I /= data_t(scalar); });
535
536 return *this;
537}
538
540 if(in_mat.is_empty()) return *this;
541
542 invalidate_sorting_flag();
543
544 const auto t_elem = n_elem + index_t(in_mat.n_elem);
545
546 reserve(t_elem);
547
548 in_mat.copy_to(row_idx, col_idx, val_idx, n_elem, 0, 0, 1);
549
550 access::rw(n_elem) = t_elem;
551
552 return *this;
553}
554
556 if(in_mat.is_empty()) return *this;
557
558 invalidate_sorting_flag();
559
560 const auto t_elem = n_elem + index_t(in_mat.n_elem);
561
562 reserve(t_elem);
563
564 in_mat.copy_to(row_idx, col_idx, val_idx, n_elem, 0, 0, -1);
565
566 access::rw(n_elem) = t_elem;
567
568 return *this;
569}
570
571template<sp_d data_t, sp_i index_t> Col<data_t> triplet_form<data_t, index_t>::diag() const {
572 Col<data_t> diag_vec(std::min(n_rows, n_cols), fill::zeros);
573 for(index_t I = 0; I < n_elem; ++I) if(row(I) == col(I)) diag_vec(row(I)) += val(I);
574 return diag_vec;
575}
576
578 auto out_mat = *this;
579
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); });
581
582 return out_mat;
583}
584
586 auto out_mat = *this;
587
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); });
589
590 return out_mat;
591}
592
594 auto out_mat = *this;
595
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); });
597
598 return out_mat;
599}
600
601template<sp_d data_t, sp_i index_t> triplet_form<data_t, index_t> triplet_form<data_t, index_t>::upper() const {
602 auto out_mat = *this;
603
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); });
605
606 return out_mat;
607}
608
609template<sp_d data_t, sp_i index_t> triplet_form<data_t, index_t> triplet_form<data_t, index_t>::lower() const {
610 auto out_mat = *this;
611
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); });
613
614 return out_mat;
615}
616
617template<sp_d data_t, sp_i index_t> triplet_form<data_t, index_t> operator+(const triplet_form<data_t, index_t>& mat_a, const triplet_form<data_t, index_t>& mat_b) {
618 auto out = mat_a;
619 out += mat_b;
620 return out;
621}
622
624 mat_a += mat_b;
625 return std::forward<triplet_form<data_t, index_t>>(mat_a);
626}
627
629 mat_b += mat_a;
630 return std::forward<triplet_form<data_t, index_t>>(mat_b);
631}
632
634 mat_a += mat_b;
635 return std::forward<triplet_form<data_t, index_t>>(mat_a);
636}
637
638#endif
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: csc_form.hpp:25
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: csr_form.hpp:25
Definition: triplet_form.hpp:62
void csr_sort()
Definition: triplet_form.hpp:404
data_t data_type
Definition: triplet_form.hpp:121
triplet_form strictly_lower() const
Definition: triplet_form.hpp:593
triplet_form operator+(const triplet_form &in_mat) const
Definition: triplet_form.hpp:258
void zeros()
Definition: triplet_form.hpp:176
data_t * val_mem()
Definition: triplet_form.hpp:157
void assemble(const triplet_form< in_dt, in_it > &, index_t, index_t, data_t)
Definition: triplet_form.hpp:474
triplet_form()=default
const index_t n_rows
Definition: triplet_form.hpp:128
triplet_form upper() const
Definition: triplet_form.hpp:601
const index_t n_alloc
Definition: triplet_form.hpp:131
data_t max() const
Definition: triplet_form.hpp:171
void csc_condense()
Definition: triplet_form.hpp:209
void csc_sort()
Definition: triplet_form.hpp:430
triplet_form strictly_upper() const
Definition: triplet_form.hpp:585
Mat< data_t > operator*(const Mat< data_t > &in_mat) const
Definition: triplet_form.hpp:245
triplet_form(triplet_form &&) noexcept
Definition: triplet_form.hpp:320
triplet_form diagonal() const
Definition: triplet_form.hpp:577
triplet_form lower() const
Definition: triplet_form.hpp:609
bool is_empty() const
Definition: triplet_form.hpp:169
data_t & at(index_t, index_t)
Definition: triplet_form.hpp:384
void full_csr_condense()
Definition: triplet_form.hpp:214
void csr_condense()
Definition: triplet_form.hpp:204
bool is_csc_sorted() const
Definition: triplet_form.hpp:167
triplet_form(const SpMat< in_dt > &)
Definition: triplet_form.hpp:357
triplet_form(const triplet_form &)
Definition: triplet_form.hpp:310
triplet_form & operator-=(const triplet_form &)
Definition: triplet_form.hpp:555
index_t * row_mem()
Definition: triplet_form.hpp:153
triplet_form & operator=(const triplet_form &)
Definition: triplet_form.hpp:331
triplet_form & operator+=(const triplet_form &)
Definition: triplet_form.hpp:539
index_t * col_mem()
Definition: triplet_form.hpp:155
data_t operator()(const index_t row, const index_t col) const
Definition: triplet_form.hpp:192
const index_t n_cols
Definition: triplet_form.hpp:129
triplet_form operator-(const triplet_form &in_mat) const
Definition: triplet_form.hpp:263
triplet_form & operator/=(T2)
const data_t * val_mem() const
Definition: triplet_form.hpp:151
triplet_form operator*(T2) const
Mat< data_t > operator*(const Col< data_t > &in_mat) const
Definition: triplet_form.hpp:237
void init(const index_t in_rows, const index_t in_cols, const index_t in_elem)
Definition: triplet_form.hpp:186
const index_t n_elem
Definition: triplet_form.hpp:130
triplet_form & operator*=(T2)
void init(const index_t in_elem)
Definition: triplet_form.hpp:181
const index_t * row_mem() const
Definition: triplet_form.hpp:147
triplet_form operator/(T2) const
void print() const
Definition: triplet_form.hpp:394
Col< data_t > diag() const
Definition: triplet_form.hpp:571
friend class triplet_form
Definition: triplet_form.hpp:126
bool is_csr_sorted() const
Definition: triplet_form.hpp:165
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)
Definition: triplet_form.hpp:229
index_t col(const index_t I) const
Definition: triplet_form.hpp:161
data_t val(const index_t I) const
Definition: triplet_form.hpp:163
void assemble(const Mat< data_t > &, const Col< uword > &)
Definition: triplet_form.hpp:456
const index_t * col_mem() const
Definition: triplet_form.hpp:149
index_t row(const index_t I) const
Definition: triplet_form.hpp:159
void full_csc_condense()
Definition: triplet_form.hpp:220
triplet_form(triplet_form< in_dt, in_it > &, SparseBase=SparseBase::ZERO, bool=false)
Definition: triplet_form.hpp:365
index_t index_type
Definition: triplet_form.hpp:122
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
SparseBase
Definition: triplet_form.hpp:27
triplet_form< data_t, index_t > operator+(const triplet_form< data_t, index_t > &mat_a, const triplet_form< data_t, index_t > &mat_b)
Definition: triplet_form.hpp:617