GCC Code Coverage Report

Directory: ./
File: Utils/CosSinDecomposition.cpp
Date: 2022-10-15 05:10:18
Exec Total Coverage
Lines: 30 41 73.2%
Functions: 1 1 100.0%
Branches: 72 158 45.6%
Decisions: 11 20 55.0%
Line Branch Decision Exec Source
1 // Copyright 2019-2022 Cambridge Quantum Computing
2 //
3 // Licensed under the Apache License, Version 2.0 (the "License");
4 // you may not use this file except in compliance with the License.
5 // You may obtain a copy of the License at
6 //
7 // http://www.apache.org/licenses/LICENSE-2.0
8 //
9 // Unless required by applicable law or agreed to in writing, software
10 // distributed under the License is distributed on an "AS IS" BASIS,
11 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12 // See the License for the specific language governing permissions and
13 // limitations under the License.
14
15 #include "CosSinDecomposition.hpp"
16
17 #include <cmath>
18 #include <tklog/TketLog.hpp>
19
20 #include "Constants.hpp"
21 #include "EigenConfig.hpp"
22 #include "MatrixAnalysis.hpp"
23
24 namespace tket {
25
26 874 csd_t CS_decomp(const Eigen::MatrixXcd &u) {
27
2/4
✓ Branch 1 taken 874 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 874 times.
1/2
✗ Decision 'true' not taken.
✓ Decision 'false' taken 874 times.
 874 if (!is_unitary(u)) {
28
0/2
✗ Decision 'true' not taken.
✗ Decision 'false' not taken.
 ✗ throw std::invalid_argument("Matrix for CS decomposition is not unitary");
29 }
30 874 unsigned N = u.rows();
31
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 874 times.
1/2
✗ Decision 'true' not taken.
✓ Decision 'false' taken 874 times.
 874 if (N % 2 != 0) {
32 ✗ throw std::invalid_argument(
33 ✗ "Matrix for CS decomposition has odd dimensions");
34 }
35 874 unsigned n = N / 2;
36
2/4
✓ Branch 1 taken 874 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 874 times.
✗ Branch 5 not taken.
 874 Eigen::MatrixXcd u00 = u.topLeftCorner(n, n);
37
2/4
✓ Branch 1 taken 874 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 874 times.
✗ Branch 5 not taken.
 874 Eigen::MatrixXcd u01 = u.topRightCorner(n, n);
38
2/4
✓ Branch 1 taken 874 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 874 times.
✗ Branch 5 not taken.
 874 Eigen::MatrixXcd u10 = u.bottomLeftCorner(n, n);
39
2/4
✓ Branch 1 taken 874 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 874 times.
✗ Branch 5 not taken.
 874 Eigen::MatrixXcd u11 = u.bottomRightCorner(n, n);
40
41 Eigen::JacobiSVD<Eigen::MatrixXcd, Eigen::NoQRPreconditioner> svd(
42
1/2
✓ Branch 1 taken 874 times.
✗ Branch 2 not taken.
 874 u00, Eigen::ComputeFullU | Eigen::ComputeFullV);
43
4/8
✓ Branch 1 taken 874 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 874 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 874 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 874 times.
✗ Branch 11 not taken.
 874 Eigen::MatrixXcd l0 = svd.matrixU().rowwise().reverse();
44
4/8
✓ Branch 1 taken 874 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 874 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 874 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 874 times.
✗ Branch 11 not taken.
 874 Eigen::MatrixXcd r0_dag = svd.matrixV().rowwise().reverse();
45
4/8
✓ Branch 1 taken 874 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 874 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 874 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 874 times.
✗ Branch 11 not taken.
 874 Eigen::MatrixXd c = svd.singularValues().reverse().asDiagonal();
46
2/4
✓ Branch 1 taken 874 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 874 times.
✗ Branch 5 not taken.
 874 Eigen::MatrixXcd r0 = r0_dag.adjoint();
47
48 // Now u00 = l0 c r0; l0 and r0 are unitary, and c is diagonal with positive
49 // non-decreasing entries. Because u00 is a submatrix of a unitary matrix, its
50 // singular values (the entries of c) are all <= 1.
51
52
2/4
✓ Branch 1 taken 874 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 874 times.
✗ Branch 5 not taken.
 874 Eigen::HouseholderQR<Eigen::MatrixXcd> qr = (u10 * r0_dag).householderQr();
53
2/4
✓ Branch 1 taken 874 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 874 times.
✗ Branch 5 not taken.
 874 Eigen::MatrixXcd l1 = qr.householderQ();
54
2/4
✓ Branch 2 taken 874 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 874 times.
✗ Branch 6 not taken.
 874 Eigen::MatrixXcd S = qr.matrixQR().triangularView<Eigen::Upper>();
55
56 // Now u10 r0* = l1 S; l1 is unitary, and S is upper triangular.
57 //
58 // Claim: S is diagonal.
59 // Proof: Since u is unitary, we have
60 // I = u00* u00 + u10* u10
61 // = (l0 c r0)* (l0 c r0) + (l1 S r0)* (l1 S r0)
62 // = r0* c l0* l0 c r0 + r0* S* l1* l1 S r0
63 // = r0* c^2 r0 + r0* S* S r0
64 // = r0* (c^2 + S* S) r0
65 // So I = c^2 + (S* S), so (S* S) = I - c^2 is a diagonal matrix with non-
66 // increasing entries in the range [0,1). As S is upper triangular, this
67 // implies that S must be diagonal. (Proof by induction over the dimension n
68 // of S: consider the two cases S_00 = 0 and S_00 != 0 and reduce to the n-1
69 // case.)
70 //
71 // We want S to be real. This is not guaranteed, though since it is diagonal
72 // it can be made so by adjusting l1. In fact, it seems that Eigen always does
73 // give us a real S. We will not assume this, but will log an informational
74 // message and do the adjustment ourselves if that behaviour does change.
75
76
3/6
✓ Branch 1 taken 874 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 874 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 874 times.
1/2
✗ Decision 'true' not taken.
✓ Decision 'false' taken 874 times.
 874 if (!S.imag().isZero()) {
77 ✗ tket_log()->info(
78 "Eigen surprisingly returned a non-real diagonal R in QR "
79 "decomposition; adjusting Q and R to make it real.");
80
0/2
✗ Decision 'true' not taken.
✗ Decision 'false' not taken.
 ✗ for (unsigned j = 0; j < n; j++) {
81 ✗ std::complex<double> z = S(j, j);
82 ✗ double r = std::abs(z);
83
0/2
✗ Decision 'true' not taken.
✗ Decision 'false' not taken.
 ✗ if (r > EPS) {
84 ✗ std::complex<double> w = std::conj(z) / r;
85 ✗ S(j, j) *= w;
86 ✗ l1.col(j) /= w;
87 }
88 }
89 }
90
91 // Now S is real and diagonal, and c^2 + S^2 = I.
92
93
2/4
✓ Branch 1 taken 874 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 874 times.
✗ Branch 5 not taken.
 874 Eigen::MatrixXd s = S.real();
94
95 // Make all entries in s non-negative.
96
2/2
✓ Branch 0 taken 2341 times.
✓ Branch 1 taken 874 times.
2/2
✓ Decision 'true' taken 2341 times.
✓ Decision 'false' taken 874 times.
 3215 for (unsigned j = 0; j < n; j++) {
97
3/4
✓ Branch 1 taken 2341 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1030 times.
✓ Branch 4 taken 1311 times.
2/2
✓ Decision 'true' taken 1030 times.
✓ Decision 'false' taken 1311 times.
 2341 if (s(j, j) < 0) {
98
2/4
✓ Branch 1 taken 1030 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1030 times.
✗ Branch 5 not taken.
 1030 s(j, j) = -s(j, j);
99
4/8
✓ Branch 1 taken 1030 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1030 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1030 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1030 times.
✗ Branch 11 not taken.
 1030 l1.col(j) = -l1.col(j);
100 }
101 }
102
103 // Finally compute r1, being careful not to divide by small things.
104
2/4
✓ Branch 1 taken 874 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 874 times.
✗ Branch 5 not taken.
 874 Eigen::MatrixXcd r1 = Eigen::MatrixXcd::Zero(n, n);
105
2/2
✓ Branch 0 taken 2341 times.
✓ Branch 1 taken 874 times.
2/2
✓ Decision 'true' taken 2341 times.
✓ Decision 'false' taken 874 times.
 3215 for (unsigned i = 0; i < n; i++) {
106
4/6
✓ Branch 1 taken 2341 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2341 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 989 times.
✓ Branch 7 taken 1352 times.
2/2
✓ Decision 'true' taken 989 times.
✓ Decision 'false' taken 1352 times.
 2341 if (s(i, i) > c(i, i)) {
107
8/16
✓ Branch 1 taken 989 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 989 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 989 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 989 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 989 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 989 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 989 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 989 times.
✗ Branch 23 not taken.
 989 r1.row(i) = -(l0.adjoint() * u01).row(i) / s(i, i);
108 } else {
109
7/14
✓ Branch 1 taken 1352 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1352 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1352 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1352 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 1352 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 1352 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 1352 times.
✗ Branch 20 not taken.
 1352 r1.row(i) = (l1.adjoint() * u11).row(i) / c(i, i);
110 }
111 }
112
1/2
✓ Branch 1 taken 874 times.
✗ Branch 2 not taken.
 1748 return {l0, l1, r0, r1, c, s};
113 874 }
114
115 } // namespace tket
116