QuBase independent and sparse support.

src/QuDynamics.jl

@@ -1,5 +1,5 @@
 module QuDynamics
-    using QuBase
+    # using QuBase
     using Compat
     using ODE
     using Expokit
@@ -11,10 +11,10 @@ module QuDynamics
     include("propodesolvers.jl")
     include("propexpmsolvers.jl")
     include("propmcwf.jl")
-    module QuTiP
-    using ..QuBase
-    using ..QuDynamics
-    include("qutipinterface.jl")
-    end
+    # module QuTiP
+    # using ..QuBase
+    # using ..QuDynamics
+    # include("qutipinterface.jl")
+    # end
     export propagate
 end

src/propexpmsolvers.jl

@@ -39,18 +39,16 @@ QuExpmV() = QuExpmV(Dict())
 function propagate(prob::QuExpokit, eq::QuEquation, t, current_t, current_qustate)
     dt = t - current_t
     dims = size(current_qustate)
-    next_state = Expokit.expmv(dt, -im*coeffs(operator(eq, t)), coeffs(vec(current_qustate)), m = get(prob.options, :m, 30), tol = get(prob.options, :tol, 1e-7))
-    CQST = QuBase.similar_type(current_qustate)
-    return CQST(reshape(next_state, dims), bases(current_qustate))
+    next_state = Expokit.expmv(dt, -im*operator(eq, t), vec(current_qustate), m = get(prob.options, :m, 30), tol = get(prob.options, :tol, 1e-7))
+    return reshape(next_state, dims)
 end
 
 function propagate(prob::QuExpmV, eq::QuEquation, t, current_t, current_qustate)
     dt = t - current_t
     dims = size(current_qustate)
-    next_state = ExpmV.expmv(dt, -im*coeffs(operator(eq, t)), coeffs(vec(current_qustate)), M = get(prob.options, :M, []), prec = get(prob.options, :prec, "double"),
+    next_state = ExpmV.expmv(dt, -im*operator(eq, t), vec(current_qustate), M = get(prob.options, :M, []), prec = get(prob.options, :prec, "double"),
                             shift = get(prob.options, :shift, false), full_term = get(prob.options, :full_term, false))
-    CQST = QuBase.similar_type(current_qustate)
-    return CQST(reshape(next_state, dims), bases(current_qustate))
+    return reshape(next_state, dims)
 end
 
 export QuExpokit,

src/propmachinery.jl

@@ -27,32 +27,65 @@ Inputs :
 Output :
 * QuPropagator construct depending on the input.
 """ ->
-immutable QuStateEvolution{QPM<:QuPropagatorMethod, QVM<:@compat(Union{QuBase.AbstractQuVector,QuBase.AbstractQuMatrix}), QE<:QuEquation}
+immutable QuStateEvolution{QPM<:QuPropagatorMethod, QVM<:@compat(Union{AbstractVector, AbstractMatrix, SparseVector, SparseMatrixCSC}), QE<:QuEquation}
     eq::QE
     init_state::QVM
     tlist
     method::QPM
     QuStateEvolution(eq, init_state, tlist, method) = new(eq, init_state, tlist, method)
 end
 
-QuStateEvolution{QPM<:QuPropagatorMethod, QV<:QuBase.AbstractQuVector, QE<: QuEquation}(eq::QE, init_state::QV, tlist, method::QPM) = QuStateEvolution{QPM,QV,QE}(eq, init_state, tlist, method)
+# QuStateEvolution{QPM<:QuPropagatorMethod, QV<:QuBase.AbstractQuVector, QE<: QuEquation}(eq::QE, init_state::QV, tlist, method::QPM) = QuStateEvolution{QPM,QV,QE}(eq, init_state, tlist, method)
 
-QuStateEvolution{QPM<:QuPropagatorMethod, QM<:QuBase.AbstractQuMatrix, QE<: QuEquation}(eq::QE, init_state::QM, tlist, method::QPM) = QuStateEvolution{QPM,QM,QE}(eq, init_state, tlist, method)
+QuStateEvolution{QPM<:QuPropagatorMethod, QV<:Union{AbstractVector, SparseVector}, QE<: QuEquation}(eq::QE, init_state::QV, tlist, method::QPM) = QuStateEvolution{QPM,QV,QE}(eq, init_state, tlist, method)
 
-QuStateEvolution{QPM<:QuPropagatorMethod, QV<:QuBase.AbstractQuVector}(eq::QuSchrodingerEq, init_state::QV, tlist, method::QPM) = QuStateEvolution{QPM,QV,QuSchrodingerEq}(eq, init_state, tlist, method)
+# QuStateEvolution{QPM<:QuPropagatorMethod, QM<:QuBase.AbstractQuMatrix, QE<: QuEquation}(eq::QE, init_state::QM, tlist, method::QPM) = QuStateEvolution{QPM,QM,QE}(eq, init_state, tlist, method)
 
-QuStateEvolution{QPM<:QuPropagatorMethod, QV<:QuBase.AbstractQuVector}(hamiltonian::QuBase.AbstractQuMatrix, init_state::QV,  tlist, method::QPM) = QuStateEvolution{QPM,QV,QuSchrodingerEq}(QuSchrodingerEq(hamiltonian),init_state, tlist, method)
+QuStateEvolution{QPM<:QuPropagatorMethod, QM<:Union{AbstractMatrix, SparseMatrixCSC}, QE<: QuEquation}(eq::QE, init_state::QM, tlist, method::QPM) = QuStateEvolution{QPM,QM,QE}(eq, init_state, tlist, method)
 
-QuStateEvolution{QPM<:QuPropagatorMethod, QM<:QuBase.AbstractQuMatrix}(eq::QuLiouvillevonNeumannEq, init_state::QM, tlist, method::QPM) = QuStateEvolution{QPM,QM,QuLiouvillevonNeumannEq}(eq, init_state, tlist, method)
+# QuStateEvolution{QPM<:QuPropagatorMethod, QV<:QuBase.AbstractQuVector}(eq::QuSchrodingerEq, init_state::QV, tlist, method::QPM) = QuStateEvolution{QPM,QV,QuSchrodingerEq}(eq, init_state, tlist, method)
 
-QuStateEvolution{QPM<:QuPropagatorMethod, QM<:QuBase.AbstractQuMatrix}(hamiltonian::QuBase.AbstractQuMatrix, init_state::QM,  tlist, method::QPM) = QuStateEvolution{QPM,QM,QuLiouvillevonNeumannEq}(QuLiouvillevonNeumannEq(liouvillian_op(hamiltonian)),init_state, tlist, method)
+QuStateEvolution{QPM<:QuPropagatorMethod, QV<:AbstractVector}(eq::QuSchrodingerEq, init_state::QV, tlist, method::QPM) = QuStateEvolution{QPM,QV,QuSchrodingerEq}(eq, init_state, tlist, method)
 
-QuStateEvolution{QPM<:QuPropagatorMethod, QM<:QuBase.AbstractQuMatrix}(eq::QuLindbladMasterEq, init_state::QM, tlist, method::QPM) = QuStateEvolution{QPM,QM,QuLindbladMasterEq}(eq, init_state, tlist, method)
+QuStateEvolution{QPM<:QuPropagatorMethod, QV<:SparseVector}(eq::QuSchrodingerEqSparse, init_state::QV, tlist, method::QPM) = QuStateEvolution{QPM,QV,QuSchrodingerEqSparse}(eq, SparseMatrixCSC(init_state), tlist, method)
 
-QuStateEvolution{QPM<:QuPropagatorMethod, QM<:QuBase.AbstractQuMatrix, COT<:QuBase.AbstractQuMatrix}(hamiltonian::QuBase.AbstractQuMatrix, collapse_ops::Vector{COT}, init_state::QM, tlist, method::QPM) = QuStateEvolution{QPM,QM,QuLindbladMasterEq}(QuLindbladMasterEq(hamiltonian,collapse_ops), init_state, tlist, method)
+# QuStateEvolution{QPM<:QuPropagatorMethod, QV<:QuBase.AbstractQuVector}(hamiltonian::QuBase.AbstractQuMatrix, init_state::QV,  tlist, method::QPM) = QuStateEvolution{QPM,QV,QuSchrodingerEq}(QuSchrodingerEq(hamiltonian),init_state, tlist, method)
+
+QuStateEvolution{QPM<:QuPropagatorMethod, QV<:AbstractVector}(hamiltonian::AbstractMatrix, init_state::QV,  tlist, method::QPM) = QuStateEvolution{QPM,QV,QuSchrodingerEq}(QuSchrodingerEq(hamiltonian),init_state, tlist, method)
+
+QuStateEvolution{QPM<:QuPropagatorMethod, QV<:SparseVector}(hamiltonian::SparseMatrixCSC, init_state::QV,  tlist, method::QPM) = QuStateEvolution{QPM,QV,QuSchrodingerEqSparse}(QuSchrodingerEqSparse(hamiltonian), SparseMatrixCSC(init_state), tlist, method)
+
+# QuStateEvolution{QPM<:QuPropagatorMethod, QM<:QuBase.AbstractQuMatrix}(eq::QuLiouvillevonNeumannEq, init_state::QM, tlist, method::QPM) = QuStateEvolution{QPM,QM,QuLiouvillevonNeumannEq}(eq, init_state, tlist, method)
+
+QuStateEvolution{QPM<:QuPropagatorMethod, QM<:AbstractMatrix}(eq::QuLiouvillevonNeumannEq, init_state::QM, tlist, method::QPM) = QuStateEvolution{QPM,QM,QuLiouvillevonNeumannEq}(eq, init_state, tlist, method)
+
+QuStateEvolution{QPM<:QuPropagatorMethod, QM<:SparseMatrixCSC}(eq::QuLiouvillevonNeumannEqSparse, init_state::QM, tlist, method::QPM) = QuStateEvolution{QPM,QM,QuLiouvillevonNeumannEqSparse}(eq, init_state, tlist, method)
+
+# QuStateEvolution{QPM<:QuPropagatorMethod, QM<:QuBase.AbstractQuMatrix}(hamiltonian::QuBase.AbstractQuMatrix, init_state::QM,  tlist, method::QPM) = QuStateEvolution{QPM,QM,QuLiouvillevonNeumannEq}(QuLiouvillevonNeumannEq(liouvillian_op(hamiltonian)),init_state, tlist, method)
+
+QuStateEvolution{QPM<:QuPropagatorMethod, QM<:AbstractMatrix}(hamiltonian::AbstractMatrix, init_state::QM,  tlist, method::QPM) = QuStateEvolution{QPM,QM,QuLiouvillevonNeumannEq}(QuLiouvillevonNeumannEq(liouvillian_op(hamiltonian)),init_state, tlist, method)
+
+QuStateEvolution{QPM<:QuPropagatorMethod, QM<:SparseMatrixCSC}(hamiltonian::SparseMatrixCSC, init_state::QM,  tlist, method::QPM) = QuStateEvolution{QPM,QM,QuLiouvillevonNeumannEqSparse}(QuLiouvillevonNeumannEqSparse(liouvillian_op(hamiltonian)),init_state, tlist, method)
+
+# QuStateEvolution{QPM<:QuPropagatorMethod, QM<:QuBase.AbstractQuMatrix}(eq::QuLindbladMasterEq, init_state::QM, tlist, method::QPM) = QuStateEvolution{QPM,QM,QuLindbladMasterEq}(eq, init_state, tlist, method)
+
+QuStateEvolution{QPM<:QuPropagatorMethod, QM<:AbstractMatrix}(eq::QuLindbladMasterEq, init_state::QM, tlist, method::QPM) = QuStateEvolution{QPM,QM,QuLindbladMasterEq}(eq, init_state, tlist, method)
+
+QuStateEvolution{QPM<:QuPropagatorMethod, QM<:SparseMatrixCSC}(eq::QuLindbladMasterEqSparse, init_state::QM, tlist, method::QPM) = QuStateEvolution{QPM,QM,QuLindbladMasterEqSparse}(eq, init_state, tlist, method)
+
+
+# QuStateEvolution{QPM<:QuPropagatorMethod, QM<:QuBase.AbstractQuMatrix, COT<:QuBase.AbstractQuMatrix}(hamiltonian::QuBase.AbstractQuMatrix, collapse_ops::Vector{COT}, init_state::QM, tlist, method::QPM) = QuStateEvolution{QPM,QM,QuLindbladMasterEq}(QuLindbladMasterEq(hamiltonian,collapse_ops), init_state, tlist, method)
+
+QuStateEvolution{QPM<:QuPropagatorMethod, QM<:AbstractMatrix, COT<:AbstractMatrix}(hamiltonian::AbstractMatrix, collapse_ops::Vector{COT}, init_state::QM, tlist, method::QPM) = QuStateEvolution{QPM,QM,QuLindbladMasterEq}(QuLindbladMasterEq(hamiltonian,collapse_ops), init_state, tlist, method)
+
+QuStateEvolution{QPM<:QuPropagatorMethod, QM<:SparseMatrixCSC, COT<:SparseMatrixCSC}(hamiltonian::SparseMatrixCSC, collapse_ops::Vector{COT}, init_state::QM, tlist, method::QPM) = QuStateEvolution{QPM,QM,QuLindbladMasterEqSparse}(QuLindbladMasterEqSparse(hamiltonian,collapse_ops), init_state, tlist, method)
 
 # QuLindbladMasterEqUncached is used as the construction of Lindblad operator is not required for every tracjectory.
-QuStateEvolution{QPM<:QuPropagatorMethod, QV<:QuBase.AbstractQuVector, COT<:QuBase.AbstractQuMatrix}(hamiltonian::QuBase.AbstractQuMatrix, collapse_ops::Vector{COT}, init_state::QV, tlist, method::QPM) = QuStateEvolution{QPM,QV,QuLindbladMasterEq}(QuLindbladMasterEqUncached(hamiltonian,collapse_ops), init_state, tlist, method)
+# QuStateEvolution{QPM<:QuPropagatorMethod, QV<:QuBase.AbstractQuVector, COT<:QuBase.AbstractQuMatrix}(hamiltonian::QuBase.AbstractQuMatrix, collapse_ops::Vector{COT}, init_state::QV, tlist, method::QPM) = QuStateEvolution{QPM,QV,QuLindbladMasterEq}(QuLindbladMasterEqUncached(hamiltonian,collapse_ops), init_state, tlist, method)
+
+QuStateEvolution{QPM<:QuPropagatorMethod, QV<:AbstractVector, COT<:AbstractMatrix}(hamiltonian::AbstractMatrix, collapse_ops::Vector{COT}, init_state::QV, tlist, method::QPM) = QuStateEvolution{QPM,QV,QuLindbladMasterEq}(QuLindbladMasterEqUncached(hamiltonian,collapse_ops), init_state, tlist, method)
+
+QuStateEvolution{QPM<:QuPropagatorMethod, QV<:SparseVector, COT<:SparseMatrixCSC}(hamiltonian::SparseMatrixCSC, collapse_ops::Vector{COT}, init_state::QV, tlist, method::QPM) = QuStateEvolution{QPM,QV,QuLindbladMasterEqSparse}(QuLindbladMasterEqUncachedSparse(hamiltonian,collapse_ops), SparseMatrixCSC(init_state), tlist, method)
 
 # QuPropagator and QuStateEvolution types are the same.
 typealias QuPropagator QuStateEvolution
@@ -169,9 +202,13 @@ Inputs :
 Output :
 * Evolved operator.
 """ ->
-QuEvolutionOp{QM<:QuBase.AbstractQuMatrix}(op::QM, dt::Float64) = expm(-im*op*dt)
+QuEvolutionOp{QM<:AbstractMatrix}(op::QM, dt::Float64) = expm(-im*op*dt)
+
+QuEvolutionOp{QM<:SparseMatrixCSC}(op::QM, dt::Float64) = expm(-im*op*dt)
+
+QuEvolutionOp{QM<:AbstractMatrix}(op::QM, tf::Float64,  ti::Float64) = QuEvolutionOp(op, tf-ti)
 
-QuEvolutionOp{QM<:QuBase.AbstractQuMatrix}(op::QM, tf::Float64,  ti::Float64) = QuEvolutionOp(op, tf-ti)
+QuEvolutionOp{QM<:SparseMatrixCSC}(op::QM, tf::Float64,  ti::Float64) = QuEvolutionOp(op, tf-ti)
 
 QuEvolutionOp{QE<:QuEquation}(eq::QE, dt::Float64) = QuEvolutionOp(operator(eq), dt)
 
@@ -182,18 +219,18 @@ function Base.show(io::IO, qprop::QuPropagator)
     println(io, "Summarizing the system :")
     if :lindblad in field_params && :hamiltonian in field_params
         println(io, "Equation type : $(typeof(qprop.eq))")
-        println(io, "Size of the Lindblad operator of the system : $(size(coeffs(qprop.eq.lindblad)))")
-        println(io, "Size of the Hamiltonian of the system : $(size(coeffs(qprop.eq.hamiltonian)))")
+        println(io, "Size of the Lindblad operator of the system : $(size(qprop.eq.lindblad))")
+        println(io, "Size of the Hamiltonian of the system : $(size(qprop.eq.hamiltonian))")
         println(io, "Number of collapse operators : $(length(qprop.eq.collapse_ops))")
-        println(io, "Size of the Density matrix : $(size(coeffs(qprop.init_state)))")
+        println(io, "Size of the Density matrix : $(size(qprop.init_state))")
     elseif :hamiltonian in field_params
         println(io, "Equation type : $(typeof(qprop.eq))")
-        println(io, "Size of the Hamiltonian of the system : $(size(coeffs(qprop.eq.hamiltonian)))")
-        println(io, "Size of the Initial state : $(size(coeffs(qprop.init_state)))")
+        println(io, "Size of the Hamiltonian of the system : $(size(qprop.eq.hamiltonian))")
+        println(io, "Size of the Initial state : $(size(qprop.init_state))")
     elseif :liouvillian in field_params
         println(io, "Equation type : $(typeof(qprop.eq))")
-        println(io, "Size of the Liouvillian of the system : $(size(coeffs(qprop.eq.liouvillian)))")
-        println(io, "Size of the Density matrix : $(size(coeffs(qprop.init_state)))")
+        println(io, "Size of the Liouvillian of the system : $(size(qprop.eq.liouvillian))")
+        println(io, "Size of the Density matrix : $(size(qprop.init_state))")
     end
     println(io, "Time steps used : $(qprop.tlist)")
     println(io, "Solver used : $(qprop.method)")

src/propmcwf.jl

@@ -37,15 +37,16 @@ Ensemble of state, number of trajectories, decomposition based on the state.
 
   Decomposition is of state `rho` if `rho` is a `QuBase.AbstractQuMatrix`.
 """ ->
-type QuMCWFEnsemble{QA<:QuBase.AbstractQuArray}
+type QuMCWFEnsemble{QA<:Union{AbstractVector, AbstractMatrix, SparseVector, SparseMatrixCSC}}
     rho::QA
     ntraj::Int
     decomp
     QuMCWFEnsemble(rho, ntraj, decomp) = new(rho, ntraj, decomp)
 end
 
-QuMCWFEnsemble{QM<:QuBase.AbstractQuMatrix}(rho::QM, ntraj=500) = QuMCWFEnsemble{QM}(rho, ntraj, eigfact(full(coeffs(rho))))
-QuMCWFEnsemble{QV<:QuBase.AbstractQuVector}(psi::QV, ntraj=500) = QuMCWFEnsemble{QV}(psi, ntraj, nothing)
+QuMCWFEnsemble{QM<:Union{AbstractMatrix, SparseMatrixCSC}}(rho::QM, ntraj=500) = QuMCWFEnsemble{QM}(rho, ntraj, eigfact(full(rho)))
+QuMCWFEnsemble{QV<:AbstractVector}(psi::QV, ntraj=500) = QuMCWFEnsemble{QV}(psi, ntraj, nothing)
+QuMCWFEnsemble{QV<:SparseVector}(psi::QV, ntraj=500) = QuMCWFEnsemble{QV}(SparseMatrixCSC(psi), ntraj, nothing)
 
 Base.length(en::QuMCWFEnsemble) = en.ntraj
 
@@ -117,7 +118,7 @@ Inputs :
 Output :
 * QuArray with eigen vector from decomposition passed as an argument.
 """ ->
-function draw{QM<:QuBase.AbstractQuMatrix}(mcwfensemble::QuMCWFEnsemble{QM})
+function draw{QM<:Union{AbstractMatrix, SparseMatrixCSC}}(mcwfensemble::QuMCWFEnsemble{QM})
     r = rand() # draw random number from [0,1)
     pacc = 0.
     for i=1:length(mcwfensemble.decomp[:values])
@@ -141,11 +142,11 @@ Inputs :
 Output :
 * Copy of the QuVector.
 """ ->
-function draw{QV<:QuBase.AbstractQuVector}(mcwfensemble::QuMCWFEnsemble{QV})
+function draw{QV<:Union{AbstractVector, SparseVector}}(mcwfensemble::QuMCWFEnsemble{QV})
     return copy(mcwfensemble.rho)
 end
 
-function propagate(prob::QuMCWF, eq::QuLindbladMasterEq, t_end, t_start, psi)
+function propagate(prob::QuMCWF, eq::Union{QuLindbladMasterEq, QuLindbladMasterEqSparse}, t_end, t_start, psi)
     const jtol = 1.e-6  # jump tolerance
     # get information of QME
     heff = eff_hamiltonian(eq)
@@ -190,5 +191,9 @@ function propagate(prob::QuMCWF, eq::QuLindbladMasterEq, t_end, t_start, psi)
     return psi
 end
 
+propagate(prob::QuMCWF, eq::QuLindbladMasterEq, t_end, t_start, psi::AbstractVector) = propagae(prob, eq, t_end, t_start, psi)
+
+propagate(prob::QuMCWF, eq::QuLindbladMasterEqSparse, t_end, t_start, psi::SparseVector) = propagate(prob, eq, t_end, t_start, SparseMatrixCSC(psi))
+
 export QuMCWFEnsemble,
       QuMCWF

src/propodesolvers.jl

@@ -31,10 +31,9 @@ for (qu_ode_type,ode_solver) in type_to_method_ode
             dims = size(current_qustate)
             # Convert the current_qustate to complex as it might result in a Inexact Error. After complex is in QuBase.jl (PR #38)
             # we could just do a complex(vec(current_qustate)) avoiding the coeffs(coeffs(vec(current_qustate))).
-            next_state = $ode_solver((t,y)-> -im*coeffs(op)*y, complex(coeffs(vec(current_qustate))), [current_t, t], points=:specified,
+            next_state = $ode_solver((t,y)-> -im*op*y, complex(vec(current_qustate)), [current_t, t], points=:specified,
                                   reltol = get(prob.options, :reltol, 1.0e-5), abstol = get(prob.options, :abstol, 1.0e-8))[2][end]
-            CQST = QuBase.similar_type(current_qustate)
-            return CQST(reshape(next_state, dims), bases(current_qustate))
+            return reshape(next_state, dims)
         end
     end
 end

src/propstepsolvers.jl

@@ -65,8 +65,7 @@ function propagate(prob::QuEuler, eq::QuEquation, t, current_t, current_qustate)
     dims = size(current_qustate)
     op = operator(eq, t)
     next_state = (eye(op)-im*op*dt)*vec(current_qustate)
-    CQST = QuBase.similar_type(current_qustate)
-    return CQST(reshape(coeffs(next_state), dims), bases(current_qustate))
+    return reshape(next_state, dims)
 end
 
 function propagate(prob::QuCrankNicolson, eq::QuEquation, t, current_t, current_qustate)
@@ -75,8 +74,7 @@ function propagate(prob::QuCrankNicolson, eq::QuEquation, t, current_t, current_
     op = operator(eq, t)
     uni = eye(op)-im*op*dt/2
     next_state = \(uni', uni*vec(current_qustate))
-    CQST = QuBase.similar_type(current_qustate)
-    return CQST(reshape(coeffs(next_state), dims), bases(current_qustate))
+    return reshape(next_state, dims)
 end
 
 function propagate(prob::QuKrylov, eq::QuEquation, t, current_t, current_qustate)
@@ -85,18 +83,23 @@ function propagate(prob::QuKrylov, eq::QuEquation, t, current_t, current_qustate
     cqs = vec(current_qustate)
     basis_size = get(prob.options,:NC, length(cqs))
     N = min(basis_size, length(cqs))
-    v = Array(typeof(cqs),0)
+    v = Array(Any,0)
     @compat sizehint!(v, N+1)
-    push!(v,zeros(cqs))
+    push!(v,zeros(Complex{Float64}, cqs))
     push!(v,cqs)
+    println(v)
     alpha = Array(Complex{Float64},0)
     @compat sizehint!(alpha, N)
+    println(alpha)
     beta = Array(Complex{Float64},0)
     @compat sizehint!(beta, N+1)
     op = operator(eq, t)
     push!(beta,0.)
     for i=2:N
+        println(i)
         w = op*v[i]
+        println(w)
+        println("alpha",alpha)
         push!(alpha, w'*v[i])
         w = w-alpha[i-1]*v[i]-beta[i-1]*v[i-1]
         push!(beta, norm(w))
@@ -111,8 +114,7 @@ function propagate(prob::QuKrylov, eq::QuEquation, t, current_t, current_qustate
     for j=1:N
         next_state = next_state + v[j]*U_k[j,1]
     end
-    CQST = QuBase.similar_type(current_qustate)
-    return CQST(reshape(coeffs(next_state), dims), bases(current_qustate))
+    return reshape(next_state, dims)
 end
 
 export QuEuler,