From 8eac5a9e2aa08bd3fe9b0ca5b9a41b6a23e4b5c7 Mon Sep 17 00:00:00 2001 From: patataman Date: Fri, 19 May 2023 15:12:38 +0900 Subject: [PATCH 01/13] add lapack svd method + test --- src/framework/lapack_protos.hpp | 36 ++++++++ .../matrix_product_state_internal.cpp | 18 +++- .../matrix_product_state_tensor.hpp | 6 +- src/simulators/matrix_product_state/svd.cpp | 89 ++++++++++++++++++- src/simulators/matrix_product_state/svd.hpp | 10 +++ .../backends/aer_simulator/test_options.py | 36 ++++++++ 6 files changed, 191 insertions(+), 4 deletions(-) create mode 100644 src/framework/lapack_protos.hpp diff --git a/src/framework/lapack_protos.hpp b/src/framework/lapack_protos.hpp new file mode 100644 index 0000000000..6ca64047f0 --- /dev/null +++ b/src/framework/lapack_protos.hpp @@ -0,0 +1,36 @@ +// 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_internal.cpp b/src/simulators/matrix_product_state/matrix_product_state_internal.cpp index b98e793ab2..d34ed843cc 100644 --- a/src/simulators/matrix_product_state/matrix_product_state_internal.cpp +++ b/src/simulators/matrix_product_state/matrix_product_state_internal.cpp @@ -542,7 +542,9 @@ void MPS::apply_swap_internal(uint_t index_A, uint_t index_B, bool swap_gate) { // to the left std::swap(qubit_ordering_.order_[index_A], qubit_ordering_.order_[index_B]); // For logging purposes: +#ifdef DEBUG print_to_log_internal_swap(index_A, index_B); +#endif // update qubit locations after all the swaps for (uint_t i = 0; i < num_qubits_; i++) @@ -664,8 +666,10 @@ void MPS::common_apply_2_qubit_gate( rvector_t lambda; double discarded_value = MPS_Tensor::Decompose(temp, left_gamma, lambda, right_gamma); +#ifdef DEBUG if (discarded_value > json_chop_threshold_) MPS::print_to_log("discarded_value=", discarded_value, ", "); +#endif if (A != 0) left_gamma.div_Gamma_by_left_Lambda(lambda_reg_[A - 1]); @@ -1786,7 +1790,12 @@ 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 (getenv("QISKIT_LAPACK_SVD")) { + 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); @@ -1811,7 +1820,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 (getenv("QISKIT_LAPACK_SVD")) { + 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_tensor.hpp b/src/simulators/matrix_product_state/matrix_product_state_tensor.hpp index 0155925843..aee7442912 100644 --- a/src/simulators/matrix_product_state/matrix_product_state_tensor.hpp +++ b/src/simulators/matrix_product_state/matrix_product_state_tensor.hpp @@ -603,7 +603,11 @@ double MPS_Tensor::Decompose(MPS_Tensor &temp, MPS_Tensor &left_gamma, left_gamma.data_ = reshape_U_after_SVD(U); lambda = S; - right_gamma.data_ = reshape_V_after_SVD(V); + if (getenv("QISKIT_LAPACK_SVD")) { + 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..a71f2b5a05 100644 --- a/src/simulators/matrix_product_state/svd.cpp +++ b/src/simulators/matrix_product_state/svd.cpp @@ -65,6 +65,12 @@ 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 @@ -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); + // V**H is not the same as V + if (getenv("QISKIT_LAPACK_SVD")) { + 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 @@ -520,6 +531,14 @@ status csvd(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) { + if (getenv("QISKIT_LAPACK_SVD")) { + 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 +571,72 @@ 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) +{ +#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[min_dim]; + complex_t work[lwork] = {0.0}; + int info; + + if (strcmp(getenv("QISKIT_LAPACK_SVD"), "DC") == 0) { + int iwork[8*min_dim] = {0}; + double rwork[5*min_dim*min_dim+5*min_dim] = {0.0}; + 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_[lwork] = {0.0}; + zgesdd_( + "A", &m, &n, lapackA, &m, lapackS, + lapackU, &m, lapackV, &n, work_, &lwork, + rwork, iwork, &info); + } else { + // Default execution follows original method + double rwork[5*min_dim] = {0.0}; + zgesvd_( + "A", "A", &m, &n, lapackA, &m, lapackS, + lapackU, &m, lapackV, &n, work, &lwork, + rwork, &info); + } + + 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 #include + #define CHOP_THRESHOLD 1e-16 namespace AER { @@ -32,11 +35,18 @@ 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); status csvd(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); +// 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); diff --git a/test/terra/backends/aer_simulator/test_options.py b/test/terra/backends/aer_simulator/test_options.py index e82c00e2d3..3f20bfe32e 100644 --- a/test/terra/backends/aer_simulator/test_options.py +++ b/test/terra/backends/aer_simulator/test_options.py @@ -282,3 +282,39 @@ def test_statevector_memory(self): "Required memory: {}".format(2 ** (n - 20) * 16) in result.results[0].status ) self.assertTrue("max memory: {}".format(max_memory_mb) in result.results[0].status) + + 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_swap = 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_swap = backend_swap.run(circuit, shots=shots).result() + original_sv = result_swap.data(0)["sv"] + + # run with lapack svd method + os.environ["QISKIT_LAPACK_SVD"] = "1" + result_swap = backend_swap.run(circuit, shots=shots).result() + lapack_sv = result_swap.data(0)["sv"] + + # run with lapack svd D&C approach + os.environ["QISKIT_LAPACK_SVD"] = "DC" + result_swap = backend_swap.run(circuit, shots=shots).result() + lapack_dc_sv = result_swap.data(0)["sv"] + + # should give the same state vector + self.assertAlmostEqual(state_fidelity(original_sv, lapack_sv), 1.0) + + # should give the same state vector + self.assertAlmostEqual(state_fidelity(original_sv, lapack_dc_sv), 1.0) From 79d05aac5ac11f559f9ed9e51004d11cec40a101 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Daniel=20Rodr=C3=ADguez=20L?= Date: Fri, 19 May 2023 15:19:36 +0900 Subject: [PATCH 02/13] forgot to unset env variable in test --- test/terra/backends/aer_simulator/test_options.py | 1 + 1 file changed, 1 insertion(+) diff --git a/test/terra/backends/aer_simulator/test_options.py b/test/terra/backends/aer_simulator/test_options.py index 3f20bfe32e..b1d0afafb8 100644 --- a/test/terra/backends/aer_simulator/test_options.py +++ b/test/terra/backends/aer_simulator/test_options.py @@ -312,6 +312,7 @@ def test_mps_svd_method(self): os.environ["QISKIT_LAPACK_SVD"] = "DC" result_swap = backend_swap.run(circuit, shots=shots).result() lapack_dc_sv = result_swap.data(0)["sv"] + os.unsetenv("QISKIT_LAPACK_SVD") # should give the same state vector self.assertAlmostEqual(state_fidelity(original_sv, lapack_sv), 1.0) From 6c3ff5cd7ec67c415c6ae4b41bd45710b8228663 Mon Sep 17 00:00:00 2001 From: patataman Date: Wed, 14 Jun 2023 14:27:47 +0900 Subject: [PATCH 03/13] fix seg fault bc arrays were too big --- src/simulators/matrix_product_state/svd.cpp | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/src/simulators/matrix_product_state/svd.cpp b/src/simulators/matrix_product_state/svd.cpp index a71f2b5a05..746d72b278 100644 --- a/src/simulators/matrix_product_state/svd.cpp +++ b/src/simulators/matrix_product_state/svd.cpp @@ -576,6 +576,7 @@ void lapack_csvd_wrapper(cmatrix_t &A, cmatrix_t &U, rvector_t &S, cmatrix_t &V) #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); @@ -589,25 +590,32 @@ void lapack_csvd_wrapper(cmatrix_t &A, cmatrix_t &U, rvector_t &S, cmatrix_t &V) *lapackV = V.move_to_buffer(); double lapackS[min_dim]; - complex_t work[lwork] = {0.0}; + complex_t work[lwork]; int info; if (strcmp(getenv("QISKIT_LAPACK_SVD"), "DC") == 0) { - int iwork[8*min_dim] = {0}; - double rwork[5*min_dim*min_dim+5*min_dim] = {0.0}; - lwork = -1; + int iwork[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_[lwork] = {0.0}; + 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); + + free(rwork); + free(work_); } else { // Default execution follows original method double rwork[5*min_dim] = {0.0}; @@ -616,7 +624,6 @@ void lapack_csvd_wrapper(cmatrix_t &A, cmatrix_t &U, rvector_t &S, cmatrix_t &V) lapackU, &m, lapackV, &n, work, &lwork, rwork, &info); } - 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); From c673b0d690518b3dcc6487d566ab39a32b2006aa Mon Sep 17 00:00:00 2001 From: patataman Date: Thu, 10 Aug 2023 11:49:15 +0900 Subject: [PATCH 04/13] code style and releasenote for PR --- ...pute-svd-with-lapack-3ee992d371d653d1.yaml | 11 +++++ src/framework/lapack_protos.hpp | 21 ++++---- .../matrix_product_state_tensor.hpp | 2 +- src/simulators/matrix_product_state/svd.cpp | 49 ++++++++----------- src/simulators/matrix_product_state/svd.hpp | 9 ++-- .../backends/aer_simulator/test_options.py | 3 +- 6 files changed, 50 insertions(+), 45 deletions(-) create mode 100644 releasenotes/notes/compute-svd-with-lapack-3ee992d371d653d1.yaml 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..1fbeba519b --- /dev/null +++ b/releasenotes/notes/compute-svd-with-lapack-3ee992d371d653d1.yaml @@ -0,0 +1,11 @@ +--- +features: + - | + Replace Qiskit SVD function with OpenBLAS/LAPACK SVD functions ``zgesvd`` + and ``zgesdd``. By default ``zgesvd`` is used, but performace decrease with + big matrices. On the other hand, ``zgesdd`` has bad performance with small + matrices, and better with big matrices. User can use ``zgesdd`` function + setting the environment variable ``QISKIT_LAPACK_SVD=DC``. +other: + - | + Performance may decrease depending on the SVD method used. diff --git a/src/framework/lapack_protos.hpp b/src/framework/lapack_protos.hpp index 6ca64047f0..18967f6496 100644 --- a/src/framework/lapack_protos.hpp +++ b/src/framework/lapack_protos.hpp @@ -3,10 +3,10 @@ #ifndef _aer_framework_lapack_protos_hpp #define _aer_framework_lapack_protos_hpp +#include #include #include #include -#include #ifdef __cplusplus extern "C" { @@ -15,19 +15,20 @@ extern "C" { // 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); + 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); + 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 } 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 aee7442912..16009ef1cd 100644 --- a/src/simulators/matrix_product_state/matrix_product_state_tensor.hpp +++ b/src/simulators/matrix_product_state/matrix_product_state_tensor.hpp @@ -604,7 +604,7 @@ double MPS_Tensor::Decompose(MPS_Tensor &temp, MPS_Tensor &left_gamma, left_gamma.data_ = reshape_U_after_SVD(U); lambda = S; if (getenv("QISKIT_LAPACK_SVD")) { - right_gamma.data_ = reshape_VH_after_SVD(V); + right_gamma.data_ = reshape_VH_after_SVD(V); } else { right_gamma.data_ = reshape_V_after_SVD(V); } diff --git a/src/simulators/matrix_product_state/svd.cpp b/src/simulators/matrix_product_state/svd.cpp index 746d72b278..f594b207c4 100644 --- a/src/simulators/matrix_product_state/svd.cpp +++ b/src/simulators/matrix_product_state/svd.cpp @@ -65,8 +65,7 @@ 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 reshape_VH_after_SVD(const cmatrix_t V) { std::vector Res(2); AER::Utils::split(V, Res[0], Res[1], 1); return Res; @@ -538,7 +537,8 @@ void csvd_wrapper(cmatrix_t &A, cmatrix_t &U, rvector_t &S, cmatrix_t &V) { } } -void qiskit_csvd_wrapper(cmatrix_t &A, cmatrix_t &U, rvector_t &S, cmatrix_t &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 @@ -571,8 +571,8 @@ void qiskit_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) -{ +void lapack_csvd_wrapper(cmatrix_t &A, cmatrix_t &U, rvector_t &S, + cmatrix_t &V) { #ifdef DEBUG cmatrix_t tempA = A; #endif @@ -585,51 +585,43 @@ void lapack_csvd_wrapper(cmatrix_t &A, cmatrix_t &U, rvector_t &S, cmatrix_t &V) 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(); + complex_t *lapackA = A.move_to_buffer(), *lapackU = U.move_to_buffer(), + *lapackV = V.move_to_buffer(); double lapackS[min_dim]; complex_t work[lwork]; int info; if (strcmp(getenv("QISKIT_LAPACK_SVD"), "DC") == 0) { - int iwork[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); + int iwork[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)); + 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); + 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)); + 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); + zgesdd_("A", &m, &n, lapackA, &m, lapackS, lapackU, &m, lapackV, &n, work_, + &lwork, rwork, iwork, &info); free(rwork); free(work_); } else { // Default execution follows original method - double rwork[5*min_dim] = {0.0}; - zgesvd_( - "A", "A", &m, &n, lapackA, &m, lapackS, - lapackU, &m, lapackV, &n, work, &lwork, - rwork, &info); + double rwork[5 * min_dim] = {0.0}; + zgesvd_("A", "A", &m, &n, lapackA, &m, lapackS, lapackU, &m, lapackV, &n, + work, &lwork, rwork, &info); } 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 #include - #define CHOP_THRESHOLD 1e-16 namespace AER { @@ -43,9 +42,11 @@ status csvd(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); // Original qiskit call -void qiskit_csvd_wrapper(cmatrix_t &C, cmatrix_t &U, rvector_t &S, cmatrix_t &V); +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 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); diff --git a/test/terra/backends/aer_simulator/test_options.py b/test/terra/backends/aer_simulator/test_options.py index d1e1112f5a..d7eae6c560 100644 --- a/test/terra/backends/aer_simulator/test_options.py +++ b/test/terra/backends/aer_simulator/test_options.py @@ -307,6 +307,7 @@ 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_swap = self.backend(method=method) @@ -338,4 +339,4 @@ def test_mps_svd_method(self): self.assertAlmostEqual(state_fidelity(original_sv, lapack_sv), 1.0) # should give the same state vector - self.assertAlmostEqual(state_fidelity(original_sv, lapack_dc_sv), 1.0) \ No newline at end of file + self.assertAlmostEqual(state_fidelity(original_sv, lapack_dc_sv), 1.0) From 180359bb9c646302b6a04fc657f398a029f106c6 Mon Sep 17 00:00:00 2001 From: patataman Date: Thu, 24 Aug 2023 01:49:43 +0000 Subject: [PATCH 05/13] address microsoft C2131? --- src/simulators/matrix_product_state/svd.cpp | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/simulators/matrix_product_state/svd.cpp b/src/simulators/matrix_product_state/svd.cpp index f594b207c4..2568fd1b09 100644 --- a/src/simulators/matrix_product_state/svd.cpp +++ b/src/simulators/matrix_product_state/svd.cpp @@ -588,12 +588,12 @@ void lapack_csvd_wrapper(cmatrix_t &A, cmatrix_t &U, rvector_t &S, complex_t *lapackA = A.move_to_buffer(), *lapackU = U.move_to_buffer(), *lapackV = V.move_to_buffer(); - double lapackS[min_dim]; - complex_t work[lwork]; + double* lapackS = new double[min_dim]; + complex_t* work = new complex_t[lwork]; int info; if (strcmp(getenv("QISKIT_LAPACK_SVD"), "DC") == 0) { - int iwork[8 * min_dim]; + 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); @@ -608,6 +608,7 @@ void lapack_csvd_wrapper(cmatrix_t &A, cmatrix_t &U, rvector_t &S, zgesdd_("A", &m, &n, lapackA, &m, lapackS, lapackU, &m, lapackV, &n, work_, &lwork, rwork, iwork, &info); + delete iwork; free(rwork); free(work_); } else { @@ -628,6 +629,9 @@ void lapack_csvd_wrapper(cmatrix_t &A, cmatrix_t &U, rvector_t &S, validate_SVD_result(tempA, U, S, V); #endif + delete lapackS; + delete work; + if (info == 0) { return; } else { From f37773354b508108a5981f3de65000626cceaa39 Mon Sep 17 00:00:00 2001 From: patataman Date: Thu, 24 Aug 2023 02:38:35 +0000 Subject: [PATCH 06/13] style + another C2131 --- src/simulators/matrix_product_state/svd.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/simulators/matrix_product_state/svd.cpp b/src/simulators/matrix_product_state/svd.cpp index 2568fd1b09..54a2f84ea7 100644 --- a/src/simulators/matrix_product_state/svd.cpp +++ b/src/simulators/matrix_product_state/svd.cpp @@ -588,12 +588,12 @@ void lapack_csvd_wrapper(cmatrix_t &A, cmatrix_t &U, rvector_t &S, 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]; + double *lapackS = new double[min_dim]; + complex_t *work = new complex_t[lwork]; int info; if (strcmp(getenv("QISKIT_LAPACK_SVD"), "DC") == 0) { - int* iwork = new int[8 * min_dim]; + 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); @@ -613,7 +613,7 @@ void lapack_csvd_wrapper(cmatrix_t &A, cmatrix_t &U, rvector_t &S, free(work_); } else { // Default execution follows original method - double rwork[5 * min_dim] = {0.0}; + 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); } From 5a91443d32cb18ef846ef23de31ba647a1725cca Mon Sep 17 00:00:00 2001 From: patataman Date: Thu, 24 Aug 2023 02:40:10 +0000 Subject: [PATCH 07/13] missing free --- src/simulators/matrix_product_state/svd.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/simulators/matrix_product_state/svd.cpp b/src/simulators/matrix_product_state/svd.cpp index 54a2f84ea7..ceb1fe98b4 100644 --- a/src/simulators/matrix_product_state/svd.cpp +++ b/src/simulators/matrix_product_state/svd.cpp @@ -616,6 +616,7 @@ void lapack_csvd_wrapper(cmatrix_t &A, cmatrix_t &U, rvector_t &S, 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); From 394298ea81c073f0b402de6cdb0320a8846dddc5 Mon Sep 17 00:00:00 2001 From: patataman Date: Thu, 24 Aug 2023 13:36:16 +0900 Subject: [PATCH 08/13] change test to support windows python3.8 --- test/terra/backends/aer_simulator/test_options.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/test/terra/backends/aer_simulator/test_options.py b/test/terra/backends/aer_simulator/test_options.py index d7eae6c560..bb7c4779b1 100644 --- a/test/terra/backends/aer_simulator/test_options.py +++ b/test/terra/backends/aer_simulator/test_options.py @@ -333,7 +333,9 @@ def test_mps_svd_method(self): os.environ["QISKIT_LAPACK_SVD"] = "DC" result_swap = backend_swap.run(circuit, shots=shots).result() lapack_dc_sv = result_swap.data(0)["sv"] - os.unsetenv("QISKIT_LAPACK_SVD") + + del os.environ["QISKIT_LAPACK_SVD"] # Python 3.8 in Windows + # os.unsetenv("QISKIT_LAPACK_SVD") # Since Python 3.9 # should give the same state vector self.assertAlmostEqual(state_fidelity(original_sv, lapack_sv), 1.0) From f7190fc0d97a53c3beaef0fd64ad25314e37083a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Daniel=20Rodr=C3=ADguez=20L?= Date: Mon, 2 Oct 2023 16:54:26 +0200 Subject: [PATCH 09/13] Update releasenotes/notes/compute-svd-with-lapack-3ee992d371d653d1.yaml Co-authored-by: merav-aharoni --- .../notes/compute-svd-with-lapack-3ee992d371d653d1.yaml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/releasenotes/notes/compute-svd-with-lapack-3ee992d371d653d1.yaml b/releasenotes/notes/compute-svd-with-lapack-3ee992d371d653d1.yaml index 1fbeba519b..aec61f076b 100644 --- a/releasenotes/notes/compute-svd-with-lapack-3ee992d371d653d1.yaml +++ b/releasenotes/notes/compute-svd-with-lapack-3ee992d371d653d1.yaml @@ -2,9 +2,9 @@ features: - | Replace Qiskit SVD function with OpenBLAS/LAPACK SVD functions ``zgesvd`` - and ``zgesdd``. By default ``zgesvd`` is used, but performace decrease with - big matrices. On the other hand, ``zgesdd`` has bad performance with small - matrices, and better with big matrices. User can use ``zgesdd`` function + 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``. other: - | From 08756076e864004243022c908e5beeb2f357c3f1 Mon Sep 17 00:00:00 2001 From: patataman Date: Mon, 2 Oct 2023 19:29:13 +0200 Subject: [PATCH 10/13] remove unnecessary comments --- .../notes/compute-svd-with-lapack-3ee992d371d653d1.yaml | 3 --- 1 file changed, 3 deletions(-) diff --git a/releasenotes/notes/compute-svd-with-lapack-3ee992d371d653d1.yaml b/releasenotes/notes/compute-svd-with-lapack-3ee992d371d653d1.yaml index aec61f076b..f861b6fbab 100644 --- a/releasenotes/notes/compute-svd-with-lapack-3ee992d371d653d1.yaml +++ b/releasenotes/notes/compute-svd-with-lapack-3ee992d371d653d1.yaml @@ -6,6 +6,3 @@ features: 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``. -other: - - | - Performance may decrease depending on the SVD method used. From 0503475dea9dfbd677e3de7ee763340ed9b37543 Mon Sep 17 00:00:00 2001 From: patataman Date: Mon, 2 Oct 2023 19:29:21 +0200 Subject: [PATCH 11/13] undo ifdef DEBUG --- .../matrix_product_state/matrix_product_state_internal.cpp | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) 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 d34ed843cc..f9346a34c1 100644 --- a/src/simulators/matrix_product_state/matrix_product_state_internal.cpp +++ b/src/simulators/matrix_product_state/matrix_product_state_internal.cpp @@ -542,9 +542,7 @@ void MPS::apply_swap_internal(uint_t index_A, uint_t index_B, bool swap_gate) { // to the left std::swap(qubit_ordering_.order_[index_A], qubit_ordering_.order_[index_B]); // For logging purposes: -#ifdef DEBUG print_to_log_internal_swap(index_A, index_B); -#endif // update qubit locations after all the swaps for (uint_t i = 0; i < num_qubits_; i++) @@ -666,10 +664,9 @@ void MPS::common_apply_2_qubit_gate( rvector_t lambda; double discarded_value = MPS_Tensor::Decompose(temp, left_gamma, lambda, right_gamma); -#ifdef DEBUG + if (discarded_value > json_chop_threshold_) MPS::print_to_log("discarded_value=", discarded_value, ", "); -#endif if (A != 0) left_gamma.div_Gamma_by_left_Lambda(lambda_reg_[A - 1]); From fda3db1084b22919335aeed84ba21cacf4925aea Mon Sep 17 00:00:00 2001 From: patataman Date: Mon, 2 Oct 2023 19:29:28 +0200 Subject: [PATCH 12/13] automatic selector for QR or D&C in LAPACK SVD --- src/simulators/matrix_product_state/svd.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/simulators/matrix_product_state/svd.cpp b/src/simulators/matrix_product_state/svd.cpp index ceb1fe98b4..154a5c12b5 100644 --- a/src/simulators/matrix_product_state/svd.cpp +++ b/src/simulators/matrix_product_state/svd.cpp @@ -112,7 +112,7 @@ 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**H is not the same as V + // When using LAPACK function, V is V dagger if (getenv("QISKIT_LAPACK_SVD")) { V.resize(new_SV_num, V.GetColumns()); } else { @@ -592,7 +592,7 @@ void lapack_csvd_wrapper(cmatrix_t &A, cmatrix_t &U, rvector_t &S, complex_t *work = new complex_t[lwork]; int info; - if (strcmp(getenv("QISKIT_LAPACK_SVD"), "DC") == 0) { + if (m >= 64 and n >= 64) { // strcmp(getenv("QISKIT_LAPACK_SVD"), "DC") == 0) { 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); From 0899baaf840df930ed5f97e536a0d1034ab63b8b Mon Sep 17 00:00:00 2001 From: patataman Date: Thu, 26 Oct 2023 17:19:02 +0200 Subject: [PATCH 13/13] codestyle, enable MPS lapack using run_options --- 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 ++ src/framework/config.hpp | 6 ++- .../matrix_product_state.hpp | 6 ++- .../matrix_product_state_internal.cpp | 13 +++--- .../matrix_product_state_internal.hpp | 5 +++ .../matrix_product_state_tensor.hpp | 14 +++--- src/simulators/matrix_product_state/svd.cpp | 45 ++++++++++++++----- src/simulators/matrix_product_state/svd.hpp | 8 +++- .../backends/aer_simulator/test_options.py | 27 +++++------ 12 files changed, 91 insertions(+), 42 deletions(-) diff --git a/qiskit_aer/backends/aer_compiler.py b/qiskit_aer/backends/aer_compiler.py index e4a3a4e9b6..0c04c1f831 100644 --- a/qiskit_aer/backends/aer_compiler.py +++ b/qiskit_aer/backends/aer_compiler.py @@ -486,6 +486,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 f845ecd6f0..6f224e49ca 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: @@ -768,6 +771,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 1901fa066f..aa0f788ec9 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 f614e4483d..a692915df5 100644 --- a/qiskit_aer/backends/wrappers/aer_controller_binding.hpp +++ b/qiskit_aer/backends/wrappers/aer_controller_binding.hpp @@ -232,6 +232,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); @@ -467,6 +468,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), @@ -561,6 +563,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/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/simulators/matrix_product_state/matrix_product_state.hpp b/src/simulators/matrix_product_state/matrix_product_state.hpp index 1c29c9bd02..64f734b96b 100644 --- a/src/simulators/matrix_product_state/matrix_product_state.hpp +++ b/src/simulators/matrix_product_state/matrix_product_state.hpp @@ -365,6 +365,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 { @@ -376,6 +379,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 { @@ -813,4 +817,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 f9346a34c1..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,8 @@ 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, ", "); @@ -1788,7 +1789,7 @@ void MPS::initialize_from_matrix(uint_t num_qubits, const cmatrix_t &mat) { remaining_matrix = mat; } else { cmatrix_t temp; - if (getenv("QISKIT_LAPACK_SVD")) { + 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); @@ -1799,9 +1800,9 @@ void MPS::initialize_from_matrix(uint_t num_qubits, const cmatrix_t &mat) { // 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 @@ -1818,7 +1819,7 @@ void MPS::initialize_from_matrix(uint_t num_qubits, const cmatrix_t &mat) { } // step 4 - create the rightmost gamma and update q_reg_ std::vector right_data; - if (getenv("QISKIT_LAPACK_SVD")) { + if (MPS::mps_lapack_) { right_data = reshape_VH_after_SVD(V); } else { right_data = reshape_V_after_SVD(V); 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 16009ef1cd..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,21 @@ 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; - if (getenv("QISKIT_LAPACK_SVD")) { + 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); diff --git a/src/simulators/matrix_product_state/svd.cpp b/src/simulators/matrix_product_state/svd.cpp index 154a5c12b5..150d97af4f 100644 --- a/src/simulators/matrix_product_state/svd.cpp +++ b/src/simulators/matrix_product_state/svd.cpp @@ -90,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; @@ -113,7 +114,7 @@ double reduce_zeros(cmatrix_t &U, rvector_t &S, cmatrix_t &V, U.resize(U.GetRows(), new_SV_num); S.resize(new_SV_num); // When using LAPACK function, V is V dagger - if (getenv("QISKIT_LAPACK_SVD")) { + if (mps_lapack) { V.resize(new_SV_num, V.GetColumns()); } else { V.resize(V.GetRows(), new_SV_num); @@ -140,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); @@ -529,8 +549,9 @@ 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) { - if (getenv("QISKIT_LAPACK_SVD")) { +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); @@ -573,9 +594,10 @@ void qiskit_csvd_wrapper(cmatrix_t &A, cmatrix_t &U, rvector_t &S, void lapack_csvd_wrapper(cmatrix_t &A, cmatrix_t &U, rvector_t &S, cmatrix_t &V) { -#ifdef DEBUG + // Activated by default as requested in the PR + // #ifdef DEBUG cmatrix_t tempA = A; -#endif + // #endif const size_t m = A.GetRows(), n = A.GetColumns(); const size_t min_dim = std::min(m, n); @@ -592,7 +614,9 @@ void lapack_csvd_wrapper(cmatrix_t &A, cmatrix_t &U, rvector_t &S, complex_t *work = new complex_t[lwork]; int info; - if (m >= 64 and n >= 64) { // strcmp(getenv("QISKIT_LAPACK_SVD"), "DC") == 0) { + 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); @@ -626,9 +650,10 @@ void lapack_csvd_wrapper(cmatrix_t &A, cmatrix_t &U, rvector_t &S, for (int i = 0; i < min_dim; i++) S.push_back(lapackS[i]); -#ifdef DEBUG - validate_SVD_result(tempA, U, S, V); -#endif + // Activated by default as requested in the PR + // #ifdef DEBUG + validate_SVdD_result(tempA, U, S, V); + // #endif delete lapackS; delete work; diff --git a/src/simulators/matrix_product_state/svd.hpp b/src/simulators/matrix_product_state/svd.hpp index 729b42dc91..fac77797fb 100644 --- a/src/simulators/matrix_product_state/svd.hpp +++ b/src/simulators/matrix_product_state/svd.hpp @@ -37,10 +37,12 @@ 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); // Entry point for the SVD calculation -void csvd_wrapper(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, + bool lapack); // Original qiskit call void qiskit_csvd_wrapper(cmatrix_t &C, cmatrix_t &U, rvector_t &S, cmatrix_t &V); @@ -50,6 +52,8 @@ void lapack_csvd_wrapper(cmatrix_t &C, cmatrix_t &U, rvector_t &S, 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 bb7c4779b1..51dff7d4be 100644 --- a/test/terra/backends/aer_simulator/test_options.py +++ b/test/terra/backends/aer_simulator/test_options.py @@ -310,10 +310,12 @@ def test_mps_svd_method(self): shots = 4000 method = "matrix_product_state" - backend_swap = self.backend(method=method) + 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]) @@ -321,24 +323,17 @@ def test_mps_svd_method(self): circuit.cx(0, i) circuit.save_statevector("sv") - result_swap = backend_swap.run(circuit, shots=shots).result() - original_sv = result_swap.data(0)["sv"] + result_og = backend_og.run(circuit, shots=shots).result() + original_sv = result_og.data(0)["sv"] # run with lapack svd method - os.environ["QISKIT_LAPACK_SVD"] = "1" - result_swap = backend_swap.run(circuit, shots=shots).result() - lapack_sv = result_swap.data(0)["sv"] - - # run with lapack svd D&C approach - os.environ["QISKIT_LAPACK_SVD"] = "DC" - result_swap = backend_swap.run(circuit, shots=shots).result() - lapack_dc_sv = result_swap.data(0)["sv"] + result_lapack = backend_lapack.run(circuit, shots=shots, mps_lapack=True).result() + lapack_sv = result_lapack.data(0)["sv"] - del os.environ["QISKIT_LAPACK_SVD"] # Python 3.8 in Windows - # os.unsetenv("QISKIT_LAPACK_SVD") # Since Python 3.9 + # 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) - - # should give the same state vector - self.assertAlmostEqual(state_fidelity(original_sv, lapack_dc_sv), 1.0)