@@ -60,11 +60,11 @@ struct sparse_lu_entry_trait<Tensor, RHSVector, XVector, enable_tensor_lu_t<Tens
6060 static constexpr Idx block_size = Tensor::RowsAtCompileTime;
6161 using Scalar = typename Tensor::Scalar;
6262 using Matrix = Eigen::Matrix<Scalar, block_size, block_size, Tensor::Options>;
63- using LUFactor = Eigen::FullPivLU<Eigen::Ref<Matrix>>;
63+ using LUFactor = Eigen::FullPivLU<Eigen::Ref<Matrix>>; // LU decomposition with full pivoting in place
6464 struct BlockPerm {
6565 typename LUFactor::PermutationPType p;
6666 typename LUFactor::PermutationQType q;
67- };
67+ }; // Extract permutation matrices p and q from LUFactor
6868 using BlockPermArray = std::vector<BlockPerm>;
6969};
7070
@@ -250,10 +250,10 @@ class SparseLUSolver {
250250 // permutation
251251 u = (block_perm.p * u.matrix ()).array ();
252252 // forward substitution, per row in u
253- for (Idx br = 0 ; br < block_size; ++br ) {
254- for (Idx bc = 0 ; bc < br ; ++bc ) {
253+ for (Idx block_row = 0 ; block_row < block_size; ++block_row ) {
254+ for (Idx block_col = 0 ; block_col < block_row ; ++block_col ) {
255255 // forward substract
256- u.row (br ) -= pivot (br, bc ) * u.row (bc );
256+ u.row (block_row ) -= pivot (block_row, block_col ) * u.row (block_col );
257257 }
258258 }
259259 }
@@ -287,12 +287,12 @@ class SparseLUSolver {
287287 // l * u = a
288288 // l0 * u00 = a0
289289 // l0 * u01 + l1 * u11 = a1
290- for (Idx bc = 0 ; bc < block_size; ++bc ) {
291- for (Idx br = 0 ; br < bc ; ++br ) {
292- l.col (bc ) -= pivot (br, bc ) * l.col (br );
290+ for (Idx block_col = 0 ; block_col < block_size; ++block_col ) {
291+ for (Idx block_row = 0 ; block_row < block_col ; ++block_row ) {
292+ l.col (block_col ) -= pivot (block_row, block_col ) * l.col (block_row );
293293 }
294294 // divide diagonal
295- l.col (bc ) = l.col (bc ) / pivot (bc, bc );
295+ l.col (block_col ) = l.col (block_col ) / pivot (block_col, block_col );
296296 }
297297 }
298298 else {
0 commit comments