suanPan
csc_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 CSC_FORM_HPP
19#define CSC_FORM_HPP
20
21#include "triplet_form.hpp"
22
23template<sp_d data_t, sp_i index_t> class csr_form;
24
25template<sp_d data_t, sp_i index_t> class csc_form final {
26 const data_t bin = data_t(0);
27
28 using index_ptr = std::unique_ptr<index_t[]>;
29 using data_ptr = std::unique_ptr<data_t[]>;
30
31 index_ptr row_idx = nullptr; // index storage
32 index_ptr col_ptr = nullptr; // index storage
33 data_ptr val_idx = nullptr; // value storage
34
35 template<sp_d in_dt, sp_i in_it> void copy_to(in_it* const new_row_idx, in_it* const new_col_ptr, in_dt* const new_val_idx) const {
36 suanpan::for_each(n_cols + 1, [&](const index_t I) { new_col_ptr[I] = in_it(col_ptr[I]); });
37 suanpan::for_each(n_elem, [&](const index_t I) {
38 new_row_idx[I] = in_it(row_idx[I]);
39 new_val_idx[I] = in_dt(val_idx[I]);
40 });
41 }
42
43 void init(const index_t in_elem) {
44 row_idx = std::move(index_ptr(new index_t[in_elem]));
45 col_ptr = std::move(index_ptr(new index_t[n_cols + 1]));
46 val_idx = std::move(data_ptr(new data_t[in_elem]));
47 }
48
49public:
50 const index_t n_rows = 0;
51 const index_t n_cols = 0;
52 const index_t n_elem = 0;
53
54 csc_form() = default;
55 csc_form(const csc_form&);
56 csc_form(csc_form&&) noexcept;
57 csc_form& operator=(const csc_form&);
58 csc_form& operator=(csc_form&&) noexcept;
59 ~csc_form() = default;
60
61 [[nodiscard]] const index_t* row_mem() const { return row_idx.get(); }
62
63 [[nodiscard]] const index_t* col_mem() const { return col_ptr.get(); }
64
65 [[nodiscard]] const data_t* val_mem() const { return val_idx.get(); }
66
67 [[nodiscard]] index_t* row_mem() { return row_idx.get(); }
68
69 [[nodiscard]] index_t* col_mem() { return col_ptr.get(); }
70
71 [[nodiscard]] data_t* val_mem() { return val_idx.get(); }
72
73 index_t row(const index_t I) const { return row_idx[I]; }
74
75 index_t col(const index_t I) const { return col_ptr[I]; }
76
77 data_t val(const index_t I) const { return val_idx[I]; }
78
79 [[nodiscard]] data_t max() const {
80 if(0 == n_elem) return data_t(0);
81 return *suanpan::max_element(val_idx.get(), val_idx.get() + n_elem);
82 }
83
84 void print() const;
85
86 template<sp_d T2> csc_form operator*(const T2 scalar) const {
87 auto copy = *this;
88 return copy *= scalar;
89 }
90
91 template<sp_d T2> csc_form operator/(const T2 scalar) const {
92 auto copy = *this;
93 return copy /= scalar;
94 }
95
96 template<sp_d T2> csc_form& operator*=(const T2 scalar) {
97 suanpan_for_each(val_idx.get(), val_idx.get() + n_elem, [=](data_t& I) { I *= data_t(scalar); });
98 return *this;
99 }
100
101 template<sp_d T2> csc_form& operator/=(const T2 scalar) {
102 suanpan_for_each(val_idx.get(), val_idx.get() + n_elem, [=](data_t& I) { I /= data_t(scalar); });
103 return *this;
104 }
105
106 template<sp_d in_dt, sp_i in_it> explicit csc_form(triplet_form<in_dt, in_it>&, SparseBase = SparseBase::ZERO, bool = false);
107
108 template<sp_d in_dt, sp_i in_it> csc_form& operator=(triplet_form<in_dt, in_it>&);
109
110 data_t operator()(const index_t in_row, const index_t in_col) const {
111 if(in_row < n_rows && in_col < n_cols) for(auto I = col_ptr[in_col]; I < col_ptr[in_col + 1]; ++I) if(in_row == row_idx[I]) return val_idx[I];
112 return access::rw(bin) = data_t(0);
113 }
114
115 Mat<data_t> operator*(const Col<data_t>& in_mat) const {
116 Mat<data_t> out_mat = arma::zeros<Mat<data_t>>(in_mat.n_rows, 1);
117 for(index_t I = 0; I < n_cols; ++I) for(auto J = col_ptr[I]; J < col_ptr[I + 1]; ++J) out_mat(row_idx[J]) += val_idx[J] * in_mat(I);
118 return out_mat;
119 }
120
121 Mat<data_t> operator*(const Mat<data_t>& in_mat) const {
122 Mat<data_t> out_mat = arma::zeros<Mat<data_t>>(in_mat.n_rows, in_mat.n_cols);
123 for(index_t I = 0; I < n_cols; ++I) for(auto J = col_ptr[I]; J < col_ptr[I + 1]; ++J) out_mat.row(row_idx[J]) += val_idx[J] * in_mat.row(I);
124 return out_mat;
125 }
126};
127
128template<sp_d data_t, sp_i index_t> csc_form<data_t, index_t>::csc_form(const csc_form& in_mat)
129 : n_rows{in_mat.n_rows}
130 , n_cols{in_mat.n_cols}
131 , n_elem{in_mat.n_elem} {
132 init(n_elem);
133 in_mat.copy_to(row_idx.get(), col_ptr.get(), val_idx.get());
134}
135
136template<sp_d data_t, sp_i index_t> csc_form<data_t, index_t>::csc_form(csc_form&& in_mat) noexcept
137 : row_idx{std::move(in_mat.row_idx)}
138 , col_ptr{std::move(in_mat.col_ptr)}
139 , val_idx{std::move(in_mat.val_idx)}
140 , n_rows{in_mat.n_rows}
141 , n_cols{in_mat.n_cols}
142 , n_elem{in_mat.n_elem} {}
143
144template<sp_d data_t, sp_i index_t> csc_form<data_t, index_t>& csc_form<data_t, index_t>::operator=(const csc_form& in_mat) {
145 if(this == &in_mat) return *this;
146 access::rw(n_rows) = in_mat.n_rows;
147 access::rw(n_cols) = in_mat.n_cols;
148 access::rw(n_elem) = in_mat.n_elem;
149 init(n_elem);
150 in_mat.copy_to(row_idx.get(), col_ptr.get(), val_idx.get());
151 return *this;
152}
153
154template<sp_d data_t, sp_i index_t> csc_form<data_t, index_t>& csc_form<data_t, index_t>::operator=(csc_form&& in_mat) noexcept {
155 if(this == &in_mat) return *this;
156 access::rw(n_rows) = in_mat.n_rows;
157 access::rw(n_cols) = in_mat.n_cols;
158 access::rw(n_elem) = in_mat.n_elem;
159 row_idx = std::move(in_mat.row_idx);
160 col_ptr = std::move(in_mat.col_ptr);
161 val_idx = std::move(in_mat.val_idx);
162 return *this;
163}
164
165template<sp_d data_t, sp_i index_t> void csc_form<data_t, index_t>::print() const {
166 suanpan_info("A sparse matrix in triplet form with size of {} by {}, the sparsity of {:.3f}%.\n", n_rows, n_cols, 1E2 - static_cast<double>(n_elem) / static_cast<double>(n_rows) / static_cast<double>(n_cols) * 1E2);
167 if(n_elem > index_t(1000)) {
168 suanpan_info("More than 1000 elements exist.\n");
169 return;
170 }
171
172 index_t c_idx = 1;
173 for(index_t I = 0; I < n_elem; ++I) {
174 if(I >= col_ptr[c_idx]) ++c_idx;
175 suanpan_info("({}, {}) ===> {:+.8E}\n", row_idx[I], c_idx - 1, val_idx[I]);
176 }
177}
178
179template<sp_d data_t, sp_i index_t> template<sp_d in_dt, sp_i in_it> csc_form<data_t, index_t>::csc_form(triplet_form<in_dt, in_it>& in_mat, const SparseBase base, const bool full)
180 : n_rows(index_t(in_mat.n_rows))
181 , n_cols(index_t(in_mat.n_cols)) {
182 if(full) in_mat.full_csc_condense();
183 else in_mat.csc_condense();
184
185 init(access::rw(n_elem) = index_t(in_mat.n_elem));
186
187 const sp_i auto shift = index_t(base);
188
189 suanpan::for_each(in_mat.n_elem, [&](const in_it I) {
190 row_idx[I] = index_t(in_mat.row_idx[I]) + shift;
191 val_idx[I] = data_t(in_mat.val_idx[I]);
192 });
193
194 in_it current_pos = 0, current_col = 0;
195
196 while(current_pos < in_mat.n_elem)
197 if(in_mat.col_idx[current_pos] < current_col) ++current_pos;
198 else col_ptr[current_col++] = index_t(current_pos) + shift;
199
200 col_ptr[0] = shift;
201 col_ptr[n_cols] = n_elem + shift;
202}
203
204template<sp_d data_t, sp_i index_t> template<sp_d in_dt, sp_i in_it> csc_form<data_t, index_t>& csc_form<data_t, index_t>::operator=(triplet_form<in_dt, in_it>& in_mat) {
205 in_mat.csc_condense();
206
207 access::rw(n_rows) = index_t(in_mat.n_rows);
208 access::rw(n_cols) = index_t(in_mat.n_cols);
209
210 init(access::rw(n_elem) = index_t(in_mat.n_elem));
211
212 suanpan::for_each(in_mat.n_elem, [&](const in_it I) {
213 row_idx[I] = index_t(in_mat.row_idx[I]);
214 val_idx[I] = data_t(in_mat.val_idx[I]);
215 });
216
217 in_it current_pos = 0, current_col = 0;
218
219 while(current_pos < in_mat.n_elem)
220 if(in_mat.col_idx[current_pos] < current_col) ++current_pos;
221 else col_ptr[current_col++] = index_t(current_pos);
222
223 col_ptr[0] = index_t(0);
224 col_ptr[n_cols] = n_elem;
225
226 return *this;
227}
228
229#endif
Definition: csc_form.hpp:25
void print() const
Definition: csc_form.hpp:165
csc_form & operator/=(const T2 scalar)
Definition: csc_form.hpp:101
const index_t n_rows
Definition: csc_form.hpp:50
csc_form & operator=(triplet_form< in_dt, in_it > &)
index_t * row_mem()
Definition: csc_form.hpp:67
csc_form()=default
csc_form & operator=(const csc_form &)
Definition: csc_form.hpp:144
Mat< data_t > operator*(const Mat< data_t > &in_mat) const
Definition: csc_form.hpp:121
data_t val(const index_t I) const
Definition: csc_form.hpp:77
index_t * col_mem()
Definition: csc_form.hpp:69
const index_t n_cols
Definition: csc_form.hpp:51
csc_form operator*(const T2 scalar) const
Definition: csc_form.hpp:86
data_t * val_mem()
Definition: csc_form.hpp:71
csc_form & operator*=(const T2 scalar)
Definition: csc_form.hpp:96
data_t max() const
Definition: csc_form.hpp:79
csc_form operator/(const T2 scalar) const
Definition: csc_form.hpp:91
Mat< data_t > operator*(const Col< data_t > &in_mat) const
Definition: csc_form.hpp:115
const index_t * col_mem() const
Definition: csc_form.hpp:63
data_t operator()(const index_t in_row, const index_t in_col) const
Definition: csc_form.hpp:110
const index_t * row_mem() const
Definition: csc_form.hpp:61
const data_t * val_mem() const
Definition: csc_form.hpp:65
index_t col(const index_t I) const
Definition: csc_form.hpp:75
const index_t n_elem
Definition: csc_form.hpp:52
index_t row(const index_t I) const
Definition: csc_form.hpp:73
Definition: csr_form.hpp:25
Definition: triplet_form.hpp:62
const index_t n_rows
Definition: triplet_form.hpp:128
void csc_condense()
Definition: triplet_form.hpp:209
const index_t n_cols
Definition: triplet_form.hpp:129
const index_t n_elem
Definition: triplet_form.hpp:130
void full_csc_condense()
Definition: triplet_form.hpp:220
Definition: suanPan.h:331
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
SparseBase
Definition: triplet_form.hpp:27