diff --git a/src/QuDynamics.jl b/src/QuDynamics.jl index a2e1840..21d7e91 100644 --- a/src/QuDynamics.jl +++ b/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 diff --git a/src/propexpmsolvers.jl b/src/propexpmsolvers.jl index 881258d..462efc2 100644 --- a/src/propexpmsolvers.jl +++ b/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, diff --git a/src/propmachinery.jl b/src/propmachinery.jl index 2b089f6..0e3cfa8 100644 --- a/src/propmachinery.jl +++ b/src/propmachinery.jl @@ -27,7 +27,7 @@ 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 @@ -35,24 +35,57 @@ immutable QuStateEvolution{QPM<:QuPropagatorMethod, QVM<:@compat(Union{QuBase.Ab 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)") diff --git a/src/propmcwf.jl b/src/propmcwf.jl index 61b99db..36b4004 100644 --- a/src/propmcwf.jl +++ b/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 diff --git a/src/propodesolvers.jl b/src/propodesolvers.jl index 082ed5b..661e5c3 100644 --- a/src/propodesolvers.jl +++ b/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 diff --git a/src/propstepsolvers.jl b/src/propstepsolvers.jl index 6100e39..8978daa 100644 --- a/src/propstepsolvers.jl +++ b/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, diff --git a/src/quequations.jl b/src/quequations.jl index 54281dc..51e9796 100644 --- a/src/quequations.jl +++ b/src/quequations.jl @@ -12,7 +12,7 @@ Schrodinger Equation type Hamiltonian of the system is the only field for the type. """ -> -immutable QuSchrodingerEq{H<:QuBase.AbstractQuMatrix} <: QuEquation{1} +immutable QuSchrodingerEq{H<:AbstractMatrix} <: QuEquation{1} hamiltonian::H QuSchrodingerEq(hamiltonian) = new(hamiltonian) end @@ -30,7 +30,14 @@ Inputs : Output : * QuSchrodingerEq type construct. """ -> -QuSchrodingerEq{H<:QuBase.AbstractQuMatrix}(hamiltonian::H) = QuSchrodingerEq{H}(hamiltonian) +QuSchrodingerEq{H<:AbstractMatrix}(hamiltonian::H) = QuSchrodingerEq{H}(hamiltonian) + +immutable QuSchrodingerEqSparse{H<:SparseMatrixCSC} <: QuEquation{1} + hamiltonian::H + QuSchrodingerEqSparse(hamiltonian) = new(hamiltonian) +end + +QuSchrodingerEqSparse{H<:SparseMatrixCSC}(hamiltonian::H) = QuSchrodingerEqSparse{H}(hamiltonian) @doc """ Liouville von Neumann Equation type @@ -41,7 +48,7 @@ Liouville von Neumann Equation type Liouvillian of the system is the only field for the type. """ -> -immutable QuLiouvillevonNeumannEq{H<:QuBase.AbstractQuMatrix} <: QuEquation{1} +immutable QuLiouvillevonNeumannEq{H<:AbstractMatrix} <: QuEquation{1} liouvillian::H QuLiouvillevonNeumannEq(liouvillian) = new(liouvillian) end @@ -59,7 +66,14 @@ Inputs : Output : * QuLiouvillevonNeumannEq type construct. """ -> -QuLiouvillevonNeumannEq{H<:QuBase.AbstractQuMatrix}(liouvillian::H) = QuLiouvillevonNeumannEq{H}(liouvillian) +QuLiouvillevonNeumannEq{H<:AbstractMatrix}(liouvillian::H) = QuLiouvillevonNeumannEq{H}(liouvillian) + +immutable QuLiouvillevonNeumannEqSparse{H<:SparseMatrixCSC} <: QuEquation{1} + liouvillian::H + QuLiouvillevonNeumannEqSparse(liouvillian) = new(liouvillian) +end + +QuLiouvillevonNeumannEqSparse{H<:SparseMatrixCSC}(liouvillian::H) = QuLiouvillevonNeumannEqSparse{H}(liouvillian) @doc """ Lindblad Master Equation type @@ -78,7 +92,7 @@ Lindblad Master Equation type """ -> # immutable QuLindbladMasterEq{CF, L<:@compat(Union{QuBase.AbstractQuMatrix, Nothing}), H<:QuBase.AbstractQuMatrix, V<:QuBase.AbstractQuMatrix} <: QuEquation{CF} # (Use this for running the benchmarks on version < 0.4) -immutable QuLindbladMasterEq{CF, L<:@compat(Union{QuBase.AbstractQuMatrix, Void}), H<:QuBase.AbstractQuMatrix, V<:QuBase.AbstractQuMatrix} <: QuEquation{CF} +immutable QuLindbladMasterEq{CF, L<:Union{AbstractMatrix, Void}, H<:AbstractMatrix, V<:AbstractMatrix} <: QuEquation{CF} lindblad::L hamiltonian::H collapse_ops::Vector{V} @@ -102,7 +116,7 @@ Inputs : Output : * QuLindbladMasterEq type construct. """ -> -function QuLindbladMasterEq{H<:QuBase.AbstractQuMatrix, V<:QuBase.AbstractQuMatrix}(hamiltonian::H, collapse_ops::Vector{V}) +function QuLindbladMasterEq{H<:AbstractMatrix, V<:AbstractMatrix}(hamiltonian::H, collapse_ops::Vector{V}) lop = lindblad_op(hamiltonian, collapse_ops) QuLindbladMasterEq{1,typeof(lop),H,V}(lop, hamiltonian, collapse_ops) end @@ -126,10 +140,27 @@ Output : * QuLindbladMasterEq type construct. """ -> # QuLindbladMasterEq{0,Nothing,H,V}(nothing, hamiltonian, collapse_ops) (Use this for running the benchmarks on version < 0.4) -function QuLindbladMasterEqUncached{H<:QuBase.AbstractQuMatrix, V<:QuBase.AbstractQuMatrix}(hamiltonian::H, collapse_ops::Vector{V}) +function QuLindbladMasterEqUncached{H<:AbstractMatrix, V<:AbstractMatrix}(hamiltonian::H, collapse_ops::Vector{V}) QuLindbladMasterEq{0,Void,H,V}(nothing, hamiltonian, collapse_ops) end +immutable QuLindbladMasterEqSparse{CF, L<:Union{SparseMatrixCSC, Void}, H<:SparseMatrixCSC, V<:SparseMatrixCSC} <: QuEquation{CF} + lindblad::L + hamiltonian::H + collapse_ops::Vector{V} + QuLindbladMasterEqSparse(lindblad, hamiltonian, collapse_ops) = new(lindblad, hamiltonian, collapse_ops) +end + +function QuLindbladMasterEqSparse{H<:SparseMatrixCSC, V<:SparseMatrixCSC}(hamiltonian::H, collapse_ops::Vector{V}) + lop = lindblad_op(hamiltonian, collapse_ops) + QuLindbladMasterEqSparse{1,typeof(lop),H,V}(lop, hamiltonian, collapse_ops) +end + +# QuLindbladMasterEq{0,Nothing,H,V}(nothing, hamiltonian, collapse_ops) (Use this for running the benchmarks on version < 0.4) +function QuLindbladMasterEqUncachedSparse{H<:SparseMatrixCSC, V<:SparseMatrixCSC}(hamiltonian::H, collapse_ops::Vector{V}) + QuLindbladMasterEqSparse{0,Void,H,V}(nothing, hamiltonian, collapse_ops) +end + @doc """ Lindblad operator construct from the `Hamiltonian` and `collapse operators` @@ -146,12 +177,12 @@ Inputs : Output : * Lindblad operator. """ -> -function lindblad_op(hamiltonian::QuBase.AbstractQuMatrix, collapse_ops::Vector) - nb = size(coeffs(hamiltonian), 1) +function lindblad_op(hamiltonian::Union{AbstractMatrix, SparseMatrixCSC}, collapse_ops::Vector) + nb = size(hamiltonian, 1) nlop = length(collapse_ops) - heff = zeros(typeof(im*one(eltype(hamiltonian))), size(coeffs(hamiltonian))) + heff = zeros(typeof(im*one(eltype(hamiltonian))), size(hamiltonian)) for l=1:nlop - heff = heff + 0.5*coeffs(collapse_ops[l])'*coeffs(collapse_ops[l]) + heff = heff + 0.5*collapse_ops[l]'*collapse_ops[l] end SI = Array(Int,0) SJ = Array(Int,0) @@ -181,7 +212,7 @@ function lindblad_op(hamiltonian::QuBase.AbstractQuMatrix, collapse_ops::Vector) end end end - return QuBase.similar_type(hamiltonian)(sparse(SI, SJ, Lvals, nb*nb, nb*nb)) + return sparse(SI, SJ, Lvals, nb*nb, nb*nb) end @doc """ @@ -197,7 +228,7 @@ Inputs : Output : * Liouvillian operator. """ -> -liouvillian_op(hamiltonian::QuBase.AbstractQuMatrix) = lindblad_op(hamiltonian, []) +liouvillian_op(hamiltonian::Union{AbstractMatrix, SparseMatrixCSC}) = lindblad_op(hamiltonian, []) @doc """ Liouvillian of the QuLiouvillevonNeumannEq type. @@ -215,7 +246,7 @@ Inputs : Output : * Liouvillian of the system """ -> -function operator(qu_eq::QuLiouvillevonNeumannEq, t=0) +function operator(qu_eq::Union{QuLiouvillevonNeumannEq, QuLiouvillevonNeumannEqSparse}, t=0) return qu_eq.liouvillian end @@ -235,7 +266,7 @@ Inputs : Output : * Hamiltonian of the system """ -> -function operator(qu_eq::QuSchrodingerEq, t=0) +function operator(qu_eq::Union{QuSchrodingerEq, QuSchrodingerEqSparse}, t=0) return qu_eq.hamiltonian end @@ -256,7 +287,7 @@ Inputs : Output : * Lindblad operator of the system """ -> -function operator(qu_eq::QuLindbladMasterEq{1}, t=0) +function operator(qu_eq::Union{QuLindbladMasterEq{1}, QuLindbladMasterEqSparse{1}}, t=0) return qu_eq.lindblad end @@ -277,7 +308,7 @@ Inputs : Output : * Lindblad operator of the system """ -> -function operator(qu_eq::QuLindbladMasterEq{0}, t=0) +function operator(qu_eq::Union{QuLindbladMasterEq{0}, QuLindbladMasterEqSparse{0}}, t=0) return lindblad_op(qu_eq.hamiltonian, qu_eq.collapse_ops) end @@ -292,7 +323,7 @@ Inputs : Output * Effective hamiltonian of the system """ -> -function eff_hamiltonian(lme::QuLindbladMasterEq) +function eff_hamiltonian(lme::Union{QuLindbladMasterEq, QuLindbladMasterEqSparse}) heff = lme.hamiltonian for c_ops in lme.collapse_ops heff = heff - im*0.5*c_ops'*c_ops @@ -302,6 +333,10 @@ end export QuEquation, QuSchrodingerEq, + QuSchrodingerEqSparse, QuLiouvillevonNeumannEq, + QuLiouvillevonNeumannEqSparse, QuLindbladMasterEq, - QuLindbladMasterEqUncached + QuLindbladMasterEqSparse, + QuLindbladMasterEqUncached, + QuLindbladMasterEqUncachedSparse diff --git a/test/runtests.jl b/test/runtests.jl index e7a7aee..858d264 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,4 +4,4 @@ using Base.Test using Compat include("propagatortest.jl") -include("qutipinterfacetest.jl") +# include("qutipinterfacetest.jl")