diff --git a/src/bicgstab.jl b/src/bicgstab.jl
index c3b914599..3e5635775 100644
--- a/src/bicgstab.jl
+++ b/src/bicgstab.jl
@@ -26,10 +26,10 @@ export bicgstab, bicgstab!
 
 Solve the square linear system Ax = b using the BICGSTAB method.
 BICGSTAB requires two initial vectors `b` and `c`.
-The relation `bᵀc ≠ 0` must be satisfied and by default `c = b`.
+The relation `bᴴc ≠ 0` must be satisfied and by default `c = b`.
 
 The Biconjugate Gradient Stabilized method is a variant of BiCG, like CGS,
-but using different updates for the Aᵀ-sequence in order to obtain smoother
+but using different updates for the Aᴴ-sequence in order to obtain smoother
 convergence than CGS.
 
 If BICGSTAB stagnates, we recommend DQGMRES and BiLQ as alternative methods for unsymmetric square systems.
@@ -157,7 +157,7 @@ function bicgstab!(solver :: BicgstabSolver{T,FC,S}, A, b :: AbstractVector{FC};
   if next_ρ == 0
     stats.niter = 0
     stats.solved, stats.inconsistent = false, false
-    stats.status = "Breakdown bᵀc = 0"
+    stats.status = "Breakdown bᴴc = 0"
     solver.warm_start = false
     return solver
   end
diff --git a/src/bilq.jl b/src/bilq.jl
index ce84d3ec1..f40538245 100644
--- a/src/bilq.jl
+++ b/src/bilq.jl
@@ -24,7 +24,7 @@ export bilq, bilq!
 Solve the square linear system Ax = b using the BiLQ method.
 
 BiLQ is based on the Lanczos biorthogonalization process and requires two initial vectors `b` and `c`.
-The relation `bᵀc ≠ 0` must be satisfied and by default `c = b`.
+The relation `bᴴc ≠ 0` must be satisfied and by default `c = b`.
 When `A` is symmetric and `b = c`, BiLQ is equivalent to SYMMLQ.
 
 An option gives the possibility of transferring to the BiCG point,
@@ -90,7 +90,7 @@ function bilq!(solver :: BilqSolver{T,FC,S}, A, b :: AbstractVector{FC}; c :: Ab
   ktypeof(c) == S || error("ktypeof(c) ≠ $S")
 
   # Compute the adjoint of A
-  Aᵀ = A'
+  Aᴴ = A'
 
   # Set up workspace.
   uₖ₋₁, uₖ, q, vₖ₋₁, vₖ = solver.uₖ₋₁, solver.uₖ, solver.q, solver.vₖ₋₁, solver.vₖ
@@ -127,25 +127,25 @@ function bilq!(solver :: BilqSolver{T,FC,S}, A, b :: AbstractVector{FC}; c :: Ab
   kdisplay(iter, verbose) && @printf("%5d  %7.1e\n", iter, bNorm)
 
   # Initialize the Lanczos biorthogonalization process.
-  cᵗb = @kdot(n, c, r₀)  # ⟨c,r₀⟩
-  if cᵗb == 0
+  cᴴb = @kdot(n, c, r₀)  # ⟨c,r₀⟩
+  if cᴴb == 0
     stats.niter = 0
     stats.solved = false
     stats.inconsistent = false
-    stats.status = "Breakdown bᵀc = 0"
+    stats.status = "Breakdown bᴴc = 0"
     solver.warm_start = false
     return solver
   end
 
-  βₖ = √(abs(cᵗb))            # β₁γ₁ = cᵀ(b - Ax₀)
-  γₖ = cᵗb / βₖ               # β₁γ₁ = cᵀ(b - Ax₀)
+  βₖ = √(abs(cᴴb))            # β₁γ₁ = cᴴ(b - Ax₀)
+  γₖ = cᴴb / βₖ               # β₁γ₁ = cᴴ(b - Ax₀)
   vₖ₋₁ .= zero(FC)            # v₀ = 0
   uₖ₋₁ .= zero(FC)            # u₀ = 0
   vₖ .= r₀ ./ βₖ              # v₁ = (b - Ax₀) / β₁
   uₖ .= c ./ conj(γₖ)         # u₁ = c / γ̄₁
   cₖ₋₁ = cₖ = -one(T)         # Givens cosines used for the LQ factorization of Tₖ
   sₖ₋₁ = sₖ = zero(FC)        # Givens sines used for the LQ factorization of Tₖ
-  d̅ .= zero(FC)               # Last column of D̅ₖ = Vₖ(Qₖ)ᵀ
+  d̅ .= zero(FC)               # Last column of D̅ₖ = Vₖ(Qₖ)ᴴ
   ζₖ₋₁ = ζbarₖ = zero(FC)     # ζₖ₋₁ and ζbarₖ are the last components of z̅ₖ = (L̅ₖ)⁻¹β₁e₁
   ζₖ₋₂ = ηₖ = zero(FC)        # ζₖ₋₂ and ηₖ are used to update ζₖ₋₁ and ζbarₖ
   δbarₖ₋₁ = δbarₖ = zero(FC)  # Coefficients of Lₖ₋₁ and L̅ₖ modified over the course of two iterations
@@ -165,10 +165,10 @@ function bilq!(solver :: BilqSolver{T,FC,S}, A, b :: AbstractVector{FC}; c :: Ab
 
     # Continue the Lanczos biorthogonalization process.
     # AVₖ  = VₖTₖ    + βₖ₊₁vₖ₊₁(eₖ)ᵀ = Vₖ₊₁Tₖ₊₁.ₖ
-    # AᵀUₖ = Uₖ(Tₖ)ᵀ + γ̄ₖ₊₁uₖ₊₁(eₖ)ᵀ = Uₖ₊₁(Tₖ.ₖ₊₁)ᵀ
+    # AᴴUₖ = Uₖ(Tₖ)ᴴ + γ̄ₖ₊₁uₖ₊₁(eₖ)ᵀ = Uₖ₊₁(Tₖ.ₖ₊₁)ᴴ
 
     mul!(q, A , vₖ)  # Forms vₖ₊₁ : q ← Avₖ
-    mul!(p, Aᵀ, uₖ)  # Forms uₖ₊₁ : p ← Aᵀuₖ
+    mul!(p, Aᴴ, uₖ)  # Forms uₖ₊₁ : p ← Aᴴuₖ
 
     @kaxpy!(n, -γₖ, vₖ₋₁, q)  # q ← q - γₖ * vₖ₋₁
     @kaxpy!(n, -βₖ, uₖ₋₁, p)  # p ← p - β̄ₖ * uₖ₋₁
@@ -178,9 +178,9 @@ function bilq!(solver :: BilqSolver{T,FC,S}, A, b :: AbstractVector{FC}; c :: Ab
     @kaxpy!(n, -     αₖ , vₖ, q)    # q ← q - αₖ * vₖ
     @kaxpy!(n, -conj(αₖ), uₖ, p)    # p ← p - ᾱₖ * uₖ
 
-    pᵗq = @kdot(n, p, q)      # pᵗq  = ⟨p,q⟩
-    βₖ₊₁ = √(abs(pᵗq))        # βₖ₊₁ = √(|pᵗq|)
-    γₖ₊₁ = pᵗq / βₖ₊₁         # γₖ₊₁ = pᵗq / βₖ₊₁
+    pᴴq = @kdot(n, p, q)      # pᴴq  = ⟨p,q⟩
+    βₖ₊₁ = √(abs(pᴴq))        # βₖ₊₁ = √(|pᴴq|)
+    γₖ₊₁ = pᴴq / βₖ₊₁         # γₖ₊₁ = pᴴq / βₖ₊₁
 
     # Update the LQ factorization of Tₖ = L̅ₖQₖ.
     # [ α₁ γ₂ 0  •  •  •  0 ]   [ δ₁   0    •   •   •    •    0   ]
@@ -235,7 +235,7 @@ function bilq!(solver :: BilqSolver{T,FC,S}, A, b :: AbstractVector{FC}; c :: Ab
       ηₖ   = -ϵₖ₋₂ * ζₖ₋₂ - λₖ₋₁ * ζₖ₋₁
     end
 
-    # Relations for the directions dₖ₋₁ and d̅ₖ, the last two columns of D̅ₖ = Vₖ(Qₖ)ᵀ.
+    # Relations for the directions dₖ₋₁ and d̅ₖ, the last two columns of D̅ₖ = Vₖ(Qₖ)ᴴ.
     # [d̅ₖ₋₁ vₖ] [cₖ  s̄ₖ] = [dₖ₋₁ d̅ₖ] ⟷ dₖ₋₁ = cₖ * d̅ₖ₋₁ + sₖ * vₖ
     #           [sₖ -cₖ]             ⟷ d̅ₖ   = s̄ₖ * d̅ₖ₋₁ - cₖ * vₖ
     if iter ≥ 2
@@ -258,13 +258,13 @@ function bilq!(solver :: BilqSolver{T,FC,S}, A, b :: AbstractVector{FC}; c :: Ab
     @. vₖ₋₁ = vₖ # vₖ₋₁ ← vₖ
     @. uₖ₋₁ = uₖ # uₖ₋₁ ← uₖ
 
-    if pᵗq ≠ 0
+    if pᴴq ≠ 0
       @. vₖ = q / βₖ₊₁        # βₖ₊₁vₖ₊₁ = q
       @. uₖ = p / conj(γₖ₊₁)  # γ̄ₖ₊₁uₖ₊₁ = p
     end
 
     # Compute ⟨vₖ,vₖ₊₁⟩ and ‖vₖ₊₁‖
-    vₖᵀvₖ₊₁ = @kdot(n, vₖ₋₁, vₖ)
+    vₖᴴvₖ₊₁ = @kdot(n, vₖ₋₁, vₖ)
     norm_vₖ₊₁ = @knrm2(n, vₖ)
 
     # Compute BiLQ residual norm
@@ -274,7 +274,7 @@ function bilq!(solver :: BilqSolver{T,FC,S}, A, b :: AbstractVector{FC}; c :: Ab
     else
       μₖ = βₖ * (sₖ₋₁ * ζₖ₋₂ - cₖ₋₁ * cₖ * ζₖ₋₁) + αₖ * sₖ * ζₖ₋₁
       ωₖ = βₖ₊₁ * sₖ * ζₖ₋₁
-      θₖ = conj(μₖ) * ωₖ * vₖᵀvₖ₊₁
+      θₖ = conj(μₖ) * ωₖ * vₖᴴvₖ₊₁
       rNorm_lq = sqrt(abs2(μₖ) * norm_vₖ^2 + abs2(ωₖ) * norm_vₖ₊₁^2 + 2 * real(θₖ))
     end
     history && push!(rNorms, rNorm_lq)
@@ -300,7 +300,7 @@ function bilq!(solver :: BilqSolver{T,FC,S}, A, b :: AbstractVector{FC}; c :: Ab
     solved_lq = rNorm_lq ≤ ε
     solved_cg = transfer_to_bicg && (abs(δbarₖ) > eps(T)) && (rNorm_cg ≤ ε)
     tired = iter ≥ itmax
-    breakdown = !solved_lq && !solved_cg && (pᵗq == 0)
+    breakdown = !solved_lq && !solved_cg && (pᴴq == 0)
     kdisplay(iter, verbose) && @printf("%5d  %7.1e\n", iter, rNorm_lq)
   end
   (verbose > 0) && @printf("\n")
diff --git a/src/bilqr.jl b/src/bilqr.jl
index 09fef1f6c..7284597dc 100644
--- a/src/bilqr.jl
+++ b/src/bilqr.jl
@@ -1,5 +1,5 @@
 # An implementation of BILQR for the solution of square
-# consistent linear adjoint systems Ax = b and Aᵀy = c.
+# consistent linear adjoint systems Ax = b and Aᴴy = c.
 #
 # This method is described in
 #
@@ -24,11 +24,11 @@ export bilqr, bilqr!
 Combine BiLQ and QMR to solve adjoint systems.
 
     [0  A] [y] = [b]
-    [Aᵀ 0] [x]   [c]
+    [Aᴴ 0] [x]   [c]
 
-The relation `bᵀc ≠ 0` must be satisfied.
+The relation `bᴴc ≠ 0` must be satisfied.
 BiLQ is used for solving primal system `Ax = b`.
-QMR is used for solving dual system `Aᵀy = c`.
+QMR is used for solving dual system `Aᴴy = c`.
 
 An option gives the possibility of transferring from the BiLQ point to the
 BiCG point, when it exists. The transfer is based on the residual norm.
@@ -94,7 +94,7 @@ function bilqr!(solver :: BilqrSolver{T,FC,S}, A, b :: AbstractVector{FC}, c ::
   ktypeof(c) == S || error("ktypeof(c) ≠ $S")
 
   # Compute the adjoint of A
-  Aᵀ = A'
+  Aᴴ = A'
 
   # Set up workspace.
   uₖ₋₁, uₖ, q, vₖ₋₁, vₖ = solver.uₖ₋₁, solver.uₖ, solver.q, solver.vₖ₋₁, solver.vₖ
@@ -109,7 +109,7 @@ function bilqr!(solver :: BilqrSolver{T,FC,S}, A, b :: AbstractVector{FC}, c ::
   if warm_start
     mul!(r₀, A, Δx)
     @kaxpby!(n, one(FC), b, -one(FC), r₀)
-    mul!(s₀, Aᵀ, Δy)
+    mul!(s₀, Aᴴ, Δy)
     @kaxpby!(n, one(FC), c, -one(FC), s₀)
   end
 
@@ -117,7 +117,7 @@ function bilqr!(solver :: BilqrSolver{T,FC,S}, A, b :: AbstractVector{FC}, c ::
   x .= zero(FC)          # x₀
   bNorm = @knrm2(n, r₀)  # rNorm = ‖r₀‖
 
-  # Initial solution t₀ and residual norm ‖s₀‖ = ‖c - Aᵀy₀‖.
+  # Initial solution t₀ and residual norm ‖s₀‖ = ‖c - Aᴴy₀‖.
   t .= zero(FC)          # t₀
   cNorm = @knrm2(n, s₀)  # sNorm = ‖s₀‖
 
@@ -132,34 +132,34 @@ function bilqr!(solver :: BilqrSolver{T,FC,S}, A, b :: AbstractVector{FC}, c ::
   kdisplay(iter, verbose) && @printf("%5d  %7.1e  %7.1e\n", iter, bNorm, cNorm)
 
   # Initialize the Lanczos biorthogonalization process.
-  cᵗb = @kdot(n, s₀, r₀)  # ⟨s₀,r₀⟩ = ⟨c - Aᵀy₀,b - Ax₀⟩
-  if cᵗb == 0
+  cᴴb = @kdot(n, s₀, r₀)  # ⟨s₀,r₀⟩ = ⟨c - Aᴴy₀,b - Ax₀⟩
+  if cᴴb == 0
     stats.niter = 0
     stats.solved_primal = false
     stats.solved_dual = false
-    stats.status = "Breakdown bᵀc = 0"
+    stats.status = "Breakdown bᴴc = 0"
     solver.warm_start = false
     return solver
   end
 
   # Set up workspace.
-  βₖ = √(abs(cᵗb))            # β₁γ₁ = (c - Aᵀy₀)ᵀ(b - Ax₀)
-  γₖ = cᵗb / βₖ               # β₁γ₁ = (c - Aᵀy₀)ᵀ(b - Ax₀)
+  βₖ = √(abs(cᴴb))            # β₁γ₁ = (c - Aᴴy₀)ᴴ(b - Ax₀)
+  γₖ = cᴴb / βₖ               # β₁γ₁ = (c - Aᴴy₀)ᴴ(b - Ax₀)
   vₖ₋₁ .= zero(FC)            # v₀ = 0
   uₖ₋₁ .= zero(FC)            # u₀ = 0
   vₖ .= r₀ ./ βₖ              # v₁ = (b - Ax₀) / β₁
-  uₖ .= s₀ ./ conj(γₖ)        # u₁ = (c - Aᵀy₀) / γ̄₁
+  uₖ .= s₀ ./ conj(γₖ)        # u₁ = (c - Aᴴy₀) / γ̄₁
   cₖ₋₁ = cₖ = -one(T)         # Givens cosines used for the LQ factorization of Tₖ
   sₖ₋₁ = sₖ = zero(FC)        # Givens sines used for the LQ factorization of Tₖ
-  d̅ .= zero(FC)               # Last column of D̅ₖ = Vₖ(Qₖ)ᵀ
+  d̅ .= zero(FC)               # Last column of D̅ₖ = Vₖ(Qₖ)ᴴ
   ζₖ₋₁ = ζbarₖ = zero(FC)     # ζₖ₋₁ and ζbarₖ are the last components of z̅ₖ = (L̅ₖ)⁻¹β₁e₁
   ζₖ₋₂ = ηₖ = zero(FC)        # ζₖ₋₂ and ηₖ are used to update ζₖ₋₁ and ζbarₖ
   δbarₖ₋₁ = δbarₖ = zero(FC)  # Coefficients of Lₖ₋₁ and L̅ₖ modified over the course of two iterations
   ψbarₖ₋₁ = ψₖ₋₁ = zero(FC)   # ψₖ₋₁ and ψbarₖ are the last components of h̅ₖ = Qₖγ̄₁e₁
   norm_vₖ = bNorm / βₖ        # ‖vₖ‖ is used for residual norm estimates
   ϵₖ₋₃ = λₖ₋₂ = zero(FC)      # Components of Lₖ₋₁
-  wₖ₋₃ .= zero(FC)            # Column k-3 of Wₖ = Uₖ(Lₖ)⁻ᵀ
-  wₖ₋₂ .= zero(FC)            # Column k-2 of Wₖ = Uₖ(Lₖ)⁻ᵀ
+  wₖ₋₃ .= zero(FC)            # Column k-3 of Wₖ = Uₖ(Lₖ)⁻ᴴ
+  wₖ₋₂ .= zero(FC)            # Column k-2 of Wₖ = Uₖ(Lₖ)⁻ᴴ
   τₖ = zero(T)                # τₖ is used for the dual residual norm estimate
 
   # Stopping criterion.
@@ -180,10 +180,10 @@ function bilqr!(solver :: BilqrSolver{T,FC,S}, A, b :: AbstractVector{FC}, c ::
 
     # Continue the Lanczos biorthogonalization process.
     # AVₖ  = VₖTₖ    + βₖ₊₁vₖ₊₁(eₖ)ᵀ = Vₖ₊₁Tₖ₊₁.ₖ
-    # AᵀUₖ = Uₖ(Tₖ)ᵀ + γ̄ₖ₊₁uₖ₊₁(eₖ)ᵀ = Uₖ₊₁(Tₖ.ₖ₊₁)ᵀ
+    # AᴴUₖ = Uₖ(Tₖ)ᴴ + γ̄ₖ₊₁uₖ₊₁(eₖ)ᵀ = Uₖ₊₁(Tₖ.ₖ₊₁)ᴴ
 
     mul!(q, A , vₖ)  # Forms vₖ₊₁ : q ← Avₖ
-    mul!(p, Aᵀ, uₖ)  # Forms uₖ₊₁ : p ← Aᵀuₖ
+    mul!(p, Aᴴ, uₖ)  # Forms uₖ₊₁ : p ← Aᴴuₖ
 
     @kaxpy!(n, -γₖ, vₖ₋₁, q)  # q ← q - γₖ * vₖ₋₁
     @kaxpy!(n, -βₖ, uₖ₋₁, p)  # p ← p - β̄ₖ * uₖ₋₁
@@ -193,9 +193,9 @@ function bilqr!(solver :: BilqrSolver{T,FC,S}, A, b :: AbstractVector{FC}, c ::
     @kaxpy!(n, -     αₖ , vₖ, q)  # q ← q - αₖ * vₖ
     @kaxpy!(n, -conj(αₖ), uₖ, p)  # p ← p - ᾱₖ * uₖ
 
-    pᵗq = @kdot(n, p, q)  # pᵗq  = ⟨p,q⟩
-    βₖ₊₁ = √(abs(pᵗq))    # βₖ₊₁ = √(|pᵗq|)
-    γₖ₊₁ = pᵗq / βₖ₊₁     # γₖ₊₁ = pᵗq / βₖ₊₁
+    pᴴq = @kdot(n, p, q)  # pᴴq  = ⟨p,q⟩
+    βₖ₊₁ = √(abs(pᴴq))    # βₖ₊₁ = √(|pᴴq|)
+    γₖ₊₁ = pᴴq / βₖ₊₁     # γₖ₊₁ = pᴴq / βₖ₊₁
 
     # Update the LQ factorization of Tₖ = L̅ₖQₖ.
     # [ α₁ γ₂ 0  •  •  •  0 ]   [ δ₁   0    •   •   •    •    0   ]
@@ -251,7 +251,7 @@ function bilqr!(solver :: BilqrSolver{T,FC,S}, A, b :: AbstractVector{FC}, c ::
         ηₖ   = -ϵₖ₋₂ * ζₖ₋₂ - λₖ₋₁ * ζₖ₋₁
       end
 
-      # Relations for the directions dₖ₋₁ and d̅ₖ, the last two columns of D̅ₖ = Vₖ(Qₖ)ᵀ.
+      # Relations for the directions dₖ₋₁ and d̅ₖ, the last two columns of D̅ₖ = Vₖ(Qₖ)ᴴ.
       # [d̅ₖ₋₁ vₖ] [cₖ  s̄ₖ] = [dₖ₋₁ d̅ₖ] ⟷ dₖ₋₁ = cₖ * d̅ₖ₋₁ + sₖ * vₖ
       #           [sₖ -cₖ]             ⟷ d̅ₖ   = s̄ₖ * d̅ₖ₋₁ - cₖ * vₖ
       if iter ≥ 2
@@ -271,7 +271,7 @@ function bilqr!(solver :: BilqrSolver{T,FC,S}, A, b :: AbstractVector{FC}, c ::
       end
 
       # Compute ⟨vₖ,vₖ₊₁⟩ and ‖vₖ₊₁‖
-      vₖᵀvₖ₊₁ = @kdot(n, vₖ, q) / βₖ₊₁
+      vₖᴴvₖ₊₁ = @kdot(n, vₖ, q) / βₖ₊₁
       norm_vₖ₊₁ = @knrm2(n, q) / βₖ₊₁
 
       # Compute BiLQ residual norm
@@ -281,7 +281,7 @@ function bilqr!(solver :: BilqrSolver{T,FC,S}, A, b :: AbstractVector{FC}, c ::
       else
         μₖ = βₖ * (sₖ₋₁ * ζₖ₋₂ - cₖ₋₁ * cₖ * ζₖ₋₁) + αₖ * sₖ * ζₖ₋₁
         ωₖ = βₖ₊₁ * sₖ * ζₖ₋₁
-        θₖ = conj(μₖ) * ωₖ * vₖᵀvₖ₊₁
+        θₖ = conj(μₖ) * ωₖ * vₖᴴvₖ₊₁
         rNorm_lq = sqrt(abs2(μₖ) * norm_vₖ^2 + abs2(ωₖ) * norm_vₖ₊₁^2 + 2 * real(θₖ))
       end
       history && push!(rNorms, rNorm_lq)
@@ -318,7 +318,7 @@ function bilqr!(solver :: BilqrSolver{T,FC,S}, A, b :: AbstractVector{FC}, c ::
         ψbarₖ = sₖ * ψbarₖ₋₁
       end
 
-      # Compute the direction wₖ₋₁, the last column of Wₖ₋₁ = (Uₖ₋₁)(Lₖ₋₁)⁻ᵀ ⟷ (L̄ₖ₋₁)(Wₖ₋₁)ᵀ = (Uₖ₋₁)ᵀ.
+      # Compute the direction wₖ₋₁, the last column of Wₖ₋₁ = (Uₖ₋₁)(Lₖ₋₁)⁻ᴴ ⟷ (L̄ₖ₋₁)(Wₖ₋₁)ᵀ = (Uₖ₋₁)ᵀ.
       # w₁ = u₁ / δ̄₁
       if iter == 2
         wₖ₋₁ = wₖ₋₂
@@ -372,7 +372,7 @@ function bilqr!(solver :: BilqrSolver{T,FC,S}, A, b :: AbstractVector{FC}, c ::
     @. vₖ₋₁ = vₖ  # vₖ₋₁ ← vₖ
     @. uₖ₋₁ = uₖ  # uₖ₋₁ ← uₖ
 
-    if pᵗq ≠ zero(FC)
+    if pᴴq ≠ zero(FC)
       @. vₖ = q / βₖ₊₁        # βₖ₊₁vₖ₊₁ = q
       @. uₖ = p / conj(γₖ₊₁)  # γ̄ₖ₊₁uₖ₊₁ = p
     end
@@ -392,7 +392,7 @@ function bilqr!(solver :: BilqrSolver{T,FC,S}, A, b :: AbstractVector{FC}, c ::
 
     user_requested_exit = callback(solver) :: Bool
     tired = iter ≥ itmax
-    breakdown = !solved_lq && !solved_cg && (pᵗq == 0)
+    breakdown = !solved_lq && !solved_cg && (pᴴq == 0)
 
     kdisplay(iter, verbose) &&  solved_primal && !solved_dual && @printf("%5d  %7s  %7.1e\n", iter, "", sNorm)
     kdisplay(iter, verbose) && !solved_primal &&  solved_dual && @printf("%5d  %7.1e  %7s\n", iter, rNorm_lq, "")
diff --git a/src/cg_lanczos.jl b/src/cg_lanczos.jl
index a8e24f02f..2f2dae16d 100644
--- a/src/cg_lanczos.jl
+++ b/src/cg_lanczos.jl
@@ -111,7 +111,7 @@ function cg_lanczos!(solver :: CgLanczosSolver{T,FC,S}, A, b :: AbstractVector{F
     Mv .= b
   end
   MisI || mulorldiv!(v, M, Mv, ldiv)  # v₁ = M⁻¹r₀
-  β = sqrt(@kdotr(n, v, Mv))          # β₁ = v₁ᵀ M v₁
+  β = sqrt(@kdotr(n, v, Mv))          # β₁ = v₁ᴴ M v₁
   σ = β
   rNorm = σ
   history && push!(rNorms, rNorm)
@@ -157,10 +157,10 @@ function cg_lanczos!(solver :: CgLanczosSolver{T,FC,S}, A, b :: AbstractVector{F
     # Form next Lanczos vector.
     # βₖ₊₁Mvₖ₊₁ = Avₖ - δₖMvₖ - βₖMvₖ₋₁
     mul!(Mv_next, A, v)        # Mvₖ₊₁ ← Avₖ
-    δ = @kdotr(n, v, Mv_next)  # δₖ = vₖᵀ A vₖ
+    δ = @kdotr(n, v, Mv_next)  # δₖ = vₖᴴ A vₖ
 
     # Check curvature. Exit fast if requested.
-    # It is possible to show that σₖ² (δₖ - ωₖ₋₁ / γₖ₋₁) = pₖᵀ A pₖ.
+    # It is possible to show that σₖ² (δₖ - ωₖ₋₁ / γₖ₋₁) = pₖᴴ A pₖ.
     γ = one(T) / (δ - ω / γ)  # γₖ = 1 / (δₖ - ωₖ₋₁ / γₖ₋₁)
     indefinite |= (γ ≤ 0)
     (check_curvature & indefinite) && continue
@@ -172,7 +172,7 @@ function cg_lanczos!(solver :: CgLanczosSolver{T,FC,S}, A, b :: AbstractVector{F
     end
     @. Mv = Mv_next                      # Mvₖ ← Mvₖ₊₁
     MisI || mulorldiv!(v, M, Mv, ldiv)   # vₖ₊₁ = M⁻¹ * Mvₖ₊₁
-    β = sqrt(@kdotr(n, v, Mv))           # βₖ₊₁ = vₖ₊₁ᵀ M vₖ₊₁
+    β = sqrt(@kdotr(n, v, Mv))           # βₖ₊₁ = vₖ₊₁ᴴ M vₖ₊₁
     @kscal!(n, one(FC) / β, v)           # vₖ₊₁  ←  vₖ₊₁ / βₖ₊₁
     MisI || @kscal!(n, one(FC) / β, Mv)  # Mvₖ₊₁ ← Mvₖ₊₁ / βₖ₊₁
     Anorm2 += β_prev^2 + β^2 + δ^2       # Use ‖Tₖ₊₁‖₂ as increasing approximation of ‖A‖₂.
diff --git a/src/cg_lanczos_shift.jl b/src/cg_lanczos_shift.jl
index 01f11e41f..ff873e5b4 100644
--- a/src/cg_lanczos_shift.jl
+++ b/src/cg_lanczos_shift.jl
@@ -92,7 +92,7 @@ function cg_lanczos_shift!(solver :: CgLanczosShiftSolver{T,FC,S}, A, b :: Abstr
   end
   Mv .= b                             # Mv₁ ← b
   MisI || mulorldiv!(v, M, Mv, ldiv)  # v₁ = M⁻¹ * Mv₁
-  β = sqrt(@kdotr(n, v, Mv))          # β₁ = v₁ᵀ M v₁
+  β = sqrt(@kdotr(n, v, Mv))          # β₁ = v₁ᴴ M v₁
   rNorms .= β
   if history
     for i = 1 : nshifts
@@ -157,7 +157,7 @@ function cg_lanczos_shift!(solver :: CgLanczosShiftSolver{T,FC,S}, A, b :: Abstr
     # Form next Lanczos vector.
     # βₖ₊₁Mvₖ₊₁ = Avₖ - δₖMvₖ - βₖMvₖ₋₁
     mul!(Mv_next, A, v)                  # Mvₖ₊₁ ← Avₖ
-    δ = @kdotr(n, v, Mv_next)            # δₖ = vₖᵀ A vₖ
+    δ = @kdotr(n, v, Mv_next)            # δₖ = vₖᴴ A vₖ
     @kaxpy!(n, -δ, Mv, Mv_next)          # Mvₖ₊₁ ← Mvₖ₊₁ - δₖMvₖ
     if iter > 0
       @kaxpy!(n, -β, Mv_prev, Mv_next)   # Mvₖ₊₁ ← Mvₖ₊₁ - βₖMvₖ₋₁
@@ -165,12 +165,12 @@ function cg_lanczos_shift!(solver :: CgLanczosShiftSolver{T,FC,S}, A, b :: Abstr
     end
     @. Mv = Mv_next                      # Mvₖ ← Mvₖ₊₁
     MisI || mulorldiv!(v, M, Mv, ldiv)   # vₖ₊₁ = M⁻¹ * Mvₖ₊₁
-    β = sqrt(@kdotr(n, v, Mv))           # βₖ₊₁ = vₖ₊₁ᵀ M vₖ₊₁
+    β = sqrt(@kdotr(n, v, Mv))           # βₖ₊₁ = vₖ₊₁ᴴ M vₖ₊₁
     @kscal!(n, one(FC) / β, v)           # vₖ₊₁  ←  vₖ₊₁ / βₖ₊₁
     MisI || @kscal!(n, one(FC) / β, Mv)  # Mvₖ₊₁ ← Mvₖ₊₁ / βₖ₊₁
 
-    # Check curvature: vₖᵀ(A + sᵢI)vₖ = vₖᵀAvₖ + sᵢ‖vₖ‖² = δₖ + ρₖ * sᵢ with ρₖ = ‖vₖ‖².
-    # It is possible to show that σₖ² (δₖ + ρₖ * sᵢ - ωₖ₋₁ / γₖ₋₁) = pₖᵀ (A + sᵢ I) pₖ.
+    # Check curvature: vₖᴴ(A + sᵢI)vₖ = vₖᴴAvₖ + sᵢ‖vₖ‖² = δₖ + ρₖ * sᵢ with ρₖ = ‖vₖ‖².
+    # It is possible to show that σₖ² (δₖ + ρₖ * sᵢ - ωₖ₋₁ / γₖ₋₁) = pₖᴴ (A + sᵢ I) pₖ.
     MisI || (ρ = @kdotr(n, v, v))
     for i = 1 : nshifts
       δhat[i] = δ + ρ * shifts[i]
diff --git a/src/cgls.jl b/src/cgls.jl
index f5529fbfb..43fa5a6b6 100644
--- a/src/cgls.jl
+++ b/src/cgls.jl
@@ -5,7 +5,7 @@
 #
 # equivalently, of the normal equations
 #
-#  AᵀAx = Aᵀb.
+#  AᴴAx = Aᴴb.
 #
 # CGLS is formally equivalent to applying the conjugate gradient method
 # to the normal equations but should be more stable. It is also formally
@@ -45,11 +45,11 @@ Solve the regularized linear least-squares problem
 using the Conjugate Gradient (CG) method, where λ ≥ 0 is a regularization
 parameter. This method is equivalent to applying CG to the normal equations
 
-    (AᵀA + λI) x = Aᵀb
+    (AᴴA + λI) x = Aᴴb
 
 but is more stable.
 
-CGLS produces monotonic residuals ‖r‖₂ but not optimality residuals ‖Aᵀr‖₂.
+CGLS produces monotonic residuals ‖r‖₂ but not optimality residuals ‖Aᴴr‖₂.
 It is formally equivalent to LSQR, though can be slightly less accurate,
 but simpler to implement.
 
@@ -95,7 +95,7 @@ function cgls!(solver :: CglsSolver{T,FC,S}, A, b :: AbstractVector{FC};
   ktypeof(b) == S || error("ktypeof(b) ≠ $S")
 
   # Compute the adjoint of A
-  Aᵀ = A'
+  Aᴴ = A'
 
   # Set up workspace.
   allocate_if(!MisI, solver, :Mr, S, m)
@@ -117,9 +117,9 @@ function cgls!(solver :: CglsSolver{T,FC,S}, A, b :: AbstractVector{FC};
     return solver
   end
   MisI || mulorldiv!(Mr, M, r, ldiv)
-  mul!(s, Aᵀ, Mr)
+  mul!(s, Aᴴ, Mr)
   p .= s
-  γ = @kdotr(n, s, s)  # γ = sᵀs
+  γ = @kdotr(n, s, s)  # γ = sᴴs
   iter = 0
   itmax == 0 && (itmax = m + n)
 
@@ -128,7 +128,7 @@ function cgls!(solver :: CglsSolver{T,FC,S}, A, b :: AbstractVector{FC};
   history && push!(rNorms, rNorm)
   history && push!(ArNorms, ArNorm)
   ε = atol + rtol * ArNorm
-  (verbose > 0) && @printf("%5s  %8s  %8s\n", "k", "‖Aᵀr‖", "‖r‖")
+  (verbose > 0) && @printf("%5s  %8s  %8s\n", "k", "‖Aᴴr‖", "‖r‖")
   kdisplay(iter, verbose) && @printf("%5d  %8.2e  %8.2e\n", iter, ArNorm, rNorm)
 
   status = "unknown"
@@ -140,8 +140,8 @@ function cgls!(solver :: CglsSolver{T,FC,S}, A, b :: AbstractVector{FC};
   while ! (solved || tired || user_requested_exit)
     mul!(q, A, p)
     MisI || mulorldiv!(Mq, M, q, ldiv)
-    δ = @kdotr(m, q, Mq)  # δ = qᵀMq
-    λ > 0 && (δ += λ * @kdotr(n, p, p))  # δ = δ + pᵀp
+    δ = @kdotr(m, q, Mq)  # δ = qᴴMq
+    λ > 0 && (δ += λ * @kdotr(n, p, p))  # δ = δ + pᴴp
     α = γ / δ
 
     # if a trust-region constraint is give, compute step to the boundary
@@ -154,9 +154,9 @@ function cgls!(solver :: CglsSolver{T,FC,S}, A, b :: AbstractVector{FC};
     @kaxpy!(n,  α, p, x)     # Faster than x = x + α * p
     @kaxpy!(m, -α, q, r)     # Faster than r = r - α * q
     MisI || mulorldiv!(Mr, M, r, ldiv)
-    mul!(s, Aᵀ, Mr)
+    mul!(s, Aᴴ, Mr)
     λ > 0 && @kaxpy!(n, -λ, x, s)   # s = A' * r - λ * x
-    γ_next = @kdotr(n, s, s)   # γ_next = sᵀs
+    γ_next = @kdotr(n, s, s)   # γ_next = sᴴs
     β = γ_next / γ
     @kaxpby!(n, one(FC), s, β, p) # p = s + βp
     γ = γ_next
diff --git a/src/cgne.jl b/src/cgne.jl
index 2f720b57c..68039d2de 100644
--- a/src/cgne.jl
+++ b/src/cgne.jl
@@ -10,7 +10,7 @@
 # and is equivalent to applying the conjugate gradient method
 # to the linear system
 #
-#  AAᵀy = b.
+#  AAᴴy = b.
 #
 # This method is also known as Craig's method, CGME, and other
 # names, and is described in
@@ -46,7 +46,7 @@ using the Conjugate Gradient (CG) method, where λ ≥ 0 is a regularization
 parameter. This method is equivalent to applying CG to the normal equations
 of the second kind
 
-    (AAᵀ + λI) y = b
+    (AAᴴ + λI) y = b
 
 but is more stable. When λ = 0, this method solves the minimum-norm problem
 
@@ -104,12 +104,12 @@ function cgne!(solver :: CgneSolver{T,FC,S}, A, b :: AbstractVector{FC};
   ktypeof(b) == S || error("ktypeof(b) ≠ $S")
 
   # Compute the adjoint of A
-  Aᵀ = A'
+  Aᴴ = A'
 
   # Set up workspace.
   allocate_if(!NisI, solver, :z, S, m)
   allocate_if(λ > 0, solver, :s, S, m)
-  x, p, Aᵀz, r, q, s, stats = solver.x, solver.p, solver.Aᵀz, solver.r, solver.q, solver.s, solver.stats
+  x, p, Aᴴz, r, q, s, stats = solver.x, solver.p, solver.Aᴴz, solver.r, solver.q, solver.s, solver.stats
   rNorms = stats.residuals
   reset!(stats)
   z = NisI ? r : solver.z
@@ -126,7 +126,7 @@ function cgne!(solver :: CgneSolver{T,FC,S}, A, b :: AbstractVector{FC};
     return solver
   end
   λ > 0 && (s .= r)
-  mul!(p, Aᵀ, z)
+  mul!(p, Aᴴ, z)
 
   # Use ‖p‖ to detect inconsistent system.
   # An inconsistent system will necessarily have AA' singular.
@@ -161,8 +161,8 @@ function cgne!(solver :: CgneSolver{T,FC,S}, A, b :: AbstractVector{FC};
     NisI || mulorldiv!(z, N, r, ldiv)
     γ_next = @kdotr(m, r, z)  # Faster than γ_next = dot(r, z)
     β = γ_next / γ
-    mul!(Aᵀz, Aᵀ, z)
-    @kaxpby!(n, one(FC), Aᵀz, β, p)  # Faster than p = Aᵀz + β * p
+    mul!(Aᴴz, Aᴴ, z)
+    @kaxpby!(n, one(FC), Aᴴz, β, p)  # Faster than p = Aᴴz + β * p
     pNorm = @knrm2(n, p)
     if λ > 0
       @kaxpby!(m, one(FC), r, β, s)  # s = r + β * s
diff --git a/src/cgs.jl b/src/cgs.jl
index c1eb1056e..592eb1b2d 100644
--- a/src/cgs.jl
+++ b/src/cgs.jl
@@ -21,7 +21,7 @@ export cgs, cgs!
 
 Solve the consistent linear system Ax = b using conjugate gradient squared algorithm.
 CGS requires two initial vectors `b` and `c`.
-The relation `bᵀc ≠ 0` must be satisfied and by default `c = b`.
+The relation `bᴴc ≠ 0` must be satisfied and by default `c = b`.
 
 From "Iterative Methods for Sparse Linear Systems (Y. Saad)" :
 
@@ -142,7 +142,7 @@ function cgs!(solver :: CgsSolver{T,FC,S}, A, b :: AbstractVector{FC}; c :: Abst
   if ρ == 0
     stats.niter = 0
     stats.solved, stats.inconsistent = false, false
-    stats.status = "Breakdown bᵀc = 0"
+    stats.status = "Breakdown bᴴc = 0"
     solver.warm_start =false
     return solver
   end
diff --git a/src/cr.jl b/src/cr.jl
index c678c7d29..4405eda76 100644
--- a/src/cr.jl
+++ b/src/cr.jl
@@ -149,7 +149,7 @@ function cr!(solver :: CrSolver{T,FC,S}, A, b :: AbstractVector{FC};
   (verbose > 0) && @printf("%5s %8s %8s %8s\n", "k", "‖x‖", "‖r‖", "quad")
   kdisplay(iter, verbose) && @printf("    %d  %8.1e %8.1e %8.1e\n", iter, xNorm, rNorm, m)
 
-  descent = pr > 0 # pᵀr > 0 means p is a descent direction
+  descent = pr > 0 # pᴴr > 0 means p is a descent direction
   solved = rNorm ≤ ε
   tired = iter ≥ itmax
   on_boundary = false
@@ -161,7 +161,7 @@ function cr!(solver :: CrSolver{T,FC,S}, A, b :: AbstractVector{FC};
     if linesearch
       if (pAp ≤ γ * pNorm²) || (ρ ≤ γ * rNorm²)
         npcurv = true
-        (verbose > 0) && @printf("nonpositive curvature detected: pᵀAp = %8.1e and rᵀAr = %8.1e\n", pAp, ρ)
+        (verbose > 0) && @printf("nonpositive curvature detected: pᴴAp = %8.1e and rᴴAr = %8.1e\n", pAp, ρ)
         stats.solved = solved
         stats.inconsistent = false
         stats.status = "nonpositive curvature"
@@ -182,16 +182,16 @@ function cr!(solver :: CrSolver{T,FC,S}, A, b :: AbstractVector{FC};
       tr = maximum(to_boundary(x, r, radius; flip = false, xNorm2 = xNorm², dNorm2 = rNorm²))
       (verbose > 0) && @printf("t1 = %8.1e, t2 = %8.1e and tr = %8.1e\n", t1, t2, tr)
 
-      if abspAp ≤ γ * pNorm * @knrm2(n, q) # pᵀAp ≃ 0
+      if abspAp ≤ γ * pNorm * @knrm2(n, q) # pᴴAp ≃ 0
         npcurv = true # nonpositive curvature
-        (verbose > 0) && @printf("pᵀAp = %8.1e ≃ 0\n", pAp)
-        if abspr ≤ γ * pNorm * rNorm # pᵀr ≃ 0
-          (verbose > 0) && @printf("pᵀr = %8.1e ≃ 0, redefining p := r\n", pr)
+        (verbose > 0) && @printf("pᴴAp = %8.1e ≃ 0\n", pAp)
+        if abspr ≤ γ * pNorm * rNorm # pᴴr ≃ 0
+          (verbose > 0) && @printf("pᴴr = %8.1e ≃ 0, redefining p := r\n", pr)
           p = r # - ∇q(x)
           q = Ar
-          # q(x + αr) = q(x) - α ‖r‖² + ½ α² rᵀAr
-          # 1) if rᵀAr > 0, the quadratic decreases from α = 0 to α = ‖r‖² / rᵀAr
-          # 2) if rᵀAr ≤ 0, the quadratic decreases to -∞ in the direction r
+          # q(x + αr) = q(x) - α ‖r‖² + ½ α² rᴴAr
+          # 1) if rᴴAr > 0, the quadratic decreases from α = 0 to α = ‖r‖² / rᴴAr
+          # 2) if rᴴAr ≤ 0, the quadratic decreases to -∞ in the direction r
           if ρ > 0 # case 1
             (verbose > 0) && @printf("quadratic is convex in direction r, curv = %8.1e\n", ρ)
             α = min(tr, rNorm² / ρ)
@@ -200,12 +200,12 @@ function cr!(solver :: CrSolver{T,FC,S}, A, b :: AbstractVector{FC};
             α = tr
           end
         else
-          # q_p = q(x + α_p * p) - q(x) = -α_p * rᵀp + ½ (α_p)² * pᵀAp
-          # q_r = q(x + α_r * r) - q(x) = -α_r * ‖r‖² + ½ (α_r)² * rᵀAr
+          # q_p = q(x + α_p * p) - q(x) = -α_p * rᴴp + ½ (α_p)² * pᴴAp
+          # q_r = q(x + α_r * r) - q(x) = -α_r * ‖r‖² + ½ (α_r)² * rᴴAr
           # Δ = q_p - q_r. If Δ > 0, r is followed, else p is followed
           α = descent ? t1 : t2
           ρ > 0 && (tr = min(tr, rNorm² / ρ))
-          Δ = -α * pr + tr * rNorm² - (tr)^2 * ρ / 2 # as pᵀAp = 0
+          Δ = -α * pr + tr * rNorm² - (tr)^2 * ρ / 2 # as pᴴAp = 0
           if Δ > 0 # direction r engenders a better decrease
             (verbose > 0) && @printf("direction r engenders a bigger decrease. q_p - q_r = %8.1e > 0\n", Δ)
             (verbose > 0) && @printf("redefining p := r\n")
@@ -218,7 +218,7 @@ function cr!(solver :: CrSolver{T,FC,S}, A, b :: AbstractVector{FC};
         end
 
       elseif pAp > 0 && ρ > 0 # no negative curvature
-        (verbose > 0) && @printf("positive curvatures along p and r. pᵀAp = %8.1e and rᵀAr = %8.1e\n", pAp, ρ)
+        (verbose > 0) && @printf("positive curvatures along p and r. pᴴAp = %8.1e and rᴴAr = %8.1e\n", pAp, ρ)
         α = ρ / @kdotr(n, q, Mq)
         if α ≥ t1
           α = t1
@@ -227,8 +227,8 @@ function cr!(solver :: CrSolver{T,FC,S}, A, b :: AbstractVector{FC};
 
       elseif pAp > 0 && ρ < 0
         npcurv = true
-        (verbose > 0) && @printf("pᵀAp = %8.1e > 0 and rᵀAr = %8.1e < 0\n", pAp, ρ)
-        # q_p is minimal for α_p = rᵀp / pᵀAp
+        (verbose > 0) && @printf("pᴴAp = %8.1e > 0 and rᴴAr = %8.1e < 0\n", pAp, ρ)
+        # q_p is minimal for α_p = rᴴp / pᴴAp
         α = descent ?  min(t1, pr / pAp) : max(t2, pr / pAp)
         Δ = -α * pr + tr * rNorm² + (α^2 * pAp - (tr)^2 * ρ) / 2
         if Δ > 0
@@ -243,7 +243,7 @@ function cr!(solver :: CrSolver{T,FC,S}, A, b :: AbstractVector{FC};
 
       elseif pAp < 0 && ρ > 0
         npcurv = true
-        (verbose > 0) && @printf("pᵀAp = %8.1e < 0 and rᵀAr = %8.1e > 0\n", pAp, ρ)
+        (verbose > 0) && @printf("pᴴAp = %8.1e < 0 and rᴴAr = %8.1e > 0\n", pAp, ρ)
         α = descent ? t1 : t2
         tr = min(tr, rNorm² / ρ)
         Δ = -α * pr + tr * rNorm² + (α^2 * pAp - (tr)^2 * ρ) / 2
@@ -259,7 +259,7 @@ function cr!(solver :: CrSolver{T,FC,S}, A, b :: AbstractVector{FC};
 
       elseif pAp < 0 && ρ < 0
         npcurv = true
-        (verbose > 0) && @printf("negative curvatures along p and r. pᵀAp = %8.1e and rᵀAr = %8.1e\n", pAp, ρ)
+        (verbose > 0) && @printf("negative curvatures along p and r. pᴴAp = %8.1e and rᴴAr = %8.1e\n", pAp, ρ)
         α = descent ? t1 : t2
         Δ = -α * pr + tr * rNorm² + (α^2 * pAp - (tr)^2 * ρ) / 2
         if Δ > 0
@@ -330,9 +330,9 @@ function cr!(solver :: CrSolver{T,FC,S}, A, b :: AbstractVector{FC};
       solver.warm_start = false
       return solver
     end
-    pr = rNorm² + β * pr - β * α * pAp # pᵀr
+    pr = rNorm² + β * pr - β * α * pAp # pᴴr
     abspr = abs(pr)
-    pAp = ρ + β^2 * pAp # pᵀq
+    pAp = ρ + β^2 * pAp # pᴴq
     abspAp = abs(pAp)
     descent = pr > 0
 
diff --git a/src/craig.jl b/src/craig.jl
index 20597ea02..5759e31df 100644
--- a/src/craig.jl
+++ b/src/craig.jl
@@ -11,7 +11,7 @@
 # and is equivalent to applying the conjugate gradient method
 # to the linear system
 #
-#  AAᵀy = b.
+#  AAᴴy = b.
 #
 # This method, sometimes known under the name CRAIG, is the
 # Golub-Kahan implementation of CGNE, and is described in
@@ -52,14 +52,14 @@ regularization parameter. This method is equivalent to CGNE but is more
 stable.
 
 For a system in the form Ax = b, Craig's method is equivalent to applying
-CG to AAᵀy = b and recovering x = Aᵀy. Note that y are the Lagrange
+CG to AAᴴy = b and recovering x = Aᴴy. Note that y are the Lagrange
 multipliers of the least-norm problem
 
     minimize ‖x‖  s.t.  Ax = b.
 
 If `λ > 0`, CRAIG solves the symmetric and quasi-definite system
 
-    [ -F     Aᵀ ] [ x ]   [ 0 ]
+    [ -F     Aᴴ ] [ x ]   [ 0 ]
     [  A   λ²E  ] [ y ] = [ b ],
 
 where E and F are symmetric and positive definite.
@@ -70,12 +70,12 @@ The system above represents the optimality conditions of
 
     min ‖x‖²_F + λ²‖y‖²_E  s.t.  Ax + λ²Ey = b.
 
-For a symmetric and positive definite matrix `K`, the K-norm of a vector `x` is `‖x‖²_K = xᵀKx`.
-CRAIG is then equivalent to applying CG to `(AF⁻¹Aᵀ + λ²E)y = b` with `Fx = Aᵀy`.
+For a symmetric and positive definite matrix `K`, the K-norm of a vector `x` is `‖x‖²_K = xᴴKx`.
+CRAIG is then equivalent to applying CG to `(AF⁻¹Aᴴ + λ²E)y = b` with `Fx = Aᴴy`.
 
 If `λ = 0`, CRAIG solves the symmetric and indefinite system
 
-    [ -F   Aᵀ ] [ x ]   [ 0 ]
+    [ -F   Aᴴ ] [ x ]   [ 0 ]
     [  A   0  ] [ y ] = [ b ].
 
 The system above represents the optimality conditions of
@@ -134,13 +134,13 @@ function craig!(solver :: CraigSolver{T,FC,S}, A, b :: AbstractVector{FC};
   ktypeof(b) == S || error("ktypeof(b) ≠ $S")
 
   # Compute the adjoint of A
-  Aᵀ = A'
+  Aᴴ = A'
 
   # Set up workspace.
   allocate_if(!MisI, solver, :u , S, m)
   allocate_if(!NisI, solver, :v , S, n)
   allocate_if(λ > 0, solver, :w2, S, n)
-  x, Nv, Aᵀu, y, w = solver.x, solver.Nv, solver.Aᵀu, solver.y, solver.w
+  x, Nv, Aᴴu, y, w = solver.x, solver.Nv, solver.Aᴴu, solver.y, solver.w
   Mu, Av, w2, stats = solver.Mu, solver.Av, solver.w2, solver.stats
   rNorms = stats.residuals
   reset!(stats)
@@ -180,7 +180,7 @@ function craig!(solver :: CraigSolver{T,FC,S}, A, b :: AbstractVector{FC};
 
   Anorm² = zero(T) # Estimate of ‖A‖²_F.
   Anorm  = zero(T)
-  Dnorm² = zero(T) # Estimate of ‖(AᵀA)⁻¹‖².
+  Dnorm² = zero(T) # Estimate of ‖(AᴴA)⁻¹‖².
   Acond  = zero(T) # Estimate of cond(A).
   xNorm² = zero(T) # Estimate of ‖x‖².
   xNorm  = zero(T)
@@ -212,9 +212,9 @@ function craig!(solver :: CraigSolver{T,FC,S}, A, b :: AbstractVector{FC};
 
   while ! (solved || inconsistent || ill_cond || tired || user_requested_exit)
     # Generate the next Golub-Kahan vectors
-    # 1. αₖ₊₁Nvₖ₊₁ = Aᵀuₖ₊₁ - βₖ₊₁Nvₖ
-    mul!(Aᵀu, Aᵀ, u)
-    @kaxpby!(n, one(FC), Aᵀu, -β, Nv)
+    # 1. αₖ₊₁Nvₖ₊₁ = Aᴴuₖ₊₁ - βₖ₊₁Nvₖ
+    mul!(Aᴴu, Aᴴ, u)
+    @kaxpby!(n, one(FC), Aᴴu, -β, Nv)
     NisI || mulorldiv!(v, N, Nv, ldiv)
     α = sqrt(@kdotr(n, v, Nv))
     if α == 0
diff --git a/src/craigmr.jl b/src/craigmr.jl
index e08bb9c36..854e3df98 100644
--- a/src/craigmr.jl
+++ b/src/craigmr.jl
@@ -10,7 +10,7 @@
 # and is equivalent to applying the conjugate residual method
 # to the linear system
 #
-#  AAᵀy = b.
+#  AAᴴy = b.
 #
 # This method is equivalent to CRMR, and is described in
 #
@@ -44,7 +44,7 @@ using the CRAIGMR method, where λ ≥ 0 is a regularization parameter.
 This method is equivalent to applying the Conjugate Residuals method
 to the normal equations of the second kind
 
-    (AAᵀ + λ²I) y = b
+    (AAᴴ + λ²I) y = b
 
 but is more stable. When λ = 0, this method solves the minimum-norm problem
 
@@ -52,7 +52,7 @@ but is more stable. When λ = 0, this method solves the minimum-norm problem
 
 If `λ > 0`, CRAIGMR solves the symmetric and quasi-definite system
 
-    [ -F    Aᵀ ] [ x ]   [ 0 ]
+    [ -F    Aᴴ ] [ x ]   [ 0 ]
     [  A  λ²E  ] [ y ] = [ b ],
 
 where E and F are symmetric and positive definite.
@@ -63,12 +63,12 @@ The system above represents the optimality conditions of
 
     min ‖x‖²_F + λ²‖y‖²_E  s.t.  Ax + λ²Ey = b.
 
-For a symmetric and positive definite matrix `K`, the K-norm of a vector `x` is `‖x‖²_K = xᵀKx`.
-CRAIGMR is then equivalent to applying MINRES to `(AF⁻¹Aᵀ + λ²E)y = b` with `Fx = Aᵀy`.
+For a symmetric and positive definite matrix `K`, the K-norm of a vector `x` is `‖x‖²_K = xᴴKx`.
+CRAIGMR is then equivalent to applying MINRES to `(AF⁻¹Aᴴ + λ²E)y = b` with `Fx = Aᴴy`.
 
 If `λ = 0`, CRAIGMR solves the symmetric and indefinite system
 
-    [ -F   Aᵀ ] [ x ]   [ 0 ]
+    [ -F   Aᴴ ] [ x ]   [ 0 ]
     [  A   0  ] [ y ] = [ b ].
 
 The system above represents the optimality conditions of
@@ -129,20 +129,20 @@ function craigmr!(solver :: CraigmrSolver{T,FC,S}, A, b :: AbstractVector{FC};
   ktypeof(b) == S || error("ktypeof(b) ≠ $S")
 
   # Compute the adjoint of A
-  Aᵀ = A'
+  Aᴴ = A'
 
   # Set up workspace.
   allocate_if(!MisI, solver, :u, S, m)
   allocate_if(!NisI, solver, :v, S, n)
   allocate_if(λ > 0, solver, :q, S, n)
-  x, Nv, Aᵀu, d, y, Mu = solver.x, solver.Nv, solver.Aᵀu, solver.d, solver.y, solver.Mu
+  x, Nv, Aᴴu, d, y, Mu = solver.x, solver.Nv, solver.Aᴴu, solver.d, solver.y, solver.Mu
   w, wbar, Av, q, stats = solver.w, solver.wbar, solver.Av, solver.q, solver.stats
   rNorms, ArNorms = stats.residuals, stats.Aresiduals
   reset!(stats)
   u = MisI ? Mu : solver.u
   v = NisI ? Nv : solver.v
 
-  # Compute y such that AAᵀy = b. Then recover x = Aᵀy.
+  # Compute y such that AAᴴy = b. Then recover x = Aᴴy.
   x .= zero(FC)
   y .= zero(FC)
   Mu .= b
@@ -161,9 +161,9 @@ function craigmr!(solver :: CraigmrSolver{T,FC,S}, A, b :: AbstractVector{FC};
   # β₁Mu₁ = b.
   @kscal!(m, one(FC)/β, u)
   MisI || @kscal!(m, one(FC)/β, Mu)
-  # α₁Nv₁ = Aᵀu₁.
-  mul!(Aᵀu, Aᵀ, u)
-  Nv .= Aᵀu
+  # α₁Nv₁ = Aᴴu₁.
+  mul!(Aᴴu, Aᴴ, u)
+  Nv .= Aᴴu
   NisI || mulorldiv!(v, N, Nv, ldiv)
   α = sqrt(@kdotr(n, v, Nv))
   Anorm² = α * α
@@ -171,10 +171,10 @@ function craigmr!(solver :: CraigmrSolver{T,FC,S}, A, b :: AbstractVector{FC};
   iter = 0
   itmax == 0 && (itmax = m + n)
 
-  (verbose > 0) && @printf("%5s  %7s  %7s  %7s  %7s  %8s  %8s  %7s\n", "k", "‖r‖", "‖Aᵀr‖", "β", "α", "cos", "sin", "‖A‖²")
+  (verbose > 0) && @printf("%5s  %7s  %7s  %7s  %7s  %8s  %8s  %7s\n", "k", "‖r‖", "‖Aᴴr‖", "β", "α", "cos", "sin", "‖A‖²")
   kdisplay(iter, verbose) && @printf("%5d  %7.1e  %7.1e  %7.1e  %7.1e  %8.1e  %8.1e  %7.1e\n", iter, β, α, β, α, 0, 1, Anorm²)
 
-  # Aᵀb = 0 so x = 0 is a minimum least-squares solution
+  # Aᴴb = 0 so x = 0 is a minimum least-squares solution
   if α == 0
     stats.niter = 0
     stats.solved, stats.inconsistent = true, false
@@ -288,9 +288,9 @@ function craigmr!(solver :: CraigmrSolver{T,FC,S}, A, b :: AbstractVector{FC};
     # xₖ = Dₖzₖ
     @kaxpy!(n, ζ, d, x)
 
-    # 2. αₖ₊₁Nvₖ₊₁ = Aᵀuₖ₊₁ - βₖ₊₁Nvₖ
-    mul!(Aᵀu, Aᵀ, u)
-    @kaxpby!(n, one(FC), Aᵀu, -β, Nv)
+    # 2. αₖ₊₁Nvₖ₊₁ = Aᴴuₖ₊₁ - βₖ₊₁Nvₖ
+    mul!(Aᴴu, Aᴴ, u)
+    @kaxpby!(n, one(FC), Aᴴu, -β, Nv)
     NisI || mulorldiv!(v, N, Nv, ldiv)
     α = sqrt(@kdotr(n, v, Nv))
     Anorm² = Anorm² + α * α  # = ‖Lₖ‖
diff --git a/src/crls.jl b/src/crls.jl
index 6410fb836..b041f8e9f 100644
--- a/src/crls.jl
+++ b/src/crls.jl
@@ -5,7 +5,7 @@
 #
 # equivalently, of the linear system
 #
-#  AᵀAx = Aᵀb.
+#  AᴴAx = Aᴴb.
 #
 # This implementation follows the formulation given in
 #
@@ -37,11 +37,11 @@ Solve the linear least-squares problem
 using the Conjugate Residuals (CR) method. This method is equivalent to
 applying MINRES to the normal equations
 
-    (AᵀA + λI) x = Aᵀb.
+    (AᴴA + λI) x = Aᴴb.
 
 This implementation recurs the residual r := b - Ax.
 
-CRLS produces monotonic residuals ‖r‖₂ and optimality residuals ‖Aᵀr‖₂.
+CRLS produces monotonic residuals ‖r‖₂ and optimality residuals ‖Aᴴr‖₂.
 It is formally equivalent to LSMR, though can be substantially less accurate,
 but simpler to implement.
 
@@ -86,7 +86,7 @@ function crls!(solver :: CrlsSolver{T,FC,S}, A, b :: AbstractVector{FC};
   ktypeof(b) == S || error("ktypeof(b) ≠ $S")
 
   # Compute the adjoint of A
-  Aᵀ = A'
+  Aᴴ = A'
 
   # Set up workspace.
   allocate_if(!MisI, solver, :Ms, S, m)
@@ -112,13 +112,13 @@ function crls!(solver :: CrlsSolver{T,FC,S}, A, b :: AbstractVector{FC};
   end
 
   MisI || mulorldiv!(Mr, M, r, ldiv)
-  mul!(Ar, Aᵀ, Mr)  # - λ * x0 if x0 ≠ 0.
+  mul!(Ar, Aᴴ, Mr)  # - λ * x0 if x0 ≠ 0.
   mul!(s, A, Ar)
   MisI || mulorldiv!(Ms, M, s, ldiv)
 
   p  .= Ar
   Ap .= s
-  mul!(q, Aᵀ, Ms)  # Ap
+  mul!(q, Aᴴ, Ms)  # Ap
   λ > 0 && @kaxpy!(n, λ, p, q)  # q = q + λ * p
   γ  = @kdotr(m, s, Ms)  # Faster than γ = dot(s, Ms)
   iter = 0
@@ -128,7 +128,7 @@ function crls!(solver :: CrlsSolver{T,FC,S}, A, b :: AbstractVector{FC};
   λ > 0 && (γ += λ * ArNorm * ArNorm)
   history && push!(ArNorms, ArNorm)
   ε = atol + rtol * ArNorm
-  (verbose > 0) && @printf("%5s  %8s  %8s\n", "k", "‖Aᵀr‖", "‖r‖")
+  (verbose > 0) && @printf("%5s  %8s  %8s\n", "k", "‖Aᴴr‖", "‖r‖")
   kdisplay(iter, verbose) && @printf("%5d  %8.2e  %8.2e\n", iter, ArNorm, rNorm)
 
   status = "unknown"
@@ -147,11 +147,11 @@ function crls!(solver :: CrlsSolver{T,FC,S}, A, b :: AbstractVector{FC};
     if radius > 0
       pNorm = @knrm2(n, p)
       if @kdotr(m, Ap, Ap) ≤ ε * sqrt(qNorm²) * pNorm # the quadratic is constant in the direction p
-        psd = true # det(AᵀA) = 0
-        p = Ar # p = Aᵀr
+        psd = true # det(AᴴA) = 0
+        p = Ar # p = Aᴴr
         pNorm² = ArNorm * ArNorm
-        mul!(q, Aᵀ, s)
-        α = min(ArNorm^2 / γ, maximum(to_boundary(x, p, radius, flip = false, dNorm2 = pNorm²))) # the quadratic is minimal in the direction Aᵀr for α = ‖Ar‖²/γ
+        mul!(q, Aᴴ, s)
+        α = min(ArNorm^2 / γ, maximum(to_boundary(x, p, radius, flip = false, dNorm2 = pNorm²))) # the quadratic is minimal in the direction Aᴴr for α = ‖Ar‖²/γ
       else
         pNorm² = pNorm * pNorm
         σ = maximum(to_boundary(x, p, radius, flip = false, dNorm2 = pNorm²))
@@ -177,7 +177,7 @@ function crls!(solver :: CrlsSolver{T,FC,S}, A, b :: AbstractVector{FC};
     @kaxpby!(n, one(FC), Ar, β, p)    # Faster than  p = Ar + β *  p
     @kaxpby!(m, one(FC), s, β, Ap)    # Faster than Ap =  s + β * Ap
     MisI || mulorldiv!(MAp, M, Ap, ldiv)
-    mul!(q, Aᵀ, MAp)
+    mul!(q, Aᴴ, MAp)
     λ > 0 && @kaxpy!(n, λ, p, q)  # q = q + λ * p
 
     γ = γ_next
diff --git a/src/crmr.jl b/src/crmr.jl
index 6ed2b3c60..3fff12b08 100644
--- a/src/crmr.jl
+++ b/src/crmr.jl
@@ -10,7 +10,7 @@
 # and is equivalent to applying the conjugate residual method
 # to the linear system
 #
-#  AAᵀy = b.
+#  AAᴴy = b.
 #
 # This method is equivalent to CRAIGMR, described in
 #
@@ -44,7 +44,7 @@ using the Conjugate Residual (CR) method, where λ ≥ 0 is a regularization
 parameter. This method is equivalent to applying CR to the normal equations
 of the second kind
 
-    (AAᵀ + λI) y = b
+    (AAᴴ + λI) y = b
 
 but is more stable. When λ = 0, this method solves the minimum-norm problem
 
@@ -102,19 +102,19 @@ function crmr!(solver :: CrmrSolver{T,FC,S}, A, b :: AbstractVector{FC};
   ktypeof(b) == S || error("ktypeof(b) ≠ $S")
 
   # Compute the adjoint of A
-  Aᵀ = A'
+  Aᴴ = A'
 
   # Set up workspace.
   allocate_if(!NisI, solver, :Nq, S, m)
   allocate_if(λ > 0, solver, :s , S, m)
-  x, p, Aᵀr, r = solver.x, solver.p, solver.Aᵀr, solver.r
+  x, p, Aᴴr, r = solver.x, solver.p, solver.Aᴴr, solver.r
   q, s, stats = solver.q, solver.s, solver.stats
   rNorms, ArNorms = stats.residuals, stats.Aresiduals
   reset!(stats)
   Nq = NisI ? q : solver.Nq
 
   x .= zero(FC)              # initial estimation x = 0
-  mulorldiv!(r, N, b, ldiv)  # initial residual r = M * (b - Ax) = M * b
+  mulorldiv!(r, N, b, ldiv)  # initial residual r = N * (b - Ax) = N * b
   bNorm = @knrm2(m, r)       # norm(b - A * x0) if x0 ≠ 0.
   rNorm = bNorm              # + λ * ‖x0‖ if x0 ≠ 0 and λ > 0.
   history && push!(rNorms, rNorm)
@@ -126,9 +126,9 @@ function crmr!(solver :: CrmrSolver{T,FC,S}, A, b :: AbstractVector{FC};
     return solver
   end
   λ > 0 && (s .= r)
-  mul!(Aᵀr, Aᵀ, r)  # - λ * x0 if x0 ≠ 0.
-  p .= Aᵀr
-  γ = @kdotr(n, Aᵀr, Aᵀr)  # Faster than γ = dot(Aᵀr, Aᵀr)
+  mul!(Aᴴr, Aᴴ, r)  # - λ * x0 if x0 ≠ 0.
+  p .= Aᴴr
+  γ = @kdotr(n, Aᴴr, Aᴴr)  # Faster than γ = dot(Aᴴr, Aᴴr)
   λ > 0 && (γ += λ * rNorm * rNorm)
   iter = 0
   itmax == 0 && (itmax = m + n)
@@ -137,7 +137,7 @@ function crmr!(solver :: CrmrSolver{T,FC,S}, A, b :: AbstractVector{FC};
   history && push!(ArNorms, ArNorm)
   ɛ_c = atol + rtol * rNorm  # Stopping tolerance for consistent systems.
   ɛ_i = atol + rtol * ArNorm  # Stopping tolerance for inconsistent systems.
-  (verbose > 0) && @printf("%5s  %8s  %8s\n", "k", "‖Aᵀr‖", "‖r‖")
+  (verbose > 0) && @printf("%5s  %8s  %8s\n", "k", "‖Aᴴr‖", "‖r‖")
   kdisplay(iter, verbose) && @printf("%5d  %8.2e  %8.2e\n", iter, ArNorm, rNorm)
 
   status = "unknown"
@@ -150,16 +150,16 @@ function crmr!(solver :: CrmrSolver{T,FC,S}, A, b :: AbstractVector{FC};
     mul!(q, A, p)
     λ > 0 && @kaxpy!(m, λ, s, q)  # q = q + λ * s
     NisI || mulorldiv!(Nq, N, q, ldiv)
-    α = γ / @kdotr(m, q, Nq)   # Compute qᵗ * M * q
+    α = γ / @kdotr(m, q, Nq)   # Compute qᴴ * N * q
     @kaxpy!(n,  α, p, x)       # Faster than  x =  x + α *  p
     @kaxpy!(m, -α, Nq, r)      # Faster than  r =  r - α * Nq
     rNorm = @knrm2(m, r)       # norm(r)
-    mul!(Aᵀr, Aᵀ, r)
-    γ_next = @kdotr(n, Aᵀr, Aᵀr)  # Faster than γ_next = dot(Aᵀr, Aᵀr)
+    mul!(Aᴴr, Aᴴ, r)
+    γ_next = @kdotr(n, Aᴴr, Aᴴr)  # Faster than γ_next = dot(Aᴴr, Aᴴr)
     λ > 0 && (γ_next += λ * rNorm * rNorm)
     β = γ_next / γ
 
-    @kaxpby!(n, one(FC), Aᵀr, β, p)  # Faster than  p = Aᵀr + β * p
+    @kaxpby!(n, one(FC), Aᴴr, β, p)  # Faster than  p = Aᴴr + β * p
     if λ > 0
       @kaxpby!(m, one(FC), r, β, s) # s = r + β * s
     end
diff --git a/src/fom.jl b/src/fom.jl
index fcae5cf62..b212129ef 100644
--- a/src/fom.jl
+++ b/src/fom.jl
@@ -211,7 +211,7 @@ function fom!(solver :: FomSolver{T,FC,S}, A, b :: AbstractVector{FC};
       mul!(w, A, p)                                  # w ← AN⁻¹vₖ
       MisI || mulorldiv!(q, M, w, ldiv)              # q ← M⁻¹AN⁻¹vₖ
       for i = 1 : inner_iter
-        U[nr+i] = @kdot(n, V[i], q)      # hᵢₖ = qᵀvᵢ
+        U[nr+i] = @kdot(n, V[i], q)      # hᵢₖ = (vᵢ)ᴴq
         @kaxpy!(n, -U[nr+i], V[i], q)    # q ← q - hᵢₖvᵢ
       end
 
diff --git a/src/gmres.jl b/src/gmres.jl
index 388a4ab96..32999aa23 100644
--- a/src/gmres.jl
+++ b/src/gmres.jl
@@ -214,7 +214,7 @@ function gmres!(solver :: GmresSolver{T,FC,S}, A, b :: AbstractVector{FC};
       mul!(w, A, p)                                  # w ← AN⁻¹vₖ
       MisI || mulorldiv!(q, M, w, ldiv)              # q ← M⁻¹AN⁻¹vₖ
       for i = 1 : inner_iter
-        R[nr+i] = @kdot(n, V[i], q)      # hᵢₖ = qᵀvᵢ
+        R[nr+i] = @kdot(n, V[i], q)      # hᵢₖ = (vᵢ)ᴴq
         @kaxpy!(n, -R[nr+i], V[i], q)    # q ← q - hᵢₖvᵢ
       end
 
@@ -245,7 +245,7 @@ function gmres!(solver :: GmresSolver{T,FC,S}, A, b :: AbstractVector{FC};
       # [s̄ₖ -cₖ] [hₖ₊₁.ₖ]   [ 0  ]
       (c[inner_iter], s[inner_iter], R[nr+inner_iter]) = sym_givens(R[nr+inner_iter], Hbis)
 
-      # Update zₖ = (Qₖ)ᵀβe₁
+      # Update zₖ = (Qₖ)ᴴβe₁
       ζₖ₊₁          = conj(s[inner_iter]) * z[inner_iter]
       z[inner_iter] =      c[inner_iter]  * z[inner_iter]
 
diff --git a/src/gpmr.jl b/src/gpmr.jl
index b10942995..82499b50e 100644
--- a/src/gpmr.jl
+++ b/src/gpmr.jl
@@ -28,7 +28,7 @@ GPMR solves the unsymmetric partitioned linear system
     [  B  μI ] [ y ]   [ c ],
 
 where λ and μ are real or complex numbers.
-`A` can have any shape and `B` has the shape of `Aᵀ`.
+`A` can have any shape and `B` has the shape of `Aᴴ`.
 `A`, `B`, `b` and `c` must be all nonzero.
 
 This implementation allows left and right block diagonal preconditioners
@@ -172,7 +172,7 @@ function gpmr!(solver :: GpmrSolver{T,FC,S}, A, B, b :: AbstractVector{FC}, c ::
   gs .= zero(FC)  # Givens sines used for the factorization QₖRₖ = Sₖ₊₁.ₖ.
   gc .= zero(T)   # Givens cosines used for the factorization QₖRₖ = Sₖ₊₁.ₖ.
   R  .= zero(FC)  # Upper triangular matrix Rₖ.
-  zt .= zero(FC)  # Rₖzₖ = tₖ with (tₖ, τbar₂ₖ₊₁, τbar₂ₖ₊₂) = (Qₖ)ᵀ(βe₁ + γe₂).
+  zt .= zero(FC)  # Rₖzₖ = tₖ with (tₖ, τbar₂ₖ₊₁, τbar₂ₖ₊₂) = (Qₖ)ᴴ(βe₁ + γe₂).
 
   # Warm-start
   # If λ ≠ 0, Cb₀ = Cb - CAΔy - λΔx because CM = Iₘ and E = Iₘ
@@ -259,8 +259,8 @@ function gpmr!(solver :: GpmrSolver{T,FC,S}, A, B, b :: AbstractVector{FC}, c ::
     DisI || mulorldiv!(p, D, dB, ldiv)        # p  = DBEvₖ
 
     for i = 1 : iter
-      hᵢₖ = @kdot(m, V[i], q)    # hᵢ.ₖ = vᵢAuₖ
-      fᵢₖ = @kdot(n, U[i], p)    # fᵢ.ₖ = uᵢBvₖ
+      hᵢₖ = @kdot(m, V[i], q)    # hᵢ.ₖ = (vᵢ)ᴴq
+      fᵢₖ = @kdot(n, U[i], p)    # fᵢ.ₖ = (uᵢ)ᴴp
       @kaxpy!(m, -hᵢₖ, V[i], q)  # q ← q - hᵢ.ₖvᵢ
       @kaxpy!(n, -fᵢₖ, U[i], p)  # p ← p - fᵢ.ₖuᵢ
       R[nr₂ₖ + 2i-1] = hᵢₖ
@@ -270,8 +270,8 @@ function gpmr!(solver :: GpmrSolver{T,FC,S}, A, B, b :: AbstractVector{FC}, c ::
     # Reorthogonalization of the Krylov basis.
     if reorthogonalization
       for i = 1 : iter
-        Htmp = @kdot(m, V[i], q)    # hₜₘₚ = qᵀvᵢ
-        Ftmp = @kdot(n, U[i], p)    # fₜₘₚ = pᵀuᵢ
+        Htmp = @kdot(m, V[i], q)    # hₜₘₚ = (vᵢ)ᴴq
+        Ftmp = @kdot(n, U[i], p)    # fₜₘₚ = (uᵢ)ᴴp
         @kaxpy!(m, -Htmp, V[i], q)  # q ← q - hₜₘₚvᵢ
         @kaxpy!(n, -Ftmp, U[i], p)  # p ← p - fₜₘₚuᵢ
         R[nr₂ₖ + 2i-1] += Htmp                            # hᵢ.ₖ = hᵢ.ₖ + hₜₘₚ
diff --git a/src/krylov_solvers.jl b/src/krylov_solvers.jl
index d557d91ae..abd0c7352 100644
--- a/src/krylov_solvers.jl
+++ b/src/krylov_solvers.jl
@@ -1092,7 +1092,7 @@ may be used in order to create these vectors.
 mutable struct CgneSolver{T,FC,S} <: KrylovSolver{T,FC,S}
   x     :: S
   p     :: S
-  Aᵀz   :: S
+  Aᴴz   :: S
   r     :: S
   q     :: S
   s     :: S
@@ -1105,13 +1105,13 @@ function CgneSolver(n, m, S)
   T   = real(FC)
   x   = S(undef, m)
   p   = S(undef, m)
-  Aᵀz = S(undef, m)
+  Aᴴz = S(undef, m)
   r   = S(undef, n)
   q   = S(undef, n)
   s   = S(undef, 0)
   z   = S(undef, 0)
   stats = SimpleStats(0, false, false, T[], T[], T[], "unknown")
-  solver = CgneSolver{T,FC,S}(x, p, Aᵀz, r, q, s, z, stats)
+  solver = CgneSolver{T,FC,S}(x, p, Aᴴz, r, q, s, z, stats)
   return solver
 end
 
@@ -1134,7 +1134,7 @@ may be used in order to create these vectors.
 mutable struct CrmrSolver{T,FC,S} <: KrylovSolver{T,FC,S}
   x     :: S
   p     :: S
-  Aᵀr   :: S
+  Aᴴr   :: S
   r     :: S
   q     :: S
   Nq    :: S
@@ -1147,13 +1147,13 @@ function CrmrSolver(n, m, S)
   T   = real(FC)
   x   = S(undef, m)
   p   = S(undef, m)
-  Aᵀr = S(undef, m)
+  Aᴴr = S(undef, m)
   r   = S(undef, n)
   q   = S(undef, n)
   Nq  = S(undef, 0)
   s   = S(undef, 0)
   stats = SimpleStats(0, false, false, T[], T[], T[], "unknown")
-  solver = CrmrSolver{T,FC,S}(x, p, Aᵀr, r, q, Nq, s, stats)
+  solver = CrmrSolver{T,FC,S}(x, p, Aᴴr, r, q, Nq, s, stats)
   return solver
 end
 
@@ -1176,7 +1176,7 @@ may be used in order to create these vectors.
 mutable struct LslqSolver{T,FC,S} <: KrylovSolver{T,FC,S}
   x       :: S
   Nv      :: S
-  Aᵀu     :: S
+  Aᴴu     :: S
   w̄       :: S
   Mu      :: S
   Av      :: S
@@ -1191,7 +1191,7 @@ function LslqSolver(n, m, S; window :: Int=5)
   T   = real(FC)
   x   = S(undef, m)
   Nv  = S(undef, m)
-  Aᵀu = S(undef, m)
+  Aᴴu = S(undef, m)
   w̄   = S(undef, m)
   Mu  = S(undef, n)
   Av  = S(undef, n)
@@ -1199,7 +1199,7 @@ function LslqSolver(n, m, S; window :: Int=5)
   v   = S(undef, 0)
   err_vec = zeros(T, window)
   stats = LSLQStats(0, false, false, T[], T[], T[], false, T[], T[], "unknown")
-  solver = LslqSolver{T,FC,S}(x, Nv, Aᵀu, w̄, Mu, Av, u, v, err_vec, stats)
+  solver = LslqSolver{T,FC,S}(x, Nv, Aᴴu, w̄, Mu, Av, u, v, err_vec, stats)
   return solver
 end
 
@@ -1222,7 +1222,7 @@ may be used in order to create these vectors.
 mutable struct LsqrSolver{T,FC,S} <: KrylovSolver{T,FC,S}
   x       :: S
   Nv      :: S
-  Aᵀu     :: S
+  Aᴴu     :: S
   w       :: S
   Mu      :: S
   Av      :: S
@@ -1237,7 +1237,7 @@ function LsqrSolver(n, m, S; window :: Int=5)
   T   = real(FC)
   x   = S(undef, m)
   Nv  = S(undef, m)
-  Aᵀu = S(undef, m)
+  Aᴴu = S(undef, m)
   w   = S(undef, m)
   Mu  = S(undef, n)
   Av  = S(undef, n)
@@ -1245,7 +1245,7 @@ function LsqrSolver(n, m, S; window :: Int=5)
   v   = S(undef, 0)
   err_vec = zeros(T, window)
   stats = SimpleStats(0, false, false, T[], T[], T[], "unknown")
-  solver = LsqrSolver{T,FC,S}(x, Nv, Aᵀu, w, Mu, Av, u, v, err_vec, stats)
+  solver = LsqrSolver{T,FC,S}(x, Nv, Aᴴu, w, Mu, Av, u, v, err_vec, stats)
   return solver
 end
 
@@ -1268,7 +1268,7 @@ may be used in order to create these vectors.
 mutable struct LsmrSolver{T,FC,S} <: KrylovSolver{T,FC,S}
   x       :: S
   Nv      :: S
-  Aᵀu     :: S
+  Aᴴu     :: S
   h       :: S
   hbar    :: S
   Mu      :: S
@@ -1284,7 +1284,7 @@ function LsmrSolver(n, m, S; window :: Int=5)
   T    = real(FC)
   x    = S(undef, m)
   Nv   = S(undef, m)
-  Aᵀu  = S(undef, m)
+  Aᴴu  = S(undef, m)
   h    = S(undef, m)
   hbar = S(undef, m)
   Mu   = S(undef, n)
@@ -1293,7 +1293,7 @@ function LsmrSolver(n, m, S; window :: Int=5)
   v    = S(undef, 0)
   err_vec = zeros(T, window)
   stats = LsmrStats(0, false, false, T[], T[], zero(T), zero(T), zero(T), zero(T), zero(T), "unknown")
-  solver = LsmrSolver{T,FC,S}(x, Nv, Aᵀu, h, hbar, Mu, Av, u, v, err_vec, stats)
+  solver = LsmrSolver{T,FC,S}(x, Nv, Aᴴu, h, hbar, Mu, Av, u, v, err_vec, stats)
   return solver
 end
 
@@ -1316,7 +1316,7 @@ may be used in order to create these vectors.
 mutable struct LnlqSolver{T,FC,S} <: KrylovSolver{T,FC,S}
   x     :: S
   Nv    :: S
-  Aᵀu   :: S
+  Aᴴu   :: S
   y     :: S
   w̄     :: S
   Mu    :: S
@@ -1332,7 +1332,7 @@ function LnlqSolver(n, m, S)
   T   = real(FC)
   x   = S(undef, m)
   Nv  = S(undef, m)
-  Aᵀu = S(undef, m)
+  Aᴴu = S(undef, m)
   y   = S(undef, n)
   w̄   = S(undef, n)
   Mu  = S(undef, n)
@@ -1341,7 +1341,7 @@ function LnlqSolver(n, m, S)
   v   = S(undef, 0)
   q   = S(undef, 0)
   stats = LNLQStats(0, false, T[], false, T[], T[], "unknown")
-  solver = LnlqSolver{T,FC,S}(x, Nv, Aᵀu, y, w̄, Mu, Av, u, v, q, stats)
+  solver = LnlqSolver{T,FC,S}(x, Nv, Aᴴu, y, w̄, Mu, Av, u, v, q, stats)
   return solver
 end
 
@@ -1364,7 +1364,7 @@ may be used in order to create these vectors.
 mutable struct CraigSolver{T,FC,S} <: KrylovSolver{T,FC,S}
   x     :: S
   Nv    :: S
-  Aᵀu   :: S
+  Aᴴu   :: S
   y     :: S
   w     :: S
   Mu    :: S
@@ -1380,7 +1380,7 @@ function CraigSolver(n, m, S)
   T   = real(FC)
   x   = S(undef, m)
   Nv  = S(undef, m)
-  Aᵀu = S(undef, m)
+  Aᴴu = S(undef, m)
   y   = S(undef, n)
   w   = S(undef, n)
   Mu  = S(undef, n)
@@ -1389,7 +1389,7 @@ function CraigSolver(n, m, S)
   v   = S(undef, 0)
   w2  = S(undef, 0)
   stats = SimpleStats(0, false, false, T[], T[], T[], "unknown")
-  solver = CraigSolver{T,FC,S}(x, Nv, Aᵀu, y, w, Mu, Av, u, v, w2, stats)
+  solver = CraigSolver{T,FC,S}(x, Nv, Aᴴu, y, w, Mu, Av, u, v, w2, stats)
   return solver
 end
 
@@ -1412,7 +1412,7 @@ may be used in order to create these vectors.
 mutable struct CraigmrSolver{T,FC,S} <: KrylovSolver{T,FC,S}
   x     :: S
   Nv    :: S
-  Aᵀu   :: S
+  Aᴴu   :: S
   d     :: S
   y     :: S
   Mu    :: S
@@ -1430,7 +1430,7 @@ function CraigmrSolver(n, m, S)
   T    = real(FC)
   x    = S(undef, m)
   Nv   = S(undef, m)
-  Aᵀu  = S(undef, m)
+  Aᴴu  = S(undef, m)
   d    = S(undef, m)
   y    = S(undef, n)
   Mu   = S(undef, n)
@@ -1441,7 +1441,7 @@ function CraigmrSolver(n, m, S)
   v    = S(undef, 0)
   q    = S(undef, 0)
   stats = SimpleStats(0, false, false, T[], T[], T[], "unknown")
-  solver = CraigmrSolver{T,FC,S}(x, Nv, Aᵀu, d, y, Mu, w, wbar, Av, u, v, q, stats)
+  solver = CraigmrSolver{T,FC,S}(x, Nv, Aᴴu, d, y, Mu, w, wbar, Av, u, v, q, stats)
   return solver
 end
 
diff --git a/src/krylov_utils.jl b/src/krylov_utils.jl
index 6f0c1c382..c61bf2e5e 100644
--- a/src/krylov_utils.jl
+++ b/src/krylov_utils.jl
@@ -164,14 +164,14 @@ function to_boundary(x :: Vector{T}, d :: Vector{T},
                      radius :: T; flip :: Bool=false, xNorm2 :: T=zero(T), dNorm2 :: T=zero(T)) where T <: Number
   radius > 0 || error("radius must be positive")
 
-  # ‖d‖² σ² + 2 xᵀd σ + (‖x‖² - radius²).
-  xd = dot(x, d)
-  flip && (xd = -xd)
+  # ‖d‖² σ² + (xᴴd + dᴴx) σ + (‖x‖² - radius²).
+  rxd = real(dot(x, d))
+  flip && (rxd = -rxd)
   dNorm2 == zero(T) && (dNorm2 = dot(d, d))
   dNorm2 == zero(T) && error("zero direction")
   xNorm2 == zero(T) && (xNorm2 = dot(x, x))
   (xNorm2 ≤ radius * radius) || error(@sprintf("outside of the trust region: ‖x‖²=%7.1e, Δ²=%7.1e", xNorm2, radius * radius))
-  roots = roots_quadratic(dNorm2, 2 * xd, xNorm2 - radius * radius)
+  roots = roots_quadratic(dNorm2, 2 * rxd, xNorm2 - radius * radius)
   return roots # `σ1` and `σ2`
 end
 
diff --git a/src/lnlq.jl b/src/lnlq.jl
index a1f890de2..db0a7c951 100644
--- a/src/lnlq.jl
+++ b/src/lnlq.jl
@@ -9,9 +9,9 @@
 # and is equivalent to applying the SYMMLQ method
 # to the linear system
 #
-#  AAᵀy = b with x = Aᵀy and can be reformulated as
+#  AAᴴy = b with x = Aᴴy and can be reformulated as
 #
-#  [ -I  Aᵀ ][ x ] = [ 0 ]
+#  [ -I  Aᴴ ][ x ] = [ 0 ]
 #  [  A     ][ y ]   [ b ].
 #
 # This method is based on the Golub-Kahan bidiagonalization process and is described in
@@ -41,14 +41,14 @@ Find the least-norm solution of the consistent linear system
 using the LNLQ method, where λ ≥ 0 is a regularization parameter.
 
 For a system in the form Ax = b, LNLQ method is equivalent to applying
-SYMMLQ to AAᵀy = b and recovering x = Aᵀy but is more stable.
+SYMMLQ to AAᴴy = b and recovering x = Aᴴy but is more stable.
 Note that y are the Lagrange multipliers of the least-norm problem
 
     minimize ‖x‖  s.t.  Ax = b.
 
 If `λ > 0`, LNLQ solves the symmetric and quasi-definite system
 
-    [ -F    Aᵀ ] [ x ]   [ 0 ]
+    [ -F    Aᴴ ] [ x ]   [ 0 ]
     [  A  λ²E  ] [ y ] = [ b ],
 
 where E and F are symmetric and positive definite.
@@ -59,12 +59,12 @@ The system above represents the optimality conditions of
 
     min ‖x‖²_F + λ²‖y‖²_E  s.t.  Ax + λ²Ey = b.
 
-For a symmetric and positive definite matrix `K`, the K-norm of a vector `x` is `‖x‖²_K = xᵀKx`.
-LNLQ is then equivalent to applying SYMMLQ to `(AF⁻¹Aᵀ + λ²E)y = b` with `Fx = Aᵀy`.
+For a symmetric and positive definite matrix `K`, the K-norm of a vector `x` is `‖x‖²_K = xᴴKx`.
+LNLQ is then equivalent to applying SYMMLQ to `(AF⁻¹Aᴴ + λ²E)y = b` with `Fx = Aᴴy`.
 
 If `λ = 0`, LNLQ solves the symmetric and indefinite system
 
-    [ -F   Aᵀ ] [ x ]   [ 0 ]
+    [ -F   Aᴴ ] [ x ]   [ 0 ]
     [  A   0  ] [ y ] = [ b ].
 
 The system above represents the optimality conditions of
@@ -126,13 +126,13 @@ function lnlq!(solver :: LnlqSolver{T,FC,S}, A, b :: AbstractVector{FC};
   ktypeof(b) == S || error("ktypeof(b) ≠ $S")
 
   # Compute the adjoint of A
-  Aᵀ = A'
+  Aᴴ = A'
 
   # Set up workspace.
   allocate_if(!MisI, solver, :u, S, m)
   allocate_if(!NisI, solver, :v, S, n)
   allocate_if(λ > 0, solver, :q, S, n)
-  x, Nv, Aᵀu, y, w̄ = solver.x, solver.Nv, solver.Aᵀu, solver.y, solver.w̄
+  x, Nv, Aᴴu, y, w̄ = solver.x, solver.Nv, solver.Aᴴu, solver.y, solver.w̄
   Mu, Av, q, stats = solver.Mu, solver.Av, solver.q, solver.stats
   rNorms, xNorms, yNorms = stats.residuals, stats.error_bnd_x, stats.error_bnd_y
   reset!(stats)
@@ -179,9 +179,9 @@ function lnlq!(solver :: LnlqSolver{T,FC,S}, A, b :: AbstractVector{FC};
     MisI || @kscal!(m, one(FC) / βₖ, Mu)
   end
 
-  # α₁Nv₁ = Aᵀu₁.
-  mul!(Aᵀu, Aᵀ, u)
-  Nv .= Aᵀu
+  # α₁Nv₁ = Aᴴu₁.
+  mul!(Aᴴu, Aᴴ, u)
+  Nv .= Aᴴu
   NisI || mulorldiv!(v, N, Nv, ldiv)  # v₁ = N⁻¹ * Nv₁
   αₖ = sqrt(@kdotr(n, v, Nv))         # α₁ = ‖v₁‖_N
   if αₖ ≠ 0
@@ -190,8 +190,8 @@ function lnlq!(solver :: LnlqSolver{T,FC,S}, A, b :: AbstractVector{FC};
   end
 
   w̄ .= u           # Direction w̄₁
-  cₖ = zero(T)     # Givens cosines used for the LQ factorization of (Lₖ)ᵀ
-  sₖ = zero(FC)    # Givens sines used for the LQ factorization of (Lₖ)ᵀ
+  cₖ = zero(T)     # Givens cosines used for the LQ factorization of (Lₖ)ᴴ
+  sₖ = zero(FC)    # Givens sines used for the LQ factorization of (Lₖ)ᴴ
   ζₖ₋₁ = zero(FC)  # ζₖ₋₁ and ζbarₖ are the last components of z̅ₖ
   ηₖ = zero(FC)    # Coefficient of M̅ₖ
 
@@ -214,7 +214,7 @@ function lnlq!(solver :: LnlqSolver{T,FC,S}, A, b :: AbstractVector{FC};
     αhatₖ = αₖ
   end
 
-  # Begin the LQ factorization of (Lₖ)ᵀ = M̅ₖQₖ.
+  # Begin the LQ factorization of (Lₖ)ᴴ = M̅ₖQₖ.
   # [ α₁ β₂ 0  •  •  •  0 ]   [ ϵ₁  0   •   •   •   •   0   ]
   # [ 0  α₂ •  •        • ]   [ η₂  ϵ₂  •               •   ]
   # [ •  •  •  •  •     • ]   [ 0   •   •   •           •   ]
@@ -225,7 +225,7 @@ function lnlq!(solver :: LnlqSolver{T,FC,S}, A, b :: AbstractVector{FC};
 
   ϵbarₖ = αhatₖ  # ϵbar₁ = αhat₁
 
-  # Hₖ = Bₖ(Lₖ)ᵀ = [   Lₖ(Lₖ)ᵀ   ] ⟹ (Hₖ₋₁)ᵀ = [Lₖ₋₁Mₖ₋₁  0] Qₖ
+  # Hₖ = Bₖ(Lₖ)ᴴ = [   Lₖ(Lₖ)ᴴ   ] ⟹ (Hₖ₋₁)ᴴ = [Lₖ₋₁Mₖ₋₁  0] Qₖ
   #                [ αₖβₖ₊₁(eₖ)ᵀ ]
   #
   # Solve Lₖtₖ = β₁e₁ and M̅ₖz̅ₖ = tₖ
@@ -273,7 +273,7 @@ function lnlq!(solver :: LnlqSolver{T,FC,S}, A, b :: AbstractVector{FC};
 
     # Continue the generalized Golub-Kahan bidiagonalization.
     # AVₖ    = MUₖ₊₁Bₖ
-    # AᵀUₖ₊₁ = NVₖ(Bₖ)ᵀ + αₖ₊₁Nvₖ₊₁(eₖ₊₁)ᵀ = NVₖ₊₁(Lₖ₊₁)ᵀ
+    # AᴴUₖ₊₁ = NVₖ(Bₖ)ᴴ + αₖ₊₁Nvₖ₊₁(eₖ₊₁)ᴴ = NVₖ₊₁(Lₖ₊₁)ᴴ
     #
     #      [ α₁ 0  •  •  •  •  0 ]
     #      [ β₂ α₂ •           • ]
@@ -296,9 +296,9 @@ function lnlq!(solver :: LnlqSolver{T,FC,S}, A, b :: AbstractVector{FC};
       MisI || @kscal!(m, one(FC) / βₖ₊₁, Mu)
     end
 
-    # αₖ₊₁Nvₖ₊₁ = Aᵀuₖ₊₁ - βₖ₊₁Nvₖ
-    mul!(Aᵀu, Aᵀ, u)
-    @kaxpby!(n, one(FC), Aᵀu, -βₖ₊₁, Nv)
+    # αₖ₊₁Nvₖ₊₁ = Aᴴuₖ₊₁ - βₖ₊₁Nvₖ
+    mul!(Aᴴu, Aᴴ, u)
+    @kaxpby!(n, one(FC), Aᴴu, -βₖ₊₁, Nv)
     NisI || mulorldiv!(v, N, Nv, ldiv)  # vₖ₊₁ = N⁻¹ * Nvₖ₊₁
     αₖ₊₁ = sqrt(@kdotr(n, v, Nv))       # αₖ₊₁ = ‖vₖ₊₁‖_N
     if αₖ₊₁ ≠ 0
@@ -353,7 +353,7 @@ function lnlq!(solver :: LnlqSolver{T,FC,S}, A, b :: AbstractVector{FC};
       ρbar = ssig * μbar + csig * σₑₛₜ
     end
 
-    # Continue the LQ factorization of (Lₖ₊₁)ᵀ.
+    # Continue the LQ factorization of (Lₖ₊₁)ᴴ.
     # [ηₖ ϵbarₖ βₖ₊₁] [1     0     0 ] = [ηₖ  ϵₖ     0    ]
     # [0    0   αₖ₊₁] [0   cₖ₊₁  sₖ₊₁]   [0  ηₖ₊₁  ϵbarₖ₊₁]
     #                 [0   sₖ₊₁ -cₖ₊₁]
diff --git a/src/lslq.jl b/src/lslq.jl
index 908de19c5..d43d4a089 100644
--- a/src/lslq.jl
+++ b/src/lslq.jl
@@ -5,7 +5,7 @@
 #
 # equivalently, of the normal equations
 #
-#  AᵀAx = Aᵀb.
+#  AᴴAx = Aᴴb.
 #
 # LSLQ is formally equivalent to applying SYMMLQ to the normal equations
 # but should be more stable.
@@ -41,7 +41,7 @@ Solve the regularized linear least-squares problem
 using the LSLQ method, where λ ≥ 0 is a regularization parameter.
 LSLQ is formally equivalent to applying SYMMLQ to the normal equations
 
-    (AᵀA + λ²I) x = Aᵀb
+    (AᴴA + λ²I) x = Aᴴb
 
 but is more stable.
 
@@ -62,7 +62,7 @@ but is more stable.
 If `λ > 0`, we solve the symmetric and quasi-definite system
 
     [ E      A ] [ r ]   [ b ]
-    [ Aᵀ  -λ²F ] [ x ] = [ 0 ],
+    [ Aᴴ  -λ²F ] [ x ] = [ 0 ],
 
 where E and F are symmetric and positive definite.
 Preconditioners M = E⁻¹ ≻ 0 and N = F⁻¹ ≻ 0 may be provided in the form of linear operators.
@@ -72,19 +72,19 @@ The system above represents the optimality conditions of
 
     minimize ‖b - Ax‖²_E⁻¹ + λ²‖x‖²_F.
 
-For a symmetric and positive definite matrix `K`, the K-norm of a vector `x` is `‖x‖²_K = xᵀKx`.
-LSLQ is then equivalent to applying SYMMLQ to `(AᵀE⁻¹A + λ²F)x = AᵀE⁻¹b` with `r = E⁻¹(b - Ax)`.
+For a symmetric and positive definite matrix `K`, the K-norm of a vector `x` is `‖x‖²_K = xᴴKx`.
+LSLQ is then equivalent to applying SYMMLQ to `(AᴴE⁻¹A + λ²F)x = AᴴE⁻¹b` with `r = E⁻¹(b - Ax)`.
 
 If `λ = 0`, we solve the symmetric and indefinite system
 
     [ E    A ] [ r ]   [ b ]
-    [ Aᵀ   0 ] [ x ] = [ 0 ].
+    [ Aᴴ   0 ] [ x ] = [ 0 ].
 
 The system above represents the optimality conditions of
 
     minimize ‖b - Ax‖²_E⁻¹.
 
-In this case, `N` can still be specified and indicates the weighted norm in which `x` and `Aᵀr` should be measured.
+In this case, `N` can still be specified and indicates the weighted norm in which `x` and `Aᴴr` should be measured.
 `r` can be recovered by computing `E⁻¹(b - Ax)`.
 
 * `λ` is a regularization parameter (see the problem statement above)
@@ -116,8 +116,8 @@ In this case, `N` can still be specified and indicates the weighted norm in whic
 The iterations stop as soon as one of the following conditions holds true:
 
 * the optimality residual is sufficiently small (`stats.status = "found approximate minimum least-squares solution"`) in the sense that either
-  * ‖Aᵀr‖ / (‖A‖ ‖r‖) ≤ atol, or
-  * 1 + ‖Aᵀr‖ / (‖A‖ ‖r‖) ≤ 1
+  * ‖Aᴴr‖ / (‖A‖ ‖r‖) ≤ atol, or
+  * 1 + ‖Aᴴr‖ / (‖A‖ ‖r‖) ≤ 1
 * an approximate zero-residual solution has been found (`stats.status = "found approximate zero-residual solution"`) in the sense that either
   * ‖r‖ / ‖b‖ ≤ btol + atol ‖A‖ * ‖xᴸ‖ / ‖b‖, or
   * 1 + ‖r‖ / ‖b‖ ≤ 1
@@ -177,12 +177,12 @@ function lslq!(solver :: LslqSolver{T,FC,S}, A, b :: AbstractVector{FC};
   ktypeof(b) == S || error("ktypeof(b) ≠ $S")
 
   # Compute the adjoint of A
-  Aᵀ = A'
+  Aᴴ = A'
 
   # Set up workspace.
   allocate_if(!MisI, solver, :u, S, m)
   allocate_if(!NisI, solver, :v, S, n)
-  x, Nv, Aᵀu, w̄ = solver.x, solver.Nv, solver.Aᵀu, solver.w̄
+  x, Nv, Aᴴu, w̄ = solver.x, solver.Nv, solver.Aᴴu, solver.w̄
   Mu, Av, err_vec, stats = solver.Mu, solver.Av, solver.err_vec, solver.stats
   rNorms, ArNorms, err_lbnds = stats.residuals, stats.Aresiduals, stats.err_lbnds
   err_ubnds_lq, err_ubnds_cg = stats.err_ubnds_lq, stats.err_ubnds_cg
@@ -213,12 +213,12 @@ function lslq!(solver :: LslqSolver{T,FC,S}, A, b :: AbstractVector{FC};
 
   @kscal!(m, one(FC)/β₁, u)
   MisI || @kscal!(m, one(FC)/β₁, Mu)
-  mul!(Aᵀu, Aᵀ, u)
-  Nv .= Aᵀu
+  mul!(Aᴴu, Aᴴ, u)
+  Nv .= Aᴴu
   NisI || mulorldiv!(v, N, Nv, ldiv)
   α = sqrt(@kdotr(n, v, Nv))  # = α₁
 
-  # Aᵀb = 0 so x = 0 is a minimum least-squares solution
+  # Aᴴb = 0 so x = 0 is a minimum least-squares solution
   if α == 0
     stats.niter = 0
     stats.solved, stats.inconsistent = true, false
@@ -274,7 +274,7 @@ function lslq!(solver :: LslqSolver{T,FC,S}, A, b :: AbstractVector{FC};
   iter = 0
   itmax == 0 && (itmax = m + n)
 
-  (verbose > 0) && @printf("%5s  %7s  %7s  %7s  %7s  %8s  %8s  %7s  %7s  %7s\n", "k", "‖r‖", "‖Aᵀr‖", "β", "α", "cos", "sin", "‖A‖²", "κ(A)", "‖xL‖")
+  (verbose > 0) && @printf("%5s  %7s  %7s  %7s  %7s  %8s  %8s  %7s  %7s  %7s\n", "k", "‖r‖", "‖Aᴴr‖", "β", "α", "cos", "sin", "‖A‖²", "κ(A)", "‖xL‖")
   kdisplay(iter, verbose) && @printf("%5d  %7.1e  %7.1e  %7.1e  %7.1e  %8.1e  %8.1e  %7.1e  %7.1e  %7.1e\n", iter, rNorm, ArNorm, β, α, c, s, Anorm², Acond, xlqNorm)
 
   status = "unknown"
@@ -298,9 +298,9 @@ function lslq!(solver :: LslqSolver{T,FC,S}, A, b :: AbstractVector{FC};
       @kscal!(m, one(FC)/β, u)
       MisI || @kscal!(m, one(FC)/β, Mu)
 
-      # 2. αₖ₊₁Nvₖ₊₁ = Aᵀuₖ₊₁ - βₖ₊₁Nvₖ
-      mul!(Aᵀu, Aᵀ, u)
-      @kaxpby!(n, one(FC), Aᵀu, -β, Nv)
+      # 2. αₖ₊₁Nvₖ₊₁ = Aᴴuₖ₊₁ - βₖ₊₁Nvₖ
+      mul!(Aᴴu, Aᴴ, u)
+      @kaxpby!(n, one(FC), Aᴴu, -β, Nv)
       NisI || mulorldiv!(v, N, Nv, ldiv)
       α = sqrt(@kdotr(n, v, Nv))
       if α ≠ 0
diff --git a/src/lsmr.jl b/src/lsmr.jl
index f4d8349d1..78db5db59 100644
--- a/src/lsmr.jl
+++ b/src/lsmr.jl
@@ -5,7 +5,7 @@
 #
 # equivalently, of the normal equations
 #
-#  AᵀAx = Aᵀb.
+#  AᴴAx = Aᴴb.
 #
 # LSMR is formally equivalent to applying MINRES to the normal equations
 # but should be more stable. It is also formally equivalent to CRLS though
@@ -46,21 +46,21 @@ Solve the regularized linear least-squares problem
 using the LSMR method, where λ ≥ 0 is a regularization parameter.
 LSMR is formally equivalent to applying MINRES to the normal equations
 
-    (AᵀA + λ²I) x = Aᵀb
+    (AᴴA + λ²I) x = Aᴴb
 
 (and therefore to CRLS) but is more stable.
 
-LSMR produces monotonic residuals ‖r‖₂ and optimality residuals ‖Aᵀr‖₂.
+LSMR produces monotonic residuals ‖r‖₂ and optimality residuals ‖Aᴴr‖₂.
 It is formally equivalent to CRLS, though can be substantially more accurate.
 
 LSMR can be also used to find a null vector of a singular matrix A
-by solving the problem `min ‖Aᵀx - b‖` with any nonzero vector `b`.
-At a minimizer, the residual vector `r = b - Aᵀx` will satisfy `Ar = 0`.
+by solving the problem `min ‖Aᴴx - b‖` with any nonzero vector `b`.
+At a minimizer, the residual vector `r = b - Aᴴx` will satisfy `Ar = 0`.
 
 If `λ > 0`, we solve the symmetric and quasi-definite system
 
     [ E      A ] [ r ]   [ b ]
-    [ Aᵀ  -λ²F ] [ x ] = [ 0 ],
+    [ Aᴴ  -λ²F ] [ x ] = [ 0 ],
 
 where E and F are symmetric and positive definite.
 Preconditioners M = E⁻¹ ≻ 0 and N = F⁻¹ ≻ 0 may be provided in the form of linear operators.
@@ -70,19 +70,19 @@ The system above represents the optimality conditions of
 
     minimize ‖b - Ax‖²_E⁻¹ + λ²‖x‖²_F.
 
-For a symmetric and positive definite matrix `K`, the K-norm of a vector `x` is `‖x‖²_K = xᵀKx`.
-LSMR is then equivalent to applying MINRES to `(AᵀE⁻¹A + λ²F)x = AᵀE⁻¹b` with `r = E⁻¹(b - Ax)`.
+For a symmetric and positive definite matrix `K`, the K-norm of a vector `x` is `‖x‖²_K = xᴴKx`.
+LSMR is then equivalent to applying MINRES to `(AᴴE⁻¹A + λ²F)x = AᴴE⁻¹b` with `r = E⁻¹(b - Ax)`.
 
 If `λ = 0`, we solve the symmetric and indefinite system
 
     [ E    A ] [ r ]   [ b ]
-    [ Aᵀ   0 ] [ x ] = [ 0 ].
+    [ Aᴴ   0 ] [ x ] = [ 0 ].
 
 The system above represents the optimality conditions of
 
     minimize ‖b - Ax‖²_E⁻¹.
 
-In this case, `N` can still be specified and indicates the weighted norm in which `x` and `Aᵀr` should be measured.
+In this case, `N` can still be specified and indicates the weighted norm in which `x` and `Aᴴr` should be measured.
 `r` can be recovered by computing `E⁻¹(b - Ax)`.
 
 The callback is called as `callback(solver)` and should return `true` if the main loop should terminate,
@@ -134,12 +134,12 @@ function lsmr!(solver :: LsmrSolver{T,FC,S}, A, b :: AbstractVector{FC};
   ktypeof(b) == S || error("ktypeof(b) ≠ $S")
 
   # Compute the adjoint of A
-  Aᵀ = A'
+  Aᴴ = A'
 
   # Set up workspace.
   allocate_if(!MisI, solver, :u, S, m)
   allocate_if(!NisI, solver, :v, S, n)
-  x, Nv, Aᵀu, h, hbar = solver.x, solver.Nv, solver.Aᵀu, solver.h, solver.hbar
+  x, Nv, Aᴴu, h, hbar = solver.x, solver.Nv, solver.Aᴴu, solver.h, solver.hbar
   Mu, Av, err_vec, stats = solver.Mu, solver.Av, solver.err_vec, solver.stats
   rNorms, ArNorms = stats.residuals, stats.Aresiduals
   reset!(stats)
@@ -166,8 +166,8 @@ function lsmr!(solver :: LsmrSolver{T,FC,S}, A, b :: AbstractVector{FC};
 
   @kscal!(m, one(FC)/β₁, u)
   MisI || @kscal!(m, one(FC)/β₁, Mu)
-  mul!(Aᵀu, Aᵀ, u)
-  Nv .= Aᵀu
+  mul!(Aᴴu, Aᴴ, u)
+  Nv .= Aᴴu
   NisI || mulorldiv!(v, N, Nv, ldiv)
   α = sqrt(@kdotr(n, v, Nv))
 
@@ -210,10 +210,10 @@ function lsmr!(solver :: LsmrSolver{T,FC,S}, A, b :: AbstractVector{FC};
   iter = 0
   itmax == 0 && (itmax = m + n)
 
-  (verbose > 0) && @printf("%5s  %7s  %7s  %7s  %7s  %8s  %8s  %7s\n", "k", "‖r‖", "‖Aᵀr‖", "β", "α", "cos", "sin", "‖A‖²")
+  (verbose > 0) && @printf("%5s  %7s  %7s  %7s  %7s  %8s  %8s  %7s\n", "k", "‖r‖", "‖Aᴴr‖", "β", "α", "cos", "sin", "‖A‖²")
   kdisplay(iter, verbose) && @printf("%5d  %7.1e  %7.1e  %7.1e  %7.1e  %8.1e  %8.1e  %7.1e\n", iter, β₁, α, β₁, α, 0, 1, Anorm²)
 
-  # Aᵀb = 0 so x = 0 is a minimum least-squares solution
+  # Aᴴb = 0 so x = 0 is a minimum least-squares solution
   if α == 0
     stats.niter = 0
     stats.solved, stats.inconsistent = true, false
@@ -248,9 +248,9 @@ function lsmr!(solver :: LsmrSolver{T,FC,S}, A, b :: AbstractVector{FC};
       @kscal!(m, one(FC)/β, u)
       MisI || @kscal!(m, one(FC)/β, Mu)
 
-      # 2. αₖ₊₁Nvₖ₊₁ = Aᵀuₖ₊₁ - βₖ₊₁Nvₖ
-      mul!(Aᵀu, Aᵀ, u)
-      @kaxpby!(n, one(FC), Aᵀu, -β, Nv)
+      # 2. αₖ₊₁Nvₖ₊₁ = Aᴴuₖ₊₁ - βₖ₊₁Nvₖ
+      mul!(Aᴴu, Aᴴ, u)
+      @kaxpby!(n, one(FC), Aᴴu, -β, Nv)
       NisI || mulorldiv!(v, N, Nv, ldiv)
       α = sqrt(@kdotr(n, v, Nv))
       if α ≠ 0
diff --git a/src/lsqr.jl b/src/lsqr.jl
index dd3779dce..083b2f9f9 100644
--- a/src/lsqr.jl
+++ b/src/lsqr.jl
@@ -5,7 +5,7 @@
 #
 # equivalently, of the normal equations
 #
-#  AᵀAx = Aᵀb.
+#  AᴴAx = Aᴴb.
 #
 # LSQR is formally equivalent to applying the conjugate gradient method
 # to the normal equations but should be more stable. It is also formally
@@ -45,17 +45,17 @@ Solve the regularized linear least-squares problem
 using the LSQR method, where λ ≥ 0 is a regularization parameter.
 LSQR is formally equivalent to applying CG to the normal equations
 
-    (AᵀA + λ²I) x = Aᵀb
+    (AᴴA + λ²I) x = Aᴴb
 
 (and therefore to CGLS) but is more stable.
 
-LSQR produces monotonic residuals ‖r‖₂ but not optimality residuals ‖Aᵀr‖₂.
+LSQR produces monotonic residuals ‖r‖₂ but not optimality residuals ‖Aᴴr‖₂.
 It is formally equivalent to CGLS, though can be slightly more accurate.
 
 If `λ > 0`, LSQR solves the symmetric and quasi-definite system
 
     [ E      A ] [ r ]   [ b ]
-    [ Aᵀ  -λ²F ] [ x ] = [ 0 ],
+    [ Aᴴ  -λ²F ] [ x ] = [ 0 ],
 
 where E and F are symmetric and positive definite.
 Preconditioners M = E⁻¹ ≻ 0 and N = F⁻¹ ≻ 0 may be provided in the form of linear operators.
@@ -65,19 +65,19 @@ The system above represents the optimality conditions of
 
     minimize ‖b - Ax‖²_E⁻¹ + λ²‖x‖²_F.
 
-For a symmetric and positive definite matrix `K`, the K-norm of a vector `x` is `‖x‖²_K = xᵀKx`.
-LSQR is then equivalent to applying CG to `(AᵀE⁻¹A + λ²F)x = AᵀE⁻¹b` with `r = E⁻¹(b - Ax)`.
+For a symmetric and positive definite matrix `K`, the K-norm of a vector `x` is `‖x‖²_K = xᴴKx`.
+LSQR is then equivalent to applying CG to `(AᴴE⁻¹A + λ²F)x = AᴴE⁻¹b` with `r = E⁻¹(b - Ax)`.
 
 If `λ = 0`, we solve the symmetric and indefinite system
 
     [ E    A ] [ r ]   [ b ]
-    [ Aᵀ   0 ] [ x ] = [ 0 ].
+    [ Aᴴ   0 ] [ x ] = [ 0 ].
 
 The system above represents the optimality conditions of
 
     minimize ‖b - Ax‖²_E⁻¹.
 
-In this case, `N` can still be specified and indicates the weighted norm in which `x` and `Aᵀr` should be measured.
+In this case, `N` can still be specified and indicates the weighted norm in which `x` and `Aᴴr` should be measured.
 `r` can be recovered by computing `E⁻¹(b - Ax)`.
 
 The callback is called as `callback(solver)` and should return `true` if the main loop should terminate,
@@ -129,12 +129,12 @@ function lsqr!(solver :: LsqrSolver{T,FC,S}, A, b :: AbstractVector{FC};
   ktypeof(b) == S || error("ktypeof(b) ≠ $S")
 
   # Compute the adjoint of A
-  Aᵀ = A'
+  Aᴴ = A'
 
   # Set up workspace.
   allocate_if(!MisI, solver, :u, S, m)
   allocate_if(!NisI, solver, :v, S, n)
-  x, Nv, Aᵀu, w = solver.x, solver.Nv, solver.Aᵀu, solver.w
+  x, Nv, Aᴴu, w = solver.x, solver.Nv, solver.Aᴴu, solver.w
   Mu, Av, err_vec, stats = solver.Mu, solver.Av, solver.err_vec, solver.stats
   rNorms, ArNorms = stats.residuals, stats.Aresiduals
   reset!(stats)
@@ -162,8 +162,8 @@ function lsqr!(solver :: LsqrSolver{T,FC,S}, A, b :: AbstractVector{FC};
 
   @kscal!(m, one(FC)/β₁, u)
   MisI || @kscal!(m, one(FC)/β₁, Mu)
-  mul!(Aᵀu, Aᵀ, u)
-  Nv .= Aᵀu
+  mul!(Aᴴu, Aᴴ, u)
+  Nv .= Aᴴu
   NisI || mulorldiv!(v, N, Nv, ldiv)
   Anorm² = @kdotr(n, v, Nv)
   Anorm = sqrt(Anorm²)
@@ -184,7 +184,7 @@ function lsqr!(solver :: LsqrSolver{T,FC,S}, A, b :: AbstractVector{FC};
   iter = 0
   itmax == 0 && (itmax = m + n)
 
-  (verbose > 0) && @printf("%5s  %7s  %7s  %7s  %7s  %7s  %7s  %7s  %7s\n", "k", "α", "β", "‖r‖", "‖Aᵀr‖", "compat", "backwrd", "‖A‖", "κ(A)")
+  (verbose > 0) && @printf("%5s  %7s  %7s  %7s  %7s  %7s  %7s  %7s  %7s\n", "k", "α", "β", "‖r‖", "‖Aᴴr‖", "compat", "backwrd", "‖A‖", "κ(A)")
   kdisplay(iter, verbose) && @printf("%5d  %7.1e  %7.1e  %7.1e  %7.1e  %7.1e  %7.1e  %7.1e  %7.1e\n", iter, β₁, α, β₁, α, 0, 1, Anorm, Acond)
 
   rNorm = β₁
@@ -194,7 +194,7 @@ function lsqr!(solver :: LsqrSolver{T,FC,S}, A, b :: AbstractVector{FC};
   history && push!(rNorms, r2Norm)
   ArNorm = ArNorm0 = α * β
   history && push!(ArNorms, ArNorm)
-  # Aᵀb = 0 so x = 0 is a minimum least-squares solution
+  # Aᴴb = 0 so x = 0 is a minimum least-squares solution
   if α == 0
     stats.niter = 0
     stats.solved, stats.inconsistent = true, false
@@ -237,9 +237,9 @@ function lsqr!(solver :: LsqrSolver{T,FC,S}, A, b :: AbstractVector{FC};
       Anorm² = Anorm² + α * α + β * β  # = ‖B_{k-1}‖²
       λ > 0 && (Anorm² += λ²)
 
-      # 2. αₖ₊₁Nvₖ₊₁ = Aᵀuₖ₊₁ - βₖ₊₁Nvₖ
-      mul!(Aᵀu, Aᵀ, u)
-      @kaxpby!(n, one(FC), Aᵀu, -β, Nv)
+      # 2. αₖ₊₁Nvₖ₊₁ = Aᴴuₖ₊₁ - βₖ₊₁Nvₖ
+      mul!(Aᴴu, Aᴴ, u)
+      @kaxpby!(n, one(FC), Aᴴu, -β, Nv)
       NisI || mulorldiv!(v, N, Nv, ldiv)
       α = sqrt(@kdotr(n, v, Nv))
       if α ≠ 0
diff --git a/src/minres.jl b/src/minres.jl
index cbaefee9f..d3b8732ee 100644
--- a/src/minres.jl
+++ b/src/minres.jl
@@ -50,7 +50,7 @@ MINRES is formally equivalent to applying CR to Ax=b when A is positive
 definite, but is typically more stable and also applies to the case where
 A is indefinite.
 
-MINRES produces monotonic residuals ‖r‖₂ and optimality residuals ‖Aᵀr‖₂.
+MINRES produces monotonic residuals ‖r‖₂ and optimality residuals ‖Aᴴr‖₂.
 
 A preconditioner M may be provided in the form of a linear operator and is
 assumed to be symmetric and positive definite.
@@ -189,7 +189,7 @@ function minres!(solver :: MinresSolver{T,FC,S}, A, b :: AbstractVector{FC};
   iter = 0
   itmax == 0 && (itmax = 2*n)
 
-  (verbose > 0) && @printf("%5s  %7s  %7s  %7s  %8s  %8s  %7s  %7s  %7s  %7s\n", "k", "‖r‖", "‖Aᵀr‖", "β", "cos", "sin", "‖A‖", "κ(A)", "test1", "test2")
+  (verbose > 0) && @printf("%5s  %7s  %7s  %7s  %8s  %8s  %7s  %7s  %7s  %7s\n", "k", "‖r‖", "‖Aᴴr‖", "β", "cos", "sin", "‖A‖", "κ(A)", "test1", "test2")
   kdisplay(iter, verbose) && @printf("%5d  %7.1e  %7.1e  %7.1e  %8.1e  %8.1e  %7.1e  %7.1e\n", iter, rNorm, ArNorm, β, cs, sn, ANorm, Acond)
 
   tol = atol + rtol * β₁
@@ -241,7 +241,7 @@ function minres!(solver :: MinresSolver{T,FC,S}, A, b :: AbstractVector{FC};
     ϵ = sn * β
     δbar = -cs * β
     root = sqrt(γbar * γbar + δbar * δbar)
-    ArNorm = ϕbar * root  # = ‖Aᵀrₖ₋₁‖
+    ArNorm = ϕbar * root  # = ‖Aᴴrₖ₋₁‖
     history && push!(ArNorms, ArNorm)
 
     # Compute the next plane rotation.
@@ -295,7 +295,7 @@ function minres!(solver :: MinresSolver{T,FC,S}, A, b :: AbstractVector{FC};
     kdisplay(iter, verbose) && @printf("%5d  %7.1e  %7.1e  %7.1e  %8.1e  %8.1e  %7.1e  %7.1e  %7.1e  %7.1e\n", iter, rNorm, ArNorm, β, cs, sn, ANorm, Acond, test1, test2)
 
     if iter == 1 && β / β₁ ≤ 10 * ϵM
-      # Aᵀb = 0 so x = 0 is a minimum least-squares solution
+      # Aᴴb = 0 so x = 0 is a minimum least-squares solution
       stats.niter = 0
       stats.solved, stats.inconsistent = true, true
       stats.status = "x is a minimum least-squares solution"
diff --git a/src/minres_qlp.jl b/src/minres_qlp.jl
index bbfbf856b..509a7ef4e 100644
--- a/src/minres_qlp.jl
+++ b/src/minres_qlp.jl
@@ -246,7 +246,7 @@ function minres_qlp!(solver :: MinresQlpSolver{T,FC,S}, A, b :: AbstractVector{F
     # [sₖ -cₖ] [βₖ₊₁ ]   [0 ]
     (cₖ, sₖ, λₖ) = sym_givens(λbarₖ, βₖ₊₁)
 
-    # Compute [   zₖ  ] = (Qₖ)ᵀβ₁e₁
+    # Compute [   zₖ  ] = (Qₖ)ᴴβ₁e₁
     #         [ζbarₖ₊₁]
     #
     # [cₖ  sₖ] [ζbarₖ] = [   ζₖ  ]
@@ -312,7 +312,7 @@ function minres_qlp!(solver :: MinresQlpSolver{T,FC,S}, A, b :: AbstractVector{F
       τₖ   = (ξₖ - ψbarₖ₋₁ * τₖ₋₁) / μbarₖ
     end
 
-    # Compute directions wₖ₋₂, ẘₖ₋₁ and w̄ₖ, last columns of Wₖ = Vₖ(Pₖ)ᵀ
+    # Compute directions wₖ₋₂, ẘₖ₋₁ and w̄ₖ, last columns of Wₖ = Vₖ(Pₖ)ᴴ
     if iter == 1
       # w̅₁ = v₁
       @. wₖ = vₖ
diff --git a/src/qmr.jl b/src/qmr.jl
index eb4a4eb46..d4b684601 100644
--- a/src/qmr.jl
+++ b/src/qmr.jl
@@ -32,7 +32,7 @@ export qmr, qmr!
 Solve the square linear system Ax = b using the QMR method.
 
 QMR is based on the Lanczos biorthogonalization process and requires two initial vectors `b` and `c`.
-The relation `bᵀc ≠ 0` must be satisfied and by default `c = b`.
+The relation `bᴴc ≠ 0` must be satisfied and by default `c = b`.
 When `A` is symmetric and `b = c`, QMR is equivalent to MINRES.
 
 QMR can be warm-started from an initial guess `x0` with the method
@@ -96,7 +96,7 @@ function qmr!(solver :: QmrSolver{T,FC,S}, A, b :: AbstractVector{FC}; c :: Abst
   ktypeof(c) == S || error("ktypeof(c) ≠ $S")
 
   # Compute the adjoint of A
-  Aᵀ = A'
+  Aᴴ = A'
 
   # Set up workspace.
   uₖ₋₁, uₖ, q, vₖ₋₁, vₖ, p = solver.uₖ₋₁, solver.uₖ, solver.q, solver.vₖ₋₁, solver.vₖ, solver.p
@@ -133,18 +133,18 @@ function qmr!(solver :: QmrSolver{T,FC,S}, A, b :: AbstractVector{FC}; c :: Abst
   kdisplay(iter, verbose) && @printf("%5d  %7.1e\n", iter, rNorm)
 
   # Initialize the Lanczos biorthogonalization process.
-  cᵗb = @kdot(n, c, r₀)  # ⟨c,r₀⟩
-  if cᵗb == 0
+  cᴴb = @kdot(n, c, r₀)  # ⟨c,r₀⟩
+  if cᴴb == 0
     stats.niter = 0
     stats.solved = false
     stats.inconsistent = false
-    stats.status = "Breakdown bᵀc = 0"
+    stats.status = "Breakdown bᴴc = 0"
     solver.warm_start = false
     return solver
   end
 
-  βₖ = √(abs(cᵗb))             # β₁γ₁ = cᵀ(b - Ax₀)
-  γₖ = cᵗb / βₖ                # β₁γ₁ = cᵀ(b - Ax₀)
+  βₖ = √(abs(cᴴb))             # β₁γ₁ = cᴴ(b - Ax₀)
+  γₖ = cᴴb / βₖ                # β₁γ₁ = cᴴ(b - Ax₀)
   vₖ₋₁ .= zero(FC)             # v₀ = 0
   uₖ₋₁ .= zero(FC)             # u₀ = 0
   vₖ .= r₀ ./ βₖ               # v₁ = (b - Ax₀) / β₁
@@ -153,7 +153,7 @@ function qmr!(solver :: QmrSolver{T,FC,S}, A, b :: AbstractVector{FC}; c :: Abst
   sₖ₋₂ = sₖ₋₁ = sₖ = zero(FC)  # Givens sines used for the QR factorization of Tₖ₊₁.ₖ
   wₖ₋₂ .= zero(FC)             # Column k-2 of Wₖ = Vₖ(Rₖ)⁻¹
   wₖ₋₁ .= zero(FC)             # Column k-1 of Wₖ = Vₖ(Rₖ)⁻¹
-  ζbarₖ = βₖ                   # ζbarₖ is the last component of z̅ₖ = (Qₖ)ᵀβ₁e₁
+  ζbarₖ = βₖ                   # ζbarₖ is the last component of z̅ₖ = (Qₖ)ᴴβ₁e₁
   τₖ = @kdotr(n, vₖ, vₖ)       # τₖ is used for the residual norm estimate
 
   # Stopping criterion.
@@ -169,10 +169,10 @@ function qmr!(solver :: QmrSolver{T,FC,S}, A, b :: AbstractVector{FC}; c :: Abst
 
     # Continue the Lanczos biorthogonalization process.
     # AVₖ  = VₖTₖ    + βₖ₊₁vₖ₊₁(eₖ)ᵀ = Vₖ₊₁Tₖ₊₁.ₖ
-    # AᵀUₖ = Uₖ(Tₖ)ᵀ + γ̄ₖ₊₁uₖ₊₁(eₖ)ᵀ = Uₖ₊₁(Tₖ.ₖ₊₁)ᵀ
+    # AᴴUₖ = Uₖ(Tₖ)ᴴ + γ̄ₖ₊₁uₖ₊₁(eₖ)ᵀ = Uₖ₊₁(Tₖ.ₖ₊₁)ᴴ
 
     mul!(q, A , vₖ)  # Forms vₖ₊₁ : q ← Avₖ
-    mul!(p, Aᵀ, uₖ)  # Forms uₖ₊₁ : p ← Aᵀuₖ
+    mul!(p, Aᴴ, uₖ)  # Forms uₖ₊₁ : p ← Aᴴuₖ
 
     @kaxpy!(n, -γₖ, vₖ₋₁, q)  # q ← q - γₖ * vₖ₋₁
     @kaxpy!(n, -βₖ, uₖ₋₁, p)  # p ← p - β̄ₖ * uₖ₋₁
@@ -182,9 +182,9 @@ function qmr!(solver :: QmrSolver{T,FC,S}, A, b :: AbstractVector{FC}; c :: Abst
     @kaxpy!(n, -     αₖ , vₖ, q)    # q ← q - αₖ * vₖ
     @kaxpy!(n, -conj(αₖ), uₖ, p)    # p ← p - ᾱₖ * uₖ
 
-    pᵗq = @kdot(n, p, q)      # pᵗq  = ⟨p,q⟩
-    βₖ₊₁ = √(abs(pᵗq))        # βₖ₊₁ = √(|pᵗq|)
-    γₖ₊₁ = pᵗq / βₖ₊₁         # γₖ₊₁ = pᵗq / βₖ₊₁
+    pᴴq = @kdot(n, p, q)      # pᴴq  = ⟨p,q⟩
+    βₖ₊₁ = √(abs(pᴴq))        # βₖ₊₁ = √(|pᴴq|)
+    γₖ₊₁ = pᴴq / βₖ₊₁         # γₖ₊₁ = pᴴq / βₖ₊₁
 
     # Update the QR factorization of Tₖ₊₁.ₖ = Qₖ [ Rₖ ].
     #                                            [ Oᵀ ]
@@ -271,7 +271,7 @@ function qmr!(solver :: QmrSolver{T,FC,S}, A, b :: AbstractVector{FC}; c :: Abst
     @. vₖ₋₁ = vₖ  # vₖ₋₁ ← vₖ
     @. uₖ₋₁ = uₖ  # uₖ₋₁ ← uₖ
 
-    if pᵗq ≠ zero(FC)
+    if pᴴq ≠ zero(FC)
       @. vₖ = q / βₖ₊₁        # βₖ₊₁vₖ₊₁ = q
       @. uₖ = p / conj(γₖ₊₁)  # γ̄ₖ₊₁uₖ₊₁ = p
     end
@@ -303,7 +303,7 @@ function qmr!(solver :: QmrSolver{T,FC,S}, A, b :: AbstractVector{FC}; c :: Abst
     resid_decrease_lim = rNorm ≤ ε
     solved = resid_decrease_lim || resid_decrease_mach
     tired = iter ≥ itmax
-    breakdown = !solved && (pᵗq == 0)
+    breakdown = !solved && (pᴴq == 0)
     kdisplay(iter, verbose) && @printf("%5d  %7.1e\n", iter, rNorm)
   end
   (verbose > 0) && @printf("\n")
diff --git a/src/tricg.jl b/src/tricg.jl
index 5acff2d52..7c140a821 100644
--- a/src/tricg.jl
+++ b/src/tricg.jl
@@ -25,7 +25,7 @@ export tricg, tricg!
 TriCG solves the symmetric linear system
 
     [ τE    A ] [ x ] = [ b ]
-    [  Aᵀ  νF ] [ y ]   [ c ],
+    [  Aᴴ  νF ] [ y ]   [ c ],
 
 where τ and ν are real numbers, E = M⁻¹ ≻ 0 and F = N⁻¹ ≻ 0.
 `b` and `c` must both be nonzero.
@@ -133,7 +133,7 @@ function tricg!(solver :: TricgSolver{T,FC,S}, A, b :: AbstractVector{FC}, c ::
   warm_start && (ν ≠ 0) && !NisI && error("Warm-start with preconditioners is not supported.")
 
   # Compute the adjoint of A
-  Aᵀ = A'
+  Aᴴ = A'
 
   # Set up workspace.
   allocate_if(!MisI, solver, :vₖ, S, m)
@@ -164,12 +164,12 @@ function tricg!(solver :: TricgSolver{T,FC,S}, A, b :: AbstractVector{FC}, c ::
   N⁻¹uₖ₋₁ .= zero(FC)  # u₀ = 0
 
   # [ τI    A ] [ xₖ ] = [ b -  τΔx - AΔy ] = [ b₀ ]
-  # [  Aᵀ  νI ] [ yₖ ]   [ c - AᵀΔx - νΔy ]   [ c₀ ]
+  # [  Aᴴ  νI ] [ yₖ ]   [ c - AᴴΔx - νΔy ]   [ c₀ ]
   if warm_start
     mul!(b₀, A, Δy)
     (τ ≠ 0) && @kaxpy!(m, τ, Δx, b₀)
     @kaxpby!(m, one(FC), b, -one(FC), b₀)
-    mul!(c₀, Aᵀ, Δx)
+    mul!(c₀, Aᴴ, Δx)
     (ν ≠ 0) && @kaxpy!(n, ν, Δy, c₀)
     @kaxpby!(n, one(FC), c, -one(FC), c₀)
   end
@@ -196,7 +196,7 @@ function tricg!(solver :: TricgSolver{T,FC,S}, A, b :: AbstractVector{FC}, c ::
     error("c must be nonzero")
   end
 
-  # Initialize directions Gₖ such that Lₖ(Gₖ)ᵀ = (Wₖ)ᵀ
+  # Initialize directions Gₖ such that L̄ₖ(Gₖ)ᵀ = (Wₖ)ᵀ
   gx₂ₖ₋₁ .= zero(FC)
   gy₂ₖ₋₁ .= zero(FC)
   gx₂ₖ   .= zero(FC)
@@ -231,10 +231,10 @@ function tricg!(solver :: TricgSolver{T,FC,S}, A, b :: AbstractVector{FC}, c ::
 
     # Continue the orthogonal tridiagonalization process.
     # AUₖ  = EVₖTₖ    + βₖ₊₁Evₖ₊₁(eₖ)ᵀ = EVₖ₊₁Tₖ₊₁.ₖ
-    # AᵀVₖ = FUₖ(Tₖ)ᵀ + γₖ₊₁Fuₖ₊₁(eₖ)ᵀ = FUₖ₊₁(Tₖ.ₖ₊₁)ᵀ
+    # AᴴVₖ = FUₖ(Tₖ)ᴴ + γₖ₊₁Fuₖ₊₁(eₖ)ᵀ = FUₖ₊₁(Tₖ.ₖ₊₁)ᴴ
 
     mul!(q, A , uₖ)  # Forms Evₖ₊₁ : q ← Auₖ
-    mul!(p, Aᵀ, vₖ)  # Forms Fuₖ₊₁ : p ← Aᵀvₖ
+    mul!(p, Aᴴ, vₖ)  # Forms Fuₖ₊₁ : p ← Aᴴvₖ
 
     if iter ≥ 2
       @kaxpy!(m, -γₖ, M⁻¹vₖ₋₁, q)  # q ← q - γₖ * M⁻¹vₖ₋₁
@@ -254,14 +254,14 @@ function tricg!(solver :: TricgSolver{T,FC,S}, A, b :: AbstractVector{FC}, c ::
     #                                [0  u₁ ••• 0  uₖ]
     #
     # rₖ = [ b ] - [ τE    A ] [ xₖ ] = [ b ] - [ τE    A ] Wₖzₖ
-    #      [ c ]   [  Aᵀ  νF ] [ yₖ ]   [ c ]   [  Aᵀ  νF ]
+    #      [ c ]   [  Aᴴ  νF ] [ yₖ ]   [ c ]   [  Aᴴ  νF ]
     #
     # block-Lanczos formulation : [ τE    A ] Wₖ = [ E   0 ] Wₖ₊₁Sₖ₊₁.ₖ
-    #                             [  Aᵀ  νF ]      [ 0   F ]
+    #                             [  Aᴴ  νF ]      [ 0   F ]
     #
-    # TriCG subproblem : (Wₖ)ᵀ * rₖ = 0 ↔ Sₖ.ₖzₖ = β₁e₁ + γ₁e₂
+    # TriCG subproblem : (Wₖ)ᴴ * rₖ = 0 ↔ Sₖ.ₖzₖ = β₁e₁ + γ₁e₂
     #
-    # Update the LDLᵀ factorization of Sₖ.ₖ.
+    # Update the LDLᴴ factorization of Sₖ.ₖ.
     #
     # [ τ  α₁    γ₂ 0  •  •  •  •  0  ]
     # [ ᾱ₁ ν  β₂       •           •  ]
@@ -306,7 +306,7 @@ function tricg!(solver :: TricgSolver{T,FC,S}, A, b :: AbstractVector{FC}, c ::
       π₂ₖ   = -(δₖ * d₂ₖ₋₁ * π₂ₖ₋₁ + λₖ * d₂ₖ₋₂ * π₂ₖ₋₂ + ηₖ * d₂ₖ₋₃ * π₂ₖ₋₃) / d₂ₖ
     end
 
-    # Solve Gₖ = Wₖ(Lₖ)⁻ᵀ ⟷ L̄ₖ(Gₖ)ᵀ = (Wₖ)ᵀ.
+    # Solve Gₖ = Wₖ(Lₖ)⁻ᴴ ⟷ L̄ₖ(Gₖ)ᵀ = (Wₖ)ᵀ.
     if iter == 1
       # [ 1  0 ] [ gx₁ gy₁ ] = [ v₁ 0  ]
       # [ δ̄₁ 1 ] [ gx₂ gy₂ ]   [ 0  u₁ ]
@@ -342,7 +342,7 @@ function tricg!(solver :: TricgSolver{T,FC,S}, A, b :: AbstractVector{FC}, c ::
 
     # Compute vₖ₊₁ and uₖ₊₁
     MisI || mulorldiv!(vₖ₊₁, M, q, ldiv)  # βₖ₊₁vₖ₊₁ = MAuₖ  - γₖvₖ₋₁ - αₖvₖ
-    NisI || mulorldiv!(uₖ₊₁, N, p, ldiv)  # γₖ₊₁uₖ₊₁ = NAᵀvₖ - βₖuₖ₋₁ - ᾱₖuₖ
+    NisI || mulorldiv!(uₖ₊₁, N, p, ldiv)  # γₖ₊₁uₖ₊₁ = NAᴴvₖ - βₖuₖ₋₁ - ᾱₖuₖ
 
     βₖ₊₁ = sqrt(@kdotr(m, vₖ₊₁, q))  # βₖ₊₁ = ‖vₖ₊₁‖_E
     γₖ₊₁ = sqrt(@kdotr(n, uₖ₊₁, p))  # γₖ₊₁ = ‖uₖ₊₁‖_F
diff --git a/src/trilqr.jl b/src/trilqr.jl
index edcb4c9b9..6b0948984 100644
--- a/src/trilqr.jl
+++ b/src/trilqr.jl
@@ -1,5 +1,5 @@
 # An implementation of TRILQR for the solution of square or
-# rectangular consistent linear adjoint systems Ax = b and Aᵀy = c.
+# rectangular consistent linear adjoint systems Ax = b and Aᴴy = c.
 #
 # This method is described in
 #
@@ -24,10 +24,10 @@ export trilqr, trilqr!
 Combine USYMLQ and USYMQR to solve adjoint systems.
 
     [0  A] [y] = [b]
-    [Aᵀ 0] [x]   [c]
+    [Aᴴ 0] [x]   [c]
 
 USYMLQ is used for solving primal system `Ax = b`.
-USYMQR is used for solving dual system `Aᵀy = c`.
+USYMQR is used for solving dual system `Aᴴy = c`.
 
 An option gives the possibility of transferring from the USYMLQ point to the
 USYMCG point, when it exists. The transfer is based on the residual norm.
@@ -93,7 +93,7 @@ function trilqr!(solver :: TrilqrSolver{T,FC,S}, A, b :: AbstractVector{FC}, c :
   ktypeof(c) == S || error("ktypeof(c) ≠ $S")
 
   # Compute the adjoint of A
-  Aᵀ = A'
+  Aᴴ = A'
 
   # Set up workspace.
   uₖ₋₁, uₖ, p, d̅, x, stats = solver.uₖ₋₁, solver.uₖ, solver.p, solver.d̅, solver.x, solver.stats
@@ -107,7 +107,7 @@ function trilqr!(solver :: TrilqrSolver{T,FC,S}, A, b :: AbstractVector{FC}, c :
   if warm_start
     mul!(r₀, A, Δx)
     @kaxpby!(n, one(FC), b, -one(FC), r₀)
-    mul!(s₀, Aᵀ, Δy)
+    mul!(s₀, Aᴴ, Δy)
     @kaxpby!(n, one(FC), c, -one(FC), s₀)
   end
 
@@ -115,7 +115,7 @@ function trilqr!(solver :: TrilqrSolver{T,FC,S}, A, b :: AbstractVector{FC}, c :
   x .= zero(FC)          # x₀
   bNorm = @knrm2(m, r₀)  # rNorm = ‖r₀‖
 
-  # Initial solution y₀ and residual s₀ = c - Aᵀy₀.
+  # Initial solution y₀ and residual s₀ = c - Aᴴy₀.
   t .= zero(FC)          # t₀
   cNorm = @knrm2(n, s₀)  # sNorm = ‖s₀‖
 
@@ -136,17 +136,17 @@ function trilqr!(solver :: TrilqrSolver{T,FC,S}, A, b :: AbstractVector{FC}, c :
   vₖ₋₁ .= zero(FC)            # v₀ = 0
   uₖ₋₁ .= zero(FC)            # u₀ = 0
   vₖ .= r₀ ./ βₖ              # v₁ = (b - Ax₀) / β₁
-  uₖ .= s₀ ./ γₖ              # u₁ = (c - Aᵀy₀) / γ₁
+  uₖ .= s₀ ./ γₖ              # u₁ = (c - Aᴴy₀) / γ₁
   cₖ₋₁ = cₖ = -one(T)         # Givens cosines used for the LQ factorization of Tₖ
   sₖ₋₁ = sₖ = zero(FC)        # Givens sines used for the LQ factorization of Tₖ
-  d̅ .= zero(FC)               # Last column of D̅ₖ = Uₖ(Qₖ)ᵀ
+  d̅ .= zero(FC)               # Last column of D̅ₖ = Uₖ(Qₖ)ᴴ
   ζₖ₋₁ = ζbarₖ = zero(FC)     # ζₖ₋₁ and ζbarₖ are the last components of z̅ₖ = (L̅ₖ)⁻¹β₁e₁
   ζₖ₋₂ = ηₖ = zero(FC)        # ζₖ₋₂ and ηₖ are used to update ζₖ₋₁ and ζbarₖ
   δbarₖ₋₁ = δbarₖ = zero(FC)  # Coefficients of Lₖ₋₁ and L̅ₖ modified over the course of two iterations
   ψbarₖ₋₁ = ψₖ₋₁ = zero(FC)   # ψₖ₋₁ and ψbarₖ are the last components of h̅ₖ = Qₖγ₁e₁
   ϵₖ₋₃ = λₖ₋₂ = zero(FC)      # Components of Lₖ₋₁
-  wₖ₋₃ .= zero(FC)            # Column k-3 of Wₖ = Vₖ(Lₖ)⁻ᵀ
-  wₖ₋₂ .= zero(FC)            # Column k-2 of Wₖ = Vₖ(Lₖ)⁻ᵀ
+  wₖ₋₃ .= zero(FC)            # Column k-3 of Wₖ = Vₖ(Lₖ)⁻ᴴ
+  wₖ₋₂ .= zero(FC)            # Column k-2 of Wₖ = Vₖ(Lₖ)⁻ᴴ
 
   # Stopping criterion.
   inconsistent = false
@@ -166,10 +166,10 @@ function trilqr!(solver :: TrilqrSolver{T,FC,S}, A, b :: AbstractVector{FC}, c :
 
     # Continue the SSY tridiagonalization process.
     # AUₖ  = VₖTₖ    + βₖ₊₁vₖ₊₁(eₖ)ᵀ = Vₖ₊₁Tₖ₊₁.ₖ
-    # AᵀVₖ = Uₖ(Tₖ)ᵀ + γₖ₊₁uₖ₊₁(eₖ)ᵀ = Uₖ₊₁(Tₖ.ₖ₊₁)ᵀ
+    # AᴴVₖ = Uₖ(Tₖ)ᴴ + γₖ₊₁uₖ₊₁(eₖ)ᵀ = Uₖ₊₁(Tₖ.ₖ₊₁)ᴴ
 
     mul!(q, A , uₖ)  # Forms vₖ₊₁ : q ← Auₖ
-    mul!(p, Aᵀ, vₖ)  # Forms uₖ₊₁ : p ← Aᵀvₖ
+    mul!(p, Aᴴ, vₖ)  # Forms uₖ₊₁ : p ← Aᴴvₖ
 
     @kaxpy!(m, -γₖ, vₖ₋₁, q)  # q ← q - γₖ * vₖ₋₁
     @kaxpy!(n, -βₖ, uₖ₋₁, p)  # p ← p - βₖ * uₖ₋₁
@@ -236,7 +236,7 @@ function trilqr!(solver :: TrilqrSolver{T,FC,S}, A, b :: AbstractVector{FC}, c :
         ηₖ   = -ϵₖ₋₂ * ζₖ₋₂ - λₖ₋₁ * ζₖ₋₁
       end
 
-      # Relations for the directions dₖ₋₁ and d̅ₖ, the last two columns of D̅ₖ = Uₖ(Qₖ)ᵀ.
+      # Relations for the directions dₖ₋₁ and d̅ₖ, the last two columns of D̅ₖ = Uₖ(Qₖ)ᴴ.
       # [d̅ₖ₋₁ uₖ] [cₖ  s̄ₖ] = [dₖ₋₁ d̅ₖ] ⟷ dₖ₋₁ = cₖ * d̅ₖ₋₁ + sₖ * uₖ
       #           [sₖ -cₖ]             ⟷ d̅ₖ   = s̄ₖ * d̅ₖ₋₁ - cₖ * uₖ
       if iter ≥ 2
@@ -295,7 +295,7 @@ function trilqr!(solver :: TrilqrSolver{T,FC,S}, A, b :: AbstractVector{FC}, c :
         ψbarₖ = sₖ * ψbarₖ₋₁
       end
 
-      # Compute the direction wₖ₋₁, the last column of Wₖ₋₁ = (Vₖ₋₁)(Lₖ₋₁)⁻ᵀ ⟷ (L̄ₖ₋₁)(Wₖ₋₁)ᵀ = (Vₖ₋₁)ᵀ.
+      # Compute the direction wₖ₋₁, the last column of Wₖ₋₁ = (Vₖ₋₁)(Lₖ₋₁)⁻ᴴ ⟷ (L̄ₖ₋₁)(Wₖ₋₁)ᵀ = (Vₖ₋₁)ᵀ.
       # w₁ = v₁ / δ̄₁
       if iter == 2
         wₖ₋₁ = wₖ₋₂
diff --git a/src/trimr.jl b/src/trimr.jl
index bc53633c2..7dd826edf 100644
--- a/src/trimr.jl
+++ b/src/trimr.jl
@@ -25,7 +25,7 @@ export trimr, trimr!
 TriMR solves the symmetric linear system
 
     [ τE    A ] [ x ] = [ b ]
-    [  Aᵀ  νF ] [ y ]   [ c ],
+    [  Aᴴ  νF ] [ y ]   [ c ],
 
 where τ and ν are real numbers, E = M⁻¹ ≻ 0, F = N⁻¹ ≻ 0.
 `b` and `c` must both be nonzero.
@@ -137,7 +137,7 @@ function trimr!(solver :: TrimrSolver{T,FC,S}, A, b :: AbstractVector{FC}, c ::
   warm_start && (ν ≠ 0) && !NisI && error("Warm-start with preconditioners is not supported.")
 
   # Compute the adjoint of A
-  Aᵀ = A'
+  Aᴴ = A'
 
   # Set up workspace.
   allocate_if(!MisI, solver, :vₖ, S, m)
@@ -169,12 +169,12 @@ function trimr!(solver :: TrimrSolver{T,FC,S}, A, b :: AbstractVector{FC}, c ::
   N⁻¹uₖ₋₁ .= zero(FC)  # u₀ = 0
 
   # [ τI    A ] [ xₖ ] = [ b -  τΔx - AΔy ] = [ b₀ ]
-  # [  Aᵀ  νI ] [ yₖ ]   [ c - AᵀΔx - νΔy ]   [ c₀ ]
+  # [  Aᴴ  νI ] [ yₖ ]   [ c - AᴴΔx - νΔy ]   [ c₀ ]
   if warm_start
     mul!(b₀, A, Δy)
     (τ ≠ 0) && @kaxpy!(m, τ, Δx, b₀)
     @kaxpby!(m, one(FC), b, -one(FC), b₀)
-    mul!(c₀, Aᵀ, Δx)
+    mul!(c₀, Aᴴ, Δx)
     (ν ≠ 0) && @kaxpy!(n, ν, Δy, c₀)
     @kaxpby!(n, one(FC), c, -one(FC), c₀)
   end
@@ -244,10 +244,10 @@ function trimr!(solver :: TrimrSolver{T,FC,S}, A, b :: AbstractVector{FC}, c ::
 
     # Continue the orthogonal tridiagonalization process.
     # AUₖ  = EVₖTₖ    + βₖ₊₁Evₖ₊₁(eₖ)ᵀ = EVₖ₊₁Tₖ₊₁.ₖ
-    # AᵀVₖ = FUₖ(Tₖ)ᵀ + γₖ₊₁Fuₖ₊₁(eₖ)ᵀ = FUₖ₊₁(Tₖ.ₖ₊₁)ᵀ
+    # AᴴVₖ = FUₖ(Tₖ)ᴴ + γₖ₊₁Fuₖ₊₁(eₖ)ᵀ = FUₖ₊₁(Tₖ.ₖ₊₁)ᴴ
 
     mul!(q, A , uₖ)  # Forms Evₖ₊₁ : q ← Auₖ
-    mul!(p, Aᵀ, vₖ)  # Forms Fuₖ₊₁ : p ← Aᵀvₖ
+    mul!(p, Aᴴ, vₖ)  # Forms Fuₖ₊₁ : p ← Aᴴvₖ
 
     if iter ≥ 2
       @kaxpy!(m, -γₖ, M⁻¹vₖ₋₁, q)  # q ← q - γₖ * M⁻¹vₖ₋₁
@@ -261,7 +261,7 @@ function trimr!(solver :: TrimrSolver{T,FC,S}, A, b :: AbstractVector{FC}, c ::
 
     # Compute vₖ₊₁ and uₖ₊₁
     MisI || mulorldiv!(vₖ₊₁, M, q, ldiv)  # βₖ₊₁vₖ₊₁ = MAuₖ  - γₖvₖ₋₁ - αₖvₖ
-    NisI || mulorldiv!(uₖ₊₁, N, p, ldiv)  # γₖ₊₁uₖ₊₁ = NAᵀvₖ - βₖuₖ₋₁ - ᾱₖuₖ
+    NisI || mulorldiv!(uₖ₊₁, N, p, ldiv)  # γₖ₊₁uₖ₊₁ = NAᴴvₖ - βₖuₖ₋₁ - ᾱₖuₖ
 
     βₖ₊₁ = sqrt(@kdotr(m, vₖ₊₁, q))  # βₖ₊₁ = ‖vₖ₊₁‖_E
     γₖ₊₁ = sqrt(@kdotr(n, uₖ₊₁, p))  # γₖ₊₁ = ‖uₖ₊₁‖_F
@@ -282,10 +282,10 @@ function trimr!(solver :: TrimrSolver{T,FC,S}, A, b :: AbstractVector{FC}, c ::
     #                                [0  u₁ ••• 0  uₖ]
     #
     # rₖ = [ b ] - [ τE    A ] [ xₖ ] = [ b ] - [ τE    A ] Wₖzₖ
-    #      [ c ]   [  Aᵀ  νF ] [ yₖ ]   [ c ]   [  Aᵀ  νF ]
+    #      [ c ]   [  Aᴴ  νF ] [ yₖ ]   [ c ]   [  Aᴴ  νF ]
     #
     # block-Lanczos formulation : [ τE    A ] Wₖ = [ E   0 ] Wₖ₊₁Sₖ₊₁.ₖ
-    #                             [  Aᵀ  νF ]      [ 0   F ]
+    #                             [  Aᴴ  νF ]      [ 0   F ]
     #
     # TriMR subproblem : min ‖ rₖ ‖ ↔ min ‖ Sₖ₊₁.ₖzₖ - β₁e₁ - γ₁e₂ ‖
     #
@@ -419,7 +419,7 @@ function trimr!(solver :: TrimrSolver{T,FC,S}, A, b :: AbstractVector{FC}, c ::
       @kswap(gy₂ₖ₋₂, gy₂ₖ)
     end
 
-    # Update p̅ₖ = (Qₖ)ᵀ * (β₁e₁ + γ₁e₂)
+    # Update p̅ₖ = (Qₖ)ᴴ * (β₁e₁ + γ₁e₂)
     πbis₂ₖ   =      c₁ₖ  * πbar₂ₖ
     πbis₂ₖ₊₂ = conj(s₁ₖ) * πbar₂ₖ
     #
diff --git a/src/usymlq.jl b/src/usymlq.jl
index 71670c80f..29cd704c7 100644
--- a/src/usymlq.jl
+++ b/src/usymlq.jl
@@ -31,7 +31,7 @@ export usymlq, usymlq!
 Solve the linear system Ax = b using the USYMLQ method.
 
 USYMLQ is based on the orthogonal tridiagonalization process and requires two initial nonzero vectors `b` and `c`.
-The vector `c` is only used to initialize the process and a default value can be `b` or `Aᵀb` depending on the shape of `A`.
+The vector `c` is only used to initialize the process and a default value can be `b` or `Aᴴb` depending on the shape of `A`.
 The error norm ‖x - x*‖ monotonously decreases in USYMLQ.
 It's considered as a generalization of SYMMLQ.
 
@@ -103,7 +103,7 @@ function usymlq!(solver :: UsymlqSolver{T,FC,S}, A, b :: AbstractVector{FC}, c :
   ktypeof(c) == S || error("ktypeof(c) ≠ $S")
 
   # Compute the adjoint of A
-  Aᵀ = A'
+  Aᴴ = A'
 
   # Set up workspace.
   uₖ₋₁, uₖ, p, Δx, x = solver.uₖ₋₁, solver.uₖ, solver.p, solver.Δx, solver.x
@@ -146,7 +146,7 @@ function usymlq!(solver :: UsymlqSolver{T,FC,S}, A, b :: AbstractVector{FC}, c :
   uₖ .= c ./ γₖ               # u₁ = c / γ₁
   cₖ₋₁ = cₖ = -one(T)         # Givens cosines used for the LQ factorization of Tₖ
   sₖ₋₁ = sₖ = zero(FC)        # Givens sines used for the LQ factorization of Tₖ
-  d̅ .= zero(FC)               # Last column of D̅ₖ = Uₖ(Qₖ)ᵀ
+  d̅ .= zero(FC)               # Last column of D̅ₖ = Uₖ(Qₖ)ᴴ
   ζₖ₋₁ = ζbarₖ = zero(FC)     # ζₖ₋₁ and ζbarₖ are the last components of z̅ₖ = (L̅ₖ)⁻¹β₁e₁
   ζₖ₋₂ = ηₖ = zero(FC)        # ζₖ₋₂ and ηₖ are used to update ζₖ₋₁ and ζbarₖ
   δbarₖ₋₁ = δbarₖ = zero(FC)  # Coefficients of Lₖ₋₁ and Lₖ modified over the course of two iterations
@@ -164,10 +164,10 @@ function usymlq!(solver :: UsymlqSolver{T,FC,S}, A, b :: AbstractVector{FC}, c :
 
     # Continue the SSY tridiagonalization process.
     # AUₖ  = VₖTₖ    + βₖ₊₁vₖ₊₁(eₖ)ᵀ = Vₖ₊₁Tₖ₊₁.ₖ
-    # AᵀVₖ = Uₖ(Tₖ)ᵀ + γₖ₊₁uₖ₊₁(eₖ)ᵀ = Uₖ₊₁(Tₖ.ₖ₊₁)ᵀ
+    # AᴴVₖ = Uₖ(Tₖ)ᴴ + γₖ₊₁uₖ₊₁(eₖ)ᵀ = Uₖ₊₁(Tₖ.ₖ₊₁)ᴴ
 
     mul!(q, A , uₖ)  # Forms vₖ₊₁ : q ← Auₖ
-    mul!(p, Aᵀ, vₖ)  # Forms uₖ₊₁ : p ← Aᵀvₖ
+    mul!(p, Aᴴ, vₖ)  # Forms uₖ₊₁ : p ← Aᴴvₖ
 
     @kaxpy!(m, -γₖ, vₖ₋₁, q)  # q ← q - γₖ * vₖ₋₁
     @kaxpy!(n, -βₖ, uₖ₋₁, p)  # p ← p - βₖ * uₖ₋₁
@@ -233,7 +233,7 @@ function usymlq!(solver :: UsymlqSolver{T,FC,S}, A, b :: AbstractVector{FC}, c :
       ηₖ   = -ϵₖ₋₂ * ζₖ₋₂ - λₖ₋₁ * ζₖ₋₁
     end
 
-    # Relations for the directions dₖ₋₁ and d̅ₖ, the last two columns of D̅ₖ = Uₖ(Qₖ)ᵀ.
+    # Relations for the directions dₖ₋₁ and d̅ₖ, the last two columns of D̅ₖ = Uₖ(Qₖ)ᴴ.
     # [d̅ₖ₋₁ uₖ] [cₖ  s̄ₖ] = [dₖ₋₁ d̅ₖ] ⟷ dₖ₋₁ = cₖ * d̅ₖ₋₁ + sₖ * uₖ
     #           [sₖ -cₖ]             ⟷ d̅ₖ   = s̄ₖ * d̅ₖ₋₁ - cₖ * uₖ
     if iter ≥ 2
diff --git a/src/usymqr.jl b/src/usymqr.jl
index 863390c3f..45c95c88d 100644
--- a/src/usymqr.jl
+++ b/src/usymqr.jl
@@ -31,7 +31,7 @@ export usymqr, usymqr!
 Solve the linear system Ax = b using the USYMQR method.
 
 USYMQR is based on the orthogonal tridiagonalization process and requires two initial nonzero vectors `b` and `c`.
-The vector `c` is only used to initialize the process and a default value can be `b` or `Aᵀb` depending on the shape of `A`.
+The vector `c` is only used to initialize the process and a default value can be `b` or `Aᴴb` depending on the shape of `A`.
 The residual norm ‖b - Ax‖ monotonously decreases in USYMQR.
 It's considered as a generalization of MINRES.
 
@@ -100,13 +100,13 @@ function usymqr!(solver :: UsymqrSolver{T,FC,S}, A, b :: AbstractVector{FC}, c :
   ktypeof(c) == S || error("ktypeof(c) ≠ $S")
 
   # Compute the adjoint of A
-  Aᵀ = A'
+  Aᴴ = A'
 
   # Set up workspace.
   vₖ₋₁, vₖ, q, Δx, x, p = solver.vₖ₋₁, solver.vₖ, solver.q, solver.Δx, solver.x, solver.p
   wₖ₋₂, wₖ₋₁, uₖ₋₁, uₖ, stats = solver.wₖ₋₂, solver.wₖ₋₁, solver.uₖ₋₁, solver.uₖ, solver.stats
   warm_start = solver.warm_start
-  rNorms, AᵀrNorms = stats.residuals, stats.Aresiduals
+  rNorms, AᴴrNorms = stats.residuals, stats.Aresiduals
   reset!(stats)
   r₀ = warm_start ? q : b
 
@@ -133,7 +133,7 @@ function usymqr!(solver :: UsymqrSolver{T,FC,S}, A, b :: AbstractVector{FC}, c :
 
   ε = atol + rtol * rNorm
   κ = zero(T)
-  (verbose > 0) && @printf("%5s  %7s  %7s\n", "k", "‖rₖ‖", "‖Aᵀrₖ₋₁‖")
+  (verbose > 0) && @printf("%5s  %7s  %7s\n", "k", "‖rₖ‖", "‖Aᴴrₖ₋₁‖")
   kdisplay(iter, verbose) && @printf("%5d  %7.1e  %7s\n", iter, rNorm, "✗ ✗ ✗ ✗")
 
   βₖ = @knrm2(m, r₀)           # β₁ = ‖v₁‖ = ‖r₀‖
@@ -146,7 +146,7 @@ function usymqr!(solver :: UsymqrSolver{T,FC,S}, A, b :: AbstractVector{FC}, c :
   sₖ₋₂ = sₖ₋₁ = sₖ = zero(FC)  # Givens sines used for the QR factorization of Tₖ₊₁.ₖ
   wₖ₋₂ .= zero(FC)             # Column k-2 of Wₖ = Uₖ(Rₖ)⁻¹
   wₖ₋₁ .= zero(FC)             # Column k-1 of Wₖ = Uₖ(Rₖ)⁻¹
-  ζbarₖ = βₖ                   # ζbarₖ is the last component of z̅ₖ = (Qₖ)ᵀβ₁e₁
+  ζbarₖ = βₖ                   # ζbarₖ is the last component of z̅ₖ = (Qₖ)ᴴβ₁e₁
 
   # Stopping criterion.
   solved = rNorm ≤ ε
@@ -161,10 +161,10 @@ function usymqr!(solver :: UsymqrSolver{T,FC,S}, A, b :: AbstractVector{FC}, c :
 
     # Continue the SSY tridiagonalization process.
     # AUₖ  = VₖTₖ    + βₖ₊₁vₖ₊₁(eₖ)ᵀ = Vₖ₊₁Tₖ₊₁.ₖ
-    # AᵀVₖ = Uₖ(Tₖ)ᵀ + γₖ₊₁uₖ₊₁(eₖ)ᵀ = Uₖ₊₁(Tₖ.ₖ₊₁)ᵀ
+    # AᴴVₖ = Uₖ(Tₖ)ᴴ + γₖ₊₁uₖ₊₁(eₖ)ᵀ = Uₖ₊₁(Tₖ.ₖ₊₁)ᴴ
 
     mul!(q, A , uₖ)  # Forms vₖ₊₁ : q ← Auₖ
-    mul!(p, Aᵀ, vₖ)  # Forms uₖ₊₁ : p ← Aᵀvₖ
+    mul!(p, Aᴴ, vₖ)  # Forms uₖ₊₁ : p ← Aᴴvₖ
 
     @kaxpy!(m, -γₖ, vₖ₋₁, q) # q ← q - γₖ * vₖ₋₁
     @kaxpy!(n, -βₖ, uₖ₋₁, p) # p ← p - βₖ * uₖ₋₁
@@ -254,9 +254,9 @@ function usymqr!(solver :: UsymqrSolver{T,FC,S}, A, b :: AbstractVector{FC}, c :
     rNorm = abs(ζbarₖ₊₁)
     history && push!(rNorms, rNorm)
 
-    # Compute ‖Aᵀrₖ₋₁‖ = |ζbarₖ| * √(|δbarₖ|² + |λbarₖ|²).
-    AᵀrNorm = abs(ζbarₖ) * √(abs2(δbarₖ) + abs2(cₖ₋₁ * γₖ₊₁))
-    history && push!(AᵀrNorms, AᵀrNorm)
+    # Compute ‖Aᴴrₖ₋₁‖ = |ζbarₖ| * √(|δbarₖ|² + |λbarₖ|²).
+    AᴴrNorm = abs(ζbarₖ) * √(abs2(δbarₖ) + abs2(cₖ₋₁ * γₖ₊₁))
+    history && push!(AᴴrNorms, AᴴrNorm)
 
     # Compute uₖ₊₁ and uₖ₊₁.
     @. vₖ₋₁ = vₖ # vₖ₋₁ ← vₖ
@@ -286,12 +286,12 @@ function usymqr!(solver :: UsymqrSolver{T,FC,S}, A, b :: AbstractVector{FC}, c :
     βₖ    = βₖ₊₁
 
     # Update stopping criterion.
-    iter == 1 && (κ = atol + rtol * AᵀrNorm)
+    iter == 1 && (κ = atol + rtol * AᴴrNorm)
     user_requested_exit = callback(solver) :: Bool
     solved = rNorm ≤ ε
-    inconsistent = !solved && AᵀrNorm ≤ κ
+    inconsistent = !solved && AᴴrNorm ≤ κ
     tired = iter ≥ itmax
-    kdisplay(iter, verbose) && @printf("%5d  %7.1e  %7.1e\n", iter, rNorm, AᵀrNorm)
+    kdisplay(iter, verbose) && @printf("%5d  %7.1e  %7.1e\n", iter, rNorm, AᴴrNorm)
   end
   (verbose > 0) && @printf("\n")
   tired               && (status = "maximum number of iterations exceeded")