private static List TwoLevelDecompose()

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;
            }