Add support for BLAS SVD functions in MPS simulation

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

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,

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

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

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``.

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/framework/lapack_protos.hpp

@@ -0,0 +1,37 @@
+// Dependencies: BLAS - LAPACK
+
+#ifndef _aer_framework_lapack_protos_hpp
+#define _aer_framework_lapack_protos_hpp
+
+#include <array>
+#include <complex>
+#include <iostream>
+#include <vector>
+
+#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<double> *a, const size_t *lda,
+             double *s, std::complex<double> *u, const size_t *ldu,
+             std::complex<double> *vt, const size_t *ldvt,
+             std::complex<double> *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<double> *a, const size_t *lda, double *s,
+             std::complex<double> *u, const size_t *ldu,
+             std::complex<double> *vt, const size_t *ldvt,
+             std::complex<double> *work, const size_t *lwork, double *rwork,
+             int *iwork, int *info);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif // end __lapack_protos_h_

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<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,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<cmatrix_t> right_data = reshape_V_after_SVD(V);
+  std::vector<cmatrix_t> 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);

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,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;
 }