scnc
/
qiskit-aer
mirror of https://github.com/Qiskit/qiskit-aer.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Add hermitian eigensolver (#889) 

Co-authored-by: Juan Gomez <atilag@gmail.com>
Co-authored-by: Christopher J. Wood <cjwood@us.ibm.com>
Co-authored-by: Eric Eastman <ereastman96@gmail.com>
This commit is contained in:
Victor Villar 2020-09-14 18:58:09 +02:00 committed by GitHub
parent 88fddc313d
commit 801e88efe9
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
14 changed files with 854 additions and 174 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

21
.github/workflows/main.yml vendored
View File

@ -64,6 +64,13 @@ jobs:
python-version: 3.7
- name: Install deps
run: pip install conan
- name: Install openblas
run: |
set -e
sudo apt-get update
sudo apt-get install -y libopenblas-dev
shell: bash
if: runner.os == 'Linux'
- name: Add msbuild to PATH
uses: microsoft/setup-msbuild@v1.0.1
if: runner.os == 'Windows'
@ -99,6 +106,13 @@ jobs:
python-version: ${{ matrix.python-version }}
- name: Install Deps
run: python -m pip install -U setuptools wheel virtualenv
- name: Install openblas
run: |
set -e
sudo apt-get update
sudo apt-get install -y libopenblas-dev
shell: bash
if: runner.os == 'Linux'
- name: Build Sdist
run: python setup.py sdist
- name: Install from sdist and test
@ -156,6 +170,13 @@ jobs:
if: runner.os == 'Windows'
- name: Install Deps
run: python -m pip install -U -r requirements-dev.txt git+https://github.com/Qiskit/qiskit-terra
- name: Install openblas
run: |
set -e
sudo apt-get update
sudo apt-get install -y libopenblas-dev
shell: bash
if: runner.os == 'Linux'
- name: Install Aer
run: python -m pip install -U .
if: runner.os != 'Windows'

5
CMakeLists.txt
View File

@ -340,8 +340,5 @@ endif()
# Tests
if(BUILD_TESTS)
# These tests were meant to be as an example, or template for C++ specific
# code, but they're not being maintained, so we are disabling them until we
# have real C++ tests in the codebase.
# add_subdirectory(test)
add_subdirectory(test)
endif()

16
CONTRIBUTING.md
View File

@ -724,6 +724,22 @@ Manual for `stestr` can be found [here](https://stestr.readthedocs.io/en/latest/
The integration tests for Qiskit python extension are included in: `test/terra`.
## C++ Tests
Our C++ unit tests use the Catch2 framework, an include-only C++ unit-testing framework.
Catch2 framework documentation can be found [here](https://github.com/catchorg/Catch2).
Then, in any case, build Aer with the extra cmake argument BUILD_TESTS set to true:
```
python ./setup.py bdist_wheel --build-type=Debug -- -DBUILD_TESTS=True -- -j4 2>&1 |tee build.log
```
The test executable will be placed into the source test directory and can be run by:
```
qiskit-aer$ ./test/unitc_tests [Catch2-options]
```
## Platform support
Bear in mind that every new feature/change needs to be compatible with all our

163
src/framework/blas_protos.hpp Executable file
View File

@ -0,0 +1,163 @@
/**
* This code is part of Qiskit.
*
* (C) Copyright IBM 2018, 2019, 2020.
*
* This code is licensed under the Apache License, Version 2.0. You may
* obtain a copy of this license in the LICENSE.txt file in the root directory
* of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
*
* Any modifications or derivative works of this code must retain this
* copyright notice, and modified files need to carry a notice indicating
* that they have been altered from the originals.
*/
// Dependencies: BLAS
// These are the declarations for the various high-performance matrix routines
// used by the matrix class. An openblas install is required.
#ifndef _aer_framework_blas_protos_hpp
#define _aer_framework_blas_protos_hpp
#include <complex>
#include <iostream>
#include <vector>
#include <array>
#ifdef __cplusplus
extern "C" {
#endif
//===========================================================================
// Prototypes for level 3 BLAS
//===========================================================================
// Single-Precison Real Matrix-Vector Multiplcation
void sgemv_(const char *TransA, const size_t *M, const size_t *N,
const float *alpha, const float *A, const size_t *lda,
const float *x, const size_t *incx, const float *beta, float *y,
const size_t *lincy);
// Double-Precison Real Matrix-Vector Multiplcation
void dgemv_(const char *TransA, const size_t *M, const size_t *N,
const double *alpha, const double *A, const size_t *lda,
const double *x, const size_t *incx, const double *beta, double *y,
const size_t *lincy);
// Single-Precison Complex Matrix-Vector Multiplcation
void cgemv_(const char *TransA, const size_t *M, const size_t *N,
const std::complex<float> *alpha, const std::complex<float> *A,
const size_t *lda, const std::complex<float> *x, const size_t *incx,
const std::complex<float> *beta, std::complex<float> *y,
const size_t *lincy);
// Double-Precison Real Matrix-Vector Multiplcation
void zgemv_(const char *TransA, const size_t *M, const size_t *N,
const std::complex<double> *alpha, const std::complex<double> *A,
const size_t *lda, const std::complex<double> *x,
const size_t *incx, const std::complex<double> *beta,
std::complex<double> *y, const size_t *lincy);
// Single-Precison Real Matrix-Matrix Multiplcation
void sgemm_(const char *TransA, const char *TransB, const size_t *M,
const size_t *N, const size_t *K, const float *alpha,
const float *A, const size_t *lda, const float *B,
const size_t *lba, const float *beta, float *C, size_t *ldc);
// Double-Precison Real Matrix-Matrix Multiplcation
void dgemm_(const char *TransA, const char *TransB, const size_t *M,
const size_t *N, const size_t *K, const double *alpha,
const double *A, const size_t *lda, const double *B,
const size_t *lba, const double *beta, double *C, size_t *ldc);
// Single-Precison Complex Matrix-Matrix Multiplcation
void cgemm_(const char *TransA, const char *TransB, const size_t *M,
const size_t *N, const size_t *K, const std::complex<float> *alpha,
const std::complex<float> *A, const size_t *lda,
const std::complex<float> *B, const size_t *ldb,
const std::complex<float> *beta, std::complex<float> *C,
size_t *ldc);
// Double-Precison Complex Matrix-Matrix Multiplcation
void zgemm_(const char *TransA, const char *TransB, const size_t *M,
const size_t *N, const size_t *K, const std::complex<double> *alpha,
const std::complex<double> *A, const size_t *lda,
const std::complex<double> *B, const size_t *ldb,
const std::complex<double> *beta, std::complex<double> *C,
size_t *ldc);
// Reduces a Single-Precison Complex Hermitian matrix A to real symmetric tridiagonal form
void chetrd_(char *TRANS, int *N, std::complex<float> *A,
int *LDA, float *d, float *e, std::complex<float> *tau,
std::complex<float> *work, int *lwork, int *info);
// Reduces a Double-Precison Complex Hermitian matrix A to real symmetric tridiagonal form T
void zhetrd_(char *TRANS, int *N, std::complex<double> *A,
int *LDA, double *d, double *e, std::complex<double> *tau,
std::complex<double> *work, int *lwork, int *info);
// Computes all eigenvalues and, optionally, eigenvectors of a
// Single-Precison Complex symmetric positive definite tridiagonal matrix
void cpteqr_(char* compz, int *n, float *d, float *e,
std::complex<float> *z, int* ldz,
std::complex<float> *work, int *info);
// Computes all eigenvalues and, optionally, eigenvectors of a
// Double-Precison Complex symmetric positive definite tridiagonal matrix
void zpteqr_(char* compz, int *n, double *d, double *e,
std::complex<double> *z, int* ldz,
std::complex<double> *work, int *info);
// Computes selected eigenvalues and, optionally, eigenvectors
// of a Single-Precison Complex Hermitian matrix A
void cheevx_(char *jobz, char *range, char *uplo, int *n,
std::complex<float> *a, int *lda, float *vl,
float *vu, int *il, int *iu, float *abstol,
int *m, float *w, std::complex<float> *z, int *ldz,
std::complex<float> *work, int *lwork, float *rwork,
int *iwork, int *ifail, int *info);
// Computes selected eigenvalues and, optionally, eigenvectors
// of a Double-Precison Complex Hermitian matrix A
void zheevx_(char *jobz, char *range, char *uplo, int *n,
std::complex<double> *a, int *lda, double *vl,
double *vu, int *il, int *iu, double *abstol,
int *m, double *w, std::complex<double> *z, int *ldz,
std::complex<double> *work, int *lwork, double *rwork,
int *iwork, int *ifail, int *info);
// Determines Single-Precision machine parameters.
float slamch_(char *cmach);
// Determines Double-Precision machine parameters.
double dlamch_(char *cmach);
#ifdef __cplusplus
}
#endif
namespace AerBlas {
std::array<char, 3> Trans = {'N', 'T', 'C'};
/* Trans (input) CHARACTER*1.
On entry, TRANSA specifies the form of op( A ) to be used in the
matrix multiplication as follows:
= 'N' no transpose;
= 'T' transpose of A;
= 'C' hermitian conjugate of A.
*/
std::array<char, 2> UpLo = {'U', 'L'};
/* UpLo (input) CHARACTER*1
= 'U': Upper triangle of A is stored;
= 'L': Lower triangle of A is stored.
*/
std::array<char, 2> Jobz = {'V', 'N'};
/* Jobz (input) CHARACTER*1
= 'N': Compute eigenvalues only;
= 'V': Compute eigenvalues and eigenvectors.
*/
std::array<char, 3> Range = {'A', 'V', 'I'};
/* Range (input) CHARACTER*1
= 'A': all eigenvalues will be found.
= 'V': all eigenvalues in the half-open interval
(VL,VU] will be found.
= 'I': the IL-th through IU-th eigenvalues will
be found.
*/
} // namespace AerBlas
#endif // end _blas_protos_h_

20
src/framework/linalg/almost_equal.hpp
View File

@ -29,14 +29,22 @@ namespace Linalg {
// If we have numbers closer to 0, then max_diff can be set to a value
// way smaller than epsilon. For numbers larger than 1.0, epsilon will
// scale (the bigger the number, the bigger the epsilon).
template <typename T, typename = enable_if_numeric_t<T>>
template <typename T, typename = typename std::enable_if<std::is_arithmetic<T>::value, T>::type >
bool almost_equal(T f1, T f2,
T max_diff = std::numeric_limits<T>::epsilon(),
T max_relative_diff = std::numeric_limits<T>::epsilon()) {
T diff = std::abs(f1 - f2);
if (diff <= max_diff) return true;
return diff <=
max_relative_diff * std::max(std::abs(f1), std::abs(f2));
T diff = std::abs(f1 - f2);
if (diff <= max_diff) return true;
return diff <= max_relative_diff * std::max(std::abs(f1), std::abs(f2));
}
template <typename T>
bool almost_equal(const std::complex<T>& f1, const std::complex<T>& f2,
T max_diff = std::numeric_limits<T>::epsilon(),
T max_relative_diff = std::numeric_limits<T>::epsilon()) {
return almost_equal<T>(f1.real(), f2.real(), max_diff, max_relative_diff)
&& almost_equal<T>(f1.imag(), f2.imag(), max_diff, max_relative_diff);
}
//------------------------------------------------------------------------------
@ -44,4 +52,4 @@ bool almost_equal(T f1, T f2,
//------------------------------------------------------------------------------
} // end namespace AER
//------------------------------------------------------------------------------
#endif
#endif

101
src/framework/linalg/eigensystem.hpp Normal file
View File

@ -0,0 +1,101 @@
/**
* This code is part of Qiskit.
*
* (C) Copyright IBM 2018, 2019, 2020.
*
* This code is licensed under the Apache License, Version 2.0. You may
* obtain a copy of this license in the LICENSE.txt file in the root directory
* of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
*
* Any modifications or derivative works of this code must retain this
* copyright notice, and modified files need to carry a notice indicating
* that they have been altered from the originals.
*/
#ifndef _aer_framework_linalg_eigensystem_hpp_
#define _aer_framework_linalg_eigensystem_hpp_
#include <type_traits>
#include "framework/blas_protos.hpp"
#include "framework/matrix.hpp"
/**
* Returns the eigenvalues and eigenvectors
* of a Hermitian matrix.
* Uses the blas function ?heevx
* @param hermitian_matrix: The Hermitian matrix.
* @param eigenvalues: On output: vector with the eignevalues of the matrix (input is overwritten)
* @param eigenvectors: On output: matrix with the eigenvectors stored as columns.
*
* @returns: void
*/
template <class T>
void eigensystem_hermitian(const matrix<std::complex<T>>& hermitian_matrix,
/* out */ std::vector<T>& eigenvalues,
/* out */ matrix<std::complex<T>>& eigenvectors);
template<typename T>
struct HeevxFuncs;
template<>
struct HeevxFuncs<double>{
HeevxFuncs() = delete;
static decltype(zheevx_)& heevx;
static decltype(dlamch_)& lamch;
};
decltype(zheevx_)& HeevxFuncs<double>::heevx = zheevx_;
decltype(dlamch_)& HeevxFuncs<double>::lamch = dlamch_;
template<>
struct HeevxFuncs<float>{
HeevxFuncs() = delete;
static decltype(cheevx_)& heevx;
static decltype(slamch_)& lamch;
};
decltype(cheevx_)& HeevxFuncs<float>::heevx = cheevx_;
decltype(slamch_)& HeevxFuncs<float>::lamch = slamch_;
template <typename T>
void eigensystem_hermitian(const matrix<std::complex<T>>& hermitian_matrix,
std::vector<T>& eigenvalues,
matrix<std::complex<T>>& eigenvectors) {
if ( hermitian_matrix.GetRows() != hermitian_matrix.GetColumns() ) {
throw std::runtime_error("Input matrix in eigensystem_hermitian "
"function is not a square matrix.");
}
int n = static_cast<int>(hermitian_matrix.GetLD());
int ldz{n}, lda{n}, lwork{2*n};
int il{0}, iu{0}; // not referenced if range='A'
T vl{0.0}, vu{0.0}; // not referenced if range='A'
char cmach{'S'};
T abstol{static_cast<T>(2.0*HeevxFuncs<T>::lamch(&cmach))};
int m{0}; // number of eigenvalues found
int info{0};
eigenvectors.resize(ldz, n);
eigenvalues.clear();
eigenvalues.resize(n);
matrix<std::complex<T>> heevx_copy{hermitian_matrix};
auto work = std::vector<std::complex<T>>(lwork, {0.0, 0.0});
auto rwork = std::vector<T>(7*n, 0.0);
auto iwork = std::vector<int>(5*n, 0);
auto ifail = std::vector<int>(n, 0);
HeevxFuncs<T>::heevx(&AerBlas::Jobz[0], &AerBlas::Range[0], &AerBlas::UpLo[0], &n,
heevx_copy.data(), &lda, &vl, &vu, &il, &iu,
&abstol, &m, eigenvalues.data(), eigenvectors.data(), &ldz, work.data(),
&lwork, rwork.data(), iwork.data(), ifail.data(), &info);
if(info){
throw std::runtime_error("Something went wrong in heevx call within eigensystem_hermitian funcion. "
"Check that input matrix is really hermitian");
}
}
#endif // _aer_framework_linalg_eigensystem_hpp_

5
src/framework/linalg/linalg.hpp
View File

@ -16,13 +16,14 @@
#define _aer_framework_linalg_hpp_
#include "framework/linalg/almost_equal.hpp"
#include "framework/linalg/eigensystem.hpp"
#include "framework/linalg/linops/linops_array.hpp"
#include "framework/linalg/linops/linops_generic.hpp"
#include "framework/linalg/linops/linops_json.hpp"
#include "framework/linalg/linops/linops_map.hpp"
#include "framework/linalg/linops/linops_matrix.hpp"
#include "framework/linalg/linops/linops_unordered_map.hpp"
#include "framework/linalg/linops/linops_vector.hpp"
#include "framework/linalg/linops/linops_generic.hpp"
#include "framework/linalg/square.hpp"
#endif
#endif

142
src/framework/matrix.hpp
View File

@ -36,78 +36,8 @@ Multiplication is done with the C wrapper of the fortran blas library.
#include <vector>
#include <array>
/*******************************************************************************
*
* BLAS headers
*
******************************************************************************/
const std::array<char, 3> Trans = {'N', 'T', 'C'};
/* Trans (input) CHARACTER*1.
On entry, TRANSA specifies the form of op( A ) to be used in the
matrix multiplication as follows:
= 'N' no transpose;
= 'T' transpose of A;
= 'C' hermitian conjugate of A.
*/
#ifdef __cplusplus
extern "C" {
#endif
//===========================================================================
// Prototypes for level 3 BLAS
//===========================================================================
// Single-Precison Real Matrix-Vector Multiplcation
void sgemv_(const char *TransA, const size_t *M, const size_t *N,
const float *alpha, const float *A, const size_t *lda,
const float *x, const size_t *incx, const float *beta, float *y,
const size_t *lincy);
// Double-Precison Real Matrix-Vector Multiplcation
void dgemv_(const char *TransA, const size_t *M, const size_t *N,
const double *alpha, const double *A, const size_t *lda,
const double *x, const size_t *incx, const double *beta, double *y,
const size_t *lincy);
// Single-Precison Complex Matrix-Vector Multiplcation
void cgemv_(const char *TransA, const size_t *M, const size_t *N,
const std::complex<float> *alpha, const std::complex<float> *A,
const size_t *lda, const std::complex<float> *x, const size_t *incx,
const std::complex<float> *beta, std::complex<float> *y,
const size_t *lincy);
// Double-Precison Real Matrix-Vector Multiplcation
void zgemv_(const char *TransA, const size_t *M, const size_t *N,
const std::complex<double> *alpha, const std::complex<double> *A,
const size_t *lda, const std::complex<double> *x,
const size_t *incx, const std::complex<double> *beta,
std::complex<double> *y, const size_t *lincy);
// Single-Precison Real Matrix-Matrix Multiplcation
void sgemm_(const char *TransA, const char *TransB, const size_t *M,
const size_t *N, const size_t *K, const float *alpha,
const float *A, const size_t *lda, const float *B,
const size_t *lba, const float *beta, float *C, size_t *ldc);
// Double-Precison Real Matrix-Matrix Multiplcation
void dgemm_(const char *TransA, const char *TransB, const size_t *M,
const size_t *N, const size_t *K, const double *alpha,
const double *A, const size_t *lda, const double *B,
const size_t *lba, const double *beta, double *C, size_t *ldc);
// Single-Precison Complex Matrix-Matrix Multiplcation
void cgemm_(const char *TransA, const char *TransB, const size_t *M,
const size_t *N, const size_t *K, const std::complex<float> *alpha,
const std::complex<float> *A, const size_t *lda,
const std::complex<float> *B, const size_t *ldb,
const std::complex<float> *beta, std::complex<float> *C,
size_t *ldc);
// Double-Precison Complex Matrix-Matrix Multiplcation
void zgemm_(const char *TransA, const char *TransB, const size_t *M,
const size_t *N, const size_t *K, const std::complex<double> *alpha,
const std::complex<double> *A, const size_t *lda,
const std::complex<double> *B, const size_t *ldb,
const std::complex<double> *beta, std::complex<double> *C,
size_t *ldc);
#ifdef __cplusplus
}
#endif
#include "framework/blas_protos.hpp"
#include "framework/linalg/enable_if_numeric.hpp"
/*******************************************************************************
*
@ -268,7 +198,7 @@ public:
// Return the size of the underlying array
size_t size() const { return size_; }
// Return True if size == 0
bool empty() const { return size_ == 0; }
@ -279,10 +209,13 @@ public:
void fill(const T& val);
// Resize the matrix and reset to zero if different size
void initialize(size_t row, size_t col);
void initialize(size_t row, size_t col);
// Resize the matrix keeping current values
void resize(size_t row, size_t col);
void resize(size_t row, size_t col);
// Addressing elements by row or column
std::vector<T> row_index(size_t row) const;
std::vector<T> col_index(size_t col) const;
// overloading functions.
matrix<T> operator+(const matrix<T> &A);
@ -307,7 +240,7 @@ protected:
// size_ = rows*colums dimensions of the vector representation
// LD is the leading dimeonsion and for Column major order is in general eqaul
// to rows
// the ptr to the vector containing the matrix
T* data_ = nullptr;
};
@ -357,7 +290,7 @@ matrix<T>& matrix<T>::operator=(matrix<T>&& other) noexcept {
template <class T>
matrix<T> &matrix<T>::operator=(const matrix<T> &other) {
if (rows_ != other.rows_ || cols_ != other.cols_) {
if (rows_ != other.rows_ || cols_ != other.cols_) {
// size delete re-construct
// the matrix
free(data_);
@ -531,6 +464,37 @@ void matrix<T>::resize(size_t rows, size_t cols) {
data_ = tempmat;
}
// Addressing elements by row or column
template <class T> inline std::vector<T> matrix<T>::row_index(size_t row) const {
#ifdef DEBUG
if (row >= rows_) {
std::cerr << "Error: matrix class operator row_index out of bounds "
<< row << " >= " << rows_ << std::endl;
exit(1);
}
#endif
std::vector<T> ret;
ret.reserve(cols_);
for(size_t i = 0; i < cols_; i++)
ret.emplace_back(data_[i * rows_ + row]);
return std::move(ret);
}
template <class T> inline std::vector<T> matrix<T>::col_index(size_t col) const {
#ifdef DEBUG
if (col >= cols_) {
std::cerr << "Error: matrix class operator col_index out of bounds "
<< col << " >= " << cols_ << std::endl;
exit(1);
}
#endif
std::vector<T> ret;
ret.reserve(rows_);
// we want the elements for all rows i..rows_ and column col
for(size_t i = 0; i < rows_; i++)
ret.emplace_back(data_[col * rows_ + i]);
return std::move(ret);
}
template <class T> inline size_t matrix<T>::GetRows() const {
// returns the rows of the matrix
return rows_;
@ -712,7 +676,7 @@ inline matrix<double> operator*(const matrix<double> &A,
// C-> alpha*op(A)*op(B) +beta C
matrix<double> C(A.rows_, B.cols_);
double alpha = 1.0, beta = 0.0;
dgemm_(&Trans[0], &Trans[0], &A.rows_, &B.cols_, &A.cols_, &alpha, A.data_,
dgemm_(&AerBlas::Trans[0], &AerBlas::Trans[0], &A.rows_, &B.cols_, &A.cols_, &alpha, A.data_,
&A.LD_, B.data_, &B.LD_, &beta, C.data_, &C.LD_);
// cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, A.rows_, B.cols_,
// A.cols_, 1.0, A.data_, A.LD_, B.data_, B.LD_, 0.0, C.data_, C.LD_);
@ -724,7 +688,7 @@ inline matrix<float> operator*(const matrix<float> &A, const matrix<float> &B) {
// C-> alpha*op(A)*op(B) +beta C
matrix<float> C(A.rows_, B.cols_);
float alpha = 1.0, beta = 0.0;
sgemm_(&Trans[0], &Trans[0], &A.rows_, &B.cols_, &A.cols_, &alpha, A.data_,
sgemm_(&AerBlas::Trans[0], &AerBlas::Trans[0], &A.rows_, &B.cols_, &A.cols_, &alpha, A.data_,
&A.LD_, B.data_, &B.LD_, &beta, C.data_, &C.LD_);
// cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, A.rows_, B.cols_,
// A.cols_, 1.0, A.data_, A.LD_, B.data_, B.LD_, 0.0, C.data_, C.LD_);
@ -738,7 +702,7 @@ operator*(const matrix<std::complex<float>> &A,
// C-> alpha*op(A)*op(B) +beta C
matrix<std::complex<float>> C(A.rows_, B.cols_);
std::complex<float> alpha = 1.0, beta = 0.0;
cgemm_(&Trans[0], &Trans[0], &A.rows_, &B.cols_, &A.cols_, &alpha, A.data_,
cgemm_(&AerBlas::Trans[0], &AerBlas::Trans[0], &A.rows_, &B.cols_, &A.cols_, &alpha, A.data_,
&A.LD_, B.data_, &B.LD_, &beta, C.data_, &C.LD_);
// cblas_zgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, A.rows_, B.cols_,
// A.cols_, &alpha, A.data_, A.LD_, B.data_, B.LD_, &beta, C.data_, C.LD_);
@ -752,7 +716,7 @@ operator*(const matrix<std::complex<double>> &A,
// C-> alpha*op(A)*op(B) +beta C
matrix<std::complex<double>> C(A.rows_, B.cols_);
std::complex<double> alpha = 1.0, beta = 0.0;
zgemm_(&Trans[0], &Trans[0], &A.rows_, &B.cols_, &A.cols_, &alpha, A.data_,
zgemm_(&AerBlas::Trans[0], &AerBlas::Trans[0], &A.rows_, &B.cols_, &A.cols_, &alpha, A.data_,
&A.LD_, B.data_, &B.LD_, &beta, C.data_, &C.LD_);
// cblas_zgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, A.rows_, B.cols_,
// A.cols_, &alpha, A.data_, A.LD_, B.data_, B.LD_, &beta, C.data_, C.LD_);
@ -766,7 +730,7 @@ operator*(const matrix<float> &A, const matrix<std::complex<float>> &B) {
matrix<std::complex<float>> C(A.rows_, B.cols_), Ac(A.rows_, A.cols_);
Ac = A;
std::complex<float> alpha = 1.0, beta = 0.0;
cgemm_(&Trans[0], &Trans[0], &Ac.rows_, &B.cols_, &Ac.cols_, &alpha, Ac.data_,
cgemm_(&AerBlas::Trans[0], &AerBlas::Trans[0], &Ac.rows_, &B.cols_, &Ac.cols_, &alpha, Ac.data_,
&Ac.LD_, B.data_, &B.LD_, &beta, C.data_, &C.LD_);
// cblas_zgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, Ac.rows_, B.cols_,
// Ac.cols_, &alpha, Ac.data_, Ac.LD_, B.data_, B.LD_, &beta, C.data_, C.LD_);
@ -780,7 +744,7 @@ operator*(const matrix<double> &A, const matrix<std::complex<double>> &B) {
matrix<std::complex<double>> C(A.rows_, B.cols_), Ac(A.rows_, A.cols_);
Ac = A;
std::complex<double> alpha = 1.0, beta = 0.0;
zgemm_(&Trans[0], &Trans[0], &Ac.rows_, &B.cols_, &Ac.cols_, &alpha, Ac.data_,
zgemm_(&AerBlas::Trans[0], &AerBlas::Trans[0], &Ac.rows_, &B.cols_, &Ac.cols_, &alpha, Ac.data_,
&Ac.LD_, B.data_, &B.LD_, &beta, C.data_, &C.LD_);
// cblas_zgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, Ac.rows_, B.cols_,
// Ac.cols_, &alpha, Ac.data_, Ac.LD_, B.data_, B.LD_, &beta, C.data_, C.LD_);
@ -794,7 +758,7 @@ operator*(const matrix<std::complex<float>> &A, const matrix<float> &B) {
matrix<std::complex<float>> C(A.rows_, B.cols_), Bc(B.rows_, B.cols_);
Bc = B;
std::complex<float> alpha = 1.0, beta = 0.0;
cgemm_(&Trans[0], &Trans[0], &A.rows_, &Bc.cols_, &A.cols_, &alpha, A.data_,
cgemm_(&AerBlas::Trans[0], &AerBlas::Trans[0], &A.rows_, &Bc.cols_, &A.cols_, &alpha, A.data_,
&A.LD_, Bc.data_, &Bc.LD_, &beta, C.data_, &C.LD_);
// cblas_zgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, A.rows_, Bc.cols_,
// A.cols_, &alpha, A.data_, A.LD_, Bc.data_, Bc.LD_, &beta, C.data_, C.LD_);
@ -808,7 +772,7 @@ operator*(const matrix<std::complex<double>> &A, const matrix<double> &B) {
matrix<std::complex<double>> C(A.rows_, B.cols_), Bc(B.rows_, B.cols_);
Bc = B;
std::complex<double> alpha = 1.0, beta = 0.0;
zgemm_(&Trans[0], &Trans[0], &A.rows_, &Bc.cols_, &A.cols_, &alpha, A.data_,
zgemm_(&AerBlas::Trans[0], &AerBlas::Trans[0], &A.rows_, &Bc.cols_, &A.cols_, &alpha, A.data_,
&A.LD_, Bc.data_, &Bc.LD_, &beta, C.data_, &C.LD_);
// cblas_zgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, A.rows_, Bc.cols_,
// A.cols_, &alpha, A.data_, A.LD_, Bc.data_, Bc.LD_, &beta, C.data_, C.LD_);
@ -822,7 +786,7 @@ inline std::vector<float> operator*(const matrix<float> &A,
std::vector<float> y(A.rows_);
float alpha = 1.0, beta = 0.0;
const size_t incx = 1, incy = 1;
sgemv_(&Trans[0], &A.rows_, &A.cols_, &alpha, A.data_, &A.LD_, x.data(), &incx,
sgemv_(&AerBlas::Trans[0], &A.rows_, &A.cols_, &alpha, A.data_, &A.LD_, x.data(), &incx,
&beta, y.data(), &incy);
return y;
}
@ -833,7 +797,7 @@ inline std::vector<double> operator*(const matrix<double> &A,
std::vector<double> y(A.rows_);
double alpha = 1.0, beta = 0.0;
const size_t incx = 1, incy = 1;
dgemv_(&Trans[0], &A.rows_, &A.cols_, &alpha, A.data_, &A.LD_, x.data(), &incx,
dgemv_(&AerBlas::Trans[0], &A.rows_, &A.cols_, &alpha, A.data_, &A.LD_, x.data(), &incx,
&beta, y.data(), &incy);
return y;
}
@ -845,7 +809,7 @@ operator*(const matrix<std::complex<float>> &A,
std::vector<std::complex<float>> y(A.rows_);
std::complex<float> alpha = 1.0, beta = 0.0;
const size_t incx = 1, incy = 1;
cgemv_(&Trans[0], &A.rows_, &A.cols_, &alpha, A.data_, &A.LD_, x.data(), &incx,
cgemv_(&AerBlas::Trans[0], &A.rows_, &A.cols_, &alpha, A.data_, &A.LD_, x.data(), &incx,
&beta, y.data(), &incy);
return y;
}
@ -857,7 +821,7 @@ operator*(const matrix<std::complex<double>> &A,
std::vector<std::complex<double>> y(A.rows_);
std::complex<double> alpha = 1.0, beta = 0.0;
const size_t incx = 1, incy = 1;
zgemv_(&Trans[0], &A.rows_, &A.cols_, &alpha, A.data_, &A.LD_, x.data(), &incx,
zgemv_(&AerBlas::Trans[0], &A.rows_, &A.cols_, &alpha, A.data_, &A.LD_, x.data(), &incx,
&beta, y.data(), &incy);
return y;
}

3
src/framework/types.hpp
View File

@ -39,8 +39,11 @@ namespace AER {
using int_t = int_fast64_t;
using uint_t = uint_fast64_t;
using complex_t = std::complex<double>;
using complexf_t = std::complex<float>;
using cvector_t = std::vector<complex_t>;
using cvectorf_t = std::vector<complexf_t>;
using cmatrix_t = matrix<complex_t>;
using cmatrixf_t = matrix<complexf_t>;
using rvector_t = std::vector<double>;
using rmatrix_t = matrix<double>;
using reg_t = std::vector<uint_t>;

2
src/simulators/matrix_product_state/matrix_product_state.hpp
View File

@ -340,7 +340,7 @@ void State::initialize_qreg(uint_t num_qubits, const matrixproductstate_t &state
throw std::invalid_argument("MatrixProductState::State::initialize: initial state does not match qubit number");
}
#ifdef DEBUG
cout << "initialize with state not supported yet";
std::cout << "initialize with state not supported yet";
#endif
}

52
test/CMakeLists.txt
View File

@ -1,42 +1,20 @@
add_executable(test_snapshot "src/test_snapshot.cpp")
set_target_properties(test_snapshot PROPERTIES
LINKER_LANGUAGE CXX
CXX_STANDARD 14)
target_include_directories(test_snapshot
PRIVATE ${AER_SIMULATOR_CPP_SRC_DIR}
PRIVATE ${AER_SIMULATOR_CPP_EXTERNAL_LIBS})
target_link_libraries(test_snapshot
PRIVATE CONAN_PKG::catch2
PRIVATE ${AER_LIBRARIES})
add_test(test_snapshot test_snapshot)
macro(add_test_executable target_name)
add_executable(${target_name} ${ARGN})
set_target_properties(${target_name} PROPERTIES
LINKER_LANGUAGE CXX
CXX_STANDARD 14)
target_include_directories(${target_name}
PRIVATE ${AER_SIMULATOR_CPP_SRC_DIR}
PRIVATE ${AER_SIMULATOR_CPP_EXTERNAL_LIBS})
target_link_libraries(${target_name}
PRIVATE CONAN_PKG::catch2
PRIVATE ${AER_LIBRARIES})
add_test(${target_name} ${target_name})
endmacro()
add_executable(test_snapshot_bdd "src/test_snapshot_bdd.cpp")
set_target_properties(test_snapshot_bdd PROPERTIES
LINKER_LANGUAGE CXX
CXX_STANDARD 14)
target_include_directories(test_snapshot_bdd
PRIVATE ${AER_SIMULATOR_CPP_SRC_DIR}
PRIVATE ${AER_SIMULATOR_CPP_EXTERNAL_LIBS})
target_link_libraries(test_snapshot_bdd
PRIVATE CONAN_PKG::catch2
PRIVATE ${AER_LIBRARIES})
add_test(test_snapshot_bdd test_snapshot_bdd)
add_executable(test_utils "src/test_utils.cpp")
set_target_properties(test_utils PROPERTIES
LINKER_LANGUAGE CXX
CXX_STANDARD 14)
target_include_directories(test_utils
PRIVATE ${AER_SIMULATOR_CPP_SRC_DIR}
PRIVATE ${AER_SIMULATOR_CPP_EXTERNAL_LIBS})
target_link_libraries(test_utils
PRIVATE CONAN_PKG::catch2
PRIVATE ${AER_LIBRARIES})
add_test(test_utils test_utils)
add_test_executable(test_linalg "src/test_linalg.cpp")
# Don't forget to add your test target here
add_custom_target(build_tests
test_snapshot
test_snapshot_bdd
test_utils)
test_linalg)

361
test/src/test_linalg.cpp Normal file
View File

@ -0,0 +1,361 @@
/**
* This code is part of Qiskit.
*
* (C) Copyright IBM 2018, 2019, 2020.
*
* This code is licensed under the Apache License, Version 2.0. You may
* obtain a copy of this license in the LICENSE.txt file in the root directory
* of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
*
* Any modifications or derivative works of this code must retain this
* copyright notice, and modified files need to carry a notice indicating
* that they have been altered from the originals.
*/
#include <framework/stl_ostream.hpp>
#include <map>
#include <type_traits>
#include <framework/linalg/linalg.hpp>
#include <framework/types.hpp>
#include <framework/utils.hpp>
#define CATCH_CONFIG_MAIN
#include <catch2/catch.hpp>
#include "utils.hpp"
using namespace AER::Test::Utilities;
namespace {
// check if polar coordinates are almost equal
// r -> (0,inf)
// angle -> (-PI, PI)
template <typename T>
T eps() {
return std::numeric_limits<T>::epsilon();
}
template<typename T>
bool check_polar_coords(T r, T angle, T r_2, T angle_2, T max_diff = eps<T>(), T max_relative_diff = eps<T>());
template<typename T>
bool check_eigenvector(const std::vector<std::complex<T>> &expected_eigen,
const std::vector<std::complex<T>> &actual_eigen,
T max_diff = eps<T>(),
T max_relative_diff = eps<T>());
// Sometimes eigenvectors differ by a factor and/or phase
// This function compares them taking this into account
template<typename T>
bool check_all_eigenvectors(const matrix<std::complex<T>> &expected,
const matrix<std::complex<T>> &actual,
T max_diff = eps<T>(),
T max_relative_diff = eps<T>());
template<typename T>
using scenarioData = std::tuple<std::string, matrix<std::complex<T>>, matrix<std::complex<T>>, std::vector<T>>;
template<typename T> matrix<std::complex<T>> herm_mat_2d();
template<typename T> matrix<std::complex<T>> herm_mat_2d_eigenvectors();
template<typename T> std::vector<T> herm_mat_2d_eigenvalues();
template<typename T> scenarioData<T> get_herm_2d_scen();
template<typename T> matrix<std::complex<T>> psd_mat_2d();
template<typename T> matrix<std::complex<T>> psd_mat_2d_eigenvectors();
template<typename T> std::vector<T> psd_mat_2d_eigenvalues();
template<typename T> scenarioData<T> get_psd_2d_scen();
template<typename T> matrix<std::complex<T>> psd_mat_2d_with_zero();
template<typename T> matrix<std::complex<T>> psd_mat_2d_wiht_zero_eigenvectors();
template<typename T> std::vector<T> psd_mat_2d_wiht_zero_eigenvalues();
template<typename T> scenarioData<T> get_psd_mat_2d_wiht_zero_scen();
}
TEST_CASE("Basic Matrix Ops", "[matrix]") {
auto mat = matrix<std::complex<double>>(2,2);
mat(0,0) = std::complex<double>(1.0, 0.0);
mat(0,1) = std::complex<double>(1.0, -2.0);
mat(1,0) = std::complex<double>(1.0, 2.0);
mat(1,1) = std::complex<double>(0.5, 0.0);
SECTION("matrix - col_index") {
std::vector<std::complex<double>> col0{ {1.0, 0.0}, {1.0, 2.0} };
std::vector<std::complex<double>> col1{ {1.0, -2.0}, {0.5, 0.0} };
REQUIRE(compare(col0, mat.col_index(0)));
REQUIRE(compare(col1, mat.col_index(1)));
}
SECTION("matrix - col_index") {
std::vector<std::complex<double>> row0{ {1.0, 0.0}, {1.0, -2.0} };
std::vector<std::complex<double>> row1{ {1.0, 2.0}, {0.5, 0.0} };
REQUIRE(compare(row0, mat.row_index(0)));
REQUIRE(compare(row1, mat.row_index(1)));
}
}
TEMPLATE_TEST_CASE("Linear Algebra utilities", "[eigen_hermitian]", float, double) {
std::string scenario_name;
matrix<std::complex<TestType>> herm_mat;
matrix<std::complex<TestType>> expected_eigenvectors;
std::vector<TestType> expected_eigenvalues;
std::tie(scenario_name, herm_mat, expected_eigenvectors, expected_eigenvalues) =
GENERATE(get_herm_2d_scen<TestType>(), get_psd_2d_scen<TestType>(), get_psd_mat_2d_wiht_zero_scen<TestType>());
// We are checking results from a numerical method so we allow for some room in comparisons
TestType eps_threshold = 5 * eps<TestType>();
SECTION(scenario_name + ": zheevx") {
SECTION("sanity check - eigenvals/vecs should recreate original") {
// sanity check
matrix<std::complex<TestType>> sanity_value(herm_mat.GetRows(), herm_mat.GetColumns());
for (size_t j=0; j < expected_eigenvalues.size(); j++) {
sanity_value += expected_eigenvalues[j] * AER::Utils::projector(expected_eigenvectors.col_index(j));
}
REQUIRE(compare(herm_mat, sanity_value, eps_threshold, eps_threshold));
}
SECTION("actual check - heevx returns correctly") {
std::vector<TestType> eigenvalues;
matrix<std::complex<TestType>> eigenvectors;
eigensystem_hermitian(herm_mat, eigenvalues, eigenvectors);
// test equality
REQUIRE(check_all_eigenvectors(expected_eigenvectors, eigenvectors,
eps_threshold, eps_threshold));
REQUIRE(compare(expected_eigenvalues, eigenvalues, eps_threshold, eps_threshold));
// test reconstruction
matrix<std::complex<TestType>> value(herm_mat.GetRows(), herm_mat.GetColumns());
for (size_t j=0; j < eigenvalues.size(); j++) {
value += AER::Utils::projector(eigenvectors.col_index(j)) * eigenvalues[j];
}
REQUIRE(compare(herm_mat, value, eps_threshold, eps_threshold));
}
}
}
TEST_CASE( "Framework Utilities", "[almost_equal]" ) {
SECTION( "The maximum difference between two scalars over 1.0 is greater than epsilon, so they are amlmost equal" ) {
double first = 1.0 + eps<double>();
double actual = 1.0;
// Because the max_diff param is bigger than epsilon, this should be almost equal
REQUIRE(AER::Linalg::almost_equal(first, actual)); //, 1e-15, 1e-15));
}
SECTION( "The difference between two scalars really close to 0 should say are almost equal" ) {
double first = 5e-323; // Really close to the min magnitude of double
double actual = 6e-323;
REQUIRE(AER::Linalg::almost_equal(first, actual)); //, 1e-323, 1e-323));
}
SECTION( "The maximum difference between two complex of doubles over 1.0 is greater than epsilon, so they are almost equal" ) {
std::complex<double> first = { eps<double>() + double(1.0), eps<double>() + double(1.0)};
std::complex<double> actual {1.0, 1.0};
// Because the max_diff param is bigger than epsilon, this should be almost equal
REQUIRE(AER::Linalg::almost_equal(first, actual)); //, 1e-15, 1e-15));
}
SECTION( "The difference between two complex of doubles really close to 0 should say are almost equal" ) {
std::complex<double> first = {5e-323, 5e-323}; // Really close to the min magnitude of double
std::complex<double> actual = {6e-323, 6e-323};
REQUIRE(AER::Linalg::almost_equal(first, actual)); // 1e-323, 1e-323));
}
}
TEST_CASE( "Test_utils", "[check_polar_coords]" ) {
auto r = 1.0;
auto angle = M_PI_2;
auto r_2 = 1.0 + eps<double>();
auto angle_2 = M_PI_2 + eps<double>();
SECTION("Check 2 numbers that are equal"){
REQUIRE(check_polar_coords(r, angle, r_2, angle_2));
}
SECTION("Check 2 numbers that differ in absolute value"){
r_2 = r_2 + 1e3 * eps<double>();
REQUIRE(!check_polar_coords(r, angle, r_2, angle_2));
}
SECTION("Check 2 numbers that differ in absolute value"){
angle_2 = angle_2 + 1e3 * eps<double>();
REQUIRE(!check_polar_coords(r, angle, r_2, angle_2));
}
SECTION("Check corner case: close to +/-0 angles"){
angle = 0.0 - eps<double>()/2.;
angle_2 = -angle;
REQUIRE(check_polar_coords(r, angle, r_2, angle_2));
}
SECTION("Check corner case: angle PI and angle -PI"){
angle = M_PI - eps<double>();
angle_2 = -angle;
REQUIRE(check_polar_coords(r, angle, r_2, angle_2));
}
}
namespace {
template<typename T>
matrix<std::complex<T>> herm_mat_2d() {
auto mat = matrix<std::complex<T>>(2,2);
mat(0,0) = std::complex<T>(1.0,0.0);
mat(0,1) = std::complex<T>(1.0, 2.0);
mat(1,0) = std::complex<T>(1.0, -2.0);
mat(1,1) = std::complex<T>(0.5, 0.0);
return mat;
}
template<typename T>
matrix<std::complex<T>> herm_mat_2d_eigenvectors() {
auto mat = matrix<std::complex<T>>(2,2);
auto den = 3. * std::sqrt(5.);
mat(0,0) = std::complex<T>(-2/den, -4/den);
mat(1,0) = std::complex<T>(5/den, 0.0);
mat(0,1) = std::complex<T>(-1./3.,-2./3.);
mat(1,1) = std::complex<T>(-2./3.,0);
return mat;
}
template<typename T>
std::vector<T> herm_mat_2d_eigenvalues(){
return { -1.5, 3.0 };
}
template <typename T>
scenarioData<T> get_herm_2d_scen(){
return {"Hermitian matrix 2x2", herm_mat_2d<T>(), herm_mat_2d_eigenvectors<T>(), herm_mat_2d_eigenvalues<T>()};
}
template<typename T>
matrix<std::complex<T>> psd_mat_2d(){
auto psd_matrix = matrix<std::complex<T>>(2, 2);
psd_matrix(0,0) = std::complex<T>(13., 0.);
psd_matrix(0,1) = std::complex<T>(0., 5.);
psd_matrix(1,0) = std::complex<T>(0., -5.);
psd_matrix(1,1) = std::complex<T>(2., 0.);
return psd_matrix;
}
template<typename T>
matrix<std::complex<T>> psd_mat_2d_eigenvectors(){
matrix<std::complex<T>> expected_eigenvectors(2, 2);
expected_eigenvectors(0,0) = std::complex<T>(0,-0.3605966767761846214491);
expected_eigenvectors(0,1) = std::complex<T>(0,-0.9327218431547380506075);
expected_eigenvectors(1,0) = std::complex<T>(0.9327218431547380506075,0);
expected_eigenvectors(1,1) = std::complex<T>(-0.3605966767761846214491,0);
return expected_eigenvectors;
}
template<typename T>
std::vector<T> psd_mat_2d_eigenvalues(){
return {0.06696562634074720854471252, 14.93303437365925212532147};
}
template <typename T>
scenarioData<T> get_psd_2d_scen(){
return {"PSD matrix 2x2", psd_mat_2d<T>(), psd_mat_2d_eigenvectors<T>(), psd_mat_2d_eigenvalues<T>()};
}
template<typename T>
matrix<std::complex<T>> psd_mat_2d_with_zero(){
auto psd_matrix = matrix<std::complex<T>>(2, 2);
psd_matrix(0,0) = std::complex<T>(1., 0.);
psd_matrix(0,1) = std::complex<T>(2., 0.);
psd_matrix(1,0) = std::complex<T>(2., 0.);
psd_matrix(1,1) = std::complex<T>(4., 0.);
return psd_matrix;
}
template<typename T>
matrix<std::complex<T>> psd_mat_2d_wiht_zero_eigenvectors(){
matrix<std::complex<T>> expected_eigenvectors(2, 2);
expected_eigenvectors(0,0) = std::complex<T>(-2./std::sqrt(5.),0);
expected_eigenvectors(0,1) = std::complex<T>(1./std::sqrt(5.),0);
expected_eigenvectors(1,0) = std::complex<T>(1./std::sqrt(5.),0);
expected_eigenvectors(1,1) = std::complex<T>(2./std::sqrt(5.),0);
return expected_eigenvectors;
}
template<typename T>
std::vector<T> psd_mat_2d_wiht_zero_eigenvalues(){
return {0.0, 5.0};
}
template <typename T>
scenarioData<T> get_psd_mat_2d_wiht_zero_scen(){
return {"PSD matrix 2x2 with a zero eigen value", psd_mat_2d_with_zero<T>(), psd_mat_2d_wiht_zero_eigenvectors<T>(), psd_mat_2d_wiht_zero_eigenvalues<T>()};
}
template <typename T>
bool check_polar_coords(T r, T angle, T r_2, T angle_2, T max_diff, T max_relative_diff){
if(!AER::Linalg::almost_equal(r, r_2, max_diff, max_relative_diff)) return false;
if(!AER::Linalg::almost_equal(angle, angle_2, max_diff, max_relative_diff)){
// May be corner case with PI and -PI
T angle_plus = angle > 0. ? angle : angle + 2*M_PI;
T angle_2_plus = angle_2 > 0. ? angle_2 : angle_2 + 2 * M_PI;
if(!AER::Linalg::almost_equal(angle_plus, angle_2_plus, max_diff, max_relative_diff)) return false;
}
return true;
}
template <typename T>
bool check_eigenvector(const std::vector<std::complex<T>>& expected_eigen,
const std::vector<std::complex<T>>& actual_eigen,
T max_diff,
T max_relative_diff){
auto div = expected_eigen[0] / actual_eigen[0];
T r = std::abs(div);
T angle = std::arg(div);
for (int j = 1; j < expected_eigen.size(); j++) {
auto div_2 = expected_eigen[j] / actual_eigen[j];
T r_2 = std::abs(div_2);
T angle_2 = std::arg(div_2);
// Check that factor is consistent across components
if (!check_polar_coords(r, angle, r_2, angle_2, max_diff, max_relative_diff)) {
return false;
}
}
return true;
}
template <typename T>
bool check_all_eigenvectors(const matrix<std::complex<T>>& expected,
const matrix<std::complex<T>>& actual,
T max_diff,
T max_relative_diff) {
auto col_num = expected.GetColumns();
if (expected.size() != actual.size() || expected.GetColumns() != expected.GetColumns()) {
return false;
}
for (int i = 0; i < col_num; i++) {
auto expected_eigen = expected.col_index(i);
auto actual_eigen = actual.col_index(i);
if(!check_eigenvector(expected_eigen, actual_eigen, max_diff, max_relative_diff)){
std::cout << "Expected: " << std::setprecision(16) << expected << std::endl;
std::cout << "Actual: " << actual << std::endl;
return false;
}
}
return true;
}
}

35
test/src/test_utils.cpp
View File

@ -1,35 +0,0 @@
#define CATCH_CONFIG_MAIN
#include <map>
#include <catch2/catch.hpp>
#include <cmath>
#include <limits>
#include "framework/linalg/almost_equal.hpp"
#include "utils.hpp"
namespace AER{
namespace Test{
TEST_CASE( "Framework Utilities", "[almost_equal]" ) {
SECTION( "The maximum difference between two numbers over 1.0 is greater than epsilon, so they are amlmost equal" ) {
double first = 1.0 + std::numeric_limits<double>::epsilon();
double actual = 1.0;
// Because the max_diff param is bigger than epsilon, this should be almost equal
REQUIRE(Linalg::almost_equal<decltype(first)>(first, actual, 1e-15, 1e-15));
}
SECTION( "The difference between two numbers really close to 0 should say are almost equal" ) {
double first = 5e-323; // Really close to the min magnitude of double
double actual = 6e-323;
REQUIRE(Linalg::almost_equal<decltype(first)>(first, actual, 1e-323, 1e-323));
}
}
//------------------------------------------------------------------------------
} // end namespace Test
//------------------------------------------------------------------------------
} // end namespace AER
//------------------------------------------------------------------------------

102
test/src/utils.hpp
View File

@ -1,6 +1,39 @@
/**
* This code is part of Qiskit.
*
* (C) Copyright IBM 2018, 2019, 2020.
*
* This code is licensed under the Apache License, Version 2.0. You may
* obtain a copy of this license in the LICENSE.txt file in the root directory
* of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
*
* Any modifications or derivative works of this code must retain this
* copyright notice, and modified files need to carry a notice indicating
* that they have been altered from the originals.
*/
#include <string>
#include "framework/json.hpp"
#include <catch2/catch.hpp>
namespace Catch {
template <typename T>
std::string convertMyTypeToString(const std::vector<T> &value) {
std::stringstream oss;
oss << value;
return oss.str();
}
template <typename T> struct StringMaker<std::vector<T>> {
static std::string convert(const std::vector<AER::complex_t> &value) {
return convertMyTypeToString(value);
}
};
} // namespace Catch
namespace AER{
namespace Test{
namespace Utilities {
@ -15,6 +48,75 @@ namespace Utilities {
return start;
}
template<typename T, typename U>
bool _compare(const matrix<T> &lhs, const matrix<T> &rhs, U max_diff = std::numeric_limits<U>::epsilon(),
U max_relative_diff = std::numeric_limits<U>::epsilon()) {
bool res = true;
std::ostringstream message;
if (lhs.size() != rhs.size()) {
res = false;
}
for (size_t i = 0; i < lhs.GetRows(); ++i) {
for (size_t j = 0; j < lhs.GetColumns(); ++j) {
if (!(AER::Linalg::almost_equal(lhs(i, j), rhs(i, j), max_diff, max_relative_diff))) {
message << "Matrices differ at element: (" << i << ", " << j << ")"
<< std::setprecision(22) << ". [" << lhs(i, j)
<< "] != [" << rhs(i, j) << "]\n";
res = false;
}
}
}
if (!res) {
message << "Matrices differ: " << lhs
<< " != " << rhs << std::endl;
std::cout << message.str();
}
return res;
}
template<typename T>
bool compare(const matrix<T> &lhs, const matrix<T> &rhs, T max_diff = std::numeric_limits<T>::epsilon(),
T max_relative_diff = std::numeric_limits<T>::epsilon()) {
return _compare(lhs, rhs, max_diff, max_relative_diff);
}
template<typename T>
bool compare(const matrix<std::complex<T>> &lhs, const matrix<std::complex<T>> &rhs,
T max_diff = std::numeric_limits<T>::epsilon(),
T max_relative_diff = std::numeric_limits<T>::epsilon()) {
return _compare(lhs, rhs, max_diff, max_relative_diff);
}
template<typename T, typename U>
bool _compare(const std::vector<T> &lhs, const std::vector<T> &rhs, U max_diff = std::numeric_limits<U>::epsilon(),
U max_relative_diff = std::numeric_limits<U>::epsilon()) {
if (lhs.size() != rhs.size())
return false;
for (size_t i = 0; i < lhs.size(); ++i) {
if (!(AER::Linalg::almost_equal(lhs[i], rhs[i], max_diff, max_relative_diff))) {
std::cout << "Vectors differ at element: " << i << std::setprecision(22) << ". [" << lhs[i]
<< "] != [" << rhs[i] << "]\n";
std::cout << "Vectors differ: " << Catch::convertMyTypeToString(lhs)
<< " != " << Catch::convertMyTypeToString(rhs) << std::endl;
return false;
}
}
return true;
}
template<typename T>
bool compare(const std::vector<T> &lhs, const std::vector<T> &rhs, T max_diff = std::numeric_limits<T>::epsilon(),
T max_relative_diff = std::numeric_limits<T>::epsilon()) {
return _compare(lhs, rhs, max_diff, max_relative_diff);
}
template<typename T>
bool compare(const std::vector<std::complex<T>> &lhs, const std::vector<std::complex<T>> &rhs,
T max_diff = std::numeric_limits<T>::epsilon(),
T max_relative_diff = std::numeric_limits<T>::epsilon()) {
return _compare(lhs, rhs, max_diff, max_relative_diff);
}
} // End of Utilities namespace
} // End of Test namespace
} // End of AER namspace