diff --git a/src/block_krylov_solvers.jl b/src/block_krylov_solvers.jl
index c6c128aeb..21d361f59 100644
--- a/src/block_krylov_solvers.jl
+++ b/src/block_krylov_solvers.jl
@@ -79,6 +79,7 @@ for (KS, fun, nsol, nA, nAt, warm_start) in [
     Atprod(solver :: $KS) = $nAt * solver.stats.niter
     nsolution(solver :: $KS) = $nsol
     if $nsol == 1
+      solution_tuple(solver :: $KS) = (solver.X,)
       solution(solver :: $KS) = solver.X
       solution(solver :: $KS, p :: Integer) = (p == 1) ? solution(solver) : error("solution(solver) has only one output.")
     end
diff --git a/src/krylov_solve.jl b/src/krylov_solve.jl
index c95a3705e..5429e831c 100644
--- a/src/krylov_solve.jl
+++ b/src/krylov_solve.jl
@@ -56,7 +56,7 @@ for (KS, fun, fun2, args, def_args, optargs, def_optargs, kwargs, def_kwargs) in
         timemax -= elapsed_time
         $(fun)(solver, $(args...); $(kwargs...))
         solver.stats.timer += elapsed_time
-        return ((solution(solver),)..., solver.stats)
+        return (solution_tuple(solver)..., solver.stats)
       end
 
       if !isempty($optargs)
@@ -68,7 +68,7 @@ for (KS, fun, fun2, args, def_args, optargs, def_optargs, kwargs, def_kwargs) in
           timemax -= elapsed_time
           $(fun)(solver, $(args...); $(kwargs...))
           solver.stats.timer += elapsed_time
-          return ((solution(solver),)..., solver.stats)
+          return (solution_tuple(solver)..., solver.stats)
         end
       end
     end
@@ -107,7 +107,7 @@ for (KS, fun, fun2, args, def_args, optargs, def_optargs, kwargs, def_kwargs) in
       timemax -= elapsed_time
       $(fun)(solver, $(args...); $(kwargs...))
       solver.stats.timer += elapsed_time
-      return ((solution(solver),)..., solver.stats)
+      return (solution_tuple(solver)..., solver.stats)
     end
 
     if !isempty($optargs)
@@ -119,7 +119,7 @@ for (KS, fun, fun2, args, def_args, optargs, def_optargs, kwargs, def_kwargs) in
         timemax -= elapsed_time
         $(fun)(solver, $(args...); $(kwargs...))
         solver.stats.timer += elapsed_time
-        return ((solution(solver),)..., solver.stats)
+        return (solution_tuple(solver)..., solver.stats)
       end
     end
 
diff --git a/src/krylov_solvers.jl b/src/krylov_solvers.jl
index d9f9fd9b1..aca001b32 100644
--- a/src/krylov_solvers.jl
+++ b/src/krylov_solvers.jl
@@ -1944,11 +1944,13 @@ for (KS, fun, nsol, nA, nAt, warm_start) in [
     end
     nsolution(solver :: $KS) = $nsol
     if $nsol == 1
+      solution_tuple(solver :: $KS) = (solver.x,)
       solution(solver :: $KS) = solver.x
       solution(solver :: $KS, p :: Integer) = (p == 1) ? solution(solver) : error("solution(solver) has only one output.")
     end
     if $nsol == 2
-      solution(solver :: $KS) = solver.x, solver.y
+      solution_tuple(solver :: $KS) = (solver.x, solver.y)
+      solution(solver :: $KS) = (solver.x, solver.y)
       solution(solver :: $KS, p :: Integer) = (1 ≤ p ≤ 2) ? solution(solver)[p] : error("solution(solver) has only two outputs.")
     end
     if $KS ∈ (BilqrSolver, TrilqrSolver)