in Standard/src/Synthesis/UnitaryDecomposition.cs [43:104]
private static List<TwoLevelUnitary> TwoLevelDecompose(Complex[,] A)
{
int n = A.GetLength(0);
var result = new List<TwoLevelUnitary>();
for (int i = 0; i < n - 2; i++)
{
for (int j = n - 1; j > i; j--)
{
// At this step we will multiply A (from the right) by such two-level
// unitary that A[i, j] becomes zero.
var a = A[i, j - 1];
var b = A[i, j];
Complex[,] mx;
if (A[i, j].Magnitude <= tol)
{
// Element is already zero. We shouldn't do anything, which equivalent
// to multiplying by identity matrix.
mx = new Complex[,] { { 1.0, 0.0 }, { 0.0, 1.0 } };
// But if it's last in a row, ensure that diagonal element will be 1.
if (j == i + 1)
{
mx = new Complex[,] { { 1 / a, 0.0 }, { 0.0, a } };
}
}
else if (A[i, j - 1].Magnitude <= tol)
{
// Just swap columns with Pauli X matrix.
mx = new Complex[,] { { 0.0, 1.0 }, { 1.0, 0.0 } };
// But if it's last in a row, ensure that diagonal element will be 1.
if (j == i + 1)
{
mx = new Complex[,] { { 0.0, b }, { 1 / b, 0.0 } };
}
}
else
{
mx = MakeEliminatingMatrix(A[i, j - 1], A[i, j]);
}
var u2x2 = new TwoLevelUnitary(mx, j - 1, j);
u2x2.MultiplyRight(A);
u2x2.ConjugateTranspose();
if (!u2x2.IsIdentity())
{
result.Add(u2x2);
}
}
// At this point all elements in row i and in column i are zeros, with exception
// of diagonal element A[i, i], which is equal to 1.
Debug.Assert((A[i, i] - 1.0).Magnitude < tol);
}
var lastMatrix = new TwoLevelUnitary(new Complex[,] {
{ A[n - 2, n - 2], A[n - 2, n - 1] },
{ A[n - 1, n - 2], A[n - 1, n - 1] } }, n - 2, n - 1);
if (!lastMatrix.IsIdentity())
{
result.Add(lastMatrix);
}
return result;
}