From 86a27e3faafd89776387deda14992553782d65ab Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Daniel=20Rodr=C3=ADguez=20L?= Date: Wed, 10 Jan 2024 01:56:12 +0100 Subject: [PATCH] Add support for BLAS SVD functions in MPS simulation (#1897) * add lapack svd method + test * forgot to unset env variable in test * fix seg fault bc arrays were too big * code style and releasenote for PR * address microsoft C2131? * style + another C2131 * missing free * change test to support windows python3.8 * Update releasenotes/notes/compute-svd-with-lapack-3ee992d371d653d1.yaml Co-authored-by: merav-aharoni * remove unnecessary comments * undo ifdef DEBUG * automatic selector for QR or D&C in LAPACK SVD * codestyle, enable MPS lapack using run_options --------- Co-authored-by: merav-aharoni Co-authored-by: Jun Doi --- qiskit_aer/backends/aer_compiler.py | 1 + qiskit_aer/backends/aer_simulator.py | 4 + qiskit_aer/backends/qasm_simulator.py | 1 + .../wrappers/aer_controller_binding.hpp | 3 + ...pute-svd-with-lapack-3ee992d371d653d1.yaml | 8 ++ src/framework/config.hpp | 6 +- src/framework/lapack_protos.hpp | 37 ++++++ .../matrix_product_state.hpp | 6 +- .../matrix_product_state_internal.cpp | 24 +++- .../matrix_product_state_internal.hpp | 5 + .../matrix_product_state_tensor.hpp | 18 ++- src/simulators/matrix_product_state/svd.cpp | 121 +++++++++++++++++- src/simulators/matrix_product_state/svd.hpp | 19 ++- .../backends/aer_simulator/test_options.py | 35 +++++ 14 files changed, 269 insertions(+), 19 deletions(-) create mode 100644 releasenotes/notes/compute-svd-with-lapack-3ee992d371d653d1.yaml create mode 100644 src/framework/lapack_protos.hpp diff --git a/qiskit_aer/backends/aer_compiler.py b/qiskit_aer/backends/aer_compiler.py index c3f7067738..bbb2e18a25 100644 --- a/qiskit_aer/backends/aer_compiler.py +++ b/qiskit_aer/backends/aer_compiler.py @@ -537,6 +537,7 @@ def compile_circuit(circuits, basis_gates=None, optypes=None): "chop_threshold": (float, np.floating), "mps_parallel_threshold": (int, np.integer), "mps_omp_threads": (int, np.integer), + "mps_lapack": (bool, np.bool_), "tensor_network_num_sampling_qubits": (int, np.integer), "use_cuTensorNet_autotuning": (bool, np.bool_), "parameterizations": (list), diff --git a/qiskit_aer/backends/aer_simulator.py b/qiskit_aer/backends/aer_simulator.py index 8ed40b170e..9638bd8a71 100644 --- a/qiskit_aer/backends/aer_simulator.py +++ b/qiskit_aer/backends/aer_simulator.py @@ -448,6 +448,9 @@ class AerSimulator(AerBackend): * ``mps_omp_threads`` (int): This option sets the number of OMP threads (Default: 1). + * ``mps_lapack`` (bool): This option indicates to compute the SVD function + using OpenBLAS/Lapack interface (Default: False). + These backend options only apply when using the ``tensor_network`` simulation method: @@ -785,6 +788,7 @@ def _default_options(cls): chop_threshold=1e-8, mps_parallel_threshold=14, mps_omp_threads=1, + mps_lapack=False, # tensor network options tensor_network_num_sampling_qubits=10, use_cuTensorNet_autotuning=False, diff --git a/qiskit_aer/backends/qasm_simulator.py b/qiskit_aer/backends/qasm_simulator.py index 06288afe95..a36731403a 100644 --- a/qiskit_aer/backends/qasm_simulator.py +++ b/qiskit_aer/backends/qasm_simulator.py @@ -527,6 +527,7 @@ def _default_options(cls): chop_threshold=1e-8, mps_parallel_threshold=14, mps_omp_threads=1, + mps_lapack=False, ) @classmethod diff --git a/qiskit_aer/backends/wrappers/aer_controller_binding.hpp b/qiskit_aer/backends/wrappers/aer_controller_binding.hpp index 9c15f08650..ea3b48dfae 100644 --- a/qiskit_aer/backends/wrappers/aer_controller_binding.hpp +++ b/qiskit_aer/backends/wrappers/aer_controller_binding.hpp @@ -242,6 +242,7 @@ void bind_aer_controller(MODULE m) { aer_config.def_readwrite("mps_parallel_threshold", &Config::mps_parallel_threshold); aer_config.def_readwrite("mps_omp_threads", &Config::mps_omp_threads); + aer_config.def_readwrite("mps_lapack", &Config::mps_lapack); // # tensor network options aer_config.def_readwrite("tensor_network_num_sampling_qubits", &Config::tensor_network_num_sampling_qubits); @@ -477,6 +478,7 @@ void bind_aer_controller(MODULE m) { write_value(30, config.chop_threshold), write_value(41, config.mps_parallel_threshold), write_value(42, config.mps_omp_threads), + write_value(101, config.mps_lapack), write_value(43, config.tensor_network_num_sampling_qubits), write_value(44, config.use_cuTensorNet_autotuning), write_value(45, config.library_dir), @@ -571,6 +573,7 @@ void bind_aer_controller(MODULE m) { read_value(t, 30, config.chop_threshold); read_value(t, 41, config.mps_parallel_threshold); read_value(t, 42, config.mps_omp_threads); + read_value(t, 101, config.mps_lapack); read_value(t, 43, config.tensor_network_num_sampling_qubits); read_value(t, 44, config.use_cuTensorNet_autotuning); read_value(t, 45, config.library_dir); diff --git a/releasenotes/notes/compute-svd-with-lapack-3ee992d371d653d1.yaml b/releasenotes/notes/compute-svd-with-lapack-3ee992d371d653d1.yaml new file mode 100644 index 0000000000..f861b6fbab --- /dev/null +++ b/releasenotes/notes/compute-svd-with-lapack-3ee992d371d653d1.yaml @@ -0,0 +1,8 @@ +--- +features: + - | + Replace Qiskit SVD function with OpenBLAS/LAPACK SVD functions ``zgesvd`` + and ``zgesdd``. By default ``zgesvd`` is used. Performance of ``zgesdd`` is better than + that of ``zgesvd`` on large matrices, whereas ``zgesvd`` performs better on small matrices. + User can use ``zgesdd`` function + setting the environment variable ``QISKIT_LAPACK_SVD=DC``. diff --git a/src/framework/config.hpp b/src/framework/config.hpp index 1074f7acdf..ae95508ef2 100644 --- a/src/framework/config.hpp +++ b/src/framework/config.hpp @@ -125,6 +125,7 @@ struct Config { double chop_threshold = 1e-8; uint_t mps_parallel_threshold = 14; uint_t mps_omp_threads = 1; + bool mps_lapack = false; // # tensor network options uint_t tensor_network_num_sampling_qubits = 10; bool use_cuTensorNet_autotuning = false; @@ -231,6 +232,7 @@ struct Config { chop_threshold = 1e-8; mps_parallel_threshold = 14; mps_omp_threads = 1; + mps_lapack = false; // # tensor network options tensor_network_num_sampling_qubits = 10; use_cuTensorNet_autotuning = false; @@ -359,6 +361,7 @@ struct Config { chop_threshold = other.chop_threshold; mps_parallel_threshold = other.mps_parallel_threshold; mps_omp_threads = other.mps_omp_threads; + mps_lapack = other.mps_lapack; // # tensor network options tensor_network_num_sampling_qubits = other.tensor_network_num_sampling_qubits; @@ -499,6 +502,7 @@ inline void from_json(const json_t &js, Config &config) { get_value(config.chop_threshold, "chop_threshold", js); get_value(config.mps_parallel_threshold, "mps_parallel_threshold", js); get_value(config.mps_omp_threads, "mps_omp_threads", js); + get_value(config.mps_lapack, "mps_lapack", js); // # tensor network options get_value(config.tensor_network_num_sampling_qubits, "tensor_network_num_sampling_qubits", js); @@ -542,4 +546,4 @@ inline void from_json(const json_t &js, Config &config) { } // namespace AER -#endif \ No newline at end of file +#endif diff --git a/src/framework/lapack_protos.hpp b/src/framework/lapack_protos.hpp new file mode 100644 index 0000000000..18967f6496 --- /dev/null +++ b/src/framework/lapack_protos.hpp @@ -0,0 +1,37 @@ +// Dependencies: BLAS - LAPACK + +#ifndef _aer_framework_lapack_protos_hpp +#define _aer_framework_lapack_protos_hpp + +#include +#include +#include +#include + +#ifdef __cplusplus +extern "C" { +#endif + +// LAPACK SVD function +// https://netlib.org/lapack/explore-html/d3/da8/group__complex16_g_esing_gad6f0c85f3cca2968e1ef901d2b6014ee.html +void zgesvd_(const char *jobu, const char *jobvt, const size_t *m, + const size_t *n, std::complex *a, const size_t *lda, + double *s, std::complex *u, const size_t *ldu, + std::complex *vt, const size_t *ldvt, + std::complex *work, const size_t *lwork, double *rwork, + int *info); + +// D&C approach +// https://netlib.org/lapack/explore-html/d3/da8/group__complex16_g_esing_gaccb06ed106ce18814ad7069dcb43aa27.html +void zgesdd_(const char *jobz, const size_t *m, const size_t *n, + std::complex *a, const size_t *lda, double *s, + std::complex *u, const size_t *ldu, + std::complex *vt, const size_t *ldvt, + std::complex *work, const size_t *lwork, double *rwork, + int *iwork, int *info); + +#ifdef __cplusplus +} +#endif + +#endif // end __lapack_protos_h_ diff --git a/src/simulators/matrix_product_state/matrix_product_state.hpp b/src/simulators/matrix_product_state/matrix_product_state.hpp index b1ae10c90f..b323d95a51 100644 --- a/src/simulators/matrix_product_state/matrix_product_state.hpp +++ b/src/simulators/matrix_product_state/matrix_product_state.hpp @@ -369,6 +369,9 @@ void State::set_config(const Config &config) { MPS::set_mps_swap_direction(MPS_swap_direction::SWAP_RIGHT); else MPS::set_mps_swap_direction(MPS_swap_direction::SWAP_LEFT); + + // Set LAPACK SVD + MPS::set_mps_lapack_svd(config.mps_lapack); } void State::add_metadata(ExperimentResult &result) const { @@ -380,6 +383,7 @@ void State::add_metadata(ExperimentResult &result) const { "matrix_product_state_sample_measure_algorithm"); if (MPS::get_mps_log_data()) result.metadata.add("{" + MPS::output_log() + "}", "MPS_log_data"); + result.metadata.add(MPS::get_mps_lapack_svd(), "matrix_product_state_lapack"); } void State::output_bond_dimensions(const Operations::Op &op) const { @@ -828,4 +832,4 @@ std::pair State::sample_measure_with_prob(const reg_t &qubits, //------------------------------------------------------------------------- } // end namespace AER //------------------------------------------------------------------------- -#endif \ No newline at end of file +#endif diff --git a/src/simulators/matrix_product_state/matrix_product_state_internal.cpp b/src/simulators/matrix_product_state/matrix_product_state_internal.cpp index b98e793ab2..e190577c7c 100644 --- a/src/simulators/matrix_product_state/matrix_product_state_internal.cpp +++ b/src/simulators/matrix_product_state/matrix_product_state_internal.cpp @@ -44,6 +44,7 @@ enum MPS_swap_direction MPS::mps_swap_direction_ = double MPS::json_chop_threshold_ = 1E-8; std::stringstream MPS::logging_str_; bool MPS::mps_log_data_ = 0; +bool MPS::mps_lapack_ = false; //------------------------------------------------------------------------ // local function declarations @@ -662,8 +663,9 @@ void MPS::common_apply_2_qubit_gate( MPS_Tensor left_gamma, right_gamma; rvector_t lambda; - double discarded_value = - MPS_Tensor::Decompose(temp, left_gamma, lambda, right_gamma); + double discarded_value = MPS_Tensor::Decompose(temp, left_gamma, lambda, + right_gamma, MPS::mps_lapack_); + if (discarded_value > json_chop_threshold_) MPS::print_to_log("discarded_value=", discarded_value, ", "); @@ -1786,16 +1788,21 @@ void MPS::initialize_from_matrix(uint_t num_qubits, const cmatrix_t &mat) { if (first_iter) { remaining_matrix = mat; } else { - cmatrix_t temp = mul_matrix_by_lambda(V, S); + cmatrix_t temp; + if (MPS::mps_lapack_) { // When using Lapack, V is V dagger + temp = mul_matrix_by_lambda(AER::Utils::dagger(V), S); + } else { + temp = mul_matrix_by_lambda(V, S); + } remaining_matrix = AER::Utils::dagger(temp); } reshaped_matrix = reshape_matrix(remaining_matrix); // step 2 - SVD S.clear(); S.resize(std::min(reshaped_matrix.GetRows(), reshaped_matrix.GetColumns())); - csvd_wrapper(reshaped_matrix, U, S, V); + csvd_wrapper(reshaped_matrix, U, S, V, MPS::mps_lapack_); reduce_zeros(U, S, V, MPS_Tensor::get_max_bond_dimension(), - MPS_Tensor::get_truncation_threshold()); + MPS_Tensor::get_truncation_threshold(), MPS::mps_lapack_); // step 3 - update q_reg_ with new gamma and new lambda // increment number of qubits in the MPS structure @@ -1811,7 +1818,12 @@ void MPS::initialize_from_matrix(uint_t num_qubits, const cmatrix_t &mat) { first_iter = false; } // step 4 - create the rightmost gamma and update q_reg_ - std::vector right_data = reshape_V_after_SVD(V); + std::vector right_data; + if (MPS::mps_lapack_) { + right_data = reshape_VH_after_SVD(V); + } else { + right_data = reshape_V_after_SVD(V); + } MPS_Tensor right_gamma(right_data[0], right_data[1]); q_reg_.push_back(right_gamma); diff --git a/src/simulators/matrix_product_state/matrix_product_state_internal.hpp b/src/simulators/matrix_product_state/matrix_product_state_internal.hpp index 19c898fa8a..e6a82c28cc 100644 --- a/src/simulators/matrix_product_state/matrix_product_state_internal.hpp +++ b/src/simulators/matrix_product_state/matrix_product_state_internal.hpp @@ -319,6 +319,8 @@ class MPS { mps_swap_direction_ = direction; } + static void set_mps_lapack_svd(bool mps_lapack) { mps_lapack_ = mps_lapack; } + static uint_t get_omp_threads() { return omp_threads_; } static uint_t get_omp_threshold() { return omp_threshold_; } static double get_json_chop_threshold() { return json_chop_threshold_; } @@ -330,6 +332,8 @@ class MPS { static bool get_mps_log_data() { return mps_log_data_; } + static bool get_mps_lapack_svd() { return mps_lapack_; } + static MPS_swap_direction get_swap_direction() { return mps_swap_direction_; } //---------------------------------------------------------------- @@ -564,6 +568,7 @@ class MPS { static std::stringstream logging_str_; static bool mps_log_data_; static MPS_swap_direction mps_swap_direction_; + static bool mps_lapack_; }; inline std::ostream &operator<<(std::ostream &out, const rvector_t &vec) { diff --git a/src/simulators/matrix_product_state/matrix_product_state_tensor.hpp b/src/simulators/matrix_product_state/matrix_product_state_tensor.hpp index 0155925843..b769e8a59f 100644 --- a/src/simulators/matrix_product_state/matrix_product_state_tensor.hpp +++ b/src/simulators/matrix_product_state/matrix_product_state_tensor.hpp @@ -158,7 +158,8 @@ class MPS_Tensor { const rvector_t &lambda, const MPS_Tensor &right_gamma, bool mul_by_lambda); static double Decompose(MPS_Tensor &temp, MPS_Tensor &left_gamma, - rvector_t &lambda, MPS_Tensor &right_gamma); + rvector_t &lambda, MPS_Tensor &right_gamma, + bool mps_lapack); static void reshape_for_3_qubits_before_SVD(const std::vector data, MPS_Tensor &reshaped_tensor); static void contract_2_dimensions(const MPS_Tensor &left_gamma, @@ -590,20 +591,25 @@ void MPS_Tensor::contract_2_dimensions(const MPS_Tensor &left_gamma, // Returns: none. //--------------------------------------------------------------- double MPS_Tensor::Decompose(MPS_Tensor &temp, MPS_Tensor &left_gamma, - rvector_t &lambda, MPS_Tensor &right_gamma) { + rvector_t &lambda, MPS_Tensor &right_gamma, + bool mps_lapack) { cmatrix_t C; C = reshape_before_SVD(temp.data_); cmatrix_t U, V; rvector_t S(std::min(C.GetRows(), C.GetColumns())); - csvd_wrapper(C, U, S, V); + csvd_wrapper(C, U, S, V, mps_lapack); double discarded_value = 0.0; - discarded_value = - reduce_zeros(U, S, V, max_bond_dimension_, truncation_threshold_); + discarded_value = reduce_zeros(U, S, V, max_bond_dimension_, + truncation_threshold_, mps_lapack); left_gamma.data_ = reshape_U_after_SVD(U); lambda = S; - right_gamma.data_ = reshape_V_after_SVD(V); + if (mps_lapack) { // When using Lapack V is V dagger + right_gamma.data_ = reshape_VH_after_SVD(V); + } else { + right_gamma.data_ = reshape_V_after_SVD(V); + } return discarded_value; } diff --git a/src/simulators/matrix_product_state/svd.cpp b/src/simulators/matrix_product_state/svd.cpp index 9a6624e60f..150d97af4f 100644 --- a/src/simulators/matrix_product_state/svd.cpp +++ b/src/simulators/matrix_product_state/svd.cpp @@ -65,6 +65,11 @@ std::vector reshape_V_after_SVD(const cmatrix_t V) { AER::Utils::split(AER::Utils::dagger(V), Res[0], Res[1], 1); return Res; } +std::vector reshape_VH_after_SVD(const cmatrix_t V) { + std::vector Res(2); + AER::Utils::split(V, Res[0], Res[1], 1); + return Res; +} //------------------------------------------------------------- // function name: num_of_SV @@ -85,7 +90,8 @@ uint_t num_of_SV(rvector_t S, double threshold) { } double reduce_zeros(cmatrix_t &U, rvector_t &S, cmatrix_t &V, - uint_t max_bond_dimension, double truncation_threshold) { + uint_t max_bond_dimension, double truncation_threshold, + bool mps_lapack) { uint_t SV_num = num_of_SV(S, CHOP_THRESHOLD); uint_t new_SV_num = SV_num; @@ -107,7 +113,12 @@ double reduce_zeros(cmatrix_t &U, rvector_t &S, cmatrix_t &V, } U.resize(U.GetRows(), new_SV_num); S.resize(new_SV_num); - V.resize(V.GetRows(), new_SV_num); + // When using LAPACK function, V is V dagger + if (mps_lapack) { + V.resize(new_SV_num, V.GetColumns()); + } else { + V.resize(V.GetRows(), new_SV_num); + } // discarded_value is the sum of the squares of the Schmidt coeffients // that were discarded by approximation @@ -130,9 +141,28 @@ double reduce_zeros(cmatrix_t &U, rvector_t &S, cmatrix_t &V, return discarded_value; } +void validate_SVdD_result(const cmatrix_t &A, const cmatrix_t &U, + const rvector_t &S, const cmatrix_t &V) { + const uint_t nrows = A.GetRows(), ncols = A.GetColumns(); + + cmatrix_t diag_S = diag(S, nrows, ncols); + cmatrix_t product = U * diag_S; + product = product * V; + + for (uint_t ii = 0; ii < nrows; ii++) + for (uint_t jj = 0; jj < ncols; jj++) + if (!Linalg::almost_equal(std::abs(A(ii, jj)), std::abs(product(ii, jj)), + THRESHOLD)) { + std::cout << std::abs(A(ii, jj)) << " vs " << std::abs(product(ii, jj)) + << std::endl; + throw std::runtime_error("Error: Wrong SVD calculations: A != USV*"); + } +} + void validate_SVD_result(const cmatrix_t &A, const cmatrix_t &U, const rvector_t &S, const cmatrix_t &V) { const uint_t nrows = A.GetRows(), ncols = A.GetColumns(); + cmatrix_t diag_S = diag(S, nrows, ncols); cmatrix_t product = U * diag_S; product = product * AER::Utils::dagger(V); @@ -519,7 +549,17 @@ status csvd(cmatrix_t &A, cmatrix_t &U, rvector_t &S, cmatrix_t &V) { return SUCCESS; } -void csvd_wrapper(cmatrix_t &A, cmatrix_t &U, rvector_t &S, cmatrix_t &V) { +void csvd_wrapper(cmatrix_t &A, cmatrix_t &U, rvector_t &S, cmatrix_t &V, + bool lapack) { + if (lapack) { + lapack_csvd_wrapper(A, U, S, V); + } else { + qiskit_csvd_wrapper(A, U, S, V); + } +} + +void qiskit_csvd_wrapper(cmatrix_t &A, cmatrix_t &U, rvector_t &S, + cmatrix_t &V) { cmatrix_t copied_A = A; int times = 0; #ifdef DEBUG @@ -552,4 +592,79 @@ void csvd_wrapper(cmatrix_t &A, cmatrix_t &U, rvector_t &S, cmatrix_t &V) { S[k] /= pow(mul_factor, times); } +void lapack_csvd_wrapper(cmatrix_t &A, cmatrix_t &U, rvector_t &S, + cmatrix_t &V) { + // Activated by default as requested in the PR + // #ifdef DEBUG + cmatrix_t tempA = A; + // #endif + + const size_t m = A.GetRows(), n = A.GetColumns(); + const size_t min_dim = std::min(m, n); + const size_t lda = std::max(m, n); + size_t lwork = 2 * min_dim + lda; + + U.resize(m, m); + V.resize(n, n); + + complex_t *lapackA = A.move_to_buffer(), *lapackU = U.move_to_buffer(), + *lapackV = V.move_to_buffer(); + + double *lapackS = new double[min_dim]; + complex_t *work = new complex_t[lwork]; + int info; + + if (m >= 64 && n >= 64) { + // From experimental results, matrices equal or bigger than this size + // perform better using Divide and Conquer approach + int *iwork = new int[8 * min_dim]; + int rwork_size = std::max(5 * min_dim * min_dim + 5 * min_dim, + 2 * m * n + 2 * min_dim * min_dim + min_dim); + + double *rwork = (double *)calloc(rwork_size, sizeof(double)); + lwork = -1; + zgesdd_("A", &m, &n, lapackA, &m, lapackS, lapackU, &m, lapackV, &n, work, + &lwork, rwork, iwork, &info); + + lwork = (int)work[0].real(); + complex_t *work_ = (complex_t *)calloc(lwork, sizeof(complex_t)); + + zgesdd_("A", &m, &n, lapackA, &m, lapackS, lapackU, &m, lapackV, &n, work_, + &lwork, rwork, iwork, &info); + + delete iwork; + free(rwork); + free(work_); + } else { + // Default execution follows original method + double *rwork = (double *)calloc(5 * min_dim, sizeof(double)); + zgesvd_("A", "A", &m, &n, lapackA, &m, lapackS, lapackU, &m, lapackV, &n, + work, &lwork, rwork, &info); + free(rwork); + } + A = cmatrix_t::move_from_buffer(m, n, lapackA); + U = cmatrix_t::move_from_buffer(m, m, lapackU); + V = cmatrix_t::move_from_buffer(n, n, lapackV); + + S.clear(); + for (int i = 0; i < min_dim; i++) + S.push_back(lapackS[i]); + + // Activated by default as requested in the PR + // #ifdef DEBUG + validate_SVdD_result(tempA, U, S, V); + // #endif + + delete lapackS; + delete work; + + if (info == 0) { + return; + } else { + std::stringstream ss; + ss << " SVD failed"; + throw std::runtime_error(ss.str()); + } +} + } // namespace AER diff --git a/src/simulators/matrix_product_state/svd.hpp b/src/simulators/matrix_product_state/svd.hpp index 2a999e40d4..fac77797fb 100644 --- a/src/simulators/matrix_product_state/svd.hpp +++ b/src/simulators/matrix_product_state/svd.hpp @@ -15,8 +15,10 @@ #ifndef SVD_HPP_ #define SVD_HPP_ +#include "framework/lapack_protos.hpp" #include "framework/types.hpp" #include "framework/utils.hpp" + #include #include @@ -32,13 +34,26 @@ cmatrix_t reshape_before_SVD(std::vector data); std::vector reshape_U_after_SVD(cmatrix_t U); rvector_t reshape_S_after_SVD(rvector_t S); std::vector reshape_V_after_SVD(const cmatrix_t V); +std::vector reshape_VH_after_SVD(const cmatrix_t V); uint_t num_of_SV(rvector_t S, double threshold); double reduce_zeros(cmatrix_t &U, rvector_t &S, cmatrix_t &V, - uint_t max_bond_dimension, double truncation_threshold); + uint_t max_bond_dimension, double truncation_threshold, + bool mps_lapack); status csvd(cmatrix_t &C, cmatrix_t &U, rvector_t &S, cmatrix_t &V); -void csvd_wrapper(cmatrix_t &C, cmatrix_t &U, rvector_t &S, cmatrix_t &V); +// Entry point for the SVD calculation +void csvd_wrapper(cmatrix_t &C, cmatrix_t &U, rvector_t &S, cmatrix_t &V, + bool lapack); +// Original qiskit call +void qiskit_csvd_wrapper(cmatrix_t &C, cmatrix_t &U, rvector_t &S, + cmatrix_t &V); +// Lapack call +void lapack_csvd_wrapper(cmatrix_t &C, cmatrix_t &U, rvector_t &S, + cmatrix_t &V); + void validate_SVD_result(const cmatrix_t &A, const cmatrix_t &U, const rvector_t &S, const cmatrix_t &V); +void validate_SVdD_result(const cmatrix_t &A, const cmatrix_t &U, + const rvector_t &S, const cmatrix_t &V); //------------------------------------------------------------------------- } // end namespace AER diff --git a/test/terra/backends/aer_simulator/test_options.py b/test/terra/backends/aer_simulator/test_options.py index e670177866..ca63207181 100644 --- a/test/terra/backends/aer_simulator/test_options.py +++ b/test/terra/backends/aer_simulator/test_options.py @@ -302,3 +302,38 @@ def test_num_qubits(self, method): num_qubits = FakeMontreal().configuration().num_qubits backend = AerSimulator.from_backend(FakeMontreal(), method=method) self.assertGreaterEqual(backend.configuration().num_qubits, num_qubits) + + def test_mps_svd_method(self): + """Test env. variabe to change MPS SVD method""" + # based on test_mps_options test + import os + + shots = 4000 + method = "matrix_product_state" + backend_og = self.backend(method=method) + backend_lapack = self.backend(method=method) + + n = 10 + circuit = QuantumCircuit(n) + + for times in range(2): + for i in range(0, n, 2): + circuit.unitary(random_unitary(4), [i, i + 1]) + for i in range(1, n - 1): + circuit.cx(0, i) + circuit.save_statevector("sv") + + result_og = backend_og.run(circuit, shots=shots).result() + original_sv = result_og.data(0)["sv"] + + # run with lapack svd method + result_lapack = backend_lapack.run(circuit, shots=shots, mps_lapack=True).result() + lapack_sv = result_lapack.data(0)["sv"] + + # result_lapack should have the metadata indicating that it used LAPACK + # for the SVD + self.assertTrue("matrix_product_state_lapack" in result_lapack._get_experiment().metadata) + self.assertTrue(result_lapack._get_experiment().metadata["matrix_product_state_lapack"]) + + # should give the same state vector + self.assertAlmostEqual(state_fidelity(original_sv, lapack_sv), 1.0)