Bodnarchuk, Kaluzhnin, Kotov, Romov’s Theorem on inclusion of Polymorphism ($Pol \rho \subseteq Pol \sigma$) Bodnarchuk, Kaluzhnin, Kotov, Romov’s paper [1]  is well-known. Anne Fearnley [2] infered from it the following theroem and used it to prove the inclusion of polymorphisms.
Theorem (Bodnarchuk, Kaluzhnin, Kotov, Romov). Let $A$ be a finite set. Let
$\rho \subseteq A^{h}$ , and let $\sigma \subseteq A^{l}$ be a relation without repetitions. Then $Pol \rho \subseteq Pol \sigma$ if and only if
there exist $m \geq l, n < m^{h}$ and an $n \times h$ matrix $X = (x_{ij})$ with $x_{ij} \in  \{1, \dots, m\}$
such that $(a_1 ,\dots , a_l) \in \sigma$ iff there exist $a_{l+1} , \dots , a_{m}$ such that for all $i = 1, \dots , n$,
$(a_{x_{i,1}}, a_{x_{i,2}}, \dots , a_{x_{i,h}}) \in \rho$.
My questions are:


*

*I have read [1] several times and am unable to find the proof for necessity part. 
Could you tell me how is the matrix $X$ constructed for the necessity part? Or could you recommend another resource for a complete proof? 

*Anne Fearnley ([2], page 8) used the matrix $X$ =  $
\begin{pmatrix}
  3 & 4 & 1\\
  5 & 3 & 2\\
\end{pmatrix}
$ in the above theorem to prove the following
$Pol\{(0, 0, 0), (1, 1, 1), (0, 1, 2)\} \subset Pol\{(0, 0), (1, 1), (1, 2), (2, 0)\}$,
but how was this matrix $X$ constructed? Is this just by trial? Or is there a general way to construct such a matrix given two relations? 
[1] V. G. Bodnarchuk, L. A. Kaluzhnin, V. N. Kotov, and B. A. Romov. Galois theory for Post algebras I–II, Kibernetika, 3 (1969), pp. 1–10 
and 5 (1969), pp. 1–9 (in Russian); Cybernetics, (1969), pp. 243–252, 531–539 (English version), 1969.
[2] Anne Fearnley, The monoidal interval for the monoid generated by two constants, Journal of Multiple-Valued Logic and Soft Computing, 15(5–6), pp. 597–609, 2009, http://www3.sympatico.ca/anathia/Anne_Fearnley/2-const.pdf‎.
[3] David Geiger, Closed systems of functions and predicates, Pacific J. Math. Volume 27,
Number 1 (1968), 95–100.
 A: $\DeclareMathOperator\Pol{Pol}$In model-theoretic terms, if $A$ is a finite set, and $\rho,\sigma$ relations on $A$, then $\Pol\rho\subseteq\Pol\sigma$ iff $\sigma$ is definable in the structure $(A,\rho)$ by a positive primitive (pp) formula, that is,
$$\tag{$*$}\sigma(x_1,\dots,x_l)\iff\exists y_1,\dots,y_k\,\bigwedge_{i\le n}\theta_i(x_1,\dots,x_l,y_1,\dots,y_k),$$
where each $\theta_i$ is an atomic formula (i.e., $\rho$ or $=$ applied to some of the variables $x_j,y_j$).
Notice that identities of the form $y_i=y_j$ or $y_i=x_j$ can be eliminated from $(*)$ by substituting $y_j$ (resp. $x_j$) for $y_i$ and removing the $y_i$ variable. Thus, we can write $\sigma$ as a conjunction of an $=$-free pp formula and some equalities of the form $x_i=x_j$ for $i\ne j$. However, if $\sigma$ is without repetitions, it cannot imply any equality of the latter form, hence we can in fact express $\sigma$ by an $=$-free pp formula, and this is the conclusion of the theorem you quote.
Now, the nontrivial direction of the characterization above can be proved as follows. Let $\sigma^+$ be the intersection of all (finitely many) pp-definable relations that include $\sigma$. Then $\sigma^+$ is itself pp-definable, hence it suffices to show that $\sigma=\sigma^+$. Assume for contradiction that $\sigma^+(\vec c)$ and $\neg\sigma(\vec c)$ for some $\vec c\in A^l$, and let $I$ be the set of all pp-definable relations $\tau\subseteq A^l$ such that $\neg\tau(\vec c)$. For every $\tau\in I$, we have $\sigma^+\nsubseteq\tau$, hence $\sigma\nsubseteq\tau$ by the definition of $\sigma^+$. Thus, we can fix $a_1^\tau,\dots,a_l^\tau\in A$ such that $\sigma(\vec a^\tau)$ and $\neg\tau(\vec a^\tau)$.
For each $i=1,\dots,l$, let $a_i$ be the element of the cartesian product $A^I$ whose $\tau$-coordinate is $a^\tau_i$ for every $\tau\in I$. By construction, we have
$$\tag{$**$}A^I\models\phi(a_1,\dots,a_l)\implies A\models\phi(c_1,\dots,c_l)$$
for every pp formula $\phi$. Enumerate $A^I=\{u_j:j=1,\dots,r\}$, and put
$$\phi(x_1,\dots,x_l)=\exists y_1,\dots,y_r\bigwedge\{\theta(x_1,\dots,x_l,y_1,\dots,y_r):A^I\models\theta(a_1,\dots,a_l,u_1,\dots,u_r)\},$$
where $\theta$ runs over all atomic formulas. Since $A^I\models\phi(\vec a)$ and $\phi$ is pp, we have $A\models\phi(\vec c)$ by $(**)$, hence we can fix $v_1,\dots,v_r\in A$ that witness the existential quantifiers in $\phi$. Then $f(u_i):=v_i$ defines a homomorphism $f\colon(A,\rho)^I\to(A,\rho)$ such that $f(a_i)=c_i$. In clone-theoretic terminology, this means that $f$ is a polymorphism of $\rho$, but since $\sigma(\vec a^\tau)$ and $\neg\sigma(\vec c)$, it is not a polymorphism of $\sigma$, which contradicts the assumption.
