Fix SolveLower loop bounds in runtime BlockCholeskyDecomposition

[CodeReclaimers]

Summary

• SolveLower inner loops iterated up to NumBlocks instead of BlockSize, accessing out-of-bounds block elements when the two differ. The compile-time specialization (line 278-280) correctly uses BlockSize.

Test plan

• Code inspection: the runtime SolveLower (line 463-465) now matches the compile-time version (line 278-280) and the sibling SolveUpper (line 488-490), all of which use BlockSize.

Test code

// Verify SolveLower produces correct results when BlockSize != NumBlocks.
// Before the fix, the inner loops iterate over NumBlocks elements instead
// of BlockSize elements, producing wrong results.

#include <Mathematics/CholeskyDecomposition.h>
#include <cmath>
#include <cstdlib>
#include <iostream>

using namespace gte;

int main()
{
    // BlockSize=2, NumBlocks=3 => 6x6 matrix.
    // These must differ to expose the bug.
    BlockCholeskyDecomposition<double> decomp(2, 3);

    // Build a 6x6 positive definite matrix A = I + ones*ones^T
    // (identity plus rank-1 update, guaranteed positive definite).
    int32_t const N = 6;
    GMatrix<double> A(N, N);
    for (int32_t r = 0; r < N; ++r)
        for (int32_t c = 0; c < N; ++c)
            A(r, c) = (r == c) ? 2.0 : 1.0;

    // Convert to block matrix and factor.
    BlockCholeskyDecomposition<double>::BlockMatrix ABlock;
    decomp.Convert(A, ABlock);
    decomp.Factor(ABlock);

    // Set up RHS: B = [1, 2, 3, 4, 5, 6]
    BlockCholeskyDecomposition<double>::BlockVector BBlock(3);
    for (int32_t b = 0; b < 3; ++b)
    {
        BBlock[b].SetSize(2);
        BBlock[b][0] = 2.0 * b + 1.0;
        BBlock[b][1] = 2.0 * b + 2.0;
    }

    // Solve L*Y = B
    BlockCholeskyDecomposition<double>::BlockVector YBlock = BBlock;
    decomp.SolveLower(ABlock, YBlock);

    // Solve L^T*X = Y
    decomp.SolveUpper(ABlock, YBlock);

    // Convert solution back and verify A*X ≈ B.
    GVector<double> X(N);
    for (int32_t b = 0; b < 3; ++b)
    {
        X[2 * b] = YBlock[b][0];
        X[2 * b + 1] = YBlock[b][1];
    }

    double maxErr = 0.0;
    for (int32_t r = 0; r < N; ++r)
    {
        double ax = 0.0;
        for (int32_t c = 0; c < N; ++c)
            ax += A(r, c) * X[c];
        double b = static_cast<double>(r + 1);
        maxErr = std::max(maxErr, std::abs(ax - b));
    }

    if (maxErr > 1e-10)
    {
        std::cerr << "FAIL: max residual " << maxErr << std::endl;
        return EXIT_FAILURE;
    }
    std::cout << "PASS: max residual " << maxErr << std::endl;
    return EXIT_SUCCESS;
}

🤖 Generated with Claude Code