Add a new file krylov_show.jl

Author: amontoison

src/Krylov.jl

@@ -12,6 +12,8 @@ include("block_krylov_utils.jl")
 include("block_krylov_processes.jl")
 include("block_krylov_workspaces.jl")
 
+include("krylov_show.jl")
+
 include("block_minres.jl")
 include("block_gmres.jl")
 

src/block_krylov_workspaces.jl

@@ -127,100 +127,3 @@ function BlockGmresWorkspace(A, B; memory::Integer = 5)
   SV = matrix_to_vector(SM)
   BlockGmresWorkspace(m, n, p, SV, SM; memory)
 end
-
-for (KS, fun, nsol, nA, nAt, warm_start) in [
-  (:BlockMinresWorkspace, :block_minres!, 1, 1, 0, true)
-  (:BlockGmresWorkspace , :block_gmres! , 1, 1, 0, true)
-]
-  @eval begin
-    krylov_elapsed_time(workspace :: $KS) = workspace.stats.timer
-    size(workspace :: $KS) = workspace.m, workspace.n
-    nrhs(workspace :: $KS) = workspace.p
-    statistics(workspace :: $KS) = workspace.stats
-    niterations(workspace :: $KS) = workspace.stats.niter
-    Aprod(workspace :: $KS) = $nA * workspace.stats.niter
-    Atprod(workspace :: $KS) = $nAt * workspace.stats.niter
-    nsolution(workspace :: $KS) = $nsol
-    if $nsol == 1
-      solution(workspace :: $KS) = workspace.X
-      solution(workspace :: $KS, p :: Integer) = (p == 1) ? solution(workspace) : error("solution(workspace) has only one output.")
-      results(workspace :: $KS) = (workspace.X, workspace.stats)
-    end
-    issolved(workspace :: $KS) = workspace.stats.solved
-    if $warm_start
-      function warm_start!(workspace :: $KS, X0)
-        n2, p2 = size(X0)
-        SM = typeof(workspace.X)
-        (workspace.n == n2 && workspace.p == p2) || error("X0 should have size ($n, $p)")
-        allocate_if(true, workspace, :ΔX, SM, workspace.n, workspace.p)
-        copyto!(workspace.ΔX, X0)
-        workspace.warm_start = true
-        return workspace
-      end
-    end
-  end
-end
-
-function ksizeof(attribute)
-  if isa(attribute, Vector{<:AbstractVector}) && !isempty(attribute)
-    # A vector of vectors is a vector of pointers in Julia.
-    # All vectors inside a vector have the same size in Krylov.jl
-    size_attribute = sizeof(attribute) + length(attribute) * ksizeof(attribute[1])
-  else
-    size_attribute = sizeof(attribute)
-  end
-  return size_attribute
-end
-
-function sizeof(stats_workspace :: Union{KrylovStats, KrylovWorkspace, BlockKrylovWorkspace})
-  type = typeof(stats_workspace)
-  nfields = fieldcount(type)
-  storage = 0
-  for i = 1:nfields
-    field_i = getfield(stats_workspace, i)
-    size_i = ksizeof(field_i)
-    storage += size_i
-  end
-  return storage
-end
-
-"""
-    show(io, workspace; show_stats=true)
-
-Statistics of `workspace` are displayed if `show_stats` is set to true.
-"""
-function show(io :: IO, workspace :: Union{KrylovWorkspace{T,FC,S}, BlockKrylovWorkspace{T,FC,S}}; show_stats :: Bool=true) where {T <: AbstractFloat, FC <: FloatOrComplex{T}, S <: AbstractVector{FC}}
-  type_workspace = typeof(workspace)
-  name_workspace = string(type_workspace.name.name)
-  name_stats = string(typeof(workspace.stats).name.name)
-  nbytes = sizeof(workspace)
-  storage = format_bytes(nbytes)
-  architecture = S <: Vector ? "CPU" : "GPU"
-  l1 = max(length(name_workspace), length(string(FC)) + 11)  # length("Precision: ") = 11
-  nchar = type_workspace <: Union{CgLanczosShiftWorkspace, FomWorkspace, DiomWorkspace, DqgmresWorkspace, GmresWorkspace, FgmresWorkspace, GpmrWorkspace, BlockGmresWorkspace} ? 8 : 0  # length("Vector{}") = 8
-  l2 = max(ndigits(workspace.m) + 7, length(architecture) + 14, length(string(S)) + nchar)  # length("nrows: ") = 7 and length("Architecture: ") = 14
-  l2 = max(l2, length(name_stats) + 2 + length(string(T)))  # length("{}") = 2
-  l3 = max(ndigits(workspace.n) + 7, length(storage) + 9)  # length("Storage: ") = 9 and length("cols: ") = 7
-  format = Printf.Format("│%$(l1)s│%$(l2)s│%$(l3)s│\n")
-  format2 = Printf.Format("│%$(l1+1)s│%$(l2)s│%$(l3)s│\n")
-  @printf(io, "┌%s┬%s┬%s┐\n", "─"^l1, "─"^l2, "─"^l3)
-  Printf.format(io, format, "$(name_workspace)", "nrows: $(workspace.m)", "ncols: $(workspace.n)")
-  @printf(io, "├%s┼%s┼%s┤\n", "─"^l1, "─"^l2, "─"^l3)
-  Printf.format(io, format, "Precision: $FC", "Architecture: $architecture","Storage: $storage")
-  @printf(io, "├%s┼%s┼%s┤\n", "─"^l1, "─"^l2, "─"^l3)
-  Printf.format(io, format, "Attribute", "Type", "Size")
-  @printf(io, "├%s┼%s┼%s┤\n", "─"^l1, "─"^l2, "─"^l3)
-  for i=1:fieldcount(type_workspace)
-    name_i = fieldname(type_workspace, i)
-    type_i = fieldtype(type_workspace, i)
-    field_i = getfield(workspace, name_i)
-    size_i = ksizeof(field_i)
-    (size_i ≠ 0) && Printf.format(io, format, string(name_i), type_i, format_bytes(size_i))
-  end
-  @printf(io, "└%s┴%s┴%s┘\n","─"^l1,"─"^l2,"─"^l3)
-  if show_stats
-    @printf(io, "\n")
-    show(io, workspace.stats)
-  end
-  return nothing
-end

src/interface.jl

@@ -26,6 +26,11 @@
 
 export krylov_workspace, krylov_solve, krylov_solve!
 
+export solution, nsolution, statistics, issolved, issolved_primal, issolved_dual
+export niterations, Aprod, Atprod, Bprod, warm_start!, krylov_elapsed_time
+
+import Base.size
+
 """
     workspace = krylov_workspace(Val(method), args...; kwargs...)
 
@@ -50,6 +55,79 @@ The argument `workspace` must be a subtype of [`KrylovWorkspace`](@ref) or [`Blo
 """
 function krylov_solve! end
 
+"""
+    solution(workspace)
+
+Return the solution(s) stored in the `workspace`.
+Optionally you can specify which solution you want to recover,
+`solution(workspace, 1)` returns `x` and `solution(workspace, 2)` returns `y`.
+"""
+function solution end
+
+"""
+    nsolution(workspace)
+
+Return the number of outputs of `solution(workspace)`.
+"""
+function nsolution end
+
+"""
+    statistics(workspace)
+
+Return the statistics stored in the `workspace`.
+"""
+function statistics end
+
+"""
+    issolved(workspace)
+
+Return a boolean that determines whether the Krylov method associated to `workspace` succeeded.
+"""
+function issolved end
+
+"""
+    niterations(workspace)
+
+Return the number of iterations performed by the Krylov method associated to `workspace`.
+"""
+function niterations end
+
+"""
+    Aprod(workspace)
+
+Return the number of operator-vector products with `A` performed by the Krylov method associated to `workspace`.
+"""
+function Aprod end
+
+"""
+    Atprod(workspace)
+
+Return the number of operator-vector products with `A'` performed by the Krylov method associated to `workspace`.
+"""
+function Atprod end
+
+"""
+    krylov_elapsed_time(workspace)
+
+Return the time elapsed (in seconds) during the last call to the Krylov solver associated with `workspace`.
+"""
+function krylov_elapsed_time end
+
+"""
+    results(workspace)
+
+Return a tuple containing the solution(s) and the statistics associated with the `workspace`.
+Allows retrieving the output arguments of an out-of-place method from the in-place method.
+
+For example, instead of `x, stats = cg(A, b)`, you can use:
+```julia
+workspace = CgWorkspace(A, b)
+cg!(workspace, A, b)
+x, stats = results(workspace)
+```
+"""
+function results end
+
 # Krylov methods
 for (workspace, krylov, args, def_args, optargs, def_optargs, kwargs, def_kwargs) in [
   (:LsmrWorkspace     , :lsmr      , args_lsmr      , def_args_lsmr      , ()                , ()                    , kwargs_lsmr      , def_kwargs_lsmr      )
@@ -333,3 +411,129 @@ for (workspace, krylov, args, def_args, optargs, def_optargs, kwargs, def_kwargs
     end
   end
 end
+
+for (KS, fun, nsol, nA, nAt, warm_start) in [
+  (:CarWorkspace      , :car!       , 1, 1, 0, true )
+  (:LsmrWorkspace     , :lsmr!      , 1, 1, 1, false)
+  (:CgsWorkspace      , :cgs!       , 1, 2, 0, true )
+  (:UsymlqWorkspace   , :usymlq!    , 1, 1, 1, true )
+  (:LnlqWorkspace     , :lnlq!      , 2, 1, 1, false)
+  (:BicgstabWorkspace , :bicgstab!  , 1, 2, 0, true )
+  (:CrlsWorkspace     , :crls!      , 1, 1, 1, false)
+  (:LsqrWorkspace     , :lsqr!      , 1, 1, 1, false)
+  (:MinresWorkspace   , :minres!    , 1, 1, 0, true )
+  (:MinaresWorkspace  , :minares!   , 1, 1, 0, true )
+  (:CgneWorkspace     , :cgne!      , 1, 1, 1, false)
+  (:DqgmresWorkspace  , :dqgmres!   , 1, 1, 0, true )
+  (:SymmlqWorkspace   , :symmlq!    , 1, 1, 0, true )
+  (:TrimrWorkspace    , :trimr!     , 2, 1, 1, true )
+  (:UsymqrWorkspace   , :usymqr!    , 1, 1, 1, true )
+  (:BilqrWorkspace    , :bilqr!     , 2, 1, 1, true )
+  (:CrWorkspace       , :cr!        , 1, 1, 0, true )
+  (:CraigmrWorkspace  , :craigmr!   , 2, 1, 1, false)
+  (:TricgWorkspace    , :tricg!     , 2, 1, 1, true )
+  (:CraigWorkspace    , :craig!     , 2, 1, 1, false)
+  (:DiomWorkspace     , :diom!      , 1, 1, 0, true )
+  (:LslqWorkspace     , :lslq!      , 1, 1, 1, false)
+  (:TrilqrWorkspace   , :trilqr!    , 2, 1, 1, true )
+  (:CrmrWorkspace     , :crmr!      , 1, 1, 1, false)
+  (:CgWorkspace       , :cg!        , 1, 1, 0, true )
+  (:CglsWorkspace     , :cgls!      , 1, 1, 1, false)
+  (:CgLanczosWorkspace, :cg_lanczos!, 1, 1, 0, true )
+  (:BilqWorkspace     , :bilq!      , 1, 1, 1, true )
+  (:MinresQlpWorkspace, :minres_qlp!, 1, 1, 0, true )
+  (:QmrWorkspace      , :qmr!       , 1, 1, 1, true )
+  (:GmresWorkspace    , :gmres!     , 1, 1, 0, true )
+  (:FgmresWorkspace   , :fgmres!    , 1, 1, 0, true )
+  (:FomWorkspace      , :fom!       , 1, 1, 0, true )
+  (:GpmrWorkspace     , :gpmr!      , 2, 1, 0, true )
+  (:CgLanczosShiftWorkspace  , :cg_lanczos_shift!  , 1, 1, 0, false)
+  (:CglsLanczosShiftWorkspace, :cgls_lanczos_shift!, 1, 1, 1, false)
+]
+  @eval begin
+    krylov_elapsed_time(workspace :: $KS) = workspace.stats.timer
+    size(workspace :: $KS) = workspace.m, workspace.n
+    statistics(workspace :: $KS) = workspace.stats
+    niterations(workspace :: $KS) = workspace.stats.niter
+    Aprod(workspace :: $KS) = $nA * workspace.stats.niter
+    Atprod(workspace :: $KS) = $nAt * workspace.stats.niter
+    if $KS == GpmrWorkspace
+      Bprod(workspace :: $KS) = workspace.stats.niter
+    end
+    nsolution(workspace :: $KS) = $nsol
+    if $nsol == 1
+      solution(workspace :: $KS) = workspace.x
+      solution(workspace :: $KS, p :: Integer) = (p == 1) ? solution(workspace) : error("solution(workspace) has only one output.")
+      results(workspace :: $KS) = (workspace.x, workspace.stats)
+    end
+    if $nsol == 2
+      solution(workspace :: $KS) = (workspace.x, workspace.y)
+      solution(workspace :: $KS, p :: Integer) = (1 ≤ p ≤ 2) ? solution(workspace)[p] : error("solution(workspace) has only two outputs.")
+      results(workspace :: $KS) = (workspace.x, workspace.y, workspace.stats)
+    end
+    if $KS ∈ (BilqrWorkspace, TrilqrWorkspace)
+      issolved_primal(workspace :: $KS) = workspace.stats.solved_primal
+      issolved_dual(workspace :: $KS) = workspace.stats.solved_dual
+      issolved(workspace :: $KS) = issolved_primal(workspace) && issolved_dual(workspace)
+    else
+      issolved(workspace :: $KS) = workspace.stats.solved
+    end
+    if $warm_start
+      if $KS in (BilqrWorkspace, TrilqrWorkspace, TricgWorkspace, TrimrWorkspace, GpmrWorkspace)
+        function warm_start!(workspace :: $KS, x0, y0)
+          length(x0) == workspace.n || error("x0 should have size $n")
+          length(y0) == workspace.m || error("y0 should have size $m")
+          S = typeof(workspace.x)
+          allocate_if(true, workspace, :Δx, S, workspace.x)  # The length of Δx is n
+          allocate_if(true, workspace, :Δy, S, workspace.y)  # The length of Δy is m
+          kcopy!(workspace.n, workspace.Δx, x0)
+          kcopy!(workspace.m, workspace.Δy, y0)
+          workspace.warm_start = true
+          return workspace
+        end
+      else
+        function warm_start!(workspace :: $KS, x0)
+          S = typeof(workspace.x)
+          length(x0) == workspace.n || error("x0 should have size $n")
+          allocate_if(true, workspace, :Δx, S, workspace.x)  # The length of Δx is n
+          kcopy!(workspace.n, workspace.Δx, x0)
+          workspace.warm_start = true
+          return workspace
+        end
+      end
+    end
+  end
+end
+
+for (KS, fun, nsol, nA, nAt, warm_start) in [
+  (:BlockMinresWorkspace, :block_minres!, 1, 1, 0, true)
+  (:BlockGmresWorkspace , :block_gmres! , 1, 1, 0, true)
+]
+  @eval begin
+    krylov_elapsed_time(workspace :: $KS) = workspace.stats.timer
+    size(workspace :: $KS) = workspace.m, workspace.n
+    nrhs(workspace :: $KS) = workspace.p
+    statistics(workspace :: $KS) = workspace.stats
+    niterations(workspace :: $KS) = workspace.stats.niter
+    Aprod(workspace :: $KS) = $nA * workspace.stats.niter
+    Atprod(workspace :: $KS) = $nAt * workspace.stats.niter
+    nsolution(workspace :: $KS) = $nsol
+    if $nsol == 1
+      solution(workspace :: $KS) = workspace.X
+      solution(workspace :: $KS, p :: Integer) = (p == 1) ? solution(workspace) : error("solution(workspace) has only one output.")
+      results(workspace :: $KS) = (workspace.X, workspace.stats)
+    end
+    issolved(workspace :: $KS) = workspace.stats.solved
+    if $warm_start
+      function warm_start!(workspace :: $KS, X0)
+        n2, p2 = size(X0)
+        SM = typeof(workspace.X)
+        (workspace.n == n2 && workspace.p == p2) || error("X0 should have size ($n, $p)")
+        allocate_if(true, workspace, :ΔX, SM, workspace.n, workspace.p)
+        copyto!(workspace.ΔX, X0)
+        workspace.warm_start = true
+        return workspace
+      end
+    end
+  end
+end