Quantum++  v1.0.0-beta1
C++11 quantum computing library
util.h
Go to the documentation of this file.
1 /*
2  * Quantum++
3  *
4  * Copyright (c) 2013 - 2016 Vlad Gheorghiu (vgheorgh@gmail.com)
5  *
6  * This file is part of Quantum++.
7  *
8  * Quantum++ is free software: you can redistribute it and/or modify
9  * it under the terms of the GNU General Public License as published by
10  * the Free Software Foundation, either version 3 of the License, or
11  * (at your option) any later version.
12  *
13  * Quantum++ is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the GNU General Public License
19  * along with Quantum++. If not, see <http://www.gnu.org/licenses/>.
20  */
21 
27 #ifndef INTERNAL_UTIL_H_
28 #define INTERNAL_UTIL_H_
29 
30 // silence g++4.9 bogus warning -Warray-bounds and -Wmaybe-uninitialized
31 // in qpp::util::_multiidx2n()
32 #if (__GNUC__ && !__clang__)
33 #pragma GCC diagnostic ignored "-Warray-bounds"
34 #pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
35 #endif
36 
37 namespace qpp
38 {
43 namespace internal
44 {
45 
46 // integer index to multi-index, use C-style array for speed
47 // standard lexicographical order, e.g. 00, 01, 10, 11
48 inline void n2multiidx(idx n, idx numdims, const idx* dims, idx* result)
49 noexcept
50 {
51  // no error checks to improve speed
52  for (idx i = 0; i < numdims; ++i)
53  {
54  result[numdims - i - 1] = n % (dims[numdims - i - 1]);
55  n /= (dims[numdims - i - 1]);
56  }
57 }
58 
59 // multi-index to integer index, use C-style array for speed,
60 // standard lexicographical order, e.g. 00->0, 01->1, 10->2, 11->3
61 inline idx multiidx2n(const idx* midx, idx numdims, const idx* dims)
62 noexcept
63 {
64  // no error checks to improve speed
65 
66  // Static allocation for speed!
67  // double the size for matrices reshaped as vectors
68  idx part_prod[2 * maxn];
69 
70  idx result = 0;
71  part_prod[numdims - 1] = 1;
72  for (idx i = 1; i < numdims; ++i)
73  {
74  part_prod[numdims - i - 1] =
75  part_prod[numdims - i] * dims[numdims - i];
76  result += midx[numdims - i - 1] * part_prod[numdims - i - 1];
77  }
78 
79  return result + midx[numdims - 1];
80 }
81 
82 // check square matrix
83 template<typename Derived>
84 bool check_square_mat(const Eigen::MatrixBase<Derived>& A)
85 {
86  return A.rows() == A.cols();
87 }
88 
89 // check whether input is a vector or not
90 template<typename Derived>
91 bool check_vector(const Eigen::MatrixBase<Derived>& A)
92 {
93  return A.rows() == 1 || A.cols() == 1;
94 }
95 
96 // check whether input is a row vector or not
97 template<typename Derived>
98 bool check_rvector(const Eigen::MatrixBase<Derived>& A)
99 {
100  return A.rows() == 1;
101 }
102 
103 // check whether input is a column vector or not
104 template<typename Derived>
105 bool check_cvector(const Eigen::MatrixBase<Derived>& A)
106 {
107  return A.cols() == 1;
108 }
109 
110 // check non-zero size of object that supports size() function
111 template<typename T>
112 bool check_nonzero_size(const T& x) noexcept
113 {
114  return x.size() != 0;
115 }
116 
117 // check that all sizes match
118 template<typename T1, typename T2>
119 bool check_matching_sizes(const T1& lhs, const T2& rhs) noexcept
120 {
121  return lhs.size() == rhs.size();
122 }
123 
124 // check that dims is a valid dimension vector
125 inline bool check_dims(const std::vector<idx>& dims)
126 {
127  if (dims.size() == 0)
128  return false;
129 
130  return std::find_if(std::begin(dims), std::end(dims),
131  [dims](idx i) -> bool
132  {
133  if (i == 0) return true;
134  else return false;
135  }) == std::end(dims);
136 }
137 
138 // check that valid dims match the dimensions
139 // of valid (non-zero sized) square matrix
140 template<typename Derived>
141 bool check_dims_match_mat(const std::vector<idx>& dims,
142  const Eigen::MatrixBase<Derived>& A)
143 {
144  idx proddim = std::accumulate(std::begin(dims), std::end(dims),
145  static_cast<idx>(1), std::multiplies<idx>());
146 
147  return proddim == static_cast<idx>(A.rows());
148 }
149 
150 // check that valid dims match the dimensions of valid column vector
151 template<typename Derived>
152 bool check_dims_match_cvect(const std::vector<idx>& dims,
153  const Eigen::MatrixBase<Derived>& V)
154 {
155  idx proddim = std::accumulate(std::begin(dims), std::end(dims),
156  static_cast<idx>(1), std::multiplies<idx>());
157 
158  return proddim == static_cast<idx>(V.rows());
159 }
160 
161 // check that valid dims match the dimensions of valid row vector
162 template<typename Derived>
163 bool check_dims_match_rvect(const std::vector<idx>& dims,
164  const Eigen::MatrixBase<Derived>& V)
165 {
166  idx proddim = std::accumulate(std::begin(dims), std::end(dims),
167  static_cast<idx>(1), std::multiplies<idx>());;
168 
169  return proddim == static_cast<idx>(V.cols());
170 }
171 
172 // check that all elements in valid dims equal to dim
173 inline bool check_eq_dims(const std::vector<idx>& dims, idx dim) noexcept
174 {
175  for (idx i : dims)
176  if (i != dim)
177  return false;
178 
179  return true;
180 }
181 
182 // check that subsys is valid with respect to valid dims
183 inline bool check_subsys_match_dims(const std::vector<idx>& subsys,
184  const std::vector<idx>& dims)
185 {
186  // // check non-zero sized subsystems
187  // if (subsys.size() == 0)
188  // return false;
189 
190  // check valid number of subsystems
191  if (subsys.size() > dims.size())
192  return false;
193 
194  // sort the subsystems
195  std::vector<idx> subsyssort = subsys;
196  std::sort(std::begin(subsyssort), std::end(subsyssort));
197 
198  // check duplicates
199  if (std::unique(std::begin(subsyssort), std::end(subsyssort))
200  != std::end(subsyssort))
201  return false;
202 
203  // check range of subsystems
204  return std::find_if(std::begin(subsyssort), std::end(subsyssort),
205  [dims](idx i) -> bool
206  {
207  return i > dims.size() - 1;
208  }) == std::end(subsyssort);
209 }
210 
211 // check matrix is 2 x 2
212 template<typename Derived>
213 bool check_qubit_matrix(const Eigen::MatrixBase<Derived>& A) noexcept
214 {
215  return A.rows() == 2 && A.cols() == 2;
216 }
217 
218 // check column vector is 2 x 1
219 template<typename Derived>
220 bool check_qubit_cvector(const Eigen::MatrixBase<Derived>& V) noexcept
221 {
222  return V.rows() == 2 && V.cols() == 1;
223 }
224 
225 // check row vector is 1 x 2
226 template<typename Derived>
227 bool check_qubit_rvector(const Eigen::MatrixBase<Derived>& V) noexcept
228 {
229  return V.rows() == 1 && V.cols() == 2;
230 }
231 
232 // check row vector is 1 x 2 or 2 x 1
233 template<typename Derived>
234 bool check_qubit_vector(const Eigen::MatrixBase<Derived>& V) noexcept
235 {
236  return (V.rows() == 1 && V.cols() == 2) ||
237  (V.rows() == 2 && V.cols() == 1);
238 }
239 
240 
241 // check valid permutation
242 inline bool check_perm(const std::vector<idx>& perm)
243 {
244  if (perm.size() == 0)
245  return false;
246 
247  std::vector<idx> ordered(perm.size());
248  std::iota(std::begin(ordered), std::end(ordered), 0);
249 
250  return std::is_permutation(std::begin(ordered), std::end(ordered),
251  std::begin(perm));
252 }
253 
254 // Kronecker product of 2 matrices, preserve return type
255 // internal function for the variadic template function wrapper kron()
256 template<typename Derived1, typename Derived2>
257 dyn_mat<typename Derived1::Scalar> kron2(const Eigen::MatrixBase<Derived1>& A,
258  const Eigen::MatrixBase<Derived2>& B)
259 {
260  const dyn_mat<typename Derived1::Scalar>& rA = A.derived();
261  const dyn_mat<typename Derived2::Scalar>& rB = B.derived();
262 
263  // EXCEPTION CHECKS
264 
265  // check types
266  if (!std::is_same<typename Derived1::Scalar,
267  typename Derived2::Scalar>::value)
268  throw Exception("qpp::kron()", Exception::Type::TYPE_MISMATCH);
269 
270  // check zero-size
272  throw Exception("qpp::kron()", Exception::Type::ZERO_SIZE);
273 
274  // check zero-size
276  throw Exception("qpp::kron()", Exception::Type::ZERO_SIZE);
277  // END EXCEPTION CHECKS
278 
279  idx Acols = static_cast<idx>(rA.cols());
280  idx Arows = static_cast<idx>(rA.rows());
281  idx Bcols = static_cast<idx>(rB.cols());
282  idx Brows = static_cast<idx>(rB.rows());
283 
285  result.resize(Arows * Brows, Acols * Bcols);
286 
287 #ifdef WITH_OPENMP_
288 #pragma omp parallel for collapse(2)
289 #endif // WITH_OPENMP_
290  for (idx j = 0; j < Acols; ++j) // column major order for speed
291  for (idx i = 0; i < Arows; ++i)
292  result.block(i * Brows, j * Bcols, Brows, Bcols) = rA(i, j) * rB;
293 
294  return result;
295 }
296 
297 // Direct sum of 2 matrices, preserve return type
298 // internal function for the variadic template function wrapper dirsum()
299 template<typename Derived1, typename Derived2>
301  const Eigen::MatrixBase<Derived1>& A,
302  const Eigen::MatrixBase<Derived2>& B)
303 {
304  const dyn_mat<typename Derived1::Scalar>& rA = A.derived();
305  const dyn_mat<typename Derived2::Scalar>& rB = B.derived();
306 
307  // EXCEPTION CHECKS
308 
309  // check types
310  if (!std::is_same<typename Derived1::Scalar,
311  typename Derived2::Scalar>::value)
312  throw Exception("qpp::dirsum()", Exception::Type::TYPE_MISMATCH);
313 
314  // check zero-size
316  throw Exception("qpp::dirsum()", Exception::Type::ZERO_SIZE);
317 
318  // check zero-size
320  throw Exception("qpp::dirsum()", Exception::Type::ZERO_SIZE);
321  // END EXCEPTION CHECKS
322 
323  idx Acols = static_cast<idx>(rA.cols());
324  idx Arows = static_cast<idx>(rA.rows());
325  idx Bcols = static_cast<idx>(rB.cols());
326  idx Brows = static_cast<idx>(rB.rows());
327 
330  Acols + Bcols);
331 
332  result.block(0, 0, Arows, Acols) = rA;
333  result.block(Arows, Acols, Brows, Bcols) = rB;
334 
335  return result;
336 }
337 
338 // may be useful, extracts variadic template argument pack into a std::vector
339 template<typename T>
340 // ends the recursion
341 void variadic_vector_emplace(std::vector<T>&)
342 {
343 }
344 
345 // may be useful, extracts variadic template argument pack into a std::vector
346 template<typename T, typename First, typename ... Args>
347 void variadic_vector_emplace(std::vector<T>& v, First&& first, Args&& ... args)
348 {
349  v.emplace_back(std::forward<First>(first));
350  variadic_vector_emplace(v, std::forward<Args>(args)...);
351 }
352 
353 // returns the number of subsystems (each subsystem assumed of the same
354 // dimension d) from an object (ket/bra/density matrix) of size sz
355 inline idx get_num_subsys(idx sz, idx d)
356 {
357  return static_cast<idx>(std::llround(std::log2(sz) / std::log2(d)));
358 }
359 
360 // returns the dimension of a subsystem (each subsystem assumed of the same
361 // dimension d) from an object (ket/bra/density matrix) of size sz consisting
362 // of N subsystems
363 inline idx get_dim_subsys(idx sz, idx N)
364 {
365  if (N == 2)
366  return static_cast<idx>(std::llround(std::sqrt(sz)));
367 
368  return static_cast<idx>(std::llround(std::pow(sz, 1./N)));
369 }
370 
371 } /* namespace internal */
372 } /* namespace qpp */
373 
374 #endif /* INTERNAL_UTIL_H_ */
bool check_nonzero_size(const T &x) noexcept
Definition: util.h:112
constexpr idx maxn
Maximum number of allowed qu(d)its (subsystems)
Definition: constants.h:74
bool check_subsys_match_dims(const std::vector< idx > &subsys, const std::vector< idx > &dims)
Definition: util.h:183
dyn_mat< typename Derived1::Scalar > kron2(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
Definition: util.h:257
bool check_cvector(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:105
Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic > dyn_mat
Dynamic Eigen matrix over the field specified by Scalar.
Definition: types.h:78
bool check_dims_match_mat(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &A)
Definition: util.h:141
bool check_qubit_cvector(const Eigen::MatrixBase< Derived > &V) noexcept
Definition: util.h:220
Quantum++ main namespace.
Definition: codes.h:30
bool check_perm(const std::vector< idx > &perm)
Definition: util.h:242
idx get_num_subsys(idx sz, idx d)
Definition: util.h:355
bool check_rvector(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:98
bool check_dims_match_cvect(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &V)
Definition: util.h:152
bool check_eq_dims(const std::vector< idx > &dims, idx dim) noexcept
Definition: util.h:173
bool check_vector(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:91
bool check_qubit_rvector(const Eigen::MatrixBase< Derived > &V) noexcept
Definition: util.h:227
bool check_qubit_vector(const Eigen::MatrixBase< Derived > &V) noexcept
Definition: util.h:234
bool check_matching_sizes(const T1 &lhs, const T2 &rhs) noexcept
Definition: util.h:119
Generates custom exceptions, used when validating function parameters.
Definition: exception.h:39
void n2multiidx(idx n, idx numdims, const idx *dims, idx *result) noexcept
Definition: util.h:48
void variadic_vector_emplace(std::vector< T > &)
Definition: util.h:341
idx get_dim_subsys(idx sz, idx N)
Definition: util.h:363
bool check_square_mat(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:84
dyn_mat< typename Derived1::Scalar > dirsum2(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
Definition: util.h:300
bool check_dims_match_rvect(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &V)
Definition: util.h:163
bool check_dims(const std::vector< idx > &dims)
Definition: util.h:125
std::size_t idx
Non-negative integer index.
Definition: types.h:36
bool check_qubit_matrix(const Eigen::MatrixBase< Derived > &A) noexcept
Definition: util.h:213
idx multiidx2n(const idx *midx, idx numdims, const idx *dims) noexcept
Definition: util.h:61