# factorization into `$XX^{T}$` may not be unique!
raise ValueError('Unable to determine extremity of ``A``.')
+
+
+def completely_positive_operators_gens(K):
+ r"""
+ Return a list of generators (matrices) for the completely-positive
+ cone of ``K``.
+
+ INPUT:
+
+ - ``K`` -- a closed convex rational polyhedral cone.
+
+ OUTPUT:
+
+ A list of matrices, the conic hull of which is the
+ completely-positive cone of ``K``.
+
+ SETUP::
+
+ sage: from mjo.cone.completely_positive import (
+ ....: completely_positive_operators_gens,
+ ....: is_completely_positive )
+ sage: from mjo.cone.nonnegative_orthant import nonnegative_orthant
+ sage: from mjo.matrix_vector import isomorphism
+
+ EXAMPLES::
+
+ sage: K = nonnegative_orthant(2)
+ sage: completely_positive_operators_gens(K)
+ [
+ [1 0] [0 0]
+ [0 0], [0 1]
+ ]
+ sage: all( is_completely_positive(M)
+ ....: for M in completely_positive_operators_gens(K) )
+ True
+
+ TESTS:
+
+ The completely-positive cone of ``K`` is subdual::
+
+ sage: K = random_cone(max_ambient_dim=8, max_rays=10)
+ sage: cp_gens = completely_positive_operators_gens(K)
+ sage: n = K.lattice_dim()
+ sage: M = MatrixSpace(QQ, n, n)
+ sage: (p, p_inv) = isomorphism(M)
+ sage: L = ToricLattice(n**2)
+ sage: cp_cone = Cone( (p(m) for m in cp_gens), lattice=L )
+ sage: copos_cone = Cone(cp_cone.dual().rays(), lattice=L )
+ sage: all( x in copos_cone for x in cp_cone )
+ True
+
+ """
+ return [ x.tensor_product(x) for x in K ]
+