Add support for BLAS SVD functions in MPS simulation

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),

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,

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

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);

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
+#endif

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<uint_t, double> State::sample_measure_with_prob(const reg_t &qubits,
 //-------------------------------------------------------------------------
 } // end namespace AER
 //-------------------------------------------------------------------------
-#endif
+#endif

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<cmatrix_t> 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);

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) {

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<cmatrix_t> 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);

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;

src/simulators/matrix_product_state/svd.hpp

@@ -37,10 +37,12 @@ std::vector<cmatrix_t> reshape_V_after_SVD(const cmatrix_t V);
 std::vector<cmatrix_t> 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