From fc8677ba1fa0a5254f3072f3c19482cecdac17f5 Mon Sep 17 00:00:00 2001
From: Alexis Montoison <alexis.montoison@polymtl.ca>
Date: Fri, 9 Sep 2022 13:43:59 -0400
Subject: [PATCH] Test Krylov macros

---
 src/krylov_utils.jl | 151 ++++++++++++-----------
 test/test_aux.jl    | 288 ++++++++++++++++++++++++++------------------
 2 files changed, 246 insertions(+), 193 deletions(-)

diff --git a/src/krylov_utils.jl b/src/krylov_utils.jl
index c61bf2e5e..46c9d6cd6 100644
--- a/src/krylov_utils.jl
+++ b/src/krylov_utils.jl
@@ -175,6 +175,49 @@ function to_boundary(x :: Vector{T}, d :: Vector{T},
   return roots # `σ1` and `σ2`
 end
 
+"""
+    s = vec2str(x; ndisp)
+
+Display an array in the form
+
+    [ -3.0e-01 -5.1e-01  1.9e-01 ... -2.3e-01 -4.4e-01  2.4e-01 ]
+
+with (ndisp - 1)/2 elements on each side.
+"""
+function vec2str(x :: AbstractVector{T}; ndisp :: Int=7) where T <: Union{AbstractFloat, Missing}
+  n = length(x)
+  if n ≤ ndisp
+    ndisp = n
+    nside = n
+  else
+    nside = max(1, div(ndisp - 1, 2))
+  end
+  s = "["
+  i = 1
+  while i ≤ nside
+    if x[i] !== missing
+      s *= @sprintf("%8.1e ", x[i])
+    else
+      s *= " ✗✗✗✗ "
+    end
+      i += 1
+  end
+  if i ≤ div(n, 2)
+    s *= "... "
+  end
+  i = max(i, n - nside + 1)
+  while i ≤ n
+    if x[i] !== missing
+      s *= @sprintf("%8.1e ", x[i])
+    else
+      s *= " ✗✗✗✗ "
+    end
+    i += 1
+  end
+  s *= "]"
+  return s
+end
+
 """
     S = ktypeof(v)
 
@@ -209,76 +252,75 @@ end
 
 Create an AbstractVector of storage type `S` of length `n` only composed of zero.
 """
-@inline kzeros(S, n) = fill!(S(undef, n), zero(eltype(S)))
+kzeros(S, n) = fill!(S(undef, n), zero(eltype(S)))
 
 """
     v = kones(S, n)
 
 Create an AbstractVector of storage type `S` of length `n` only composed of one.
 """
-@inline kones(S, n) = fill!(S(undef, n), one(eltype(S)))
+kones(S, n) = fill!(S(undef, n), one(eltype(S)))
 
-@inline allocate_if(bool, solver, v, S, n) = bool && isempty(solver.:($v)) && (solver.:($v) = S(undef, n))
+allocate_if(bool, solver, v, S, n) = bool && isempty(solver.:($v)) && (solver.:($v) = S(undef, n))
 
-@inline kdisplay(iter, verbose) = (verbose > 0) && (mod(iter, verbose) == 0)
+kdisplay(iter, verbose) = (verbose > 0) && (mod(iter, verbose) == 0)
 
-@inline mulorldiv!(y, P, x, ldiv::Bool) = ldiv ? ldiv!(y, P, x) : mul!(y, P, x)
+mulorldiv!(y, P, x, ldiv::Bool) = ldiv ? ldiv!(y, P, x) : mul!(y, P, x)
 
-@inline krylov_dot(n :: Integer, x :: Vector{T}, dx :: Integer, y :: Vector{T}, dy :: Integer) where T <: BLAS.BlasReal = BLAS.dot(n, x, dx, y, dy)
-@inline krylov_dot(n :: Integer, x :: Vector{T}, dx :: Integer, y :: Vector{T}, dy :: Integer) where T <: BLAS.BlasComplex = BLAS.dotc(n, x, dx, y, dy)
-@inline krylov_dot(n :: Integer, x :: AbstractVector{T}, dx :: Integer, y :: AbstractVector{T}, dy :: Integer) where T <: Number = dot(x, y)
+kdot(n :: Integer, x :: Vector{T}, dx :: Integer, y :: Vector{T}, dy :: Integer) where T <: BLAS.BlasReal = BLAS.dot(n, x, dx, y, dy)
+kdot(n :: Integer, x :: Vector{T}, dx :: Integer, y :: Vector{T}, dy :: Integer) where T <: BLAS.BlasComplex = BLAS.dotc(n, x, dx, y, dy)
+kdot(n :: Integer, x :: AbstractVector{T}, dx :: Integer, y :: AbstractVector{T}, dy :: Integer) where T <: FloatOrComplex = dot(x, y)
 
-@inline krylov_dotr(n :: Integer, x :: AbstractVector{T}, dx :: Integer, y :: AbstractVector{T}, dy :: Integer) where T <: AbstractFloat = krylov_dot(n, x, dx, y, dy)
-@inline krylov_dotr(n :: Integer, x :: AbstractVector{Complex{T}}, dx :: Integer, y :: AbstractVector{Complex{T}}, dy :: Integer) where T <: AbstractFloat = real(krylov_dot(n, x, dx, y, dy))
+kdotr(n :: Integer, x :: AbstractVector{T}, dx :: Integer, y :: AbstractVector{T}, dy :: Integer) where T <: AbstractFloat = kdot(n, x, dx, y, dy)
+kdotr(n :: Integer, x :: AbstractVector{Complex{T}}, dx :: Integer, y :: AbstractVector{Complex{T}}, dy :: Integer) where T <: AbstractFloat = real(kdot(n, x, dx, y, dy))
 
-@inline krylov_norm2(n :: Integer, x :: Vector{T}, dx :: Integer) where T <: BLAS.BlasFloat = BLAS.nrm2(n, x, dx)
-@inline krylov_norm2(n :: Integer, x :: AbstractVector{T}, dx :: Integer) where T <: Number = norm(x)
+knrm2(n :: Integer, x :: Vector{T}, dx :: Integer) where T <: BLAS.BlasFloat = BLAS.nrm2(n, x, dx)
+knrm2(n :: Integer, x :: AbstractVector{T}, dx :: Integer) where T <: FloatOrComplex = norm(x)
 
-@inline krylov_scal!(n :: Integer, s :: T, x :: Vector{T}, dx :: Integer) where T <: BLAS.BlasFloat = BLAS.scal!(n, s, x, dx)
-@inline krylov_scal!(n :: Integer, s :: T, x :: AbstractVector{T}, dx :: Integer) where T <: Number = (x .*= s)
-@inline krylov_scal!(n :: Integer, s :: T, x :: AbstractVector{Complex{T}}, dx :: Integer) where T <: AbstractFloat = krylov_scal!(n, Complex{T}(s), x, dx)
+kscal!(n :: Integer, s :: T, x :: Vector{T}, dx :: Integer) where T <: BLAS.BlasFloat = BLAS.scal!(n, s, x, dx)
+kscal!(n :: Integer, s :: T, x :: AbstractVector{T}, dx :: Integer) where T <: FloatOrComplex = (x .*= s)
+kscal!(n :: Integer, s :: T, x :: AbstractVector{Complex{T}}, dx :: Integer) where T <: AbstractFloat = kscal!(n, Complex{T}(s), x, dx)
 
-@inline krylov_axpy!(n :: Integer, s :: T, x :: Vector{T}, dx :: Integer, y :: Vector{T}, dy :: Integer) where T <: BLAS.BlasFloat = BLAS.axpy!(n, s, x, dx, y, dy)
-@inline krylov_axpy!(n :: Integer, s :: T, x :: AbstractVector{T}, dx :: Integer, y :: AbstractVector{T}, dy :: Integer) where T <: Number = axpy!(s, x, y)
-@inline krylov_axpy!(n :: Integer, s :: T, x :: AbstractVector{Complex{T}}, dx :: Integer, y :: AbstractVector{Complex{T}}, dy :: Integer) where T <: AbstractFloat = krylov_axpy!(n, Complex{T}(s), x, dx, y, dy)
+kaxpy!(n :: Integer, s :: T, x :: Vector{T}, dx :: Integer, y :: Vector{T}, dy :: Integer) where T <: BLAS.BlasFloat = BLAS.axpy!(n, s, x, dx, y, dy)
+kaxpy!(n :: Integer, s :: T, x :: AbstractVector{T}, dx :: Integer, y :: AbstractVector{T}, dy :: Integer) where T <: FloatOrComplex = axpy!(s, x, y)
+kaxpy!(n :: Integer, s :: T, x :: AbstractVector{Complex{T}}, dx :: Integer, y :: AbstractVector{Complex{T}}, dy :: Integer) where T <: AbstractFloat = kaxpy!(n, Complex{T}(s), x, dx, y, dy)
 
-@inline krylov_axpby!(n :: Integer, s :: T, x :: Vector{T}, dx :: Integer, t :: T, y :: Vector{T}, dy :: Integer) where T <: BLAS.BlasFloat = BLAS.axpby!(n, s, x, dx, t, y, dy)
-@inline krylov_axpby!(n :: Integer, s :: T, x :: AbstractVector{T}, dx :: Integer, t :: T, y :: AbstractVector{T}, dy :: Integer) where T <: Number = axpby!(s, x, t, y)
-@inline krylov_axpby!(n :: Integer, s :: T, x :: AbstractVector{Complex{T}}, dx :: Integer, t :: Complex{T}, y :: AbstractVector{Complex{T}}, dy :: Integer) where T <: AbstractFloat = krylov_axpby!(n, Complex{T}(s), x, dx, t, y, dy)
-@inline krylov_axpby!(n :: Integer, s :: Complex{T}, x :: AbstractVector{Complex{T}}, dx :: Integer, t :: T, y :: AbstractVector{Complex{T}}, dy :: Integer) where T <: AbstractFloat = krylov_axpby!(n, s, x, dx, Complex{T}(t), y, dy)
-@inline krylov_axpby!(n :: Integer, s :: T, x :: AbstractVector{Complex{T}}, dx :: Integer, t :: T, y :: AbstractVector{Complex{T}}, dy :: Integer) where T <: AbstractFloat = krylov_axpby!(n, Complex{T}(s), x, dx, Complex{T}(t), y, dy)
+kaxpby!(n :: Integer, s :: T, x :: Vector{T}, dx :: Integer, t :: T, y :: Vector{T}, dy :: Integer) where T <: BLAS.BlasFloat = BLAS.axpby!(n, s, x, dx, t, y, dy)
+kaxpby!(n :: Integer, s :: T, x :: AbstractVector{T}, dx :: Integer, t :: T, y :: AbstractVector{T}, dy :: Integer) where T <: FloatOrComplex = axpby!(s, x, t, y)
+kaxpby!(n :: Integer, s :: T, x :: AbstractVector{Complex{T}}, dx :: Integer, t :: Complex{T}, y :: AbstractVector{Complex{T}}, dy :: Integer) where T <: AbstractFloat = kaxpby!(n, Complex{T}(s), x, dx, t, y, dy)
+kaxpby!(n :: Integer, s :: Complex{T}, x :: AbstractVector{Complex{T}}, dx :: Integer, t :: T, y :: AbstractVector{Complex{T}}, dy :: Integer) where T <: AbstractFloat = kaxpby!(n, s, x, dx, Complex{T}(t), y, dy)
+kaxpby!(n :: Integer, s :: T, x :: AbstractVector{Complex{T}}, dx :: Integer, t :: T, y :: AbstractVector{Complex{T}}, dy :: Integer) where T <: AbstractFloat = kaxpby!(n, Complex{T}(s), x, dx, Complex{T}(t), y, dy)
 
-@inline krylov_copy!(n :: Integer, x :: Vector{T}, dx :: Integer, y :: Vector{T}, dy :: Integer) where T <: BLAS.BlasFloat = BLAS.blascopy!(n, x, dx, y, dy)
-@inline krylov_copy!(n :: Integer, x :: AbstractVector{T}, dx :: Integer, y :: AbstractVector{T}, dy :: Integer) where T <: Number = copyto!(y, x)
+kcopy!(n :: Integer, x :: Vector{T}, dx :: Integer, y :: Vector{T}, dy :: Integer) where T <: BLAS.BlasFloat = BLAS.blascopy!(n, x, dx, y, dy)
+kcopy!(n :: Integer, x :: AbstractVector{T}, dx :: Integer, y :: AbstractVector{T}, dy :: Integer) where T <: FloatOrComplex = copyto!(y, x)
 
 # the macros are just for readability, so we don't have to write the increments (always equal to 1)
-
 macro kdot(n, x, y)
-  return esc(:(krylov_dot($n, $x, 1, $y, 1)))
+  return esc(:(Krylov.kdot($n, $x, 1, $y, 1)))
 end
 
 macro kdotr(n, x, y)
-  return esc(:(krylov_dotr($n, $x, 1, $y, 1)))
+  return esc(:(Krylov.kdotr($n, $x, 1, $y, 1)))
 end
 
 macro knrm2(n, x)
-  return esc(:(krylov_norm2($n, $x, 1)))
+  return esc(:(Krylov.knrm2($n, $x, 1)))
 end
 
 macro kscal!(n, s, x)
-  return esc(:(krylov_scal!($n, $s, $x, 1)))
+  return esc(:(Krylov.kscal!($n, $s, $x, 1)))
 end
 
 macro kaxpy!(n, s, x, y)
-  return esc(:(krylov_axpy!($n, $s, $x, 1, $y, 1)))
+  return esc(:(Krylov.kaxpy!($n, $s, $x, 1, $y, 1)))
 end
 
 macro kaxpby!(n, s, x, t, y)
-  return esc(:(krylov_axpby!($n, $s, $x, 1, $t, $y, 1)))
+  return esc(:(Krylov.kaxpby!($n, $s, $x, 1, $t, $y, 1)))
 end
 
 macro kcopy!(n, x, y)
-  return esc(:(krylov_copy!($n, $x, 1, $y, 1)))
+  return esc(:(Krylov.kcopy!($n, $x, 1, $y, 1)))
 end
 
 macro kswap(x, y)
@@ -292,46 +334,3 @@ end
 macro kref!(n, x, y, c, s)
   return esc(:(reflect!($x, $y, $c, $s)))
 end
-
-"""
-    s = vec2str(x; ndisp)
-
-Display an array in the form
-
-    [ -3.0e-01 -5.1e-01  1.9e-01 ... -2.3e-01 -4.4e-01  2.4e-01 ]
-
-with (ndisp - 1)/2 elements on each side.
-"""
-function vec2str(x :: AbstractVector{T}; ndisp :: Int=7) where T <: Union{AbstractFloat, Missing}
-  n = length(x)
-  if n ≤ ndisp
-    ndisp = n
-    nside = n
-  else
-    nside = max(1, div(ndisp - 1, 2))
-  end
-  s = "["
-  i = 1
-  while i ≤ nside
-    if x[i] !== missing
-      s *= @sprintf("%8.1e ", x[i])
-    else
-      s *= " ✗✗✗✗ "
-    end
-      i += 1
-  end
-  if i ≤ div(n, 2)
-    s *= "... "
-  end
-  i = max(i, n - nside + 1)
-  while i ≤ n
-    if x[i] !== missing
-      s *= @sprintf("%8.1e ", x[i])
-    else
-      s *= " ✗✗✗✗ "
-    end
-    i += 1
-  end
-  s *= "]"
-  return s
-end
diff --git a/test/test_aux.jl b/test/test_aux.jl
index 11bdb7c2d..215ffc4b8 100644
--- a/test/test_aux.jl
+++ b/test/test_aux.jl
@@ -1,119 +1,173 @@
 @testset "aux" begin
-  # test Givens reflector corner cases
-  (c, s, ρ) = Krylov.sym_givens(0.0, 0.0)
-  @test (c == 1.0) && (s == 0.0) && (ρ == 0.0)
-
-  a = 3.14
-  (c, s, ρ) = Krylov.sym_givens(a, 0.0)
-  @test (c == 1.0) && (s == 0.0) && (ρ == a)
-  (c, s, ρ) = Krylov.sym_givens(-a, 0.0)
-  @test (c == -1.0) && (s == 0.0) && (ρ == a)
-
-  b = 3.14
-  (c, s, ρ) = Krylov.sym_givens(0.0, b)
-  @test (c == 0.0) && (s == 1.0) && (ρ == b)
-  (c, s, ρ) = Krylov.sym_givens(0.0, -b)
-  @test (c == 0.0) && (s == -1.0) && (ρ == b)
-
-  (c, s, ρ) = Krylov.sym_givens(Complex(0.0), Complex(0.0))
-  @test (c == 1.0) && (s == Complex(0.0)) && (ρ == Complex(0.0))
-
-  a = Complex(1.0, 1.0)
-  (c, s, ρ) = Krylov.sym_givens(a, Complex(0.0))
-  @test (c == 1.0) && (s == Complex(0.0)) && (ρ == a)
-  (c, s, ρ) = Krylov.sym_givens(-a, Complex(0.0))
-  @test (c == 1.0) && (s == Complex(0.0)) && (ρ == -a)
-
-  b = Complex(1.0, 1.0)
-  (c, s, ρ) = Krylov.sym_givens(Complex(0.0), b)
-  @test (c == 0.0) && (s == Complex(1.0)) && (ρ == b)
-  (c, s, ρ) = Krylov.sym_givens(Complex(0.0), -b)
-  @test (c == 0.0) && (s == Complex(1.0)) && (ρ == -b)
-
-  # test roots of a quadratic
-  roots = Krylov.roots_quadratic(0.0, 0.0, 0.0)
-  @test length(roots) == 1
-  @test roots[1] == 0.0
-
-  roots = Krylov.roots_quadratic(0.0, 0.0, 1.0)
-  @test length(roots) == 0
-
-  roots = Krylov.roots_quadratic(0.0, 3.14, -1.0)
-  @test length(roots) == 1
-  @test roots[1] == 1.0 / 3.14
-
-  roots = Krylov.roots_quadratic(1.0, 0.0, 1.0)
-  @test length(roots) == 0
-
-  roots = Krylov.roots_quadratic(1.0, 0.0, 0.0)
-  @test length(roots) == 2
-  @test roots[1] == 0.0
-  @test roots[2] == 0.0
-
-  roots = Krylov.roots_quadratic(1.0, 3.0, 2.0)
-  @test length(roots) == 2
-  @test roots[1] ≈ -2.0
-  @test roots[2] ≈ -1.0
-
-  roots = Krylov.roots_quadratic(1.0e+8, 1.0, 1.0)
-  @test length(roots) == 0
-
-  # ill-conditioned quadratic
-  roots = Krylov.roots_quadratic(-1.0e-8, 1.0e+5, 1.0, nitref=0)
-  @test length(roots) == 2
-  @test roots[1] == 1.0e+13
-  @test roots[2] == 0.0
-
-  # iterative refinement is crucial!
-  roots = Krylov.roots_quadratic(-1.0e-8, 1.0e+5, 1.0, nitref=1)
-  @test length(roots) == 2
-  @test roots[1] == 1.0e+13
-  @test roots[2] == -1.0e-05
-
-  # not ill-conditioned quadratic
-  roots = Krylov.roots_quadratic(-1.0e-7, 1.0, 1.0, nitref=0)
-  @test length(roots) == 2
-  @test isapprox(roots[1],  1.0e+7, rtol=1.0e-6)
-  @test isapprox(roots[2], -1.0, rtol=1.0e-6)
-
-  roots = Krylov.roots_quadratic(-1.0e-7, 1.0, 1.0, nitref=1)
-  @test length(roots) == 2
-  @test isapprox(roots[1], 1.0e+7, rtol=1.0e-6)
-  @test isapprox(roots[2], -1.0, rtol=1.0e-6)
-
-  # test trust-region boundary
-  x = ones(5)
-  d = ones(5); d[1:2:5] .= -1
-  @test_throws ErrorException Krylov.to_boundary(x, d, -1.0)
-  @test_throws ErrorException Krylov.to_boundary(x, d, 0.5)
-  @test_throws ErrorException Krylov.to_boundary(x, zeros(5), 1.0)
-  @test maximum(Krylov.to_boundary(x, d, 5.0)) ≈ 2.209975124224178
-  @test minimum(Krylov.to_boundary(x, d, 5.0)) ≈ -1.8099751242241782
-  @test maximum(Krylov.to_boundary(x, d, 5.0, flip=true)) ≈ 1.8099751242241782
-  @test minimum(Krylov.to_boundary(x, d, 5.0, flip=true)) ≈ -2.209975124224178
-
-  # test kzeros and kones
-  @test Krylov.kzeros(Vector{Float64}, 10) == zeros(10)
-  @test Krylov.kones(Vector{Float64}, 10) == ones(10)
-
-  # test ktypeof
-  a = rand(Float32, 10)
-  b = view(a, 4:8)
-  @test Krylov.ktypeof(a) == Vector{Float32}
-  @test Krylov.ktypeof(b) == Vector{Float32}
-
-  a = rand(Float64, 10)
-  b = view(a, 4:8)
-  @test Krylov.ktypeof(a) == Vector{Float64}
-  @test Krylov.ktypeof(b) == Vector{Float64}
-
-  a = sprand(Float32, 10, 0.5)
-  b = view(a, 4:8)
-  @test Krylov.ktypeof(a) == Vector{Float32}
-  @test Krylov.ktypeof(b) == Vector{Float32}
-
-  a = sprand(Float64, 10, 0.5)
-  b = view(a, 4:8)
-  @test Krylov.ktypeof(a) == Vector{Float64}
-  @test Krylov.ktypeof(b) == Vector{Float64}
+
+  @testset "sym_givens" begin
+    # test Givens reflector corner cases
+    (c, s, ρ) = Krylov.sym_givens(0.0, 0.0)
+    @test (c == 1.0) && (s == 0.0) && (ρ == 0.0)
+
+    a = 3.14
+    (c, s, ρ) = Krylov.sym_givens(a, 0.0)
+    @test (c == 1.0) && (s == 0.0) && (ρ == a)
+    (c, s, ρ) = Krylov.sym_givens(-a, 0.0)
+    @test (c == -1.0) && (s == 0.0) && (ρ == a)
+
+    b = 3.14
+    (c, s, ρ) = Krylov.sym_givens(0.0, b)
+    @test (c == 0.0) && (s == 1.0) && (ρ == b)
+    (c, s, ρ) = Krylov.sym_givens(0.0, -b)
+    @test (c == 0.0) && (s == -1.0) && (ρ == b)
+
+    (c, s, ρ) = Krylov.sym_givens(Complex(0.0), Complex(0.0))
+    @test (c == 1.0) && (s == Complex(0.0)) && (ρ == Complex(0.0))
+
+    a = Complex(1.0, 1.0)
+    (c, s, ρ) = Krylov.sym_givens(a, Complex(0.0))
+    @test (c == 1.0) && (s == Complex(0.0)) && (ρ == a)
+    (c, s, ρ) = Krylov.sym_givens(-a, Complex(0.0))
+    @test (c == 1.0) && (s == Complex(0.0)) && (ρ == -a)
+
+    b = Complex(1.0, 1.0)
+    (c, s, ρ) = Krylov.sym_givens(Complex(0.0), b)
+    @test (c == 0.0) && (s == Complex(1.0)) && (ρ == b)
+    (c, s, ρ) = Krylov.sym_givens(Complex(0.0), -b)
+    @test (c == 0.0) && (s == Complex(1.0)) && (ρ == -b)
+  end
+
+  @testset "roots_quadratic" begin
+    # test roots of a quadratic
+    roots = Krylov.roots_quadratic(0.0, 0.0, 0.0)
+    @test length(roots) == 1
+    @test roots[1] == 0.0
+
+    roots = Krylov.roots_quadratic(0.0, 0.0, 1.0)
+    @test length(roots) == 0
+
+    roots = Krylov.roots_quadratic(0.0, 3.14, -1.0)
+    @test length(roots) == 1
+    @test roots[1] == 1.0 / 3.14
+
+    roots = Krylov.roots_quadratic(1.0, 0.0, 1.0)
+    @test length(roots) == 0
+
+    roots = Krylov.roots_quadratic(1.0, 0.0, 0.0)
+    @test length(roots) == 2
+    @test roots[1] == 0.0
+    @test roots[2] == 0.0
+
+    roots = Krylov.roots_quadratic(1.0, 3.0, 2.0)
+    @test length(roots) == 2
+    @test roots[1] ≈ -2.0
+    @test roots[2] ≈ -1.0
+
+    roots = Krylov.roots_quadratic(1.0e+8, 1.0, 1.0)
+    @test length(roots) == 0
+
+    # ill-conditioned quadratic
+    roots = Krylov.roots_quadratic(-1.0e-8, 1.0e+5, 1.0, nitref=0)
+    @test length(roots) == 2
+    @test roots[1] == 1.0e+13
+    @test roots[2] == 0.0
+
+    # iterative refinement is crucial!
+    roots = Krylov.roots_quadratic(-1.0e-8, 1.0e+5, 1.0, nitref=1)
+    @test length(roots) == 2
+    @test roots[1] == 1.0e+13
+    @test roots[2] == -1.0e-05
+
+    # not ill-conditioned quadratic
+    roots = Krylov.roots_quadratic(-1.0e-7, 1.0, 1.0, nitref=0)
+    @test length(roots) == 2
+    @test isapprox(roots[1],  1.0e+7, rtol=1.0e-6)
+    @test isapprox(roots[2], -1.0, rtol=1.0e-6)
+
+    roots = Krylov.roots_quadratic(-1.0e-7, 1.0, 1.0, nitref=1)
+    @test length(roots) == 2
+    @test isapprox(roots[1], 1.0e+7, rtol=1.0e-6)
+    @test isapprox(roots[2], -1.0, rtol=1.0e-6)
+  end
+
+  @testset "to_boundary" begin
+    # test trust-region boundary
+    x = ones(5)
+    d = ones(5); d[1:2:5] .= -1
+    @test_throws ErrorException Krylov.to_boundary(x, d, -1.0)
+    @test_throws ErrorException Krylov.to_boundary(x, d, 0.5)
+    @test_throws ErrorException Krylov.to_boundary(x, zeros(5), 1.0)
+    @test maximum(Krylov.to_boundary(x, d, 5.0)) ≈ 2.209975124224178
+    @test minimum(Krylov.to_boundary(x, d, 5.0)) ≈ -1.8099751242241782
+    @test maximum(Krylov.to_boundary(x, d, 5.0, flip=true)) ≈ 1.8099751242241782
+    @test minimum(Krylov.to_boundary(x, d, 5.0, flip=true)) ≈ -2.209975124224178
+  end
+
+  @testset "kzeros" begin
+    # test kzeros
+    @test Krylov.kzeros(Vector{Float64}, 10) == zeros(Float64, 10)
+    @test Krylov.kzeros(Vector{ComplexF32}, 10) == zeros(ComplexF32, 10)
+  end
+
+  @testset "kones" begin
+    # test kones
+    @test Krylov.kones(Vector{Float64}, 10) == ones(Float64, 10)
+    @test Krylov.kones(Vector{ComplexF32}, 10) == ones(ComplexF32, 10)
+  end
+
+  @testset "ktypeof" begin
+    # test ktypeof
+    a = rand(Float32, 10)
+    b = view(a, 4:8)
+    @test Krylov.ktypeof(a) == Vector{Float32}
+    @test Krylov.ktypeof(b) == Vector{Float32}
+
+    a = rand(Float64, 10)
+    b = view(a, 4:8)
+    @test Krylov.ktypeof(a) == Vector{Float64}
+    @test Krylov.ktypeof(b) == Vector{Float64}
+
+    a = sprand(Float32, 10, 0.5)
+    b = view(a, 4:8)
+    @test Krylov.ktypeof(a) == Vector{Float32}
+    @test Krylov.ktypeof(b) == Vector{Float32}
+
+    a = sprand(Float64, 10, 0.5)
+    b = view(a, 4:8)
+    @test Krylov.ktypeof(a) == Vector{Float64}
+    @test Krylov.ktypeof(b) == Vector{Float64}
+  end
+
+  @testset "macros" begin
+    # test macros
+    for FC ∈ (Float16, Float32, Float64, Complex{Float16}, Complex{Float32}, Complex{Float64})
+      n = 10
+      x = rand(FC, n)
+      y = rand(FC, n)
+      a = rand(FC)
+      b = rand(FC)
+      c = rand(FC)
+      s = rand(FC)
+
+      T = real(FC)
+      a2 = rand(T)
+      b2 = rand(T)
+
+      Krylov.@kdot(n, x, y)
+
+      Krylov.@kdotr(n, x, y)
+
+      Krylov.@knrm2(n, x)
+
+      Krylov.@kaxpy!(n, a, x, y)
+      Krylov.@kaxpy!(n, a2, x, y)
+
+      Krylov.@kaxpby!(n, a, x, b, y)
+      Krylov.@kaxpby!(n, a2, x, b, y)
+      Krylov.@kaxpby!(n, a, x, b2, y)
+      Krylov.@kaxpby!(n, a2, x, b2, y)
+
+      Krylov.@kcopy!(n, x, y)
+
+      Krylov.@kswap(x, y)
+
+      Krylov.@kref!(n, x, y, c, s)
+    end
+  end
 end