Skip to content

Commit

Permalink
MPS: Added thresholds to loops that use parallelization (Qiskit#638)
Browse files Browse the repository at this point in the history 
* Adding conditions on invoking parallelization
  • Loading branch information
@merav-aharoni
merav-aharoni authored Mar 2, 2020
1 parent a965674 commit 7fb9147
Show file tree
Hide file tree
Showing 4 changed files with 100 additions and 101 deletions.
43 changes: 13 additions & 30 deletions src/simulators/matrix_product_state/matrix_product_state.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -271,15 +271,6 @@ class State : public Base::State<matrixproductstate_t> {
// Config Settings
//-----------------------------------------------------------------------

// OpenMP qubit threshold
int omp_qubit_threshold_ = 14;

// QubitVector sample measure index size
int sample_measure_index_size_ = 10;

// Threshold for chopping small values to zero in JSON
double json_chop_threshold_ = 1e-15;

// Table of allowed gate names to gate enum class members
const static stringmap_t<Gates> gateset_;

Expand Down Expand Up @@ -343,7 +334,6 @@ const stringmap_t<Snapshots> State::snapshotset_({
//-------------------------------------------------------------------------

void State::initialize_qreg(uint_t num_qubits) {
initialize_omp();
qreg_.initialize((uint_t)num_qubits);
}

Expand All @@ -352,8 +342,6 @@ void State::initialize_qreg(uint_t num_qubits, const matrixproductstate_t &state
if (qreg_.num_qubits() != num_qubits) {
throw std::invalid_argument("MatrixProductState::State::initialize: initial state does not match qubit number");
}
initialize_omp();
//qreg_.initialize((uint_t)num_qubits, state);
#ifdef DEBUG
cout << "initialize with state not supported yet";
#endif
Expand All @@ -364,7 +352,6 @@ void State::initialize_qreg(uint_t num_qubits, const cvector_t &statevector) {
if (qreg_.num_qubits() != num_qubits) {
throw std::invalid_argument("MatrixProductState::State::initialize: initial state does not match qubit number");
}
initialize_omp();

// internal bit ordering is the opposite of ordering in Qasm, so must
// reverse order before starting
Expand All @@ -373,12 +360,6 @@ void State::initialize_qreg(uint_t num_qubits, const cvector_t &statevector) {
qreg_.initialize_from_statevector(num_qubits, mps_format_state_vector);
}

void State::initialize_omp() {
qreg_.set_omp_threshold(omp_qubit_threshold_);
if (BaseState::threads_ > 0)
qreg_.set_omp_threads(BaseState::threads_); // set allowed OMP threads in MPS
}

size_t State::required_memory_mb(uint_t num_qubits,
const std::vector<Operations::Op> &ops) const {
// for each qubit we have a tensor structure.
Expand All @@ -394,23 +375,25 @@ size_t State::required_memory_mb(uint_t num_qubits,
void State::set_config(const json_t &config) {

// Set threshold for truncating snapshots
JSON::get_value(json_chop_threshold_, "chop_threshold", config);
qreg_.set_json_chop_threshold(json_chop_threshold_);
uint_t json_chop_threshold;
if (JSON::get_value(json_chop_threshold, "chop_threshold", config))
MPS::set_json_chop_threshold(json_chop_threshold);

// Set OMP threshold for state update functions
JSON::get_value(omp_qubit_threshold_, "statevector_parallel_threshold", config);
uint_t omp_qubit_threshold;
if (JSON::get_value(omp_qubit_threshold, "mps_parallel_threshold", config))
MPS::set_omp_threshold(omp_qubit_threshold);

// Set the sample measure indexing size
int index_size;
if (JSON::get_value(index_size, "statevector_sample_measure_opt", config)) {
qreg_.set_sample_measure_index_size(index_size);
if (JSON::get_value(index_size, "mps_sample_measure_opt", config)) {
MPS::set_sample_measure_index_size(index_size);
};

// Enable sorted gate optimzations
bool gate_opt = false;
JSON::get_value(gate_opt, "statevector_gate_opt", config);
if (gate_opt)
qreg_.enable_gate_opt();
// if (JSON::get_value(gate_opt, "mps_gate_opt", config))
// MPS::set_enable_gate_opt(gate_opt);
}

//=========================================================================
Expand Down Expand Up @@ -481,7 +464,7 @@ void State::snapshot_pauli_expval(const Operations::Op &op,
}

// add to snapshot
Utils::chop_inplace(expval, json_chop_threshold_);
Utils::chop_inplace(expval, MPS::get_json_chop_threshold());
switch (type) {
case SnapshotDataType::average:
data.add_average_snapshot("expectation_value", op.string_params[0],
Expand Down Expand Up @@ -521,7 +504,7 @@ void State::snapshot_matrix_expval(const Operations::Op &op,
}

// add to snapshot
Utils::chop_inplace(expval, json_chop_threshold_);
Utils::chop_inplace(expval, MPS::get_json_chop_threshold());
switch (type) {
case SnapshotDataType::average:
data.add_average_snapshot("expectation_value", op.string_params[0],
Expand Down Expand Up @@ -550,7 +533,7 @@ void State::snapshot_probabilities(const Operations::Op &op,
SnapshotDataType type) {
rvector_t prob_vector;
qreg_.get_probabilities_vector(prob_vector, op.qubits);
auto probs = Utils::vec2ket(prob_vector, json_chop_threshold_, 16);
auto probs = Utils::vec2ket(prob_vector, MPS::get_json_chop_threshold(), 16);
bool variance = type == SnapshotDataType::average_var;
data.add_average_snapshot("probabilities", op.string_params[0],
BaseState::creg_.memory_hex(), probs, variance);
Expand Down
50 changes: 27 additions & 23 deletions src/simulators/matrix_product_state/matrix_product_state_internal.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,10 @@ static const cmatrix_t zero_measure =
static const cmatrix_t one_measure =
AER::Utils::make_matrix<complex_t>({{{0, 0}, {0, 0}},
{{0, 0}, {1, 0}}});

uint_t MPS::omp_threads_ = 1;
uint_t MPS::omp_threshold_ = 14;
int MPS::sample_measure_index_size_ = 10;
double MPS::json_chop_threshold_ = 1E-8;
//------------------------------------------------------------------------
// local function declarations
//------------------------------------------------------------------------
Expand Down Expand Up @@ -162,7 +165,6 @@ uint_t reorder_qubits(const reg_t qubits, uint_t index) {

uint_t reverse_bits(uint_t num, uint_t len) {
uint_t sum = 0;
// std::assert(num < pow(2, len));
for (uint_t i=0; i<len; ++i) {
if ((num & 0x1) == 1) {
sum += 1ULL << (len-1-i); // adding pow(2, len-1-i)
Expand All @@ -180,7 +182,8 @@ std::vector<T> reverse_all_bits(const std::vector<T>& statevector, uint_t num_qu
{
uint_t length = statevector.size(); // length = pow(2, num_qubits_)
std::vector<T> output_vector(length);
#pragma omp parallel for

#pragma omp parallel for if (length > MPS::get_omp_threshold() && MPS::get_omp_threads() > 1) num_threads(MPS::get_omp_threads())
for (int_t i = 0; i < static_cast<int_t>(length); i++) {
output_vector[i] = statevector[reverse_bits(i, num_qubits)];
}
Expand All @@ -199,17 +202,16 @@ std::vector<uint_t> calc_new_indices(const reg_t &indices) {
}

cmatrix_t mul_matrix_by_lambda(const cmatrix_t &mat,
const rvector_t &lambda)
{
const rvector_t &lambda) {
if (lambda == rvector_t {1.0}) return mat;
cmatrix_t res_mat(mat);
uint_t num_rows = mat.GetRows(), num_cols = mat.GetColumns();

#ifdef _WIN32
#pragma omp parallel for
#else
#pragma omp parallel for collapse(2)
#endif
#ifdef _WIN32
#pragma omp parallel for if (num_rows*num_cols > 64 && MPS::get_omp_threads() > 1) num_threads(MPS::get_omp_threads())
#else
#pragma omp parallel for collapse(2) if (num_rows*num_cols > 64 && MPS::get_omp_threads() > 1) num_threads(MPS::get_omp_threads())
#endif
for(int_t row = 0; row < static_cast<int_t>(num_rows); row++) {
for(int_t col = 0; col < static_cast<int_t>(num_cols); col++) {
res_mat(row, col) = mat(row, col) * lambda[col];
Expand Down Expand Up @@ -752,11 +754,12 @@ cmatrix_t MPS::density_matrix(const reg_t &qubits) const
MPS_Tensor psi = temp_MPS.state_vec_as_MPS(new_qubits.front(), new_qubits.back());
uint_t size = psi.get_dim();
cmatrix_t rho(size,size);
#ifdef _WIN32
#pragma omp parallel for
#else
#pragma omp parallel for collapse(2)
#endif

#ifdef _WIN32
#pragma omp parallel for if (size > omp_threshold_ && omp_threads_ > 1) num_threads(omp_threads_)
#else
#pragma omp parallel for collapse(2) if (size > omp_threshold_ && omp_threads_ > 1) num_threads(omp_threads_)
#endif
for(int_t i = 0; i < static_cast<int_t>(size); i++) {
for(int_t j = 0; j < static_cast<int_t>(size); j++) {
rho(i,j) = AER::Utils::sum( AER::Utils::elementwise_multiplication(psi.get_data(i), AER::Utils::conjugate(psi.get_data(j))) );
Expand Down Expand Up @@ -911,6 +914,7 @@ complex_t MPS::expectation_value_pauli(const reg_t &qubits, const std::string &m
left_tensor.mul_Gamma_by_right_Lambda(temp_MPS.lambda_reg_[first_index]);
}


// Step 2 - prepare the dagger of left_tensor
MPS_Tensor left_tensor_dagger(AER::Utils::dagger(left_tensor.get_data(0)),
AER::Utils::dagger(left_tensor.get_data(1)));
Expand All @@ -922,11 +926,10 @@ complex_t MPS::expectation_value_pauli(const reg_t &qubits, const std::string &m
// Before contraction, Gamma0' has size a1 x a0 x i, Gamma0 has size i x a0 x a1
// result = left_contract is a matrix of size a1 x a1
cmatrix_t final_contract;
MPS_Tensor::contract_2_dimensions(left_tensor_dagger, left_tensor,
MPS_Tensor::contract_2_dimensions(left_tensor_dagger, left_tensor, omp_threads_,
final_contract);

for (uint_t qubit_num=first_index+1; qubit_num<=last_index; qubit_num++) {

// Step 5 - multiply next Gamma by its left lambda (same as Step 1)
// next gamma has dimensions a0 x a1 x i
MPS_Tensor next_gamma = temp_MPS.q_reg_[qubit_num];
Expand Down Expand Up @@ -955,8 +958,9 @@ complex_t MPS::expectation_value_pauli(const reg_t &qubits, const std::string &m
// with next_gamma_dagger (i x a2 x a1) (same as Step 4)
// here we need to contract across two dimensions: a1 and i
// result is a matrix of size a2 x a2
MPS_Tensor::contract_2_dimensions(next_gamma_dagger, next_contract,
final_contract);
MPS_Tensor::contract_2_dimensions(next_gamma_dagger, next_contract, omp_threads_,
final_contract);

}

// Step 10 - contract over final matrix of size aN x aN
Expand Down Expand Up @@ -1030,10 +1034,12 @@ void MPS::full_state_vector(cvector_t& statevector) const
MPS_Tensor mps_vec = state_vec_as_MPS(0, num_qubits_-1);
uint_t length = 1ULL << num_qubits_; // length = pow(2, num_qubits_)
statevector.resize(length);
#pragma omp parallel for

#pragma omp parallel for if (num_qubits_ > omp_threshold_ && omp_threads_ > 1) num_threads(omp_threads_)
for (int_t i = 0; i < static_cast<int_t>(length); i++) {
statevector[i] = mps_vec.get_data(reverse_bits(i, num_qubits_))(0,0);
}

#ifdef DEBUG
std::cout << *this;
#endif
Expand Down Expand Up @@ -1076,7 +1082,7 @@ uint_t MPS::apply_measure(uint_t qubit,

// step 1 - measure qubit 0 in Z basis
double exp_val = real(expectation_value_pauli(qubits_to_update, "Z"));

// step 2 - compute probability for 0 or 1 result
double prob0 = (1 + exp_val ) / 2;
double prob1 = 1 - prob0;
Expand All @@ -1095,7 +1101,6 @@ uint_t MPS::apply_measure(uint_t qubit,
measurement_matrix = one_measure;
measurement_matrix = measurement_matrix * (1 / sqrt(prob1));
}

apply_matrix(qubits_to_update, measurement_matrix);

// step 4 - propagate the changes to all qubits to the right
Expand All @@ -1111,7 +1116,6 @@ uint_t MPS::apply_measure(uint_t qubit,
break; // no need to propagate if no entanglement
apply_2_qubit_gate(i-1, i, id, cmatrix_t(1));
}

return measurement;
}

Expand Down
51 changes: 33 additions & 18 deletions src/simulators/matrix_product_state/matrix_product_state_internal.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -196,27 +196,39 @@ class MPS{
void full_state_vector(cvector_t &state_vector) const;
void get_probabilities_vector(rvector_t& probvector, const reg_t &qubits) const;

//methods from qasm_controller that are not supported yet
void set_omp_threads(int threads) {
static void set_omp_threads(uint_t threads) {
if (threads > 0)
omp_threads_ = threads;
}

void set_omp_threshold(int omp_qubit_threshold) {
static void set_omp_threshold(uint_t omp_qubit_threshold) {
if (omp_qubit_threshold > 0)
omp_threshold_ = omp_qubit_threshold;
}

void set_json_chop_threshold(double json_chop_threshold) {
static void set_json_chop_threshold(double json_chop_threshold) {
json_chop_threshold_ = json_chop_threshold;
}

void set_sample_measure_index_size(int index_size){
static void set_sample_measure_index_size(uint_t index_size){
sample_measure_index_size_ = index_size;
}
static void set_enable_gate_opt(bool enable_gate_opt) {
enable_gate_opt_ = enable_gate_opt;
}

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_;
}
static uint_t get_sample_measure_index_size(){
return sample_measure_index_size_;
}

void enable_gate_opt() {
std::cout << "enable_gate_opt not supported yet" <<std::endl;
static bool get_enable_gate_opt() {
return enable_gate_opt_;
}

// void store_measure(const AER::reg_t outcome, const AER::reg_t &cmemory, const AER::reg_t &cregister) const{
Expand Down Expand Up @@ -250,10 +262,12 @@ class MPS{

void initialize_from_statevector(uint_t num_qubits, cvector_t state_vector) {
cmatrix_t statevector_as_matrix(1, state_vector.size());
#pragma omp parallel for
for (int_t i=0; i<static_cast<int_t>(state_vector.size()); i++) {
statevector_as_matrix(0, i) = state_vector[i];
}

#pragma omp parallel for if (num_qubits_ > MPS::get_omp_threshold() && MPS::get_omp_threads() > 1) num_threads(MPS::get_omp_threads())
for (int_t i=0; i<static_cast<int_t>(state_vector.size()); i++) {
statevector_as_matrix(0, i) = state_vector[i];
}

initialize_from_matrix(num_qubits, statevector_as_matrix);
}
void initialize_from_matrix(uint_t num_qubits, cmatrix_t mat);
Expand Down Expand Up @@ -358,11 +372,12 @@ class MPS{
//-----------------------------------------------------------------------
// Config settings
//-----------------------------------------------------------------------
uint_t omp_threads_ = 1; // Disable multithreading by default
uint_t omp_threshold_ = 14; // Qubit threshold for multithreading when enabled
int sample_measure_index_size_ = 10; // Sample measure indexing qubit size
double json_chop_threshold_ = 1E-8; // Threshold for choping small values
static uint_t omp_threads_; // Disable multithreading by default
static uint_t omp_threshold_; // Qubit threshold for multithreading when enabled
static int sample_measure_index_size_; // Sample measure indexing qubit size
static double json_chop_threshold_; // Threshold for choping small values
// in JSON serialization
static bool enable_gate_opt_; // allow optimizations on gates
};

inline std::ostream &operator<<(std::ostream &out, const rvector_t &vec) {
Expand Down
Loading

0 comments on commit 7fb9147

Please sign in to comment.