TDLS 0.4.0
Tiny Device-callable Linear Solvers
Loading...
Searching...
No Matches
solver_static.hpp
Go to the documentation of this file.
1#ifndef TDLS_SOLVERS_ARROWHEAD_LUPP_SOLVER_STATIC_HPP
2#define TDLS_SOLVERS_ARROWHEAD_LUPP_SOLVER_STATIC_HPP
3
4
5
106
107
108
109#include <cmath>
110#include <type_traits>
111#include <utility>
112
113#include <tdls/core/macros.hpp>
114#include <tdls/core/math.hpp>
115#include <tdls/core/target.hpp>
116#include <tdls/solvers/arrowhead_lupp/config.hpp>
117#include <tdls/solvers/arrowhead_lupp/storage_map.hpp>
118#include <tdls/solvers/tiled_lupp/solver_static.hpp>
119
120
121
122// The unrolling requested by TDLS_UNROLL_FORCE is a performance hint: a
123// failed or misplaced one does not affect correctness. clang reports a
124// forced unrolling that the optimizer could not perform through
125// -Wpass-failed (clang 12 does, on the multiplier pass). nvcc 12.0 places
126// the _Pragma of a loop nested in the body of another TDLS_ARROW_LUPP_FOR,
127// a macro argument, before the enclosing statement, and reports #20058
128// ("the #pragma unroll must precede a loop"); nvcc 12.9 and 13.0 place it
129// on the loop. Both suppressions are scoped to this header and to those
130// warnings only.
131#if defined(TDLS_DETAIL_COMPILER_CLANG)
132#pragma clang diagnostic push
133#pragma clang diagnostic ignored "-Wpass-failed"
134#endif
135#if defined(TDLS_DETAIL_COMPILER_NVCC)
136#pragma nv_diagnostic push
137#pragma nv_diag_suppress 20058
138#endif
139
140namespace tdls {
141
142
143
144/* Addressing macros. They expand inside member functions where the
145 parameter names (A, A_stride, y, rhs_stride, work, work_stride), the
146 residency template booleans and the class constants are in scope.
147 #undef'd at the end of this header. */
148
152#define TDLS_ARROW_LUPP_A(idx) A[internal_matrix ? unsigned(idx) : TDLS_STRIDED(idx, A_stride)]
156#define TDLS_ARROW_LUPP_Y(i) y[internal_rhs ? unsigned(i) : TDLS_STRIDED(i, rhs_stride)]
163#define TDLS_ARROW_LUPP_YW(m, i) \
164 y[internal_rhs ? unsigned((ycols == 1 ? 0 : (m)) * dimension + (i)) \
165 : TDLS_STRIDED(i, rhs_stride) + \
166 (ycols == 1 ? std::size_t(0) : TDLS_STRIDED_WIDE(m, xcol_stride))]
170#define TDLS_ARROW_LUPP_W(e) work[internal_work ? unsigned(e) : TDLS_STRIDED(e, work_stride)]
174#define TDLS_ARROW_LUPP_F(r, c) TDLS_ARROW_LUPP_W((r) * fcols + (c))
181#if defined(TDLS_DETAIL_SIDE_DISPATCH)
182#define TDLS_ARROW_LUPP_SIDED(...) \
183 if constexpr (sided) { \
184 if (!std::is_constant_evaluated()) { \
185 if (detail::on_device_side()) return SolverGpu::__VA_ARGS__; \
186 return SolverCpu::__VA_ARGS__; \
187 } \
188 }
189#else
190#define TDLS_ARROW_LUPP_SIDED(...)
191#endif
197#define TDLS_ARROW_LUPP_FOR(i, first, last, ...) \
198 if constexpr (Config.unroll_loops) { \
199 TDLS_UNROLL_FORCE \
200 for (int i = (first); i < (last); ++i) { \
201 __VA_ARGS__ \
202 } \
203 } else { \
204 for (int i = (first); i < (last); ++i) { \
205 __VA_ARGS__ \
206 } \
207 }
208
209
210
211namespace detail {
212
228template<typename T, int S, int B, bool col, bool diag>
229struct ArrowheadStreamedBlock {
230 T Dl[S * S];
231 T W[(B + 1) * S];
232 T Abk[diag ? S : B * S];
236 TDLS_HOST_DEVICE TDLS_FORCE_INLINE_ALL constexpr T& d(const int i, const int j) noexcept {
237 return Dl[col ? j * S + i : i * S + j];
238 }
242 TDLS_HOST_DEVICE TDLS_FORCE_INLINE_ALL constexpr T& akb(const int i, const int c) noexcept {
243 return W[c * S + i];
244 }
248 TDLS_HOST_DEVICE TDLS_FORCE_INLINE_ALL constexpr T& y(const int i) noexcept {
249 return W[B * S + i];
250 }
253 TDLS_HOST_DEVICE TDLS_FORCE_INLINE_ALL constexpr T& coupling(const int i) noexcept {
254 return Abk[i];
255 }
259 TDLS_HOST_DEVICE TDLS_FORCE_INLINE_ALL constexpr T& abk(const int r, const int i) noexcept {
260 return Abk[r * S + i];
261 }
262};
263
264
277template<typename T, int S, int B, bool col, bool diag>
278struct ArrowheadRecomputedBlock : ArrowheadStreamedBlock<T, S, B, col, diag> {
282 bool recovery = false;
283};
284
286enum class ArrowheadStreamedOutcome {
287 solved,
288 tangent,
290 singular
291};
292
295enum class ArrowheadRecomputedOutcome {
296 solved,
297 tangent,
299 singular,
300 deferred
303};
304
305} // namespace detail
306
307
308
314template<typename T, ArrowheadShape Shape, ArrowheadLUppConfig<T> Config>
315struct ArrowheadLUppSolverStatic {
317 using storage_map = ArrowheadStorageMap<Config.storage, Config.coupling, Config.layout>;
319 static constexpr ArrowheadShape shape = Shape;
321 static constexpr int dimension = Shape.n();
323 static constexpr int storage_size = storage_map::size(Shape);
324
325 private:
326 static constexpr int B = Shape.border;
327 static constexpr int S = Shape.block;
328 static constexpr int NB = Shape.blocks;
329 static constexpr bool col = Config.layout == MatrixLayout::ColMajor;
330 static constexpr bool diag = Config.coupling == ArrowheadCoupling::Diagonal;
333 static constexpr bool reg = detail::arrowhead_register_blocks(Config, S);
336 static constexpr int rhs_pass = detail::arrowhead_rhs_pass_width(Config, reg);
337#if defined(TDLS_DETAIL_SIDE_DISPATCH)
340 static constexpr ArrowheadLUppConfig<T> config_gpu =
341 detail::arrowhead_resolved_config<T, true>(Config, S);
343 static constexpr ArrowheadLUppConfig<T> config_cpu =
344 detail::arrowhead_resolved_config<T, false>(Config, S);
347 static constexpr bool sided =
348 config_gpu.block_row_interchange != config_cpu.block_row_interchange ||
349 config_gpu.rhs_pass_width != config_cpu.rhs_pass_width;
351 using SolverGpu = ArrowheadLUppSolverStatic<T, Shape, config_gpu>;
353 using SolverCpu = ArrowheadLUppSolverStatic<T, Shape, config_cpu>;
354#endif
357 template<typename U, ArrowheadShape OtherShape, ArrowheadLUppConfig<U> OtherConfig>
358 friend struct ArrowheadLUppSolverStatic;
363 static constexpr int R = diag ? S : B;
364 static constexpr int R0 = diag ? Shape.coupling_row : 0;
365
371 [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr ArrowheadShape
372 shape_value() noexcept {
373 return {B, S, NB, R0};
374 }
377 static constexpr int front_rows = S + B;
380 template<int nc>
381 static constexpr int front_cols_for = S + B + R + nc;
384 static constexpr int extra_cols = R > S ? R - S : 0;
385
386 public:
390 template<int nrhs>
391 static constexpr int multirhs_workspace_size =
392 front_rows * front_cols_for<nrhs> + NB * S * extra_cols;
395 static constexpr int workspace_size = multirhs_workspace_size<1>;
397 static constexpr int pass_width = rhs_pass == 0 ? B + 1 : rhs_pass;
398
399 private:
402 template<int ny>
403 static constexpr int pass_width_for = rhs_pass == 0 || rhs_pass > B + ny ? B + ny : rhs_pass;
408 template<int ny>
409 static constexpr bool forward_check_for = !reg && col && pass_width_for<ny> == B + ny;
414 template<int ny>
415 static constexpr bool stores_recovery_for = ny > 0;
423 static constexpr bool multiplier_update =
424 !col && !diag && Config.condensation == ArrowheadCondensation::MinFlops;
427 static constexpr bool in_place = Config.schur_in_place;
428
429 static_assert(Shape.border >= 1 && Shape.block >= 1 && Shape.blocks >= 1,
430 "ArrowheadLUppSolverStatic: border, block and blocks must be >= 1");
431 static_assert(diag || Shape.coupling_row == 0,
432 "ArrowheadLUppSolverStatic: coupling_row applies to the diagonal coupling only");
433 static_assert(!diag || (Shape.coupling_row >= 0 && Shape.coupling_row + S <= B),
434 "ArrowheadLUppSolverStatic: the diagonal coupling rows exceed the border");
435 static_assert(Config.storage != ArrowheadStorage::Bsr ||
436 (B % S == 0 && (!diag || Shape.coupling_row % S == 0)),
437 "ArrowheadLUppSolverStatic: the Bsr storage needs a border (and a coupling_row) "
438 "multiple of the block");
439 static_assert(storage_map::wide_size(Shape) < (1LL << 31),
440 "ArrowheadLUppSolverStatic: the storage exceeds 2^31 entries");
441 static_assert(Config.condensation_pivot_threshold.is_finite() &&
442 T(Config.condensation_pivot_threshold) >= T(0) &&
443 T(Config.condensation_pivot_threshold) <= T(1),
444 "ArrowheadLUppSolverStatic: condensation_pivot_threshold must be in [0, 1]");
445
449 static constexpr T condensation_pivot_threshold = Config.condensation_pivot_threshold;
450 static_assert(Config.rhs_pass_width >= -1 && Config.rhs_pass_width <= Shape.border + 1,
451 "ArrowheadLUppSolverStatic: rhs_pass_width must be in [-1, border + 1]");
452
453 using BlockSolver = TiledLUppSolverStatic<T, S, detail::arrowhead_block_config(Config)>;
454 using BorderSolver = TiledLUppSolverStatic<T, B, detail::arrowhead_border_config(Config)>;
457 using PhysicalBlockSolver =
458 TiledLUppSolverStatic<T, S, detail::arrowhead_physical_block_config(Config)>;
459
460 // The recomputed solve runs the second pass of the solver of its store
461 // of refused blocks, another instantiation.
462 template<typename U, ArrowheadShape OtherShape, ArrowheadLUppConfig<U> OtherConfig>
463 friend struct ArrowheadLUppSolverStatic;
466 static constexpr ArrowheadLUppConfig<T> store_config = [] {
467 ArrowheadLUppConfig<T> c = Config;
468 c.storage = ArrowheadStorage::Packed;
469 c.schur_in_place = false;
470 c.streamed_out_of_line = false;
471 return c;
472 }();
477 template<int capacity>
478 [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr ArrowheadShape
479 store_shape() noexcept {
480 return {B, S, capacity, R0};
481 }
489 template<int capacity, int nrhs>
490 [[nodiscard]] TDLS_HOST_DEVICE static constexpr int recomputed_size() noexcept {
491 if constexpr (capacity == 0) {
492 return 0;
493 } else {
494 using Store =
495 ArrowheadLUppSolverStatic<T, ArrowheadShape{B, S, capacity, R0}, store_config>;
496 return Store::template multirhs_workspace_size<nrhs> + Store::storage_size +
497 Store::dimension;
498 }
499 }
502 template<int capacity>
503 using RefusedStore =
504 ArrowheadLUppSolverStatic<T, ArrowheadShape{B, S, capacity, R0}, store_config>;
505
506 public:
512 using StreamedBlock = detail::ArrowheadStreamedBlock<T, S, B, col, diag>;
515 using RecomputedBlock = detail::ArrowheadRecomputedBlock<T, S, B, col, diag>;
516
546 template<bool internal_rhs, bool internal_work, bool internal_matrix,
547 int A_static_stride = dynamic_stride, int work_static_stride = dynamic_stride,
548 int rhs_static_stride = dynamic_stride, typename Provider>
549 [[nodiscard]] TDLS_HOST_DEVICE static constexpr bool
550 solve_inplace_streamed(Provider&& provider, T* TDLS_RESTRICT A, const int A_stride,
551 T* TDLS_RESTRICT work, const int work_stride, T* TDLS_RESTRICT y,
552 const int rhs_stride, int& refused_count) noexcept {
553 if constexpr (Config.streamed_out_of_line)
554 return solve_core_out_of_line<1, false, internal_rhs, internal_work, internal_matrix,
555 A_static_stride, work_static_stride, rhs_static_stride>(
556 A, A_stride, work, work_stride, y, rhs_stride, 0, 0, refused_count, &provider);
557 else
558 return solve_core<1, false, internal_rhs, internal_work, internal_matrix,
559 A_static_stride, work_static_stride, rhs_static_stride>(
560 A, A_stride, work, work_stride, y, rhs_stride, 0, 0, refused_count, &provider);
561 }
562
565 using StreamedOutcome = detail::ArrowheadStreamedOutcome;
566
618 template<int ntan, bool internal_rhs, bool internal_work, bool internal_matrix,
619 int A_static_stride = dynamic_stride, int work_static_stride = dynamic_stride,
620 int rhs_static_stride = dynamic_stride, int x_static_stride = dynamic_stride,
621 typename Provider, typename Converged, bool outlined = false>
622 [[nodiscard]] TDLS_HOST_DEVICE static constexpr StreamedOutcome
623 solve_inplace_streamed_or_tangent(Provider&& provider, Converged&& converged,
624 T* TDLS_RESTRICT A, const int A_stride, T* TDLS_RESTRICT work,
625 const int work_stride, T* TDLS_RESTRICT y,
626 const int rhs_stride, const int col0, T* TDLS_RESTRICT x,
627 const int x_stride, const std::size_t xcol_stride,
628 int& refused_count) noexcept {
629 static_assert(ntan >= 1 && ntan <= Shape.border,
630 "ArrowheadLUppSolverStatic: ntan must be in [1, border]");
631 TDLS_ARROW_LUPP_SIDED(template solve_inplace_streamed_or_tangent<
632 ntan, internal_rhs, internal_work, internal_matrix, A_static_stride,
633 work_static_stride, rhs_static_stride, x_static_stride, Provider,
634 Converged, outlined>(
635 std::forward<Provider>(provider), std::forward<Converged>(converged), A, A_stride, work,
636 work_stride, y, rhs_stride, col0, x, x_stride, xcol_stride, refused_count))
637 if constexpr (Config.streamed_out_of_line && !outlined)
638 return streamed_or_tangent_out_of_line<
639 ntan, internal_rhs, internal_work, internal_matrix, A_static_stride,
640 work_static_stride, rhs_static_stride, x_static_stride>(
641 provider, converged, A, A_stride, work, work_stride, y, rhs_stride, col0, x,
642 x_stride, xcol_stride, refused_count);
643 else
644 return streamed_or_tangent_core<ntan, internal_rhs, internal_work, internal_matrix,
645 A_static_stride, work_static_stride, rhs_static_stride,
646 x_static_stride>(provider, converged, A, A_stride, work,
647 work_stride, y, rhs_stride, col0, x,
648 x_stride, xcol_stride, refused_count);
649 }
650
681 template<int nrhs, bool internal_rhs, bool internal_work, bool internal_matrix,
682 int A_static_stride = dynamic_stride, int work_static_stride = dynamic_stride,
683 int rhs_static_stride = dynamic_stride, typename Provider>
684 [[nodiscard]] TDLS_HOST_DEVICE static constexpr bool solve_border_canonical_multirhs_streamed(
685 Provider&& provider, T* TDLS_RESTRICT A, const int A_stride, T* TDLS_RESTRICT work,
686 const int work_stride, const int col0, T* TDLS_RESTRICT x, const int rhs_stride,
687 const std::size_t xcol_stride, int& refused_count) noexcept {
688 static_assert(nrhs >= 1 && nrhs <= Shape.border,
689 "ArrowheadLUppSolverStatic: nrhs must be in [1, border]");
690 if constexpr (Config.streamed_out_of_line)
691 return solve_core_out_of_line<nrhs, true, internal_rhs, internal_work, internal_matrix,
692 A_static_stride, work_static_stride, rhs_static_stride>(
693 A, A_stride, work, work_stride, x, rhs_stride, xcol_stride, col0, refused_count,
694 &provider);
695 else
696 return solve_core<nrhs, true, internal_rhs, internal_work, internal_matrix,
697 A_static_stride, work_static_stride, rhs_static_stride>(
698 A, A_stride, work, work_stride, x, rhs_stride, xcol_stride, col0, refused_count,
699 &provider);
700 }
701
723 template<bool internal_rhs, bool internal_work, bool internal_matrix,
724 int A_static_stride = dynamic_stride, int work_static_stride = dynamic_stride,
725 int rhs_static_stride = dynamic_stride>
726 [[nodiscard]] TDLS_HOST_DEVICE static constexpr bool
727 solve_inplace(T* TDLS_RESTRICT A, const int A_stride, T* TDLS_RESTRICT work,
728 const int work_stride, T* TDLS_RESTRICT y, const int rhs_stride,
729 int& refused_count) noexcept {
730 return solve_core<1, false, internal_rhs, internal_work, internal_matrix, A_static_stride,
731 work_static_stride, rhs_static_stride>(A, A_stride, work, work_stride, y,
732 rhs_stride, 0, 0, refused_count);
733 }
734
752 template<bool internal_rhs, bool internal_work, bool internal_matrix,
753 int A_static_stride = dynamic_stride, int work_static_stride = dynamic_stride,
754 int rhs_static_stride = dynamic_stride>
755 [[nodiscard]] TDLS_HOST_DEVICE static constexpr bool
756 solve_inplace(T* TDLS_RESTRICT A, const int A_stride, T* TDLS_RESTRICT work,
757 const int work_stride, T* TDLS_RESTRICT y, const int rhs_stride) noexcept {
758 int refused_count = 0;
759 return solve_inplace<internal_rhs, internal_work, internal_matrix, A_static_stride,
760 work_static_stride, rhs_static_stride>(A, A_stride, work, work_stride,
761 y, rhs_stride, refused_count);
762 }
763
787 template<bool internal_rhs, bool internal_work, bool internal_matrix,
788 int A_static_stride = dynamic_stride, int work_static_stride = dynamic_stride,
789 int rhs_static_stride = dynamic_stride>
790 [[nodiscard]] TDLS_HOST_DEVICE static constexpr bool
791 solve(T* TDLS_RESTRICT A, const int A_stride, T* TDLS_RESTRICT work, const int work_stride,
792 const T* TDLS_RESTRICT b, T* TDLS_RESTRICT x, const int rhs_stride,
793 int& refused_count) noexcept {
794 const int stride = detail::stride_or<rhs_static_stride>(rhs_stride);
795 for (int i = 0; i < dimension; ++i) {
796 if constexpr (internal_rhs)
797 x[i] = b[i];
798 else
799 x[TDLS_STRIDED(i, stride)] = b[TDLS_STRIDED(i, stride)];
800 }
801 return solve_inplace<internal_rhs, internal_work, internal_matrix, A_static_stride,
802 work_static_stride, rhs_static_stride>(A, A_stride, work, work_stride,
803 x, rhs_stride, refused_count);
804 }
805
825 template<bool internal_rhs, bool internal_work, bool internal_matrix,
826 int A_static_stride = dynamic_stride, int work_static_stride = dynamic_stride,
827 int rhs_static_stride = dynamic_stride>
828 [[nodiscard]] TDLS_HOST_DEVICE static constexpr bool
829 solve(T* TDLS_RESTRICT A, const int A_stride, T* TDLS_RESTRICT work, const int work_stride,
830 const T* TDLS_RESTRICT b, T* TDLS_RESTRICT x, const int rhs_stride) noexcept {
831 int refused_count = 0;
832 return solve<internal_rhs, internal_work, internal_matrix, A_static_stride,
833 work_static_stride, rhs_static_stride>(A, A_stride, work, work_stride, b, x,
834 rhs_stride, refused_count);
835 }
836
869 template<int nrhs, bool internal_rhs, bool internal_work, bool internal_matrix,
870 int A_static_stride = dynamic_stride, int work_static_stride = dynamic_stride,
871 int rhs_static_stride = dynamic_stride>
872 [[nodiscard]] TDLS_HOST_DEVICE static constexpr bool
873 solve_inplace_multirhs(T* TDLS_RESTRICT A, const int A_stride, T* TDLS_RESTRICT work,
874 const int work_stride, T* TDLS_RESTRICT y, const int rhs_stride,
875 const std::size_t xcol_stride, int& refused_count) noexcept {
876 static_assert(nrhs >= 1, "ArrowheadLUppSolverStatic: nrhs must be >= 1");
877 return solve_core<nrhs, false, internal_rhs, internal_work, internal_matrix,
878 A_static_stride, work_static_stride, rhs_static_stride>(
879 A, A_stride, work, work_stride, y, rhs_stride, xcol_stride, 0, refused_count);
880 }
881
904 template<int nrhs, bool internal_rhs, bool internal_work, bool internal_matrix,
905 int A_static_stride = dynamic_stride, int work_static_stride = dynamic_stride,
906 int rhs_static_stride = dynamic_stride>
907 [[nodiscard]] TDLS_HOST_DEVICE static constexpr bool
908 solve_inplace_multirhs(T* TDLS_RESTRICT A, const int A_stride, T* TDLS_RESTRICT work,
909 const int work_stride, T* TDLS_RESTRICT y, const int rhs_stride,
910 const std::size_t xcol_stride) noexcept {
911 int refused_count = 0;
912 return solve_inplace_multirhs<nrhs, internal_rhs, internal_work, internal_matrix,
913 A_static_stride, work_static_stride, rhs_static_stride>(
914 A, A_stride, work, work_stride, y, rhs_stride, xcol_stride, refused_count);
915 }
916
944 template<int nrhs, bool internal_rhs, bool internal_work, bool internal_matrix,
945 int A_static_stride = dynamic_stride, int work_static_stride = dynamic_stride,
946 int rhs_static_stride = dynamic_stride>
947 [[nodiscard]] TDLS_HOST_DEVICE static constexpr bool
948 solve_multirhs(T* TDLS_RESTRICT A, const int A_stride, T* TDLS_RESTRICT work,
949 const int work_stride, const T* TDLS_RESTRICT b, T* TDLS_RESTRICT x,
950 const int rhs_stride, const std::size_t xcol_stride,
951 int& refused_count) noexcept {
952 const int stride = detail::stride_or<rhs_static_stride>(rhs_stride);
953 for (int m = 0; m < nrhs; ++m)
954 for (int i = 0; i < dimension; ++i) {
955 if constexpr (internal_rhs)
956 x[m * dimension + i] = b[m * dimension + i];
957 else
958 x[TDLS_STRIDED(i, stride) + TDLS_STRIDED_WIDE(m, xcol_stride)] =
959 b[TDLS_STRIDED(i, stride) + TDLS_STRIDED_WIDE(m, xcol_stride)];
960 }
961 return solve_inplace_multirhs<nrhs, internal_rhs, internal_work, internal_matrix,
962 A_static_stride, work_static_stride, rhs_static_stride>(
963 A, A_stride, work, work_stride, x, rhs_stride, xcol_stride, refused_count);
964 }
965
988 template<int nrhs, bool internal_rhs, bool internal_work, bool internal_matrix,
989 int A_static_stride = dynamic_stride, int work_static_stride = dynamic_stride,
990 int rhs_static_stride = dynamic_stride>
991 [[nodiscard]] TDLS_HOST_DEVICE static constexpr bool
992 solve_multirhs(T* TDLS_RESTRICT A, const int A_stride, T* TDLS_RESTRICT work,
993 const int work_stride, const T* TDLS_RESTRICT b, T* TDLS_RESTRICT x,
994 const int rhs_stride, const std::size_t xcol_stride) noexcept {
995 int refused_count = 0;
996 return solve_multirhs<nrhs, internal_rhs, internal_work, internal_matrix, A_static_stride,
997 work_static_stride, rhs_static_stride>(
998 A, A_stride, work, work_stride, b, x, rhs_stride, xcol_stride, refused_count);
999 }
1000
1043 template<int nrhs, bool internal_rhs, bool internal_work, bool internal_matrix,
1044 int A_static_stride = dynamic_stride, int work_static_stride = dynamic_stride,
1045 int rhs_static_stride = dynamic_stride>
1046 [[nodiscard]] TDLS_HOST_DEVICE static constexpr bool
1047 solve_border_canonical_multirhs(T* TDLS_RESTRICT A, const int A_stride, T* TDLS_RESTRICT work,
1048 const int work_stride, const int col0, T* TDLS_RESTRICT x,
1049 const int rhs_stride, const std::size_t xcol_stride,
1050 int& refused_count) noexcept {
1051 static_assert(nrhs >= 1 && nrhs <= Shape.border,
1052 "ArrowheadLUppSolverStatic: nrhs must be in [1, border]");
1053 return solve_core<nrhs, true, internal_rhs, internal_work, internal_matrix, A_static_stride,
1054 work_static_stride, rhs_static_stride>(
1055 A, A_stride, work, work_stride, x, rhs_stride, xcol_stride, col0, refused_count);
1056 }
1057
1082 template<int nrhs, bool internal_rhs, bool internal_work, bool internal_matrix,
1083 int A_static_stride = dynamic_stride, int work_static_stride = dynamic_stride,
1084 int rhs_static_stride = dynamic_stride>
1085 [[nodiscard]] TDLS_HOST_DEVICE static constexpr bool
1086 solve_border_canonical_multirhs(T* TDLS_RESTRICT A, const int A_stride, T* TDLS_RESTRICT work,
1087 const int work_stride, const int col0, T* TDLS_RESTRICT x,
1088 const int rhs_stride, const std::size_t xcol_stride) noexcept {
1089 int refused_count = 0;
1090 return solve_border_canonical_multirhs<nrhs, internal_rhs, internal_work, internal_matrix,
1091 A_static_stride, work_static_stride,
1092 rhs_static_stride>(
1093 A, A_stride, work, work_stride, col0, x, rhs_stride, xcol_stride, refused_count);
1094 }
1095
1098 using RecomputedOutcome = detail::ArrowheadRecomputedOutcome;
1099
1106 template<int capacity, int nrhs = 1>
1107 static constexpr int recomputed_workspace_size = recomputed_size<capacity, nrhs>();
1108
1166 template<int capacity, bool internal_rhs, bool internal_work, bool internal_matrix,
1167 int Abb_static_stride = dynamic_stride, int work_static_stride = dynamic_stride,
1168 int rhs_static_stride = dynamic_stride, typename Provider>
1169 [[nodiscard]] TDLS_HOST_DEVICE static constexpr RecomputedOutcome
1170 solve_inplace_recomputed(Provider&& provider, const T* TDLS_RESTRICT Abb, const int Abb_stride,
1171 T* TDLS_RESTRICT work, const int work_stride, T* TDLS_RESTRICT y,
1172 const int rhs_stride, int& refused_count) noexcept {
1173 return recomputed_core<0, capacity, internal_rhs, internal_work, internal_matrix,
1174 Abb_static_stride, work_static_stride, rhs_static_stride,
1175 dynamic_stride>(provider, Abb, Abb_stride, work, work_stride, y,
1176 rhs_stride, 0, nullptr, 1, 0, refused_count);
1177 }
1178
1230 template<int ntan, int capacity, bool internal_rhs, bool internal_work, bool internal_matrix,
1231 int Abb_static_stride = dynamic_stride, int work_static_stride = dynamic_stride,
1232 int rhs_static_stride = dynamic_stride, int x_static_stride = dynamic_stride,
1233 typename Provider, typename Converged>
1234 [[nodiscard]] TDLS_HOST_DEVICE static constexpr RecomputedOutcome
1235 solve_inplace_recomputed_or_tangent(Provider&& provider, Converged&& converged,
1236 const T* TDLS_RESTRICT Abb, const int Abb_stride,
1237 T* TDLS_RESTRICT work, const int work_stride,
1238 T* TDLS_RESTRICT y, const int rhs_stride, const int col0,
1239 T* TDLS_RESTRICT x, const int x_stride,
1240 const std::size_t xcol_stride,
1241 int& refused_count) noexcept {
1242 static_assert(ntan >= 1 && ntan <= Shape.border,
1243 "ArrowheadLUppSolverStatic: ntan must be in [1, border]");
1244 return recomputed_core<ntan, capacity, internal_rhs, internal_work, internal_matrix,
1245 Abb_static_stride, work_static_stride, rhs_static_stride,
1246 x_static_stride>(provider, Abb, Abb_stride, work, work_stride, y,
1247 rhs_stride, col0, x, x_stride, xcol_stride,
1248 refused_count, &converged);
1249 }
1250
1251 private:
1252// gcc cannot always prove that the provider writes every entry of the
1253// streamed block before it is read, and reports the reads through
1254// -Wmaybe-uninitialized at -O3 (see condense_or_store_streamed): the report
1255// is spurious. The suppression is scoped to this function and to that
1256// warning only.
1257#if defined(TDLS_DETAIL_COMPILER_GCC)
1258#pragma GCC diagnostic push
1259#pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
1260#endif
1294 template<int ntan, int capacity, bool internal_rhs, bool internal_work, bool internal_matrix,
1295 int Abb_static_stride = dynamic_stride, int work_static_stride = dynamic_stride,
1296 int rhs_static_stride = dynamic_stride, int x_static_stride = dynamic_stride,
1297 typename Provider, typename... Converged>
1298 [[nodiscard]] TDLS_HOST_DEVICE static constexpr RecomputedOutcome
1299 recomputed_core(Provider& provider, const T* TDLS_RESTRICT Abb, const int Abb_stride_arg,
1300 [[maybe_unused]] T* TDLS_RESTRICT work,
1301 [[maybe_unused]] const int work_stride_arg, T* TDLS_RESTRICT y,
1302 const int rhs_stride_arg, [[maybe_unused]] const int col0,
1303 [[maybe_unused]] T* TDLS_RESTRICT x, [[maybe_unused]] const int x_stride_arg,
1304 [[maybe_unused]] const std::size_t xcol_stride, int& refused_count,
1305 Converged*... converged) noexcept {
1306 const int Abb_stride = detail::stride_or<Abb_static_stride>(Abb_stride_arg);
1307 [[maybe_unused]] const int work_stride =
1308 detail::stride_or<work_static_stride>(work_stride_arg);
1309 const int rhs_stride = detail::stride_or<rhs_static_stride>(rhs_stride_arg);
1310 [[maybe_unused]] const int x_stride = detail::stride_or<x_static_stride>(x_stride_arg);
1311 TDLS_ARROW_LUPP_SIDED(
1312 template recomputed_core<ntan, capacity, internal_rhs, internal_work, internal_matrix,
1313 Abb_static_stride, work_static_stride, rhs_static_stride,
1314 x_static_stride>(provider, Abb, Abb_stride, work, work_stride,
1315 y, rhs_stride, col0, x, x_stride, xcol_stride,
1316 refused_count, converged...))
1317 static_assert(capacity >= 0 && capacity <= NB,
1318 "ArrowheadLUppSolverStatic: the capacity of the store of refused blocks "
1319 "must be in [0, blocks]");
1320 static_assert(!in_place, "ArrowheadLUppSolverStatic: the recomputed solve accumulates the "
1321 "Schur complement apart (schur_in_place off)");
1322 // Columns of the front reserved in the workspace: the tangent
1323 // columns of the canonical second pass, one at least.
1324 constexpr int fnc = ntan > 1 ? ntan : 1;
1325 // Rows R0 .. R0 + R of the Schur complement and of its right-hand
1326 // side.
1327 T Sacc[R * B];
1328 T racc[R];
1329 TDLS_ARROW_LUPP_FOR(r, 0, R,
1330 TDLS_ARROW_LUPP_FOR(c, 0, B,
1331 Sacc[r * B + c] = border_entry<internal_matrix>(
1332 Abb, Abb_stride, R0 + r, c);) racc[r] =
1333 TDLS_ARROW_LUPP_Y(R0 + r);)
1334 constexpr int words = (NB + 31) / 32;
1335 unsigned refused[words];
1336 TDLS_ARROW_LUPP_FOR(e, 0, words, refused[e] = 0u;)
1337 int count = 0;
1338 refused_count = 0;
1339#if defined(TDLS_DETAIL_CODEGEN_NVCXX)
1340 // One block for the condensation and the recovery, safe since the
1341 // provider writes every entry of the block. With one block per loop,
1342 // nvc++ 25.11 at -O3 generates device code that reads the local
1343 // memory out of bounds with the dense coupling.
1344 RecomputedBlock shared_blk;
1345#endif
1346 for (int k = 0; k < NB; ++k) {
1347#if defined(TDLS_DETAIL_CODEGEN_NVCXX)
1348 RecomputedBlock& blk = shared_blk;
1349#else
1350 RecomputedBlock blk;
1351#endif
1352 blk.recovery = false;
1353 provider(k, blk);
1354 // A save of the block, to keep it as given if it is refused: D_k,
1355 // and [A_kb | y_k] when the check reads the forward stage.
1356 constexpr bool forward_check = forward_check_for<1>;
1357 T Dsave[S * S];
1358 T Wsave[forward_check ? (B + 1) * S : 1];
1359 TDLS_ARROW_LUPP_FOR(e, 0, S * S, Dsave[e] = blk.Dl[e];)
1360 if constexpr (forward_check) {
1361 TDLS_ARROW_LUPP_FOR(e, 0, (B + 1) * S, Wsave[e] = blk.W[e];)
1362 } else {
1363 Wsave[0] = T(0);
1364 }
1365 if (condense_streamed<1, true, true, true>(nullptr, 1, nullptr, 1, 0, k, Sacc, racc,
1366 blk))
1367 continue;
1368 // Refused: kept in the store as the provider gave it, or the
1369 // solve is deferred.
1370 ++count;
1371 refused_count = count;
1372 refused[k >> 5] |= 1u << (k & 31);
1373 if constexpr (capacity == 0) {
1374 return RecomputedOutcome::deferred;
1375 } else {
1376 if (count > capacity) return RecomputedOutcome::deferred;
1377 TDLS_ARROW_LUPP_FOR(e, 0, S * S, blk.Dl[e] = Dsave[e];)
1378 if constexpr (forward_check) {
1379 TDLS_ARROW_LUPP_FOR(e, 0, (B + 1) * S, blk.W[e] = Wsave[e];)
1380 }
1381 keep_refused<capacity, fnc, internal_work>(work, work_stride, count - 1, blk);
1382 }
1383 }
1384 // The border system: the accumulated rows, the others as given.
1385 T Sbb[B * B];
1386 TDLS_ARROW_LUPP_FOR(
1387 r, 0, B, const bool acc = r >= R0 && r < R0 + R; TDLS_ARROW_LUPP_FOR(
1388 c, 0, B,
1389 Sbb[r * B + c] = acc ? Sacc[(acc ? r - R0 : 0) * B + c]
1390 : border_entry<internal_matrix>(Abb, Abb_stride, r, c);))
1391 if constexpr (sizeof...(Converged) > 0) {
1392 if ((*streamed_arg(converged...))()) {
1393 // The tangent: unit columns, the refused blocks of the store
1394 // by the second pass of the canonical solve.
1395 T rt[B * fnc];
1396 TDLS_ARROW_LUPP_FOR(
1397 m, 0, fnc,
1398 TDLS_ARROW_LUPP_FOR(r, 0, B, rt[m * B + r] = r == col0 + m ? T(1) : T(0);))
1399 if constexpr (capacity > 0) {
1400 if (count > 0 &&
1401 !second_pass_refused<capacity, fnc, fnc, true, internal_work,
1402 work_static_stride>(work, work_stride, count, Sbb, rt))
1403 return RecomputedOutcome::singular;
1404 }
1405 int pt[B];
1406 bool ok;
1407 if constexpr (fnc == 1)
1408 ok = BorderSolver::template solve_inplace<true, true, true>(Sbb, 1, pt, 1, rt,
1409 1);
1410 else
1411 ok = BorderSolver::template solve_inplace_multirhs<fnc, true, true, true>(
1412 Sbb, 1, pt, 1, rt, 1, B);
1413 if (!ok) return RecomputedOutcome::singular;
1414 TDLS_ARROW_LUPP_FOR(
1415 m, 0, fnc,
1416 TDLS_ARROW_LUPP_FOR(r, 0, B,
1417 x[internal_rhs ? unsigned(m * B + r)
1418 : TDLS_STRIDED(r, x_stride) +
1419 TDLS_STRIDED_WIDE(m, xcol_stride)] =
1420 rt[m * B + r];))
1421 return RecomputedOutcome::tangent;
1422 }
1423 }
1424 T rb[B];
1425 TDLS_ARROW_LUPP_FOR(r, 0, B, const bool acc = r >= R0 && r < R0 + R;
1426 rb[r] = acc ? racc[acc ? r - R0 : 0] : TDLS_ARROW_LUPP_Y(r);)
1427 if constexpr (capacity > 0) {
1428 if (count > 0 &&
1429 !second_pass_refused<capacity, fnc, 1, false, internal_work, work_static_stride>(
1430 work, work_stride, count, Sbb, rb))
1431 return RecomputedOutcome::singular;
1432 }
1433 int pb[B];
1434 if (!BorderSolver::template solve_inplace<true, true, true>(Sbb, 1, pb, 1, rb, 1))
1435 return RecomputedOutcome::singular;
1436 // The accepted blocks, recomputed, the last one first; then the
1437 // refused ones from the store; then the border.
1438 for (int k = NB - 1; k >= 0; --k) {
1439 if ((refused[k >> 5] >> (k & 31)) & 1u) continue;
1440#if defined(TDLS_DETAIL_CODEGEN_NVCXX)
1441 RecomputedBlock& blk = shared_blk;
1442#else
1443 RecomputedBlock blk;
1444#endif
1445 blk.recovery = true;
1446 provider(k, blk);
1447 T xk[S];
1448 TDLS_ARROW_LUPP_FOR(i, 0, S, xk[i] = blk.W[B * S + i];
1449 TDLS_ARROW_LUPP_FOR(c, 0, B, xk[i] -= blk.W[c * S + i] * rb[c];))
1450 recover_block(blk.Dl, xk);
1451 TDLS_ARROW_LUPP_FOR(i, 0, S, TDLS_ARROW_LUPP_Y(B + k * S + i) = xk[i];)
1452 }
1453 if constexpr (capacity > 0) {
1454 if (count > 0)
1455 recover_refused<capacity, fnc, internal_rhs, internal_work, work_static_stride,
1456 rhs_static_stride>(work, work_stride, y, rhs_stride, refused, count,
1457 rb);
1458 }
1459 TDLS_ARROW_LUPP_FOR(r, 0, B, TDLS_ARROW_LUPP_Y(r) = rb[r];)
1460 return RecomputedOutcome::solved;
1461 }
1462#if defined(TDLS_DETAIL_COMPILER_GCC)
1463#pragma GCC diagnostic pop
1464#endif
1465
1474 template<bool internal_matrix>
1475 [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr T
1476 border_entry(const T* TDLS_RESTRICT Abb, const int Abb_stride, const int r,
1477 const int c) noexcept {
1478 const int e = col ? c * B + r : r * B + c;
1479 return Abb[internal_matrix ? unsigned(e) : TDLS_STRIDED(e, Abb_stride)];
1480 }
1481
1493 template<int capacity, int fnc, bool internal_work>
1494 TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
1495 keep_refused(T* TDLS_RESTRICT work, const int work_stride, const int j,
1496 const StreamedBlock& blk) noexcept {
1497 using Store = RefusedStore<capacity>;
1498 using Map = typename Store::storage_map;
1499 constexpr int off = Store::template multirhs_workspace_size<fnc>;
1500 constexpr int oy = off + Store::storage_size;
1501 TDLS_ARROW_LUPP_FOR(
1502 i, 0, S,
1503 TDLS_ARROW_LUPP_FOR(j2, 0, S,
1504 TDLS_ARROW_LUPP_W(off + Map::d(store_shape<capacity>(), j, i, j2)) =
1505 blk.Dl[col ? j2 * S + i : i * S + j2];)
1506 TDLS_ARROW_LUPP_FOR(c, 0, B,
1507 TDLS_ARROW_LUPP_W(off + Map::akb(store_shape<capacity>(), j, i,
1508 c)) = blk.W[c * S + i];)
1509 TDLS_ARROW_LUPP_W(oy + B + j * S + i) = blk.W[B * S + i];)
1510 if constexpr (diag) {
1511 TDLS_ARROW_LUPP_FOR(
1512 i, 0, S,
1513 TDLS_ARROW_LUPP_W(off + Map::abk_diag(store_shape<capacity>(), j, i)) = blk.Abk[i];)
1514 } else {
1515 TDLS_ARROW_LUPP_FOR(
1516 r, 0, B,
1517 TDLS_ARROW_LUPP_FOR(i, 0, S,
1518 TDLS_ARROW_LUPP_W(off + Map::abk(store_shape<capacity>(), j, r,
1519 i)) = blk.Abk[r * S + i];))
1520 }
1521 }
1522
1542 template<int capacity, int fnc, int nc, bool canonical, bool internal_work,
1543 int work_static_stride = dynamic_stride>
1544 TDLS_HOST_DEVICE static constexpr bool
1545 second_pass_refused(T* TDLS_RESTRICT work, const int work_stride_arg, const int count,
1546 T* TDLS_RESTRICT Sbb, T* TDLS_RESTRICT rb) noexcept {
1547 const int work_stride = detail::stride_or<work_static_stride>(work_stride_arg);
1548 using Store = RefusedStore<capacity>;
1549 constexpr int off = Store::template multirhs_workspace_size<fnc>;
1550 constexpr int oy = off + Store::storage_size;
1551 constexpr int words = (capacity + 31) / 32;
1552 unsigned held[words];
1553 for (int e = 0; e < words; ++e)
1554 held[e] = 0u;
1555 for (int j = 0; j < count; ++j)
1556 held[j >> 5] |= 1u << (j & 31);
1557 T* const store = work + (internal_work ? off : TDLS_STRIDED(off, work_stride));
1558 T* const ystore = work + (internal_work ? oy : TDLS_STRIDED(oy, work_stride));
1559 return Store::template second_pass<nc, canonical, internal_work, internal_work,
1560 internal_work, work_static_stride, work_static_stride,
1561 work_static_stride>(
1562 store, work_stride, work, work_stride, ystore, work_stride, 0, held, Sbb, rb,
1563 BlockSolver::singular_floor);
1564 }
1565
1566// The block the provider wrote for the recovery (see the suppression
1567// before condense_streamed): scoped to recover_block and to that warning.
1568#if defined(TDLS_DETAIL_COMPILER_GCC)
1569#pragma GCC diagnostic push
1570#pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
1571#endif
1578 TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
1579 recover_block(T* TDLS_RESTRICT Dl, T* TDLS_RESTRICT x) noexcept {
1580 int pl[S];
1581 if constexpr (reg) {
1582 (void)register_factor(Dl, pl);
1583 register_solve<1, true>(Dl, pl, x);
1584 } else {
1585 (void)BlockSolver::template factorize<true, true>(Dl, 1, pl, 1);
1586 BlockSolver::template substitute_inplace<true, true, true>(Dl, 1, pl, 1, x, 1);
1587 }
1588 }
1589
1590#if defined(TDLS_DETAIL_COMPILER_GCC)
1591#pragma GCC diagnostic pop
1592#endif
1593
1613 template<int capacity, int fnc, bool internal_rhs, bool internal_work,
1614 int work_static_stride = dynamic_stride, int rhs_static_stride = dynamic_stride>
1615 TDLS_HOST_DEVICE static constexpr void
1616 recover_refused(const T* TDLS_RESTRICT work, const int work_stride_arg, T* TDLS_RESTRICT y,
1617 const int rhs_stride_arg, const unsigned* TDLS_RESTRICT refused,
1618 const int count, const T* TDLS_RESTRICT rb) noexcept {
1619 const int work_stride = detail::stride_or<work_static_stride>(work_stride_arg);
1620 const int rhs_stride = detail::stride_or<rhs_static_stride>(rhs_stride_arg);
1621 using Store = RefusedStore<capacity>;
1622 using Map = typename Store::storage_map;
1623 constexpr int off = Store::template multirhs_workspace_size<fnc>;
1624 constexpr int oy = off + Store::storage_size;
1625 constexpr int fcols = front_cols_for<1>;
1626 T w[R];
1627 for (int a = 0; a < R; ++a)
1628 w[a] = T(0);
1629 int j = count;
1630 for (int k = NB - 1; k >= 0; --k) {
1631 if (!((refused[k >> 5] >> (k & 31)) & 1u)) continue;
1632 --j;
1633 T s[S];
1634 for (int i = 0; i < S; ++i) {
1635 s[i] = TDLS_ARROW_LUPP_W(oy + B + j * S + i);
1636 for (int c = 0; c < B; ++c)
1637 s[i] -=
1638 TDLS_ARROW_LUPP_W(off + Map::akb(store_shape<capacity>(), j, i, c)) * rb[c];
1639 for (int a = 0; a < R; ++a) {
1640 const T h =
1641 a < S ? TDLS_ARROW_LUPP_W(off + Map::d(store_shape<capacity>(), j, i, a))
1642 : TDLS_ARROW_LUPP_W(front_rows * fcols + (j * S + i) * extra_cols +
1643 a - S);
1644 s[i] -= h * w[a];
1645 }
1646 }
1647 if constexpr (diag) {
1648 for (int i = 0; i < S; ++i)
1649 w[i] += TDLS_ARROW_LUPP_W(off + Map::abk_diag(store_shape<capacity>(), j, i)) *
1650 s[i];
1651 } else {
1652 for (int r = 0; r < B; ++r)
1653 for (int i = 0; i < S; ++i)
1654 w[r] +=
1655 TDLS_ARROW_LUPP_W(off + Map::abk(store_shape<capacity>(), j, r, i)) *
1656 s[i];
1657 }
1658 for (int i = 0; i < S; ++i)
1659 TDLS_ARROW_LUPP_Y(B + k * S + i) = s[i];
1660 }
1661 }
1662
1663 private:
1696 template<int ntan, bool internal_rhs, bool internal_work, bool internal_matrix,
1697 int A_static_stride = dynamic_stride, int work_static_stride = dynamic_stride,
1698 int rhs_static_stride = dynamic_stride, int x_static_stride = dynamic_stride,
1699 typename Provider, typename Converged>
1700 TDLS_HOST_DEVICE TDLS_FORCE_INLINE_ALL static constexpr StreamedOutcome
1701 streamed_or_tangent_core(Provider& provider, Converged& converged, T* TDLS_RESTRICT A,
1702 const int A_stride_arg, T* TDLS_RESTRICT work,
1703 const int work_stride_arg, T* TDLS_RESTRICT y,
1704 const int rhs_stride_arg, const int col0, T* TDLS_RESTRICT x,
1705 const int x_stride_arg, const std::size_t xcol_stride,
1706 int& refused_count) noexcept {
1707 const int A_stride = detail::stride_or<A_static_stride>(A_stride_arg);
1708 const int work_stride = detail::stride_or<work_static_stride>(work_stride_arg);
1709 const int rhs_stride = detail::stride_or<rhs_static_stride>(rhs_stride_arg);
1710 const int x_stride = detail::stride_or<x_static_stride>(x_stride_arg);
1711 constexpr int ny = 1;
1712 [[maybe_unused]] constexpr int ycols = 1;
1713 T Sacc[in_place ? 1 : R * B];
1714 T racc[in_place ? 1 : R];
1715 if constexpr (in_place) {
1716 Sacc[0] = T(0);
1717 racc[0] = T(0);
1718 } else {
1719 TDLS_ARROW_LUPP_FOR(
1720 r, 0, R,
1721 TDLS_ARROW_LUPP_FOR(c, 0, B,
1722 Sacc[r * B + c] = TDLS_ARROW_LUPP_A(
1723 storage_map::abb(shape_value(), R0 + r, c));) racc[r] =
1724 TDLS_ARROW_LUPP_YW(0, R0 + r);)
1725 }
1726 constexpr int words = (NB + 31) / 32;
1727 unsigned refused[words];
1728 TDLS_ARROW_LUPP_FOR(e, 0, words, refused[e] = 0u;)
1729 refused_count = 0;
1730 for (int k = 0; k < NB; ++k) {
1731 StreamedBlock blk;
1732 provider(k, blk);
1733 if (!condense_or_store_streamed<ny, internal_rhs, internal_matrix>(
1734 A, A_stride, y, rhs_stride, 0, k, Sacc, racc, blk)) {
1735 refused[k >> 5] |= 1u << (k & 31);
1736 ++refused_count;
1737 }
1738 }
1739 T Sbb[B * B];
1740 int pb[B];
1741 // The rows as stored and the accumulated rows in separate loops, as
1742 // in solve_core.
1743 constexpr int acc0 = in_place ? B : R0;
1744 constexpr int acc1 = in_place ? B : R0 + R;
1745 TDLS_ARROW_LUPP_FOR(r, 0, acc0,
1746 TDLS_ARROW_LUPP_FOR(c, 0, B,
1747 Sbb[r * B + c] = TDLS_ARROW_LUPP_A(
1748 storage_map::abb(shape_value(), r, c));))
1749 TDLS_ARROW_LUPP_FOR(r, acc0, acc1,
1750 TDLS_ARROW_LUPP_FOR(c, 0, B, Sbb[r * B + c] = Sacc[(r - R0) * B + c];))
1751 TDLS_ARROW_LUPP_FOR(r, acc1, B,
1752 TDLS_ARROW_LUPP_FOR(c, 0, B,
1753 Sbb[r * B + c] = TDLS_ARROW_LUPP_A(
1754 storage_map::abb(shape_value(), r, c));))
1755 if (converged()) {
1756 T rt[B * ntan];
1757 TDLS_ARROW_LUPP_FOR(
1758 m, 0, ntan,
1759 TDLS_ARROW_LUPP_FOR(r, 0, B, rt[m * B + r] = r == col0 + m ? T(1) : T(0);))
1760 // the refused blocks, stored as is, by the second pass of the
1761 // canonical solve, which writes neither x nor y
1762 if (refused_count > 0 &&
1763 !second_pass<ntan, true, internal_rhs, internal_work, internal_matrix,
1764 A_static_stride, work_static_stride, rhs_static_stride>(
1765 A, A_stride, work, work_stride, y, rhs_stride, 0, refused, Sbb, rt,
1766 BlockSolver::singular_floor))
1767 return StreamedOutcome::singular;
1768 bool ok;
1769 if constexpr (ntan == 1)
1770 ok = BorderSolver::template solve_inplace<true, true, true>(Sbb, 1, pb, 1, rt, 1);
1771 else
1772 ok = BorderSolver::template solve_inplace_multirhs<ntan, true, true, true>(
1773 Sbb, 1, pb, 1, rt, 1, B);
1774 if (!ok) return StreamedOutcome::singular;
1775 TDLS_ARROW_LUPP_FOR(
1776 m, 0, ntan,
1777 TDLS_ARROW_LUPP_FOR(r, 0, B,
1778 x[internal_rhs ? unsigned(m * B + r)
1779 : TDLS_STRIDED(r, x_stride) +
1780 TDLS_STRIDED_WIDE(m, xcol_stride)] =
1781 rt[m * B + r];))
1782 return StreamedOutcome::tangent;
1783 }
1784 T rb[B];
1785 TDLS_ARROW_LUPP_FOR(r, 0, acc0, rb[r] = TDLS_ARROW_LUPP_YW(0, r);)
1786 TDLS_ARROW_LUPP_FOR(r, acc0, acc1, rb[r] = racc[r - R0];)
1787 TDLS_ARROW_LUPP_FOR(r, acc1, B, rb[r] = TDLS_ARROW_LUPP_YW(0, r);)
1788 if (refused_count > 0 &&
1789 !second_pass<1, false, internal_rhs, internal_work, internal_matrix, A_static_stride,
1790 work_static_stride, rhs_static_stride>(A, A_stride, work, work_stride, y,
1791 rhs_stride, 0, refused, Sbb, rb,
1792 BlockSolver::singular_floor))
1793 return StreamedOutcome::singular;
1794 if (!BorderSolver::template solve_inplace<true, true, true>(Sbb, 1, pb, 1, rb, 1))
1795 return StreamedOutcome::singular;
1796 TDLS_ARROW_LUPP_FOR(r, 0, B, TDLS_ARROW_LUPP_YW(0, r) = rb[r];)
1797 recover<1, internal_rhs, internal_work, internal_matrix>(A, A_stride, work, work_stride, y,
1798 rhs_stride, 0, refused, rb);
1799 return StreamedOutcome::solved;
1800 }
1801
1805 template<int ntan, bool internal_rhs, bool internal_work, bool internal_matrix,
1806 int A_static_stride, int work_static_stride, int rhs_static_stride,
1807 int x_static_stride, typename Provider, typename Converged>
1808 TDLS_HOST_DEVICE TDLS_NO_INLINE_GPU static constexpr StreamedOutcome
1809 streamed_or_tangent_out_of_line(Provider& provider, Converged& converged, T* TDLS_RESTRICT A,
1810 const int A_stride, T* TDLS_RESTRICT work,
1811 const int work_stride, T* TDLS_RESTRICT y, const int rhs_stride,
1812 const int col0, T* TDLS_RESTRICT x, const int x_stride,
1813 const std::size_t xcol_stride, int& refused_count) noexcept {
1814 return solve_inplace_streamed_or_tangent<
1815 ntan, internal_rhs, internal_work, internal_matrix, A_static_stride, work_static_stride,
1816 rhs_static_stride, x_static_stride, Provider&, Converged&, true>(
1817 provider, converged, A, A_stride, work, work_stride, y, rhs_stride, col0, x, x_stride,
1818 xcol_stride, refused_count);
1819 }
1820
1824 template<int nc, bool canonical, bool internal_rhs, bool internal_work, bool internal_matrix,
1825 int A_static_stride, int work_static_stride, int rhs_static_stride, typename Provider>
1826 TDLS_HOST_DEVICE TDLS_NO_INLINE_GPU static constexpr bool
1827 solve_core_out_of_line(T* TDLS_RESTRICT A, const int A_stride, T* TDLS_RESTRICT work,
1828 const int work_stride, T* TDLS_RESTRICT y, const int rhs_stride,
1829 const std::size_t xcol_stride, const int col0, int& refused_count,
1830 Provider* provider) noexcept {
1831 return solve_core<nc, canonical, internal_rhs, internal_work, internal_matrix,
1832 A_static_stride, work_static_stride, rhs_static_stride>(
1833 A, A_stride, work, work_stride, y, rhs_stride, xcol_stride, col0, refused_count,
1834 provider);
1835 }
1836
1844 template<typename P>
1845 [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_ALL static constexpr P*
1846 streamed_arg(P* p) noexcept {
1847 return p;
1848 }
1849
1882 template<int nc, bool canonical, bool internal_rhs, bool internal_work, bool internal_matrix,
1883 int A_static_stride = dynamic_stride, int work_static_stride = dynamic_stride,
1884 int rhs_static_stride = dynamic_stride, typename... Provider>
1885 TDLS_HOST_DEVICE TDLS_FORCE_INLINE_ALL static constexpr bool
1886 solve_core(T* TDLS_RESTRICT A, const int A_stride_arg, T* TDLS_RESTRICT work,
1887 const int work_stride_arg, T* TDLS_RESTRICT y, const int rhs_stride_arg,
1888 [[maybe_unused]] const std::size_t xcol_stride, [[maybe_unused]] const int col0,
1889 int& refused_count, Provider*... provider) noexcept {
1890 const int A_stride = detail::stride_or<A_static_stride>(A_stride_arg);
1891 const int work_stride = detail::stride_or<work_static_stride>(work_stride_arg);
1892 const int rhs_stride = detail::stride_or<rhs_static_stride>(rhs_stride_arg);
1893 TDLS_ARROW_LUPP_SIDED(
1894 template solve_core<nc, canonical, internal_rhs, internal_work, internal_matrix,
1895 A_static_stride, work_static_stride, rhs_static_stride>(
1896 A, A_stride, work, work_stride, y, rhs_stride, xcol_stride, col0, refused_count,
1897 provider...))
1898 constexpr bool streamed = sizeof...(Provider) > 0;
1899 // Right-hand sides carried through the blocks.
1900 constexpr int ny = canonical ? 0 : nc;
1901 // Columns of y, for TDLS_ARROW_LUPP_YW.
1902 [[maybe_unused]] constexpr int ycols = nc;
1903 // Rows R0 .. R0 + R of the Schur complement and of its right-hand
1904 // sides (column m at m * R), unless accumulated in place.
1905 T Sacc[in_place ? 1 : R * B];
1906 T racc[in_place || ny == 0 ? 1 : R * ny];
1907 if constexpr (in_place) {
1908 // Placeholders, never read: initialized for the flow analysis
1909 // of the compilers.
1910 Sacc[0] = T(0);
1911 racc[0] = T(0);
1912 } else {
1913 if constexpr (ny == 0) racc[0] = T(0);
1914 TDLS_ARROW_LUPP_FOR(
1915 r, 0, R,
1916 TDLS_ARROW_LUPP_FOR(c, 0, B,
1917 Sacc[r * B + c] = TDLS_ARROW_LUPP_A(
1918 storage_map::abb(shape_value(), R0 + r, c));)
1919 TDLS_ARROW_LUPP_FOR(m, 0, ny, racc[m * R + r] = TDLS_ARROW_LUPP_YW(m, R0 + r);))
1920 }
1921 constexpr int words = (NB + 31) / 32;
1922 unsigned refused[words];
1923 TDLS_ARROW_LUPP_FOR(e, 0, words, refused[e] = 0u;)
1924 refused_count = 0;
1925 for (int k = 0; k < NB; ++k) {
1926 if constexpr (streamed) {
1927 static_assert(
1928 ny <= 1,
1929 "ArrowheadLUppSolverStatic: streamed blocks carry one right-hand side");
1930 StreamedBlock blk;
1931 (*streamed_arg(provider...))(k, blk);
1932 if (!condense_or_store_streamed<ny, internal_rhs, internal_matrix>(
1933 A, A_stride, y, rhs_stride, xcol_stride, k, Sacc, racc, blk)) {
1934 refused[k >> 5] |= 1u << (k & 31);
1935 ++refused_count;
1936 }
1937 } else {
1938 if (!condense<ny, internal_rhs, internal_matrix>(A, A_stride, y, rhs_stride,
1939 xcol_stride, k, Sacc, racc)) {
1940 refused[k >> 5] |= 1u << (k & 31);
1941 ++refused_count;
1942 }
1943 }
1944 }
1945 // The border system: the accumulated rows, the others as stored;
1946 // column m of its right-hand sides at m * B.
1947 T Sbb[B * B];
1948 T rb[B * nc];
1949 // The rows as stored and the accumulated rows in separate loops:
1950 // GCC 11 and 12 vectorize a choice of source per entry wrongly (-O3
1951 // with AVX2 or AVX-512, unroll_loops off).
1952 constexpr int acc0 = in_place ? B : R0;
1953 constexpr int acc1 = in_place ? B : R0 + R;
1954 TDLS_ARROW_LUPP_FOR(
1955 r, 0, acc0,
1956 TDLS_ARROW_LUPP_FOR(
1957 c, 0, B, Sbb[r * B + c] = TDLS_ARROW_LUPP_A(storage_map::abb(shape_value(), r, c));)
1958 TDLS_ARROW_LUPP_FOR(
1959 m, 0, nc, if constexpr (canonical) rb[m * B + r] = r == col0 + m ? T(1) : T(0);
1960 else rb[m * B + r] = TDLS_ARROW_LUPP_YW(m, r);))
1961 TDLS_ARROW_LUPP_FOR(r, acc0, acc1,
1962 TDLS_ARROW_LUPP_FOR(c, 0, B, Sbb[r * B + c] = Sacc[(r - R0) * B + c];)
1963 TDLS_ARROW_LUPP_FOR(m, 0, nc,
1964 if constexpr (canonical) rb[m * B + r] =
1965 r == col0 + m ? T(1) : T(0);
1966 else rb[m * B + r] = racc[m * R + r - R0];))
1967 TDLS_ARROW_LUPP_FOR(
1968 r, acc1, B,
1969 TDLS_ARROW_LUPP_FOR(
1970 c, 0, B, Sbb[r * B + c] = TDLS_ARROW_LUPP_A(storage_map::abb(shape_value(), r, c));)
1971 TDLS_ARROW_LUPP_FOR(
1972 m, 0, nc, if constexpr (canonical) rb[m * B + r] = r == col0 + m ? T(1) : T(0);
1973 else rb[m * B + r] = TDLS_ARROW_LUPP_YW(m, r);))
1974 if (refused_count > 0 &&
1975 !second_pass<nc, canonical, internal_rhs, internal_work, internal_matrix,
1976 A_static_stride, work_static_stride, rhs_static_stride>(
1977 A, A_stride, work, work_stride, y, rhs_stride, xcol_stride, refused, Sbb, rb,
1978 BlockSolver::singular_floor))
1979 return false;
1980 int pb[B];
1981 if constexpr (nc == 1) {
1982 if (!BorderSolver::template solve_inplace<true, true, true>(Sbb, 1, pb, 1, rb, 1))
1983 return false;
1984 } else {
1985 if (!BorderSolver::template solve_inplace_multirhs<nc, true, true, true>(Sbb, 1, pb, 1,
1986 rb, 1, B))
1987 return false;
1988 }
1989 if constexpr (canonical) {
1990 TDLS_ARROW_LUPP_FOR(
1991 m, 0, nc,
1992 TDLS_ARROW_LUPP_FOR(r, 0, B,
1993 y[internal_rhs ? unsigned(m * B + r)
1994 : TDLS_STRIDED(r, rhs_stride) +
1995 TDLS_STRIDED_WIDE(m, xcol_stride)] =
1996 rb[m * B + r];))
1997 } else {
1998 TDLS_ARROW_LUPP_FOR(
1999 m, 0, nc, TDLS_ARROW_LUPP_FOR(r, 0, B, TDLS_ARROW_LUPP_YW(m, r) = rb[m * B + r];))
2000 recover<nc, internal_rhs, internal_work, internal_matrix>(
2001 A, A_stride, work, work_stride, y, rhs_stride, xcol_stride, refused, rb);
2002 }
2003 return true;
2004 }
2005
2032 template<int ny, bool internal_rhs, bool internal_matrix>
2033 TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
2034 condense(T* TDLS_RESTRICT A, const int A_stride, T* TDLS_RESTRICT y, const int rhs_stride,
2035 [[maybe_unused]] const std::size_t xcol_stride, const int k, T* TDLS_RESTRICT Sacc,
2036 T* TDLS_RESTRICT racc) noexcept {
2037 constexpr bool forward_check = forward_check_for<ny>;
2038 [[maybe_unused]] constexpr int ycols = ny;
2039 // The block is copied in the order of its storage into a local
2040 // array, and TiledLUpp reads it in that layout. Measured on the RTX
2041 // 5080, a copy in the workspace instead (global memory, SoA) costs
2042 // 3 to 56 % more time although it shrinks the stack: the repeated
2043 // accesses of the factorization stay in the L1 cache here.
2044 T Dl[S * S];
2045 TDLS_ARROW_LUPP_FOR(i, 0, S,
2046 TDLS_ARROW_LUPP_FOR(j, 0, S,
2047 Dl[col ? j * S + i : i * S + j] = TDLS_ARROW_LUPP_A(
2048 storage_map::d(shape_value(), k, i, j));))
2049 int pl[S];
2050 // Factored by TiledLUpp, or in registers: then Dl holds U in
2051 // logical row order and the check of L is made.
2052 bool ok = true;
2053 if constexpr (reg) {
2054 const int verdict = register_factor(Dl, pl);
2055 if (verdict == 2) return false;
2056 ok = verdict == 0;
2057 } else {
2058 if (!BlockSolver::template factorize<true, true>(Dl, 1, pl, 1)) return false;
2059 }
2060 // The couplings of the block under the diagonal coupling; a
2061 // placeholder, never read, otherwise (initialized for the flow
2062 // analysis of the compilers, as the other placeholders below).
2063 T acoef[diag ? S : 1];
2064 if constexpr (diag) {
2065 TDLS_ARROW_LUPP_FOR(
2066 i, 0, S, acoef[i] = TDLS_ARROW_LUPP_A(storage_map::abk_diag(shape_value(), k, i));)
2067 } else {
2068 acoef[0] = T(0);
2069 }
2070 // U in logical row order, unless the forward stage of the solve
2071 // gives the border multipliers, or the block is factored in
2072 // registers (U is then Dl).
2073 T U[forward_check || reg ? 1 : S * S];
2074 if constexpr (forward_check || reg) U[0] = T(0);
2075 if constexpr (!reg) ok = lower_pass<!forward_check>(Dl, pl, U);
2076 const T* const Uf = reg ? Dl : U;
2077 if constexpr (forward_check) {
2078 T W[(B + ny) * S];
2079 TDLS_ARROW_LUPP_FOR(
2080 i, 0, S,
2081 TDLS_ARROW_LUPP_FOR(
2082 c, 0, B,
2083 W[c * S + i] = TDLS_ARROW_LUPP_A(storage_map::akb(shape_value(), k, i, c));)
2084 TDLS_ARROW_LUPP_FOR(m, 0, ny,
2085 W[(B + m) * S + i] = TDLS_ARROW_LUPP_YW(m, B + k * S + i);))
2086 BlockSolver::template forward_inplace_multirhs<B + ny, true, true, true>(Dl, 1, pl, 1,
2087 W, 1, S);
2088 // Column c of U^-T A_kb is row c of the border multipliers.
2089 const T u = condensation_pivot_threshold;
2090 TDLS_ARROW_LUPP_FOR(e, 0, B * S, ok = ok && u * detail::abs(W[e]) <= T(1);)
2091 if (!ok) return false;
2092 BlockSolver::template backward_inplace_multirhs<B + ny, true, true, true>(Dl, 1, pl, 1,
2093 W, 1, S);
2094 condense_pass<0, true, ny, internal_rhs, internal_matrix>(
2095 A, A_stride, y, rhs_stride, xcol_stride, k, Dl, pl, acoef, nullptr, Sacc, racc, W);
2096 } else {
2097 T Mb[multiplier_update ? R * S : 1];
2098 if constexpr (!multiplier_update) Mb[0] = T(0);
2099 ok = border_multipliers<internal_matrix>(A, A_stride, k, Uf, acoef, Mb) && ok;
2100 if (!ok) return false;
2101 if constexpr (multiplier_update && stores_recovery_for<ny>)
2102 store_upper<internal_matrix>(A, A_stride, k, Uf);
2103 condense_pass<0, false, ny, internal_rhs, internal_matrix>(
2104 A, A_stride, y, rhs_stride, xcol_stride, k, Dl, pl, acoef, Mb, Sacc, racc, nullptr);
2105 }
2106 return true;
2107 }
2108
2109// gcc cannot always prove that the provider writes every entry of the
2110// streamed block before it is read or copied, and reports the reads through
2111// -Wmaybe-uninitialized at -O3 (gcc 11; gcc 13 with -march=native; gcc 10
2112// on the reads of the block factored in registers). The provider writes
2113// the block (StreamedBlock), and the copies only keep it as given: the
2114// report is spurious. The suppression is scoped to the three functions
2115// below and to that warning only.
2116#if defined(TDLS_DETAIL_COMPILER_GCC)
2117#pragma GCC diagnostic push
2118#pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
2119#endif
2147 template<int ny, bool internal_rhs, bool internal_matrix, bool recompute = false>
2148 TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
2149 condense_streamed(T* TDLS_RESTRICT A, const int A_stride, T* TDLS_RESTRICT y,
2150 const int rhs_stride, [[maybe_unused]] const std::size_t xcol_stride,
2151 const int k, T* TDLS_RESTRICT Sacc, T* TDLS_RESTRICT racc,
2152 StreamedBlock& blk) noexcept {
2153 constexpr bool forward_check = forward_check_for<ny>;
2154 constexpr int cols = B + ny;
2155 [[maybe_unused]] constexpr int ycols = ny == 0 ? 1 : ny;
2156 const StreamedBlock* const bp = &blk;
2157 T* const Dl = blk.Dl;
2158 T* const W = blk.W;
2159 int pl[S];
2160 bool ok = true;
2161 if constexpr (reg) {
2162 const int verdict = register_factor(Dl, pl);
2163 if (verdict == 2) return false;
2164 ok = verdict == 0;
2165 } else {
2166 if (!BlockSolver::template factorize<true, true>(Dl, 1, pl, 1)) return false;
2167 }
2168 const T* const acoef = blk.Abk;
2169 T U[forward_check || reg ? 1 : S * S];
2170 if constexpr (forward_check || reg) U[0] = T(0);
2171 if constexpr (!reg) ok = lower_pass<!forward_check>(Dl, pl, U);
2172 const T* const Uf = reg ? Dl : U;
2173 T Mb[multiplier_update ? R * S : 1];
2174 if constexpr (!multiplier_update) Mb[0] = T(0);
2175 if constexpr (forward_check) {
2176 BlockSolver::template forward_inplace_multirhs<cols, true, true, true>(Dl, 1, pl, 1, W,
2177 1, S);
2178 const T u = condensation_pivot_threshold;
2179 TDLS_ARROW_LUPP_FOR(e, 0, B * S, ok = ok && u * detail::abs(W[e]) <= T(1);)
2180 if (!ok) return false;
2181 BlockSolver::template backward_inplace_multirhs<cols, true, true, true>(Dl, 1, pl, 1, W,
2182 1, S);
2183 } else {
2184 ok = border_multipliers<internal_matrix>(A, A_stride, k, Uf, acoef, Mb, bp) && ok;
2185 if (!ok) return false;
2186 if constexpr (pass_width_for<ny> < cols) {
2187 // rhs_pass_width columns at a time, read from the block as
2188 // the stored path reads them from the storage (only one
2189 // slice alive at a time); all of them at once otherwise
2190 if constexpr (multiplier_update && stores_recovery_for<ny> && !recompute)
2191 store_upper<internal_matrix>(A, A_stride, k, Uf);
2192 condense_pass<0, false, ny, internal_rhs, internal_matrix, recompute>(
2193 A, A_stride, y, rhs_stride, xcol_stride, k, Dl, pl, acoef, Mb, Sacc, racc,
2194 nullptr, bp);
2195 return true;
2196 }
2197 if constexpr (multiplier_update) {
2198 if constexpr (stores_recovery_for<ny> && !recompute)
2199 store_upper<internal_matrix>(A, A_stride, k, Uf);
2200 if constexpr (reg)
2201 register_solve<cols>(Dl, pl, W);
2202 else
2203 BlockSolver::template forward_inplace_multirhs<cols, true, true, true>(
2204 Dl, 1, pl, 1, W, 1, S);
2205 } else if constexpr (reg) {
2206 register_solve<cols>(Dl, pl, W);
2207 } else if constexpr (Config.gather_rhs) {
2208 T X[cols * S];
2209 BlockSolver::template substitute_multirhs<cols, true, true, true>(Dl, 1, pl, 1, W,
2210 X, 1, S);
2211 TDLS_ARROW_LUPP_FOR(e, 0, cols * S, W[e] = X[e];)
2212 } else {
2213 BlockSolver::template substitute_inplace_multirhs<cols, true, true, true>(
2214 Dl, 1, pl, 1, W, 1, S);
2215 }
2216 }
2217 // the recovery data (not for the canonical solve, nor for the
2218 // recomputed one) and the Schur complement, all the columns (the
2219 // entries are distinct per column: the order of the passes of the
2220 // stored path does not matter)
2221 if constexpr (stores_recovery_for<ny> && !recompute) {
2222 TDLS_ARROW_LUPP_FOR(
2223 i, 0, S,
2224 TDLS_ARROW_LUPP_FOR(j, 0, cols,
2225 if (j < B) TDLS_ARROW_LUPP_A(
2226 storage_map::akb(shape_value(), k, i, j)) = W[j * S + i];
2227 else TDLS_ARROW_LUPP_YW(j - B, B + k * S + i) = W[j * S + i];))
2228 }
2229 if constexpr (multiplier_update) {
2230 TDLS_ARROW_LUPP_FOR(
2231 r, 0, R,
2232 TDLS_ARROW_LUPP_FOR(
2233 j, 0, cols, T sum = T(0);
2234 TDLS_ARROW_LUPP_FOR(l, 0, S, sum += Mb[r * S + l] * W[j * S + l];)
2235 schur_sub<ny, internal_rhs, internal_matrix>(
2236 A, A_stride, y, rhs_stride, xcol_stride, Sacc, racc, r, j, sum);))
2237 } else if constexpr (diag) {
2238 TDLS_ARROW_LUPP_FOR(
2239 i, 0, S,
2240 TDLS_ARROW_LUPP_FOR(j, 0, cols,
2241 schur_sub<ny, internal_rhs, internal_matrix>(
2242 A, A_stride, y, rhs_stride, xcol_stride, Sacc, racc, i, j,
2243 acoef[i] * W[j * S + i]);))
2244 } else {
2245 TDLS_ARROW_LUPP_FOR(
2246 r, 0, B, T a[S];
2247 TDLS_ARROW_LUPP_FOR(i, 0, S, a[i] = blk.Abk[r * S + i];) TDLS_ARROW_LUPP_FOR(
2248 j, 0, cols, T sum = T(0);
2249 TDLS_ARROW_LUPP_FOR(i, 0, S, sum += a[i] * W[j * S + i];)
2250 schur_sub<ny, internal_rhs, internal_matrix>(
2251 A, A_stride, y, rhs_stride, xcol_stride, Sacc, racc, r, j, sum);))
2252 }
2253 return true;
2254 }
2255
2277 template<int ny, bool internal_rhs, bool internal_matrix>
2278 TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
2279 condense_or_store_streamed(T* TDLS_RESTRICT A, const int A_stride, T* TDLS_RESTRICT y,
2280 const int rhs_stride, const std::size_t xcol_stride, const int k,
2281 T* TDLS_RESTRICT Sacc, T* TDLS_RESTRICT racc,
2282 StreamedBlock& blk) noexcept {
2283 constexpr bool forward_check = forward_check_for<ny>;
2284 constexpr int wsave = forward_check ? (B + ny) * S : 1;
2285 T Dsave[S * S];
2286 T Wsave[wsave];
2287 TDLS_ARROW_LUPP_FOR(e, 0, S * S, Dsave[e] = blk.Dl[e];)
2288 if constexpr (forward_check) {
2289 TDLS_ARROW_LUPP_FOR(e, 0, wsave, Wsave[e] = blk.W[e];)
2290 } else {
2291 Wsave[0] = T(0);
2292 }
2293 if (condense_streamed<ny, internal_rhs, internal_matrix>(A, A_stride, y, rhs_stride,
2294 xcol_stride, k, Sacc, racc, blk))
2295 return true;
2296 TDLS_ARROW_LUPP_FOR(e, 0, S * S, blk.Dl[e] = Dsave[e];)
2297 if constexpr (forward_check) { TDLS_ARROW_LUPP_FOR(e, 0, wsave, blk.W[e] = Wsave[e];) }
2298 store_block<ny, internal_rhs, internal_matrix>(A, A_stride, y, rhs_stride, xcol_stride, k,
2299 blk);
2300 return false;
2301 }
2302
2317 template<int ny, bool internal_rhs, bool internal_matrix>
2318 TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
2319 store_block(T* TDLS_RESTRICT A, const int A_stride, T* TDLS_RESTRICT y, const int rhs_stride,
2320 [[maybe_unused]] const std::size_t xcol_stride, const int k,
2321 const StreamedBlock& blk) noexcept {
2322 [[maybe_unused]] constexpr int ycols = ny == 0 ? 1 : ny;
2323 TDLS_ARROW_LUPP_FOR(
2324 i, 0, S,
2325 TDLS_ARROW_LUPP_FOR(j, 0, S,
2326 TDLS_ARROW_LUPP_A(storage_map::d(shape_value(), k, i, j)) =
2327 blk.Dl[col ? j * S + i : i * S + j];)
2328 TDLS_ARROW_LUPP_FOR(c, 0, B,
2329 TDLS_ARROW_LUPP_A(storage_map::akb(shape_value(), k, i, c)) =
2330 blk.W[c * S + i];))
2331 if constexpr (diag) {
2332 TDLS_ARROW_LUPP_FOR(i, 0, S,
2333 TDLS_ARROW_LUPP_A(storage_map::abk_diag(shape_value(), k, i)) =
2334 blk.Abk[i];)
2335 } else {
2336 TDLS_ARROW_LUPP_FOR(
2337 r, 0, B,
2338 TDLS_ARROW_LUPP_FOR(i, 0, S,
2339 TDLS_ARROW_LUPP_A(storage_map::abk(shape_value(), k, r, i)) =
2340 blk.Abk[r * S + i];))
2341 }
2342 if constexpr (ny == 1) {
2343 TDLS_ARROW_LUPP_FOR(i, 0, S, TDLS_ARROW_LUPP_YW(0, B + k * S + i) = blk.W[B * S + i];)
2344 }
2345 }
2346#if defined(TDLS_DETAIL_COMPILER_GCC)
2347#pragma GCC diagnostic pop
2348#endif
2349
2358 template<bool keep_u>
2359 TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
2360 lower_pass(const T* TDLS_RESTRICT Dl, const int* TDLS_RESTRICT pl,
2361 T* TDLS_RESTRICT U) noexcept {
2362 const T u = condensation_pivot_threshold;
2363 bool ok = true;
2364 TDLS_ARROW_LUPP_FOR(l, 0, S, const int p = pl[l]; TDLS_ARROW_LUPP_FOR(
2365 j, 0, S,
2366 if (j < l) {
2367 ok = ok && u * detail::abs(Dl[p * S + j]) <= T(1);
2368 } else if constexpr (keep_u) { U[l * S + j] = Dl[p * S + j]; }))
2369 return ok;
2370 }
2371
2388 template<bool internal_matrix, typename... Streamed>
2389 TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
2390 border_multipliers(const T* TDLS_RESTRICT A, const int A_stride, const int k,
2391 const T* TDLS_RESTRICT U, const T* TDLS_RESTRICT acoef, T* TDLS_RESTRICT Mb,
2392 const Streamed* TDLS_RESTRICT... blk) noexcept {
2393 constexpr bool streamed = sizeof...(Streamed) > 0;
2394 const T u = condensation_pivot_threshold;
2395 bool ok = true;
2396 if constexpr (col) {
2397 TDLS_ARROW_LUPP_FOR(
2398 c, 0, B, T m[S];
2399 TDLS_ARROW_LUPP_FOR(
2400 i, 0, S, if constexpr (streamed) m[i] = streamed_arg(blk...)->W[c * S + i];
2401 else m[i] = TDLS_ARROW_LUPP_A(storage_map::akb(shape_value(), k, i, c));) ok =
2402 multipliers_pass(m, 0, U, u) && ok;)
2403 } else if constexpr (diag) {
2404 // Border row coupling_row + i holds one entry, in column i.
2405 TDLS_ARROW_LUPP_FOR(i, 0, S, T m[S];
2406 TDLS_ARROW_LUPP_FOR(j, 0, S, m[j] = j == i ? acoef[i] : T(0);) ok =
2407 multipliers_pass(m, i, U, u) && ok;)
2408 } else if constexpr (multiplier_update) {
2409 TDLS_ARROW_LUPP_FOR(r, 0, B,
2410 TDLS_ARROW_LUPP_FOR(j, 0, S,
2411 if constexpr (streamed) Mb[r * S + j] =
2412 streamed_arg(blk...)->Abk[r * S + j];
2413 else Mb[r * S + j] = TDLS_ARROW_LUPP_A(
2414 storage_map::abk(shape_value(), k, r, j));)
2415 ok = multipliers_pass(Mb + r * S, 0, U, u) && ok;)
2416 } else {
2417 // MinTraffic: each row checked, then dropped.
2418 TDLS_ARROW_LUPP_FOR(
2419 r, 0, B, T m[S];
2420 TDLS_ARROW_LUPP_FOR(
2421 j, 0, S, if constexpr (streamed) m[j] = streamed_arg(blk...)->Abk[r * S + j];
2422 else m[j] = TDLS_ARROW_LUPP_A(storage_map::abk(shape_value(), k, r, j));) ok =
2423 multipliers_pass(m, 0, U, u) && ok;)
2424 }
2425 return ok;
2426 }
2427
2428// The block factored in registers is, on the streamed and recomputed
2429// paths, the block the provider wrote (see the suppression before
2430// condense_streamed): gcc 10 reports its reads through
2431// -Wmaybe-uninitialized, spuriously. Scoped to the functions below and to
2432// that warning only.
2433#if defined(TDLS_DETAIL_COMPILER_GCC)
2434#pragma GCC diagnostic push
2435#pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
2436#endif
2445 [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr T
2446 magnitude(const T x) noexcept {
2447#if defined(__has_builtin)
2448#if __has_builtin(__builtin_is_constant_evaluated)
2449 if (__builtin_is_constant_evaluated()) return x < T(0) ? -x : x;
2450#endif
2451#elif defined(TDLS_DETAIL_COMPILER_MSVC) && _MSC_VER >= 1925
2452 if (__builtin_is_constant_evaluated()) return x < T(0) ? -x : x;
2453#endif
2454 return std::fabs(x);
2455 }
2456
2464 TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr bool
2465 multipliers_pass(T* TDLS_RESTRICT m, const int first, const T* TDLS_RESTRICT U,
2466 const T u) noexcept {
2467 bool ok = true;
2468 if constexpr (reg) {
2469 // The same operations over full ranges, the triangle selected
2470 // inside, so that U stays in registers (see register_factor).
2471 TDLS_ARROW_LUPP_FOR(
2472 l, 0, S, if (l >= first) {
2473 m[l] *= U[l * S + l];
2474 ok = ok && u * magnitude(m[l]) <= T(1);
2475 TDLS_ARROW_LUPP_FOR(j, 0, S, if (j > l) m[j] -= m[l] * U[l * S + j];)
2476 })
2477 } else {
2478 TDLS_ARROW_LUPP_FOR(l, first, S, m[l] *= U[l * S + l];
2479 ok = ok && u * detail::abs(m[l]) <= T(1);
2480 TDLS_ARROW_LUPP_FOR(j, l + 1, S, m[j] -= m[l] * U[l * S + j];))
2481 }
2482 return ok;
2483 }
2484
2493 template<bool internal_matrix>
2494 TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
2495 store_upper(T* TDLS_RESTRICT A, const int A_stride, const int k,
2496 const T* TDLS_RESTRICT U) noexcept {
2497 if constexpr (reg) {
2498 TDLS_ARROW_LUPP_FOR(l, 0, S,
2499 TDLS_ARROW_LUPP_FOR(j, 0, S,
2500 if (j >= l) TDLS_ARROW_LUPP_A(storage_map::d(
2501 shape_value(), k, l, j)) = U[l * S + j];))
2502 } else {
2503 TDLS_ARROW_LUPP_FOR(l, 0, S,
2504 TDLS_ARROW_LUPP_FOR(j, l, S,
2505 TDLS_ARROW_LUPP_A(storage_map::d(
2506 shape_value(), k, l, j)) = U[l * S + j];))
2507 }
2508 }
2509
2527 TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr int
2528 register_factor(T* TDLS_RESTRICT F, int* TDLS_RESTRICT p) noexcept {
2529 return PhysicalBlockSolver::template physical_factorize_checked<true, true>(
2530 F, 1, p, 1, condensation_pivot_threshold);
2531 }
2532
2547 template<int w, bool backward = !multiplier_update>
2548 TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
2549 register_solve(const T* TDLS_RESTRICT F, const int* TDLS_RESTRICT p,
2550 T* TDLS_RESTRICT W) noexcept {
2551 PhysicalBlockSolver::template forward_inplace_multirhs<w, true, true, true>(F, 1, p, 1, W,
2552 1, S);
2553 if constexpr (col || backward)
2554 PhysicalBlockSolver::template backward_inplace_multirhs<w, true, true, true>(F, 1, p, 1,
2555 W, 1, S);
2556 }
2557
2558#if defined(TDLS_DETAIL_COMPILER_GCC)
2559#pragma GCC diagnostic pop
2560#endif
2561
2577 template<int ny, bool internal_rhs, bool internal_matrix>
2578 TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
2579 schur_sub(T* TDLS_RESTRICT A, const int A_stride, T* TDLS_RESTRICT y, const int rhs_stride,
2580 [[maybe_unused]] const std::size_t xcol_stride, T* TDLS_RESTRICT Sacc,
2581 T* TDLS_RESTRICT racc, const int r, const int c, const T v) noexcept {
2582 [[maybe_unused]] constexpr int ycols = ny;
2583 if constexpr (in_place) {
2584 if (c < B)
2585 TDLS_ARROW_LUPP_A(storage_map::abb(shape_value(), R0 + r, c)) -= v;
2586 else
2587 TDLS_ARROW_LUPP_YW(c - B, R0 + r) -= v;
2588 } else {
2589 if (c < B)
2590 Sacc[r * B + c] -= v;
2591 else
2592 racc[(ny == 1 ? 0 : c - B) * R + r] -= v;
2593 }
2594 }
2595
2626 template<int c0, bool solved, int ny, bool internal_rhs, bool internal_matrix,
2627 bool recompute = false, typename... Streamed>
2628 TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
2629 condense_pass(T* TDLS_RESTRICT A, const int A_stride, T* TDLS_RESTRICT y, const int rhs_stride,
2630 [[maybe_unused]] const std::size_t xcol_stride, const int k,
2631 const T* TDLS_RESTRICT Dl, const int* TDLS_RESTRICT pl,
2632 const T* TDLS_RESTRICT acoef, const T* TDLS_RESTRICT Mb, T* TDLS_RESTRICT Sacc,
2633 T* TDLS_RESTRICT racc, T* TDLS_RESTRICT Wsolved,
2634 const Streamed* TDLS_RESTRICT... blk) noexcept {
2635 constexpr bool streamed = sizeof...(Streamed) > 0;
2636 constexpr int cols = B + ny;
2637 constexpr int w = cols - c0 < pass_width_for<ny> ? cols - c0 : pass_width_for<ny>;
2638 [[maybe_unused]] constexpr int ycols = ny;
2639 T Wl[solved ? 1 : w * S];
2640 T* W = solved ? Wsolved : Wl;
2641 if constexpr (!solved) {
2642 if constexpr (streamed) {
2643 // the slice of the streamed block: [A_kb | y_k] by columns
2644 TDLS_ARROW_LUPP_FOR(
2645 i, 0, S,
2646 TDLS_ARROW_LUPP_FOR(j, 0, w,
2647 W[j * S + i] = streamed_arg(blk...)->W[(c0 + j) * S + i];))
2648 } else {
2649 TDLS_ARROW_LUPP_FOR(
2650 i, 0, S,
2651 TDLS_ARROW_LUPP_FOR(
2652 j, 0, w,
2653 W[j * S + i] = c0 + j < B ? TDLS_ARROW_LUPP_A(storage_map::akb(
2654 shape_value(), k, i, c0 + j))
2655 : TDLS_ARROW_LUPP_YW(c0 + j - B, B + k * S + i);))
2656 }
2657 if constexpr (reg) {
2658 register_solve<w>(Dl, pl, W);
2659 } else if constexpr (multiplier_update) {
2660 BlockSolver::template forward_inplace_multirhs<w, true, true, true>(Dl, 1, pl, 1, W,
2661 1, S);
2662 } else if constexpr (Config.gather_rhs) {
2663 T X[w * S];
2664 BlockSolver::template substitute_multirhs<w, true, true, true>(Dl, 1, pl, 1, W, X,
2665 1, S);
2666 TDLS_ARROW_LUPP_FOR(e, 0, w * S, W[e] = X[e];)
2667 } else {
2668 BlockSolver::template substitute_inplace_multirhs<w, true, true, true>(Dl, 1, pl, 1,
2669 W, 1, S);
2670 }
2671 }
2672 if constexpr (stores_recovery_for<ny> && !recompute) {
2673 TDLS_ARROW_LUPP_FOR(
2674 i, 0, S,
2675 TDLS_ARROW_LUPP_FOR(
2676 j, 0, w,
2677 if (c0 + j < B) TDLS_ARROW_LUPP_A(
2678 storage_map::akb(shape_value(), k, i, c0 + j)) = W[j * S + i];
2679 else TDLS_ARROW_LUPP_YW(c0 + j - B, B + k * S + i) = W[j * S + i];))
2680 }
2681 if constexpr (multiplier_update) {
2682 TDLS_ARROW_LUPP_FOR(
2683 r, 0, R,
2684 TDLS_ARROW_LUPP_FOR(
2685 j, 0, w, T sum = T(0);
2686 TDLS_ARROW_LUPP_FOR(l, 0, S, sum += Mb[r * S + l] * W[j * S + l];)
2687 schur_sub<ny, internal_rhs, internal_matrix>(
2688 A, A_stride, y, rhs_stride, xcol_stride, Sacc, racc, r, c0 + j, sum);))
2689 } else if constexpr (diag) {
2690 TDLS_ARROW_LUPP_FOR(
2691 i, 0, S,
2692 TDLS_ARROW_LUPP_FOR(j, 0, w,
2693 schur_sub<ny, internal_rhs, internal_matrix>(
2694 A, A_stride, y, rhs_stride, xcol_stride, Sacc, racc, i,
2695 c0 + j, acoef[i] * W[j * S + i]);))
2696 } else {
2697 TDLS_ARROW_LUPP_FOR(
2698 r, 0, B, T a[S]; TDLS_ARROW_LUPP_FOR(
2699 i, 0, S, if constexpr (streamed) a[i] = streamed_arg(blk...)->Abk[r * S + i];
2700 else a[i] = TDLS_ARROW_LUPP_A(storage_map::abk(shape_value(), k, r, i));)
2701 TDLS_ARROW_LUPP_FOR(j, 0, w, T sum = T(0);
2702 TDLS_ARROW_LUPP_FOR(i, 0, S, sum += a[i] * W[j * S + i];)
2703 schur_sub<ny, internal_rhs, internal_matrix>(
2704 A, A_stride, y, rhs_stride, xcol_stride, Sacc, racc,
2705 r, c0 + j, sum);))
2706 }
2707 if constexpr (c0 + w < cols)
2708 condense_pass<c0 + w, false, ny, internal_rhs, internal_matrix, recompute>(
2709 A, A_stride, y, rhs_stride, xcol_stride, k, Dl, pl, acoef, Mb, Sacc, racc, nullptr,
2710 blk...);
2711 }
2712
2733 template<int nc, bool internal_rhs, bool internal_work, bool internal_matrix>
2734 TDLS_HOST_DEVICE TDLS_FORCE_INLINE_GPU static constexpr void
2735 recover(const T* TDLS_RESTRICT A, const int A_stride, const T* TDLS_RESTRICT work,
2736 const int work_stride, T* TDLS_RESTRICT y, const int rhs_stride,
2737 [[maybe_unused]] const std::size_t xcol_stride, const unsigned* TDLS_RESTRICT refused,
2738 const T* TDLS_RESTRICT rb) noexcept {
2739 constexpr int fcols = front_cols_for<nc>;
2740 [[maybe_unused]] constexpr int ycols = nc;
2741 // w and s: column m at m * R and m * S.
2742 T w[R * nc];
2743 TDLS_ARROW_LUPP_FOR(e, 0, R * nc, w[e] = T(0);)
2744 for (int k = NB - 1; k >= 0; --k) {
2745 T s[S * nc];
2746 TDLS_ARROW_LUPP_FOR(
2747 i, 0, S,
2748 TDLS_ARROW_LUPP_FOR(m, 0, nc, s[m * S + i] = TDLS_ARROW_LUPP_YW(m, B + k * S + i);)
2749 TDLS_ARROW_LUPP_FOR(
2750 c, 0, B,
2751 const T a = TDLS_ARROW_LUPP_A(storage_map::akb(shape_value(), k, i, c));
2752 TDLS_ARROW_LUPP_FOR(m, 0, nc, s[m * S + i] -= a * rb[m * B + c];)))
2753 if ((refused[k >> 5] >> (k & 31)) & 1u) {
2754 // x_k = c_k - G_k xb - H_k w, then w += A_bk x_k.
2755 TDLS_ARROW_LUPP_FOR(
2756 i, 0, S,
2757 TDLS_ARROW_LUPP_FOR(
2758 a, 0, R,
2759 const T h = a < S
2760 ? TDLS_ARROW_LUPP_A(storage_map::d(shape_value(), k, i, a))
2761 : TDLS_ARROW_LUPP_W(front_rows * fcols +
2762 (k * S + i) * extra_cols + a - S);
2763 TDLS_ARROW_LUPP_FOR(m, 0, nc, s[m * S + i] -= h * w[m * R + a];)))
2764 if constexpr (diag) {
2765 TDLS_ARROW_LUPP_FOR(
2766 i, 0, S,
2767 const T coef =
2768 TDLS_ARROW_LUPP_A(storage_map::abk_diag(shape_value(), k, i));
2769 TDLS_ARROW_LUPP_FOR(m, 0, nc, w[m * R + i] += coef * s[m * S + i];))
2770 } else {
2771 TDLS_ARROW_LUPP_FOR(
2772 r, 0, B,
2773 TDLS_ARROW_LUPP_FOR(
2774 i, 0, S,
2775 const T coef =
2776 TDLS_ARROW_LUPP_A(storage_map::abk(shape_value(), k, r, i));
2777 TDLS_ARROW_LUPP_FOR(m, 0, nc, w[m * R + r] += coef * s[m * S + i];)))
2778 }
2779 } else if constexpr (multiplier_update) {
2780 // x_k = U^-1 s, U in the upper triangle of D_k, 1 / pivot on
2781 // its diagonal. Columns j = l + 1 + u, u < t: gcc -O2 takes
2782 // the empty range j in [l + 1, S) of S = 1 for a nonempty
2783 // one and warns of an out-of-bounds s.
2784 TDLS_ARROW_LUPP_FOR(
2785 t, 0, S, const int l = S - 1 - t;
2786 TDLS_ARROW_LUPP_FOR(
2787 u, 0, t, const int j = l + 1 + u;
2788 const T ulj = TDLS_ARROW_LUPP_A(storage_map::d(shape_value(), k, l, j));
2789 TDLS_ARROW_LUPP_FOR(m, 0, nc, s[m * S + l] -= ulj * s[m * S + j];))
2790 const T ull = TDLS_ARROW_LUPP_A(storage_map::d(shape_value(), k, l, l));
2791 TDLS_ARROW_LUPP_FOR(m, 0, nc, s[m * S + l] *= ull;))
2792 }
2793 TDLS_ARROW_LUPP_FOR(
2794 i, 0, S,
2795 TDLS_ARROW_LUPP_FOR(m, 0, nc, TDLS_ARROW_LUPP_YW(m, B + k * S + i) = s[m * S + i];))
2796 }
2797 }
2798
2839 template<int nc, bool canonical, bool internal_rhs, bool internal_work, bool internal_matrix,
2840 int A_static_stride = dynamic_stride, int work_static_stride = dynamic_stride,
2841 int rhs_static_stride = dynamic_stride>
2842 TDLS_HOST_DEVICE TDLS_NO_INLINE static constexpr bool
2843 second_pass(T* TDLS_RESTRICT A, const int A_stride_arg, T* TDLS_RESTRICT work,
2844 const int work_stride_arg, T* TDLS_RESTRICT y, const int rhs_stride_arg,
2845 [[maybe_unused]] const std::size_t xcol_stride,
2846 const unsigned* TDLS_RESTRICT refused, T* TDLS_RESTRICT Sbb, T* TDLS_RESTRICT rb,
2847 const T singular_floor) noexcept {
2848 const int A_stride = detail::stride_or<A_static_stride>(A_stride_arg);
2849 const int work_stride = detail::stride_or<work_static_stride>(work_stride_arg);
2850 const int rhs_stride = detail::stride_or<rhs_static_stride>(rhs_stride_arg);
2851 constexpr int fcols = front_cols_for<nc>;
2852 [[maybe_unused]] constexpr int ycols = nc;
2853 // First column of the right-hand sides in the front.
2854 constexpr int rhs0 = S + B + R;
2855 const T u = condensation_pivot_threshold;
2856 // The border rows [T A_bk | Sbb | T | rb], T starting as the
2857 // injection of the R coupled rows.
2858 for (int r = 0; r < B; ++r) {
2859 for (int c = 0; c < B; ++c)
2860 TDLS_ARROW_LUPP_F(S + r, S + c) = Sbb[r * B + c];
2861 for (int a = 0; a < R; ++a)
2862 TDLS_ARROW_LUPP_F(S + r, S + B + a) = r == R0 + a ? T(1) : T(0);
2863 for (int m = 0; m < nc; ++m)
2864 TDLS_ARROW_LUPP_F(S + r, rhs0 + m) = rb[m * B + r];
2865 }
2866 for (int k = 0; k < NB; ++k) {
2867 if (!((refused[k >> 5] >> (k & 31)) & 1u)) continue;
2868 // The rows of the block: [D_k | A_kb | 0 | Y_k].
2869 for (int i = 0; i < S; ++i) {
2870 for (int j = 0; j < S; ++j)
2871 TDLS_ARROW_LUPP_F(i, j) =
2872 TDLS_ARROW_LUPP_A(storage_map::d(shape_value(), k, i, j));
2873 for (int c = 0; c < B; ++c)
2874 TDLS_ARROW_LUPP_F(i, S + c) =
2875 TDLS_ARROW_LUPP_A(storage_map::akb(shape_value(), k, i, c));
2876 for (int a = 0; a < R; ++a)
2877 TDLS_ARROW_LUPP_F(i, S + B + a) = T(0);
2878 for (int m = 0; m < nc; ++m) {
2879 if constexpr (canonical)
2880 TDLS_ARROW_LUPP_F(i, rhs0 + m) = T(0);
2881 else
2882 TDLS_ARROW_LUPP_F(i, rhs0 + m) = TDLS_ARROW_LUPP_YW(m, B + k * S + i);
2883 }
2884 }
2885 // The border rows meet the columns of the block through T A_bk.
2886 for (int r = 0; r < B; ++r) {
2887 for (int j = 0; j < S; ++j) {
2888 T v = T(0);
2889 if constexpr (diag) {
2890 v = TDLS_ARROW_LUPP_F(S + r, S + B + j) *
2891 TDLS_ARROW_LUPP_A(storage_map::abk_diag(shape_value(), k, j));
2892 } else {
2893 for (int a = 0; a < B; ++a)
2894 v += TDLS_ARROW_LUPP_F(S + r, S + B + a) *
2895 TDLS_ARROW_LUPP_A(storage_map::abk(shape_value(), k, a, j));
2896 }
2897 TDLS_ARROW_LUPP_F(S + r, j) = v;
2898 }
2899 }
2900 // S pivots over the S + B rows: a row of the block when it
2901 // reaches u times the largest candidate, the largest otherwise.
2902 bool block_row[front_rows];
2903 for (int e = 0; e < front_rows; ++e)
2904 block_row[e] = e < S;
2905 for (int j = 0; j < S; ++j) {
2906 int pmax = j, pblk = -1;
2907 T vmax = T(0), vblk = T(0);
2908 for (int e = j; e < front_rows; ++e) {
2909 const T v = detail::abs(TDLS_ARROW_LUPP_F(e, j));
2910 if (v > vmax) {
2911 vmax = v;
2912 pmax = e;
2913 }
2914 if (block_row[e] && v > vblk) {
2915 vblk = v;
2916 pblk = e;
2917 }
2918 }
2919 if (!(vmax >= singular_floor)) return false;
2920 const int p = pblk >= 0 && vblk >= u * vmax && vblk >= singular_floor ? pblk : pmax;
2921 if (p != j) {
2922 for (int c = j; c < fcols; ++c) {
2923 const T t = TDLS_ARROW_LUPP_F(p, c);
2924 TDLS_ARROW_LUPP_F(p, c) = TDLS_ARROW_LUPP_F(j, c);
2925 TDLS_ARROW_LUPP_F(j, c) = t;
2926 }
2927 const bool t = block_row[p];
2928 block_row[p] = block_row[j];
2929 block_row[j] = t;
2930 }
2931 const T inv = T(1) / TDLS_ARROW_LUPP_F(j, j);
2932 TDLS_ARROW_LUPP_F(j, j) = inv;
2933 for (int e = j + 1; e < front_rows; ++e) {
2934 const T l = TDLS_ARROW_LUPP_F(e, j) * inv;
2935 if (l == T(0)) continue;
2936 for (int c = j + 1; c < fcols; ++c)
2937 TDLS_ARROW_LUPP_F(e, c) -= l * TDLS_ARROW_LUPP_F(j, c);
2938 }
2939 }
2940 if constexpr (!canonical) {
2941 // [G_k | H_k | c_k] = U^-1 [V | Z | h] on the S pivot rows,
2942 // then stored in place of A_kb, D_k (and after the front)
2943 // and Y_k.
2944 for (int c = S; c < fcols; ++c) {
2945 for (int i = S - 1; i >= 0; --i) {
2946 T t = TDLS_ARROW_LUPP_F(i, c);
2947 for (int m = i + 1; m < S; ++m)
2948 t -= TDLS_ARROW_LUPP_F(i, m) * TDLS_ARROW_LUPP_F(m, c);
2949 TDLS_ARROW_LUPP_F(i, c) = t * TDLS_ARROW_LUPP_F(i, i);
2950 }
2951 }
2952 for (int i = 0; i < S; ++i) {
2953 for (int c = 0; c < B; ++c)
2954 TDLS_ARROW_LUPP_A(storage_map::akb(shape_value(), k, i, c)) =
2955 TDLS_ARROW_LUPP_F(i, S + c);
2956 for (int a = 0; a < R; ++a) {
2957 if (a < S)
2958 TDLS_ARROW_LUPP_A(storage_map::d(shape_value(), k, i, a)) =
2959 TDLS_ARROW_LUPP_F(i, S + B + a);
2960 else
2961 TDLS_ARROW_LUPP_W(front_rows * fcols + (k * S + i) * extra_cols + a -
2962 S) = TDLS_ARROW_LUPP_F(i, S + B + a);
2963 }
2964 for (int m = 0; m < nc; ++m)
2965 TDLS_ARROW_LUPP_YW(m, B + k * S + i) = TDLS_ARROW_LUPP_F(i, rhs0 + m);
2966 }
2967 }
2968 }
2969 for (int r = 0; r < B; ++r) {
2970 for (int c = 0; c < B; ++c)
2971 Sbb[r * B + c] = TDLS_ARROW_LUPP_F(S + r, S + c);
2972 for (int m = 0; m < nc; ++m)
2973 rb[m * B + r] = TDLS_ARROW_LUPP_F(S + r, rhs0 + m);
2974 }
2975 return true;
2976 }
2977};
2978
2979
2980
2981#undef TDLS_ARROW_LUPP_A
2982#undef TDLS_ARROW_LUPP_Y
2983#undef TDLS_ARROW_LUPP_YW
2984#undef TDLS_ARROW_LUPP_W
2985#undef TDLS_ARROW_LUPP_F
2986#undef TDLS_ARROW_LUPP_FOR
2987#undef TDLS_ARROW_LUPP_SIDED
2988
2989
2990
2991} // namespace tdls
2992
2993#if defined(TDLS_DETAIL_COMPILER_NVCC)
2994#pragma nv_diagnostic pop
2995#endif
2996#if defined(TDLS_DETAIL_COMPILER_CLANG)
2997#pragma clang diagnostic pop
2998#endif
2999
3000
3001
3002#endif // TDLS_SOLVERS_ARROWHEAD_LUPP_SOLVER_STATIC_HPP
config.hpp
Shape and compile-time configuration of the ArrowheadLUpp solver family.
tdls::detail::ArrowheadRecomputedOutcome
ArrowheadRecomputedOutcome
Outcome of ArrowheadLUppSolverStatic::solve_inplace_recomputed and solve_inplace_recomputed_or_tangen...
Definition solver_static.hpp:295
tdls::detail::ArrowheadRecomputedOutcome::deferred
@ deferred
TDLS_ARROW_LUPP_A
#define TDLS_ARROW_LUPP_A(idx)
Stored entry at index idx of the matrix: contiguous under internal residency, strided otherwise.
Definition solver_static.hpp:152
TDLS_ARROW_LUPP_W
#define TDLS_ARROW_LUPP_W(e)
Entry e of the workspace: contiguous under internal residency, strided otherwise.
Definition solver_static.hpp:170
TDLS_ARROW_LUPP_F
#define TDLS_ARROW_LUPP_F(r, c)
Entry (r, c) of the front of the second pass, row-major at the start of the workspace,...
Definition solver_static.hpp:174
TDLS_ARROW_LUPP_SIDED
#define TDLS_ARROW_LUPP_SIDED(...)
Under the side dispatch of nvc++ (TDLS_TARGET_GPU), when the host and the device choose differently (...
Definition solver_static.hpp:190
TDLS_ARROW_LUPP_YW
#define TDLS_ARROW_LUPP_YW(m, i)
Entry i of column m of the right-hand sides and solutions: columns of dimension contiguous entries un...
Definition solver_static.hpp:163
TDLS_ARROW_LUPP_FOR
#define TDLS_ARROW_LUPP_FOR(i, first, last,...)
for (int i = first; i < last; ++i) { body } under the unroll policy of the configuration: the two-bra...
Definition solver_static.hpp:197
TDLS_ARROW_LUPP_Y
#define TDLS_ARROW_LUPP_Y(i)
Entry i of the right-hand side and solution: contiguous under internal residency, strided otherwise.
Definition solver_static.hpp:156
tdls::detail::ArrowheadStreamedOutcome
ArrowheadStreamedOutcome
Outcome of ArrowheadLUppSolverStatic::solve_inplace_streamed_or_tangent.
Definition solver_static.hpp:286
tdls::detail::ArrowheadStreamedOutcome::singular
@ singular
the matrix is singular
tdls::detail::ArrowheadStreamedOutcome::tangent
@ tangent
tdls::detail::ArrowheadStreamedOutcome::solved
@ solved
y holds the solution (converged() returned false)
macros.hpp
Toolchain detection and portability macros.
TDLS_NO_INLINE_GPU
#define TDLS_NO_INLINE_GPU
Keeps a function out of line on GPU only: noinline under CUDA, HIP and in the device pass of SYCL,...
Definition macros.hpp:125
TDLS_FORCE_INLINE_GPU
#define TDLS_FORCE_INLINE_GPU
Inlining forced on GPU only: __forceinline__ under CUDA, the equivalent attribute under HIP (see TDLS...
Definition macros.hpp:77
TDLS_FORCE_INLINE_ALL
#define TDLS_FORCE_INLINE_ALL
Inlining forced on every backend, CPU included: TDLS_FORCE_INLINE_GPU on GPU, __forceinline under MSV...
Definition macros.hpp:101
TDLS_RESTRICT
#define TDLS_RESTRICT
Non-aliasing pointer qualifier (__restrict__; __restrict on MSVC; nothing in the code that nvc++ gene...
Definition macros.hpp:174
TDLS_STRIDED_WIDE
#define TDLS_STRIDED_WIDE(i, s)
Offset of entry i of a strided object, i times stride s, for a 64-bit stride s (std::size_t); i is a ...
Definition macros.hpp:227
TDLS_HOST_DEVICE
#define TDLS_HOST_DEVICE
__host__ __device__ under CUDA, the equivalent attributes under HIP, empty elsewhere.
Definition macros.hpp:50
TDLS_STRIDED
#define TDLS_STRIDED(i, s)
Offset of entry i of a strided object, i times stride s, in 64-bit arithmetic (std::size_t); i and s ...
Definition macros.hpp:206
TDLS_NO_INLINE
#define TDLS_NO_INLINE
Keeps a function out of line on every backend.
Definition macros.hpp:146
math.hpp
Scalar math helpers usable in constant expressions.
tdls::dynamic_stride
constexpr int dynamic_stride
The compile-time stride of an operand whose stride is known at run time only: the default of the stri...
Definition options.hpp:181
storage_map.hpp
Position of every entry of an arrowhead matrix in its storage.
tdls::ArrowheadLUppConfig
Compile-time knobs of the ArrowheadLUpp solvers, carrying the tuned defaults.
Definition config.hpp:136
tdls::ArrowheadLUppConfig::schur_in_place
bool schur_in_place
Definition config.hpp:293
tdls::ArrowheadLUppConfig::streamed_out_of_line
bool streamed_out_of_line
Definition config.hpp:311
tdls::ArrowheadLUppConfig::storage
ArrowheadStorage storage
Definition config.hpp:147
tdls::ArrowheadLUppConfig::rhs_pass_width
int rhs_pass_width
Definition config.hpp:274
tdls::ArrowheadLUppConfig::block_row_interchange
RowInterchange block_row_interchange
Definition config.hpp:231
tdls::ArrowheadLUppSolverStatic
Arrowhead solver with a compile-time shape.
Definition solver_static.hpp:315
tdls::ArrowheadLUppSolverStatic::solve
static TDLS_HOST_DEVICE constexpr bool solve(T *TDLS_RESTRICT A, const int A_stride, T *TDLS_RESTRICT work, const int work_stride, const T *TDLS_RESTRICT b, T *TDLS_RESTRICT x, const int rhs_stride) noexcept
Solves A x = b, b left untouched.
Definition solver_static.hpp:829
tdls::ArrowheadLUppSolverStatic::solve_border_canonical_multirhs
static TDLS_HOST_DEVICE constexpr bool solve_border_canonical_multirhs(T *TDLS_RESTRICT A, const int A_stride, T *TDLS_RESTRICT work, const int work_stride, const int col0, T *TDLS_RESTRICT x, const int rhs_stride, const std::size_t xcol_stride) noexcept
Columns col0 .. col0 + nrhs - 1 of the border block of A^-1 (col0 + nrhs <= border).
Definition solver_static.hpp:1086
tdls::ArrowheadLUppSolverStatic::solve_inplace
static TDLS_HOST_DEVICE constexpr bool solve_inplace(T *TDLS_RESTRICT A, const int A_stride, T *TDLS_RESTRICT work, const int work_stride, T *TDLS_RESTRICT y, const int rhs_stride) noexcept
Solves A x = y in place.
Definition solver_static.hpp:756
tdls::ArrowheadLUppSolverStatic::ArrowheadLUppSolverStatic
friend struct ArrowheadLUppSolverStatic
Definition solver_static.hpp:358
tdls::ArrowheadLUppSolverStatic::dimension
static constexpr int dimension
Dimension of the matrix, border + block * blocks.
Definition solver_static.hpp:321
tdls::ArrowheadLUppSolverStatic::storage_map
ArrowheadStorageMap< Config.storage, Config.coupling, Config.layout > storage_map
Index of every entry of the matrix in its storage.
Definition solver_static.hpp:317
tdls::ArrowheadLUppSolverStatic::shape
static constexpr ArrowheadShape shape
Shape of the matrix.
Definition solver_static.hpp:319
tdls::ArrowheadLUppSolverStatic::solve_multirhs
static TDLS_HOST_DEVICE constexpr bool solve_multirhs(T *TDLS_RESTRICT A, const int A_stride, T *TDLS_RESTRICT work, const int work_stride, const T *TDLS_RESTRICT b, T *TDLS_RESTRICT x, const int rhs_stride, const std::size_t xcol_stride, int &refused_count) noexcept
Solves A X = B for nrhs right-hand sides, B left untouched, counting the refused blocks: B is copied ...
Definition solver_static.hpp:948
tdls::ArrowheadLUppSolverStatic::solve_inplace
static TDLS_HOST_DEVICE constexpr bool solve_inplace(T *TDLS_RESTRICT A, const int A_stride, T *TDLS_RESTRICT work, const int work_stride, T *TDLS_RESTRICT y, const int rhs_stride, int &refused_count) noexcept
Solves A x = y in place, counting the refused blocks.
Definition solver_static.hpp:727
tdls::ArrowheadLUppSolverStatic::workspace_size
static constexpr int workspace_size
Definition solver_static.hpp:395
tdls::ArrowheadLUppSolverStatic::solve_border_canonical_multirhs_streamed
static TDLS_HOST_DEVICE constexpr bool solve_border_canonical_multirhs_streamed(Provider &&provider, T *TDLS_RESTRICT A, const int A_stride, T *TDLS_RESTRICT work, const int work_stride, const int col0, T *TDLS_RESTRICT x, const int rhs_stride, const std::size_t xcol_stride, int &refused_count) noexcept
solve_border_canonical_multirhs with the blocks streamed (see solve_inplace_streamed; blk....
Definition solver_static.hpp:684
tdls::ArrowheadLUppSolverStatic::StreamedBlock
detail::ArrowheadStreamedBlock< T, S, B, col, diag > StreamedBlock
Block k of a streamed solve (detail::ArrowheadStreamedBlock): the working arrays of the solver,...
Definition solver_static.hpp:512
tdls::ArrowheadLUppSolverStatic::RecomputedBlock
detail::ArrowheadRecomputedBlock< T, S, B, col, diag > RecomputedBlock
Block k of a recomputed solve (detail::ArrowheadRecomputedBlock): a StreamedBlock and the recovery fl...
Definition solver_static.hpp:515
tdls::ArrowheadLUppSolverStatic::RecomputedOutcome
detail::ArrowheadRecomputedOutcome RecomputedOutcome
Outcome of solve_inplace_recomputed and solve_inplace_recomputed_or_tangent (detail::ArrowheadRecompu...
Definition solver_static.hpp:1098
tdls::ArrowheadLUppSolverStatic::solve_inplace_recomputed_or_tangent
static TDLS_HOST_DEVICE constexpr RecomputedOutcome solve_inplace_recomputed_or_tangent(Provider &&provider, Converged &&converged, const T *TDLS_RESTRICT Abb, const int Abb_stride, T *TDLS_RESTRICT work, const int work_stride, T *TDLS_RESTRICT y, const int rhs_stride, const int col0, T *TDLS_RESTRICT x, const int x_stride, const std::size_t xcol_stride, int &refused_count) noexcept
One recomputed pass of a Newton iteration: the blocks condensed as solve_inplace_recomputed condenses...
Definition solver_static.hpp:1235
tdls::ArrowheadLUppSolverStatic::solve_border_canonical_multirhs
static TDLS_HOST_DEVICE constexpr bool solve_border_canonical_multirhs(T *TDLS_RESTRICT A, const int A_stride, T *TDLS_RESTRICT work, const int work_stride, const int col0, T *TDLS_RESTRICT x, const int rhs_stride, const std::size_t xcol_stride, int &refused_count) noexcept
Columns col0 .. col0 + nrhs - 1 of the border block of A^-1, the inverse of the Schur complement of t...
Definition solver_static.hpp:1047
tdls::ArrowheadLUppSolverStatic::solve_multirhs
static TDLS_HOST_DEVICE constexpr bool solve_multirhs(T *TDLS_RESTRICT A, const int A_stride, T *TDLS_RESTRICT work, const int work_stride, const T *TDLS_RESTRICT b, T *TDLS_RESTRICT x, const int rhs_stride, const std::size_t xcol_stride) noexcept
Solves A X = B for nrhs right-hand sides, B left untouched.
Definition solver_static.hpp:992
tdls::ArrowheadLUppSolverStatic::StreamedOutcome
detail::ArrowheadStreamedOutcome StreamedOutcome
Outcome of solve_inplace_streamed_or_tangent (detail::ArrowheadStreamedOutcome).
Definition solver_static.hpp:565
tdls::ArrowheadLUppSolverStatic::pass_width
static constexpr int pass_width
Columns [A_kb | y_k] solved per pass, one right-hand side.
Definition solver_static.hpp:397
tdls::ArrowheadLUppSolverStatic::solve_inplace_streamed_or_tangent
static TDLS_HOST_DEVICE constexpr StreamedOutcome solve_inplace_streamed_or_tangent(Provider &&provider, Converged &&converged, T *TDLS_RESTRICT A, const int A_stride, T *TDLS_RESTRICT work, const int work_stride, T *TDLS_RESTRICT y, const int rhs_stride, const int col0, T *TDLS_RESTRICT x, const int x_stride, const std::size_t xcol_stride, int &refused_count) noexcept
One streamed pass of a Newton iteration: the blocks condensed as solve_inplace_streamed does,...
Definition solver_static.hpp:623
tdls::ArrowheadLUppSolverStatic::solve
static TDLS_HOST_DEVICE constexpr bool solve(T *TDLS_RESTRICT A, const int A_stride, T *TDLS_RESTRICT work, const int work_stride, const T *TDLS_RESTRICT b, T *TDLS_RESTRICT x, const int rhs_stride, int &refused_count) noexcept
Solves A x = b, b left untouched, counting the refused blocks: b is copied into x,...
Definition solver_static.hpp:791
tdls::ArrowheadLUppSolverStatic::recomputed_workspace_size
static constexpr int recomputed_workspace_size
Entries of the workspace of solve_inplace_recomputed with a store of capacity refused blocks,...
Definition solver_static.hpp:1107
tdls::ArrowheadLUppSolverStatic::solve_inplace_recomputed
static TDLS_HOST_DEVICE constexpr RecomputedOutcome solve_inplace_recomputed(Provider &&provider, const T *TDLS_RESTRICT Abb, const int Abb_stride, T *TDLS_RESTRICT work, const int work_stride, T *TDLS_RESTRICT y, const int rhs_stride, int &refused_count) noexcept
Solves A x = y in place, the blocks recomputed: provider(k,blk) fills block k (RecomputedBlock) once ...
Definition solver_static.hpp:1170
tdls::ArrowheadLUppSolverStatic::multirhs_workspace_size
static constexpr int multirhs_workspace_size
Definition solver_static.hpp:391
tdls::ArrowheadLUppSolverStatic::solve_inplace_multirhs
static TDLS_HOST_DEVICE constexpr bool solve_inplace_multirhs(T *TDLS_RESTRICT A, const int A_stride, T *TDLS_RESTRICT work, const int work_stride, T *TDLS_RESTRICT y, const int rhs_stride, const std::size_t xcol_stride, int &refused_count) noexcept
Solves A X = Y in place for nrhs right-hand sides, counting the refused blocks. The right-hand sides ...
Definition solver_static.hpp:873
tdls::ArrowheadLUppSolverStatic::storage_size
static constexpr int storage_size
Number of stored entries of the matrix.
Definition solver_static.hpp:323
tdls::ArrowheadLUppSolverStatic::solve_inplace_multirhs
static TDLS_HOST_DEVICE constexpr bool solve_inplace_multirhs(T *TDLS_RESTRICT A, const int A_stride, T *TDLS_RESTRICT work, const int work_stride, T *TDLS_RESTRICT y, const int rhs_stride, const std::size_t xcol_stride) noexcept
Solves A X = Y in place for nrhs right-hand sides.
Definition solver_static.hpp:908
tdls::ArrowheadLUppSolverStatic::solve_inplace_streamed
static TDLS_HOST_DEVICE constexpr bool solve_inplace_streamed(Provider &&provider, T *TDLS_RESTRICT A, const int A_stride, T *TDLS_RESTRICT work, const int work_stride, T *TDLS_RESTRICT y, const int rhs_stride, int &refused_count) noexcept
Solves A x = y in place, the blocks streamed: provider(k, blk) fills block k (StreamedBlock) when the...
Definition solver_static.hpp:550
tdls::ArrowheadShape
Shape of an arrowhead matrix: a template argument of the compile-time solver, a function argument of ...
Definition config.hpp:111
tdls::ArrowheadShape::n
TDLS_HOST_DEVICE constexpr int n() const noexcept
Dimension of the matrix.
Definition config.hpp:125
tdls::ArrowheadStorageMap
Index of every entry of an arrowhead matrix in its storage.
Definition storage_map.hpp:56
tdls::ArrowheadStorageMap::wide_size
static TDLS_HOST_DEVICE constexpr long long wide_size(const ArrowheadShape &s) noexcept
Number of stored entries, in 64-bit arithmetic (the overflow check of valid()).
Definition storage_map.hpp:87
tdls::ArrowheadStorageMap::abb
TDLS_HOST_DEVICE static TDLS_FORCE_INLINE_GPU constexpr int abb(const ArrowheadShape &s, const int r, const int c) noexcept
Border-border entry (r, c), r and c < border.
Definition storage_map.hpp:121
tdls::ArrowheadStorageMap::size
TDLS_HOST_DEVICE static TDLS_FORCE_INLINE_GPU constexpr int size(const ArrowheadShape &s) noexcept
Number of stored entries.
Definition storage_map.hpp:111
tdls::ArrowheadStorageMap::abk
TDLS_HOST_DEVICE static TDLS_FORCE_INLINE_GPU constexpr int abk(const ArrowheadShape &s, const int k, const int r, const int i) noexcept
Coupling entry of border row r and unknown i of block k (dense coupling).
Definition storage_map.hpp:145
tdls::ArrowheadStorageMap::akb
TDLS_HOST_DEVICE static TDLS_FORCE_INLINE_GPU constexpr int akb(const ArrowheadShape &s, const int k, const int i, const int c) noexcept
Entry of row i of block k in border column c.
Definition storage_map.hpp:189
tdls::ArrowheadStorageMap::d
TDLS_HOST_DEVICE static TDLS_FORCE_INLINE_GPU constexpr int d(const ArrowheadShape &s, const int k, const int i, const int j) noexcept
Entry (i, j) of the diagonal block k.
Definition storage_map.hpp:213
tdls::ArrowheadStorageMap::abk_diag
TDLS_HOST_DEVICE static TDLS_FORCE_INLINE_GPU constexpr int abk_diag(const ArrowheadShape &s, const int k, const int i) noexcept
Coupling entry of border row coupling_row + i and unknown i of block k (diagonal coupling).
Definition storage_map.hpp:167
tdls::TiledLUppSolverStatic
Tiled dense LU factorization with logical partial pivoting and out-of-tile pivot recovery,...
Definition solver_static.hpp:255
tdls::TiledLUppSolverStatic::singular_floor
static constexpr T singular_floor
Singularity floor of the out-of-tile recovery, read once from the configuration (see TiledLUppConfig:...
Definition solver_static.hpp:269
tdls::detail::ArrowheadRecomputedBlock
Block k of a recomputed solve (solve_inplace_recomputed, solve_inplace_recomputed_or_tangent): a Stre...
Definition solver_static.hpp:278
tdls::detail::ArrowheadRecomputedBlock::recovery
bool recovery
Definition solver_static.hpp:282
tdls::detail::ArrowheadStreamedBlock
Block k of a streamed solve of ArrowheadLUppSolverStatic (ArrowheadLUppSolverStatic::StreamedBlock): ...
Definition solver_static.hpp:229
tdls::detail::ArrowheadStreamedBlock::Abk
T Abk[diag ? S :B *S]
Definition solver_static.hpp:232
tdls::detail::ArrowheadStreamedBlock::y
TDLS_HOST_DEVICE TDLS_FORCE_INLINE_ALL constexpr T & y(const int i) noexcept
Definition solver_static.hpp:248
tdls::detail::ArrowheadStreamedBlock::Dl
T Dl[S *S]
D_k in the layout of the factorization.
Definition solver_static.hpp:230
tdls::detail::ArrowheadStreamedBlock::W
T W[(B+1) *S]
[A_kb | y_k], column c at c * S
Definition solver_static.hpp:231
tdls::detail::ArrowheadStreamedBlock::d
TDLS_HOST_DEVICE TDLS_FORCE_INLINE_ALL constexpr T & d(const int i, const int j) noexcept
Definition solver_static.hpp:236
tdls::detail::ArrowheadStreamedBlock::akb
TDLS_HOST_DEVICE TDLS_FORCE_INLINE_ALL constexpr T & akb(const int i, const int c) noexcept
Definition solver_static.hpp:242
tdls::detail::ArrowheadStreamedBlock::coupling
TDLS_HOST_DEVICE TDLS_FORCE_INLINE_ALL constexpr T & coupling(const int i) noexcept
Definition solver_static.hpp:253
tdls::detail::ArrowheadStreamedBlock::abk
TDLS_HOST_DEVICE TDLS_FORCE_INLINE_ALL constexpr T & abk(const int r, const int i) noexcept
Definition solver_static.hpp:259
target.hpp
Detection of the compiler, of the language and of the side of the compilation pass: the one place whe...
solver_static.hpp
TiledLUpp solver: LU with partial pivoting on tile grids, one thread / work-item per system.
Generated by doxygen 1.12.0