From 3f69bd00a53d764eab91aafd5bb089627272168d Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@gmx.de>
Date: Tue, 30 Jul 2019 18:39:10 +0200
Subject: [PATCH 01/35] add generalized dot product

---
 stdlib/LinearAlgebra/src/bidiag.jl     | 23 ++++++++
 stdlib/LinearAlgebra/src/diagonal.jl   |  5 +-
 stdlib/LinearAlgebra/src/generic.jl    | 51 ++++++++++++++++
 stdlib/LinearAlgebra/src/matmul.jl     |  3 +
 stdlib/LinearAlgebra/src/triangular.jl | 80 ++++++++++++++++++++++++++
 stdlib/LinearAlgebra/src/tridiag.jl    | 28 +++++++++
 stdlib/SparseArrays/src/linalg.jl      | 38 ++++++++++++
 7 files changed, 226 insertions(+), 2 deletions(-)

diff --git a/stdlib/LinearAlgebra/src/bidiag.jl b/stdlib/LinearAlgebra/src/bidiag.jl
index 318c0e5ca6dfd..7eee900d07ab1 100644
--- a/stdlib/LinearAlgebra/src/bidiag.jl
+++ b/stdlib/LinearAlgebra/src/bidiag.jl
@@ -647,6 +647,29 @@ function *(A::SymTridiagonal, B::Diagonal)
     A_mul_B_td!(Tridiagonal(zeros(TS, size(A, 1)-1), zeros(TS, size(A, 1)), zeros(TS, size(A, 1)-1)), A, B)
 end
 
+function dot(x::AbstractVector, B::Bidiagonal, y::AbstractVector)
+    nx, ny = length(x), length(y)
+    (nx == size(B, 1) == ny) || throw(DimensionMismatch())
+    if iszero(nx)
+        return dot(zero(eltype(x)), zero(eltype(B)), zero(eltype(y)))
+    end
+    if A.uplo == 'U'
+        r = adjoint(x[1]) * B.dv[1] * y[1]
+        @inbounds for j in 2:m-1
+            r += (adjoint(x[j-1])*B.ev[j] + adjoint(x[j])*B.dv[j]) * y[j]
+        end
+        r += (adjoint(x[nx-1]) * B.ev[nx-1] + adjoint(x[nx]) * B.dv[nx]) * y[nx]
+        return r
+    else # A.uplo = 'L'
+        r = (adjoint(x[1]) * B.dv[1] + adjoint(x[2]) * B.ev[1]) * y[1]
+        @inbounds for j in 2:m-1
+            r += (adjoint(x[j])*B.dv[j] + adjoint(x[j+1])*B.ev[j]) * y[j]
+        end
+        r += adjoint(x[nx]) * A.d[nx] * y[nx]
+        return r
+    end
+end
+
 #Linear solvers
 ldiv!(A::Union{Bidiagonal, AbstractTriangular}, b::AbstractVector) = naivesub!(A, b)
 ldiv!(A::Transpose{<:Any,<:Bidiagonal}, b::AbstractVector) = ldiv!(copy(A), b)
diff --git a/stdlib/LinearAlgebra/src/diagonal.jl b/stdlib/LinearAlgebra/src/diagonal.jl
index cd17474407d6f..7e1c5142737b7 100644
--- a/stdlib/LinearAlgebra/src/diagonal.jl
+++ b/stdlib/LinearAlgebra/src/diagonal.jl
@@ -637,11 +637,12 @@ end
 
 # disambiguation methods: * of Diagonal and Adj/Trans AbsVec
 *(x::Adjoint{<:Any,<:AbstractVector}, D::Diagonal) = Adjoint(map((t,s) -> t'*s, D.diag, parent(x)))
-*(x::Adjoint{<:Any,<:AbstractVector}, D::Diagonal, y::AbstractVector) =
-    mapreduce(t -> t[1]*t[2]*t[3], +, zip(x, D.diag, y))
 *(x::Transpose{<:Any,<:AbstractVector}, D::Diagonal) = Transpose(map((t,s) -> transpose(t)*s, D.diag, parent(x)))
 *(x::Transpose{<:Any,<:AbstractVector}, D::Diagonal, y::AbstractVector) =
     mapreduce(t -> t[1]*t[2]*t[3], +, zip(x, D.diag, y))
+function dot(x::AbstractVector, D::Diagonal, y::AbstractVector)
+    mapreduce(t -> adjoint(t[1])*t[2]*t[3], +, zip(x, D.diag, y))
+end
 
 function cholesky!(A::Diagonal, ::Val{false} = Val(false); check::Bool = true)
     info = 0
diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl
index b2267463aab45..4e12ea99e9c54 100644
--- a/stdlib/LinearAlgebra/src/generic.jl
+++ b/stdlib/LinearAlgebra/src/generic.jl
@@ -874,6 +874,57 @@ function dot(x::AbstractArray, y::AbstractArray)
     s
 end
 
+"""
+    dot(x, A, y)
+
+Compute the generalized dot product `x' * A * y` between two vectors `x` and `y`.
+For complex vectors, the first vector is conjugated.
+"""
+dot(x::Number, A::Number, y::Number) = conj(x) * A * y
+
+function dot(x::AbstractVector, A::AbstractMatrix, y::AbstractVector)
+    m, n = size(A)
+    (length(x) == m && n == length(y)) || throw(DimensionMismatch())
+    if iszero(m) || iszero(n)
+        return dot(zero(eltype(x)), zero(eltype(A)), zero(eltype(y)))
+    end
+
+    if m == n
+        if istril(A)
+            if istriu(A)
+                return dot(x, Diagonal(A), y)
+            elseif isbanded(A, -1, 0)
+                return dot(x, Bidiagonal(A, 'L'), y)
+            else
+                return dot(x, LowerTriangular(A), y)
+            end
+        end
+        if istriu(A)
+            if isbanded(A, 0, 1)
+                return dot(x, Bidiagonal(A, 'U'), y)
+            else
+                return dot(x, UpperTriangular(A), y)
+            end
+        end
+        if isbanded(A, -1, 1)
+            return dot(x, Tridiagonal(A), y)
+        end
+    end
+
+    T = typeof(dot(first(x), first(A), first(y)))
+    s = zero(T)
+    @inbounds for j in 1:n
+        yj = y[j]
+        if !iszero(yj)
+            temp = zero(T)
+            @simd for i in 1:m
+                temp += adjoint(x[i]) * A[i,j]
+            end
+            s += temp * yj
+        end
+    end
+    return s
+end
 
 ###########################################################################################
 
diff --git a/stdlib/LinearAlgebra/src/matmul.jl b/stdlib/LinearAlgebra/src/matmul.jl
index 0d6ac7a62e1ed..77e5c4f224195 100644
--- a/stdlib/LinearAlgebra/src/matmul.jl
+++ b/stdlib/LinearAlgebra/src/matmul.jl
@@ -3,6 +3,7 @@
 # matmul.jl: Everything to do with dense matrix multiplication
 
 matprod(x, y) = x*y + x*y
+gendot(x, A, y) = adjoint(x)*A*y + adjoint(x)*A*y
 
 # dot products
 
@@ -40,6 +41,8 @@ function *(transx::Transpose{<:Any,<:StridedVector{T}}, y::StridedVector{T}) whe
     return BLAS.dotu(x, y)
 end
 
+(*)(x::Adjoint{<:Any,<:AbstractVector}, A::AbstractMatrix, y::AbstractVector) = dot(parent(x), A, y)
+
 # Matrix-vector multiplication
 function (*)(A::StridedMatrix{T}, x::StridedVector{S}) where {T<:BlasFloat,S<:Real}
     TS = promote_op(matprod, T, S)
diff --git a/stdlib/LinearAlgebra/src/triangular.jl b/stdlib/LinearAlgebra/src/triangular.jl
index 903e8954319ee..4235e19bfbc47 100644
--- a/stdlib/LinearAlgebra/src/triangular.jl
+++ b/stdlib/LinearAlgebra/src/triangular.jl
@@ -545,6 +545,86 @@ end
 rmul!(A::Union{UpperTriangular,LowerTriangular}, c::Number) = mul!(A, A, c)
 lmul!(c::Number, A::Union{UpperTriangular,LowerTriangular}) = mul!(A, c, A)
 
+function dot(x::AbstractVector, A::UpperTriangular, y::AbstractVector)
+    m = size(A, 1)
+    (length(x) == m == length(y)) || throw(DimensionMismatch())
+    if iszero(m)
+        return dot(zero(eltype(x)), zero(eltype(A)), zero(eltype(y)))
+    end
+    T = typeof(dot(first(x), first(A), first(y)))
+    r = zero(T)
+    @inbounds for j in 1:m
+        yj = y[j]
+        if !iszero(yj)
+            temp = zero(T)
+            @simd for i in 1:j
+                temp += adjoint(x[i]) * A[i,j]
+            end
+            r += temp * yj
+        end
+    end
+    return r
+end
+function dot(x::AbstractVector, A::UnitUpperTriangular, y::AbstractVector)
+    m = size(A, 1)
+    (length(x) == m == length(y)) || throw(DimensionMismatch())
+    if iszero(m)
+        return dot(zero(eltype(x)), zero(eltype(A)), zero(eltype(y)))
+    end
+    T = typeof(dot(first(x), first(A), first(y)))
+    r = zero(T)
+    @inbounds for j in 1:m
+        yj = y[j]
+        if !iszero(yj)
+            temp = zero(T)
+            @simd for i in 1:j-1
+                temp += adjoint(x[i]) * A[i,j]
+            end
+            temp += A[j,j]
+            r += temp * yj
+        end
+    end
+    return r
+end
+function dot(x::AbstractVector, A::LowerTriangular, y::AbstractVector)
+    m = size(A, 1)
+    (length(x) == m == length(y)) || throw(DimensionMismatch())
+    if iszero(m)
+        return dot(zero(eltype(x)), zero(eltype(A)), zero(eltype(y)))
+    end
+    T = typeof(dot(first(x), first(A), first(y)))
+    r = zero(T)
+    @inbounds for j in 1:n
+        yj = y[j]
+        if !iszero(yj)
+            temp = zero(T)
+            @simd for i in j:m
+                temp += adjoint(x[i]) * A[i,j]
+            end
+            r += temp * yj
+        end
+    end
+    return r
+end
+function dot(x::AbstractVector, A::UnitLowerTriangular, y::AbstractVector)
+    m = size(A, 1)
+    (length(x) == m == length(y)) || throw(DimensionMismatch())
+
+    T = typeof(dot(first(x), first(A), first(y)))
+    r = zero(T)
+    @inbounds for j in 1:n
+        yj = y[j]
+        if !iszero(yj)
+            temp = convert(T, A[j,j])
+            @simd for i in j+1:m
+                temp += adjoint(x[i]) * A[i,j]
+            end
+            r += temp * yj
+        end
+    end
+    return r
+end
+
 fillstored!(A::LowerTriangular, x)     = (fillband!(A.data, x, 1-size(A,1), 0); A)
 fillstored!(A::UnitLowerTriangular, x) = (fillband!(A.data, x, 1-size(A,1), -1); A)
 fillstored!(A::UpperTriangular, x)     = (fillband!(A.data, x, 0, size(A,2)-1); A)
diff --git a/stdlib/LinearAlgebra/src/tridiag.jl b/stdlib/LinearAlgebra/src/tridiag.jl
index 72d8226f476e9..fac0ee0fb3074 100644
--- a/stdlib/LinearAlgebra/src/tridiag.jl
+++ b/stdlib/LinearAlgebra/src/tridiag.jl
@@ -202,6 +202,20 @@ end
     return C
 end
 
+function dot(x::AbstractVector, S::SymTridiagonal, y::AbstractVector)
+    nx, ny = length(x), length(y)
+    (nx == size(S, 1) == ny) || throw(DimensionMismatch())
+    if iszero(nx)
+        return dot(zero(eltype(x)), zero(eltype(A)), zero(eltype(y)))
+    end
+    r = (adjoint(x[1]) * S.dv[1] + adjoint(x[2]) * adjoint(S.ev[1])) * y[1]
+    @inbounds for j in 2:nx-1
+        r += (adjoint(x[j-1])*S.ev[j-1] + adjoint(x[j])*S.dv[j] + adjoint(x[j+1])*adjoint(S.ev[j])) * yj
+    end
+    r += (adjoint(x[nx-1]) * S.ev[nx-1] + adjoint(x[nx]) * adjoint(S.dv[nx])) * y[nx]
+    return r
+end
+
 (\)(T::SymTridiagonal, B::StridedVecOrMat) = ldlt(T)\B
 
 # division with optional shift for use in shifted-Hessenberg solvers (hessenberg.jl):
@@ -657,3 +671,17 @@ end
 
 Base._sum(A::Tridiagonal, ::Colon) = sum(A.d) + sum(A.dl) + sum(A.du)
 Base._sum(A::SymTridiagonal, ::Colon) = sum(A.dv) + 2sum(A.ev)
+
+function dot(x::AbstractVector, A::Tridiagonal, y::AbstractVector)
+    nx, ny = length(x), length(y)
+    (nx == size(A, 1) == ny) || throw(DimensionMismatch())
+    if iszero(nx)
+        return dot(zero(eltype(x)), zero(eltype(A)), zero(eltype(y)))
+    end
+    r = (adjoint(x[1]) * A.d[1] + adjoint(x[2]) * adjoint(A.dl[1])) * y[1]
+    @inbounds for j in 2:nx-1
+        r += (adjoint(x[j-1])*A.du[j-1] + adjoint(x[j])*A.d[j] + adjoint(x[j+1])*A.dl[j]) * yj
+    end
+    r += (adjoint(x[nx-1]) * A.du[nx-1] + adjoint(x[nx]) * A.d[nx]) * y[nx]
+    return r
+end
diff --git a/stdlib/SparseArrays/src/linalg.jl b/stdlib/SparseArrays/src/linalg.jl
index 437c5a8683eea..a5990e18b76fc 100644
--- a/stdlib/SparseArrays/src/linalg.jl
+++ b/stdlib/SparseArrays/src/linalg.jl
@@ -321,6 +321,44 @@ function dot(A::AbstractSparseMatrixCSC{T1,S1},B::AbstractSparseMatrixCSC{T2,S2}
     return r
 end
 
+function dot(x::AbstractVector, A::SparseMatrixCSC, y::AbstractVector)
+    (length(x) == A.m && A.n == length(y)) || throw(DimensionMismatch())
+    if iszero(A.m) || iszero(A.n)
+        return dot(zero(eltype(x)), zero(eltype(A)), zero(eltype(y)))
+    end
+    T = promote_type(eltype(x), eltype(A), eltype(y))
+    r = zero(T)
+    rvals = rowvals(A)
+    nzvals = nonzeros(A)
+    @inbounds for col in 1:A.n
+        ycol = y[col]
+        if !iszero(ycol)
+            temp = zero(T)
+            for k in nzrange(A, col)
+                temp += x[rvals[k]] * nzvals[k]
+            end
+            r += temp * ycol
+        end
+    end
+    return r
+end
+function dot(x::SparseVector, A::SparseMatrixCSC, y::SparseVector)
+    x.n == A.m && A.n == y.n || throw(DimensionMismatch())
+    if iszero(A.m) || iszero(A.n)
+        return dot(zero(eltype(x)), zero(eltype(A)), zero(eltype(y)))
+    end
+    r = zero(promote_type(eltype(x), eltype(A), eltype(y)))
+    for (yi, yv) in zip(y.nzind, y.nzval)
+        A_ptr_lo = A.colptr[yi]
+        A_ptr_hi = A.colptr[yi+1] - 1
+        if A_ptr_lo <= A_ptr_hi
+            r += _spdot(dot, 1, length(x.nzind), x.nzind, x.nzval,
+                                            A_ptr_lo, A_ptr_hi, A.rowval, A.nzval) * yv
+        end
+    end
+    r
+end
+
 ## triangular sparse handling
 
 possible_adjoint(adj::Bool, a::Real ) = a

From 271781b514a9cba75f24a19114a7d90625f148fa Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@gmx.de>
Date: Wed, 31 Jul 2019 15:26:45 +0200
Subject: [PATCH 02/35] add generalized dot for Adjoint and Transpose

---
 stdlib/LinearAlgebra/src/generic.jl | 64 +++++++++++++++++++----------
 1 file changed, 43 insertions(+), 21 deletions(-)

diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl
index 4e12ea99e9c54..e1555e2cfdf25 100644
--- a/stdlib/LinearAlgebra/src/generic.jl
+++ b/stdlib/LinearAlgebra/src/generic.jl
@@ -889,38 +889,60 @@ function dot(x::AbstractVector, A::AbstractMatrix, y::AbstractVector)
         return dot(zero(eltype(x)), zero(eltype(A)), zero(eltype(y)))
     end
 
-    if m == n
-        if istril(A)
-            if istriu(A)
-                return dot(x, Diagonal(A), y)
-            elseif isbanded(A, -1, 0)
-                return dot(x, Bidiagonal(A, 'L'), y)
-            else
-                return dot(x, LowerTriangular(A), y)
+    T = typeof(dot(first(x), first(A), first(y)))
+    s = zero(T)
+    @inbounds for j in 1:n
+        yj = y[j]
+        if !iszero(yj)
+            temp = zero(T)
+            @simd for i in 1:m
+                temp += adjoint(x[i]) * A[i,j]
             end
+            s += temp * yj
         end
-        if istriu(A)
-            if isbanded(A, 0, 1)
-                return dot(x, Bidiagonal(A, 'U'), y)
-            else
-                return dot(x, UpperTriangular(A), y)
+    end
+    return s
+end
+function dot(x::AbstractVector, adjA::Adjoint, y::AbstractVector)
+    A = adjA.parent
+    m, n = size(A)
+    (length(x) == n && m == length(y)) || throw(DimensionMismatch())
+    if iszero(m) || iszero(n)
+        return dot(zero(eltype(y)), zero(eltype(A)), zero(eltype(x)))
+    end
+
+    T = typeof(dot(first(y), first(A), first(x)))
+    s = zero(T)
+    @inbounds for j in 1:n
+        xj = x[j]
+        if !iszero(xj)
+            temp = zero(T)
+            @simd for i in 1:m
+                temp += adjoint(y[i]) * A[i,j]
             end
-        end
-        if isbanded(A, -1, 1)
-            return dot(x, Tridiagonal(A), y)
+            s += temp * xj
         end
     end
+    return adjoint(s)
+end
+function dot(x::AbstractVector, transA::Transpose, y::AbstractVector)
+    A = transA.parent
+    m, n = size(A)
+    (length(x) == n && m == length(y)) || throw(DimensionMismatch())
+    if iszero(m) || iszero(n)
+        return dot(zero(eltype(y)), zero(eltype(A)), zero(eltype(x)))
+    end
 
-    T = typeof(dot(first(x), first(A), first(y)))
+    T = typeof(dot(first(y), first(A), first(x)))
     s = zero(T)
     @inbounds for j in 1:n
-        yj = y[j]
-        if !iszero(yj)
+        xj = x[j]
+        if !iszero(xj)
             temp = zero(T)
             @simd for i in 1:m
-                temp += adjoint(x[i]) * A[i,j]
+                temp += transpose(y[i]) * A[i,j]
             end
-            s += temp * yj
+            s += temp * conj(xj)
         end
     end
     return s

From 819c601fc982ecafe883e2132b080bf255abf3c2 Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@gmx.de>
Date: Wed, 31 Jul 2019 15:27:14 +0200
Subject: [PATCH 03/35] add "generalized" dot for UniformScalings

---
 stdlib/LinearAlgebra/src/uniformscaling.jl | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/stdlib/LinearAlgebra/src/uniformscaling.jl b/stdlib/LinearAlgebra/src/uniformscaling.jl
index 7b6c3a58000dc..f5836bbdf7c07 100644
--- a/stdlib/LinearAlgebra/src/uniformscaling.jl
+++ b/stdlib/LinearAlgebra/src/uniformscaling.jl
@@ -400,3 +400,5 @@ Array(s::UniformScaling, dims::Dims{2}) = Matrix(s, dims)
 ## Diagonal construction from UniformScaling
 Diagonal{T}(s::UniformScaling, m::Integer) where {T} = Diagonal{T}(fill(T(s.λ), m))
 Diagonal(s::UniformScaling, m::Integer) = Diagonal{eltype(s)}(s, m)
+
+dot(x::AbstractVector, J::UniformScaling, y::AbstractVector) = J.λ * dot(x, y)

From 4ea06def68da238870e8ac296d43dd3234d1d3f1 Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@gmx.de>
Date: Wed, 31 Jul 2019 15:27:48 +0200
Subject: [PATCH 04/35] fix adjoint/transpose in tridiags

---
 stdlib/LinearAlgebra/src/tridiag.jl | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/stdlib/LinearAlgebra/src/tridiag.jl b/stdlib/LinearAlgebra/src/tridiag.jl
index fac0ee0fb3074..97ca8167ffb63 100644
--- a/stdlib/LinearAlgebra/src/tridiag.jl
+++ b/stdlib/LinearAlgebra/src/tridiag.jl
@@ -208,11 +208,11 @@ function dot(x::AbstractVector, S::SymTridiagonal, y::AbstractVector)
     if iszero(nx)
         return dot(zero(eltype(x)), zero(eltype(A)), zero(eltype(y)))
     end
-    r = (adjoint(x[1]) * S.dv[1] + adjoint(x[2]) * adjoint(S.ev[1])) * y[1]
+    r = (adjoint(x[1]) * S.dv[1] + adjoint(x[2]) * transpose(S.ev[1])) * y[1]
     @inbounds for j in 2:nx-1
-        r += (adjoint(x[j-1])*S.ev[j-1] + adjoint(x[j])*S.dv[j] + adjoint(x[j+1])*adjoint(S.ev[j])) * yj
+        r += (adjoint(x[j-1])*S.ev[j-1] + adjoint(x[j])*S.dv[j] + adjoint(x[j+1])*transpose(S.ev[j])) * yj
     end
-    r += (adjoint(x[nx-1]) * S.ev[nx-1] + adjoint(x[nx]) * adjoint(S.dv[nx])) * y[nx]
+    r += (adjoint(x[nx-1]) * S.ev[nx-1] + adjoint(x[nx]) * transpose(S.dv[nx])) * y[nx]
     return r
 end
 
@@ -678,7 +678,7 @@ function dot(x::AbstractVector, A::Tridiagonal, y::AbstractVector)
     if iszero(nx)
         return dot(zero(eltype(x)), zero(eltype(A)), zero(eltype(y)))
     end
-    r = (adjoint(x[1]) * A.d[1] + adjoint(x[2]) * adjoint(A.dl[1])) * y[1]
+    r = (adjoint(x[1]) * A.d[1] + adjoint(x[2]) * A.dl[1]) * y[1]
     @inbounds for j in 2:nx-1
         r += (adjoint(x[j-1])*A.du[j-1] + adjoint(x[j])*A.d[j] + adjoint(x[j+1])*A.dl[j]) * yj
     end

From 512c122d8a65886f9811a5031c0bd1a5728ecbfb Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@gmx.de>
Date: Wed, 31 Jul 2019 18:31:16 +0200
Subject: [PATCH 05/35] improve generic dot, add tests

---
 stdlib/LinearAlgebra/src/generic.jl  | 60 ++++------------------------
 stdlib/LinearAlgebra/src/matmul.jl   |  1 -
 stdlib/LinearAlgebra/test/generic.jl | 28 +++++++++++++
 3 files changed, 35 insertions(+), 54 deletions(-)

diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl
index e1555e2cfdf25..ce6f486c8811a 100644
--- a/stdlib/LinearAlgebra/src/generic.jl
+++ b/stdlib/LinearAlgebra/src/generic.jl
@@ -883,19 +883,15 @@ For complex vectors, the first vector is conjugated.
 dot(x::Number, A::Number, y::Number) = conj(x) * A * y
 
 function dot(x::AbstractVector, A::AbstractMatrix, y::AbstractVector)
-    m, n = size(A)
-    (length(x) == m && n == length(y)) || throw(DimensionMismatch())
-    if iszero(m) || iszero(n)
-        return dot(zero(eltype(x)), zero(eltype(A)), zero(eltype(y)))
-    end
-
+    (axes(x)..., axes(y)...) == axes(A) || throw(DimensionMismatch())
     T = typeof(dot(first(x), first(A), first(y)))
+    zeroxA = zero(adjoint(first(x))*first(A))
     s = zero(T)
-    @inbounds for j in 1:n
+    @inbounds for j in eachindex(y)
         yj = y[j]
         if !iszero(yj)
-            temp = zero(T)
-            @simd for i in 1:m
+            temp = zeroxA
+            @simd for i in eachindex(x)
                 temp += adjoint(x[i]) * A[i,j]
             end
             s += temp * yj
@@ -903,50 +899,8 @@ function dot(x::AbstractVector, A::AbstractMatrix, y::AbstractVector)
     end
     return s
 end
-function dot(x::AbstractVector, adjA::Adjoint, y::AbstractVector)
-    A = adjA.parent
-    m, n = size(A)
-    (length(x) == n && m == length(y)) || throw(DimensionMismatch())
-    if iszero(m) || iszero(n)
-        return dot(zero(eltype(y)), zero(eltype(A)), zero(eltype(x)))
-    end
-
-    T = typeof(dot(first(y), first(A), first(x)))
-    s = zero(T)
-    @inbounds for j in 1:n
-        xj = x[j]
-        if !iszero(xj)
-            temp = zero(T)
-            @simd for i in 1:m
-                temp += adjoint(y[i]) * A[i,j]
-            end
-            s += temp * xj
-        end
-    end
-    return adjoint(s)
-end
-function dot(x::AbstractVector, transA::Transpose, y::AbstractVector)
-    A = transA.parent
-    m, n = size(A)
-    (length(x) == n && m == length(y)) || throw(DimensionMismatch())
-    if iszero(m) || iszero(n)
-        return dot(zero(eltype(y)), zero(eltype(A)), zero(eltype(x)))
-    end
-
-    T = typeof(dot(first(y), first(A), first(x)))
-    s = zero(T)
-    @inbounds for j in 1:n
-        xj = x[j]
-        if !iszero(xj)
-            temp = zero(T)
-            @simd for i in 1:m
-                temp += transpose(y[i]) * A[i,j]
-            end
-            s += temp * conj(xj)
-        end
-    end
-    return s
-end
+dot(x::AbstractVector, adjA::Adjoint, y::AbstractVector) = adjoint(dot(y, adjA.parent, x))
+dot(x::AbstractVector, transA::Transpose{<:Real}, y::AbstractVector) = adjoint(dot(y, transA.parent, x))
 
 ###########################################################################################
 
diff --git a/stdlib/LinearAlgebra/src/matmul.jl b/stdlib/LinearAlgebra/src/matmul.jl
index 77e5c4f224195..2debe360d7238 100644
--- a/stdlib/LinearAlgebra/src/matmul.jl
+++ b/stdlib/LinearAlgebra/src/matmul.jl
@@ -3,7 +3,6 @@
 # matmul.jl: Everything to do with dense matrix multiplication
 
 matprod(x, y) = x*y + x*y
-gendot(x, A, y) = adjoint(x)*A*y + adjoint(x)*A*y
 
 # dot products
 
diff --git a/stdlib/LinearAlgebra/test/generic.jl b/stdlib/LinearAlgebra/test/generic.jl
index 8880325d21e03..37dc3906925aa 100644
--- a/stdlib/LinearAlgebra/test/generic.jl
+++ b/stdlib/LinearAlgebra/test/generic.jl
@@ -409,4 +409,32 @@ end
     @test all(!isnan, lmul!(false, Any[NaN]))
 end
 
+@testset "generalized dot #32739" begin
+    for elty in (Int, Float32, Float64, BigFloat, Complex{Float32}, Complex{Float64}, Complex{BigFloat})
+        n = 10
+        if elty <: Int
+            A = rand(-n:n, n, n)
+            x = rand(-n:n, n)
+            y = rand(-n:n, n)
+        elseif elty <: Real
+            A = convert(Matrix{elty}, randn(n,n))
+            x = rand(elty, n)
+            y = rand(elty, n)
+        else
+            A = convert(Matrix{elty}, complex.(randn(n,n), randn(n,n)))
+            x = rand(elty, n)
+            y = rand(elty, n)
+        end
+        @test dot(x, A, y) ≈ x' * A * y ≈ *(x', A, y) ≈ (x'A)*y
+        @test dot(x, A', y) ≈ x' * A' * y ≈ *(x', A', y) ≈ (x'A')*y
+        elty <: Real && @test dot(x, transpose(A), y) ≈ x' * transpose(A) * y ≈ *(x', transpose(A), y) ≈ (x'*transpose(A))*y
+        B = reshape([A], 1, 1)
+        x = [x]
+        y = [y]
+        @test dot(x, B, y) ≈ x' * B * y ≈ *(x', B, y) ≈ (x'B)*y
+        @test dot(x, B', y) ≈ x' * B' * y ≈ *(x', B', y) ≈ (x'B')*y
+        elty <: Real && @test dot(x, transpose(B), y) ≈ x' * transpose(B) * y ≈ *(x', transpose(B), y) ≈ (x'*transpose(B))*y
+    end
+end
+
 end # module TestGeneric

From b068ed1afe120881a6c08539ab4f7dc3c671b94f Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@gmx.de>
Date: Thu, 1 Aug 2019 10:56:20 +0200
Subject: [PATCH 06/35] fix typos, optimize *diag, require_one_based_indexing

---
 stdlib/LinearAlgebra/src/bidiag.jl     | 27 ++++++++++++++++----------
 stdlib/LinearAlgebra/src/triangular.jl |  4 ++++
 stdlib/LinearAlgebra/src/tridiag.jl    | 26 ++++++++++++++++++-------
 stdlib/SparseArrays/src/linalg.jl      |  1 +
 4 files changed, 41 insertions(+), 17 deletions(-)

diff --git a/stdlib/LinearAlgebra/src/bidiag.jl b/stdlib/LinearAlgebra/src/bidiag.jl
index 7eee900d07ab1..e108796f9aa43 100644
--- a/stdlib/LinearAlgebra/src/bidiag.jl
+++ b/stdlib/LinearAlgebra/src/bidiag.jl
@@ -648,24 +648,31 @@ function *(A::SymTridiagonal, B::Diagonal)
 end
 
 function dot(x::AbstractVector, B::Bidiagonal, y::AbstractVector)
+    require_one_based_indexing(x, y)
     nx, ny = length(x), length(y)
     (nx == size(B, 1) == ny) || throw(DimensionMismatch())
     if iszero(nx)
         return dot(zero(eltype(x)), zero(eltype(B)), zero(eltype(y)))
     end
-    if A.uplo == 'U'
-        r = adjoint(x[1]) * B.dv[1] * y[1]
-        @inbounds for j in 2:m-1
-            r += (adjoint(x[j-1])*B.ev[j] + adjoint(x[j])*B.dv[j]) * y[j]
+    ev, dv = B.ev, B.dv
+    if B.uplo == 'U'
+        x₀ = adjoint(x[1])
+        r = x₀ * dv[1] * y[1]
+        @inbounds for j in 2:nx-1
+            x₋, x₀ = x₀, adjoint(x[j])
+            r += (x₋*ev[j-1] + x₀*dv[j]) * y[j]
         end
-        r += (adjoint(x[nx-1]) * B.ev[nx-1] + adjoint(x[nx]) * B.dv[nx]) * y[nx]
+        r += (x₀*ev[nx-1] + adjoint(x[nx])*dv[nx]) * y[nx]
         return r
-    else # A.uplo = 'L'
-        r = (adjoint(x[1]) * B.dv[1] + adjoint(x[2]) * B.ev[1]) * y[1]
-        @inbounds for j in 2:m-1
-            r += (adjoint(x[j])*B.dv[j] + adjoint(x[j+1])*B.ev[j]) * y[j]
+    else # B.uplo == 'L'
+        x₀ = adjoint(x[1])
+        x₊ = adjoint(x[2])
+        r = (x₀*dv[1] + x₊*ev[1]) * y[1]
+        @inbounds for j in 2:nx-1
+            x₀, x₊ = x₊, adjoint(x[j+1])
+            r += (x₀*dv[j] + x₊*ev[j]) * y[j]
         end
-        r += adjoint(x[nx]) * A.d[nx] * y[nx]
+        r += x₊ * dv[nx] * y[nx]
         return r
     end
 end
diff --git a/stdlib/LinearAlgebra/src/triangular.jl b/stdlib/LinearAlgebra/src/triangular.jl
index 4235e19bfbc47..d83fb414e5802 100644
--- a/stdlib/LinearAlgebra/src/triangular.jl
+++ b/stdlib/LinearAlgebra/src/triangular.jl
@@ -546,6 +546,7 @@ rmul!(A::Union{UpperTriangular,LowerTriangular}, c::Number) = mul!(A, A, c)
 lmul!(c::Number, A::Union{UpperTriangular,LowerTriangular}) = mul!(A, c, A)
 
 function dot(x::AbstractVector, A::UpperTriangular, y::AbstractVector)
+    require_one_based_indexing(x, y)
     m = size(A, 1)
     (length(x) == m == length(y)) || throw(DimensionMismatch())
     if iszero(m)
@@ -566,6 +567,7 @@ function dot(x::AbstractVector, A::UpperTriangular, y::AbstractVector)
     return r
 end
 function dot(x::AbstractVector, A::UnitUpperTriangular, y::AbstractVector)
+    require_one_based_indexing(x, y)
     m = size(A, 1)
     (length(x) == m == length(y)) || throw(DimensionMismatch())
     if iszero(m)
@@ -587,6 +589,7 @@ function dot(x::AbstractVector, A::UnitUpperTriangular, y::AbstractVector)
     return r
 end
 function dot(x::AbstractVector, A::LowerTriangular, y::AbstractVector)
+    require_one_based_indexing(x, y)
     m = size(A, 1)
     (length(x) == m == length(y)) || throw(DimensionMismatch())
     if iszero(m)
@@ -607,6 +610,7 @@ function dot(x::AbstractVector, A::LowerTriangular, y::AbstractVector)
     return r
 end
 function dot(x::AbstractVector, A::UnitLowerTriangular, y::AbstractVector)
+    require_one_based_indexing(x, y)
     m = size(A, 1)
     (length(x) == m == length(y)) || throw(DimensionMismatch())
 
diff --git a/stdlib/LinearAlgebra/src/tridiag.jl b/stdlib/LinearAlgebra/src/tridiag.jl
index 97ca8167ffb63..356e6c532ec31 100644
--- a/stdlib/LinearAlgebra/src/tridiag.jl
+++ b/stdlib/LinearAlgebra/src/tridiag.jl
@@ -203,16 +203,23 @@ end
 end
 
 function dot(x::AbstractVector, S::SymTridiagonal, y::AbstractVector)
+    require_one_based_indexing(x, y)
     nx, ny = length(x), length(y)
     (nx == size(S, 1) == ny) || throw(DimensionMismatch())
     if iszero(nx)
-        return dot(zero(eltype(x)), zero(eltype(A)), zero(eltype(y)))
+        return dot(zero(eltype(x)), zero(eltype(S)), zero(eltype(y)))
     end
-    r = (adjoint(x[1]) * S.dv[1] + adjoint(x[2]) * transpose(S.ev[1])) * y[1]
+    dv, ev = S.dv, S.ev
+    x₀ = adjoint(x[1])
+    x₊ = adjoint(x[2])
+    sub = transpose(ev[1])
+    r = (x₀*dv[1] + x₊*sub) * y[1]
     @inbounds for j in 2:nx-1
-        r += (adjoint(x[j-1])*S.ev[j-1] + adjoint(x[j])*S.dv[j] + adjoint(x[j+1])*transpose(S.ev[j])) * yj
+        x₋, x₀, x₊ = x₀, x₊, adjoint(x[j+1])
+        sup, sub = transpose(sub), transpose(ev[j])
+        r += (x₋*sup + x₀*dv[j] + x₊*sub) * y[j]
     end
-    r += (adjoint(x[nx-1]) * S.ev[nx-1] + adjoint(x[nx]) * transpose(S.dv[nx])) * y[nx]
+    r += (x₀*transpose(sub) + x_₊*dv[nx]) * y[nx]
     return r
 end
 
@@ -673,15 +680,20 @@ Base._sum(A::Tridiagonal, ::Colon) = sum(A.d) + sum(A.dl) + sum(A.du)
 Base._sum(A::SymTridiagonal, ::Colon) = sum(A.dv) + 2sum(A.ev)
 
 function dot(x::AbstractVector, A::Tridiagonal, y::AbstractVector)
+    require_one_based_indexing(x, y)
     nx, ny = length(x), length(y)
     (nx == size(A, 1) == ny) || throw(DimensionMismatch())
     if iszero(nx)
         return dot(zero(eltype(x)), zero(eltype(A)), zero(eltype(y)))
     end
-    r = (adjoint(x[1]) * A.d[1] + adjoint(x[2]) * A.dl[1]) * y[1]
+    x₀ = adjoint(x[1])
+    x₊ = adjoint(x[2])
+    dl, d, du = A.dl, A.d, A.du
+    r = (x₀*d[1] + x₊*dl[1]) * y[1]
     @inbounds for j in 2:nx-1
-        r += (adjoint(x[j-1])*A.du[j-1] + adjoint(x[j])*A.d[j] + adjoint(x[j+1])*A.dl[j]) * yj
+        x₋, x₀, x₊ = x₀, x₊, adjoint(x[j+1])
+        r += (x₋*du[j-1] + x₀*d[j] + x₊*dl[j]) * y[j]
     end
-    r += (adjoint(x[nx-1]) * A.du[nx-1] + adjoint(x[nx]) * A.d[nx]) * y[nx]
+    r += (x₀*du[nx-1] + x₊*d[nx]) * y[nx]
     return r
 end
diff --git a/stdlib/SparseArrays/src/linalg.jl b/stdlib/SparseArrays/src/linalg.jl
index a5990e18b76fc..1a8d961a2901e 100644
--- a/stdlib/SparseArrays/src/linalg.jl
+++ b/stdlib/SparseArrays/src/linalg.jl
@@ -322,6 +322,7 @@ function dot(A::AbstractSparseMatrixCSC{T1,S1},B::AbstractSparseMatrixCSC{T2,S2}
 end
 
 function dot(x::AbstractVector, A::SparseMatrixCSC, y::AbstractVector)
+    require_one_based_indexing(x, y)
     (length(x) == A.m && A.n == length(y)) || throw(DimensionMismatch())
     if iszero(A.m) || iszero(A.n)
         return dot(zero(eltype(x)), zero(eltype(A)), zero(eltype(y)))

From 96b3a8bf793a7506c7c8b9bbbc265c1140f31bda Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@gmx.de>
Date: Thu, 1 Aug 2019 10:56:35 +0200
Subject: [PATCH 07/35] add tests

---
 stdlib/LinearAlgebra/test/bidiag.jl         | 13 +++++++++++++
 stdlib/LinearAlgebra/test/triangular.jl     |  7 +++++++
 stdlib/LinearAlgebra/test/tridiag.jl        |  5 +++++
 stdlib/LinearAlgebra/test/uniformscaling.jl |  8 ++++++++
 4 files changed, 33 insertions(+)

diff --git a/stdlib/LinearAlgebra/test/bidiag.jl b/stdlib/LinearAlgebra/test/bidiag.jl
index 166449018d1f0..b0578687d5271 100644
--- a/stdlib/LinearAlgebra/test/bidiag.jl
+++ b/stdlib/LinearAlgebra/test/bidiag.jl
@@ -455,4 +455,17 @@ end
     @test A * Tridiagonal(ones(1, 1)) == A
 end
 
+@testset "generalized dot" begin
+    for elty in (Float64, ComplexF64)
+        dv = randn(elty, 5)
+        ev = randn(elty, 4)
+        x = randn(elty, 5)
+        y = randn(elty, 5)
+        for uplo in (:U, :L)
+            B = Bidiagonal(dv, ev, uplo)
+            @test dot(x, B, y) ≈ x' * B * y ≈ *(x', B, y) ≈ (x'B)*y
+        end
+    end
+end
+
 end # module TestBidiagonal
diff --git a/stdlib/LinearAlgebra/test/triangular.jl b/stdlib/LinearAlgebra/test/triangular.jl
index 612fe05374dcd..94492c7e4b898 100644
--- a/stdlib/LinearAlgebra/test/triangular.jl
+++ b/stdlib/LinearAlgebra/test/triangular.jl
@@ -236,6 +236,13 @@ for elty1 in (Float32, Float64, BigFloat, ComplexF32, ComplexF64, Complex{BigFlo
             end
         end
 
+        # generalized dot
+        for eltyb in (Float32, Float64, BigFloat, ComplexF32, ComplexF64, Complex{BigFloat})
+            b1 = convert(Vector{eltyb}, (elty1 <: Complex ? real(A1) : A1)*fill(1., n))
+            b2 = convert(Vector{eltyb}, (elty1 <: Complex ? real(A1) : A1)*randn(n))
+            @test dot(b1, A1, b2) ≈ b1' * A1 * b2 ≈ *(b1', A1, b2) ≈ (b1'A1)*b2
+        end
+
         # Binary operations
         @test A1*0.5 == Matrix(A1)*0.5
         @test 0.5*A1 == 0.5*Matrix(A1)
diff --git a/stdlib/LinearAlgebra/test/tridiag.jl b/stdlib/LinearAlgebra/test/tridiag.jl
index 686a1010914ab..bc3eb4c63e99e 100644
--- a/stdlib/LinearAlgebra/test/tridiag.jl
+++ b/stdlib/LinearAlgebra/test/tridiag.jl
@@ -390,6 +390,11 @@ end
                 end
             end
         end
+        @testset "generalized dot" begin
+            x = fill(convert(elty, 1), n)
+            y = fill(convert(elty, 1), n)
+            @test dot(x, A, y) ≈ x' * A * y ≈ *(x', A, y) ≈ (x'A)*y
+        end
     end
 end
 
diff --git a/stdlib/LinearAlgebra/test/uniformscaling.jl b/stdlib/LinearAlgebra/test/uniformscaling.jl
index 21a8cf0bb542c..fd0c80e90326d 100644
--- a/stdlib/LinearAlgebra/test/uniformscaling.jl
+++ b/stdlib/LinearAlgebra/test/uniformscaling.jl
@@ -330,4 +330,12 @@ end
     @test I(3) == [1 0 0; 0 1 0; 0 0 1]
 end
 
+@testset "generalized dot" begin
+    x = rand(-10:10, 3)
+    y = rand(-10:10, 3)
+    λ = rand(-10:10)
+    J = UniformScaling(λ)
+    @test dot(x, J, y) == λ*dot(x, y)
+end
+
 end # module TestUniformscaling

From 6012d59a1ee9afb05430547dbd1565be5db7ca62 Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@gmx.de>
Date: Thu, 1 Aug 2019 12:04:34 +0200
Subject: [PATCH 08/35] fix typos in triangular and tridiag

---
 stdlib/LinearAlgebra/src/triangular.jl | 12 +++++++-----
 stdlib/LinearAlgebra/src/tridiag.jl    |  2 +-
 2 files changed, 8 insertions(+), 6 deletions(-)

diff --git a/stdlib/LinearAlgebra/src/triangular.jl b/stdlib/LinearAlgebra/src/triangular.jl
index d83fb414e5802..65934c10dbba5 100644
--- a/stdlib/LinearAlgebra/src/triangular.jl
+++ b/stdlib/LinearAlgebra/src/triangular.jl
@@ -582,7 +582,7 @@ function dot(x::AbstractVector, A::UnitUpperTriangular, y::AbstractVector)
             @simd for i in 1:j-1
                 temp += adjoint(x[i]) * A[i,j]
             end
-            temp += A[j,j]
+            temp += adjoint(x[j])
             r += temp * yj
         end
     end
@@ -597,7 +597,7 @@ function dot(x::AbstractVector, A::LowerTriangular, y::AbstractVector)
     end
     T = typeof(dot(first(x), first(A), first(y)))
     r = zero(T)
-    @inbounds for j in 1:n
+    @inbounds for j in 1:m
         yj = y[j]
         if !iszero(yj)
             temp = zero(T)
@@ -613,13 +613,15 @@ function dot(x::AbstractVector, A::UnitLowerTriangular, y::AbstractVector)
     require_one_based_indexing(x, y)
     m = size(A, 1)
     (length(x) == m == length(y)) || throw(DimensionMismatch())
-
+    if iszero(m)
+        return dot(zero(eltype(x)), zero(eltype(A)), zero(eltype(y)))
+    end
     T = typeof(dot(first(x), first(A), first(y)))
     r = zero(T)
-    @inbounds for j in 1:n
+    @inbounds for j in 1:m
         yj = y[j]
         if !iszero(yj)
-            temp = convert(T, A[j,j])
+            temp = adjoint(x[j])
             @simd for i in j+1:m
                 temp += adjoint(x[i]) * A[i,j]
             end
diff --git a/stdlib/LinearAlgebra/src/tridiag.jl b/stdlib/LinearAlgebra/src/tridiag.jl
index 356e6c532ec31..6fc1d3625dd95 100644
--- a/stdlib/LinearAlgebra/src/tridiag.jl
+++ b/stdlib/LinearAlgebra/src/tridiag.jl
@@ -219,7 +219,7 @@ function dot(x::AbstractVector, S::SymTridiagonal, y::AbstractVector)
         sup, sub = transpose(sub), transpose(ev[j])
         r += (x₋*sup + x₀*dv[j] + x₊*sub) * y[j]
     end
-    r += (x₀*transpose(sub) + x_₊*dv[nx]) * y[nx]
+    r += (x₀*transpose(sub) + x₊*dv[nx]) * y[nx]
     return r
 end
 

From 430320b982d9000044d43df8acc0205017cf6eb9 Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@gmx.de>
Date: Thu, 1 Aug 2019 13:06:27 +0200
Subject: [PATCH 09/35] fix BigFloat tests in triangular

---
 stdlib/LinearAlgebra/test/triangular.jl | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/stdlib/LinearAlgebra/test/triangular.jl b/stdlib/LinearAlgebra/test/triangular.jl
index 94492c7e4b898..d207d3ea1cdbf 100644
--- a/stdlib/LinearAlgebra/test/triangular.jl
+++ b/stdlib/LinearAlgebra/test/triangular.jl
@@ -240,7 +240,11 @@ for elty1 in (Float32, Float64, BigFloat, ComplexF32, ComplexF64, Complex{BigFlo
         for eltyb in (Float32, Float64, BigFloat, ComplexF32, ComplexF64, Complex{BigFloat})
             b1 = convert(Vector{eltyb}, (elty1 <: Complex ? real(A1) : A1)*fill(1., n))
             b2 = convert(Vector{eltyb}, (elty1 <: Complex ? real(A1) : A1)*randn(n))
-            @test dot(b1, A1, b2) ≈ b1' * A1 * b2 ≈ *(b1', A1, b2) ≈ (b1'A1)*b2
+            if elty1 in (BigFloat, Complex{BigFloat}) || eltyb in (BigFloat, Complex{BigFloat})
+                @test dot(b1, A1, b2) ≈ (b1'A1)*b2  atol=sqrt(max(eps(real(float(one(elty1)))),eps(real(float(one(eltyb))))))*n*n
+            else
+                @test dot(b1, A1, b2) ≈ (b1'A1)*b2
+            end
         end
 
         # Binary operations

From 5008546028ab95e7555672a4e02495d13194ebd2 Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@gmx.de>
Date: Thu, 1 Aug 2019 17:11:52 +0200
Subject: [PATCH 10/35] add sparse tests (and minor fix)

---
 stdlib/SparseArrays/src/linalg.jl  | 2 +-
 stdlib/SparseArrays/test/sparse.jl | 9 +++++++++
 2 files changed, 10 insertions(+), 1 deletion(-)

diff --git a/stdlib/SparseArrays/src/linalg.jl b/stdlib/SparseArrays/src/linalg.jl
index 1a8d961a2901e..9393bfab17358 100644
--- a/stdlib/SparseArrays/src/linalg.jl
+++ b/stdlib/SparseArrays/src/linalg.jl
@@ -336,7 +336,7 @@ function dot(x::AbstractVector, A::SparseMatrixCSC, y::AbstractVector)
         if !iszero(ycol)
             temp = zero(T)
             for k in nzrange(A, col)
-                temp += x[rvals[k]] * nzvals[k]
+                temp += adjoint(x[rvals[k]]) * nzvals[k]
             end
             r += temp * ycol
         end
diff --git a/stdlib/SparseArrays/test/sparse.jl b/stdlib/SparseArrays/test/sparse.jl
index bd85f28609c44..45099f399fb51 100644
--- a/stdlib/SparseArrays/test/sparse.jl
+++ b/stdlib/SparseArrays/test/sparse.jl
@@ -424,6 +424,15 @@ end
     @test_throws DimensionMismatch dot(sprand(5,5,0.2),sprand(5,6,0.2))
 end
 
+@testset "generalized dot product" begin
+    for i = 1:5
+        A = sprand(ComplexF64, 10, 15, 0.4)
+        x = sprand(ComplexF64, 10, 0.5)
+        y = sprand(ComplexF64, 15, 0.5)
+        @test dot(x, A, y) ≈ dot(Vector(x), A, Vector(y)) ≈ (Vector(x)' * Matrix(A)) * Vector(y)
+    end
+end
+
 const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test")
 isdefined(Main, :Quaternions) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "Quaternions.jl"))
 using .Main.Quaternions

From a28280fe60b613492c6f76ef1b821a1b6e45ba0b Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@gmx.de>
Date: Tue, 6 Aug 2019 13:14:12 +0200
Subject: [PATCH 11/35] handle block arrays of varying lengths

---
 stdlib/LinearAlgebra/src/generic.jl | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl
index ce6f486c8811a..aeb25f830c47a 100644
--- a/stdlib/LinearAlgebra/src/generic.jl
+++ b/stdlib/LinearAlgebra/src/generic.jl
@@ -885,12 +885,11 @@ dot(x::Number, A::Number, y::Number) = conj(x) * A * y
 function dot(x::AbstractVector, A::AbstractMatrix, y::AbstractVector)
     (axes(x)..., axes(y)...) == axes(A) || throw(DimensionMismatch())
     T = typeof(dot(first(x), first(A), first(y)))
-    zeroxA = zero(adjoint(first(x))*first(A))
     s = zero(T)
     @inbounds for j in eachindex(y)
         yj = y[j]
         if !iszero(yj)
-            temp = zeroxA
+            temp = zero(adjoint(yj))
             @simd for i in eachindex(x)
                 temp += adjoint(x[i]) * A[i,j]
             end

From 93bb8c77877d9e2b31a692225363e8682dd63402 Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@gmx.de>
Date: Tue, 6 Aug 2019 13:52:06 +0200
Subject: [PATCH 12/35] make generalized dot act recursively

---
 stdlib/LinearAlgebra/src/bidiag.jl     | 12 ++++----
 stdlib/LinearAlgebra/src/diagonal.jl   |  4 ++-
 stdlib/LinearAlgebra/src/generic.jl    |  6 ++--
 stdlib/LinearAlgebra/src/matmul.jl     |  3 +-
 stdlib/LinearAlgebra/src/triangular.jl | 42 ++++++++++++--------------
 stdlib/LinearAlgebra/src/tridiag.jl    | 12 ++++----
 6 files changed, 41 insertions(+), 38 deletions(-)

diff --git a/stdlib/LinearAlgebra/src/bidiag.jl b/stdlib/LinearAlgebra/src/bidiag.jl
index e108796f9aa43..758e5cbd1dd02 100644
--- a/stdlib/LinearAlgebra/src/bidiag.jl
+++ b/stdlib/LinearAlgebra/src/bidiag.jl
@@ -657,22 +657,22 @@ function dot(x::AbstractVector, B::Bidiagonal, y::AbstractVector)
     ev, dv = B.ev, B.dv
     if B.uplo == 'U'
         x₀ = adjoint(x[1])
-        r = x₀ * dv[1] * y[1]
+        r = dot(x[1], dv[1], y[1])
         @inbounds for j in 2:nx-1
             x₋, x₀ = x₀, adjoint(x[j])
-            r += (x₋*ev[j-1] + x₀*dv[j]) * y[j]
+            r += dot(adjoint(x₋*ev[j-1] + x₀*dv[j]), y[j])
         end
-        r += (x₀*ev[nx-1] + adjoint(x[nx])*dv[nx]) * y[nx]
+        r += dot(adjoint(x₀*ev[nx-1] + adjoint(x[nx])*dv[nx]), y[nx])
         return r
     else # B.uplo == 'L'
         x₀ = adjoint(x[1])
         x₊ = adjoint(x[2])
-        r = (x₀*dv[1] + x₊*ev[1]) * y[1]
+        r = dot(adjoint(x₀*dv[1] + x₊*ev[1]), y[1])
         @inbounds for j in 2:nx-1
             x₀, x₊ = x₊, adjoint(x[j+1])
-            r += (x₀*dv[j] + x₊*ev[j]) * y[j]
+            r += dot(adjoint(x₀*dv[j] + x₊*ev[j]), y[j])
         end
-        r += x₊ * dv[nx] * y[nx]
+        r += dot(adjoint(x₊), dv[nx], y[nx])
         return r
     end
 end
diff --git a/stdlib/LinearAlgebra/src/diagonal.jl b/stdlib/LinearAlgebra/src/diagonal.jl
index 7e1c5142737b7..bc531d9ac4801 100644
--- a/stdlib/LinearAlgebra/src/diagonal.jl
+++ b/stdlib/LinearAlgebra/src/diagonal.jl
@@ -638,10 +638,12 @@ end
 # disambiguation methods: * of Diagonal and Adj/Trans AbsVec
 *(x::Adjoint{<:Any,<:AbstractVector}, D::Diagonal) = Adjoint(map((t,s) -> t'*s, D.diag, parent(x)))
 *(x::Transpose{<:Any,<:AbstractVector}, D::Diagonal) = Transpose(map((t,s) -> transpose(t)*s, D.diag, parent(x)))
+*(x::Adjoint{<:Any,<:AbstractVector}, D::Diagonal, y::AbstractVector) =
+    mapreduce(t -> t[1]*t[2]*t[3], +, zip(x, D.diag, y))
 *(x::Transpose{<:Any,<:AbstractVector}, D::Diagonal, y::AbstractVector) =
     mapreduce(t -> t[1]*t[2]*t[3], +, zip(x, D.diag, y))
 function dot(x::AbstractVector, D::Diagonal, y::AbstractVector)
-    mapreduce(t -> adjoint(t[1])*t[2]*t[3], +, zip(x, D.diag, y))
+    mapreduce(t -> dot(t[1], t[2], t[3]), +, zip(x, D.diag, y))
 end
 
 function cholesky!(A::Diagonal, ::Val{false} = Val(false); check::Bool = true)
diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl
index aeb25f830c47a..880a75c5745ed 100644
--- a/stdlib/LinearAlgebra/src/generic.jl
+++ b/stdlib/LinearAlgebra/src/generic.jl
@@ -886,14 +886,16 @@ function dot(x::AbstractVector, A::AbstractMatrix, y::AbstractVector)
     (axes(x)..., axes(y)...) == axes(A) || throw(DimensionMismatch())
     T = typeof(dot(first(x), first(A), first(y)))
     s = zero(T)
+    i₁ = first(eachindex(x))
+    x₁ = first(x)
     @inbounds for j in eachindex(y)
         yj = y[j]
         if !iszero(yj)
-            temp = zero(adjoint(yj))
+            temp = zero(adjoint(x₁) * A[i₁,j])
             @simd for i in eachindex(x)
                 temp += adjoint(x[i]) * A[i,j]
             end
-            s += temp * yj
+            s += dot(adjoint(temp), yj)
         end
     end
     return s
diff --git a/stdlib/LinearAlgebra/src/matmul.jl b/stdlib/LinearAlgebra/src/matmul.jl
index 2debe360d7238..1bafd979ea2f1 100644
--- a/stdlib/LinearAlgebra/src/matmul.jl
+++ b/stdlib/LinearAlgebra/src/matmul.jl
@@ -40,7 +40,8 @@ function *(transx::Transpose{<:Any,<:StridedVector{T}}, y::StridedVector{T}) whe
     return BLAS.dotu(x, y)
 end
 
-(*)(x::Adjoint{<:Any,<:AbstractVector}, A::AbstractMatrix, y::AbstractVector) = dot(parent(x), A, y)
+# I'm no longer sure we want ternary multiplication to be generally interpreted recursively
+# (*)(x::Adjoint{<:Any,<:AbstractVector}, A::AbstractMatrix, y::AbstractVector) = dot(parent(x), A, y)
 
 # Matrix-vector multiplication
 function (*)(A::StridedMatrix{T}, x::StridedVector{S}) where {T<:BlasFloat,S<:Real}
diff --git a/stdlib/LinearAlgebra/src/triangular.jl b/stdlib/LinearAlgebra/src/triangular.jl
index 65934c10dbba5..1d3b16d7544de 100644
--- a/stdlib/LinearAlgebra/src/triangular.jl
+++ b/stdlib/LinearAlgebra/src/triangular.jl
@@ -552,16 +552,16 @@ function dot(x::AbstractVector, A::UpperTriangular, y::AbstractVector)
     if iszero(m)
         return dot(zero(eltype(x)), zero(eltype(A)), zero(eltype(y)))
     end
-    T = typeof(dot(first(x), first(A), first(y)))
-    r = zero(T)
+    x₁ = first(x)
+    r = zero(typeof(dot(x₁, first(A), first(y))))
     @inbounds for j in 1:m
         yj = y[j]
         if !iszero(yj)
-            temp = zero(T)
-            @simd for i in 1:j
+            temp = adjoint(x₁) * A[1,j]
+            @simd for i in 2:j
                 temp += adjoint(x[i]) * A[i,j]
             end
-            r += temp * yj
+            r += dot(adjoint(temp), yj)
         end
     end
     return r
@@ -573,17 +573,17 @@ function dot(x::AbstractVector, A::UnitUpperTriangular, y::AbstractVector)
     if iszero(m)
         return dot(zero(eltype(x)), zero(eltype(A)), zero(eltype(y)))
     end
-    T = typeof(dot(first(x), first(A), first(y)))
-    r = zero(T)
-    @inbounds for j in 1:m
+    x₁ = first(x)
+    r = dot(x₁, y[1])
+    @inbounds for j in 2:m
         yj = y[j]
         if !iszero(yj)
-            temp = zero(T)
-            @simd for i in 1:j-1
+            temp = adjoint(x₁) * A[1,j]
+            @simd for i in 2:j-1
                 temp += adjoint(x[i]) * A[i,j]
             end
-            temp += adjoint(x[j])
-            r += temp * yj
+            r += dot(adjoint(temp), yj)
+            r += dot(x[j], yj)
         end
     end
     return r
@@ -595,16 +595,15 @@ function dot(x::AbstractVector, A::LowerTriangular, y::AbstractVector)
     if iszero(m)
         return dot(zero(eltype(x)), zero(eltype(A)), zero(eltype(y)))
     end
-    T = typeof(dot(first(x), first(A), first(y)))
-    r = zero(T)
+    r = zero(typeof(dot(first(x), first(A), first(y))))
     @inbounds for j in 1:m
         yj = y[j]
         if !iszero(yj)
-            temp = zero(T)
-            @simd for i in j:m
+            temp = adjoint(x[j]) * A[j,j]
+            @simd for i in j+1:m
                 temp += adjoint(x[i]) * A[i,j]
             end
-            r += temp * yj
+            r += dot(adjoint(temp), yj)
         end
     end
     return r
@@ -616,16 +615,15 @@ function dot(x::AbstractVector, A::UnitLowerTriangular, y::AbstractVector)
     if iszero(m)
         return dot(zero(eltype(x)), zero(eltype(A)), zero(eltype(y)))
     end
-    T = typeof(dot(first(x), first(A), first(y)))
-    r = zero(T)
-    @inbounds for j in 1:m
+    r = dot(x[1], y[1])
+    @inbounds for j in 2:m
         yj = y[j]
         if !iszero(yj)
-            temp = adjoint(x[j])
+            temp = dot(x[j], yj)
             @simd for i in j+1:m
                 temp += adjoint(x[i]) * A[i,j]
             end
-            r += temp * yj
+            r += dot(adjoint(temp), yj)
         end
     end
     return r
diff --git a/stdlib/LinearAlgebra/src/tridiag.jl b/stdlib/LinearAlgebra/src/tridiag.jl
index 6fc1d3625dd95..413a2b0fe9a92 100644
--- a/stdlib/LinearAlgebra/src/tridiag.jl
+++ b/stdlib/LinearAlgebra/src/tridiag.jl
@@ -213,13 +213,13 @@ function dot(x::AbstractVector, S::SymTridiagonal, y::AbstractVector)
     x₀ = adjoint(x[1])
     x₊ = adjoint(x[2])
     sub = transpose(ev[1])
-    r = (x₀*dv[1] + x₊*sub) * y[1]
+    r = dot(adjoint(x₀*dv[1] + x₊*sub), y[1])
     @inbounds for j in 2:nx-1
         x₋, x₀, x₊ = x₀, x₊, adjoint(x[j+1])
         sup, sub = transpose(sub), transpose(ev[j])
-        r += (x₋*sup + x₀*dv[j] + x₊*sub) * y[j]
+        r += dot(adjoint(x₋*sup + x₀*dv[j] + x₊*sub), y[j])
     end
-    r += (x₀*transpose(sub) + x₊*dv[nx]) * y[nx]
+    r += dot(adjoint(x₀*transpose(sub) + x₊*dv[nx]), y[nx])
     return r
 end
 
@@ -689,11 +689,11 @@ function dot(x::AbstractVector, A::Tridiagonal, y::AbstractVector)
     x₀ = adjoint(x[1])
     x₊ = adjoint(x[2])
     dl, d, du = A.dl, A.d, A.du
-    r = (x₀*d[1] + x₊*dl[1]) * y[1]
+    r = dot(adjoint(x₀*d[1] + x₊*dl[1]), y[1])
     @inbounds for j in 2:nx-1
         x₋, x₀, x₊ = x₀, x₊, adjoint(x[j+1])
-        r += (x₋*du[j-1] + x₀*d[j] + x₊*dl[j]) * y[j]
+        r += dot(adjoint(x₋*du[j-1] + x₀*d[j] + x₊*dl[j]), y[j])
     end
-    r += (x₀*du[nx-1] + x₊*d[nx]) * y[nx]
+    r += dot(adjoint(x₀*du[nx-1] + x₊*d[nx]), y[nx])
     return r
 end

From 6aea6bb09294b0f767f4b8ae30171c3dff6b73f6 Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@gmx.de>
Date: Tue, 6 Aug 2019 22:42:14 +0200
Subject: [PATCH 13/35] add generalized dot for symmetric/Hermitian matrices

---
 stdlib/LinearAlgebra/src/symmetric.jl  | 10 ++++++++++
 stdlib/LinearAlgebra/test/symmetric.jl |  4 ++++
 2 files changed, 14 insertions(+)

diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl
index 918e5747e3c69..04d7123023e4d 100644
--- a/stdlib/LinearAlgebra/src/symmetric.jl
+++ b/stdlib/LinearAlgebra/src/symmetric.jl
@@ -457,6 +457,16 @@ end
 
 *(A::HermOrSym, B::HermOrSym) = A * copyto!(similar(parent(B)), B)
 
+function dot(x::AbstractVector, A::RealHermSymComplexHerm, y::AbstractVector)
+    Δ = A.uplo == 'U' ? UpperTriangular(A) : LowerTriangular(A)
+    D = Diagonal(A)
+    if x === y
+        return 2dot(x, Δ, x) - dot(x, D, x)
+    else
+        return dot(x, Δ, y) + dot(y, Δ, x) - dot(x, D, x)
+    end
+end
+
 # Fallbacks to avoid generic_matvecmul!/generic_matmatmul!
 ## Symmetric{<:Number} and Hermitian{<:Real} are invariant to transpose; peel off the t
 *(transA::Transpose{<:Any,<:RealHermSymComplexSym}, B::AbstractVector) = transA.parent * B
diff --git a/stdlib/LinearAlgebra/test/symmetric.jl b/stdlib/LinearAlgebra/test/symmetric.jl
index a7538fcd12382..e8a8f078c118d 100644
--- a/stdlib/LinearAlgebra/test/symmetric.jl
+++ b/stdlib/LinearAlgebra/test/symmetric.jl
@@ -365,6 +365,10 @@ end
                 @test Symmetric(asym)\b  ≈ asym\b
             end
         end
+        @testset "generalized dot product" begin
+            @test dot(x, aherm, y) ≈ x'aherm*y ≈ x'*aherm*y
+            @test dot(x, aherm, x) ≈ x'aherm*x ≈ x'*aherm*x
+        end
     end
 end
 

From 4855aa54a81f6e939e7581a2c0fb4b50c3e690d6 Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@gmx.de>
Date: Wed, 7 Aug 2019 08:55:35 +0200
Subject: [PATCH 14/35] fix triangular case

---
 stdlib/LinearAlgebra/src/triangular.jl | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/stdlib/LinearAlgebra/src/triangular.jl b/stdlib/LinearAlgebra/src/triangular.jl
index 1d3b16d7544de..304bfaa474678 100644
--- a/stdlib/LinearAlgebra/src/triangular.jl
+++ b/stdlib/LinearAlgebra/src/triangular.jl
@@ -552,9 +552,9 @@ function dot(x::AbstractVector, A::UpperTriangular, y::AbstractVector)
     if iszero(m)
         return dot(zero(eltype(x)), zero(eltype(A)), zero(eltype(y)))
     end
-    x₁ = first(x)
-    r = zero(typeof(dot(x₁, first(A), first(y))))
-    @inbounds for j in 1:m
+    x₁ = x[1]
+    r = dot(x₁, A[1,1], y[1])
+    @inbounds for j in 2:m
         yj = y[j]
         if !iszero(yj)
             temp = adjoint(x₁) * A[1,j]
@@ -615,8 +615,8 @@ function dot(x::AbstractVector, A::UnitLowerTriangular, y::AbstractVector)
     if iszero(m)
         return dot(zero(eltype(x)), zero(eltype(A)), zero(eltype(y)))
     end
-    r = dot(x[1], y[1])
-    @inbounds for j in 2:m
+    r = zero(typeof(dot(first(x), first(y))))
+    @inbounds for j in 1:m
         yj = y[j]
         if !iszero(yj)
             temp = dot(x[j], yj)

From 427b849c0ef6ff17f6629e8fd503e9467227055d Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@gmx.de>
Date: Wed, 7 Aug 2019 10:28:35 +0200
Subject: [PATCH 15/35] more complete tests for Symmetric/Hermitian

---
 stdlib/LinearAlgebra/src/symmetric.jl  |  2 +-
 stdlib/LinearAlgebra/test/symmetric.jl | 12 ++++++++++--
 2 files changed, 11 insertions(+), 3 deletions(-)

diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl
index 04d7123023e4d..a96540cf47658 100644
--- a/stdlib/LinearAlgebra/src/symmetric.jl
+++ b/stdlib/LinearAlgebra/src/symmetric.jl
@@ -463,7 +463,7 @@ function dot(x::AbstractVector, A::RealHermSymComplexHerm, y::AbstractVector)
     if x === y
         return 2dot(x, Δ, x) - dot(x, D, x)
     else
-        return dot(x, Δ, y) + dot(y, Δ, x) - dot(x, D, x)
+        return dot(x, Δ, y) + dot(y, Δ, x) - dot(x, D, y)
     end
 end
 
diff --git a/stdlib/LinearAlgebra/test/symmetric.jl b/stdlib/LinearAlgebra/test/symmetric.jl
index e8a8f078c118d..dfe75563113ea 100644
--- a/stdlib/LinearAlgebra/test/symmetric.jl
+++ b/stdlib/LinearAlgebra/test/symmetric.jl
@@ -366,8 +366,16 @@ end
             end
         end
         @testset "generalized dot product" begin
-            @test dot(x, aherm, y) ≈ x'aherm*y ≈ x'*aherm*y
-            @test dot(x, aherm, x) ≈ x'aherm*x ≈ x'*aherm*x
+            for uplo in (:U, :L)
+                @test dot(x, Hermitian(aherm, uplo), y) ≈ x'aherm*y ≈ x'*Hermitian(aherm, uplo)*y
+                @test dot(x, Hermitian(aherm, uplo), x) ≈ x'aherm*x ≈ x'*Hermitian(aherm, uplo)*x
+            end
+            if eltya <: Real
+                for uplo in (:U, :L)
+                    @test dot(x, Symmetric(aherm, uplo), y) ≈ x'aherm*y ≈ x'*Symmetric(aherm, uplo)*y
+                    @test dot(x, Symmetric(aherm, uplo), x) ≈ x'aherm*x ≈ x'*Symmetric(aherm, uplo)*x
+                end
+            end
         end
     end
 end

From 93b59be3ca88d9d86ebab41be671c524316a3d8d Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@gmx.de>
Date: Wed, 7 Aug 2019 10:29:00 +0200
Subject: [PATCH 16/35] fix UnitLowerTriangular case

---
 stdlib/LinearAlgebra/src/triangular.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/stdlib/LinearAlgebra/src/triangular.jl b/stdlib/LinearAlgebra/src/triangular.jl
index 304bfaa474678..d93f4f7cfe3f4 100644
--- a/stdlib/LinearAlgebra/src/triangular.jl
+++ b/stdlib/LinearAlgebra/src/triangular.jl
@@ -619,7 +619,7 @@ function dot(x::AbstractVector, A::UnitLowerTriangular, y::AbstractVector)
     @inbounds for j in 1:m
         yj = y[j]
         if !iszero(yj)
-            temp = dot(x[j], yj)
+            temp = adjoint(x[j])
             @simd for i in j+1:m
                 temp += adjoint(x[i]) * A[i,j]
             end

From 791be8b23fdb633e97e66ff8d39c22a54cc2e9c4 Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@gmx.de>
Date: Wed, 7 Aug 2019 12:07:35 +0200
Subject: [PATCH 17/35] fix complex case in symmetric gendot

---
 stdlib/LinearAlgebra/src/symmetric.jl  | 15 +++++++++++++--
 stdlib/LinearAlgebra/test/symmetric.jl |  8 ++++----
 2 files changed, 17 insertions(+), 6 deletions(-)

diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl
index a96540cf47658..c2dfa75c9ea83 100644
--- a/stdlib/LinearAlgebra/src/symmetric.jl
+++ b/stdlib/LinearAlgebra/src/symmetric.jl
@@ -458,12 +458,23 @@ end
 *(A::HermOrSym, B::HermOrSym) = A * copyto!(similar(parent(B)), B)
 
 function dot(x::AbstractVector, A::RealHermSymComplexHerm, y::AbstractVector)
+    require_one_based_indexing(x, y)
+    (length(x) == length(y) == size(A, 1)) || throw(DimensionMismatch())
     Δ = A.uplo == 'U' ? UpperTriangular(A) : LowerTriangular(A)
     D = Diagonal(A)
     if x === y
-        return 2dot(x, Δ, x) - dot(x, D, x)
+        r = 2real(dot(x, Δ, x))
+        @inbounds @simd for i in eachindex(x)
+            r -= dot(x[i], A[i,i], x[i])
+        end
+        return r
     else
-        return dot(x, Δ, y) + dot(y, Δ, x) - dot(x, D, y)
+        r = dot(x, Δ, y)
+        r += conj(dot(y, Δ, x))
+        @inbounds @simd for i in eachindex(x)
+            r -= dot(x[i], A[i,i], y[i])
+        end
+        return r
     end
 end
 
diff --git a/stdlib/LinearAlgebra/test/symmetric.jl b/stdlib/LinearAlgebra/test/symmetric.jl
index dfe75563113ea..f1c93c5f1e857 100644
--- a/stdlib/LinearAlgebra/test/symmetric.jl
+++ b/stdlib/LinearAlgebra/test/symmetric.jl
@@ -367,13 +367,13 @@ end
         end
         @testset "generalized dot product" begin
             for uplo in (:U, :L)
-                @test dot(x, Hermitian(aherm, uplo), y) ≈ x'aherm*y ≈ x'*Hermitian(aherm, uplo)*y
-                @test dot(x, Hermitian(aherm, uplo), x) ≈ x'aherm*x ≈ x'*Hermitian(aherm, uplo)*x
+                @test dot(x, Hermitian(aherm, uplo), y) ≈ dot(x, Hermitian(aherm, uplo)*y) ≈ dot(x, Matrix(Hermitian(aherm, uplo))*y)
+                @test dot(x, Hermitian(aherm, uplo), x) ≈ dot(x, Hermitian(aherm, uplo)*x) ≈ dot(x, Matrix(Hermitian(aherm, uplo))*x)
             end
             if eltya <: Real
                 for uplo in (:U, :L)
-                    @test dot(x, Symmetric(aherm, uplo), y) ≈ x'aherm*y ≈ x'*Symmetric(aherm, uplo)*y
-                    @test dot(x, Symmetric(aherm, uplo), x) ≈ x'aherm*x ≈ x'*Symmetric(aherm, uplo)*x
+                    @test dot(x, Symmetric(aherm, uplo), y) ≈ dot(x, Symmetric(aherm, uplo)*y) ≈ dot(x, Matrix(Symmetric(aherm, uplo))*y)
+                    @test dot(x, Symmetric(aherm, uplo), x) ≈ dot(x, Symmetric(aherm, uplo)*x) ≈ dot(x, Matrix(Symmetric(aherm, uplo))*x)
                 end
             end
         end

From 5d6bbbbba9dc6a8c6efaf9a4e2fc10f44e7611d8 Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@gmx.de>
Date: Wed, 7 Aug 2019 14:02:40 +0200
Subject: [PATCH 18/35] interpret dot(x, A, y) as dot(A'x, y), test accordingly

---
 stdlib/LinearAlgebra/src/bidiag.jl      | 20 ++++++++++----------
 stdlib/LinearAlgebra/src/generic.jl     | 12 +++++++-----
 stdlib/LinearAlgebra/src/matmul.jl      |  3 ---
 stdlib/LinearAlgebra/src/triangular.jl  | 24 ++++++++++++------------
 stdlib/LinearAlgebra/src/tridiag.jl     | 24 ++++++++++++------------
 stdlib/LinearAlgebra/test/bidiag.jl     |  2 +-
 stdlib/LinearAlgebra/test/generic.jl    | 12 ++++++------
 stdlib/LinearAlgebra/test/symmetric.jl  |  8 ++++----
 stdlib/LinearAlgebra/test/triangular.jl |  4 ++--
 stdlib/LinearAlgebra/test/tridiag.jl    |  2 +-
 10 files changed, 55 insertions(+), 56 deletions(-)

diff --git a/stdlib/LinearAlgebra/src/bidiag.jl b/stdlib/LinearAlgebra/src/bidiag.jl
index 758e5cbd1dd02..8e707c4516308 100644
--- a/stdlib/LinearAlgebra/src/bidiag.jl
+++ b/stdlib/LinearAlgebra/src/bidiag.jl
@@ -656,23 +656,23 @@ function dot(x::AbstractVector, B::Bidiagonal, y::AbstractVector)
     end
     ev, dv = B.ev, B.dv
     if B.uplo == 'U'
-        x₀ = adjoint(x[1])
+        x₀ = x[1]
         r = dot(x[1], dv[1], y[1])
         @inbounds for j in 2:nx-1
-            x₋, x₀ = x₀, adjoint(x[j])
-            r += dot(adjoint(x₋*ev[j-1] + x₀*dv[j]), y[j])
+            x₋, x₀ = x₀, x[j]
+            r += dot(adjoint(ev[j-1])*x₋ + adjoint(dv[j])*x₀, y[j])
         end
-        r += dot(adjoint(x₀*ev[nx-1] + adjoint(x[nx])*dv[nx]), y[nx])
+        r += dot(adjoint(ev[nx-1])*x₀ + adjoint(dv[nx])*x[nx], y[nx])
         return r
     else # B.uplo == 'L'
-        x₀ = adjoint(x[1])
-        x₊ = adjoint(x[2])
-        r = dot(adjoint(x₀*dv[1] + x₊*ev[1]), y[1])
+        x₀ = x[1]
+        x₊ = x[2]
+        r = dot(adjoint(dv[1])*x₀ + adjoint(ev[1])*x₊, y[1])
         @inbounds for j in 2:nx-1
-            x₀, x₊ = x₊, adjoint(x[j+1])
-            r += dot(adjoint(x₀*dv[j] + x₊*ev[j]), y[j])
+            x₀, x₊ = x₊, x[j+1]
+            r += dot(adjoint(dv[j])*x₀ + adjoint(ev[j])*x₊, y[j])
         end
-        r += dot(adjoint(x₊), dv[nx], y[nx])
+        r += dot(x₊, dv[nx], y[nx])
         return r
     end
 end
diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl
index 880a75c5745ed..a5c4f4cba5182 100644
--- a/stdlib/LinearAlgebra/src/generic.jl
+++ b/stdlib/LinearAlgebra/src/generic.jl
@@ -877,8 +877,10 @@ end
 """
     dot(x, A, y)
 
-Compute the generalized dot product `x' * A * y` between two vectors `x` and `y`.
-For complex vectors, the first vector is conjugated.
+Compute the generalized dot product `dot(A'x, y)` between two vectors `x` and `y`,
+without storing the intermediate result of `A'x`. As for the two-argument
+[`dot(_,_)`](@ref), this acts recursively. Moreover, for complex vectors, the
+first vector is conjugated.
 """
 dot(x::Number, A::Number, y::Number) = conj(x) * A * y
 
@@ -891,11 +893,11 @@ function dot(x::AbstractVector, A::AbstractMatrix, y::AbstractVector)
     @inbounds for j in eachindex(y)
         yj = y[j]
         if !iszero(yj)
-            temp = zero(adjoint(x₁) * A[i₁,j])
+            temp = zero(adjoint(A[i₁,j]) * x₁)
             @simd for i in eachindex(x)
-                temp += adjoint(x[i]) * A[i,j]
+                temp += adjoint(A[i,j]) * x[i]
             end
-            s += dot(adjoint(temp), yj)
+            s += dot(temp, yj)
         end
     end
     return s
diff --git a/stdlib/LinearAlgebra/src/matmul.jl b/stdlib/LinearAlgebra/src/matmul.jl
index 1bafd979ea2f1..0d6ac7a62e1ed 100644
--- a/stdlib/LinearAlgebra/src/matmul.jl
+++ b/stdlib/LinearAlgebra/src/matmul.jl
@@ -40,9 +40,6 @@ function *(transx::Transpose{<:Any,<:StridedVector{T}}, y::StridedVector{T}) whe
     return BLAS.dotu(x, y)
 end
 
-# I'm no longer sure we want ternary multiplication to be generally interpreted recursively
-# (*)(x::Adjoint{<:Any,<:AbstractVector}, A::AbstractMatrix, y::AbstractVector) = dot(parent(x), A, y)
-
 # Matrix-vector multiplication
 function (*)(A::StridedMatrix{T}, x::StridedVector{S}) where {T<:BlasFloat,S<:Real}
     TS = promote_op(matprod, T, S)
diff --git a/stdlib/LinearAlgebra/src/triangular.jl b/stdlib/LinearAlgebra/src/triangular.jl
index d93f4f7cfe3f4..1c156fa6bf326 100644
--- a/stdlib/LinearAlgebra/src/triangular.jl
+++ b/stdlib/LinearAlgebra/src/triangular.jl
@@ -557,11 +557,11 @@ function dot(x::AbstractVector, A::UpperTriangular, y::AbstractVector)
     @inbounds for j in 2:m
         yj = y[j]
         if !iszero(yj)
-            temp = adjoint(x₁) * A[1,j]
+            temp = adjoint(A[1,j]) * x₁
             @simd for i in 2:j
-                temp += adjoint(x[i]) * A[i,j]
+                temp += adjoint(A[i,j]) * x[i]
             end
-            r += dot(adjoint(temp), yj)
+            r += dot(temp, yj)
         end
     end
     return r
@@ -578,11 +578,11 @@ function dot(x::AbstractVector, A::UnitUpperTriangular, y::AbstractVector)
     @inbounds for j in 2:m
         yj = y[j]
         if !iszero(yj)
-            temp = adjoint(x₁) * A[1,j]
+            temp = adjoint(A[1,j]) * x₁
             @simd for i in 2:j-1
-                temp += adjoint(x[i]) * A[i,j]
+                temp += adjoint(A[i,j]) * x[i]
             end
-            r += dot(adjoint(temp), yj)
+            r += dot(temp, yj)
             r += dot(x[j], yj)
         end
     end
@@ -599,11 +599,11 @@ function dot(x::AbstractVector, A::LowerTriangular, y::AbstractVector)
     @inbounds for j in 1:m
         yj = y[j]
         if !iszero(yj)
-            temp = adjoint(x[j]) * A[j,j]
+            temp = adjoint(A[j,j]) * x[j]
             @simd for i in j+1:m
-                temp += adjoint(x[i]) * A[i,j]
+                temp += adjoint(A[i,j]) * x[i]
             end
-            r += dot(adjoint(temp), yj)
+            r += dot(temp, yj)
         end
     end
     return r
@@ -619,11 +619,11 @@ function dot(x::AbstractVector, A::UnitLowerTriangular, y::AbstractVector)
     @inbounds for j in 1:m
         yj = y[j]
         if !iszero(yj)
-            temp = adjoint(x[j])
+            temp = x[j]
             @simd for i in j+1:m
-                temp += adjoint(x[i]) * A[i,j]
+                temp += adjoint(A[i,j]) * x[i]
             end
-            r += dot(adjoint(temp), yj)
+            r += dot(temp, yj)
         end
     end
     return r
diff --git a/stdlib/LinearAlgebra/src/tridiag.jl b/stdlib/LinearAlgebra/src/tridiag.jl
index 413a2b0fe9a92..670ef8dcd1fe0 100644
--- a/stdlib/LinearAlgebra/src/tridiag.jl
+++ b/stdlib/LinearAlgebra/src/tridiag.jl
@@ -210,16 +210,16 @@ function dot(x::AbstractVector, S::SymTridiagonal, y::AbstractVector)
         return dot(zero(eltype(x)), zero(eltype(S)), zero(eltype(y)))
     end
     dv, ev = S.dv, S.ev
-    x₀ = adjoint(x[1])
-    x₊ = adjoint(x[2])
+    x₀ = x[1]
+    x₊ = x[2]
     sub = transpose(ev[1])
-    r = dot(adjoint(x₀*dv[1] + x₊*sub), y[1])
+    r = dot(adjoint(dv[1])*x₀ + adjoint(sub)*x₊, y[1])
     @inbounds for j in 2:nx-1
-        x₋, x₀, x₊ = x₀, x₊, adjoint(x[j+1])
+        x₋, x₀, x₊ = x₀, x₊, x[j+1]
         sup, sub = transpose(sub), transpose(ev[j])
-        r += dot(adjoint(x₋*sup + x₀*dv[j] + x₊*sub), y[j])
+        r += dot(adjoint(sup)*x₋ + adjoint(dv[j])*x₀ + adjoint(sub)*x₊, y[j])
     end
-    r += dot(adjoint(x₀*transpose(sub) + x₊*dv[nx]), y[nx])
+    r += dot(adjoint(transpose(sub))*x₀ + adjoint(dv[nx])*x₊, y[nx])
     return r
 end
 
@@ -686,14 +686,14 @@ function dot(x::AbstractVector, A::Tridiagonal, y::AbstractVector)
     if iszero(nx)
         return dot(zero(eltype(x)), zero(eltype(A)), zero(eltype(y)))
     end
-    x₀ = adjoint(x[1])
-    x₊ = adjoint(x[2])
+    x₀ = x[1]
+    x₊ = x[2]
     dl, d, du = A.dl, A.d, A.du
-    r = dot(adjoint(x₀*d[1] + x₊*dl[1]), y[1])
+    r = dot(adjoint(d[1])*x₀ + adjoint(dl[1])*x₊, y[1])
     @inbounds for j in 2:nx-1
-        x₋, x₀, x₊ = x₀, x₊, adjoint(x[j+1])
-        r += dot(adjoint(x₋*du[j-1] + x₀*d[j] + x₊*dl[j]), y[j])
+        x₋, x₀, x₊ = x₀, x₊, x[j+1]
+        r += dot(adjoint(du[j-1])*x₋ + adjoint(d[j])*x₀ + adjoint(dl[j])*x₊, y[j])
     end
-    r += dot(adjoint(x₀*du[nx-1] + x₊*d[nx]), y[nx])
+    r += dot(adjoint(du[nx-1])*x₀ + adjoint(d[nx])*x₊, y[nx])
     return r
 end
diff --git a/stdlib/LinearAlgebra/test/bidiag.jl b/stdlib/LinearAlgebra/test/bidiag.jl
index b0578687d5271..d55fbc32a7242 100644
--- a/stdlib/LinearAlgebra/test/bidiag.jl
+++ b/stdlib/LinearAlgebra/test/bidiag.jl
@@ -463,7 +463,7 @@ end
         y = randn(elty, 5)
         for uplo in (:U, :L)
             B = Bidiagonal(dv, ev, uplo)
-            @test dot(x, B, y) ≈ x' * B * y ≈ *(x', B, y) ≈ (x'B)*y
+            @test dot(x, B, y) ≈ dot(B'x, y) ≈ dot(x, Matrix(B), y)
         end
     end
 end
diff --git a/stdlib/LinearAlgebra/test/generic.jl b/stdlib/LinearAlgebra/test/generic.jl
index 37dc3906925aa..113fa15db3366 100644
--- a/stdlib/LinearAlgebra/test/generic.jl
+++ b/stdlib/LinearAlgebra/test/generic.jl
@@ -425,15 +425,15 @@ end
             x = rand(elty, n)
             y = rand(elty, n)
         end
-        @test dot(x, A, y) ≈ x' * A * y ≈ *(x', A, y) ≈ (x'A)*y
-        @test dot(x, A', y) ≈ x' * A' * y ≈ *(x', A', y) ≈ (x'A')*y
-        elty <: Real && @test dot(x, transpose(A), y) ≈ x' * transpose(A) * y ≈ *(x', transpose(A), y) ≈ (x'*transpose(A))*y
+        @test dot(x, A, y) ≈ dot(A'x, y) ≈ *(x', A, y) ≈ (x'A)*y
+        @test dot(x, A', y) ≈ dot(A*x, y) ≈ *(x', A', y) ≈ (x'A')*y
+        elty <: Real && @test dot(x, transpose(A), y) ≈ dot(x, transpose(A)*y) ≈ *(x', transpose(A), y) ≈ (x'*transpose(A))*y
         B = reshape([A], 1, 1)
         x = [x]
         y = [y]
-        @test dot(x, B, y) ≈ x' * B * y ≈ *(x', B, y) ≈ (x'B)*y
-        @test dot(x, B', y) ≈ x' * B' * y ≈ *(x', B', y) ≈ (x'B')*y
-        elty <: Real && @test dot(x, transpose(B), y) ≈ x' * transpose(B) * y ≈ *(x', transpose(B), y) ≈ (x'*transpose(B))*y
+        @test dot(x, B, y) ≈ dot(B'x, y)
+        @test dot(x, B', y) ≈ dot(B*x, y)
+        elty <: Real && @test dot(x, transpose(B), y) ≈ dot(x, transpose(B)*y)
     end
 end
 
diff --git a/stdlib/LinearAlgebra/test/symmetric.jl b/stdlib/LinearAlgebra/test/symmetric.jl
index f1c93c5f1e857..ebdfc9c93c72f 100644
--- a/stdlib/LinearAlgebra/test/symmetric.jl
+++ b/stdlib/LinearAlgebra/test/symmetric.jl
@@ -367,13 +367,13 @@ end
         end
         @testset "generalized dot product" begin
             for uplo in (:U, :L)
-                @test dot(x, Hermitian(aherm, uplo), y) ≈ dot(x, Hermitian(aherm, uplo)*y) ≈ dot(x, Matrix(Hermitian(aherm, uplo))*y)
-                @test dot(x, Hermitian(aherm, uplo), x) ≈ dot(x, Hermitian(aherm, uplo)*x) ≈ dot(x, Matrix(Hermitian(aherm, uplo))*x)
+                @test dot(x, Hermitian(aherm, uplo), y) ≈ dot(x, Hermitian(aherm, uplo)*y) ≈ dot(x, Matrix(Hermitian(aherm, uplo)), y)
+                @test dot(x, Hermitian(aherm, uplo), x) ≈ dot(x, Hermitian(aherm, uplo)*x) ≈ dot(x, Matrix(Hermitian(aherm, uplo)), x)
             end
             if eltya <: Real
                 for uplo in (:U, :L)
-                    @test dot(x, Symmetric(aherm, uplo), y) ≈ dot(x, Symmetric(aherm, uplo)*y) ≈ dot(x, Matrix(Symmetric(aherm, uplo))*y)
-                    @test dot(x, Symmetric(aherm, uplo), x) ≈ dot(x, Symmetric(aherm, uplo)*x) ≈ dot(x, Matrix(Symmetric(aherm, uplo))*x)
+                    @test dot(x, Symmetric(aherm, uplo), y) ≈ dot(x, Symmetric(aherm, uplo)*y) ≈ dot(x, Matrix(Symmetric(aherm, uplo)), y)
+                    @test dot(x, Symmetric(aherm, uplo), x) ≈ dot(x, Symmetric(aherm, uplo)*x) ≈ dot(x, Matrix(Symmetric(aherm, uplo)), x)
                 end
             end
         end
diff --git a/stdlib/LinearAlgebra/test/triangular.jl b/stdlib/LinearAlgebra/test/triangular.jl
index d207d3ea1cdbf..fc48b840b8d16 100644
--- a/stdlib/LinearAlgebra/test/triangular.jl
+++ b/stdlib/LinearAlgebra/test/triangular.jl
@@ -241,9 +241,9 @@ for elty1 in (Float32, Float64, BigFloat, ComplexF32, ComplexF64, Complex{BigFlo
             b1 = convert(Vector{eltyb}, (elty1 <: Complex ? real(A1) : A1)*fill(1., n))
             b2 = convert(Vector{eltyb}, (elty1 <: Complex ? real(A1) : A1)*randn(n))
             if elty1 in (BigFloat, Complex{BigFloat}) || eltyb in (BigFloat, Complex{BigFloat})
-                @test dot(b1, A1, b2) ≈ (b1'A1)*b2  atol=sqrt(max(eps(real(float(one(elty1)))),eps(real(float(one(eltyb))))))*n*n
+                @test dot(b1, A1, b2) ≈ dot(A1'b1, b2)  atol=sqrt(max(eps(real(float(one(elty1)))),eps(real(float(one(eltyb))))))*n*n
             else
-                @test dot(b1, A1, b2) ≈ (b1'A1)*b2
+                @test dot(b1, A1, b2) ≈ dot(A1'b1, b2)
             end
         end
 
diff --git a/stdlib/LinearAlgebra/test/tridiag.jl b/stdlib/LinearAlgebra/test/tridiag.jl
index bc3eb4c63e99e..8431c81fd3ef8 100644
--- a/stdlib/LinearAlgebra/test/tridiag.jl
+++ b/stdlib/LinearAlgebra/test/tridiag.jl
@@ -393,7 +393,7 @@ end
         @testset "generalized dot" begin
             x = fill(convert(elty, 1), n)
             y = fill(convert(elty, 1), n)
-            @test dot(x, A, y) ≈ x' * A * y ≈ *(x', A, y) ≈ (x'A)*y
+            @test dot(x, A, y) ≈ dot(A'x, y)
         end
     end
 end

From c176bfce3b433eed136922cda4fd9fa9dcb7d6bc Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@gmx.de>
Date: Wed, 7 Aug 2019 16:05:38 +0200
Subject: [PATCH 19/35] use correct tolerance in triangular tests

---
 stdlib/LinearAlgebra/test/triangular.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/stdlib/LinearAlgebra/test/triangular.jl b/stdlib/LinearAlgebra/test/triangular.jl
index fc48b840b8d16..5f6ec8771a030 100644
--- a/stdlib/LinearAlgebra/test/triangular.jl
+++ b/stdlib/LinearAlgebra/test/triangular.jl
@@ -243,7 +243,7 @@ for elty1 in (Float32, Float64, BigFloat, ComplexF32, ComplexF64, Complex{BigFlo
             if elty1 in (BigFloat, Complex{BigFloat}) || eltyb in (BigFloat, Complex{BigFloat})
                 @test dot(b1, A1, b2) ≈ dot(A1'b1, b2)  atol=sqrt(max(eps(real(float(one(elty1)))),eps(real(float(one(eltyb))))))*n*n
             else
-                @test dot(b1, A1, b2) ≈ dot(A1'b1, b2)
+                @test dot(b1, A1, b2) ≈ dot(A1'b1, b2)  atol=sqrt(max(eps(real(float(one(elty1)))),eps(real(float(one(eltyb))))))*n*n
             end
         end
 

From 9510dfd6b0a7c8821d0d073d59265fc6a07431fd Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@gmx.de>
Date: Thu, 8 Aug 2019 10:51:51 +0200
Subject: [PATCH 20/35] add gendot for UpperHessenberg, and tests

---
 stdlib/LinearAlgebra/src/hessenberg.jl  | 31 +++++++++++++++++++++++++
 stdlib/LinearAlgebra/test/hessenberg.jl |  5 ++++
 2 files changed, 36 insertions(+)

diff --git a/stdlib/LinearAlgebra/src/hessenberg.jl b/stdlib/LinearAlgebra/src/hessenberg.jl
index 341d2c71bf332..5f9b2be1d34a2 100644
--- a/stdlib/LinearAlgebra/src/hessenberg.jl
+++ b/stdlib/LinearAlgebra/src/hessenberg.jl
@@ -284,6 +284,37 @@ function logabsdet(F::UpperHessenberg; shift::Number=false)
     return (logdeterminant, P)
 end
 
+function dot(x::AbstractVector, H::UpperHessenberg, y::AbstractVector)
+    require_one_based_indexing(x, y)
+    m = size(H, 1)
+    (length(x) == m == length(y)) || throw(DimensionMismatch())
+    if iszero(m)
+        return dot(zero(eltype(x)), zero(eltype(H)), zero(eltype(y)))
+    end
+    x₁ = x[1]
+    r = dot(x₁, H[1,1], y[1])
+    r += dot(x[2], H[2,1], y[1])
+    @inbounds for j in 2:m-1
+        yj = y[j]
+        if !iszero(yj)
+            temp = adjoint(H[1,j]) * x₁
+            @simd for i in 2:j+1
+                temp += adjoint(H[i,j]) * x[i]
+            end
+            r += dot(temp, yj)
+        end
+    end
+    ym = y[m]
+    if !iszero(ym)
+        temp = adjoint(H[1,m]) * x₁
+        @simd for i in 2:m
+            temp += adjoint(H[i,m]) * x[i]
+        end
+        r += dot(temp, ym)
+    end
+    return r
+end
+
 ######################################################################################
 # Hessenberg factorizations Q(H+μI)Q' of A+μI:
 
diff --git a/stdlib/LinearAlgebra/test/hessenberg.jl b/stdlib/LinearAlgebra/test/hessenberg.jl
index 4db168dfe9189..b7cc07d09446f 100644
--- a/stdlib/LinearAlgebra/test/hessenberg.jl
+++ b/stdlib/LinearAlgebra/test/hessenberg.jl
@@ -88,6 +88,11 @@ let n = 10
             @test det(H + shift*I) ≈ det(A + shift*I)
             @test logabsdet(H + shift*I) ≅ logabsdet(A + shift*I)
         end
+
+        HM = Matrix(h)
+        @test dot(b, h, b) ≈ dot(h'b, b) ≈ dot(b, HM, b) ≈ dot(HM'b, b)
+        c = b .+ 1
+        @test dot(b, h, c) ≈ dot(h'b, c) ≈ dot(b, HM, c) ≈ dot(HM'b, c)
     end
 end
 

From a6bbb45e5535b86583d5765c3646f72dbbb38ab5 Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@gmx.de>
Date: Fri, 16 Aug 2019 20:19:05 +0200
Subject: [PATCH 21/35] fix docstring of 3-arg dot

---
 stdlib/LinearAlgebra/src/generic.jl | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl
index a5c4f4cba5182..bf4d725b02371 100644
--- a/stdlib/LinearAlgebra/src/generic.jl
+++ b/stdlib/LinearAlgebra/src/generic.jl
@@ -877,11 +877,13 @@ end
 """
     dot(x, A, y)
 
-Compute the generalized dot product `dot(A'x, y)` between two vectors `x` and `y`,
-without storing the intermediate result of `A'x`. As for the two-argument
+Compute the generalized dot product `dot(x, A*y)` between two vectors `x` and `y`,
+without storing the intermediate result of `A*y`. As for the two-argument
 [`dot(_,_)`](@ref), this acts recursively. Moreover, for complex vectors, the
 first vector is conjugated.
 """
+function dot end
+
 dot(x::Number, A::Number, y::Number) = conj(x) * A * y
 
 function dot(x::AbstractVector, A::AbstractMatrix, y::AbstractVector)

From e4668ec0939025c3e83d74eac1de3b801589533b Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@gmx.de>
Date: Fri, 16 Aug 2019 20:27:25 +0200
Subject: [PATCH 22/35] add generic 3-arg dot for UniformScaling

---
 stdlib/LinearAlgebra/src/uniformscaling.jl | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/stdlib/LinearAlgebra/src/uniformscaling.jl b/stdlib/LinearAlgebra/src/uniformscaling.jl
index f5836bbdf7c07..11f6b860163f2 100644
--- a/stdlib/LinearAlgebra/src/uniformscaling.jl
+++ b/stdlib/LinearAlgebra/src/uniformscaling.jl
@@ -401,4 +401,7 @@ Array(s::UniformScaling, dims::Dims{2}) = Matrix(s, dims)
 Diagonal{T}(s::UniformScaling, m::Integer) where {T} = Diagonal{T}(fill(T(s.λ), m))
 Diagonal(s::UniformScaling, m::Integer) = Diagonal{eltype(s)}(s, m)
 
-dot(x::AbstractVector, J::UniformScaling, y::AbstractVector) = J.λ * dot(x, y)
+function dot(x::AbstractVector, J::UniformScaling, y::AbstractVector)
+    return mapreduce(t -> dot(t[1], J.λ, t[2]), +, zip(x, y))
+end
+dot(x::AbstractVector, J::UniformScaling{<:Union{Real,Complex}}, y::AbstractVector) = J.λ*dot(x, y)

From 09322387eb47dd2aced8621c6d01efcd41af470b Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@gmx.de>
Date: Fri, 16 Aug 2019 20:48:45 +0200
Subject: [PATCH 23/35] add generic fallback

This should be only relevant to cases like `dot(x, J, y)`, where `x` and `y` are vectors of quaternion vectors, and `J` is a quaternion `UniformScaling`.
---
 stdlib/LinearAlgebra/src/generic.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl
index bf4d725b02371..d6bb7eef2ab7a 100644
--- a/stdlib/LinearAlgebra/src/generic.jl
+++ b/stdlib/LinearAlgebra/src/generic.jl
@@ -884,6 +884,7 @@ first vector is conjugated.
 """
 function dot end
 
+dot(x, A, y) = dot(x, A*y) # generic fallback for cases that are not covered by specialized methods
 dot(x::Number, A::Number, y::Number) = conj(x) * A * y
 
 function dot(x::AbstractVector, A::AbstractMatrix, y::AbstractVector)

From 6310fa65a2e96a1d78848cd464eaf13432e17dda Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@gmx.de>
Date: Fri, 16 Aug 2019 21:01:21 +0200
Subject: [PATCH 24/35] add gendot with middle argument Number

---
 stdlib/LinearAlgebra/src/uniformscaling.jl | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/stdlib/LinearAlgebra/src/uniformscaling.jl b/stdlib/LinearAlgebra/src/uniformscaling.jl
index 11f6b860163f2..f4e9d90432ed4 100644
--- a/stdlib/LinearAlgebra/src/uniformscaling.jl
+++ b/stdlib/LinearAlgebra/src/uniformscaling.jl
@@ -404,4 +404,8 @@ Diagonal(s::UniformScaling, m::Integer) = Diagonal{eltype(s)}(s, m)
 function dot(x::AbstractVector, J::UniformScaling, y::AbstractVector)
     return mapreduce(t -> dot(t[1], J.λ, t[2]), +, zip(x, y))
 end
+function dot(x::AbstractVector, a::Number, y::AbstractVector)
+    return mapreduce(t -> dot(t[1], a, t[2]), +, zip(x, y))
+end
 dot(x::AbstractVector, J::UniformScaling{<:Union{Real,Complex}}, y::AbstractVector) = J.λ*dot(x, y)
+dot(x::AbstractVector, a::Union{Real,Complex}, y::AbstractVector) = a*dot(x, y)

From a66a6a1181bb3da6c17ca3fa6ed01424e483f942 Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@gmx.de>
Date: Wed, 28 Aug 2019 10:08:24 +0200
Subject: [PATCH 25/35] merge NEWS

---
 NEWS.md | 1 +
 1 file changed, 1 insertion(+)

diff --git a/NEWS.md b/NEWS.md
index 2ff018ce96ce4..012b3b788650e 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -39,6 +39,7 @@ Standard library changes
 
 * `qr` and `qr!` functions support `blocksize` keyword argument ([#33053]).
 
+* `dot` now admits a 3-argument method `dot(x, A, y)` to compute generalized dot products `dot(x, A*y)`, but without computing and storing the intermediate result `A*y` ([#32739]).
 
 #### SparseArrays
 

From fb97cc3bfdc50d67aa3d3bf3041b1fd516e7495d Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@gmx.de>
Date: Sun, 18 Aug 2019 19:09:12 +0200
Subject: [PATCH 26/35] attach docstring to generic fallback

---
 stdlib/LinearAlgebra/src/generic.jl | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl
index d6bb7eef2ab7a..47d0d8c9a3732 100644
--- a/stdlib/LinearAlgebra/src/generic.jl
+++ b/stdlib/LinearAlgebra/src/generic.jl
@@ -882,9 +882,8 @@ without storing the intermediate result of `A*y`. As for the two-argument
 [`dot(_,_)`](@ref), this acts recursively. Moreover, for complex vectors, the
 first vector is conjugated.
 """
-function dot end
-
 dot(x, A, y) = dot(x, A*y) # generic fallback for cases that are not covered by specialized methods
+
 dot(x::Number, A::Number, y::Number) = conj(x) * A * y
 
 function dot(x::AbstractVector, A::AbstractMatrix, y::AbstractVector)

From c9112fc159944d399cd8ea5abc281bc5606c840d Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@gmx.de>
Date: Sun, 18 Aug 2019 19:12:44 +0200
Subject: [PATCH 27/35] simplify scalar/uniform scaling gendot

---
 stdlib/LinearAlgebra/src/uniformscaling.jl | 9 ++-------
 1 file changed, 2 insertions(+), 7 deletions(-)

diff --git a/stdlib/LinearAlgebra/src/uniformscaling.jl b/stdlib/LinearAlgebra/src/uniformscaling.jl
index f4e9d90432ed4..d268120262e94 100644
--- a/stdlib/LinearAlgebra/src/uniformscaling.jl
+++ b/stdlib/LinearAlgebra/src/uniformscaling.jl
@@ -401,11 +401,6 @@ Array(s::UniformScaling, dims::Dims{2}) = Matrix(s, dims)
 Diagonal{T}(s::UniformScaling, m::Integer) where {T} = Diagonal{T}(fill(T(s.λ), m))
 Diagonal(s::UniformScaling, m::Integer) = Diagonal{eltype(s)}(s, m)
 
-function dot(x::AbstractVector, J::UniformScaling, y::AbstractVector)
-    return mapreduce(t -> dot(t[1], J.λ, t[2]), +, zip(x, y))
-end
-function dot(x::AbstractVector, a::Number, y::AbstractVector)
-    return mapreduce(t -> dot(t[1], a, t[2]), +, zip(x, y))
-end
-dot(x::AbstractVector, J::UniformScaling{<:Union{Real,Complex}}, y::AbstractVector) = J.λ*dot(x, y)
+dot(x::AbstractVector, J::UniformScaling, y::AbstractVector) = dot(x, J.λ, y)
+dot(x::AbstractVector, a::Number, y::AbstractVector) = sum(t -> dot(t[1], a, t[2]), zip(x, y))
 dot(x::AbstractVector, a::Union{Real,Complex}, y::AbstractVector) = a*dot(x, y)

From cffa4aa40afa454ca575386f9ace24ac14578f84 Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@gmx.de>
Date: Thu, 29 Aug 2019 13:38:09 +0200
Subject: [PATCH 28/35] use dot(A'x,y) for fallback

---
 stdlib/LinearAlgebra/src/generic.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl
index 47d0d8c9a3732..257708ac00774 100644
--- a/stdlib/LinearAlgebra/src/generic.jl
+++ b/stdlib/LinearAlgebra/src/generic.jl
@@ -882,7 +882,7 @@ without storing the intermediate result of `A*y`. As for the two-argument
 [`dot(_,_)`](@ref), this acts recursively. Moreover, for complex vectors, the
 first vector is conjugated.
 """
-dot(x, A, y) = dot(x, A*y) # generic fallback for cases that are not covered by specialized methods
+dot(x, A, y) = dot(A'x, y) # generic fallback for cases that are not covered by specialized methods
 
 dot(x::Number, A::Number, y::Number) = conj(x) * A * y
 

From 48bdbc1d6147e311a595187dd0075db24a5b8dd1 Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@gmx.de>
Date: Fri, 30 Aug 2019 10:25:24 +0200
Subject: [PATCH 29/35] use accessor functions in sparse code, generalize to
 Abstract..., tests

---
 stdlib/SparseArrays/src/linalg.jl  | 37 +++++++++++++++++++-----------
 stdlib/SparseArrays/test/sparse.jl |  2 ++
 2 files changed, 25 insertions(+), 14 deletions(-)

diff --git a/stdlib/SparseArrays/src/linalg.jl b/stdlib/SparseArrays/src/linalg.jl
index 9393bfab17358..2576a0752086e 100644
--- a/stdlib/SparseArrays/src/linalg.jl
+++ b/stdlib/SparseArrays/src/linalg.jl
@@ -321,17 +321,18 @@ function dot(A::AbstractSparseMatrixCSC{T1,S1},B::AbstractSparseMatrixCSC{T2,S2}
     return r
 end
 
-function dot(x::AbstractVector, A::SparseMatrixCSC, y::AbstractVector)
+function dot(x::AbstractVector, A::AbstractSparseMatrixCSC, y::AbstractVector)
     require_one_based_indexing(x, y)
-    (length(x) == A.m && A.n == length(y)) || throw(DimensionMismatch())
-    if iszero(A.m) || iszero(A.n)
+    m, n = size(A)
+    (length(x) == m && n == length(y)) || throw(DimensionMismatch())
+    if iszero(m) || iszero(n)
         return dot(zero(eltype(x)), zero(eltype(A)), zero(eltype(y)))
     end
     T = promote_type(eltype(x), eltype(A), eltype(y))
     r = zero(T)
-    rvals = rowvals(A)
-    nzvals = nonzeros(A)
-    @inbounds for col in 1:A.n
+    rvals = getrowval(A)
+    nzvals = getnzval(A)
+    @inbounds for col in 1:n
         ycol = y[col]
         if !iszero(ycol)
             temp = zero(T)
@@ -343,18 +344,26 @@ function dot(x::AbstractVector, A::SparseMatrixCSC, y::AbstractVector)
     end
     return r
 end
-function dot(x::SparseVector, A::SparseMatrixCSC, y::SparseVector)
-    x.n == A.m && A.n == y.n || throw(DimensionMismatch())
-    if iszero(A.m) || iszero(A.n)
+function dot(x::SparseVector, A::AbstractSparseMatrixCSC, y::SparseVector)
+    m, n = size(A)
+    length(x) == m && n == length(y) || throw(DimensionMismatch())
+    if iszero(m) || iszero(n)
         return dot(zero(eltype(x)), zero(eltype(A)), zero(eltype(y)))
     end
     r = zero(promote_type(eltype(x), eltype(A), eltype(y)))
-    for (yi, yv) in zip(y.nzind, y.nzval)
-        A_ptr_lo = A.colptr[yi]
-        A_ptr_hi = A.colptr[yi+1] - 1
+    xnzind = nonzeroinds(x)
+    xnzval = nonzeros(x)
+    ynzind = nonzeroinds(y)
+    ynzval = nonzeros(y)
+    Acolptr = getcolptr(A)
+    Arowval = getrowval(A)
+    Anzval = getnzval(A)
+    for (yi, yv) in zip(ynzind, ynzval)
+        A_ptr_lo = Acolptr[yi]
+        A_ptr_hi = Acolptr[yi+1] - 1
         if A_ptr_lo <= A_ptr_hi
-            r += _spdot(dot, 1, length(x.nzind), x.nzind, x.nzval,
-                                            A_ptr_lo, A_ptr_hi, A.rowval, A.nzval) * yv
+            r += _spdot(dot, 1, length(xnzind), xnzind, xnzval,
+                                            A_ptr_lo, A_ptr_hi, Arowval, Anzval) * yv
         end
     end
     r
diff --git a/stdlib/SparseArrays/test/sparse.jl b/stdlib/SparseArrays/test/sparse.jl
index 45099f399fb51..6f21e18ae0bc3 100644
--- a/stdlib/SparseArrays/test/sparse.jl
+++ b/stdlib/SparseArrays/test/sparse.jl
@@ -427,9 +427,11 @@ end
 @testset "generalized dot product" begin
     for i = 1:5
         A = sprand(ComplexF64, 10, 15, 0.4)
+        Av = view(A, :, :)
         x = sprand(ComplexF64, 10, 0.5)
         y = sprand(ComplexF64, 15, 0.5)
         @test dot(x, A, y) ≈ dot(Vector(x), A, Vector(y)) ≈ (Vector(x)' * Matrix(A)) * Vector(y)
+        @test dot(x, A, y) ≈ dot(x, Av, y)
     end
 end
 

From ff377e677600e1b7c92b0bca83461ea488c7efa7 Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@gmx.de>
Date: Fri, 30 Aug 2019 10:26:32 +0200
Subject: [PATCH 30/35] revert fallback definition

---
 stdlib/LinearAlgebra/src/generic.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl
index 257708ac00774..47d0d8c9a3732 100644
--- a/stdlib/LinearAlgebra/src/generic.jl
+++ b/stdlib/LinearAlgebra/src/generic.jl
@@ -882,7 +882,7 @@ without storing the intermediate result of `A*y`. As for the two-argument
 [`dot(_,_)`](@ref), this acts recursively. Moreover, for complex vectors, the
 first vector is conjugated.
 """
-dot(x, A, y) = dot(A'x, y) # generic fallback for cases that are not covered by specialized methods
+dot(x, A, y) = dot(x, A*y) # generic fallback for cases that are not covered by specialized methods
 
 dot(x::Number, A::Number, y::Number) = conj(x) * A * y
 

From 48874e22cbc9565a137d71520aab70ae6abffc8d Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@gmx.de>
Date: Sat, 31 Aug 2019 11:54:49 +0200
Subject: [PATCH 31/35] add compat note and jldoctest

---
 stdlib/LinearAlgebra/src/generic.jl | 15 +++++++++++++++
 1 file changed, 15 insertions(+)

diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl
index 47d0d8c9a3732..4dad22160b166 100644
--- a/stdlib/LinearAlgebra/src/generic.jl
+++ b/stdlib/LinearAlgebra/src/generic.jl
@@ -881,6 +881,21 @@ Compute the generalized dot product `dot(x, A*y)` between two vectors `x` and `y
 without storing the intermediate result of `A*y`. As for the two-argument
 [`dot(_,_)`](@ref), this acts recursively. Moreover, for complex vectors, the
 first vector is conjugated.
+
+!!! compat "Julia 1.4"
+    Three-argument `dot` requires at least Julia 1.4.
+
+# Examples
+```jldoctest
+julia> dot([1; 1], [1 2; 3 4], [2; 3])
+26
+
+julia> dot(1:5, reshape(1:25, 5, 5), 2:6)
+4850
+
+julia> ⋅(1:5, reshape(1:25, 5, 5), 2:6) == dot(1:5, reshape(1:25, 5, 5), 2:6)
+true
+```
 """
 dot(x, A, y) = dot(x, A*y) # generic fallback for cases that are not covered by specialized methods
 

From 55687957d1031212e5d51db9ff4e549d0c43683c Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@gmx.de>
Date: Fri, 30 Aug 2019 16:01:40 +0200
Subject: [PATCH 32/35] remove redundant Number version

---
 stdlib/LinearAlgebra/src/generic.jl | 2 --
 1 file changed, 2 deletions(-)

diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl
index 4dad22160b166..1ffdf95ba6978 100644
--- a/stdlib/LinearAlgebra/src/generic.jl
+++ b/stdlib/LinearAlgebra/src/generic.jl
@@ -899,8 +899,6 @@ true
 """
 dot(x, A, y) = dot(x, A*y) # generic fallback for cases that are not covered by specialized methods
 
-dot(x::Number, A::Number, y::Number) = conj(x) * A * y
-
 function dot(x::AbstractVector, A::AbstractMatrix, y::AbstractVector)
     (axes(x)..., axes(y)...) == axes(A) || throw(DimensionMismatch())
     T = typeof(dot(first(x), first(A), first(y)))

From 8051614bf5a0d54e1b508c2c02f09ac06e0af587 Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@gmx.de>
Date: Fri, 30 Aug 2019 16:34:32 +0200
Subject: [PATCH 33/35] write out loops in symmetric/hermitian case

---
 stdlib/LinearAlgebra/src/symmetric.jl | 30 +++++++++++++++------------
 1 file changed, 17 insertions(+), 13 deletions(-)

diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl
index c2dfa75c9ea83..1f3829d3727bf 100644
--- a/stdlib/LinearAlgebra/src/symmetric.jl
+++ b/stdlib/LinearAlgebra/src/symmetric.jl
@@ -460,22 +460,26 @@ end
 function dot(x::AbstractVector, A::RealHermSymComplexHerm, y::AbstractVector)
     require_one_based_indexing(x, y)
     (length(x) == length(y) == size(A, 1)) || throw(DimensionMismatch())
-    Δ = A.uplo == 'U' ? UpperTriangular(A) : LowerTriangular(A)
-    D = Diagonal(A)
-    if x === y
-        r = 2real(dot(x, Δ, x))
-        @inbounds @simd for i in eachindex(x)
-            r -= dot(x[i], A[i,i], x[i])
+    data = A.data
+    r = zero(eltype(x)) * zero(eltype(A)) * zero(eltype(y))
+    if A.uplo == 'U'
+        @inbounds for j = 1:length(y)
+            r += dot(x[j], real(data[j,j]), y[j])
+            @simd for i = 1:j-1
+                Aij = data[i,j]
+                r += dot(x[i], Aij, y[j]) + dot(x[j], adjoint(Aij), y[i])
+            end
         end
-        return r
-    else
-        r = dot(x, Δ, y)
-        r += conj(dot(y, Δ, x))
-        @inbounds @simd for i in eachindex(x)
-            r -= dot(x[i], A[i,i], y[i])
+    else # A.uplo == 'L'
+        @inbounds for j = 1:length(y)
+            r += dot(x[j], real(data[j,j]), y[j])
+            @simd for i = j+1:length(y)
+                Aij = data[i,j]
+                r += dot(x[i], Aij, y[j]) + dot(x[j], adjoint(Aij), y[i])
+            end
         end
-        return r
     end
+    return r
 end
 
 # Fallbacks to avoid generic_matvecmul!/generic_matmatmul!

From 200732c90aca493fcf941d0472bdebb1df151175 Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@gmx.de>
Date: Fri, 30 Aug 2019 16:44:52 +0200
Subject: [PATCH 34/35] test quaternions in uniformscaling gendot

---
 stdlib/LinearAlgebra/test/uniformscaling.jl | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/stdlib/LinearAlgebra/test/uniformscaling.jl b/stdlib/LinearAlgebra/test/uniformscaling.jl
index fd0c80e90326d..141594f3ef247 100644
--- a/stdlib/LinearAlgebra/test/uniformscaling.jl
+++ b/stdlib/LinearAlgebra/test/uniformscaling.jl
@@ -4,6 +4,10 @@ module TestUniformscaling
 
 using Test, LinearAlgebra, Random, SparseArrays
 
+const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test")
+isdefined(Main, :Quaternions) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "Quaternions.jl"))
+using .Main.Quaternions
+
 Random.seed!(123)
 
 @testset "basic functions" begin
@@ -336,6 +340,9 @@ end
     λ = rand(-10:10)
     J = UniformScaling(λ)
     @test dot(x, J, y) == λ*dot(x, y)
+    λ = Quaternion(0.44567, 0.755871, 0.882548, 0.423612)
+    x, y = Quaternion(rand(4)...), Quaternion(rand(4)...)
+    @test dot([x], λ*I, [y]) ≈ dot(x, λ, y) ≈ dot(x, λ*I, y) ≈ dot(x, λ*y)
 end
 
 end # module TestUniformscaling

From d30426458f6de246326b8b450ac416f95e7819b1 Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <Daniel.Karrasch@gmx.de>
Date: Fri, 30 Aug 2019 17:55:31 +0200
Subject: [PATCH 35/35] fix uniformscaling test

---
 stdlib/LinearAlgebra/test/uniformscaling.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/stdlib/LinearAlgebra/test/uniformscaling.jl b/stdlib/LinearAlgebra/test/uniformscaling.jl
index 141594f3ef247..63096fac75124 100644
--- a/stdlib/LinearAlgebra/test/uniformscaling.jl
+++ b/stdlib/LinearAlgebra/test/uniformscaling.jl
@@ -342,7 +342,7 @@ end
     @test dot(x, J, y) == λ*dot(x, y)
     λ = Quaternion(0.44567, 0.755871, 0.882548, 0.423612)
     x, y = Quaternion(rand(4)...), Quaternion(rand(4)...)
-    @test dot([x], λ*I, [y]) ≈ dot(x, λ, y) ≈ dot(x, λ*I, y) ≈ dot(x, λ*y)
+    @test dot([x], λ*I, [y]) ≈ dot(x, λ, y) ≈ dot(x, λ*y)
 end
 
 end # module TestUniformscaling