update bounds and pf format

quantum_simulation_recipe/bounds.py

@@ -3,6 +3,7 @@
 from scipy.linalg import expm
 from numpy.linalg import matrix_power
 import scipy.sparse.linalg as ssla
+from scipy.sparse import csr_matrix, csc_matrix
 
 import numpy as np
 from quantum_simulation_recipe.measure import commutator, norm
@@ -65,6 +66,8 @@ def tight_bound(h_list: list, order: int, t: float, r: int, type='spectral', ver
     elif isinstance(h_list[0], SparsePauliOp):
         n = h_list[0].num_qubits
         d = 2**n
+    # elif isinstance(h_list[0], csr_matrix):
+    #     d = h_list[0].todense().shape[0]
     else:
         raise ValueError('Hamiltonian type is not defined')
 

quantum_simulation_recipe/trotter.py

@@ -42,16 +42,17 @@ def pf(h_list, t, r: int, order: int=2, use_jax=False, return_exact=False, verbo
 # def pf_r(h_list, t, r: int, order: int=2, use_jax=False, return_exact=False, verbose=False):
     if order == 1:
         list_U = [expH(herm, t/r, use_jax=use_jax) for herm in h_list]
-        appro_U_dt = np.linalg.multi_dot(list_U)
+        # appro_U_dt = np.linalg.multi_dot(list_U)
+        appro_U_dt = sparse_multi_dot(list_U)
         if isinstance(appro_U_dt, csr_matrix):
             appro_U = appro_U_dt**r
         else:
             appro_U = np.linalg.matrix_power(appro_U_dt, r)
     elif order == 2:
         list_U = [expH(herm, t/(2*r), use_jax=use_jax) for herm in h_list]
         if verbose: print('----expm Herm finished----')
-        appro_U_dt_forward = np.linalg.multi_dot(list_U)
-        appro_U_dt_reverse = np.linalg.multi_dot(list_U[::-1])
+        appro_U_dt_forward = sparse_multi_dot(list_U)
+        appro_U_dt_reverse = sparse_multi_dot(list_U[::-1])
         # appro_U_dt = list_U[0] @ list_U[1]
         if verbose: print('----matrix product finished----')
         if isinstance(appro_U_dt_forward, csr_matrix):