add implementation of QuadraticForm

src/mstructures.jl

@@ -61,14 +61,22 @@ function MA.operate!(::UnsafeAddMul, res, c)
     return res
 end
 
-function MA.operate!(op::UnsafeAddMul, res, b, c, args::Vararg{Any, N}) where {N}
+function MA.operate!(
+    op::UnsafeAddMul,
+    res,
+    b,
+    c,
+    args::Vararg{Any,N};
+    cfs = nothing,
+) where {N}
     for (kb, vb) in nonzero_pairs(b)
         for (kc, vc) in nonzero_pairs(c)
             for (k, v) in nonzero_pairs(op.structure(kb, kc))
+                _cfs = isnothing(cfs) ? vb * vc * v : vb * vc * v * cfs
                 MA.operate!(
                     op,
                     res,
-                    SparseCoefficients((_key(op.structure, k),), (vb * vc * v,)),
+                    SparseCoefficients((_key(op.structure, k),), (_cfs,)),
                     args...,
                 )
             end
@@ -85,3 +93,34 @@ function (mstr::DiracMStructure)(x::T, y::T) where {T}
     xy = mstr.op(x, y)
     return SparseCoefficients((xy,), (1,))
 end
+
+"""
+    QuadraticForm(Q)
+A simple wrapper for representing a quadratic form.
+
+A `QuadraticForm(Q)` represents actually a **sesquilinear form** with respect to
+its basis `b`, i.e. `star.(b)·Q·b = ∑ᵢⱼ Q[i,j]·star(b[i])·b[j]`.
+
+# Necessary methods:
+ * `basis(Q)` - a **finite** basis `b` of the quadratic form;
+ * `Base.getindex(Q, i::T, j::T)` - the value at `i`-th and `j`-th basis elements
+   where `T` is the type of indicies of `basis(Q)`;
+ * `Base.eltype(Q)` - the type of `Q[i,j]`.
+"""
+struct QuadraticForm{T}
+    Q::T
+end
+
+Base.eltype(qf::QuadraticForm) = eltype(qf.Q)
+basis(qf::QuadraticForm) = basis(qf.Q)
+Base.getindex(qf::QuadraticForm, i::T, j::T) where {T} = qf.Q[i, j]
+
+function MA.operate!(op::UnsafeAddMul, res, Q::QuadraticForm)
+    for (i, b1) in pairs(basis(Q))
+        b1★ = star(b1)
+        for (j, b2) in pairs(basis(Q))
+            MA.operate!(op, res, coeffs(b1★), coeffs(b2); cfs = Q[i, j])
+        end
+    end
+    return res
+end

src/types.jl

@@ -131,3 +131,12 @@ function AlgebraElement{T}(X::AlgebraElement) where {T}
     end
     return AlgebraElement(w, parent(X))
 end
+
+function AlgebraElement(qf::QuadraticForm, A::AbstractStarAlgebra)
+    @assert all(b -> parent(b) == A, basis(qf))
+    res = zero(eltype(qf), A)
+    MA.operate_to!(coeffs(res), mstructure(A), qf)
+    return res
+end
+
+(A::AbstractStarAlgebra)(qf::QuadraticForm) = AlgebraElement(qf, A)

test/quadratic_form.jl

@@ -0,0 +1,67 @@
+struct Gramm{T,B}
+    matrix::T
+    basis::B
+end
+
+function Base.eltype(g::Gramm)
+    return promote_type(eltype(g.matrix), eltype(eltype(basis(g))))
+end
+SA.basis(g::Gramm) = g.basis
+Base.getindex(g::Gramm, i, j) = g.matrix[i, j]
+
+@testset "QuadraticForm" begin
+    A = let alph = [:a, :b, :c]
+        fw = FreeWords(alph)
+        SA.StarAlgebra(fw, SA.DiracBasis(fw))
+    end
+
+    gbasis = let (id, a, b, c) = A.(Iterators.take(SA.object(A), 4))
+        # basis has to be star-invariant:
+        bas = 1.0 * [one(A), (a + b) / 2, (a + c) / 2, (b + c) / 2]
+        SA.FixedBasis(bas, SA.DiracMStructure(*))
+    end
+
+    m = [
+        π 0 0 0
+        0 0 0 0
+        0 0 0 0
+        0 0 0 0
+    ]
+    Q = SA.QuadraticForm(Gramm(m, gbasis))
+    b = basis(Q)
+    @test A(Q) == π * b[1]' * b[1]
+
+    m = [
+        0 1//2 0 0
+        0 0 0 0
+        0 0 0 0
+        0 0 0 0
+    ]
+    Q = SA.QuadraticForm(Gramm(m, gbasis))
+    b = basis(Q)
+    @test A(Q) == 1 // 2 * b[1]' * b[2]
+
+    m = [
+        0 0 0 0
+        0 π -1 0
+        0 0 0 0
+        0 0 0 0
+    ]
+    Q = SA.QuadraticForm(Gramm(m, gbasis))
+    b = basis(Q)
+    @test A(Q) == π * b[2]' * b[2] - b[2]' * b[3]
+
+    m = [
+        0 0 0 0
+        0 0 0 0
+        0 π 0 0
+        0 0 1 0
+    ]
+    Q = SA.QuadraticForm(Gramm(m, gbasis))
+    b = basis(Q)
+    @test A(Q) == π * b[3]' * b[2] + b[4]' * b[3]
+
+    m = ones(Int, 4, 4)
+    Q = SA.QuadraticForm(Gramm(m, gbasis))
+    @test A(Q) == sum(bi' * bj for bi in gbasis for bj in gbasis)
+end

test/runtests.jl

@@ -37,6 +37,7 @@ include("test_example_acoeffs.jl")
 
     # free monoid algebra
     include("monoid_algebra.jl")
+    include("quadratic_form.jl")
 
     include("caching_allocations.jl")