Quantum++  v0.7
C++11 quantum computing library
util.h
Go to the documentation of this file.
1 /*
2  * Quantum++
3  *
4  * Copyright (c) 2013 - 2015 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++ bogus warning -Warray-bounds in _multiidx2n() in line
31 // part_prod[numdims - 1] = 1;
32 // #if (__GNUC__)
33 // #pragma GCC diagnostic ignored "-Warray-bounds"
34 // #endif
35 
36 namespace qpp
37 {
42 namespace internal
43 {
44 
45 // integer index to multi-index, use C-style array for speed
46 // standard lexicographical order, e.g. 00, 01, 10, 11
47 inline void _n2multiidx(idx n, idx numdims, const idx* dims, idx* result)
48 noexcept
49 {
50  // no error checks to improve speed
51  for (idx i = 0; i < numdims; ++i)
52  {
53  result[numdims - i - 1] = n % (dims[numdims - i - 1]);
54  n /= (dims[numdims - i - 1]);
55  }
56 }
57 
58 // multi-index to integer index, use C-style array for speed,
59 // standard lexicographical order, e.g. 00->0, 01->1, 10->2, 11->3
60 inline idx _multiidx2n(const idx* midx, idx numdims, const idx* dims)
61 noexcept
62 {
63  // no error checks to improve speed
64 
65  // Static allocation for speed!
66  // double the size for matrices reshaped as vectors
67  idx part_prod[2 * maxn];
68 
69  idx result = 0;
70  part_prod[numdims - 1] = 1;
71  for (idx i = 1; i < numdims; ++i)
72  {
73  part_prod[numdims - i - 1] =
74  part_prod[numdims - i] * dims[numdims - i];
75  result += midx[numdims - i - 1] * part_prod[numdims - i - 1];
76  }
77 
78  return result + midx[numdims - 1];
79 }
80 
81 // check square matrix
82 template<typename Derived>
83 bool _check_square_mat(const Eigen::MatrixBase<Derived>& A)
84 {
86 
87  return rA.rows() == rA.cols();
88 }
89 
90 // check whether input is a vector or not
91 template<typename Derived>
92 bool _check_vector(const Eigen::MatrixBase<Derived>& A)
93 {
95 
96  return rA.rows() == 1 || rA.cols() == 1;
97 }
98 
99 // check whether input is a row vector or not
100 template<typename Derived>
101 bool _check_rvector(const Eigen::MatrixBase<Derived>& A)
102 {
103  const dyn_mat<typename Derived::Scalar>& rA = A;
104 
105  return rA.rows() == 1;
106 }
107 
108 // check whether input is a column vector or not
109 template<typename Derived>
110 bool _check_cvector(const Eigen::MatrixBase<Derived>& A)
111 {
112  const dyn_mat<typename Derived::Scalar>& rA = A;
113 
114  return rA.cols() == 1;
115 }
116 
117 // check non-zero size of object that supports size() function
118 template<typename T>
119 bool _check_nonzero_size(const T& x) noexcept
120 {
121  return x.size() != 0;
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) quare matrix
140 template<typename Derived>
141 bool _check_dims_match_mat(const std::vector<idx>& dims,
142  const Eigen::MatrixBase<Derived>& A)
143 {
144  const dyn_mat<typename Derived::Scalar>& rA = A;
145 
146  idx proddim = std::accumulate(std::begin(dims), std::end(dims),
147  static_cast<idx>(1), std::multiplies<idx>());
148 
149  return proddim == static_cast<idx>(rA.rows());
150 }
151 
152 // check that valid dims match the dimensions of valid column vector
153 template<typename Derived>
154 bool _check_dims_match_cvect(const std::vector<idx>& dims,
155  const Eigen::MatrixBase<Derived>& V)
156 {
157  const dyn_mat<typename Derived::Scalar>& rV = V;
158 
159  idx proddim = std::accumulate(std::begin(dims), std::end(dims),
160  static_cast<idx>(1), std::multiplies<idx>());
161 
162  return proddim == static_cast<idx>(rV.rows());
163 }
164 
165 // check that valid dims match the dimensions of valid row vector
166 template<typename Derived>
167 bool _check_dims_match_rvect(const std::vector<idx>& dims,
168  const Eigen::MatrixBase<Derived>& V)
169 {
170  const dyn_mat<typename Derived::Scalar>& rV = V;
171 
172  idx proddim = std::accumulate(std::begin(dims), std::end(dims),
173  static_cast<idx>(1), std::multiplies<idx>());;
174 
175  return proddim == static_cast<idx>(rV.cols());
176 }
177 
178 // check that all elements in valid dims equal to dim
179 inline bool _check_eq_dims(const std::vector<idx>& dims, idx dim) noexcept
180 {
181  for (idx i : dims)
182  if (i != dim)
183  return false;
184 
185  return true;
186 }
187 
188 // check that subsys is valid with respect to valid dims
189 inline bool _check_subsys_match_dims(const std::vector<idx>& subsys,
190  const std::vector<idx>& dims)
191 {
192  // // check non-zero sized subsystems
193  // if (subsys.size() == 0)
194  // return false;
195 
196  // check valid number of subsystems
197  if (subsys.size() > dims.size())
198  return false;
199 
200  // sort the subsystems
201  std::vector<idx> subsyssort = subsys;
202  std::sort(std::begin(subsyssort), std::end(subsyssort));
203 
204  // check duplicates
205  if (std::unique(std::begin(subsyssort), std::end(subsyssort))
206  != std::end(subsyssort))
207  return false;
208 
209  // check range of subsystems
210  return std::find_if(std::begin(subsyssort), std::end(subsyssort),
211  [dims](idx i) -> bool
212  {
213  return i > dims.size() - 1;
214  }) == std::end(subsyssort);
215 }
216 
217 // check matrix is 2 x 2
218 template<typename Derived>
219 bool _check_qubit_matrix(const Eigen::MatrixBase<Derived>& A) noexcept
220 {
221  const dyn_mat<typename Derived::Scalar>& rA = A;
222 
223  return rA.rows() == 2 && rA.cols() == 2;
224 }
225 
226 // check column vector is 2 x 1
227 template<typename Derived>
228 bool _check_qubit_cvector(const Eigen::MatrixBase<Derived>& V) noexcept
229 {
230  const dyn_mat<typename Derived::Scalar>& rV = V;
231 
232  return rV.rows() == 2 && rV.cols() == 1;
233 }
234 
235 // check row vector is 1 x 2
236 template<typename Derived>
237 bool _check_qubit_rvector(const Eigen::MatrixBase<Derived>& V) noexcept
238 {
239  const dyn_mat<typename Derived::Scalar>& rV = V;
240 
241  return rV.rows() == 1 && rV.cols() == 2;
242 }
243 
244 // check row vector is 1 x 2 or 2 x 1
245 template<typename Derived>
246 bool _check_qubit_vector(const Eigen::MatrixBase<Derived>& V) noexcept
247 {
248  const dyn_mat<typename Derived::Scalar>& rV = V;
249 
250  return (rV.rows() == 1 && rV.cols() == 2) ||
251  (rV.rows() == 2 && rV.cols() == 1);
252 }
253 
254 
255 // check valid permutation
256 inline bool _check_perm(const std::vector<idx>& perm)
257 {
258  if (perm.size() == 0)
259  return false;
260 
261  std::vector<idx> ordered(perm.size());
262  std::iota(std::begin(ordered), std::end(ordered), 0);
263 
264  return std::is_permutation(std::begin(ordered), std::end(ordered),
265  std::begin(perm));
266 }
267 
268 // Kronecker product of 2 matrices, preserve return type
269 // internal function for the variadic template function wrapper kron()
270 template<typename Derived1, typename Derived2>
271 dyn_mat<typename Derived1::Scalar> _kron2(const Eigen::MatrixBase<Derived1>& A,
272  const Eigen::MatrixBase<Derived2>& B)
273 {
276 
277  // EXCEPTION CHECKS
278 
279  // check types
280  if (!std::is_same<typename Derived1::Scalar,
281  typename Derived2::Scalar>::value)
282  throw Exception("qpp::kron()", Exception::Type::TYPE_MISMATCH);
283 
284  // check zero-size
286  throw Exception("qpp::kron()", Exception::Type::ZERO_SIZE);
287 
288  // check zero-size
290  throw Exception("qpp::kron()", Exception::Type::ZERO_SIZE);
291 
292  idx Acols = static_cast<idx>(rA.cols());
293  idx Arows = static_cast<idx>(rA.rows());
294  idx Bcols = static_cast<idx>(rB.cols());
295  idx Brows = static_cast<idx>(rB.rows());
296 
298  result.resize(Arows * Brows, Acols * Bcols);
299 
300 #pragma omp parallel for collapse(2)
301  for (idx j = 0; j < Acols; ++j) // column major order for speed
302  for (idx i = 0; i < Arows; ++i)
303  result.block(i * Brows, j * Bcols, Brows, Bcols) = rA(i, j) * rB;
304 
305  return result;
306 }
307 
308 // Direct sum of 2 matrices, preserve return type
309 // internal function for the variadic template function wrapper dirsum()
310 template<typename Derived1, typename Derived2>
312  const Eigen::MatrixBase<Derived1>& A,
313  const Eigen::MatrixBase<Derived2>& B)
314 {
317 
318  // EXCEPTION CHECKS
319 
320  // check types
321  if (!std::is_same<typename Derived1::Scalar,
322  typename Derived2::Scalar>::value)
323  throw Exception("qpp::dirsum()", Exception::Type::TYPE_MISMATCH);
324 
325  // check zero-size
327  throw Exception("qpp::dirsum()", Exception::Type::ZERO_SIZE);
328 
329  // check zero-size
331  throw Exception("qpp::dirsum()", Exception::Type::ZERO_SIZE);
332 
333  idx Acols = static_cast<idx>(rA.cols());
334  idx Arows = static_cast<idx>(rA.rows());
335  idx Bcols = static_cast<idx>(rB.cols());
336  idx Brows = static_cast<idx>(rB.rows());
337 
340  Acols + Bcols);
341 
342  result.block(0, 0, Arows, Acols) = rA;
343  result.block(Arows, Acols, Brows, Bcols) = rB;
344 
345  return result;
346 }
347 
348 // may be useful, extracts variadic template argument pack into a std::vector
349 template<typename T>
350 // ends the recursion
351 void variadic_vector_emplace(std::vector<T>&)
352 {
353 }
354 
355 // may be useful, extracts variadic template argument pack into a std::vector
356 template<typename T, typename First, typename ... Args>
357 void variadic_vector_emplace(std::vector<T>& v, First&& first, Args&& ... args)
358 {
359  v.emplace_back(std::forward<First>(first));
360  variadic_vector_emplace(v, std::forward<Args>(args)...);
361 }
362 
363 } /* namespace internal */
364 } /* namespace qpp */
365 
366 #endif /* INTERNAL_UTIL_H_ */
bool _check_cvector(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:110
bool _check_qubit_cvector(const Eigen::MatrixBase< Derived > &V) noexcept
Definition: util.h:228
constexpr idx maxn
Maximum number of allowed qu(d)its (subsystems)
Definition: constants.h:82
bool _check_subsys_match_dims(const std::vector< idx > &subsys, const std::vector< idx > &dims)
Definition: util.h:189
bool _check_dims_match_mat(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &A)
Definition: util.h:141
bool _check_dims_match_cvect(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &V)
Definition: util.h:154
dyn_mat< typename Derived1::Scalar > _kron2(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
Definition: util.h:271
Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic > dyn_mat
Dynamic Eigen matrix over the field specified by Scalar.
Definition: types.h:73
bool _check_vector(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:92
Quantum++ main namespace.
Definition: codes.h:30
bool _check_dims_match_rvect(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &V)
Definition: util.h:167
bool _check_square_mat(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:83
bool _check_rvector(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:101
bool _check_qubit_matrix(const Eigen::MatrixBase< Derived > &A) noexcept
Definition: util.h:219
bool _check_perm(const std::vector< idx > &perm)
Definition: util.h:256
idx _multiidx2n(const idx *midx, idx numdims, const idx *dims) noexcept
Definition: util.h:60
Generates custom exceptions, used when validating function parameters.
Definition: exception.h:39
void variadic_vector_emplace(std::vector< T > &)
Definition: util.h:351
bool _check_nonzero_size(const T &x) noexcept
Definition: util.h:119
dyn_mat< typename Derived1::Scalar > _dirsum2(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
Definition: util.h:311
bool _check_qubit_vector(const Eigen::MatrixBase< Derived > &V) noexcept
Definition: util.h:246
void _n2multiidx(idx n, idx numdims, const idx *dims, idx *result) noexcept
Definition: util.h:47
std::size_t idx
Non-negative integer index.
Definition: types.h:36
bool _check_dims(const std::vector< idx > &dims)
Definition: util.h:125
bool _check_eq_dims(const std::vector< idx > &dims, idx dim) noexcept
Definition: util.h:179
bool _check_qubit_rvector(const Eigen::MatrixBase< Derived > &V) noexcept
Definition: util.h:237