diff --git a/stdlib/SparseArrays/src/linalg.jl b/stdlib/SparseArrays/src/linalg.jl
index c54a6c760c0fd..75587b203f6af 100644
--- a/stdlib/SparseArrays/src/linalg.jl
+++ b/stdlib/SparseArrays/src/linalg.jl
@@ -44,9 +44,9 @@ function mul!(C::StridedVecOrMat, A::AbstractSparseMatrixCSC, B::Union{StridedVe
     end
     C
 end
-*(A::AbstractSparseMatrixCSC{TA}, x::StridedVector{Tx}) where {TA,Tx} =
+*(A::SparseMatrixCSCUnion{TA}, x::StridedVector{Tx}) where {TA,Tx} =
     (T = promote_op(matprod, TA, Tx); mul!(similar(x, T, size(A, 1)), A, x, true, false))
-*(A::AbstractSparseMatrixCSC{TA}, B::AdjOrTransStridedOrTriangularMatrix{Tx}) where {TA,Tx} =
+*(A::SparseMatrixCSCUnion{TA}, B::AdjOrTransStridedOrTriangularMatrix{Tx}) where {TA,Tx} =
     (T = promote_op(matprod, TA, Tx); mul!(similar(B, T, (size(A, 1), size(B, 2))), A, B, true, false))
 
 function mul!(C::StridedVecOrMat, adjA::Adjoint{<:Any,<:AbstractSparseMatrixCSC}, B::Union{StridedVector,AdjOrTransStridedOrTriangularMatrix}, α::Number, β::Number)
@@ -125,7 +125,7 @@ function mul!(C::StridedVecOrMat, X::AdjOrTransStridedOrTriangularMatrix, A::Abs
     end
     C
 end
-*(X::AdjOrTransStridedOrTriangularMatrix, A::AbstractSparseMatrixCSC{TvA}) where {TvA} =
+*(X::AdjOrTransStridedOrTriangularMatrix, A::SparseMatrixCSCUnion{TvA}) where {TvA} =
     (T = promote_op(matprod, eltype(X), TvA); mul!(similar(X, T, (size(X, 1), size(A, 2))), X, A, true, false))
 
 function mul!(C::StridedVecOrMat, X::AdjOrTransStridedOrTriangularMatrix, adjA::Adjoint{<:Any,<:AbstractSparseMatrixCSC}, α::Number, β::Number)
@@ -182,11 +182,20 @@ end
 # Sparse matrix multiplication as described in [Gustavson, 1978]:
 # http://dl.acm.org/citation.cfm?id=355796
 
-*(A::AbstractSparseMatrixCSC, B::AbstractSparseMatrixCSC) = spmatmul(A,B)
-*(A::AbstractSparseMatrixCSC, B::Adjoint{<:Any,<:AbstractSparseMatrixCSC}) = spmatmul(A, copy(B))
-*(A::AbstractSparseMatrixCSC, B::Transpose{<:Any,<:AbstractSparseMatrixCSC}) = spmatmul(A, copy(B))
-*(A::Transpose{<:Any,<:AbstractSparseMatrixCSC}, B::AbstractSparseMatrixCSC) = spmatmul(copy(A), B)
-*(A::Adjoint{<:Any,<:AbstractSparseMatrixCSC}, B::AbstractSparseMatrixCSC) = spmatmul(copy(A), B)
+const SparseTriangular{Tv,Ti} = Union{UpperTriangular{Tv,<:SparseMatrixCSCUnion{Tv,Ti}},LowerTriangular{Tv,<:SparseMatrixCSCUnion{Tv,Ti}}}
+const SparseOrTri{Tv,Ti} = Union{SparseMatrixCSCUnion{Tv,Ti},SparseTriangular{Tv,Ti}}
+
+*(A::SparseOrTri, B::AbstractSparseVector) = spmatmulv(A, B)
+*(A::SparseOrTri, B::SparseColumnView) = spmatmulv(A, B)
+*(A::SparseOrTri, B::SparseVectorView) = spmatmulv(A, B)
+*(A::SparseMatrixCSCUnion, B::SparseMatrixCSCUnion) = spmatmul(A,B)
+*(A::SparseTriangular, B::SparseMatrixCSCUnion) = spmatmul(A,B)
+*(A::SparseMatrixCSCUnion, B::SparseTriangular) = spmatmul(A,B)
+*(A::SparseTriangular, B::SparseTriangular) = spmatmul1(A,B)
+*(A::SparseOrTri, B::Adjoint{<:Any,<:AbstractSparseMatrixCSC}) = spmatmul(A, copy(B))
+*(A::SparseOrTri, B::Transpose{<:Any,<:AbstractSparseMatrixCSC}) = spmatmul(A, copy(B))
+*(A::Transpose{<:Any,<:AbstractSparseMatrixCSC}, B::SparseOrTri) = spmatmul(copy(A), B)
+*(A::Adjoint{<:Any,<:AbstractSparseMatrixCSC}, B::SparseOrTri) = spmatmul(copy(A), B)
 *(A::Adjoint{<:Any,<:AbstractSparseMatrixCSC}, B::Adjoint{<:Any,<:AbstractSparseMatrixCSC}) = spmatmul(copy(A), copy(B))
 *(A::Transpose{<:Any,<:AbstractSparseMatrixCSC}, B::Transpose{<:Any,<:AbstractSparseMatrixCSC}) = spmatmul(copy(A), copy(B))
 
@@ -198,16 +207,15 @@ end
 # done by a quicksort of the row indices or by a full scan of the dense result vector.
 # The last is faster, if more than ≈ 1/32 of the result column is nonzero.
 # TODO: extend to SparseMatrixCSCUnion to allow for SubArrays (view(X, :, r)).
-function spmatmul(A::AbstractSparseMatrixCSC{TvA,TiA}, B::AbstractSparseMatrixCSC{TvB,TiB}) where {TvA,TiA,TvB,TiB}
+function spmatmul(A::SparseOrTri{TvA,TiA}, B::Union{SparseOrTri{TvB,TiB},SparseVectorUnion{TvB,TiB},SubArray{TvB,<:Any,<:AbstractSparseArray{TvB,TiB}}}) where {TvA,TiA,TvB,TiB}
+
     Tv = promote_op(matprod, TvA, TvB)
     Ti = promote_type(TiA, TiB)
     mA, nA = size(A)
     nB = size(B, 2)
     nA == size(B, 1) || throw(DimensionMismatch())
 
-    rowvalA = rowvals(A); nzvalA = nonzeros(A)
-    rowvalB = rowvals(B); nzvalB = nonzeros(B)
-    nnzC = max(estimate_mulsize(mA, nnz(A), nA, nnz(B), nB) * 11 ÷ 10, mA)
+    nnzC = min(estimate_mulsize(mA, nnz(A), nA, nnz(B), nB) * 11 ÷ 10 + mA, mA*nB)
     colptrC = Vector{Ti}(undef, nB+1)
     rowvalC = Vector{Ti}(undef, nnzC)
     nzvalC = Vector{Tv}(undef, nnzC)
@@ -221,45 +229,8 @@ function spmatmul(A::AbstractSparseMatrixCSC{TvA,TiA}, B::AbstractSparseMatrixCS
                 resize!(rowvalC, nnzC)
                 resize!(nzvalC, nnzC)
             end
-            colptrC[i] = ip0 = ip
-            k0 = ip - 1
-            for jp in nzrange(B, i)
-                nzB = nzvalB[jp]
-                j = rowvalB[jp]
-                for kp in nzrange(A, j)
-                    nzC = nzvalA[kp] * nzB
-                    k = rowvalA[kp]
-                    if xb[k]
-                        nzvalC[k+k0] += nzC
-                    else
-                        nzvalC[k+k0] = nzC
-                        xb[k] = true
-                        rowvalC[ip] = k
-                        ip += 1
-                    end
-                end
-            end
-            if ip > ip0
-                if prefer_sort(ip-k0, mA)
-                    # in-place sort of indices. Effort: O(nnz*ln(nnz)).
-                    sort!(rowvalC, ip0, ip-1, QuickSort, Base.Order.Forward)
-                    for vp = ip0:ip-1
-                        k = rowvalC[vp]
-                        xb[k] = false
-                        nzvalC[vp] = nzvalC[k+k0]
-                    end
-                else
-                    # scan result vector (effort O(mA))
-                    for k = 1:mA
-                        if xb[k]
-                            xb[k] = false
-                            rowvalC[ip0] = k
-                            nzvalC[ip0] = nzvalC[k+k0]
-                            ip0 += 1
-                        end
-                    end
-                end
-            end
+            colptrC[i] = ip
+            ip = spcolmul!(rowvalC, nzvalC, xb, i, ip, A, B)
         end
         colptrC[nB+1] = ip
     end
@@ -272,6 +243,68 @@ function spmatmul(A::AbstractSparseMatrixCSC{TvA,TiA}, B::AbstractSparseMatrixCS
     return C
 end
 
+# process single rhs column
+function spcolmul!(rowvalC, nzvalC, xb, i, ip, A, B)
+    rowvalA = rowvals(A); nzvalA = nonzeros(A)
+    rowvalB = rowvals(B); nzvalB = nonzeros(B)
+    mA = size(A, 1)
+    ip0 = ip
+    k0 = ip - 1
+    @inbounds begin
+        for jp in nzrange(B, i)
+            nzB = nzvalB[jp]
+            j = rowvalB[jp]
+            for kp in nzrange(A, j)
+                nzC = nzvalA[kp] * nzB
+                k = rowvalA[kp]
+                if xb[k]
+                    nzvalC[k+k0] += nzC
+                else
+                    nzvalC[k+k0] = nzC
+                    xb[k] = true
+                    rowvalC[ip] = k
+                    ip += 1
+                end
+            end
+        end
+        if ip > ip0
+            if prefer_sort(ip-k0, mA)
+                # in-place sort of indices. Effort: O(nnz*ln(nnz)).
+                sort!(rowvalC, ip0, ip-1, QuickSort, Base.Order.Forward)
+                for vp = ip0:ip-1
+                    k = rowvalC[vp]
+                    xb[k] = false
+                    nzvalC[vp] = nzvalC[k+k0]
+                end
+            else
+                # scan result vector (effort O(mA))
+                for k = 1:mA
+                    if xb[k]
+                        xb[k] = false
+                        rowvalC[ip0] = k
+                        nzvalC[ip0] = nzvalC[k+k0]
+                        ip0 += 1
+                    end
+                end
+            end
+        end
+    end
+    return ip
+end
+
+# special cases of same twin Upper/LowerTriangular
+spmatmul1(A, B) = spmatmul(A, B)
+function spmatmul1(A::UpperTriangular, B::UpperTriangular)
+    UpperTriangular(spmatmul(A, B))
+end
+function spmatmul1(A::LowerTriangular, B::LowerTriangular)
+    LowerTriangular(spmatmul(A, B))
+end
+# exploit spmatmul for sparse vectors and column views
+function spmatmulv(A, B)
+    spmatmul(A, B)[:,1]
+end
+
 # estimated number of non-zeros in matrix product
 # it is assumed, that the non-zero indices are distributed independently and uniformly
 # in both matrices. Over-estimation is possible if that is not the case.
@@ -765,7 +798,7 @@ function _ldiv!(U::UpperTriangularWrapped, B::StridedVecOrMat)
 end
 
 (\)(L::TriangularSparse, B::AbstractSparseMatrixCSC) = ldiv!(L, Array(B))
-(*)(L::TriangularSparse, B::AbstractSparseMatrixCSC) = lmul!(L, Array(B))
+#(*)(L::TriangularSparse, B::AbstractSparseMatrixCSC) = lmul!(L, Array(B))
 
 ## end of triangular
 
diff --git a/stdlib/SparseArrays/src/sparsematrix.jl b/stdlib/SparseArrays/src/sparsematrix.jl
index 6cb72428247e3..78ad80101dd8e 100644
--- a/stdlib/SparseArrays/src/sparsematrix.jl
+++ b/stdlib/SparseArrays/src/sparsematrix.jl
@@ -109,8 +109,15 @@ julia> nnz(A)
 ```
 """
 nnz(S::AbstractSparseMatrixCSC) = Int(getcolptr(S)[size(S, 2) + 1] - 1)
-nnz(S::ReshapedArray{T,1,<:AbstractSparseMatrixCSC}) where T = nnz(parent(S))
-count(pred, S::AbstractSparseMatrixCSC) = count(pred, nzvalview(S)) + pred(zero(eltype(S)))*(prod(size(S)) - nnz(S))
+nnz(S::ReshapedArray{<:Any,1,<:AbstractSparseMatrixCSC}) = nnz(parent(S))
+nnz(S::UpperTriangular{<:Any,<:AbstractSparseMatrixCSC}) = nnz1(S)
+nnz(S::LowerTriangular{<:Any,<:AbstractSparseMatrixCSC}) = nnz1(S)
+nnz(S::SparseMatrixCSCView) = nnz1(S)
+nnz1(S) = sum(length.(nzrange.(Ref(S), axes(S, 2))))
+
+function count(pred, S::AbstractSparseMatrixCSC)
+    count(pred, nzvalview(S)) + pred(zero(eltype(S)))*(prod(size(S)) - nnz(S))
+end
 
 """
     nonzeros(A)
@@ -138,6 +145,8 @@ julia> nonzeros(A)
 """
 nonzeros(S::SparseMatrixCSC) = getfield(S, :nzval)
 nonzeros(S::SparseMatrixCSCView)  = nonzeros(S.parent)
+nonzeros(S::UpperTriangular{<:Any,<:SparseMatrixCSCUnion}) = nonzeros(S.data)
+nonzeros(S::LowerTriangular{<:Any,<:SparseMatrixCSCUnion}) = nonzeros(S.data)
 
 """
     rowvals(A::AbstractSparseMatrixCSC)
@@ -164,6 +173,8 @@ julia> rowvals(A)
 """
 rowvals(S::SparseMatrixCSC) = getfield(S, :rowval)
 rowvals(S::SparseMatrixCSCView) = rowvals(S.parent)
+rowvals(S::UpperTriangular{<:Any,<:SparseMatrixCSCUnion}) = rowvals(S.data)
+rowvals(S::LowerTriangular{<:Any,<:SparseMatrixCSCUnion}) = rowvals(S.data)
 
 """
     nzrange(A::AbstractSparseMatrixCSC, col::Integer)
@@ -186,6 +197,8 @@ column. In conjunction with [`nonzeros`](@ref) and
 """
 nzrange(S::AbstractSparseMatrixCSC, col::Integer) = getcolptr(S)[col]:(getcolptr(S)[col+1]-1)
 nzrange(S::SparseMatrixCSCView, col::Integer) = nzrange(S.parent, S.indices[2][col])
+nzrange(S::UpperTriangular{<:Any,<:SparseMatrixCSCUnion}, i::Integer) = nzrangeup(S.data, i)
+nzrange(S::LowerTriangular{<:Any,<:SparseMatrixCSCUnion}, i::Integer) = nzrangelo(S.data, i)
 
 function Base.isstored(A::AbstractSparseMatrixCSC, i::Integer, j::Integer)
     @boundscheck checkbounds(A, i, j)
diff --git a/stdlib/SparseArrays/src/sparsevector.jl b/stdlib/SparseArrays/src/sparsevector.jl
index 2cc9268e6fea8..077bb11583330 100644
--- a/stdlib/SparseArrays/src/sparsevector.jl
+++ b/stdlib/SparseArrays/src/sparsevector.jl
@@ -32,16 +32,35 @@ SparseVector(n::Integer, nzind::Vector{Ti}, nzval::Vector{Tv}) where {Tv,Ti} =
 
 # Define an alias for a view of a whole column of a SparseMatrixCSC. Many methods can be written for the
 # union of such a view and a SparseVector so we define an alias for such a union as well
-const SparseColumnView{T}  = SubArray{T,1,<:AbstractSparseMatrixCSC,Tuple{Base.Slice{Base.OneTo{Int}},Int},false}
-const SparseVectorUnion{T} = Union{SparseVector{T}, SparseColumnView{T}}
-const AdjOrTransSparseVectorUnion{T} = LinearAlgebra.AdjOrTrans{T, <:SparseVectorUnion{T}}
+const SparseColumnView{Tv,Ti}  = SubArray{Tv,1,<:AbstractSparseMatrixCSC{Tv,Ti},Tuple{Base.Slice{Base.OneTo{Int}},Int},false}
+const SparseVectorView{Tv,Ti}  = SubArray{Tv,1,<:AbstractSparseVector{Tv,Ti},Tuple{Base.Slice{Base.OneTo{Int}}},false}
+const SparseVectorUnion{Tv,Ti} = Union{SparseVector{Tv,Ti}, SparseColumnView{Tv,Ti}, SparseVectorView{Tv,Ti}}
+const AdjOrTransSparseVectorUnion{Tv,Ti} = LinearAlgebra.AdjOrTrans{Tv, <:SparseVectorUnion{Tv,Ti}}
 
 ### Basic properties
 
 size(x::SparseVector)     = (getfield(x, :n),)
-nnz(x::SparseVector)      = length(nonzeros(x))
 count(f, x::SparseVector) = count(f, nonzeros(x)) + f(zero(eltype(x)))*(length(x) - nnz(x))
 
+# implement the nnz - nzrange - nonzeros - rowvals interface for sparse vectors
+
+nnz(x::SparseVector)      = length(nonzeros(x))
+function nnz(x::SparseColumnView)
+    rowidx, colidx = parentindices(x)
+    return length(nzrange(parent(x), colidx))
+end
+nnz(x::SparseVectorView) = nnz(x.parent)
+
+"""
+    nzrange(x::SparseVectorUnion, col)
+
+Give the range of indices to the structural nonzero values of a sparse vector.
+The column index `col` is ignored (assumed to be `1`).
+"""
+function nzrange(x::SparseVectorUnion, j::Integer)
+    j == 1 ? (1:nnz(x)) : throw(BoundsError(x, (":", j)))
+end
+
 nonzeros(x::SparseVector) = getfield(x, :nzval)
 function nonzeros(x::SparseColumnView)
     rowidx, colidx = parentindices(x)
@@ -49,6 +68,7 @@ function nonzeros(x::SparseColumnView)
     @inbounds y = view(nonzeros(A), nzrange(A, colidx))
     return y
 end
+nonzeros(x::SparseVectorView) = nonzeros(parent(x))
 
 nonzeroinds(x::SparseVector) = getfield(x, :nzind)
 function nonzeroinds(x::SparseColumnView)
@@ -57,12 +77,12 @@ function nonzeroinds(x::SparseColumnView)
     @inbounds y = view(rowvals(A), nzrange(A, colidx))
     return y
 end
+nonzeroinds(x::SparseVectorView) = nonzeroinds(parent(x))
+
+rowvals(x::SparseVectorUnion) = nonzeroinds(x)
 
 indtype(x::SparseColumnView) = indtype(parent(x))
-function nnz(x::SparseColumnView)
-    rowidx, colidx = parentindices(x)
-    return length(nzrange(parent(x), colidx))
-end
+indtype(x::SparseVectorView) = indtype(parent(x))
 
 ## similar
 #
diff --git a/stdlib/SparseArrays/test/sparse.jl b/stdlib/SparseArrays/test/sparse.jl
index a5d71f0bcce7a..b6a747c74b4b5 100644
--- a/stdlib/SparseArrays/test/sparse.jl
+++ b/stdlib/SparseArrays/test/sparse.jl
@@ -2996,4 +2996,84 @@ end
     end
 end
 
+@testset "Multiplying with triangular sparse matrices #35609 #35610" begin
+    n = 10
+    A = sprand(n, n, 5/n)
+    U = UpperTriangular(A)
+    L = LowerTriangular(A)
+    AM = Matrix(A)
+    UM = Matrix(U)
+    LM = Matrix(L)
+    Y = A * U
+    @test Y ≈ AM * UM
+    @test typeof(Y) == typeof(A)
+    Y = A * L
+    @test Y ≈ AM * LM
+    @test typeof(Y) == typeof(A)
+    Y = U * A
+    @test Y ≈ UM * AM
+    @test typeof(Y) == typeof(A)
+    Y = L * A
+    @test Y ≈ LM * AM
+    @test typeof(Y) == typeof(A)
+    Y = U * U
+    @test Y ≈ UM * UM
+    @test typeof(Y) == typeof(U)
+    Y = L * L
+    @test Y ≈ LM * LM
+    @test typeof(Y) == typeof(L)
+    Y = L * U
+    @test Y ≈ LM * UM
+    @test typeof(Y) == typeof(A)
+    Y = U * L
+    @test Y ≈ UM * LM
+    @test typeof(Y) == typeof(A)
+end
+
+#testing the sparse matrix/vector access functions nnz, nzrange, rowvals, nonzeros
+@testset "generic sparse matrix access functions" begin
+    I = [1,3,4,5, 1,3,4,5, 1,3,4,5];
+    J = [4,4,4,4, 5,5,5,5, 6,6,6,6];
+    V = [14,34,44,54, 15,35,45,55, 16,36,46,56];
+    A = sparse(I, J, V, 9, 9);
+    AU = UpperTriangular(A)
+    AL = LowerTriangular(A)
+    b = SparseVector(9, I[1:4], V[1:4])
+    c = view(A, :, 5)
+    d = view(b, :)
+
+    @testset "nnz $n" for (n, M, nz) in (("A", A, 12), ("AU", AU, 11), ("AL", AL, 3),
+                                         ("b", b, 4), ("c", c, 4), ("d", d, 4))
+        @test nnz(M) == nz
+        @test_throws BoundsError nzrange(M, 0)
+        @test_throws BoundsError nzrange(M, size(M, 2) + 1)
+    end
+    @testset "nzrange(A, $i)" for (i, nzr) in ((1,1:0),(4,1:4),(5,5:8),(6,9:12),(9,13:12))
+        @test nzrange(A, i) == nzr
+    end
+    @testset "nzrange(AU, $i)" for (i, nzr) in ((2,1:0),(4,1:3),(5,5:8),(6,9:12),(8,13:12))
+        @test nzrange(AU, i) == nzr
+    end
+    @testset "nzrange(AL, $i)" for (i, nzr) in ((3,1:0),(4,3:4),(5,8:8),(6,13:12),(7,13:12))
+        @test nzrange(AL, i) == nzr
+    end
+    @test nzrange(b, 1) == 1:4
+    @test nzrange(c, 1) == 1:4
+    @test nzrange(d, 1) == 1:4
+
+    @test rowvals(A) == I
+    @test rowvals(AL) == I
+    @test rowvals(AL) == I
+    @test rowvals(b) == I[1:4]
+    @test rowvals(c) == I[5:8]
+    @test rowvals(d) == I[1:4]
+
+    @test nonzeros(A) == V
+    @test nonzeros(AU) == V
+    @test nonzeros(AL) == V
+    @test nonzeros(b) == V[1:4]
+    @test nonzeros(c) == V[5:8]
+    @test nonzeros(d) == V[1:4]
+end
+
 end # module
diff --git a/stdlib/SparseArrays/test/sparsevector.jl b/stdlib/SparseArrays/test/sparsevector.jl
index 975a9f2885476..b69db8fdb1e2c 100644
--- a/stdlib/SparseArrays/test/sparsevector.jl
+++ b/stdlib/SparseArrays/test/sparsevector.jl
@@ -1442,4 +1442,19 @@ end
     @test nonzeros(A) !== nonzeros(B)
 end
 
+@testset "multiplication of Triangular sparse matrices with sparse vectors #35642" begin
+    n = 10
+    A = sprand(n, n, 5/n)
+    U = UpperTriangular(A)
+    L = LowerTriangular(A)
+    x = sprand(n, 5/n)
+    y = view(A, :, 6)
+    z = view(x, :)
+    ty = typeof
+    @testset "matvec multiplication $(ty(X)) * $(ty(v))" for X in (U, L), v in (x, y, z)
+        @test X * v ≈ Matrix(X) * Vector(v)
+        @test typeof(X * v) == typeof(x)
+    end
+end
+
 end # module