Repository: rakytap/sequential-quantum-gate-decomposer (language: python)

Issue:
Fixes



Discussion / hints:

Fix patch:
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 123a6f60..e2886ae7 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -188,6 +188,7 @@ elseif (CMAKE_CXX_COMPILER_ID STREQUAL "GNU")
     list(APPEND CXX_FLAGS_RELEASE "-mavx512f" "-DUSE_AVX")
   elseif (${HAVE_AVX2_EXTENSIONS})
     list(APPEND CXX_FLAGS_RELEASE "-mavx2" "-DUSE_AVX")
+    list(APPEND EXTRA_INCLUDES "./gates/kernels/include/")
   elseif (${HAVE_AVX_EXTENSIONS})
     list(APPEND CXX_FLAGS_RELEASE "-mavx" "-DUSE_AVX")
   endif()
@@ -366,6 +367,7 @@ list(APPEND qgd_files
     ${PROJECT_SOURCE_DIR}/gates/RX.cpp
     ${PROJECT_SOURCE_DIR}/gates/RZ.cpp
     ${PROJECT_SOURCE_DIR}/gates/Composite.cpp
+    ${PROJECT_SOURCE_DIR}/gates/kernels/apply_kernel_to_input_AVX.cpp
     ${PROJECT_SOURCE_DIR}/nn/NN.cpp
     ${PROJECT_SOURCE_DIR}/decomposition/Decomposition_Base.cpp
     ${PROJECT_SOURCE_DIR}/decomposition/N_Qubit_Decomposition_Base.cpp
@@ -397,9 +399,6 @@ if (DEFINED ENV{QGD_DFE})
   list(APPEND CXX_FLAGS_DEBUG "-D__DFE__=1")
   list(APPEND CXX_FLAGS_RELEASE "-D__DFE__=1")
 
-  find_library(DFE_LIB qgdDFE
-               PATHS $ENV{DFE_LIB_DIR}
-               REQUIRED)
 
 endif()
 
@@ -425,7 +424,6 @@ target_link_libraries( qgd PRIVATE
     ${TBBMALLOC_PROXY_LIB}
     ${TBB_LIB}
     ${GSL_LIB}
-    ${DFE_LIB}
     ${MPI_C_LIBRARIES}
     )
 
@@ -437,6 +435,7 @@ target_include_directories(qgd PRIVATE
                             ./decomposition/include
                             ./random_unitary/include
                             ./nn/include
+                            ./gates/kernels/include
                             ${EXTRA_INCLUDES})
 
 
diff --git a/common/common_DFE.cpp b/common/common_DFE.cpp
index 3f9ee3a3..0b29c6dc 100644
--- a/common/common_DFE.cpp
+++ b/common/common_DFE.cpp
@@ -27,6 +27,36 @@ along with this program.  If not, see http://www.gnu.org/licenses/.
 
 #include "common_DFE.h"
 
+#include <atomic>
+#include <dlfcn.h>
+#include <unistd.h>
+#include <mutex>
+
+
+// pointer to the dynamically loaded DFE library
+void* handle = NULL;
+
+/// reference counting of locking-unlocking the DFE accelerators
+std::atomic_size_t read_count(0); //readers-writer problem semaphore
+
+/// mutex to guard DFE lib locking and unlocking
+std::recursive_mutex libmutex; //writing mutex
+std::mutex libreadmutex; //reader mutex
+
+
+size_t (*get_accelerator_avail_num_dll)() = NULL;
+size_t (*get_accelerator_free_num_dll)() = NULL;
+int (*calcqgdKernelDFE_dll)(size_t rows, size_t cols, DFEgate_kernel_type* gates, int gatesNum, int gateSetNum, double* trace) = NULL;
+int (*load2LMEM_dll)(QGD_Complex16* data, size_t rows, size_t cols) = NULL;
+void (*releive_DFE_dll)() = NULL;
+int (*initialize_DFE_dll)( int accelerator_num ) = NULL;
+int (*get_chained_gates_num_dll)() = NULL;
+
+// The ID of the class initializing the DFE lib
+int initialize_id = -1;
+
+
+
 /**
 @brief ????????????
 @return ??????????
@@ -36,7 +66,153 @@ void uploadMatrix2DFE( Matrix& input ) {
     std::cout << "size in bytes of uploading to DFE: " << input.size()*2*sizeof(float) << std::endl;    
 
     // load the data to LMEM
-    load2LMEM( input.get_data(), input.rows, input.cols );
+    load2LMEM_dll( input.get_data(), input.rows, input.cols );
+
+}
+
+
+/**
+@brief ????????????
+@return ??????????
+*/
+void unload_dfe_lib()
+{
+    const std::lock_guard<std::recursive_mutex> lock(libmutex);
+    if (handle) {
+        releive_DFE_dll();
+        dlclose(handle);
+        handle = NULL;
+    }
+}
+
+
+
+
+
+
+
+/**
+@brief ????????????
+@return ??????????
+*/
+int init_dfe_lib( const int accelerator_num, int qbit_num, int initialize_id_in )  {
+
+    const std::lock_guard<std::recursive_mutex> lock(libmutex);
+
+    initialize_id = initialize_id_in;
+
+    
+    unload_dfe_lib();
+
+
+    std::string lib_name_DFE = qbit_num > 9 ? DFE_LIB_10QUBITS : DFE_LIB_9QUBITS;
+    std::string lib_name     = getenv("SLIC_CONF") ? DFE_LIB_SIM : lib_name_DFE;
+
+    // dynamic-loading the correct DFE permanent calculator (Simulator/DFE/single or dual) from shared libararies
+    handle = dlopen(lib_name.c_str(), RTLD_NOW); //"MAXELEROSDIR"
+    if (handle == NULL && qbit_num == 10 && !getenv("SLIC_CONF")) {
+        handle = dlopen(DFE_LIB_9QUBITS, RTLD_NOW);
+    }
+    if (handle == NULL) {
+
+        std::string err("init_dfe_lib: failed to load library " + lib_name);
+        throw err;
+
+    } 
+    else {
+
+        get_accelerator_avail_num_dll = (size_t (*)())dlsym(handle, "get_accelerator_avail_num");
+        get_accelerator_free_num_dll  = (size_t (*)())dlsym(handle, "get_accelerator_free_num");
+        calcqgdKernelDFE_dll          = (int (*)(size_t, size_t, DFEgate_kernel_type*, int, int, double*))dlsym(handle, "calcqgdKernelDFE");
+        load2LMEM_dll                 = (int (*)(QGD_Complex16*, size_t, size_t))dlsym(handle, "load2LMEM");
+        releive_DFE_dll               = (void (*)())dlsym(handle, "releive_DFE");
+        initialize_DFE_dll            = (int (*)(int))dlsym(handle, "initialize_DFE");
+        get_chained_gates_num_dll     = (int (*)())dlsym(handle, "get_chained_gates_num");
+
+        return initialize_DFE_dll(accelerator_num);
+
+    }
+
+
+}
+
+
+/**
+@brief ????????????
+@return ??????????
+*/
+void lock_lib()
+{
+    const std::lock_guard<std::mutex> lock(libreadmutex);
+    if (++read_count == 1) libmutex.lock();
+}
+
+
+/**
+@brief ????????????
+@return ??????????
+*/
+void unlock_lib()
+{
+    const std::lock_guard<std::mutex> lock(libreadmutex);
+    if (--read_count == 0) libmutex.unlock();
+}
+
+
+
+
+
+/**
+@brief ????????????
+@return ??????????
+*/
+size_t get_accelerator_avail_num() {
+
+    return get_accelerator_avail_num_dll();
+
+}
+
+
+/**
+@brief ????????????
+@return ??????????
+*/
+size_t get_accelerator_free_num() {
+
+    return get_accelerator_free_num_dll();
+
+}
+
+/**
+@brief ????????????
+@return ??????????
+*/
+int get_initialize_id() {
+
+    return initialize_id;
+
+}
+
+/**
+@brief ????????????
+@return ??????????
+*/
+int calcqgdKernelDFE(size_t rows, size_t cols, DFEgate_kernel_type* gates, int gatesNum, int gateSetNum, double* trace) {
+
+
+    return calcqgdKernelDFE_dll(rows, cols, gates, gatesNum, gateSetNum, trace);
+
+}
+
+
+
+/**
+@brief ????????????
+@return ??????????
+*/
+int get_chained_gates_num() {
+
+    return get_chained_gates_num_dll();
 
 }
 
diff --git a/common/include/common_DFE.h b/common/include/common_DFE.h
index 34d13b4e..1c5ffda8 100644
--- a/common/include/common_DFE.h
+++ b/common/include/common_DFE.h
@@ -38,11 +38,16 @@ along with this program.  If not, see http://www.gnu.org/licenses/.
 #include <cstring>
 #include <sstream>
 
-extern "C"
-{
+
+#define DFE_LIB_9QUBITS "libqgdDFE.so"
+#define DFE_LIB_10QUBITS "libqgdDFE_10qubits.so"
+#define DFE_LIB_SIM "libqgdDFE_SIM.so"
 
 
 
+extern "C"
+{
+
 
 /**
 @brief Fixed point data related to a gate operation
@@ -65,29 +70,43 @@ typedef struct {
 } DFEgate_kernel_type;
 
 
+}
+
+
+
 /**
 @brief ????????????
 @return ??????????
 */
-int calcqgdKernelDFE(size_t rows, size_t cols, DFEgate_kernel_type* gates, int gatesNum, int gateSetNum, double* trace);
+size_t get_accelerator_avail_num();
+
 
 /**
 @brief ????????????
 @return ??????????
 */
-int load2LMEM( QGD_Complex16* data, size_t rows, size_t cols );
+size_t get_accelerator_free_num();
 
 /**
 @brief ????????????
 @return ??????????
 */
-void releive_DFE();
+int get_initialize_id();
 
 /**
 @brief ????????????
 @return ??????????
 */
-int initialize_DFE();
+int calcqgdKernelDFE(size_t rows, size_t cols, DFEgate_kernel_type* gates, int gatesNum, int gateSetNum, double* trace);
+
+
+/**
+@brief ????????????
+@return ??????????
+*/
+int get_chained_gates_num();
+
+
 
 
 /**
@@ -97,13 +116,36 @@ int initialize_DFE();
 void uploadMatrix2DFE( Matrix& input );
 
 
+
 /**
 @brief ????????????
 @return ??????????
 */
-int get_chained_gates_num();
+void unload_dfe_lib();
+
+
+/**
+@brief ????????????
+@return ??????????
+*/
+void lock_lib();
+
+
+/**
+@brief ????????????
+@return ??????????
+*/
+void unlock_lib();
+
+
+
+
+/**
+@brief ????????????
+@return ??????????
+*/
+int init_dfe_lib( const int accelerator_num, int qbit_num, int initialize_id_in);
 
-}
 
 
 #endif
diff --git a/decomposition/Decomposition_Base.cpp b/decomposition/Decomposition_Base.cpp
index 16085ed0..9d11f1fc 100644
--- a/decomposition/Decomposition_Base.cpp
+++ b/decomposition/Decomposition_Base.cpp
@@ -1339,17 +1339,17 @@ std::vector<Gate*> Decomposition_Base::prepare_gates_to_export( std::vector<Gate
         else if (gate->get_type() == CRY_OPERATION) {
 
             // definig the parameter of the rotational angle
-            double phi;
+            double Theta;
 
             // get the inverse parameters of the RZ rotation
 
             CRY* cry_gate = static_cast<CRY*>(gate);
 
-            phi = std::fmod( parameters[parameter_idx], 2*M_PI);
+            Theta = std::fmod( 2.0*parameters[parameter_idx], 4*M_PI);
             parameter_idx = parameter_idx + 1;
 
 
-            cry_gate->set_optimized_parameters( phi );
+            cry_gate->set_optimized_parameters( Theta );
             ops_ret.push_back( static_cast<Gate*>(cry_gate) );
 
 
@@ -1417,17 +1417,17 @@ std::vector<Gate*> Decomposition_Base::prepare_gates_to_export( std::vector<Gate
         else if (gate->get_type() == ADAPTIVE_OPERATION) {
 
             // definig the parameter of the rotational angle
-            double phi;
+            double Theta;
 
             // get the inverse parameters of the RZ rotation
 
             Adaptive* ad_gate = static_cast<Adaptive*>(gate);
 
-            phi = std::fmod( parameters[parameter_idx], 2*M_PI);
+            Theta = std::fmod( 2*parameters[parameter_idx], 4*M_PI);
             parameter_idx = parameter_idx + 1;
 
 
-            ad_gate->set_optimized_parameters( phi );
+            ad_gate->set_optimized_parameters( Theta );
             ops_ret.push_back( static_cast<Gate*>(ad_gate) );
 
 
diff --git a/decomposition/N_Qubit_Decomposition_Base.cpp b/decomposition/N_Qubit_Decomposition_Base.cpp
index 8233d95a..54764527 100644
--- a/decomposition/N_Qubit_Decomposition_Base.cpp
+++ b/decomposition/N_Qubit_Decomposition_Base.cpp
@@ -72,11 +72,13 @@ N_Qubit_Decomposition_Base::N_Qubit_Decomposition_Base() {
 
     iteration_threshold_of_randomization = 2500000;
 
-#ifdef __DFE__
-    if( qbit_num >= 5 ) {
-        uploadMatrix2DFE( Umtx );
-    }
-#endif
+    // number of utilized accelerators
+    accelerator_num = 0;
+
+    // unique id indentifying the instance of the class
+    std::uniform_int_distribution<> distrib_int(0, INT_MAX);  
+    int id = distrib_int(gen);
+
 }
 
 /**
@@ -87,7 +89,7 @@ N_Qubit_Decomposition_Base::N_Qubit_Decomposition_Base() {
 @param initial_guess_in Enumeration element indicating the method to guess initial values for the optimization. Possible values: 'zeros=0' ,'random=1', 'close_to_zero=2'
 @return An instance of the class
 */
-N_Qubit_Decomposition_Base::N_Qubit_Decomposition_Base( Matrix Umtx_in, int qbit_num_in, bool optimize_layer_num_in, guess_type initial_guess_in= CLOSE_TO_ZERO ) : Decomposition_Base(Umtx_in, qbit_num_in, initial_guess_in) {
+N_Qubit_Decomposition_Base::N_Qubit_Decomposition_Base( Matrix Umtx_in, int qbit_num_in, bool optimize_layer_num_in, guess_type initial_guess_in= CLOSE_TO_ZERO, int accelerator_num_in ) : Decomposition_Base(Umtx_in, qbit_num_in, initial_guess_in) {
 
     // logical value. Set true if finding the minimum number of gate layers is required (default), or false when the maximal number of two-qubit gates is used (ideal for general unitaries).
     optimize_layer_num  = optimize_layer_num_in;
@@ -125,10 +127,18 @@ N_Qubit_Decomposition_Base::N_Qubit_Decomposition_Base( Matrix Umtx_in, int qbit
 
     iteration_threshold_of_randomization = 2500000;
 
+
+    // unique id indentifying the instance of the class
+    std::uniform_int_distribution<> distrib_int(0, INT_MAX);  
+    id = distrib_int(gen);
+
+
 #ifdef __DFE__
-    if( qbit_num >= 5 ) {
-        uploadMatrix2DFE( Umtx );
-    }
+
+    // number of utilized accelerators
+    accelerator_num = accelerator_num_in;
+#else
+    accelerator_num = 0;
 #endif
 
 }
@@ -142,7 +152,7 @@ N_Qubit_Decomposition_Base::~N_Qubit_Decomposition_Base() {
 
 
 #ifdef __DFE__
-    releive_DFE();
+    unload_dfe_lib();//releive_DFE();
 #endif      
 
 
@@ -304,6 +314,14 @@ void N_Qubit_Decomposition_Base::solve_layer_optimization_problem( int num_of_pa
 */
 void N_Qubit_Decomposition_Base::solve_layer_optimization_problem_ADAM( int num_of_parameters, gsl_vector *solution_guess_gsl) {
 
+#ifdef __DFE__
+        if ( qbit_num >= 5 ) {
+            upload_Umtx_to_DFE();
+        }
+#endif
+
+
+
         if (gates.size() == 0 ) {
             return;
         }
@@ -504,6 +522,12 @@ pure_DFE_time = 0.0;
 void N_Qubit_Decomposition_Base::solve_layer_optimization_problem_BFGS( int num_of_parameters, gsl_vector *solution_guess_gsl) {
 
 
+#ifdef __DFE__
+        if ( qbit_num >= 5 ) {
+            upload_Umtx_to_DFE();
+        }
+#endif
+
         if (gates.size() == 0 ) {
             return;
         }
@@ -609,6 +633,13 @@ void N_Qubit_Decomposition_Base::solve_layer_optimization_problem_BFGS( int num_
 void N_Qubit_Decomposition_Base::solve_layer_optimization_problem_BFGS2( int num_of_parameters, gsl_vector *solution_guess_gsl) {
 
 
+#ifdef __DFE__
+        if ( qbit_num >= 5 ) {
+            upload_Umtx_to_DFE();
+        }
+#endif
+
+
         if (gates.size() == 0 ) {
             return;
         }
@@ -1008,28 +1039,11 @@ void N_Qubit_Decomposition_Base::optimization_problem_grad( const gsl_vector* pa
 @param void_instance A void pointer pointing to the instance of the current class.
 @param f0 The value of the cost function at x0.
 @param grad A GNU Scientific Library vector containing the calculated gradient components.
-@param onlyCPU Set true to use only CPU in the calculations (has effect if compiled to use accelerator devices)
 */
 void N_Qubit_Decomposition_Base::optimization_problem_combined( const gsl_vector* parameters, void* void_instance, double* f0, gsl_vector* grad ) {
 
     N_Qubit_Decomposition_Base* instance = reinterpret_cast<N_Qubit_Decomposition_Base*>(void_instance);
 
-    instance->optimization_problem_combined( parameters, void_instance, f0, grad, false );
-
-}
-
-/**
-@brief Call to calculate both the cost function and the its gradient components.
-@param parameters A GNU Scientific Library vector containing the free parameters to be optimized.
-@param void_instance A void pointer pointing to the instance of the current class.
-@param f0 The value of the cost function at x0.
-@param grad A GNU Scientific Library vector containing the calculated gradient components.
-@param onlyCPU Set true to use only CPU in the calculations (has effect if compiled to use accelerator devices)
-*/
-void N_Qubit_Decomposition_Base::optimization_problem_combined( const gsl_vector* parameters, void* void_instance, double* f0, gsl_vector* grad, bool onlyCPU ) {
-
-    N_Qubit_Decomposition_Base* instance = reinterpret_cast<N_Qubit_Decomposition_Base*>(void_instance);
-
     // the number of free parameters
     int parameter_num_loc = instance->get_parameter_num();
 
@@ -1048,7 +1062,7 @@ void N_Qubit_Decomposition_Base::optimization_problem_combined( const gsl_vector
 ///////////////////////////////////////
 //std::cout << "number of qubits: " << instance->qbit_num << std::endl;
 //tbb::tick_count t0_DFE = tbb::tick_count::now();/////////////////////////////////    
-if ( instance->qbit_num >= 5 && !onlyCPU ) {
+if ( instance->qbit_num >= 5 && instance->get_accelerator_num() > 0 ) {
     Matrix_real parameters_mtx(parameters->data, 1, parameters->size);
 
     int gatesNum, redundantGateSets, gateSetNum;
@@ -1066,14 +1080,18 @@ if ( instance->qbit_num >= 5 && !onlyCPU ) {
 
     Matrix_real mpi_trace_DFE_mtx(mpi_gateSetNum, 3);
 
+    lock_lib();
     calcqgdKernelDFE( Umtx_loc.rows, Umtx_loc.cols, DFEgates+mpi_starting_gateSetIdx*gatesNum, gatesNum, mpi_gateSetNum, mpi_trace_DFE_mtx.get_data() );
+    unlock_lib();
 
     int bytes = mpi_trace_DFE_mtx.size()*sizeof(double);
     MPI_Allgather(mpi_trace_DFE_mtx.get_data(), bytes, MPI_BYTE, trace_DFE_mtx.get_data(), bytes, MPI_BYTE, MPI_COMM_WORLD);
 
 #else
 
+    lock_lib();
     calcqgdKernelDFE( Umtx_loc.rows, Umtx_loc.cols, DFEgates, gatesNum, gateSetNum, trace_DFE_mtx.get_data() );
+    unlock_lib();
 
 #endif  
 
@@ -1216,9 +1234,19 @@ std::string error("N_Qubit_Decomposition_Base::optimization_problem_combined");
 @param parameters The parameters for which the cost fuction shoule be calculated
 @param f0 The value of the cost function at x0.
 @param grad An array storing the calculated gradient components
-@param onlyCPU Set true to use only CPU in the calculations (has effect if compiled to use accelerator devices)
 */
-void N_Qubit_Decomposition_Base::optimization_problem_combined( const Matrix_real& parameters, double* f0, Matrix_real& grad, bool onlyCPU ) {
+void N_Qubit_Decomposition_Base::optimization_problem_combined( const Matrix_real& parameters, double* f0, Matrix_real& grad ) {
+
+#ifdef __DFE__
+
+    lock_lib();
+
+    if ( get_initialize_id() != id ) {
+        std::string err("The uploaded unitary to the DFE might not be identical to the unitary stored by this specific class instance. Please upload the unitary to DFE by the Upload_Umtx_to_DFE() method.");
+        throw err;
+    }
+
+#endif
 
     // create GSL wrappers around the pointers
     gsl_block block_tmp;
@@ -1245,7 +1273,11 @@ void N_Qubit_Decomposition_Base::optimization_problem_combined( const Matrix_rea
     grad_gsl.owner = 0;    
 
     // call the method to calculate the cost function and the gradients
-    optimization_problem_combined( &parameters_gsl, this, f0, &grad_gsl, onlyCPU );
+    optimization_problem_combined( &parameters_gsl, this, f0, &grad_gsl );
+
+#ifdef __DFE__
+    unlock_lib();
+#endif
 
 }
 
@@ -1450,3 +1482,37 @@ N_Qubit_Decomposition_Base::set_iteration_threshold_of_randomization( const unsi
 
 }
 
+
+#ifdef __DFE__
+
+void 
+N_Qubit_Decomposition_Base::upload_Umtx_to_DFE() {
+
+    lock_lib();
+
+    if ( get_initialize_id() != id ) {
+        // initialize DFE library
+        init_dfe_lib( accelerator_num, qbit_num, id );
+    }
+
+    uploadMatrix2DFE( Umtx );
+
+
+    unlock_lib();
+
+}
+
+
+/**
+@brief Get the number of accelerators to be reserved on DFEs on users demand. 
+*/
+int 
+N_Qubit_Decomposition_Base::get_accelerator_num() {
+
+    return accelerator_num;
+
+}
+
+
+#endif
+
diff --git a/decomposition/N_Qubit_Decomposition_adaptive.cpp b/decomposition/N_Qubit_Decomposition_adaptive.cpp
index 4d693312..7e4fbbff 100644
--- a/decomposition/N_Qubit_Decomposition_adaptive.cpp
+++ b/decomposition/N_Qubit_Decomposition_adaptive.cpp
@@ -83,7 +83,7 @@ N_Qubit_Decomposition_adaptive::N_Qubit_Decomposition_adaptive() : N_Qubit_Decom
 @param initial_guess_in Enumeration element indicating the method to guess initial values for the optimization. Possible values: 'zeros=0' ,'random=1', 'close_to_zero=2'
 @return An instance of the class
 */
-N_Qubit_Decomposition_adaptive::N_Qubit_Decomposition_adaptive( Matrix Umtx_in, int qbit_num_in, int level_limit_in, int level_limit_min_in ) : N_Qubit_Decomposition_Base(Umtx_in, qbit_num_in, false, RANDOM) {
+N_Qubit_Decomposition_adaptive::N_Qubit_Decomposition_adaptive( Matrix Umtx_in, int qbit_num_in, int level_limit_in, int level_limit_min_in, int accelerator_num ) : N_Qubit_Decomposition_Base(Umtx_in, qbit_num_in, false, RANDOM, accelerator_num) {
 
 
     // set the level limit
@@ -128,7 +128,7 @@ N_Qubit_Decomposition_adaptive::N_Qubit_Decomposition_adaptive( Matrix Umtx_in,
 @param initial_guess_in Enumeration element indicating the method to guess initial values for the optimization. Possible values: 'zeros=0' ,'random=1', 'close_to_zero=2'
 @return An instance of the class
 */
-N_Qubit_Decomposition_adaptive::N_Qubit_Decomposition_adaptive( Matrix Umtx_in, int qbit_num_in, int level_limit_in, int level_limit_min_in, std::vector<matrix_base<int>> topology_in ) : N_Qubit_Decomposition_Base(Umtx_in, qbit_num_in, false, RANDOM) {
+N_Qubit_Decomposition_adaptive::N_Qubit_Decomposition_adaptive( Matrix Umtx_in, int qbit_num_in, int level_limit_in, int level_limit_min_in, std::vector<matrix_base<int>> topology_in, int accelerator_num ) : N_Qubit_Decomposition_Base(Umtx_in, qbit_num_in, false, RANDOM, accelerator_num) {
 
 
 
@@ -218,10 +218,6 @@ N_Qubit_Decomposition_adaptive::start_decomposition(bool prepare_export) {
     }
 
 
-#ifdef __DFE__
-    uploadMatrix2DFE( Umtx );
-#endif
-
 
     double optimization_tolerance_orig = optimization_tolerance;
 
@@ -421,7 +417,7 @@ N_Qubit_Decomposition_adaptive::optimize_imported_gate_structure(Matrix_real& op
     // solve the optimization problem
     N_Qubit_Decomposition_custom cDecomp_custom;
     // solve the optimization problem in isolated optimization process
-    cDecomp_custom = N_Qubit_Decomposition_custom( Umtx.copy(), qbit_num, false, initial_guess);
+    cDecomp_custom = N_Qubit_Decomposition_custom( Umtx.copy(), qbit_num, false, initial_guess, accelerator_num);
     cDecomp_custom.set_custom_gate_structure( gate_structure_loc );
     cDecomp_custom.set_optimized_parameters( optimized_parameters_mtx_loc.get_data(), optimized_parameters_mtx_loc.size() );
     cDecomp_custom.set_optimization_blocks( gate_structure_loc->get_gate_num() );
@@ -525,7 +521,7 @@ N_Qubit_Decomposition_adaptive::determine_initial_gate_structure(Matrix_real& op
             [&]{            
 #endif
                 // solve the optimization problem in isolated optimization process
-                cDecomp_custom_random = N_Qubit_Decomposition_custom( Umtx.copy(), qbit_num, false, RANDOM);
+                cDecomp_custom_random = N_Qubit_Decomposition_custom( Umtx.copy(), qbit_num, false, RANDOM, accelerator_num);
                 cDecomp_custom_random.set_custom_gate_structure( gate_structure_loc );
                 cDecomp_custom_random.set_optimization_blocks( gate_structure_loc->get_gate_num() );
                 cDecomp_custom_random.set_max_iteration( max_iterations );
@@ -847,7 +843,7 @@ N_Qubit_Decomposition_adaptive::compress_gate_structure( Gates_block* gate_struc
     N_Qubit_Decomposition_custom cDecomp_custom;
        
     // solve the optimization problem in isolated optimization process
-    cDecomp_custom = N_Qubit_Decomposition_custom( Umtx.copy(), qbit_num, false, initial_guess);
+    cDecomp_custom = N_Qubit_Decomposition_custom( Umtx.copy(), qbit_num, false, initial_guess, accelerator_num);
     cDecomp_custom.set_custom_gate_structure( gate_structure_reduced );
     cDecomp_custom.set_optimized_parameters( parameters_reduced.get_data(), parameters_reduced.size() );
     cDecomp_custom.set_verbose(0);
@@ -1490,7 +1486,7 @@ N_Qubit_Decomposition_adaptive::set_unitary_from_file( std::string filename ) {
 
 #ifdef __DFE__
     if( qbit_num >= 5 ) {
-        uploadMatrix2DFE( Umtx );
+        upload_Umtx_to_DFE();
     }
 #endif
 
@@ -1506,7 +1502,7 @@ N_Qubit_Decomposition_adaptive::set_unitary( Matrix& Umtx_new ) {
 
 #ifdef __DFE__
     if( qbit_num >= 5 ) {
-        uploadMatrix2DFE( Umtx );
+        upload_Umtx_to_DFE();
     }
 #endif
 
diff --git a/decomposition/N_Qubit_Decomposition_custom.cpp b/decomposition/N_Qubit_Decomposition_custom.cpp
index 2ce98071..27d40c5c 100644
--- a/decomposition/N_Qubit_Decomposition_custom.cpp
+++ b/decomposition/N_Qubit_Decomposition_custom.cpp
@@ -61,7 +61,7 @@ N_Qubit_Decomposition_custom::N_Qubit_Decomposition_custom() : N_Qubit_Decomposi
 @param initial_guess_in Enumeration element indicating the method to guess initial values for the optimization. Possible values: 'zeros=0' ,'random=1', 'close_to_zero=2'
 @return An instance of the class
 */
-N_Qubit_Decomposition_custom::N_Qubit_Decomposition_custom( Matrix Umtx_in, int qbit_num_in, bool optimize_layer_num_in, guess_type initial_guess_in= CLOSE_TO_ZERO ) : N_Qubit_Decomposition_Base(Umtx_in, qbit_num_in, optimize_layer_num_in, initial_guess_in) {
+N_Qubit_Decomposition_custom::N_Qubit_Decomposition_custom( Matrix Umtx_in, int qbit_num_in, bool optimize_layer_num_in, guess_type initial_guess_in= CLOSE_TO_ZERO, int accelerator_num ) : N_Qubit_Decomposition_Base(Umtx_in, qbit_num_in, optimize_layer_num_in, initial_guess_in, accelerator_num) {
 
     // initialize custom gate structure
     gate_structure = NULL;
@@ -134,9 +134,6 @@ N_Qubit_Decomposition_custom::start_decomposition(bool prepare_export) {
     // setting the gate structure for optimization
     add_gate_layers();
 
-#ifdef __DFE__
-    uploadMatrix2DFE( Umtx );
-#endif
 
 /*
 if (optimized_parameters_mtx.size() > 0 ) {
@@ -148,9 +145,7 @@ std::cout << "ooooooooooooo " <<  optimized_parameters_mtx.size() << std::endl;
     // final tuning of the decomposition parameters
     final_optimization();
 
-#ifdef __DFE__
-//return;
-#endif
+
     // prepare gates to export
     if (prepare_export) {
         prepare_gates_to_export();
diff --git a/decomposition/include/N_Qubit_Decomposition_Base.h b/decomposition/include/N_Qubit_Decomposition_Base.h
index 7d4e3332..08f9aa11 100644
--- a/decomposition/include/N_Qubit_Decomposition_Base.h
+++ b/decomposition/include/N_Qubit_Decomposition_Base.h
@@ -79,7 +79,8 @@ class N_Qubit_Decomposition_Base : public Decomposition_Base {
     int random_shift_count_max;
     ///    
     double eta;
-
+    /// unique id indentifying the instance of the class
+    int id;
 protected:
 
 
@@ -107,6 +108,10 @@ class N_Qubit_Decomposition_Base : public Decomposition_Base {
     double randomization_rate;
     /// threashold of count of iterations after what the parameters are randomized if the cost function does not deacrese fast enough
     unsigned long long iteration_threshold_of_randomization;
+    /// number of utilized accelerators
+    int accelerator_num;
+
+
 
     Matrix_real randomization_probs;
     matrix_base<int> randomized_probs;
@@ -134,7 +139,7 @@ N_Qubit_Decomposition_Base();
 @param initial_guess_in Enumeration element indicating the method to guess initial values for the optimization. Possible values: 'zeros=0' ,'random=1', 'close_to_zero=2'
 @return An instance of the class
 */
-N_Qubit_Decomposition_Base( Matrix Umtx_in, int qbit_num_in, bool optimize_layer_num_in, guess_type initial_guess_in );
+N_Qubit_Decomposition_Base( Matrix Umtx_in, int qbit_num_in, bool optimize_layer_num_in, guess_type initial_guess_in, int accelerator_num_in=0 );
 
 
 
@@ -247,18 +252,6 @@ static void optimization_problem_grad( const gsl_vector* parameters, void* void_
 static void optimization_problem_combined( const gsl_vector* parameters, void* void_instance, double* f0, gsl_vector* grad );
 
 
-/**
-@brief Call to calculate both the cost function and the its gradient components.
-@param parameters A GNU Scientific Library vector containing the free parameters to be optimized.
-@param void_instance A void pointer pointing to the instance of the current class.
-@param f0 The value of the cost function at x0.
-@param grad A GNU Scientific Library vector containing the calculated gradient components.
-@param onlyCPU Set true to use only CPU in the calculations (has effect if compiled to use accelerator devices)
-*/
-static void optimization_problem_combined( const gsl_vector* parameters, void* void_instance, double* f0, gsl_vector* grad, bool onlyCPU  );
-
-
-
 /**
 @brief Call to calculate both the cost function and the its gradient components.
 @param parameters The parameters for which the cost fuction shoule be calculated
@@ -266,7 +259,7 @@ static void optimization_problem_combined( const gsl_vector* parameters, void* v
 @param grad An array storing the calculated gradient components
 @param onlyCPU Set true to use only CPU in the calculations (has effect if compiled to use accelerator devices)
 */
-void optimization_problem_combined( const Matrix_real& parameters, double* f0, Matrix_real& grad, bool onlyCPU );
+void optimization_problem_combined( const Matrix_real& parameters, double* f0, Matrix_real& grad );
 
 
 /**
@@ -367,6 +360,22 @@ void set_adaptive_eta( bool adaptive_eta_in  );
 */
 void set_randomized_radius( double radius_in  );
 
+#ifdef __DFE__
+
+/**
+@brief ?????????????
+*/
+void upload_Umtx_to_DFE();
+
+
+/**
+@brief Get the number of accelerators to be reserved on DFEs on users demand.
+*/
+int get_accelerator_num();
+
+
+#endif
+
 };
 
 
diff --git a/decomposition/include/N_Qubit_Decomposition_adaptive.h b/decomposition/include/N_Qubit_Decomposition_adaptive.h
index c078f342..ae7511a7 100644
--- a/decomposition/include/N_Qubit_Decomposition_adaptive.h
+++ b/decomposition/include/N_Qubit_Decomposition_adaptive.h
@@ -93,7 +93,7 @@ N_Qubit_Decomposition_adaptive();
 @param initial_guess_in Enumeration element indicating the method to guess initial values for the optimization. Possible values: 'zeros=0' ,'random=1', 'close_to_zero=2'
 @return An instance of the class
 */
-N_Qubit_Decomposition_adaptive( Matrix Umtx_in, int qbit_num_in, int level_limit_in, int level_limit_min_in );
+N_Qubit_Decomposition_adaptive( Matrix Umtx_in, int qbit_num_in, int level_limit_in, int level_limit_min_in, int accelerator_num=0 );
 
 
 /**
@@ -104,7 +104,7 @@ N_Qubit_Decomposition_adaptive( Matrix Umtx_in, int qbit_num_in, int level_limit
 @param initial_guess_in Enumeration element indicating the method to guess initial values for the optimization. Possible values: 'zeros=0' ,'random=1', 'close_to_zero=2'
 @return An instance of the class
 */
-N_Qubit_Decomposition_adaptive( Matrix Umtx_in, int qbit_num_in, int level_limit_in, int level_limit_min_in, std::vector<matrix_base<int>> topology_in );
+N_Qubit_Decomposition_adaptive( Matrix Umtx_in, int qbit_num_in, int level_limit_in, int level_limit_min_in, std::vector<matrix_base<int>> topology_in, int accelerator_num=0 );
 
 
 
diff --git a/decomposition/include/N_Qubit_Decomposition_custom.h b/decomposition/include/N_Qubit_Decomposition_custom.h
index a0c31745..09f279b5 100644
--- a/decomposition/include/N_Qubit_Decomposition_custom.h
+++ b/decomposition/include/N_Qubit_Decomposition_custom.h
@@ -60,7 +60,7 @@ N_Qubit_Decomposition_custom();
 @param initial_guess_in Enumeration element indicating the method to guess initial values for the optimization. Possible values: 'zeros=0' ,'random=1', 'close_to_zero=2'
 @return An instance of the class
 */
-N_Qubit_Decomposition_custom( Matrix Umtx_in, int qbit_num_in, bool optimize_layer_num_in, guess_type initial_guess_in );
+N_Qubit_Decomposition_custom( Matrix Umtx_in, int qbit_num_in, bool optimize_layer_num_in, guess_type initial_guess_in, int accelerator_num=0 );
 
 
 
diff --git a/examples/vqe/vqe_utils.py b/examples/vqe/vqe_utils.py
index 0a902cf4..92a37105 100644
--- a/examples/vqe/vqe_utils.py
+++ b/examples/vqe/vqe_utils.py
@@ -80,8 +80,8 @@ def get_optimized_circuit( alpha, optimized_parameters_in=None ):
 		# importing the quantum circuit
 		cDecompose.import_Qiskit_Circuit(qc_trial)
 
-		if isinstance(optimized_parameters_in, (np.ndarray, np.generic) ) :
-			cDecompose.set_Optimized_Parameters( optimized_parameters_in )
+		#if isinstance(optimized_parameters_in, (np.ndarray, np.generic) ) :
+		#	cDecompose.set_Optimized_Parameters( optimized_parameters_in )
 
 		# set the number of successive identical blocks in the optimalization of disentanglement of the n-th qubits
 		cDecompose.set_Optimization_Blocks( 200 )
diff --git a/gates/Gate.cpp b/gates/Gate.cpp
index a4d08b05..ac7b1282 100644
--- a/gates/Gate.cpp
+++ b/gates/Gate.cpp
@@ -25,6 +25,9 @@ along with this program.  If not, see http://www.gnu.org/licenses/.
 #include "Gate.h"
 #include "common.h"
 
+#ifdef USE_AVX 
+#include "apply_kernel_to_input_AVX.h"
+#endif
 
 
 /**
@@ -253,12 +256,17 @@ void
 Gate::apply_kernel_to(Matrix& u3_1qbit, Matrix& input, bool deriv) {
 
 
+#ifdef USE_AVX
+
+    apply_kernel_to_input_AVX(u3_1qbit, input, deriv, target_qbit, control_qbit, matrix_size);
+    return;
+
+#else
     
     int index_step_target = 1 << target_qbit;
     int current_idx = 0;
     int current_idx_pair = current_idx+index_step_target;
 
-//std::cout << "target qbit: " << target_qbit << std::endl;
 
     while ( current_idx_pair < matrix_size ) {
 
@@ -295,6 +303,7 @@ Gate::apply_kernel_to(Matrix& u3_1qbit, Matrix& input, bool deriv) {
 
                 }
 
+
             }
             else if (deriv) {
                 // when calculating derivatives, the constant element should be zeros
@@ -307,8 +316,6 @@ Gate::apply_kernel_to(Matrix& u3_1qbit, Matrix& input, bool deriv) {
             }
 
 
-//std::cout << current_idx_target << " " << current_idx_target_pair << std::endl;
-
         
         //});
         }
@@ -321,6 +328,9 @@ Gate::apply_kernel_to(Matrix& u3_1qbit, Matrix& input, bool deriv) {
     }
 
 
+#endif // USE_AVX
+
+
 }
 
 
diff --git a/gates/kernels/apply_kernel_to_input_AVX.cpp b/gates/kernels/apply_kernel_to_input_AVX.cpp
new file mode 100644
index 00000000..2f7cc5ca
--- /dev/null
+++ b/gates/kernels/apply_kernel_to_input_AVX.cpp
@@ -0,0 +1,213 @@
+/*
+Created on Fri Jun 26 14:13:26 2020
+Copyright (C) 2020 Peter Rakyta, Ph.D.
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this program.  If not, see http://www.gnu.org/licenses/.
+
+@author: Peter Rakyta, Ph.D.
+*/
+/*! \file apply_kerel_to_input_AVX.cpp
+    \brief ????????????????
+*/
+
+
+#include "apply_kernel_to_input_AVX.h"
+#include <immintrin.h>
+
+
+/**
+@brief AVX kernel to apply single qubit gate kernel on an input matrix
+@param A matrix on which the householder transformation is applied. (The output is returned via this matrix)
+@param v A matrix instance of the reflection vector
+*/
+void
+apply_kernel_to_input_AVX(Matrix& u3_1qbit, Matrix& input, const bool& deriv, const int& target_qbit, const int& control_qbit, const int& matrix_size) {
+
+
+    int index_step_target = 1 << target_qbit;
+    int current_idx = 0;
+    int current_idx_pair = current_idx+index_step_target;
+
+
+    // load elements of the U3 unitary into 256bit registers (4 registers)
+    __m128d* u3_1qubit_tmp  = (__m128d*)&u3_1qbit[0];
+    __m256d u3_1qbit_00_vec = _mm256_broadcast_pd( u3_1qubit_tmp );
+
+    u3_1qubit_tmp = (__m128d*)&u3_1qbit[1];
+    __m256d u3_1qbit_01_vec = _mm256_broadcast_pd( u3_1qubit_tmp );
+
+    u3_1qubit_tmp = (__m128d*)&u3_1qbit[2];
+    __m256d u3_1qbit_10_vec = _mm256_broadcast_pd( u3_1qubit_tmp );
+
+    u3_1qubit_tmp = (__m128d*)&u3_1qbit[3];
+    __m256d u3_1qbit_11_vec = _mm256_broadcast_pd( u3_1qubit_tmp );
+
+
+
+    while ( current_idx_pair < matrix_size ) {
+
+        for(int idx=0; idx<index_step_target; idx++) {  
+        //tbb::parallel_for(0, index_step_target, 1, [&](int idx) {  
+
+            int current_idx_loc = current_idx + idx;
+            int current_idx_pair_loc = current_idx_pair + idx;
+
+            int row_offset = current_idx_loc*input.stride;
+            int row_offset_pair = current_idx_pair_loc*input.stride;
+
+           if ( control_qbit<0 || ((current_idx_loc >> control_qbit) & 1) ) {
+
+
+                double* element      = (double*)input.get_data() + 2*row_offset;
+                double* element_pair = (double*)input.get_data() + 2*row_offset_pair;
+
+
+                __m256d neg = _mm256_setr_pd(1.0, -1.0, 1.0, -1.0); // 5th register
+
+
+                for ( int col_idx=0; col_idx<2*(input.cols-1); col_idx=col_idx+4) {
+
+                    // extract successive elements from arrays element, element_pair
+                    __m256d element_vec      = _mm256_loadu_pd( element+col_idx ); // 6th register
+                    __m256d element_pair_vec = _mm256_loadu_pd( element_pair+col_idx ); // 7th register
+
+                    //// u3_1qbit_00*element_vec ////
+
+                    // 1 calculate the multiplications  u3_1qbit_00*element_vec
+                    __m256d vec3 = _mm256_mul_pd(u3_1qbit_00_vec, element_vec); // 8th register
+
+                    // 2 Switch the real and imaginary elements of element_vec
+                    __m256d element_vec_permuted = _mm256_permute_pd(element_vec, 0x5);   // 9th register
+
+                    // 3 Negate the imaginary elements of element_vec_permuted
+                    element_vec_permuted = _mm256_mul_pd(element_vec_permuted, neg);
+
+                    // 4 Multiply elements of u3_1qbit_00*element_vec_permuted
+                    __m256d vec4 = _mm256_mul_pd(u3_1qbit_00_vec, element_vec_permuted); 
+
+                    // 5 Horizontally subtract the elements in vec3 and vec4
+                    vec3  = _mm256_hsub_pd(vec3, vec4); 
+
+
+                    //// u3_1qbit_01*element_vec_pair ////
+
+                    // 1 calculate the multiplications  u3_1qbit_01*element_pair_vec
+                    __m256d vec5 = _mm256_mul_pd(u3_1qbit_01_vec, element_pair_vec); // 10th register
+
+                    // 2 Switch the real and imaginary elements of element_vec
+                    __m256d element_pair_vec_permuted = _mm256_permute_pd(element_pair_vec, 0x5);   // 11th register
+
+                    // 3 Negate the imaginary elements of element_vec_permuted
+                    element_pair_vec_permuted = _mm256_mul_pd(element_pair_vec_permuted, neg);
+
+                    // 4 Multiply elements of u3_1qbit_01*element_vec_pair_permuted
+                    vec4 = _mm256_mul_pd(u3_1qbit_01_vec, element_pair_vec_permuted); 
+
+                    // 5 Horizontally subtract the elements in vec5 and vec4
+                    vec5  = _mm256_hsub_pd(vec5, vec4); 
+
+                    //// u3_1qbit_00*element_vec + u3_1qbit_01*element_vec_pair ////
+                    vec3 = _mm256_add_pd(vec3, vec5);
+
+
+                    // 6 store the transformed elements in vec3
+                    _mm256_storeu_pd(element+col_idx, vec3);
+
+
+                    //// u3_1qbit_10*element_vec ////
+
+                    // 1 calculate the multiplications  u3_1qbit_10*element_vec
+                    vec3 = _mm256_mul_pd(u3_1qbit_10_vec, element_vec); 
+
+                    // 4 Multiply elements of u3_1qbit_10*element_vec_permuted
+                    vec4 = _mm256_mul_pd(u3_1qbit_10_vec, element_vec_permuted); 
+
+                    // 5 Horizontally subtract the elements in vec3 and vec4
+                    vec3  = _mm256_hsub_pd(vec3, vec4); 
+
+
+                    //// u3_1qbit_01*element_vec_pair ////
+
+                    // 1 calculate the multiplications  u3_1qbit_01*element_pair_vec
+                    vec5 = _mm256_mul_pd(u3_1qbit_11_vec, element_pair_vec); 
+
+                    // 4 Multiply elements of u3_1qbit_01*element_vec_pair_permuted
+                    vec4 = _mm256_mul_pd(u3_1qbit_11_vec, element_pair_vec_permuted); 
+
+                    // 5 Horizontally subtract the elements in vec5 and vec4
+                    vec5  = _mm256_hsub_pd(vec5, vec4); 
+
+                    //// u3_1qbit_10*element_vec + u3_1qbit_11*element_vec_pair ////
+                    vec3 = _mm256_add_pd(vec3, vec5);
+
+                    // 6 store the transformed elements in vec3
+                    _mm256_storeu_pd(element_pair+col_idx, vec3);
+
+                }   
+
+                if (input.cols % 2 == 1) {     
+
+                    int col_idx = input.cols-1;       
+
+                    int index      = row_offset+col_idx;
+                    int index_pair = row_offset_pair+col_idx;                
+
+                    QGD_Complex16 element      = input[index];
+                    QGD_Complex16 element_pair = input[index_pair];              
+
+                    QGD_Complex16 tmp1 = mult(u3_1qbit[0], element);
+                    QGD_Complex16 tmp2 = mult(u3_1qbit[1], element_pair);
+ 
+                    input[index].real = tmp1.real + tmp2.real;
+                    input[index].imag = tmp1.imag + tmp2.imag;
+
+                    tmp1 = mult(u3_1qbit[2], element);
+                    tmp2 = mult(u3_1qbit[3], element_pair);
+
+                    input[index_pair].real = tmp1.real + tmp2.real;
+                    input[index_pair].imag = tmp1.imag + tmp2.imag;
+
+
+                }
+
+            }
+            else if (deriv) {
+                // when calculating derivatives, the constant element should be zeros
+                memset( input.get_data()+row_offset, 0.0, input.cols*sizeof(QGD_Complex16));
+                memset( input.get_data()+row_offset_pair, 0.0, input.cols*sizeof(QGD_Complex16));
+            }
+            else {
+                // leave the state as it is
+                continue; 
+            }
+
+
+//std::cout << current_idx_target << " " << current_idx_target_pair << std::endl;
+
+        
+        //});
+        }
+
+
+        current_idx = current_idx + (index_step_target << 1);
+        current_idx_pair = current_idx_pair + (index_step_target << 1);
+
+
+    }
+
+
+
+
+}
+
diff --git a/gates/kernels/include/apply_kernel_to_input_AVX.h b/gates/kernels/include/apply_kernel_to_input_AVX.h
new file mode 100644
index 00000000..7982ec26
--- /dev/null
+++ b/gates/kernels/include/apply_kernel_to_input_AVX.h
@@ -0,0 +1,38 @@
+/*
+Created on Fri Jun 26 14:13:26 2020
+Copyright (C) 2020 Peter Rakyta, Ph.D.
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this program.  If not, see http://www.gnu.org/licenses/.
+
+@author: Peter Rakyta, Ph.D.
+*/
+/*! \file apply_kerel_to_input_AVX.cpp
+    \brief ????????????????
+*/
+
+
+#ifndef apply_kerel_to_input_AVX_H
+#define apply_kerel_to_input_AVX_H
+
+#include "matrix.h"
+#include "common.h"
+/**
+@brief AVX kernel to apply single qubit gate kernel on an input matrix
+@param ????????
+@param ?????????
+*/
+void apply_kernel_to_input_AVX(Matrix& u3_1qbit, Matrix& input, const bool& deriv, const int& target_qbit, const int& control_qbit, const int& matrix_size);
+
+
+#endif
diff --git a/qgd_python/decomposition/qgd_N_Qubit_Decomposition_adaptive.py b/qgd_python/decomposition/qgd_N_Qubit_Decomposition_adaptive.py
index 6de4fdce..49d72221 100644
--- a/qgd_python/decomposition/qgd_N_Qubit_Decomposition_adaptive.py
+++ b/qgd_python/decomposition/qgd_N_Qubit_Decomposition_adaptive.py
@@ -41,7 +41,7 @@ class qgd_N_Qubit_Decomposition_adaptive(qgd_N_Qubit_Decomposition_adaptive_Wrap
 # @param optimize_layer_num Set true to optimize the minimum number of operation layers required in the decomposition, or false when the predefined maximal number of layer gates is used (ideal for general unitaries).
 # @param initial_guess String indicating the method to guess initial values for the optimalization. Possible values: "zeros" ,"random", "close_to_zero".
 # @return An instance of the class
-    def __init__( self, Umtx, level_limit_max=8, level_limit_min=0, topology=None ):
+    def __init__( self, Umtx, level_limit_max=8, level_limit_min=0, topology=None, accelerator_num=0 ):
 
         ## the number of qubits
         self.qbit_num = int(round( np.log2( len(Umtx) ) ))
@@ -65,7 +65,7 @@ def __init__( self, Umtx, level_limit_max=8, level_limit_min=0, topology=None ):
         
 
         # call the constructor of the wrapper class
-        super(qgd_N_Qubit_Decomposition_adaptive, self).__init__(Umtx, self.qbit_num, level_limit_max, level_limit_min, topology=topology_validated)
+        super(qgd_N_Qubit_Decomposition_adaptive, self).__init__(Umtx, self.qbit_num, level_limit_max, level_limit_min, topology=topology_validated, accelerator_num=accelerator_num)
 
 
 ##
@@ -115,6 +115,10 @@ def get_Quantum_Circuit( self ):
                 # adding CNOT gate to the quantum circuit
                 circuit.cx(gate.get("control_qbit"), gate.get("target_qbit"))
 
+            elif gate.get("type") == "CRY":
+                # adding CNOT gate to the quantum circuit
+                circuit.cry(gate.get("Theta"), gate.get("control_qbit"), gate.get("target_qbit"))
+
             elif gate.get("type") == "CZ":
                 # adding CZ gate to the quantum circuit
                 circuit.cz(gate.get("control_qbit"), gate.get("target_qbit"))
@@ -145,11 +149,19 @@ def get_Quantum_Circuit( self ):
                 circuit.rz(gate.get("Phi"), gate.get("target_qbit"))
 
             elif gate.get("type") == "X":
-                # RZ gate
+                # X gate
+                circuit.x(gate.get("target_qbit"))
+
+            elif gate.get("type") == "Y":
+                # Y gate
+                circuit.x(gate.get("target_qbit"))
+
+            elif gate.get("type") == "Z":
+                # Z gate
                 circuit.x(gate.get("target_qbit"))
 
             elif gate.get("type") == "SX":
-                # RZ gate
+                # SX gate
                 circuit.sx(gate.get("target_qbit"))
 
 
@@ -184,6 +196,10 @@ def get_Cirq_Circuit( self ):
                 # adding CNOT gate to the quantum circuit
                 circuit.append(cirq.CNOT(q[self.qbit_num-1-gate.get("control_qbit")], q[self.qbit_num-1-gate.get("target_qbit")]))
 
+            if gate.get("type") == "CRY":
+                # adding CRY gate to the quantum circuit
+                print("CRY gate needs to be implemented")
+
             elif gate.get("type") == "CZ":
                 # adding CZ gate to the quantum circuit
                 circuit.append(cirq.CZ(q[self.qbit_num-1-gate.get("control_qbit")], q[self.qbit_num-1-gate.get("target_qbit")]))
@@ -213,9 +229,18 @@ def get_Cirq_Circuit( self ):
                 circuit.append(cirq.rz(gate.get("Phi")).on(q[self.qbit_num-1-gate.get("target_qbit")]))
 
             elif gate.get("type") == "X":
-                # RZ gate
+                # X gate
                 circuit.append(cirq.x(q[self.qbit_num-1-gate.get("target_qbit")]))
 
+            elif gate.get("type") == "Y":
+                # Y gate
+                circuit.append(cirq.y(q[self.qbit_num-1-gate.get("target_qbit")]))
+
+            elif gate.get("type") == "Z":
+                # Z gate
+                circuit.append(cirq.z(q[self.qbit_num-1-gate.get("target_qbit")]))
+
+
             elif gate.get("type") == "SX":
                 # RZ gate
                 circuit.append(cirq.sx(q[self.qbit_num-1-gate.get("target_qbit")]))
@@ -346,7 +371,8 @@ def import_Qiskit_Circuit( self, qc_in ):
 
         # setting gate structure and optimized initial parameters
         self.set_Gate_Structure(Gates_Block_ret)
-        self.set_Optimized_Parameters( optimized_parameters )
+        self.set_Optimized_Parameters( np.flip(optimized_parameters,0) )
+        #self.set_Optimized_Parameters( optimized_parameters )
           
 ##
 # @brief Call to set custom layers to the gate structure that are intended to be used in the decomposition from a binary file created from SQUANDER
@@ -478,7 +504,7 @@ def get_Matrix( self, parameters = None ):
 ## 
 # @brief Call to set the optimizer used in the gate synthesis process
 # @param costfnc Variant of the cost function. Input argument 0 stands for FROBENIUS_NORM, 1 for FROBENIUS_NORM_CORRECTION1, 2 for FROBENIUS_NORM_CORRECTION2
-    def set_Cost_Function_Variant( self, costfnc="1" ):
+    def set_Cost_Function_Variant( self, costfnc=0 ):
 
         # Set the optimizer
         super(qgd_N_Qubit_Decomposition_adaptive, self).set_Cost_Function_Variant(costfnc=costfnc)  
@@ -497,17 +523,25 @@ def set_Iteration_Threshold_of_Randomization( self, threshold=2500 ):
 ## 
 # @brief Call to evaluate the cost function and the gradient components.
 # @param parameters A float64 numpy array
-# @param onlyCPU Perform calculation on CPU. Has effect if SQUANDER is compiled with DFE/Groq support.
-    def Optimization_Problem_Combined( self, parameters=None, onlyCPU=False ):
+    def Optimization_Problem_Combined( self, parameters=None ):
 
         if parameters is None:
             print( "Optimization_Problem_Combined: arary of input parameters is None")
             return None
 
         # evaluate the cost function and gradients
-        cost_function, grad = super(qgd_N_Qubit_Decomposition_adaptive, self).Optimization_Problem_Combined(parameters, onlyCPU)  
+        cost_function, grad = super(qgd_N_Qubit_Decomposition_adaptive, self).Optimization_Problem_Combined(parameters)  
 
         grad = grad.reshape( (-1,))
 
         return cost_function, grad
 
+
+
+## 
+# @brief Call to prepare the circuit to be exported into Qiskit format. (parameters and gates gets bound together, gate block structure is converted to plain structure).
+    def Prepare_Gates_To_Export(self):
+
+        # Set the optimizer
+        super(qgd_N_Qubit_Decomposition_adaptive, self).Prepare_Gates_To_Export()  
+
diff --git a/qgd_python/decomposition/qgd_N_Qubit_Decomposition_adaptive_Wrapper.cpp b/qgd_python/decomposition/qgd_N_Qubit_Decomposition_adaptive_Wrapper.cpp
index c6a433c6..5d901477 100644
--- a/qgd_python/decomposition/qgd_N_Qubit_Decomposition_adaptive_Wrapper.cpp
+++ b/qgd_python/decomposition/qgd_N_Qubit_Decomposition_adaptive_Wrapper.cpp
@@ -70,9 +70,9 @@ typedef struct qgd_N_Qubit_Decomposition_adaptive_Wrapper {
 @return Return with a void pointer pointing to an instance of N_Qubit_Decomposition class.
 */
 N_Qubit_Decomposition_adaptive* 
-create_N_Qubit_Decomposition_adaptive( Matrix& Umtx, int qbit_num, int level_limit, int level_limit_min, std::vector<matrix_base<int>> topology_in ) {
+create_N_Qubit_Decomposition_adaptive( Matrix& Umtx, int qbit_num, int level_limit, int level_limit_min, std::vector<matrix_base<int>> topology_in, int accelerator_num ) {
 
-    return new N_Qubit_Decomposition_adaptive( Umtx, qbit_num, level_limit, level_limit_min, topology_in );
+    return new N_Qubit_Decomposition_adaptive( Umtx, qbit_num, level_limit, level_limit_min, topology_in, accelerator_num );
 }
 
 
@@ -155,7 +155,7 @@ static int
 qgd_N_Qubit_Decomposition_adaptive_Wrapper_init(qgd_N_Qubit_Decomposition_adaptive_Wrapper *self, PyObject *args, PyObject *kwds)
 {
     // The tuple of expected keywords
-    static char *kwlist[] = {(char*)"Umtx", (char*)"qbit_num", (char*)"level_limit_min", (char*)"method", (char*)"topology", NULL};
+    static char *kwlist[] = {(char*)"Umtx", (char*)"qbit_num", (char*)"level_limit_min", (char*)"method", (char*)"topology", (char*)"accelerator_num", NULL};
  
     // initiate variables for input arguments
     PyObject *Umtx_arg = NULL;
@@ -163,10 +163,11 @@ qgd_N_Qubit_Decomposition_adaptive_Wrapper_init(qgd_N_Qubit_Decomposition_adapti
     int level_limit = 0;
     int level_limit_min = 0;
     PyObject *topology = NULL;
+    int accelerator_num = 0;
 
     // parsing input arguments
-    if (!PyArg_ParseTupleAndKeywords(args, kwds, "|OiiiO", kwlist,
-                                     &Umtx_arg, &qbit_num, &level_limit, &level_limit_min, &topology))
+    if (!PyArg_ParseTupleAndKeywords(args, kwds, "|OiiiOi", kwlist,
+                                     &Umtx_arg, &qbit_num, &level_limit, &level_limit_min, &topology, &accelerator_num))
         return -1;
 
     // convert python object array to numpy C API array
@@ -220,7 +221,13 @@ qgd_N_Qubit_Decomposition_adaptive_Wrapper_init(qgd_N_Qubit_Decomposition_adapti
 
     // create an instance of the class N_Qubit_Decomposition
     if (qbit_num > 0 ) {
-        self->decomp = create_N_Qubit_Decomposition_adaptive( Umtx_mtx, qbit_num, level_limit, level_limit_min, topology_Cpp);
+        try {
+            self->decomp = create_N_Qubit_Decomposition_adaptive( Umtx_mtx, qbit_num, level_limit, level_limit_min, topology_Cpp, accelerator_num);
+        }
+        catch (std::string err ) {
+            PyErr_SetString(PyExc_Exception, err.c_str());
+            return -1;
+        }
     }
     else {
         std::cout << "The number of qubits should be given as a positive integer, " << qbit_num << "  was given" << std::endl;
@@ -310,6 +317,7 @@ static PyObject *
 get_gate( N_Qubit_Decomposition_adaptive* decomp, int &idx ) {
 
 
+
     // create dictionary conatining the gate data
     PyObject* py_gate = PyDict_New();
 
@@ -477,6 +485,30 @@ get_gate( N_Qubit_Decomposition_adaptive* decomp, int &idx ) {
         Py_XDECREF(Phi);
 
 
+    }
+    else if (gate->get_type() == CRY_OPERATION ||  gate->get_type() == ADAPTIVE_OPERATION ) {
+
+        // get U3 parameters
+        CRY* cry_gate = static_cast<CRY*>(gate);
+        Matrix_real&& parameters = cry_gate->get_optimized_parameters();
+ 
+
+        // create gate parameters
+        PyObject* type = Py_BuildValue("s",  "CRY" );
+        PyObject* target_qbit = Py_BuildValue("i",  gate->get_target_qbit() );
+        PyObject* control_qbit = Py_BuildValue("i",  gate->get_control_qbit() );
+        PyObject* Theta = Py_BuildValue("f",  parameters[0] );
+
+        PyDict_SetItemString(py_gate, "type", type );
+        PyDict_SetItemString(py_gate, "target_qbit", target_qbit );
+        PyDict_SetItemString(py_gate, "control_qbit", control_qbit );            
+        PyDict_SetItemString(py_gate, "Theta", Theta );
+
+        Py_XDECREF(type);
+        Py_XDECREF(target_qbit);
+        Py_XDECREF(control_qbit);
+        Py_XDECREF(Theta);
+
     }
     else if (gate->get_type() == X_OPERATION) {
 
@@ -490,6 +522,32 @@ get_gate( N_Qubit_Decomposition_adaptive* decomp, int &idx ) {
         Py_XDECREF(type);
         Py_XDECREF(target_qbit);
 
+    }
+    else if (gate->get_type() == Y_OPERATION) {
+
+        // create gate parameters
+        PyObject* type = Py_BuildValue("s",  "Y" );
+        PyObject* target_qbit = Py_BuildValue("i",  gate->get_target_qbit() );
+
+        PyDict_SetItemString(py_gate, "type", type );
+        PyDict_SetItemString(py_gate, "target_qbit", target_qbit );
+
+        Py_XDECREF(type);
+        Py_XDECREF(target_qbit);
+
+    }
+    else if (gate->get_type() == Z_OPERATION) {
+
+        // create gate parameters
+        PyObject* type = Py_BuildValue("s",  "Z" );
+        PyObject* target_qbit = Py_BuildValue("i",  gate->get_target_qbit() );
+
+        PyDict_SetItemString(py_gate, "type", type );
+        PyDict_SetItemString(py_gate, "target_qbit", target_qbit );
+
+        Py_XDECREF(type);
+        Py_XDECREF(target_qbit);
+
     }
     else if (gate->get_type() == SX_OPERATION) {
 
@@ -641,23 +699,16 @@ static PyObject *
 qgd_N_Qubit_Decomposition_adaptive_Wrapper_get_Optimized_Parameters( qgd_N_Qubit_Decomposition_adaptive_Wrapper *self ) {
 
     int parameter_num = self->decomp->get_parameter_num();
-
-    matrix_base<double> parameters_mtx(1, parameter_num);
+    Matrix_real parameters_mtx(1, parameter_num);
     double* parameters = parameters_mtx.get_data();
     self->decomp->get_optimized_parameters(parameters);
 
-    // reversing the order
-    Matrix_real parameters_mtx_reversed(1, parameter_num);
-    double* parameters_reversed = parameters_mtx_reversed.get_data();
-    for (int idx=0; idx<parameter_num; idx++ ) {
-        parameters_reversed[idx] = parameters[parameter_num-1-idx];
-    }
-
     // convert to numpy array
-    parameters_mtx_reversed.set_owner(false);
-    PyObject * parameter_arr = matrix_real_to_numpy( parameters_mtx_reversed );
+    parameters_mtx.set_owner(false);
+    PyObject * parameter_arr = matrix_real_to_numpy( parameters_mtx );
 
     return parameter_arr;
+
 }
 
 
@@ -699,20 +750,10 @@ qgd_N_Qubit_Decomposition_adaptive_Wrapper_set_Optimized_Parameters( qgd_N_Qubit
     }
 
 
-    // get the pointer to the data stored in the input matrices
-    double* parameters = (double*)PyArray_DATA(parameters_arr);
-
-
-    npy_intp param_num = PyArray_Size( parameters_arr );
+    Matrix_real parameters_mtx = numpy2matrix_real( parameters_arr );
 
-    // reversing the order
-    matrix_base<double> parameters_mtx_reversed(param_num, 1);
-    double* parameters_reversed = parameters_mtx_reversed.get_data();
-    for (int idx=0; idx<param_num; idx++ ) {
-        parameters_reversed[idx] = parameters[param_num-1-idx];
-    }
 
-    self->decomp->set_optimized_parameters(parameters_reversed, param_num);
+    self->decomp->set_optimized_parameters(parameters_mtx.get_data(), parameters_mtx.size());
 
 
     Py_DECREF(parameters_arr);
@@ -1270,11 +1311,10 @@ qgd_N_Qubit_Decomposition_adaptive_Wrapper_Optimization_Problem_Combined( qgd_N_
 
 
     PyObject* parameters_arg = NULL;
-    bool onlyCPU = 0;
 
 
     // parsing input arguments
-    if (!PyArg_ParseTuple(args, "|Ob", &parameters_arg, &onlyCPU )) {
+    if (!PyArg_ParseTuple(args, "|Oi", &parameters_arg )) {
 
         std::string err( "Unsuccessful argument parsing not ");
         PyErr_SetString(PyExc_Exception, err.c_str());
@@ -1301,7 +1341,7 @@ qgd_N_Qubit_Decomposition_adaptive_Wrapper_Optimization_Problem_Combined( qgd_N_
     double f0;
 
     try {
-        self->decomp->optimization_problem_combined(parameters_mtx, &f0, grad_mtx, onlyCPU);
+        self->decomp->optimization_problem_combined(parameters_mtx, &f0, grad_mtx );
     }
     catch (std::string err ) {
         PyErr_SetString(PyExc_Exception, err.c_str());
@@ -1686,6 +1726,62 @@ qgd_N_Qubit_Decomposition_adaptive_Wrapper_set_Iteration_Threshold_of_Randomizat
 }
 
 
+/**
+@brief Call to upload the unitary to the DFE. (Has no effect for non-DFE builds)
+*/
+static PyObject *
+qgd_N_Qubit_Decomposition_adaptive_Wrapper_Upload_Umtx_to_DFE(qgd_N_Qubit_Decomposition_adaptive_Wrapper *self ) {
+
+#ifdef __DFE__
+
+    try {
+        self->decomp->upload_Umtx_to_DFE();
+    }
+    catch (std::string err) {
+        PyErr_SetString(PyExc_Exception, err.c_str());
+        std::cout << err << std::endl;
+        return NULL;
+    }
+    catch(...) {
+        std::string err( "Invalid pointer to decomposition class");
+        PyErr_SetString(PyExc_Exception, err.c_str());
+        return NULL;
+    }
+
+#endif
+
+    return Py_BuildValue("i", 0);
+    
+}
+
+
+
+/**
+@brief Call to prepare the circuit to be exported into Qiskit format. (parameters and gates gets bound together, gate block structure is converted to plain structure).
+*/
+static PyObject *
+qgd_N_Qubit_Decomposition_adaptive_Wrapper_Prepare_Gates_To_Export(qgd_N_Qubit_Decomposition_adaptive_Wrapper *self ) {
+
+    try {
+        self->decomp->prepare_gates_to_export();
+    }
+    catch (std::string err) {
+        PyErr_SetString(PyExc_Exception, err.c_str());
+        std::cout << err << std::endl;
+        return NULL;
+    }
+    catch(...) {
+        std::string err( "Invalid pointer to decomposition class");
+        PyErr_SetString(PyExc_Exception, err.c_str());
+        return NULL;
+    }
+
+
+    return Py_BuildValue("i", 0);
+    
+}
+
+
 /**
 @brief Structure containing metadata about the members of class qgd_N_Qubit_Decomposition_adaptive_Wrapper.
 */
@@ -1811,6 +1907,12 @@ static PyMethodDef qgd_N_Qubit_Decomposition_adaptive_Wrapper_methods[] = {
     {"Optimization_Problem_Combined", (PyCFunction) qgd_N_Qubit_Decomposition_adaptive_Wrapper_Optimization_Problem_Combined, METH_VARARGS,
      "Wrapper function to evaluate the cost function and the gradient components."
     },
+    {"Upload_Umtx_to_DFE", (PyCFunction) qgd_N_Qubit_Decomposition_adaptive_Wrapper_Upload_Umtx_to_DFE, METH_NOARGS,
+     "Call to upload the unitary to the DFE. (Has no effect for non-DFE builds)"
+    },
+    {"Prepare_Gates_To_Export", (PyCFunction) qgd_N_Qubit_Decomposition_adaptive_Wrapper_Prepare_Gates_To_Export, METH_NOARGS,
+     "Call to prepare the circuit to be exported into Qiskit format. (parameters and gates gets bound together, gate block structure is converted to plain structure)."
+    },
     {NULL}  /* Sentinel */
 };
 
diff --git a/qgd_python/decomposition/qgd_N_Qubit_Decomposition_custom.py b/qgd_python/decomposition/qgd_N_Qubit_Decomposition_custom.py
index d0ab8c66..0250ce5f 100644
--- a/qgd_python/decomposition/qgd_N_Qubit_Decomposition_custom.py
+++ b/qgd_python/decomposition/qgd_N_Qubit_Decomposition_custom.py
@@ -41,13 +41,13 @@ class qgd_N_Qubit_Decomposition_custom(qgd_N_Qubit_Decomposition_custom_Wrapper)
 # @param optimize_layer_num Set true to optimize the minimum number of operation layers required in the decomposition, or false when the predefined maximal number of layer gates is used (ideal for general unitaries).
 # @param initial_guess String indicating the method to guess initial values for the optimalization. Possible values: "zeros" ,"random", "close_to_zero".
 # @return An instance of the class
-    def __init__( self, Umtx, initial_guess="RANDOM" ):
+    def __init__( self, Umtx, initial_guess="RANDOM", accelerator_num=0 ):
 
         ## the number of qubits
         self.qbit_num = int(round( np.log2( len(Umtx) ) ))
 
         # call the constructor of the wrapper class
-        super(qgd_N_Qubit_Decomposition_custom, self).__init__(Umtx, self.qbit_num, initial_guess)
+        super(qgd_N_Qubit_Decomposition_custom, self).__init__(Umtx, self.qbit_num, initial_guess, accelerator_num)
 
 
 
@@ -98,6 +98,10 @@ def get_Quantum_Circuit( self ):
                 # adding CNOT gate to the quantum circuit
                 circuit.cx(gate.get("control_qbit"), gate.get("target_qbit"))
 
+            elif gate.get("type") == "CRY":
+                # adding CNOT gate to the quantum circuit
+                circuit.cry(gate.get("Theta"), gate.get("control_qbit"), gate.get("target_qbit"))
+
             elif gate.get("type") == "CZ":
                 # adding CZ gate to the quantum circuit
                 circuit.cz(gate.get("control_qbit"), gate.get("target_qbit"))
@@ -108,7 +112,7 @@ def get_Quantum_Circuit( self ):
 
             elif gate.get("type") == "SYC":
                 # Sycamore gate
-                print("Unsupported gate in the Qiskit circuit export: Sycamore gate")
+                print("Unsupported gate in the circuit export: Sycamore gate")
                 return None;
 
             elif gate.get("type") == "U3":
@@ -128,11 +132,19 @@ def get_Quantum_Circuit( self ):
                 circuit.rz(gate.get("Phi"), gate.get("target_qbit"))
 
             elif gate.get("type") == "X":
-                # RZ gate
+                # X gate
+                circuit.x(gate.get("target_qbit"))
+
+            elif gate.get("type") == "Y":
+                # Y gate
+                circuit.x(gate.get("target_qbit"))
+
+            elif gate.get("type") == "Z":
+                # Z gate
                 circuit.x(gate.get("target_qbit"))
 
             elif gate.get("type") == "SX":
-                # RZ gate
+                # SX gate
                 circuit.sx(gate.get("target_qbit"))
 
 
@@ -167,6 +179,10 @@ def get_Cirq_Circuit( self ):
                 # adding CNOT gate to the quantum circuit
                 circuit.append(cirq.CNOT(q[self.qbit_num-1-gate.get("control_qbit")], q[self.qbit_num-1-gate.get("target_qbit")]))
 
+            if gate.get("type") == "CRY":
+                # adding CRY gate to the quantum circuit
+                print("CRY gate needs to be implemented")
+
             elif gate.get("type") == "CZ":
                 # adding CZ gate to the quantum circuit
                 circuit.append(cirq.CZ(q[self.qbit_num-1-gate.get("control_qbit")], q[self.qbit_num-1-gate.get("target_qbit")]))
@@ -199,15 +215,22 @@ def get_Cirq_Circuit( self ):
                 # X gate
                 circuit.append(cirq.x(q[self.qbit_num-1-gate.get("target_qbit")]))
 
+            elif gate.get("type") == "Y":
+                # Y gate
+                circuit.append(cirq.y(q[self.qbit_num-1-gate.get("target_qbit")]))
+
+            elif gate.get("type") == "Z":
+                # Z gate
+                circuit.append(cirq.z(q[self.qbit_num-1-gate.get("target_qbit")]))
+
+
             elif gate.get("type") == "SX":
-                # SX gate
-                print('SX!!!!!!!!!!!!!!!!!!!')
+                # RZ gate
                 circuit.append(cirq.sx(q[self.qbit_num-1-gate.get("target_qbit")]))
 
 
         return circuit
 
-
         
 ##
 # @brief Call to import initial quantum circuit in QISKIT format to be further comporessed
@@ -326,7 +349,7 @@ def import_Qiskit_Circuit( self, qc_in ):
 
         # setting gate structure and optimized initial parameters
         self.set_Gate_Structure(Gates_Block_ret)
-        self.set_Optimized_Parameters( optimized_parameters )      
+        self.set_Optimized_Parameters( np.flip(optimized_parameters,0) )      
         
 
 
@@ -341,3 +364,10 @@ def set_Optimizer( self, optimizer="BFGS" ):
 
 
 
+## 
+# @brief Call to prepare the circuit to be exported into Qiskit format. (parameters and gates gets bound together, gate block structure is converted to plain structure).
+    def Prepare_Gates_To_Export(self):
+
+        # Set the optimizer
+        super(qgd_N_Qubit_Decomposition_custom, self).Prepare_Gates_To_Export()  
+
diff --git a/qgd_python/decomposition/qgd_N_Qubit_Decomposition_custom_Wrapper.cpp b/qgd_python/decomposition/qgd_N_Qubit_Decomposition_custom_Wrapper.cpp
index e5a74ed6..182f8b63 100644
--- a/qgd_python/decomposition/qgd_N_Qubit_Decomposition_custom_Wrapper.cpp
+++ b/qgd_python/decomposition/qgd_N_Qubit_Decomposition_custom_Wrapper.cpp
@@ -68,9 +68,9 @@ typedef struct qgd_N_Qubit_Decomposition_custom_Wrapper {
 @return Return with a void pointer pointing to an instance of N_Qubit_Decomposition class.
 */
 N_Qubit_Decomposition_custom* 
-create_N_Qubit_Decomposition_custom( Matrix& Umtx, int qbit_num, bool optimize_layer_num, guess_type initial_guess ) {
+create_N_Qubit_Decomposition_custom( Matrix& Umtx, int qbit_num, bool optimize_layer_num, guess_type initial_guess, int accelerator_num ) {
 
-    return new N_Qubit_Decomposition_custom( Umtx, qbit_num, optimize_layer_num, initial_guess );
+    return new N_Qubit_Decomposition_custom( Umtx, qbit_num, optimize_layer_num, initial_guess, accelerator_num );
 }
 
 
@@ -144,16 +144,17 @@ static int
 qgd_N_Qubit_Decomposition_custom_Wrapper_init(qgd_N_Qubit_Decomposition_custom_Wrapper *self, PyObject *args, PyObject *kwds)
 {
     // The tuple of expected keywords
-    static char *kwlist[] = {(char*)"Umtx", (char*)"qbit_num", (char*)"initial_guess", NULL};
+    static char *kwlist[] = {(char*)"Umtx", (char*)"qbit_num", (char*)"initial_guess", (char*)"accelerator_num", NULL};
  
     // initiate variables for input arguments
     PyObject *Umtx_arg = NULL;
     int  qbit_num = -1; 
     PyObject *initial_guess = NULL;
+    int accelerator_num = 0;
 
     // parsing input arguments
-    if (!PyArg_ParseTupleAndKeywords(args, kwds, "|OiO", kwlist,
-                                     &Umtx_arg, &qbit_num, &initial_guess))
+    if (!PyArg_ParseTupleAndKeywords(args, kwds, "|OiOi", kwlist,
+                                     &Umtx_arg, &qbit_num, &initial_guess, &accelerator_num))
         return -1;
 
     // convert python object array to numpy C API array
@@ -192,7 +193,13 @@ qgd_N_Qubit_Decomposition_custom_Wrapper_init(qgd_N_Qubit_Decomposition_custom_W
   
     // create an instance of the class N_Qubit_Decomposition_custom
     if (qbit_num > 0 ) {
-        self->decomp =  create_N_Qubit_Decomposition_custom( Umtx_mtx, qbit_num, false, qgd_initial_guess);
+        try {
+            self->decomp =  create_N_Qubit_Decomposition_custom( Umtx_mtx, qbit_num, false, qgd_initial_guess, accelerator_num);
+        }
+        catch (std::string err ) {
+            PyErr_SetString(PyExc_Exception, err.c_str());
+            return -1;
+        }
     }
     else {
         std::cout << "The number of qubits should be given as a positive integer, " << qbit_num << "  was given" << std::endl;
@@ -564,16 +571,9 @@ qgd_N_Qubit_Decomposition_custom_Wrapper_get_Optimized_Parameters( qgd_N_Qubit_D
     double* parameters = parameters_mtx.get_data();
     self->decomp->get_optimized_parameters(parameters);
 
-    // reversing the order
-    Matrix_real parameters_mtx_reversed(1, parameter_num);
-    double* parameters_reversed = parameters_mtx_reversed.get_data();
-    for (int idx=0; idx<parameter_num; idx++ ) {
-        parameters_reversed[idx] = parameters[parameter_num-1-idx];
-    }
-
     // convert to numpy array
-    parameters_mtx_reversed.set_owner(false);
-    PyObject * parameter_arr = matrix_real_to_numpy( parameters_mtx_reversed );
+    parameters_mtx.set_owner(false);
+    PyObject * parameter_arr = matrix_real_to_numpy( parameters_mtx );
 
     return parameter_arr;
 }
@@ -605,21 +605,10 @@ qgd_N_Qubit_Decomposition_custom_Wrapper_set_Optimized_Parameters( qgd_N_Qubit_D
     }
 
 
-    // get the pointer to the data stored in the input matrices
-    double* parameters = (double*)PyArray_DATA(parameters_arr);
-
-
-    npy_intp param_num = PyArray_Size( parameters_arr );
-
-    // reversing the order
-    Matrix_real parameters_mtx_reversed(param_num, 1);
-    double* parameters_reversed = parameters_mtx_reversed.get_data();
-    for (int idx=0; idx<param_num; idx++ ) {
-        parameters_reversed[idx] = parameters[param_num-1-idx];
-    }
+    Matrix_real parameters_mtx = numpy2matrix_real( parameters_arr );
 
 
-    self->decomp->set_optimized_parameters(parameters_reversed, param_num);
+    self->decomp->set_optimized_parameters(parameters_mtx.get_data(), parameters_mtx.size());
 
 
     Py_DECREF(parameters_arr);
@@ -1012,6 +1001,34 @@ qgd_N_Qubit_Decomposition_custom_Wrapper_set_Optimizer( qgd_N_Qubit_Decompositio
 }
 
 
+/**
+@brief Call to upload the unitary to the DFE. (Has no effect for non-DFE builds)
+*/
+static PyObject *
+qgd_N_Qubit_Decomposition_custom_Wrapper_Upload_Umtx_to_DFE(qgd_N_Qubit_Decomposition_custom_Wrapper *self ) {
+
+#ifdef __DFE__
+
+    try {
+        self->decomp->upload_Umtx_to_DFE();
+    }
+    catch (std::string err) {
+        PyErr_SetString(PyExc_Exception, err.c_str());
+        std::cout << err << std::endl;
+        return NULL;
+    }
+    catch(...) {
+        std::string err( "Invalid pointer to decomposition class");
+        PyErr_SetString(PyExc_Exception, err.c_str());
+        return NULL;
+    }
+
+#endif
+
+    return Py_BuildValue("i", 0);
+    
+}
+
 
 
 /**
@@ -1076,6 +1093,9 @@ static PyMethodDef qgd_N_Qubit_Decomposition_custom_Wrapper_methods[] = {
     {"set_Optimizer", (PyCFunction) qgd_N_Qubit_Decomposition_custom_Wrapper_set_Optimizer, METH_VARARGS | METH_KEYWORDS,
      "Wrapper method to to set the optimizer method for the gate synthesis."
     },
+    {"Upload_Umtx_to_DFE", (PyCFunction) qgd_N_Qubit_Decomposition_custom_Wrapper_Upload_Umtx_to_DFE, METH_NOARGS,
+     "Call to upload the unitary to the DFE. (Has no effect for non-DFE builds)"
+    },
     {NULL}  /* Sentinel */
 };
 


Test patch:
diff --git a/qgd_python/decomposition/test/test_IBM.py b/qgd_python/decomposition/test/test_IBM.py
index fc7f9eac..56b98e5d 100644
--- a/qgd_python/decomposition/test/test_IBM.py
+++ b/qgd_python/decomposition/test/test_IBM.py
@@ -109,7 +109,7 @@ def test_IBM_Chellenge(self):
 
 
 
-    def ctest_IBM_Chellenge_adaptive(self):
+    def test_IBM_Chellenge_adaptive(self):
         r"""
         This method is called by pytest. 
         Test to decompose a 4-qubit unitary of the IBM chellenge
diff --git a/qgd_python/decomposition/test/test_parametric_circuit.py b/qgd_python/decomposition/test/test_parametric_circuit.py
index 101b18bd..7e13d890 100644
--- a/qgd_python/decomposition/test/test_parametric_circuit.py
+++ b/qgd_python/decomposition/test/test_parametric_circuit.py
@@ -114,7 +114,7 @@ def get_unitary_distance( Umtx1, Umtx2 ):
 ##
 # @brief ???????????
 # @return ???????????
-def get_optimized_circuit( alpha, optimizer='BFGS' ):
+def get_optimized_circuit( alpha, optimizer='BFGS', optimized_parameters=None ):
 	
 	filename = 'qgd_python/decomposition/test/19CNOT.qasm'
 	qc_trial = QuantumCircuit.from_qasm_file( filename )
@@ -146,6 +146,10 @@ def get_optimized_circuit( alpha, optimizer='BFGS' ):
                 # set the optimizer
 		cDecompose.set_Optimizer( optimizer )
 
+		# set the initial parameters if given
+		if not ( optimized_parameters is None ):
+			cDecompose.set_Optimized_Parameters( optimized_parameters )
+
 		# set verbosity
 		cDecompose.set_Verbose( 4 )
 
@@ -186,8 +190,14 @@ def test_optimizer(self):
         # determine random parameter value alpha
         alpha = 1.823631161607293
 
+        # determine the quantum circuit at parameter value alpha with BFGS2 optimizer
+        qc, optimized_parameters = get_optimized_circuit( alpha, optimizer='BFGS2' )
+
         # determine the quantum circuit at parameter value alpha with BFGS optimizer
-        qc = get_optimized_circuit( alpha, optimizer='BFGS2' )
+        qc, optimized_parameters = get_optimized_circuit( alpha, optimizer='BFGS', optimized_parameters=optimized_parameters )
+
+        # determine the quantum circuit at parameter value alpha with ADAM optimizer
+        qc, optimized_parameters = get_optimized_circuit( alpha+0.005, optimizer='ADAM', optimized_parameters=optimized_parameters )