From 9b8d1bf1d192a0773531bfa00543d88ff7caf9d6 Mon Sep 17 00:00:00 2001
From: David Hong <hong@udel.edu>
Date: Thu, 22 Aug 2024 17:19:46 -0400
Subject: [PATCH 1/6] Add sortcomps/sortcomps!

---
 src/GCPDecompositions.jl |  4 +++-
 src/cpd.jl               | 19 +++++++++++++++++++
 2 files changed, 22 insertions(+), 1 deletion(-)

diff --git a/src/GCPDecompositions.jl b/src/GCPDecompositions.jl
index e04d575..17ad342 100644
--- a/src/GCPDecompositions.jl
+++ b/src/GCPDecompositions.jl
@@ -8,12 +8,14 @@ import Base: ndims, size, show, summary
 import Base: getindex
 import Base: AbstractArray, Array
 import LinearAlgebra: norm
+using Base.Order: Ordering, Reverse
 using IntervalSets: Interval
 using Random: default_rng
 
 # Exports
 export CPD
-export ncomps, normalizecomps, normalizecomps!, permutecomps, permutecomps!
+export ncomps,
+    normalizecomps, normalizecomps!, permutecomps, permutecomps!, sortcomps, sortcomps!
 export gcp
 export GCPLosses, GCPConstraints, GCPAlgorithms
 
diff --git a/src/cpd.jl b/src/cpd.jl
index bbb1c0f..1d2c672 100644
--- a/src/cpd.jl
+++ b/src/cpd.jl
@@ -236,3 +236,22 @@ function permutecomps!(M::CPD, perm::Vector)
     # Return CPD with permuted components
     return M
 end
+
+sortcomps(M::CPD; dims = :λ, order::Ordering = Reverse, kwargs...) =
+    permutecomps(M, sortperm(_sortvals(M, dims); order, kwargs...))
+sortcomps!(M::CPD; dims = :λ, order::Ordering = Reverse, kwargs...) =
+    permutecomps!(M, sortperm(_sortvals(M, dims); order, kwargs...))
+
+function _sortvals(M::CPD, dims)
+    # Check dims
+    dims_iterable = dims isa Symbol ? (dims,) : dims
+    all(d -> d === :λ || (d isa Integer && d in 1:ndims(M)), dims_iterable) || throw(
+        ArgumentError(
+            "`dims` must be `:λ`, an integer specifying a mode, or a collection, got $dims",
+        ),
+    )
+
+    # Return vector of values to sort by
+    return dims === :λ ? M.λ :
+           [map(d -> d === :λ ? M.λ[j] : view(M.U[d], :, j), dims) for j in 1:ncomps(M)]
+end

From e28b3b24e884cc17c6f1f717b67df06bfdd59c6f Mon Sep 17 00:00:00 2001
From: David Hong <hong@udel.edu>
Date: Thu, 22 Aug 2024 19:00:41 -0400
Subject: [PATCH 2/6] Add docstrings

---
 docs/src/man/main.md |  2 ++
 src/cpd.jl           | 25 +++++++++++++++++++++++--
 2 files changed, 25 insertions(+), 2 deletions(-)

diff --git a/docs/src/man/main.md b/docs/src/man/main.md
index 08be824..d99ee69 100644
--- a/docs/src/man/main.md
+++ b/docs/src/man/main.md
@@ -13,6 +13,8 @@ normalizecomps
 normalizecomps!
 permutecomps
 permutecomps!
+sortcomps
+sortcomps!
 GCPDecompositions.default_constraints
 GCPDecompositions.default_algorithm
 GCPDecompositions.default_init
diff --git a/src/cpd.jl b/src/cpd.jl
index 1d2c672..e2b2672 100644
--- a/src/cpd.jl
+++ b/src/cpd.jl
@@ -209,7 +209,7 @@ end
 Permute the components of `M`.
 `perm` is a vector or a tuple of length `ncomps(M)` specifying the permutation.
 
-See also: `permutecomps!`.
+See also: `permutecomps!`, `sortcomps`, `sortcomps!`.
 """
 permutecomps(M::CPD, perm) = permutecomps!(deepcopy(M), perm)
 
@@ -219,7 +219,7 @@ permutecomps(M::CPD, perm) = permutecomps!(deepcopy(M), perm)
 Permute the components of `M` in-place.
 `perm` is a vector or a tuple of length `ncomps(M)` specifying the permutation.
 
-See also: `permutecomps`.
+See also: `permutecomps`, `sortcomps`, `sortcomps!`.
 """
 permutecomps!(M::CPD, perm) = permutecomps!(M, collect(perm))
 function permutecomps!(M::CPD, perm::Vector)
@@ -237,8 +237,29 @@ function permutecomps!(M::CPD, perm::Vector)
     return M
 end
 
+"""
+    sortcomps(M::CPD; dims=:λ, alg::Algorithm=DEFAULT_UNSTABLE, lt=isless, by=identity, rev::Bool=false, order::Ordering=Reverse)
+
+Sort the components of `M`. `dims` specifies what part to sort by;
+it must be the symbol `:λ`, an integer in `1:ndims(M)`, or a collection of these.
+
+For the remaining keyword arguments, see the documentation of `sort!`.
+
+See also: `sortcomps!`, `sort`, `sort!`.
+"""
 sortcomps(M::CPD; dims = :λ, order::Ordering = Reverse, kwargs...) =
     permutecomps(M, sortperm(_sortvals(M, dims); order, kwargs...))
+
+"""
+    sortcomps!(M::CPD; dims=:λ, alg::Algorithm=DEFAULT_UNSTABLE, lt=isless, by=identity, rev::Bool=false, order::Ordering=Reverse)
+
+Sort the components of `M` in-place. `dims` specifies what part to sort by;
+it must be the symbol `:λ`, an integer in `1:ndims(M)`, or a collection of these.
+
+For the remaining keyword arguments, see the documentation of `sort!`.
+
+See also: `sortcomps`, `sort`, `sort!`.
+"""
 sortcomps!(M::CPD; dims = :λ, order::Ordering = Reverse, kwargs...) =
     permutecomps!(M, sortperm(_sortvals(M, dims); order, kwargs...))
 

From 6cebeac11c4691dc4695bf38f76f83a38ca21ee3 Mon Sep 17 00:00:00 2001
From: David Hong <hong@udel.edu>
Date: Wed, 28 Aug 2024 08:48:23 -0400
Subject: [PATCH 3/6] Add CPDComp type

---
 src/cpdcomp.jl | 166 +++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 166 insertions(+)
 create mode 100644 src/cpdcomp.jl

diff --git a/src/cpdcomp.jl b/src/cpdcomp.jl
new file mode 100644
index 0000000..2ac5031
--- /dev/null
+++ b/src/cpdcomp.jl
@@ -0,0 +1,166 @@
+## CPD component type
+
+"""
+    CPDComp
+
+Type for a single component of a canonical polyadic decompositions (CPD).
+
+If `M::CPDComp` is the component object,
+the scalar weight `λ` and the factor vectors `u = (u[1],...,u[N])`
+can be obtained via `M.λ` and `M.u`.
+"""
+struct CPDComp{T,N,Tu<:AbstractVector{T}}
+    λ::T
+    u::NTuple{N,Tu}
+    function CPDComp{T,N,Tu}(λ, u) where {T,N,Tu<:AbstractVector{T}}
+        Base.require_one_based_indexing(u...)
+        return new{T,N,Tu}(λ, u)
+    end
+end
+CPDComp(λ::T, u::NTuple{N,Tu}) where {T,N,Tu<:AbstractVector{T}} = CPDComp{T,N,Tu}(λ, u)
+
+ndims(::CPDComp{T,N}) where {T,N} = N
+size(M::CPDComp{T,N}, dim::Integer) where {T,N} = dim <= N ? length(M.u[dim]) : 1
+size(M::CPDComp{T,N}) where {T,N} = ntuple(d -> size(M, d), N)
+
+function show(io::IO, mime::MIME{Symbol("text/plain")}, M::CPDComp{T,N}) where {T,N}
+    # Compute displaysize for showing fields
+    LINES, COLUMNS = displaysize(io)
+    LINES_FIELD = max(LINES - 2 - N, 0) ÷ (1 + N)
+    io_field = IOContext(io, :displaysize => (LINES_FIELD, COLUMNS))
+
+    # Show summary and fields
+    summary(io, M)
+    println(io)
+    println(io, "λ weight:")
+    show(io_field, mime, M.λ)
+    for k in Base.OneTo(N)
+        println(io, "\nu[$k] factor vector:")
+        show(io_field, mime, M.u[k])
+    end
+end
+
+function summary(io::IO, M::CPDComp)
+    dimstring =
+        ndims(M) == 0 ? "0-dimensional" :
+        ndims(M) == 1 ? "$(size(M,1))-element" : join(map(string, size(M)), '×')
+    return print(io, dimstring, " ", typeof(M))
+end
+
+function getindex(M::CPDComp{T,N}, I::Vararg{Int,N}) where {T,N}
+    @boundscheck Base.checkbounds_indices(Bool, axes(M), I) || Base.throw_boundserror(M, I)
+    return M.λ * prod(M.u[k][I[k]] for k in Base.OneTo(ndims(M)))
+end
+getindex(M::CPDComp{T,N}, I::CartesianIndex{N}) where {T,N} = getindex(M, Tuple(I)...)
+
+AbstractArray(A::CPDComp) =
+    reshape(TensorKernels.khatrirao(reverse(reshape.(A.u, :, 1))...) * A.λ, size(A))
+Array(A::CPDComp) = Array(AbstractArray(A))
+
+norm(M::CPDComp, p::Real = 2) =
+    p == 2 ? norm2(M) : norm((M[I] for I in CartesianIndices(size(M))), p)
+norm2(M::CPDComp{T,N}) where {T,N} = sqrt(abs2(M.λ) * prod(sum(abs2, M.u[i]) for i in 1:N))
+
+"""
+    normalizecomps(M::CPDComp, p::Real = 2)
+
+Normalize `M` so that all its factor vectors have `p`-norm equal to unity,
+i.e., `norm(M.u[k], p) == 1` for all `k ∈ 1:ndims(M)`. The excess weight is absorbed into `M.λ`.
+Norms equal to zero are ignored (i.e., treated as though they were equal to one).
+
+The following keyword arguments can be used to modify this behavior:
+- `dims` specifies what to normalize (default: `[:λ; 1:ndims(M)]`)
+- `distribute_to` specifies where to distribute the excess weight (default: `:λ`)
+Valid options for these arguments are the symbol `:λ`, an integer in `1:ndims(M)`,
+or a collection of these.
+
+See also: `normalizecomps!`, `norm`.
+"""
+normalizecomps(M::CPDComp, p::Real = 2; dims = [:λ; 1:ndims(M)], distribute_to = :λ) =
+    normalizecomps!(deepcopy(M), p; dims, distribute_to)
+
+"""
+    normalizecomps!(M::CPDComp, p::Real = 2)
+
+Normalize `M` in-place so that all its factor vectors have `p`-norm equal to unity,
+i.e., `norm(M.u[k], p) == 1` for all `k ∈ 1:ndims(M)`. The excess weight is absorbed into `M.λ`.
+Norms equal to zero are ignored (i.e., treated as though they were equal to one).
+
+The following keyword arguments can be used to modify this behavior:
+- `dims` specifies what to normalize (default: `[:λ; 1:ndims(M)]`)
+- `distribute_to` specifies where to distribute the excess weight (default: `:λ`)
+Valid options for these arguments are the symbol `:λ`, an integer in `1:ndims(M)`,
+or a collection of these.
+
+See also: `normalizecomps`, `norm`.
+"""
+function normalizecomps!(
+    M::CPDComp{T,N},
+    p::Real = 2;
+    dims = [:λ; 1:N],
+    distribute_to = :λ,
+) where {T,N}
+    # Check dims and put into standard (mask) form
+    dims_iterable = dims isa Symbol ? (dims,) : dims
+    all(d -> d === :λ || (d isa Integer && d in 1:N), dims_iterable) || throw(
+        ArgumentError(
+            "`dims` must be `:λ`, an integer specifying a mode, or a collection, got $dims",
+        ),
+    )
+    dims_λ = :λ in dims_iterable
+    dims_u = ntuple(in(dims_iterable), N)
+
+    # Check distribute_to and put into standard (mask) form
+    dist_iterable = distribute_to isa Symbol ? (distribute_to,) : distribute_to
+    all(d -> d === :λ || (d isa Integer && d in 1:N), dist_iterable) || throw(
+        ArgumentError(
+            "`distribute_to` must be `:λ`, an integer specifying a mode, or a collection, got $distribute_to",
+        ),
+    )
+    dist_λ = :λ in dist_iterable
+    dist_u = ntuple(in(dist_iterable), N)
+
+    # Call inner function
+    return _normalizecomps!(M, p, dims_λ, dims_u, dist_λ, dist_u)
+end
+
+function _normalizecomps!(
+    M::CPDComp{T,N},
+    p::Real,
+    dims_λ::Bool,
+    dims_u::NTuple{N,Bool},
+    dist_λ::Bool,
+    dist_u::NTuple{N,Bool},
+) where {T,N}
+    # Utility function to handle zero weights and norms
+    zero_to_one(x) = iszero(x) ? oneunit(x) : x
+
+    # Normalize components and collect excess weight
+    excess = oneunit(T)
+    if dims_λ
+        _norm = zero_to_one(abs(M.λ))
+        M.λ /= _norm
+        excess *= _norm
+    end
+    for k in Base.OneTo(N)
+        if dims_u[k]
+            _norm = zero_to_one(norm(M.u[k], p))
+            M.u[k] ./= _norm
+            excess *= _norm
+        end
+    end
+
+    # Distribute excess weight (uniformly across specified parts)
+    excess = excess^(1 / count((dist_λ, dist_u...)))
+    if dist_λ
+        M.λ *= excess
+    end
+    for k in Base.OneTo(N)
+        if dist_u[k]
+            M.u[k] .*= excess
+        end
+    end
+
+    # Return normalized CPDComp
+    return M
+end

From 034d6c1201f7490fd226012b5c8775e7f0f46b7e Mon Sep 17 00:00:00 2001
From: David Hong <hong@udel.edu>
Date: Thu, 29 Aug 2024 18:04:19 -0400
Subject: [PATCH 4/6] Include and export CPDComp

---
 src/GCPDecompositions.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/GCPDecompositions.jl b/src/GCPDecompositions.jl
index 17ad342..ec0cb14 100644
--- a/src/GCPDecompositions.jl
+++ b/src/GCPDecompositions.jl
@@ -13,7 +13,7 @@ using IntervalSets: Interval
 using Random: default_rng
 
 # Exports
-export CPD
+export CPD, CPDComp
 export ncomps,
     normalizecomps, normalizecomps!, permutecomps, permutecomps!, sortcomps, sortcomps!
 export gcp
@@ -21,6 +21,7 @@ export GCPLosses, GCPConstraints, GCPAlgorithms
 
 include("tensor-kernels.jl")
 include("cpd.jl")
+include("cpdcomp.jl")
 include("gcp-losses.jl")
 include("gcp-constraints.jl")
 include("gcp-algorithms.jl")

From 34a3f26a0a52c9220e8f4fa869e69bf5c7e16ee9 Mon Sep 17 00:00:00 2001
From: David Hong <hong@udel.edu>
Date: Thu, 29 Aug 2024 18:04:37 -0400
Subject: [PATCH 5/6] Drop normalizecomps since CPDComp is immutable

---
 src/cpdcomp.jl | 104 -------------------------------------------------
 1 file changed, 104 deletions(-)

diff --git a/src/cpdcomp.jl b/src/cpdcomp.jl
index 2ac5031..0449e3c 100644
--- a/src/cpdcomp.jl
+++ b/src/cpdcomp.jl
@@ -60,107 +60,3 @@ Array(A::CPDComp) = Array(AbstractArray(A))
 norm(M::CPDComp, p::Real = 2) =
     p == 2 ? norm2(M) : norm((M[I] for I in CartesianIndices(size(M))), p)
 norm2(M::CPDComp{T,N}) where {T,N} = sqrt(abs2(M.λ) * prod(sum(abs2, M.u[i]) for i in 1:N))
-
-"""
-    normalizecomps(M::CPDComp, p::Real = 2)
-
-Normalize `M` so that all its factor vectors have `p`-norm equal to unity,
-i.e., `norm(M.u[k], p) == 1` for all `k ∈ 1:ndims(M)`. The excess weight is absorbed into `M.λ`.
-Norms equal to zero are ignored (i.e., treated as though they were equal to one).
-
-The following keyword arguments can be used to modify this behavior:
-- `dims` specifies what to normalize (default: `[:λ; 1:ndims(M)]`)
-- `distribute_to` specifies where to distribute the excess weight (default: `:λ`)
-Valid options for these arguments are the symbol `:λ`, an integer in `1:ndims(M)`,
-or a collection of these.
-
-See also: `normalizecomps!`, `norm`.
-"""
-normalizecomps(M::CPDComp, p::Real = 2; dims = [:λ; 1:ndims(M)], distribute_to = :λ) =
-    normalizecomps!(deepcopy(M), p; dims, distribute_to)
-
-"""
-    normalizecomps!(M::CPDComp, p::Real = 2)
-
-Normalize `M` in-place so that all its factor vectors have `p`-norm equal to unity,
-i.e., `norm(M.u[k], p) == 1` for all `k ∈ 1:ndims(M)`. The excess weight is absorbed into `M.λ`.
-Norms equal to zero are ignored (i.e., treated as though they were equal to one).
-
-The following keyword arguments can be used to modify this behavior:
-- `dims` specifies what to normalize (default: `[:λ; 1:ndims(M)]`)
-- `distribute_to` specifies where to distribute the excess weight (default: `:λ`)
-Valid options for these arguments are the symbol `:λ`, an integer in `1:ndims(M)`,
-or a collection of these.
-
-See also: `normalizecomps`, `norm`.
-"""
-function normalizecomps!(
-    M::CPDComp{T,N},
-    p::Real = 2;
-    dims = [:λ; 1:N],
-    distribute_to = :λ,
-) where {T,N}
-    # Check dims and put into standard (mask) form
-    dims_iterable = dims isa Symbol ? (dims,) : dims
-    all(d -> d === :λ || (d isa Integer && d in 1:N), dims_iterable) || throw(
-        ArgumentError(
-            "`dims` must be `:λ`, an integer specifying a mode, or a collection, got $dims",
-        ),
-    )
-    dims_λ = :λ in dims_iterable
-    dims_u = ntuple(in(dims_iterable), N)
-
-    # Check distribute_to and put into standard (mask) form
-    dist_iterable = distribute_to isa Symbol ? (distribute_to,) : distribute_to
-    all(d -> d === :λ || (d isa Integer && d in 1:N), dist_iterable) || throw(
-        ArgumentError(
-            "`distribute_to` must be `:λ`, an integer specifying a mode, or a collection, got $distribute_to",
-        ),
-    )
-    dist_λ = :λ in dist_iterable
-    dist_u = ntuple(in(dist_iterable), N)
-
-    # Call inner function
-    return _normalizecomps!(M, p, dims_λ, dims_u, dist_λ, dist_u)
-end
-
-function _normalizecomps!(
-    M::CPDComp{T,N},
-    p::Real,
-    dims_λ::Bool,
-    dims_u::NTuple{N,Bool},
-    dist_λ::Bool,
-    dist_u::NTuple{N,Bool},
-) where {T,N}
-    # Utility function to handle zero weights and norms
-    zero_to_one(x) = iszero(x) ? oneunit(x) : x
-
-    # Normalize components and collect excess weight
-    excess = oneunit(T)
-    if dims_λ
-        _norm = zero_to_one(abs(M.λ))
-        M.λ /= _norm
-        excess *= _norm
-    end
-    for k in Base.OneTo(N)
-        if dims_u[k]
-            _norm = zero_to_one(norm(M.u[k], p))
-            M.u[k] ./= _norm
-            excess *= _norm
-        end
-    end
-
-    # Distribute excess weight (uniformly across specified parts)
-    excess = excess^(1 / count((dist_λ, dist_u...)))
-    if dist_λ
-        M.λ *= excess
-    end
-    for k in Base.OneTo(N)
-        if dist_u[k]
-            M.u[k] .*= excess
-        end
-    end
-
-    # Return normalized CPDComp
-    return M
-end

From 0655cd17558bc2279fbe60950fa38ba7c4e25caa Mon Sep 17 00:00:00 2001
From: David Hong <hong@udel.edu>
Date: Thu, 29 Aug 2024 18:04:51 -0400
Subject: [PATCH 6/6] Add tests for CPDComp

---
 test/items/cpdcomp.jl | 127 ++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 127 insertions(+)
 create mode 100644 test/items/cpdcomp.jl

diff --git a/test/items/cpdcomp.jl b/test/items/cpdcomp.jl
new file mode 100644
index 0000000..ed83792
--- /dev/null
+++ b/test/items/cpdcomp.jl
@@ -0,0 +1,127 @@
+## CPD component type
+
+@testitem "constructors" begin
+    using OffsetArrays
+
+    @testset "T=$T" for T in [Float64, Float16]
+        λ = T(100)
+        u1, u2, u3 = T[1, 4], T[-1], T[2, 5, 8]
+
+        # Check type for various orders
+        @test CPDComp{T,0,Vector{T}}(λ, ()) isa CPDComp{T,0,Vector{T}}
+        @test CPDComp(λ, (u1,)) isa CPDComp{T,1,Vector{T}}
+        @test CPDComp(λ, (u1, u2)) isa CPDComp{T,2,Vector{T}}
+        @test CPDComp(λ, (u1, u2, u3)) isa CPDComp{T,3,Vector{T}}
+
+        # Check requirement of one-based indexing
+        O1, O2 = OffsetArray(u1, 0:1), OffsetArray(u2, 0:0)
+        @test_throws ArgumentError CPDComp(λ, (O1, O2))
+    end
+end
+
+@testitem "ndims" begin
+    λ = 100
+    u1, u2, u3 = [1, 4], [-1], [2, 5, 8]
+
+    @test ndims(CPDComp{Int,0,Vector{Int}}(λ, ())) == 0
+    @test ndims(CPDComp(λ, (u1,))) == 1
+    @test ndims(CPDComp(λ, (u1, u2))) == 2
+    @test ndims(CPDComp(λ, (u1, u2, u3))) == 3
+end
+
+@testitem "size" begin
+    λ = 100
+    u1, u2, u3 = [1, 4], [-1], [2, 5, 8]
+
+    @test size(CPDComp(λ, (u1,))) == (length(u1),)
+    @test size(CPDComp(λ, (u1, u2))) == (length(u1), length(u2))
+    @test size(CPDComp(λ, (u1, u2, u3))) == (length(u1), length(u2), length(u3))
+
+    M = CPDComp(λ, (u1, u2, u3))
+    @test size(M, 1) == 2
+    @test size(M, 2) == 1
+    @test size(M, 3) == 3
+    @test size(M, 4) == 1
+end
+
+@testitem "show / summary" begin
+    M = CPDComp(rand(), rand.((3, 4, 5)))
+    Mstring = sprint((t, s) -> show(t, "text/plain", s), M)
+    λstring = sprint((t, s) -> show(t, "text/plain", s), M.λ)
+    ustrings = sprint.((t, s) -> show(t, "text/plain", s), M.u)
+    @test Mstring == string(
+        "$(summary(M))\nλ weight:\n$λstring",
+        ["\nu[$k] factor vector:\n$ustring" for (k, ustring) in enumerate(ustrings)]...,
+    )
+end
+
+@testitem "getindex" begin
+    T = Float64
+    λ = T(100)
+    u1, u2, u3 = T[1, 4], T[-1], T[2, 5, 8]
+
+    M = CPDComp(λ, (u1, u2, u3))
+    for i1 in axes(u1, 1), i2 in axes(u2, 1), i3 in axes(u3, 1)
+        Mi = λ * u1[i1] * u2[i2] * u3[i3]
+        @test Mi == M[i1, i2, i3]
+        @test Mi == M[CartesianIndex((i1, i2, i3))]
+    end
+    @test_throws BoundsError M[length(u1)+1, 1, 1]
+    @test_throws BoundsError M[1, length(u2)+1, 1]
+    @test_throws BoundsError M[1, 1, length(u3)+1]
+
+    M = CPDComp(λ, (u1, u2))
+    for i1 in axes(u1, 1), i2 in axes(u2, 1)
+        Mi = λ * u1[i1] * u2[i2]
+        @test Mi == M[i1, i2]
+        @test Mi == M[CartesianIndex((i1, i2))]
+    end
+    @test_throws BoundsError M[length(u1)+1, 1]
+    @test_throws BoundsError M[1, length(u2)+1]
+
+    M = CPDComp(λ, (u1,))
+    for i1 in axes(u1, 1)
+        Mi = λ * u1[i1]
+        @test Mi == M[i1]
+        @test Mi == M[CartesianIndex((i1,))]
+    end
+    @test_throws BoundsError M[length(u1)+1]
+end
+
+@testitem "Array" begin
+    @testset "N=$N" for N in 1:3
+        T = Float64
+        λ = T(100)
+        u1, u2, u3 = T[1, 4], T[-1], T[2, 5, 8]
+        M = CPDComp(λ, (u1, u2, u3))
+
+        X = Array(M)
+        @test all(I -> M[I] == X[I], CartesianIndices(X))
+    end
+end
+
+@testitem "norm" begin
+    using LinearAlgebra
+
+    T = Float64
+    λ = T(100)
+    u1, u2, u3 = T[1, 4], T[-1], T[2, 5, 8]
+
+    M = CPDComp(λ, (u1, u2, u3))
+    @test norm(M) == norm(M, 2) == sqrt(sum(abs2, M[I] for I in CartesianIndices(size(M))))
+    @test norm(M, 1) == sum(abs, M[I] for I in CartesianIndices(size(M)))
+    @test norm(M, 3) ==
+          (sum(m -> abs(m)^3, M[I] for I in CartesianIndices(size(M))))^(1 / 3)
+
+    M = CPDComp(λ, (u1, u2))
+    @test norm(M) == norm(M, 2) == sqrt(sum(abs2, M[I] for I in CartesianIndices(size(M))))
+    @test norm(M, 1) == sum(abs, M[I] for I in CartesianIndices(size(M)))
+    @test norm(M, 3) ==
+          (sum(m -> abs(m)^3, M[I] for I in CartesianIndices(size(M))))^(1 / 3)
+
+    M = CPDComp(λ, (u1,))
+    @test norm(M) == norm(M, 2) == sqrt(sum(abs2, M[I] for I in CartesianIndices(size(M))))
+    @test norm(M, 1) == sum(abs, M[I] for I in CartesianIndices(size(M)))
+    @test norm(M, 3) ==
+          (sum(m -> abs(m)^3, M[I] for I in CartesianIndices(size(M))))^(1 / 3)
+end