Merge pull request #121 from ytdHuang/dev/SciML

Make HEOMLS data be `SciMLOperator`

Project.toml

@@ -1,7 +1,7 @@
 name = "HierarchicalEOM"
 uuid = "a62dbcb7-80f5-4d31-9a88-8b19fd92b128"
 authors = ["Yi-Te Huang <[email protected]>"]
-version = "2.2.4"
+version = "2.3.0"
 
 [deps]
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
@@ -37,7 +37,7 @@ LinearSolve = "2.4.2 - 2"
 OrdinaryDiffEqCore = "1"
 OrdinaryDiffEqLowOrderRK = "1"
 Pkg = "1"
-QuantumToolbox = "0.15 - 0.21"
+QuantumToolbox = "0.22"
 Reexport = "1"
 SciMLBase = "2"
 SciMLOperators = "0.3"

docs/src/libraryAPI.md

@@ -60,9 +60,10 @@ ODD
 ## HEOM Liouvillian superoperator matrices
 ```@docs
 HEOMSuperOp
-HEOMSuperOp(op, opParity::AbstractParity, refHEOMLS::AbstractHEOMLSMatrix, mul_basis::AbstractString="L")
-HEOMSuperOp(op, opParity::AbstractParity, refADOs::ADOs, mul_basis::AbstractString="L")
-HEOMSuperOp(op, opParity::AbstractParity, dims::SVector, N::Int, mul_basis::AbstractString)
+HEOMSuperOp(op, opParity::AbstractParity, refHEOMLS::AbstractHEOMLSMatrix)
+HEOMSuperOp(op, opParity::AbstractParity, refADOs::ADOs)
+HEOMSuperOp(op, opParity::AbstractParity, dims::SVector, N::Int)
+AbstractHEOMLSMatrix
 M_S
 M_S(Hsys::QuantumObject, parity::AbstractParity=EVEN; verbose::Bool=true)
 M_Boson

docs/src/spectrum.md

@@ -71,10 +71,7 @@ Finially, one can obtain the value of the power spectrum for specific ``\omega``
 !!! note "Odd-Parity for Power Spectrum"
     When ``Q`` is an operator acting on fermionic systems and has `ODD`-parity, the HEOMLS matrix ``\hat{\mathcal{M}}`` is acting on the `ODD`-parity space because ``\textbf{b}=Q\rho^{(m,n,+)}_{\textbf{j} \vert \textbf{q}}``. Therefore, remember to construct ``\hat{\mathcal{M}}`` with `ODD` [parity](@ref doc-Parity) in this kind of cases.
 
-See also the docstring : 
-```@docs
-PowerSpectrum(M::AbstractHEOMLSMatrix, ρ::Union{QuantumObject,ADOs}, P_op::Union{QuantumObject,HEOMSuperOp}, Q_op::Union{QuantumObject,HEOMSuperOp}, ωlist::AbstractVector, reverse::Bool = false; solver = UMFPACKFactorization(), verbose::Bool = true, filename::String = "", SOLVEROptions...)
-```
+See also the docstring [`PowerSpectrum`](@ref).
 
 ```julia
 M::AbstractHEOMLSMatrix
@@ -121,10 +118,7 @@ Finially, one can obtain the density of states for specific ``\omega``, namely
 !!! note "Odd-Parity for Density of States"
     As shown above, the HEOMLS matrix ``\hat{\mathcal{M}}`` acts on the `ODD`-parity space, compatibly with the parity of both the operators ``\textbf{b}_-=d\rho^{(m,n,+)}_{\textbf{j} \vert \textbf{q}}`` and ``\textbf{b}_+=d^\dagger\rho^{(m,n,+)}_{\textbf{j} \vert \textbf{q}}``. Therefore, remember to construct ``\hat{\mathcal{M}}`` with `ODD` [parity](@ref doc-Parity) for solving spectrum of fermionic systems.
 
-See also the docstring : 
-```@docs
-DensityOfStates(M::AbstractHEOMLSMatrix, ρ::QuantumObject, d_op::QuantumObject, ωlist::AbstractVector; solver=UMFPACKFactorization(), verbose::Bool = true, filename::String = "", SOLVEROptions...)
-```
+See also the docstring [`DensityOfStates`](@ref).
 
 ```julia
 Hs::QuantumObject  # system Hamiltonian

docs/src/stationary_state.md

@@ -11,12 +11,6 @@ The first method is implemented by solving the linear problem
 
 `HierarchicalEOM.jl` wraps some of the functions in [LinearSolve.jl](http://linearsolve.sciml.ai/stable/), which is a very rich numerical library for solving the linear problems and provides many solvers. It offers quite a few options for the user to tailor the solver to their specific needs. The default solver (and its corresponding settings) are chosen to suit commonly encountered problems and should work fine for most of the cases. If you require more specialized methods, such as the choice of algorithm, please refer to [LinearSolve solvers](@ref LS-solvers) and also the documentation of [LinearSolve.jl](http://linearsolve.sciml.ai/stable/).
 
-See the docstring of this method:  
-
-```@docs
-steadystate(M::AbstractHEOMLSMatrix; solver=UMFPACKFactorization(), verbose::Bool=true, SOLVEROptions...)
-```
-
 ```julia
 # the HEOMLS matrix
 M::AbstractHEOMLSMatrix  
@@ -34,13 +28,7 @@ until finding a stationary solution.
 
 `HierarchicalEOM.jl` wraps some of the functions in [DifferentialEquations.jl](https://diffeq.sciml.ai/stable/), which is a very rich numerical library for solving the differential equations and provides many ODE solvers. It offers quite a few options for the user to tailor the solver to their specific needs. The default solver (and its corresponding settings) are chosen to suit commonly encountered problems and should work fine for most of the cases. If you require more specialized methods, such as the choice of algorithm, please refer to the documentation of [DifferentialEquations.jl](https://diffeq.sciml.ai/stable/).
 
-
 ### Given the initial state as Density Operator (`QuantumObject` type)
-See the docstring of this method:  
-
-```@docs
-steadystate(M::AbstractHEOMLSMatrix, ρ0::QuantumObject, tspan::Number = Inf; solver = DP5(), termination_condition = NormTerminationMode(), verbose::Bool = true, SOLVEROptions...)
-```
 
 ```julia
 # the HEOMLS matrix
@@ -53,10 +41,6 @@ ados_steady = steadystate(M, ρ0)
 ```
 
 ### Given the initial state as Auxiliary Density Operators
-See the docstring of this method:  
-```@docs
-steadystate(M::AbstractHEOMLSMatrix, ados::ADOs, tspan::Number = Inf; solver = DP5(), termination_condition = NormTerminationMode(), verbose::Bool = true, SOLVEROptions...)
-```
 
 ```julia
 # the HEOMLS matrix

docs/src/time_evolution.md

@@ -116,15 +116,15 @@ M = M_Fermion(H0, ...)
 M = M_BosonFermion(H0, ...)
 ```
 
-To solve the dynamics characterized by ``\hat{\mathcal{M}}`` together with the time-dependent part of system Hamiltonian ``H_1(t)``, you can specify keyword arguments `H_t` and `params` while calling [`HEOMsolve`](@ref). Here, the definition of user-defined function `H_1` must be in the form `H_1(t, params::NamedTuple)` and returns the time-dependent part of system Hamiltonian or Liouvillian (in `QuantumObject` type) at any given time point `t`. The parameters `params` should be a `NamedTuple` which contains all the extra parameters you need for the function `H_1`. For example:
+To solve the dynamics characterized by ``\hat{\mathcal{M}}`` together with the time-dependent part of system Hamiltonian ``H_1(t)``, you can specify keyword arguments `H_t` and `params` while calling [`HEOMsolve`](@ref). Here, `H_t` must be specified as a `QuantumToolbox.QuantumObjectEvolution` (or `QobjEvo`), and `params` should contain all the extra parameters you need for `QobjEvo`, for example:
 ```julia
-# in this case, p should be passed in as a NamedTuple: (p0 = p0, p1 = p1, p2 = p2)
-function H_1(t, p::NamedTuple)  
-    σx = [0 1; 1 0] # Pauli-X matrix
-    return (sin(p.p0 * t) + sin(p.p1 * t) + sin(p.p2 * t)) * σx
-end
+# in this case, p will be passed in as a NamedTuple: (p0 = p0, p1 = p1, p2 = p2)
+coef(p::NamedTuple, t) = sin(p.p0 * t) + sin(p.p1 * t) + sin(p.p2 * t)
+
+σx = sigmax() # Pauli-X matrix
+H_1 = QobjEvo(σx, coef)
 ```
-The `p` will be passed to your function `H_1` directly from the keyword argument in [`HEOMsolve`](@ref) called `params`:
+The `p` can be passed to `H_1` directly from the keyword argument in [`HEOMsolve`](@ref) called `params`:
 ```julia
 M::AbstractHEOMLSMatrix
 ρ0::QuantumObject
@@ -134,24 +134,6 @@ p = (p0 = 0.1, p1 = 1, p2 = 10)
 sol = HEOMsolve(M, ρ0, tlist; H_t = H_1, params = p)
 ```
 
-!!! warning "Warning"
-    If you don't need any extra `param` in your case, you still need to put a redundant one in the definition of `H_1`, for example:
-
-```julia
-function H_1(t, p::NamedTuple)
-    σx = [0 1; 1 0] # Pauli-X matrix
-    return sin(0.1 * t) * σx
-end
-
-M::AbstractHEOMLSMatrix
-ρ0::QuantumObject
-tlist = 0:0.1:10
-
-sol = HEOMsolve(M, ρ0, tlist; H_t = H_1)
-```
-!!! note "Note"
-    The default value for `params` in `HEOMsolve` is an empty `NamedTuple()`.
-
 ## Propagator Method
 The second method is implemented by directly construct the propagator of a given [HEOMLS matrix](@ref doc-HEOMLS-Matrix) ``\hat{\mathcal{M}}``. Because ``\hat{\mathcal{M}}`` is time-independent, the equation above can be solved analytically as
 ```math

examples/dynamical_decoupling.jl

@@ -52,12 +52,12 @@ H0 = 0.5 * ω0 * σz
 # Here, we set the period of the pulses to be $\tau V = \pi/2$. Therefore, we consider two scenarios with fast and slow pulses:
 
 ## a function which returns the amplitude of the pulse at time t
-function pulse(V, Δ, t)
-    τ = 0.5 * π / V
-    period = τ + Δ
+function pulse(p, t)
+    τ = 0.5 * π / p.V
+    period = τ + p.Δ
 
     if (t % period) < τ
-        return V
+        return p.V
     else
         return 0
     end
@@ -68,9 +68,12 @@ amp_fast = 0.50
 amp_slow = 0.01
 delay = 20
 
+fastTuple = (V = amp_fast, Δ = delay)
+slowTuple = (V = amp_slow, Δ = delay)
+
 Plots.plot(
     tlist,
-    [[pulse(amp_fast, delay, t) for t in tlist], [pulse(amp_slow, delay, t) for t in tlist]],
+    [[pulse(fastTuple, t) for t in tlist], [pulse(slowTuple, t) for t in tlist]],
     label = ["Fast Pulse" "Slow Pulse"],
     linestyle = [:solid :dash],
 )
@@ -101,13 +104,10 @@ noPulseSol = HEOMsolve(M, ψ0, tlist; e_ops = [ρ01]);
 # ## Solve time evolution with time-dependent Hamiltonian
 # (see also [Time Evolution](@ref doc-Time-Evolution))
 #   
-# We need to provide a user-defined function (named as `H_D` in this case), which must be in the form `H_D(t, p::NamedTuple)` and returns the time-dependent part of system Hamiltonian (in `QuantumObject` type) at any given time point `t`. The parameter `p` should be a `NamedTuple` which contains all the extra parameters [`V` (amplitude), `Δ` (delay), and `σx` (operator) in this case] for the function `H_D`:
-H_D(t, p::NamedTuple) = pulse(p.V, p.Δ, t) * p.σx;
-
-# The parameter `p` will be passed to your function `H_D` directly from the **last required** parameter in `HEOMsolve`:
-fastTuple = (V = amp_fast, Δ = delay, σx = σx)
-slowTuple = (V = amp_slow, Δ = delay, σx = σx)
+# We need to provide a `QuantumToolbox.QuantumObjectEvolution` (named as `H_D` in this case)
+H_D = QobjEvo(σx, pulse)
 
+# The keyword argument `params` in `HEOMsolve` will be passed to the argument `p` in user-defined function (`pulse` in this case) directly:
 fastPulseSol = HEOMsolve(M, ψ0, tlist; e_ops = [ρ01], H_t = H_D, params = fastTuple)
 slowPulseSol = HEOMsolve(M, ψ0, tlist; e_ops = [ρ01], H_t = H_D, params = slowTuple)
 

ext/HierarchicalEOM_CUDAExt.jl

@@ -7,6 +7,7 @@ import CUDA
 import CUDA: cu, CuArray
 import CUDA.CUSPARSE: CuSparseVector, CuSparseMatrixCSC
 import SparseArrays: sparse, SparseVector, SparseMatrixCSC
+import SciMLOperators: MatrixOperator, ScaledOperator, AddedOperator
 
 @doc raw"""
     cu(M::AbstractHEOMLSMatrix)
@@ -19,63 +20,42 @@ cu(M::AbstractHEOMLSMatrix) = CuSparseMatrixCSC(M)
 Return a new HEOMLS-matrix-type object with `M.data` is in the type of `CuSparseMatrixCSC{ComplexF32, Int32}` for gpu calculations.
 """
 function CuSparseMatrixCSC(M::T) where {T<:AbstractHEOMLSMatrix}
-    A = M.data
-    if typeof(A) <: CuSparseMatrixCSC
-        return M
+    A_gpu = _convert_to_gpu_matrix(M.data)
+    if T <: M_S
+        return M_S(A_gpu, M.tier, M.dims, M.N, M.sup_dim, M.parity)
+    elseif T <: M_Boson
+        return M_Boson(A_gpu, M.tier, M.dims, M.N, M.sup_dim, M.parity, M.bath, M.hierarchy)
+    elseif T <: M_Fermion
+        return M_Fermion(A_gpu, M.tier, M.dims, M.N, M.sup_dim, M.parity, M.bath, M.hierarchy)
     else
-        colptr = CuArray{Int32}(A.colptr)
-        rowval = CuArray{Int32}(A.rowval)
-        nzval = CuArray{ComplexF32}(A.nzval)
-        A_gpu = CuSparseMatrixCSC{ComplexF32,Int32}(colptr, rowval, nzval, size(A))
-        if T <: M_S
-            return M_S(A_gpu, M.tier, M.dims, M.N, M.sup_dim, M.parity)
-        elseif T <: M_Boson
-            return M_Boson(A_gpu, M.tier, M.dims, M.N, M.sup_dim, M.parity, M.bath, M.hierarchy)
-        elseif T <: M_Fermion
-            return M_Fermion(A_gpu, M.tier, M.dims, M.N, M.sup_dim, M.parity, M.bath, M.hierarchy)
-        else
-            return M_Boson_Fermion(
-                A_gpu,
-                M.Btier,
-                M.Ftier,
-                M.dims,
-                M.N,
-                M.sup_dim,
-                M.parity,
-                M.Bbath,
-                M.Fbath,
-                M.hierarchy,
-            )
-        end
+        return M_Boson_Fermion(A_gpu, M.Btier, M.Ftier, M.dims, M.N, M.sup_dim, M.parity, M.Bbath, M.Fbath, M.hierarchy)
     end
 end
-function CuSparseMatrixCSC{ComplexF32}(M::HEOMSuperOp)
-    A = M.data
-    AType = typeof(A)
-    if AType == CuSparseMatrixCSC{ComplexF32,Int32}
-        return M
-    elseif AType <: CuSparseMatrixCSC
+
+CuSparseMatrixCSC{ComplexF32}(M::HEOMSuperOp) = HEOMSuperOp(_convert_to_gpu_matrix(M.data), M.dims, M.N, M.parity)
+
+function _convert_to_gpu_matrix(A::AbstractSparseMatrix)
+    if A isa CuSparseMatrixCSC{ComplexF32,Int32}
+        return A
+    elseif A isa CuSparseMatrixCSC
         colptr = CuArray{Int32}(A.colPtr)
         rowval = CuArray{Int32}(A.rowVal)
         nzval = CuArray{ComplexF32}(A.nzVal)
-        A_gpu = CuSparseMatrixCSC{ComplexF32,Int32}(colptr, rowval, nzval, size(A))
-        return HEOMSuperOp(A_gpu, M.dims, M.N, M.parity)
+        return CuSparseMatrixCSC{ComplexF32,Int32}(colptr, rowval, nzval, size(A))
     else
         colptr = CuArray{Int32}(A.colptr)
         rowval = CuArray{Int32}(A.rowval)
         nzval = CuArray{ComplexF32}(A.nzval)
-        A_gpu = CuSparseMatrixCSC{ComplexF32,Int32}(colptr, rowval, nzval, size(A))
-        return HEOMSuperOp(A_gpu, M.dims, M.N, M.parity)
+        return CuSparseMatrixCSC{ComplexF32,Int32}(colptr, rowval, nzval, size(A))
     end
 end
+_convert_to_gpu_matrix(A::MatrixOperator) = MatrixOperator(_convert_to_gpu_matrix(A.A))
+_convert_to_gpu_matrix(A::ScaledOperator) = ScaledOperator(A.λ, _convert_to_gpu_matrix(A.L))
+_convert_to_gpu_matrix(A::AddedOperator) = AddedOperator(map(op -> _convert_to_gpu_matrix(op), A.ops))
 
-function _Tr(M::AbstractHEOMLSMatrix{T}) where {T<:CuSparseMatrixCSC}
-    D = prod(M.dims)
-    return CuSparseVector(SparseVector(M.N * D^2, [1 + n * (D + 1) for n in 0:(D-1)], ones(eltype(M), D)))
-end
+_Tr(M::Type{<:CuSparseMatrixCSC}, dims::SVector, N::Int) = CuSparseVector(_Tr(eltype(M), dims, N))
 
 # change the type of `ADOs` to match the type of HEOMLS matrix
-_HandleVectorType(::AbstractHEOMLSMatrix{<:CuSparseMatrixCSC{T}}, V::SparseVector) where {T<:Number} =
-    CuArray{_CType(T)}(V)
+_HandleVectorType(M::Type{<:CuSparseMatrixCSC}, V::SparseVector) = CuArray{_CType(eltype(M))}(V)
 
 end

src/ADOs.jl

@@ -166,9 +166,9 @@ where ``O`` is the operator and ``\rho`` is the reduced density operator in the
 - `exp_val` : The expectation value
 """
 function QuantumToolbox.expect(op, ados::ADOs; take_real::Bool = true)
-    if typeof(op) <: HEOMSuperOp
+    if op isa HEOMSuperOp
         _check_sys_dim_and_ADOs_num(op, ados)
-        exp_val = dot(transpose(_Tr(ados.dims, ados.N)), (SparseMatrixCSC(op) * ados).data)
+        exp_val = dot(transpose(_Tr(eltype(ados), ados.dims, ados.N)), (SparseMatrixCSC(op) * ados).data)
     else
         _op = HandleMatrixType(op, ados.dims, "op (observable)"; type = Operator)
         exp_val = tr(_op.data * getRho(ados).data)
@@ -204,13 +204,13 @@ function QuantumToolbox.expect(op, ados_list::Vector{ADOs}; take_real::Bool = tr
         _check_sys_dim_and_ADOs_num(ados_list[1], ados_list[i])
     end
 
-    if typeof(op) <: HEOMSuperOp
+    if op isa HEOMSuperOp
         _check_sys_dim_and_ADOs_num(op, ados_list[1])
         _op = op
     else
-        _op = HEOMSuperOp(op, EVEN, dims, N, "L")
+        _op = HEOMSuperOp(spre(op), EVEN, dims, N)
     end
-    tr_op = transpose(_Tr(dims, N)) * SparseMatrixCSC(_op).data
+    tr_op = transpose(_Tr(eltype(op), dims, N)) * SparseMatrixCSC(_op).data
 
     exp_val = [dot(tr_op, ados.data) for ados in ados_list]
 

src/HeomBase.jl

@@ -2,27 +2,56 @@ export AbstractHEOMLSMatrix
 
 abstract type AbstractHEOMLSMatrix{T} end
 
-QuantumToolbox._FType(::AbstractHEOMLSMatrix{<:AbstractArray{T}}) where {T<:Number} = _FType(T)
-QuantumToolbox._CType(::AbstractHEOMLSMatrix{<:AbstractArray{T}}) where {T<:Number} = _CType(T)
+@doc raw"""
+    (M::AbstractHEOMLSMatrix)(p, t)
+
+Calculate the time-dependent AbstractHEOMLSMatrix at time `t` with parameters `p`.
+
+# Arguments
+- `p`: The parameters of the time-dependent coefficients.
+- `t`: The time at which the coefficients are evaluated.
+
+# Returns
+- The output HEOMLS matrix.
+"""
+function (M::AbstractHEOMLSMatrix)(p, t)
+    # We put 0 in the place of `u` because the time-dependence doesn't depend on the state
+    update_coefficients!(M.data, 0, p, t)
+    return concretize(M.data)
+end
+
+(M::AbstractHEOMLSMatrix)(t) = M(nothing, t)
+
+QuantumToolbox._FType(M::AbstractHEOMLSMatrix) = _FType(eltype(M))
+QuantumToolbox._CType(M::AbstractHEOMLSMatrix) = _CType(eltype(M))
+
+_get_SciML_matrix_wrapper(M::AbstractArray) = QuantumToolbox.get_typename_wrapper(M){eltype(M)}
+_get_SciML_matrix_wrapper(M::MatrixOperator) = _get_SciML_matrix_wrapper(M.A)
+_get_SciML_matrix_wrapper(M::ScaledOperator) = _get_SciML_matrix_wrapper(M.L)
+_get_SciML_matrix_wrapper(M::AddedOperator) = _get_SciML_matrix_wrapper(M.ops[1])
+_get_SciML_matrix_wrapper(M::AbstractHEOMLSMatrix) = _get_SciML_matrix_wrapper(M.data)
 
 # equal to : sparse(vec(system_identity_matrix))
-function _Tr(dims::SVector, N::Int)
+function _Tr(T::Type{<:Number}, dims::SVector, N::Int)
     D = prod(dims)
-    return SparseVector(N * D^2, [1 + n * (D + 1) for n in 0:(D-1)], ones(ComplexF64, D))
-end
-function _Tr(M::AbstractHEOMLSMatrix{T}) where {T<:SparseMatrixCSC}
-    D = prod(M.dims)
-    return SparseVector(M.N * D^2, [1 + n * (D + 1) for n in 0:(D-1)], ones(eltype(M), D))
+    return SparseVector(N * D^2, [1 + n * (D + 1) for n in 0:(D-1)], ones(T, D))
 end
+_Tr(M::AbstractHEOMLSMatrix) = _Tr(_get_SciML_matrix_wrapper(M), M.dims, M.N)
+_Tr(M::Type{<:SparseMatrixCSC}, dims::SVector, N::Int) = _Tr(eltype(M), dims, N)
 
-function HandleMatrixType(M::QuantumObject, MatrixName::String = ""; type::QuantumObjectType = Operator)
+function HandleMatrixType(M::AbstractQuantumObject, MatrixName::String = ""; type::QuantumObjectType = Operator)
     if M.type == type
         return M
     else
         error("The matrix $(MatrixName) should be an $(type).")
     end
 end
-function HandleMatrixType(M::QuantumObject, dims::SVector, MatrixName::String = ""; type::QuantumObjectType = Operator)
+function HandleMatrixType(
+    M::AbstractQuantumObject,
+    dims::SVector,
+    MatrixName::String = "";
+    type::QuantumObjectType = Operator,
+)
     if M.dims == dims
         return HandleMatrixType(M, MatrixName; type = type)
     else
@@ -35,13 +64,14 @@ HandleMatrixType(M, MatrixName::String = ""; type::QuantumObjectType = Operator)
     error("HierarchicalEOM doesn't support matrix $(MatrixName) with type : $(typeof(M))")
 
 # change the type of `ADOs` to match the type of HEOMLS matrix
-_HandleVectorType(::AbstractHEOMLSMatrix{<:AbstractSparseMatrix{T}}, V::SparseVector) where {T<:Number} =
-    Vector{_CType(T)}(V)
+_HandleVectorType(M::AbstractHEOMLSMatrix, V::SparseVector) = _HandleVectorType(_get_SciML_matrix_wrapper(M), V)
+_HandleVectorType(M::Type{<:SparseMatrixCSC}, V::SparseVector) = Vector{_CType(eltype(M))}(V)
 
-function _HandleSteadyStateMatrix(M::AbstractHEOMLSMatrix, S::Int)
+function _HandleSteadyStateMatrix(M::AbstractHEOMLSMatrix{<:MatrixOperator})
+    S = size(M, 1)
     ElType = eltype(M)
     D = prod(M.dims)
-    A = copy(M.data)
+    A = copy(M.data.A)
     A[1, 1:S] .= 0
 
     # sparse(row_idx, col_idx, values, row_dims, col_dims)