diff --git a/lib/LinearSolveAutotune/Project.toml b/lib/LinearSolveAutotune/Project.toml
index 92084ca27..83cc222db 100644
--- a/lib/LinearSolveAutotune/Project.toml
+++ b/lib/LinearSolveAutotune/Project.toml
@@ -16,6 +16,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 MKL_jll = "856f044c-d86e-5d09-b602-aeab76dc8ba7"
 Metal = "dde4c033-4e86-420c-a63e-0dd931031962"
+Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Preferences = "21216c6a-2e73-6563-6e65-726566657250"
 PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
@@ -41,6 +42,7 @@ LinearAlgebra = "1"
 LinearSolve = "3"
 MKL_jll = "2025.2.0"
 Metal = "1"
+Pkg = "1"
 Plots = "1"
 Preferences = "1"
 PrettyTables = "2"
diff --git a/lib/LinearSolveAutotune/src/benchmarking.jl b/lib/LinearSolveAutotune/src/benchmarking.jl
index f70be4557..cbfb38ef9 100644
--- a/lib/LinearSolveAutotune/src/benchmarking.jl
+++ b/lib/LinearSolveAutotune/src/benchmarking.jl
@@ -1,6 +1,7 @@
 # Core benchmarking functionality
 
 using ProgressMeter
+using LinearAlgebra
 
 """
     test_algorithm_compatibility(alg, eltype::Type, test_size::Int=4)
@@ -78,7 +79,8 @@ Benchmark the given algorithms across different matrix sizes and element types.
 Returns a DataFrame with results including element type information.
 """
 function benchmark_algorithms(matrix_sizes, algorithms, alg_names, eltypes;
-        samples = 5, seconds = 0.5, sizes = [:tiny, :small, :medium, :large])
+        samples = 5, seconds = 0.5, sizes = [:tiny, :small, :medium, :large],
+        check_correctness = true, correctness_tol = 1e-2)
 
     # Set benchmark parameters
     old_params = BenchmarkTools.DEFAULT_PARAMETERS
@@ -116,6 +118,18 @@ function benchmark_algorithms(matrix_sizes, algorithms, alg_names, eltypes;
                 A = rand(rng, eltype, n, n)
                 b = rand(rng, eltype, n)
                 u0 = rand(rng, eltype, n)
+                
+                # Compute reference solution with LUFactorization if correctness check is enabled
+                reference_solution = nothing
+                if check_correctness
+                    try
+                        ref_prob = LinearProblem(copy(A), copy(b); u0 = copy(u0))
+                        reference_solution = solve(ref_prob, LinearSolve.LUFactorization())
+                    catch e
+                        @warn "Failed to compute reference solution with LUFactorization for size $n, eltype $eltype: $e"
+                        check_correctness = false  # Disable for this size/type combination
+                    end
+                end
 
                 for (alg, name) in zip(compatible_algs, compatible_names)
                     # Update progress description
@@ -125,6 +139,7 @@ function benchmark_algorithms(matrix_sizes, algorithms, alg_names, eltypes;
                     gflops = 0.0
                     success = true
                     error_msg = ""
+                    passed_correctness = true
 
                     try
                         # Create the linear problem for this test
@@ -132,19 +147,37 @@ function benchmark_algorithms(matrix_sizes, algorithms, alg_names, eltypes;
                             u0 = copy(u0),
                             alias = LinearAliasSpecifier(alias_A = true, alias_b = true))
 
-                        # Warmup run
-                        solve(prob, alg)
-
-                        # Actual benchmark
-                        bench = @benchmark solve($prob, $alg) setup=(prob = LinearProblem(
-                            copy($A), copy($b);
-                            u0 = copy($u0),
-                            alias = LinearAliasSpecifier(alias_A = true, alias_b = true)))
-
-                        # Calculate GFLOPs
-                        min_time_sec = minimum(bench.times) / 1e9
-                        flops = luflop(n, n)
-                        gflops = flops / min_time_sec / 1e9
+                        # Warmup run and correctness check
+                        warmup_sol = solve(prob, alg)
+                        
+                        # Check correctness if reference solution is available
+                        if check_correctness && reference_solution !== nothing
+                            # Compute relative error
+                            rel_error = norm(warmup_sol.u - reference_solution.u) / norm(reference_solution.u)
+                            
+                            if rel_error > correctness_tol
+                                passed_correctness = false
+                                @warn "Algorithm $name failed correctness check for size $n, eltype $eltype. " *
+                                      "Relative error: $(round(rel_error, sigdigits=3)) > tolerance: $correctness_tol. " *
+                                      "Algorithm will be excluded from results."
+                                success = false
+                                error_msg = "Failed correctness check (rel_error = $(round(rel_error, sigdigits=3)))"
+                            end
+                        end
+                        
+                        # Only benchmark if correctness check passed
+                        if passed_correctness
+                            # Actual benchmark
+                            bench = @benchmark solve($prob, $alg) setup=(prob = LinearProblem(
+                                copy($A), copy($b);
+                                u0 = copy($u0),
+                                alias = LinearAliasSpecifier(alias_A = true, alias_b = true)))
+
+                            # Calculate GFLOPs
+                            min_time_sec = minimum(bench.times) / 1e9
+                            flops = luflop(n, n)
+                            gflops = flops / min_time_sec / 1e9
+                        end
 
                     catch e
                         success = false
diff --git a/lib/LinearSolveAutotune/src/gpu_detection.jl b/lib/LinearSolveAutotune/src/gpu_detection.jl
index 8a04fc8e0..40a52e307 100644
--- a/lib/LinearSolveAutotune/src/gpu_detection.jl
+++ b/lib/LinearSolveAutotune/src/gpu_detection.jl
@@ -1,6 +1,7 @@
 # GPU hardware and package detection
 
 using CPUSummary
+using Pkg
 
 """
     is_cuda_available()
@@ -116,9 +117,107 @@ function get_system_info()
         info["apple_accelerate_available"] = false
     end
 
+    # Add package versions
+    info["package_versions"] = get_package_versions()
+    
     return info
 end
 
+"""
+    get_package_versions()
+
+Get versions of LinearSolve-related packages and their dependencies.
+Returns a Dict with package names and versions.
+"""
+function get_package_versions()
+    versions = Dict{String, String}()
+    
+    # Get the current project's dependencies
+    deps = Pkg.dependencies()
+    
+    # List of packages we're interested in tracking
+    important_packages = [
+        "LinearSolve",
+        "LinearSolveAutotune",
+        "RecursiveFactorization",
+        "CUDA",
+        "Metal",
+        "MKL_jll",
+        "BLISBLAS",
+        "AppleAccelerate",
+        "SparseArrays",
+        "KLU",
+        "Pardiso",
+        "MKLPardiso",
+        "BandedMatrices",
+        "FastLapackInterface",
+        "HYPRE",
+        "IterativeSolvers",
+        "Krylov",
+        "KrylovKit",
+        "LinearAlgebra"
+    ]
+    
+    # Also track JLL packages for BLAS libraries
+    jll_packages = [
+        "MKL_jll",
+        "OpenBLAS_jll",
+        "OpenBLAS32_jll",
+        "blis_jll",
+        "LAPACK_jll",
+        "CompilerSupportLibraries_jll"
+    ]
+    
+    all_packages = union(important_packages, jll_packages)
+    
+    # Iterate through dependencies and collect versions
+    for (uuid, dep) in deps
+        if dep.name in all_packages
+            if dep.version !== nothing
+                versions[dep.name] = string(dep.version)
+            else
+                # Try to get version from the package itself if loaded
+                try
+                    pkg_module = Base.loaded_modules[Base.PkgId(uuid, dep.name)]
+                    if isdefined(pkg_module, :version)
+                        versions[dep.name] = string(pkg_module.version)
+                    else
+                        versions[dep.name] = "unknown"
+                    end
+                catch
+                    versions[dep.name] = "unknown"
+                end
+            end
+        end
+    end
+    
+    # Try to get Julia's LinearAlgebra stdlib version
+    try
+        versions["LinearAlgebra"] = string(VERSION)  # Stdlib version matches Julia
+    catch
+        versions["LinearAlgebra"] = "stdlib"
+    end
+    
+    # Get BLAS configuration info
+    try
+        blas_config = LinearAlgebra.BLAS.get_config()
+        if hasfield(typeof(blas_config), :loaded_libs)
+            for lib in blas_config.loaded_libs
+                if hasfield(typeof(lib), :libname)
+                    lib_name = basename(string(lib.libname))
+                    # Extract version info if available
+                    versions["BLAS_lib"] = lib_name
+                end
+            end
+        end
+    catch
+        # Fallback for older Julia versions
+        versions["BLAS_vendor"] = string(LinearAlgebra.BLAS.vendor())
+    end
+    
+    return versions
+end
+
 """
     get_detailed_system_info()
 
diff --git a/lib/LinearSolveAutotune/src/telemetry.jl b/lib/LinearSolveAutotune/src/telemetry.jl
index dd49abf53..c1e16c79a 100644
--- a/lib/LinearSolveAutotune/src/telemetry.jl
+++ b/lib/LinearSolveAutotune/src/telemetry.jl
@@ -172,6 +172,21 @@ function format_system_info_markdown(system_info::Dict)
     push!(lines, "- **CUDA Available**: $(get(system_info, "cuda_available", get(system_info, "has_cuda", false)))")
     # Handle both "has_metal" and "metal_available" keys
     push!(lines, "- **Metal Available**: $(get(system_info, "metal_available", get(system_info, "has_metal", false)))")
+    
+    # Add package versions section
+    if haskey(system_info, "package_versions")
+        push!(lines, "")
+        push!(lines, "### Package Versions")
+        pkg_versions = system_info["package_versions"]
+        
+        # Sort packages for consistent display
+        sorted_packages = sort(collect(keys(pkg_versions)))
+        
+        for pkg_name in sorted_packages
+            version = pkg_versions[pkg_name]
+            push!(lines, "- **$pkg_name**: $version")
+        end
+    end
 
     return join(lines, "\n")
 end