5#ifndef PROXSUITE_LINALG_SPARSE_LDLT_FACTORIZE_HPP
6#define PROXSUITE_LINALG_SPARSE_LDLT_FACTORIZE_HPP
9#include <Eigen/OrderingMethods>
20 return { nrows *
isize(
sizeof(I)),
isize(
alignof(I)) };
24template<
typename T,
typename I>
31 using namespace _detail;
35 at.nrows() == a.ncols(),
36 at.ncols() == a.nrows(),
39 auto pai = a.row_indices();
40 auto pax = a.values();
42 auto patp = at.col_ptrs_mut();
43 auto pati = at.row_indices_mut();
44 auto patx = at.values_mut();
46 auto _work = stack.make_new(proxsuite::linalg::veg::Tag<I>{}, at.ncols());
47 auto work = _work.ptr_mut();
50 if (a.is_compressed()) {
52 util::wrapping_inc(mut(work[util::zero_extend(pai[p])]));
55 for (
usize j = 0; j < a.ncols(); ++j) {
56 isize col_start = a.col_start(j);
57 isize col_end = a.col_end(j);
58 for (
isize p = col_start; p < col_end; ++p) {
59 util::wrapping_inc(mut(work[util::zero_extend(pai[p])]));
65 for (
usize j = 0; j <
usize(at.ncols()); ++j) {
66 patp[j + 1] = util::checked_non_negative_plus(patp[j], work[j]);
71 auto col_start = a.col_start(j);
72 auto col_end = a.col_end(j);
73 for (
usize p = col_start; p < col_end; ++p) {
74 auto i = util::zero_extend(pai[p]);
75 auto q = util::zero_extend(work[i]);
79 util::wrapping_inc(mut(work[i]));
90 return { nrows *
isize(
sizeof(I)),
isize(
alignof(I)) };
100 using namespace _detail;
104 at.nrows() == a.ncols(),
105 at.ncols() == a.nrows(),
106 at.nnz() == a.nnz());
108 auto pai = a.row_indices();
110 auto patp = at.col_ptrs_mut();
111 auto pati = at.row_indices_mut();
113 auto _work = stack.make_new(proxsuite::linalg::veg::Tag<I>{}, at.ncols());
114 auto work = _work.ptr_mut();
118 util::wrapping_inc(mut(work[util::zero_extend(pai[p])]));
122 for (
usize j = 0; j <
usize(at.ncols()); ++j) {
123 patp[j + 1] = util::checked_non_negative_plus(patp[j], work[j]);
127 for (
usize j = 0; j <
usize(a.ncols()); ++j) {
128 auto col_start = a.col_start(j);
129 auto col_end = a.col_end(j);
130 for (
usize p = col_start; p < col_end; ++p) {
131 auto i = util::zero_extend(pai[p]);
132 auto q = util::zero_extend(work[i]);
135 util::wrapping_inc(mut(work[i]));
147template<
typename T,
typename I>
151 using namespace _detail;
154 l.nrows() == l.ncols(),
155 x.nrows() == l.nrows()
161 auto pli = l.row_indices();
162 auto plx = l.values();
164 auto px = x.as_slice_mut().ptr_mut();
166 for (
usize j = 0; j < n; ++j) {
167 auto const xj = px[j];
168 auto col_start = l.col_start(j);
169 auto col_end = l.col_end(j);
172 for (
usize p = col_start + 1; p < col_end; ++p) {
173 auto i = util::zero_extend(pli[p]);
174 px[i] -= plx[p] * xj;
186template<
typename T,
typename I>
190 using namespace _detail;
193 l.nrows() == l.ncols(),
194 x.nrows() == l.nrows()
200 auto pli = l.row_indices();
201 auto plx = l.values();
203 auto px = x.as_slice_mut().ptr_mut();
212 auto col_start = l.col_start(j);
213 auto col_end = l.col_end(j);
220 usize pstart = col_start + 1;
221 usize pcount = col_end - pstart;
224 for (; p < pstart + pcount / 4 * 4; p += 4) {
225 auto i0 = util::zero_extend(pli[p + 0]);
226 auto i1 = util::zero_extend(pli[p + 1]);
227 auto i2 = util::zero_extend(pli[p + 2]);
228 auto i3 = util::zero_extend(pli[p + 3]);
229 acc0 += plx[p + 0] * px[i0];
230 acc1 += plx[p + 1] * px[i1];
231 acc2 += plx[p + 2] * px[i2];
232 acc3 += plx[p + 3] * px[i3];
234 for (; p < pstart + pcount; ++p) {
235 auto i0 = util::zero_extend(pli[p + 0]);
236 acc0 += plx[p + 0] * px[i0];
239 acc0 = (acc0 + acc1) + (acc2 + acc3);
255 return { n *
isize{
sizeof(I) },
alignof(I) };
274 using namespace _detail;
277 auto pai = a.row_indices();
280 stack.make_new_for_overwrite(proxsuite::linalg::veg::Tag<I>{},
isize(n));
281 auto pancestors = _work.ptr_mut();
284 for (
usize k = 0; k < n; ++k) {
286 pancestors[k] = I(-1);
289 auto col_start = a.col_start(k);
290 auto col_end = a.col_end(k);
293 for (
usize p = col_start; p < col_end; ++p) {
295 auto i = util::zero_extend(pai[p]);
303 auto next =
usize(-1);
305 if (node ==
usize(-1) || node >= k) {
310 next = util::sign_extend(pancestors[node]);
314 pancestors[node] = I(k);
319 if (next ==
usize(-1)) {
334 return { (k + 1) *
isize{
sizeof(bool) },
alignof(bool) };
346 bool* pmarked)
noexcept -> I*
352 auto pai = a.row_indices();
353 auto pparent = parent;
355 auto col_start = a.col_start(k_);
356 auto col_end = a.col_end(k_);
363 for (
usize p = col_start; p < col_end; ++p) {
364 auto i = util::zero_extend(pai[p]);
380 s[
isize(len)] = I(i);
381 util::wrapping_inc(mut(len));
385 i = util::sign_extend(pparent[i]);
392 usize(len) *
sizeof(I));
395 top = util::wrapping_plus(top, -len);
398 for (
usize q = top; q < n; ++q) {
399 pmarked[s[q]] =
false;
419 I* pnext_child)
noexcept ->
usize
421 using namespace _detail;
427 while (top !=
usize(-1)) {
428 auto current_node = util::zero_extend(pstack[top]);
429 auto current_child = util::sign_extend(pfirst_child[current_node]);
432 if (current_child ==
usize(-1)) {
433 post[start_index] = I(current_node);
437 util::wrapping_dec(mut(top));
440 util::wrapping_inc(mut(top));
441 pstack[top] = I(current_child);
444 pfirst_child[current_node] = pnext_child[current_child];
462 return { (3 * n) *
isize(
sizeof(I)),
alignof(I) };
477 using namespace _detail;
479 auto _work = stack.make_new_for_overwrite(proxsuite::linalg::veg::Tag<I>{},
481 I* pwork = _work.ptr_mut();
484 I* pfirst_child = pstack + n;
485 I* pnext_child = pfirst_child + n;
489 pfirst_child[j] = I(-1);
498 if (parent[
isize(j)] != I(-1)) {
500 pnext_child[j] = pfirst_child[util::zero_extend(parent[
isize(j)])];
502 pfirst_child[util::zero_extend(parent[
isize(j)])] = I(j);
506 usize start_index = 0;
507 for (
usize root = 0; root <
usize(n); ++root) {
508 if (parent[
isize(root)] == I(-1)) {
510 post, root, start_index, pstack, pfirst_child, pnext_child);
526 I* pancestor)
noexcept -> I
528 using namespace _detail;
533 if (i <= j || util::wrapping_plus(pfirst[j], I(1)) <=
534 util::wrapping_plus(pmax_first[i], I(1))) {
541 pmax_first[i] = pfirst[j];
544 usize j_prev = util::sign_extend(pprev_leaf[i]);
547 pprev_leaf[i] = I(j);
550 if (j_prev ==
usize(-1)) {
558 if (lca == util::zero_extend(pancestor[lca])) {
561 lca = util::zero_extend(pancestor[lca]);
571 usize next = util::zero_extend(pancestor[node]);
572 pancestor[node] = I(lca);
589 isize{
sizeof(I) } * (1 + 5 * n + nnz),
603 using namespace _detail;
605 auto _at_work = stack.make_new_for_overwrite(proxsuite::linalg::veg::Tag<I>{},
606 1 + 5 *
isize(n) + a.nnz());
607 auto pat_work = _at_work.ptr_mut();
609 pat_work[n] = I(a.nnz());
613 pat_work,
nullptr, pat_work + n + 1,
617 auto patp = at.col_ptrs();
618 auto pati = at.row_indices();
620 auto pwork = pat_work + n + 1 + a.nnz();
622 auto pdelta = counts;
625 auto pmax_first = pwork + n;
626 auto pprev_leaf = pwork + 2 * n;
627 auto pancestor = pwork + 3 * n;
629 auto pcounts = counts;
631 auto pparent = parent;
633 for (
usize i = 0; i < 3 * n; ++i) {
636 for (
usize i = 0; i < n; ++i) {
641 for (
usize k = 0; k < n; ++k) {
643 auto j = util::zero_extend(ppost[k]);
650 pdelta[j] = (pfirst[j] == I(-1)) ? I(1) : I(0);
657 if (j ==
usize(-1) || pfirst[j] != I(-1)) {
661 j = util::sign_extend(pparent[j]);
666 for (
usize k = 0; k < n; ++k) {
668 auto j = util::zero_extend(ppost[k]);
671 if (pparent[j] != I(-1)) {
674 util::wrapping_dec(mut(pdelta[util::zero_extend(pparent[j])]));
677 auto col_start = util::zero_extend(patp[j]);
678 auto col_end = util::zero_extend(patp[j + 1]);
681 for (
usize p = col_start; p < col_end; ++p) {
682 auto i = util::zero_extend(pati[p]);
693 util::wrapping_inc(mut(pdelta[j]));
697 util::wrapping_dec(mut(pdelta[util::zero_extend(lca)]));
702 if (pparent[j] != -1) {
704 pancestor[j] = pparent[j];
709 for (
usize j = 0; j < n; ++j) {
710 if (parent[
isize(j)] != I(-1)) {
711 pcounts[util::zero_extend(parent[
isize(j)])] = util::wrapping_plus(
712 pcounts[util::zero_extend(parent[
isize(j)])], pcounts[j]);
722 return { nnz *
isize{
sizeof(char) },
alignof(char) };
732 isize n = mat.nrows();
733 isize nnz = mat.nnz();
735 Eigen::PermutationMatrix<-1, -1, I> perm_eigen;
736 auto _ = stack.make_new(proxsuite::linalg::veg::Tag<char>{}, nnz);
738 Eigen::AMDOrdering<I>{}(
739 Eigen::Map<Eigen::SparseMatrix<char, Eigen::ColMajor, I>
const>{
748 .template selfadjointView<Eigen::Upper>(),
753 perm_eigen.indices().data(),
754 usize(n) *
sizeof(I));
763 perm_inv[util::zero_extend(perm[i])] = I(i);
773 return { n *
isize{
sizeof(I) },
alignof(I) };
780 return { n *
isize{
sizeof(I) },
alignof(I) };
791 for (
usize old_j = 0; old_j < n; ++old_j) {
792 usize new_j = util::zero_extend(pperm_inv[old_j]);
795 auto col_end = old_a.
col_end(old_j);
797 for (
usize p = col_start; p < col_end; ++p) {
800 if (old_i <= old_j) {
801 usize new_i = util::zero_extend(pperm_inv[old_i]);
802 util::wrapping_inc(mut(pcol_counts[new_i > new_j ? new_i : new_j]));
808 for (
usize i = 0; i < n; ++i) {
810 util::checked_non_negative_plus(pnew_ap[i], pcol_counts[i]);
811 pcol_counts[i] = pnew_ap[i];
825 auto _work = stack.make_new(proxsuite::linalg::veg::Tag<I>{},
isize(n));
826 I* pcol_counts = _work.ptr_mut();
829 auto pold_ai = old_a.row_indices();
831 auto pnew_ap = new_a.col_ptrs_mut();
832 auto pnew_ai = new_a.row_indices_mut();
834 auto pperm_inv = perm_inv;
838 auto pcurrent_row_index = pcol_counts;
840 for (
usize old_j = 0; old_j < n; ++old_j) {
841 usize new_j = util::zero_extend(pperm_inv[old_j]);
843 auto col_start = old_a.col_start(old_j);
844 auto col_end = old_a.col_end(old_j);
846 for (
usize p = col_start; p < col_end; ++p) {
847 usize old_i = util::zero_extend(pold_ai[p]);
849 if (old_i <= old_j) {
850 usize new_i = util::zero_extend(pperm_inv[old_i]);
852 usize new_max = new_i > new_j ? new_i : new_j;
853 usize new_min = new_i < new_j ? new_i : new_j;
855 auto row_idx = pcurrent_row_index[new_max];
856 pnew_ai[row_idx] = I(new_min);
857 pcurrent_row_index[new_max] = util::wrapping_plus(row_idx, I(1));
863template<
typename T,
typename I>
871 auto _work = stack.make_new(proxsuite::linalg::veg::Tag<I>{},
isize(n));
872 I* pcol_counts = _work.ptr_mut();
875 auto pold_ai = old_a.row_indices();
877 auto pnew_ap = new_a.col_ptrs_mut();
878 auto pnew_ai = new_a.row_indices_mut();
880 auto pperm_inv = perm_inv;
883 n, pperm_inv, old_a.symbolic(), pnew_ap, pcol_counts);
885 auto pcurrent_row_index = pcol_counts;
887 auto pold_ax = old_a.values();
888 auto pnew_ax = new_a.values_mut();
889 for (
usize old_j = 0; old_j < n; ++old_j) {
890 usize new_j = util::zero_extend(pperm_inv[old_j]);
892 auto col_start = old_a.col_start(old_j);
893 auto col_end = old_a.col_end(old_j);
895 for (
usize p = col_start; p < col_end; ++p) {
896 usize old_i = util::zero_extend(pold_ai[p]);
898 if (old_i <= old_j) {
899 usize new_i = util::zero_extend(pperm_inv[old_i]);
901 usize new_max = new_i > new_j ? new_i : new_j;
902 usize new_min = new_i < new_j ? new_i : new_j;
904 auto row_idx = pcurrent_row_index[new_max];
905 pnew_ai[row_idx] = I(new_min);
906 pnew_ax[row_idx] = pold_ax[p];
907 pcurrent_row_index[new_max] = util::wrapping_plus(row_idx, I(1));
939 constexpr isize sz{
sizeof(I) };
940 constexpr isize al{
alignof(I) };
942 StackReq perm_req{ 0, al };
952 perm_req = perm_req & StackReq{ (n + 1 + nnz) * sz, al };
958 StackReq parent_req = { n * sz, al };
959 StackReq post_req = { n * sz, al };
994 bool id_perm = perm_inv ==
nullptr;
995 bool user_perm = perm !=
nullptr;
1001 proxsuite::linalg::veg::Tag<I> tag{};
1010 auto amd_perm = stack.make_new_for_overwrite(tag,
isize(n));
1012 perm = amd_perm.ptr();
1022 auto _permuted_a_col_ptrs =
1024 .make_new_for_overwrite(tag, id_perm ? 0 : (a.ncols() + 1));
1025 auto _permuted_a_row_indices =
1027 .make_new_for_overwrite(tag, id_perm ? 0 : (a.nnz()));
1030 _permuted_a_col_ptrs.as_mut()[0] = 0;
1031 _permuted_a_col_ptrs.as_mut()[
isize(n)] = I(a.nnz());
1037 _permuted_a_col_ptrs.ptr_mut(),
1039 _permuted_a_row_indices.ptr_mut(),
1050 _permuted_a_col_ptrs.ptr(),
1052 _permuted_a_row_indices.ptr(),
1057 auto _post = stack.make_new_for_overwrite(tag,
isize(n));
1092 auto pcol_ptrs = col_ptrs;
1093 pcol_ptrs[0] = I(0);
1094 for (
usize i = 0; i < n; ++i) {
1096 util::checked_non_negative_plus(pcol_ptrs[i + 1], pcol_ptrs[i]);
1108template<
typename T,
typename I>
1111 proxsuite::linalg::veg::Tag<I> ,
1119 constexpr isize sz{
sizeof(I) };
1120 constexpr isize al{
alignof(I) };
1122 constexpr isize tsz{
sizeof(T) };
1123 constexpr isize tal{
alignof(T) };
1127 auto symb_perm_req = StackReq{ sz * (id_perm ? 0 : (n + 1 + a_nnz)), al };
1128 auto num_perm_req = StackReq{ tsz * (id_perm ? 0 : a_nnz), tal };
1130 & (StackReq{ tsz * n, tal }
1132 & (StackReq{ 2 * n * sz, al }
1133 & (StackReq{ n * tsz, tal }
1134 & StackReq{ n *
isize{
sizeof(bool) },
alignof(bool) }))));
1156template<
typename T,
typename I>
1169 using namespace _detail;
1170 isize n = a.nrows();
1172 bool id_perm = perm_inv ==
nullptr;
1174 proxsuite::linalg::veg::Tag<I> tag{};
1176 auto _permuted_a_values = stack.make_new_for_overwrite(
1177 proxsuite::linalg::veg::Tag<T>{}, id_perm ? 0 : a.nnz());
1179 auto _x = stack.make_new_for_overwrite(proxsuite::linalg::veg::Tag<T>{}, n);
1181 auto _permuted_a_col_ptrs =
1182 stack.make_new_for_overwrite(tag, id_perm ? 0 : (a.ncols() + 1));
1183 auto _permuted_a_row_indices =
1184 stack.make_new_for_overwrite(tag, id_perm ? 0 : a.nnz());
1187 _permuted_a_col_ptrs.as_mut()[0] = 0;
1188 _permuted_a_col_ptrs.as_mut()[n] = I(a.nnz());
1194 _permuted_a_col_ptrs.ptr_mut(),
1196 _permuted_a_row_indices.ptr_mut(),
1197 _permuted_a_values.ptr_mut(),
1208 _permuted_a_col_ptrs.ptr(),
1210 _permuted_a_row_indices.ptr(),
1211 _permuted_a_values.ptr(),
1214 auto _current_row_index = stack.make_new_for_overwrite(tag, n);
1215 auto _ereach_stack_storage = stack.make_new_for_overwrite(tag, n);
1217 I* pcurrent_row_index = _current_row_index.ptr_mut();
1218 T* px = _x.ptr_mut();
1223 usize(n) *
sizeof(I));
1230 I
const* plp = col_ptrs;
1232 auto _marked = stack.make_new(proxsuite::linalg::veg::Tag<bool>{}, n);
1233 for (
usize iter = 0; iter <
usize(n); ++iter) {
1234 usize ereach_count = 0;
1236 _ereach_stack_storage.ptr_mut(),
1242 auto pereach_stack = ereach_stack;
1245 T
const* pax = permuted_a.
values();
1247 I* pli = row_indices;
1251 auto col_start = permuted_a.
col_start(iter);
1252 auto col_end = permuted_a.
col_end(iter);
1257 for (
usize p = col_start; p < col_end; ++p) {
1258 auto i = util::zero_extend(pai[p]);
1262 T d = px[iter] + ((diag_to_add ==
nullptr || perm ==
nullptr)
1264 : diag_to_add[util::zero_extend(perm[iter])]);
1269 for (
usize q = 0; q < ereach_count; ++q) {
1270 usize j = util::zero_extend(pereach_stack[q]);
1271 auto col_start = util::zero_extend(plp[j]);
1272 auto row_idx = util::zero_extend(pcurrent_row_index[j]) + 1;
1275 T
const dj = plx[col_start];
1276 T
const lkj = xj / dj;
1282 for (
usize p = col_start + 1; p < row_idx; ++p) {
1283 auto i = util::zero_extend(pli[p]);
1284 px[i] -= plx[p] * xj;
1289 pli[row_idx] = I(iter);
1291 pcurrent_row_index[j] = I(row_idx);
1294 auto col_start = util::zero_extend(plp[iter]);
1295 pli[col_start] = I(iter);
#define VEG_ASSERT_ALL_OF(...)
#define HEDLEY_FALL_THROUGH
VEG_INLINE auto postorder_depth_first_search(I *post, usize root, usize start_index, I *pstack, I *pfirst_child, I *pnext_child) noexcept -> usize
void inv_perm(I *perm_inv, I const *perm, isize n) noexcept
void symmetric_permute(MatMut< T, I > new_a, MatRef< T, I > old_a, I const *perm_inv, DynStackMut stack) noexcept(VEG_CONCEPT(nothrow_copyable< T >))
void symmetric_permute_common(usize n, I const *pperm_inv, SymbolicMatRef< I > old_a, I *pnew_ap, I *pcol_counts)
VEG_INLINE auto least_common_ancestor(usize i, usize j, I const *pfirst, I *pmax_first, I *pprev_leaf, I *pancestor) noexcept -> I
auto symmetric_permute_req(proxsuite::linalg::veg::Tag< I >, isize n) noexcept -> proxsuite::linalg::veg::dynstack::StackReq
VEG_NODISCARD VEG_INLINE auto ereach(usize &count, I *s, SymbolicMatRef< I > a, I const *parent, isize k, bool *pmarked) noexcept -> I *
auto ereach_req(isize k) noexcept -> proxsuite::linalg::veg::dynstack::StackReq
void symmetric_permute_symbolic(SymbolicMatMut< I > new_a, SymbolicMatRef< I > old_a, I const *perm_inv, DynStackMut stack) noexcept
auto symmetric_permute_symbolic_req(proxsuite::linalg::veg::Tag< I >, isize n) noexcept -> proxsuite::linalg::veg::dynstack::StackReq
void factorize_symbolic_non_zeros(I *nnz_per_col, I *etree, I *perm_inv, I const *perm, SymbolicMatRef< I > a, DynStackMut stack) noexcept
void factorize_numeric(T *values, I *row_indices, proxsuite::linalg::veg::DoNotDeduce< T const * > diag_to_add, proxsuite::linalg::veg::DoNotDeduce< I const * > perm, I const *col_ptrs, I const *etree, I const *perm_inv, MatRef< T, I > a, DynStackMut stack) noexcept(false)
auto amd_req(proxsuite::linalg::veg::Tag< I >, isize, isize nnz) noexcept -> proxsuite::linalg::veg::dynstack::StackReq
auto postorder_req(proxsuite::linalg::veg::Tag< I >, isize n) noexcept -> proxsuite::linalg::veg::dynstack::StackReq
void postorder(I *post, I const *parent, isize n, DynStackMut stack) noexcept
auto transpose_req(proxsuite::linalg::veg::Tag< I >, isize nrows) noexcept -> proxsuite::linalg::veg::dynstack::StackReq
void transpose(MatMut< T, I > at, MatRef< T, I > a, DynStackMut stack) noexcept(VEG_CONCEPT(nothrow_copyable< T >))
VEG_INLINE void etree(I *parent, SymbolicMatRef< I > a, DynStackMut stack) noexcept
auto factorize_symbolic_req(proxsuite::linalg::veg::Tag< I > tag, isize n, isize nnz, Ordering o) noexcept -> proxsuite::linalg::veg::dynstack::StackReq
void factorize_symbolic_col_counts(I *col_ptrs, I *etree, I *perm_inv, I const *perm, SymbolicMatRef< I > a, DynStackMut stack) noexcept
void transpose_symbolic(SymbolicMatMut< I > at, SymbolicMatRef< I > a, DynStackMut stack) noexcept
auto transpose_symbolic_req(proxsuite::linalg::veg::Tag< I >, isize nrows) noexcept -> proxsuite::linalg::veg::dynstack::StackReq
auto factorize_numeric_req(proxsuite::linalg::veg::Tag< T >, proxsuite::linalg::veg::Tag< I >, isize n, isize a_nnz, Ordering o) noexcept -> proxsuite::linalg::veg::dynstack::StackReq
void dense_ltsolve(DenseVecMut< T > x, MatRef< T, I > l) noexcept(false)
void column_counts(I *counts, SymbolicMatRef< I > a, I const *parent, I const *post, DynStackMut stack) noexcept
auto etree_req(proxsuite::linalg::veg::Tag< I >, isize n) noexcept -> proxsuite::linalg::veg::dynstack::StackReq
void dense_lsolve(DenseVecMut< T > x, MatRef< T, I > l) noexcept(false)
auto column_counts_req(proxsuite::linalg::veg::Tag< I > tag, isize n, isize nnz) noexcept -> proxsuite::linalg::veg::dynstack::StackReq
meta::type_identity_t< T > DoNotDeduce
_detail::_meta::make_signed< usize >::Type isize
decltype(sizeof(0)) usize
auto symbolic() const noexcept -> SymbolicMatRef< I >
auto values() const noexcept -> T const *
auto col_start(usize j) const noexcept -> usize
auto col_end(usize j) const noexcept -> usize
auto row_indices() const noexcept -> I const *