diff --git a/src/SummationByPartsOperators.jl b/src/SummationByPartsOperators.jl
index cd07ccae..7dc9e015 100644
--- a/src/SummationByPartsOperators.jl
+++ b/src/SummationByPartsOperators.jl
@@ -137,7 +137,7 @@ export FilterCallback, ConstantFilter, ExponentialFilter
 export SafeMode, FastMode, ThreadedMode
 export derivative_order, accuracy_order, source_of_coefficients, grid, semidiscretize
 export mass_matrix
-export integrate, left_boundary_weight, right_boundary_weight,
+export integrate, left_boundary_weight, right_boundary_weight, weight,
        derivative_left, derivative_right,
        mul_transpose_derivative_left!, mul_transpose_derivative_right!,
        evaluate_coefficients, evaluate_coefficients!,
diff --git a/src/general_operators.jl b/src/general_operators.jl
index d02a0d35..cc8d9878 100644
--- a/src/general_operators.jl
+++ b/src/general_operators.jl
@@ -56,7 +56,10 @@ source_of_coefficients(D) = SummationByPartsOperators
     left_boundary_weight(D)
 
 Return the left-boundary weight of the (diagonal) mass matrix `M` associated to
-the derivative operator `D`.
+the SBP derivative operator `D`.
+
+See also [`right_boundary_weight`](@ref), [`weight`](@ref), and
+[`mass_matrix`](@ref).
 """
 function left_boundary_weight end
 
@@ -64,10 +67,35 @@ function left_boundary_weight end
     right_boundary_weight(D)
 
 Return the left-boundary weight of the (diagonal) mass matrix `M` associated to
-the derivative operator `D`.
+the SBP derivative operator `D`.
+
+See also [`left_boundary_weight`](@ref), [`weight`](@ref), and
+[`mass_matrix`](@ref).
 """
 function right_boundary_weight end
 
+"""
+    weight(D, i::Int)
+
+Return the `i`th weight of the (diagonal) mass matrix `M` associated to
+the SBP derivative operator `D` (starting to count at `1`).
+
+See also [`left_boundary_weight`](@ref), [`right_boundary_weight`](@ref), and
+[`mass_matrix`](@ref).
+"""
+function weight end
+
+"""
+    mass_matrix(D)
+
+Return the (diagonal) mass matrix `M` associated to the SBP derivative operator
+`D`.
+
+See also [`left_boundary_weight`](@ref), [`right_boundary_weight`](@ref), and
+[`weight`](@ref).
+"""
+function mass_matrix end
+
 function Base.summary(io::IO, D::AbstractDerivativeOperator)
     print(io, nameof(typeof(D)), "(derivative:", derivative_order(D),
               ", accuracy:", accuracy_order(D), ")")
diff --git a/src/var_coef_operators.jl b/src/var_coef_operators.jl
index aa58240a..8fb6628c 100644
--- a/src/var_coef_operators.jl
+++ b/src/var_coef_operators.jl
@@ -94,11 +94,6 @@ end
 
 
 
-"""
-    mass_matrix(D::Union{DerivativeOperator,VarCoefDerivativeOperator})
-
-Create the diagonal mass matrix for the SBP derivative operator `D`.
-"""
 function mass_matrix(D::Union{DerivativeOperator,VarCoefDerivativeOperator})
     m = fill(one(eltype(D)), length(grid(D)))
     @unpack left_weights, right_weights = D.coefficients
diff --git a/test/coupling_test.jl b/test/coupling_test.jl
index 5ce95e06..ef330e6b 100644
--- a/test/coupling_test.jl
+++ b/test/coupling_test.jl
@@ -35,8 +35,8 @@ using SummationByPartsOperators
           SummationByPartsOperators.scale_by_mass_matrix!(u, cD)
           SummationByPartsOperators.scale_by_inverse_mass_matrix!(u, cD)
           @test u ≈ v
-          @test integrate(u, cD) ≈ sum(mass_matrix(cD) * u)
-          @test integrate(u->u^2, u, cD) ≈ sum(u' * mass_matrix(cD) * u)
+          @test integrate(u, cD) ≈ sum(mass_matrix(cD) * u) ≈ sum(weight(cD, i) * u[i] for i in eachindex(u))
+          @test integrate(u->u^2, u, cD) ≈ sum(u' * mass_matrix(cD) * u) ≈ sum(weight(cD, i) * u[i]^2 for i in eachindex(u))
           @test cD * u ≈ BandedMatrix(cD) * u atol=eps(float(T)) rtol=sqrt(eps(float(T)))
           @test cD * u ≈ Matrix(cD) * u atol=eps(float(T)) rtol=sqrt(eps(float(T)))
           @test cD * u ≈ sparse(cD) * u atol=eps(float(T)) rtol=sqrt(eps(float(T)))
@@ -55,6 +55,9 @@ using SummationByPartsOperators
         M = mass_matrix(cD_central)
         @test M ≈ mass_matrix(cD_plus)
         @test M ≈ mass_matrix(cD_minus)
+        for i in 1:size(cD_central, 1)
+          @test @inferred(weight(cD_central, i)) ≈ @inferred(weight(cD_plus, i)) ≈ @inferred(weight(cD_minus, i))
+        end
         @test sum(M) ≈ xmax - xmin
         Dp = Matrix(cD_plus)
         Dm = Matrix(cD_minus)
@@ -75,9 +78,13 @@ using SummationByPartsOperators
         cD2 = couple_discontinuously(D, UniformMesh1D(xmin, xmax, N^2))
         @test grid(cD1) ≈ grid(cD2)
         @test mass_matrix(cD1) ≈ mass_matrix(cD2)
+        for i in 1:size(cD1, 1)
+          @test @inferred(weight(cD1, i)) ≈ @inferred(weight(cD2, i))
+        end
         @test Matrix(cD1) ≈ Matrix(cD2)
 
         Mcont = mass_matrix(cD_continuous)
+        @test Mcont ≈ Diagonal([weight(cD_continuous, i) for i in 1:size(cD_continuous, 1)])
         @test sum(Mcont) ≈ xmax - xmin
         Dcont = Matrix(cD_continuous)
         res = Mcont * Dcont + Dcont' * Mcont
@@ -109,6 +116,7 @@ using SummationByPartsOperators
           @test u ≈ v
           @test integrate(u, cD) ≈ sum(mass_matrix(cD) * u)
           @test integrate(u->u^2, u, cD) ≈ sum(u' * mass_matrix(cD) * u)
+          @test mass_matrix(cD) ≈ Diagonal([weight(cD, i) for i in 1:size(cD, 1)])
           @test cD * u ≈ Matrix(cD) * u atol=eps(float(T)) rtol=sqrt(eps(float(T)))
           @test cD * u ≈ sparse(cD) * u atol=eps(float(T)) rtol=sqrt(eps(float(T)))
           @test SummationByPartsOperators.xmin(cD) ≈ xmin
@@ -117,6 +125,9 @@ using SummationByPartsOperators
         M = mass_matrix(cD_central)
         @test M ≈ mass_matrix(cD_plus)
         @test M ≈ mass_matrix(cD_minus)
+        for i in 1:size(cD_central, 1)
+          @test @inferred(weight(cD_central, i)) ≈ @inferred(weight(cD_plus, i)) ≈ @inferred(weight(cD_minus, i))
+        end
         @test sum(M) ≈ xmax - xmin
         Dp = Matrix(cD_plus)
         Dm = Matrix(cD_minus)
@@ -130,6 +141,7 @@ using SummationByPartsOperators
         @test norm(res) < 100N * eps(float(T))
 
         Mcont = mass_matrix(cD_continuous)
+        @test Mcont ≈ Diagonal([weight(cD_continuous, i) for i in 1:size(cD_continuous, 1)])
         @test sum(Mcont) ≈ xmax - xmin
         Dcont = Matrix(cD_continuous)
         res = Mcont * Dcont + Dcont' * Mcont
@@ -172,6 +184,7 @@ end
             @test u ≈ v
             @test integrate(u, cD) ≈ sum(mass_matrix(cD) * u)
             @test integrate(u->u^2, u, cD) ≈ sum(u' * mass_matrix(cD) * u)
+            @test mass_matrix(cD) ≈ Diagonal([weight(cD, i) for i in 1:size(cD, 1)])
             @test cD * u ≈ BandedMatrix(cD) * u atol=eps(float(T)) rtol=sqrt(eps(float(T)))
             @test cD * u ≈ Matrix(cD) * u atol=eps(float(T)) rtol=sqrt(eps(float(T)))
             @test cD * u ≈ sparse(cD) * u atol=eps(float(T)) rtol=sqrt(eps(float(T)))
@@ -190,6 +203,9 @@ end
           M = mass_matrix(cD_central)
           @test M ≈ mass_matrix(cD_plus)
           @test M ≈ mass_matrix(cD_minus)
+          for i in 1:size(cD_central, 1)
+            @test @inferred(weight(cD_central, i)) ≈ @inferred(weight(cD_plus, i)) ≈ @inferred(weight(cD_minus, i))
+          end
           @test sum(M) ≈ xmax - xmin
           Dp = Matrix(cD_plus)
           Dm = Matrix(cD_minus)
@@ -213,6 +229,7 @@ end
           @test Matrix(cD1) ≈ Matrix(cD2)
 
           Mcont = mass_matrix(cD_continuous)
+          @test Mcont ≈ Diagonal([weight(cD_continuous, i) for i in 1:size(cD_continuous, 1)])
           @test sum(Mcont) ≈ xmax - xmin
           Dcont = Matrix(cD_continuous)
           res = Mcont * Dcont + Dcont' * Mcont
@@ -242,6 +259,7 @@ end
             @test u ≈ v
             @test integrate(u, cD) ≈ sum(mass_matrix(cD) * u)
             @test integrate(u->u^2, u, cD) ≈ sum(u' * mass_matrix(cD) * u)
+            @test mass_matrix(cD) ≈ Diagonal([weight(cD, i) for i in 1:size(cD, 1)])
             @test cD * u ≈ Matrix(cD) * u atol=eps(float(T)) rtol=sqrt(eps(float(T)))
             @test cD * u ≈ sparse(cD) * u atol=eps(float(T)) rtol=sqrt(eps(float(T)))
             @test SummationByPartsOperators.xmin(cD) ≈ xmin
@@ -250,6 +268,9 @@ end
           M = mass_matrix(cD_central)
           @test M ≈ mass_matrix(cD_plus)
           @test M ≈ mass_matrix(cD_minus)
+          for i in 1:size(cD_central, 1)
+            @test @inferred(weight(cD_central, i)) ≈ @inferred(weight(cD_plus, i)) ≈ @inferred(weight(cD_minus, i))
+          end
           @test sum(M) ≈ xmax - xmin
           Dp = Matrix(cD_plus)
           Dm = Matrix(cD_minus)
@@ -263,6 +284,7 @@ end
           @test norm(res) < 100N * eps(float(T))
 
           Mcont = mass_matrix(cD_continuous)
+          @test Mcont ≈ Diagonal([weight(cD_continuous, i) for i in 1:size(cD_continuous, 1)])
           @test sum(Mcont) ≈ xmax - xmin
           Dcont = Matrix(cD_continuous)
           res = Mcont * Dcont + Dcont' * Mcont
@@ -306,6 +328,7 @@ end
           @test u ≈ v
           @test integrate(u, cD) ≈ sum(mass_matrix(cD) * u)
           @test integrate(u->u^2, u, cD) ≈ sum(u' * mass_matrix(cD) * u)
+          @test mass_matrix(cD) ≈ Diagonal([weight(cD, i) for i in 1:size(cD, 1)])
           @test cD * u ≈ BandedMatrix(cD) * u atol=eps(float(T)) rtol=sqrt(eps(float(T)))
           @test cD * u ≈ Matrix(cD) * u atol=eps(float(T)) rtol=sqrt(eps(float(T)))
           @test cD * u ≈ sparse(cD) * u atol=eps(float(T)) rtol=sqrt(eps(float(T)))
@@ -317,6 +340,9 @@ end
           cDm_dense = Matrix(cDm)
           M = mass_matrix(cDp)
           @test M ≈ mass_matrix(cDm)
+          for i in 1:size(cD_central, 1)
+            @test @inferred(weight(cDp, i)) ≈ @inferred(weight(cDm, i))
+          end
           @test sum(M) ≈ xmax - xmin
 
           diss = M * (cDp_dense - cDm_dense)
@@ -353,6 +379,7 @@ end
           @test u ≈ v
           @test integrate(u, cD) ≈ sum(mass_matrix(cD) * u)
           @test integrate(u->u^2, u, cD) ≈ sum(u' * mass_matrix(cD) * u)
+          @test mass_matrix(cD) ≈ Diagonal([weight(cD, i) for i in 1:size(cD, 1)])
           @test cD * u ≈ Matrix(cD) * u atol=eps(float(T)) rtol=sqrt(eps(float(T)))
           @test cD * u ≈ sparse(cD) * u atol=eps(float(T)) rtol=sqrt(eps(float(T)))
           @test SummationByPartsOperators.xmin(cD) ≈ xmin
@@ -363,6 +390,9 @@ end
           cDm_dense = Matrix(cDm)
           M = mass_matrix(cDp)
           @test M ≈ mass_matrix(cDm)
+          for i in 1:size(cD_central, 1)
+            @test @inferred(weight(cDp, i)) ≈ @inferred(weight(cDm, i))
+          end
           @test sum(M) ≈ xmax - xmin
 
           diss = M * (cDp_dense - cDm_dense)
diff --git a/test/dissipation_operators_test.jl b/test/dissipation_operators_test.jl
index 9cfa7c31..d553ab9d 100644
--- a/test/dissipation_operators_test.jl
+++ b/test/dissipation_operators_test.jl
@@ -23,8 +23,11 @@ for source_D in D_test_list, source_Di in Di_test_list, acc_order in 2:2:8, T in
     end
     D === nothing && continue
 
-    @inferred mass_matrix(D)
-    H = mass_matrix(D)
+    M = @inferred mass_matrix(D)
+    for i in 1:size(D, 2)
+        @inferred weight(D, i)
+    end
+    @test M ≈ Diagonal([weight(D, i) for i in 1:size(D, 1)])
 
     for order in 2:2:8
         Di = try
@@ -37,9 +40,9 @@ for source_D in D_test_list, source_Di in Di_test_list, acc_order in 2:2:8, T in
 
         println(devnull, Di)
         println(devnull, Di.coefficients)
-        HDi = H*Matrix(Di)
-        @test norm(HDi - HDi') < 10*eps(T)
-        @test minimum(real, eigvals(HDi)) < 10*eps(T)
+        MDi = M * Matrix(Di)
+        @test norm(MDi - MDi') < 10*eps(T)
+        @test minimum(real, eigvals(MDi)) < 10*eps(T)
     end
 end
 
diff --git a/test/fourier_operators_test.jl b/test/fourier_operators_test.jl
index 59422e13..0f892c2f 100644
--- a/test/fourier_operators_test.jl
+++ b/test/fourier_operators_test.jl
@@ -83,6 +83,7 @@ for T in (Float32, Float64)
 
         @test integrate(u, D) ≈ sum(mass_matrix(D) * u)
         @test integrate(u->u^2, u, D) ≈ dot(u, mass_matrix(D), u)
+        @test mass_matrix(D) ≈ Diagonal([weight(D, i) for i in 1:size(D, 1)])
     end
 end
 
diff --git a/test/legendre_operators_test.jl b/test/legendre_operators_test.jl
index c1b91485..71ec0b86 100644
--- a/test/legendre_operators_test.jl
+++ b/test/legendre_operators_test.jl
@@ -51,6 +51,7 @@ for T in (Float32, Float64), filter_type in (ExponentialFilter(),)
             @test integrate(u->u^2, u, D) <= norm2_u
             @test integrate(u, D) ≈ sum(mass_matrix(D) * u)
             @test integrate(u->u^2, u, D) ≈ dot(u, mass_matrix(D), u)
+            @test mass_matrix(D) ≈ Diagonal([weight(D, i) for i in 1:size(D, 1)])
         end
     end
 end
diff --git a/test/periodic_operators_test.jl b/test/periodic_operators_test.jl
index f961f08b..1ee87b3d 100644
--- a/test/periodic_operators_test.jl
+++ b/test/periodic_operators_test.jl
@@ -208,6 +208,7 @@ let T=Float32
 
     @test integrate(x7, D) ≈ sum(mass_matrix(D) * x7)
     @test integrate(u->u^2, x7, D) ≈ dot(x7, mass_matrix(D), x7)
+    @test mass_matrix(D) ≈ Diagonal([weight(D, i) for i in 1:size(D, 1)])
 end
 
 # Accuracy tests with Float64.
@@ -393,6 +394,7 @@ let T = Float64
 
     @test integrate(x7, D) ≈ sum(mass_matrix(D) * x7)
     @test integrate(u->u^2, x7, D) ≈ dot(x7, mass_matrix(D), x7)
+    @test mass_matrix(D) ≈ Diagonal([weight(D, i) for i in 1:size(D, 1)])
 end
 
 # Compare Fornberg algorithm with exact representation of central derivative coefficients.
@@ -868,6 +870,7 @@ for T in (Float32, Float64)
 
         @test integrate(u, D) ≈ sum(mass_matrix(D) * u)
         @test integrate(u->u^2, u, D) ≈ dot(u, mass_matrix(D), u)
+        @test mass_matrix(D) ≈ Diagonal([weight(D, i) for i in 1:size(D, 1)])
     end
 
     for N in (8, 9), acc_order in (2, 3, 4)
diff --git a/test/upwind_operators_test.jl b/test/upwind_operators_test.jl
index 6a03c703..263a5c2f 100644
--- a/test/upwind_operators_test.jl
+++ b/test/upwind_operators_test.jl
@@ -25,6 +25,9 @@ using SummationByPartsOperators
     M = mass_matrix(Dp_bounded)
     @test M == mass_matrix(Dm_bounded)
     @test M == mass_matrix(Dc_bounded)
+    for i in 1:size(Dp_bounded, 1)
+      @test @inferred(weight(Dp_bounded, i)) ≈ @inferred(weight(Dm_bounded, i)) ≈ @inferred(weight(Dc_bounded, i))
+    end
     Dp_periodic = periodic_derivative_operator(1, acc_order, xmin, xmax, N-1, -(acc_order - 1) ÷ 2)
     Dm_periodic = periodic_derivative_operator(1, acc_order, xmin, xmax, N-1, -acc_order + (acc_order - 1) ÷ 2)
     Dp = Matrix(Dp_bounded)
@@ -83,6 +86,9 @@ end
   @test SummationByPartsOperators.xmax(D) == SummationByPartsOperators.xmax(Dp)
 
   @test mass_matrix(D) == mass_matrix(Dp)
+  for i in 1:size(D, 1)
+    @test @inferred(weight(D, i)) ≈ @inferred(weight(Dp, i))
+  end
   @test left_boundary_weight(D) == left_boundary_weight(Dp)
   @test right_boundary_weight(D) == right_boundary_weight(Dp)