diff --git a/CHANGELOG.md b/CHANGELOG.md
index e06ec247b..98a8398b3 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -9,6 +9,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 - Make CUDA conversion more general using Adapt.jl. ([#436], [#437])
 - Make the generation of `fock` states non-mutating to support Zygote.jl. ([#438])
+- Remove Reexport.jl from the dependencies. ([#443])
 
 ## [v0.29.1]
 Release date: 2025-03-07
@@ -193,3 +194,4 @@ Release date: 2024-11-13
 [#436]: https://github.com/qutip/QuantumToolbox.jl/issues/436
 [#437]: https://github.com/qutip/QuantumToolbox.jl/issues/437
 [#438]: https://github.com/qutip/QuantumToolbox.jl/issues/438
+[#443]: https://github.com/qutip/QuantumToolbox.jl/issues/443
diff --git a/Project.toml b/Project.toml
index fc5d42037..dbf8fe0ca 100644
--- a/Project.toml
+++ b/Project.toml
@@ -18,7 +18,6 @@ OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
 OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
-Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
@@ -56,7 +55,6 @@ OrdinaryDiffEqCore = "1"
 OrdinaryDiffEqTsit5 = "1"
 Pkg = "1"
 Random = "1"
-Reexport = "1"
 SciMLBase = "2"
 SciMLOperators = "0.3"
 SparseArrays = "1"
diff --git a/benchmarks/runbenchmarks.jl b/benchmarks/runbenchmarks.jl
index e921c81a4..a8419c516 100644
--- a/benchmarks/runbenchmarks.jl
+++ b/benchmarks/runbenchmarks.jl
@@ -1,4 +1,6 @@
 using BenchmarkTools
+using LinearAlgebra
+using SparseArrays
 using QuantumToolbox
 using OrdinaryDiffEq
 using LinearSolve
diff --git a/docs/make.jl b/docs/make.jl
index e624cbee3..2a6548979 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -10,7 +10,12 @@ using Changelog
 # Load of packages required to compile the extension documentation
 using CairoMakie
 
-DocMeta.setdocmeta!(QuantumToolbox, :DocTestSetup, :(using QuantumToolbox); recursive = true)
+doctest_setup = quote
+    using LinearAlgebra
+    using SparseArrays
+    using QuantumToolbox
+end
+DocMeta.setdocmeta!(QuantumToolbox, :DocTestSetup, doctest_setup; recursive = true)
 
 # some options for `makedocs`
 const DRAFT = false  # set `true`  to disable cell evaluation
diff --git a/docs/src/getting_started/type_stability.md b/docs/src/getting_started/type_stability.md
index 1daabb063..2b5153665 100644
--- a/docs/src/getting_started/type_stability.md
+++ b/docs/src/getting_started/type_stability.md
@@ -199,7 +199,7 @@ Which returns a tensor of size `2x2x2x2x2x2`. Let's check the `@code_warntype`:
 @code_warntype reshape_operator_data([2, 2, 2])
 ```
 
-We got a `Any` type, because the compiler doesn't know the size of the `dims` vector. We can fix this by using a `Tuple` (or `SVector`):
+We got a `Any` type, because the compiler doesn't know the size of the `dims` vector. We can fix this by using a `Tuple` (or `SVector` from [StaticArrays.jl](https://github.com/JuliaArrays/StaticArrays.jl)):
 
 ```@example type-stability
 typeof(reshape_operator_data((2, 2, 2)))
@@ -219,13 +219,15 @@ Finally, let's look at the benchmarks
 @benchmark reshape_operator_data($((2, 2, 2)))
 ```
 
-Which is an innocuous but huge difference in terms of performance. Hence, we highly recommend using `Tuple` or `SVector` when defining the dimensions of a user-defined [`QuantumObject`](@ref).
+Which is an innocuous but huge difference in terms of performance. Hence, we highly recommend using `Tuple` or `SVector` from [StaticArrays.jl](https://github.com/JuliaArrays/StaticArrays.jl) when defining the dimensions of a user-defined [`QuantumObject`](@ref).
 
 ## The use of `Val` in some `QuantumToolbox.jl` functions
 
 In some functions of `QuantumToolbox.jl`, you may find the use of the [`Val`](https://docs.julialang.org/en/v1/base/base/#Base.Val) type in the arguments. This is a trick to pass a value at compile time, and it is very useful to avoid type instabilities. Let's make a very simple example, where we want to create a Fock state ``|j\rangle`` of a given dimension `N`, and we give the possibility to create it as a sparse or dense vector. At first, we can write the function without using `Val`:
 
 ```@example type-stability
+using SparseArrays
+
 function my_fock(N::Int, j::Int = 0; sparse::Bool = false)
     if sparse
         array = sparsevec([j + 1], [1.0 + 0im], N)
diff --git a/docs/src/resources/api.md b/docs/src/resources/api.md
index eb768caca..fc1d9fc7c 100644
--- a/docs/src/resources/api.md
+++ b/docs/src/resources/api.md
@@ -1,5 +1,11 @@
 ```@meta
 CurrentModule = QuantumToolbox
+
+DocTestSetup = quote
+    using LinearAlgebra
+    using SparseArrays
+    using QuantumToolbox
+end
 ```
 
 # [API](@id doc-API)
diff --git a/docs/src/users_guide/QuantumObject/QuantumObject.md b/docs/src/users_guide/QuantumObject/QuantumObject.md
index 1a1aead59..0ac2a6fb7 100644
--- a/docs/src/users_guide/QuantumObject/QuantumObject.md
+++ b/docs/src/users_guide/QuantumObject/QuantumObject.md
@@ -45,6 +45,8 @@ Qobj(rand(4, 4))
 M = rand(ComplexF64, 4, 4)
 Qobj(M, dims = [2, 2])  # dims as Vector
 Qobj(M, dims = (2, 2))  # dims as Tuple (recommended)
+
+import QuantumToolbox: SVector # or using StaticArrays
 Qobj(M, dims = SVector(2, 2)) # dims as StaticArrays.SVector (recommended)
 ```
 
@@ -195,6 +197,8 @@ Vector{Int64}(v_d)
 ```
 
 ```@example Qobj
+using SparseArrays
+
 v_s = SparseVector(v_d)
 ```
 
diff --git a/docs/src/users_guide/QuantumObject/QuantumObject_functions.md b/docs/src/users_guide/QuantumObject/QuantumObject_functions.md
index ac2c62698..be76b0c06 100644
--- a/docs/src/users_guide/QuantumObject/QuantumObject_functions.md
+++ b/docs/src/users_guide/QuantumObject/QuantumObject_functions.md
@@ -43,8 +43,6 @@ Here is a table that summarizes all the supported linear algebra functions and a
 
 - [`eigenenergies`](@ref): return eigenenergies (eigenvalues)
 - [`eigenstates`](@ref): return [`EigsolveResult`](@ref) (contains eigenvalues and eigenvectors)
-- [`eigvals`](@ref): return eigenvalues
-- [`eigen`](@ref): using dense eigen solver and return [`EigsolveResult`](@ref) (contains eigenvalues and eigenvectors)
 - [`eigsolve`](@ref): using sparse eigen solver and return [`EigsolveResult`](@ref) (contains eigenvalues and eigenvectors)
 - [`eigsolve_al`](@ref): using the Arnoldi-Lindblad eigen solver and return [`EigsolveResult`](@ref) (contains eigenvalues and eigenvectors)
 
diff --git a/docs/src/users_guide/tensor.md b/docs/src/users_guide/tensor.md
index 61e125358..fc0e6cdd0 100644
--- a/docs/src/users_guide/tensor.md
+++ b/docs/src/users_guide/tensor.md
@@ -108,7 +108,7 @@ H = 0.5 * ωa * σz + ωc * a' * a + g * (a' * σm + a * σm')
 
 The partial trace is an operation that reduces the dimension of a Hilbert space by eliminating some degrees of freedom by averaging (tracing). In this sense it is therefore the converse of the tensor product. It is useful when one is interested in only a part of a coupled quantum system. For open quantum systems, this typically involves tracing over the environment leaving only the system of interest. In `QuantumToolbox` the function [`ptrace`](@ref) is used to take partial traces. [`ptrace`](@ref) takes one [`QuantumObject`](@ref) as an input, and also one argument `sel`, which marks the component systems that should be kept, and all the other components are traced out. 
 
-Remember that the index of `Julia` starts from `1`, and all the elements in `sel` should be positive `Integer`. Therefore, the type of `sel` can be either `Integer`, `Tuple`, `SVector`, or `Vector`.
+Remember that the index of `Julia` starts from `1`, and all the elements in `sel` should be positive `Integer`. Therefore, the type of `sel` can be either `Integer`, `Tuple`, `SVector` ([StaticArrays.jl](https://github.com/JuliaArrays/StaticArrays.jl)), or `Vector`.
 
 !!! warning "Beware of type-stability!"
     Although it supports also `Vector` type, it is recommended to use `Tuple` or `SVector` from [`StaticArrays.jl`](https://github.com/JuliaArrays/StaticArrays.jl) to improve performance. For a brief explanation on the impact of the type of `sel`, see the section [The Importance of Type-Stability](@ref doc:Type-Stability).
diff --git a/ext/QuantumToolboxCUDAExt.jl b/ext/QuantumToolboxCUDAExt.jl
index d05a65648..e676b4ea4 100644
--- a/ext/QuantumToolboxCUDAExt.jl
+++ b/ext/QuantumToolboxCUDAExt.jl
@@ -5,7 +5,7 @@ using QuantumToolbox: makeVal, getVal
 import QuantumToolbox: _sparse_similar, _convert_eltype_wordsize
 import CUDA: cu, CuArray, allowscalar
 import CUDA.CUSPARSE: CuSparseVector, CuSparseMatrixCSC, CuSparseMatrixCSR, AbstractCuSparseArray
-import SparseArrays: SparseVector, SparseMatrixCSC
+import SparseArrays: SparseVector, SparseMatrixCSC, sparse
 import CUDA.Adapt: adapt
 
 allowscalar(false)
diff --git a/src/QuantumToolbox.jl b/src/QuantumToolbox.jl
index 44b1aa84e..409bcfa54 100644
--- a/src/QuantumToolbox.jl
+++ b/src/QuantumToolbox.jl
@@ -1,19 +1,10 @@
 module QuantumToolbox
 
-# Re-export:
-#   1. StaticArraysCore.SVector for the type of dims
-#   2. basic functions in LinearAlgebra and SparseArrays
-#   3. some functions in SciMLOperators
-import Reexport: @reexport
-@reexport import StaticArraysCore: SVector
-@reexport using LinearAlgebra
-@reexport using SparseArrays
-@reexport import SciMLOperators: cache_operator, iscached, isconstant
-
-# other functions in LinearAlgebra
-import LinearAlgebra: BlasReal, BlasInt, BlasFloat, BlasComplex, checksquare
-import LinearAlgebra.BLAS: @blasfunc
+# Standard Julia libraries
+using LinearAlgebra
+import LinearAlgebra: BlasInt, BlasFloat, checksquare
 import LinearAlgebra.LAPACK: hseqr!
+using SparseArrays
 
 # SciML packages (for QobjEvo, OrdinaryDiffEq, and LinearSolve)
 import SciMLBase:
@@ -42,6 +33,9 @@ import SciMLBase:
     AbstractODESolution
 import StochasticDiffEq: StochasticDiffEqAlgorithm, SRA2, SRIW1
 import SciMLOperators:
+    cache_operator,
+    iscached,
+    isconstant,
     SciMLOperators,
     AbstractSciMLOperator,
     MatrixOperator,
@@ -67,11 +61,18 @@ import IncompleteLU: ilu
 import Pkg
 import Random: AbstractRNG, default_rng, seed!
 import SpecialFunctions: loggamma
-import StaticArraysCore: MVector
+import StaticArraysCore: SVector, MVector
+
+# Export functions from the other modules
+
+# LinearAlgebra
+export ishermitian, issymmetric, isposdef, dot, tr, svdvals, norm, normalize, normalize!, diag, Hermitian, Symmetric
+
+# SparseArrays
+export permute
 
-# Setting the number of threads to 1 allows
-# to achieve better performances for more massive parallelizations
-BLAS.set_num_threads(1)
+# SciMLOperators
+export cache_operator, iscached, isconstant
 
 # Utility
 include("utilities.jl")
diff --git a/src/qobj/arithmetic_and_attributes.jl b/src/qobj/arithmetic_and_attributes.jl
index 2cf9503cc..a18a383c8 100644
--- a/src/qobj/arithmetic_and_attributes.jl
+++ b/src/qobj/arithmetic_and_attributes.jl
@@ -36,16 +36,16 @@ end
 
 for op in (:(+), :(-), :(*))
     @eval begin
-        function LinearAlgebra.$op(A::AbstractQuantumObject, B::AbstractQuantumObject)
+        function Base.$op(A::AbstractQuantumObject, B::AbstractQuantumObject)
             check_dimensions(A, B)
             QType = promote_op_type(A, B)
             return QType($(op)(A.data, B.data), A.type, A.dimensions)
         end
-        LinearAlgebra.$op(A::AbstractQuantumObject) = get_typename_wrapper(A)($(op)(A.data), A.type, A.dimensions)
+        Base.$op(A::AbstractQuantumObject) = get_typename_wrapper(A)($(op)(A.data), A.type, A.dimensions)
 
-        LinearAlgebra.$op(n::T, A::AbstractQuantumObject) where {T<:Number} =
+        Base.$op(n::T, A::AbstractQuantumObject) where {T<:Number} =
             get_typename_wrapper(A)($(op)(n * I, A.data), A.type, A.dimensions)
-        LinearAlgebra.$op(A::AbstractQuantumObject, n::T) where {T<:Number} =
+        Base.$op(A::AbstractQuantumObject, n::T) where {T<:Number} =
             get_typename_wrapper(A)($(op)(A.data, n * I), A.type, A.dimensions)
     end
 end
@@ -63,7 +63,7 @@ for ADimType in (:Dimensions, :GeneralDimensions)
     for BDimType in (:Dimensions, :GeneralDimensions)
         if ADimType == BDimType == :Dimensions
             @eval begin
-                function LinearAlgebra.:(*)(
+                function Base.:(*)(
                     A::AbstractQuantumObject{OperatorQuantumObject,<:$ADimType},
                     B::AbstractQuantumObject{OperatorQuantumObject,<:$BDimType},
                 )
@@ -74,7 +74,7 @@ for ADimType in (:Dimensions, :GeneralDimensions)
             end
         else
             @eval begin
-                function LinearAlgebra.:(*)(
+                function Base.:(*)(
                     A::AbstractQuantumObject{OperatorQuantumObject,<:$ADimType},
                     B::AbstractQuantumObject{OperatorQuantumObject,<:$BDimType},
                 )
@@ -91,57 +91,41 @@ for ADimType in (:Dimensions, :GeneralDimensions)
     end
 end
 
-function LinearAlgebra.:(*)(
-    A::AbstractQuantumObject{OperatorQuantumObject},
-    B::QuantumObject{KetQuantumObject,<:Dimensions},
-)
+function Base.:(*)(A::AbstractQuantumObject{OperatorQuantumObject}, B::QuantumObject{KetQuantumObject,<:Dimensions})
     check_mul_dimensions(get_dimensions_from(A), get_dimensions_to(B))
     return QuantumObject(A.data * B.data, Ket, Dimensions(get_dimensions_to(A)))
 end
-function LinearAlgebra.:(*)(
-    A::QuantumObject{BraQuantumObject,<:Dimensions},
-    B::AbstractQuantumObject{OperatorQuantumObject},
-)
+function Base.:(*)(A::QuantumObject{BraQuantumObject,<:Dimensions}, B::AbstractQuantumObject{OperatorQuantumObject})
     check_mul_dimensions(get_dimensions_from(A), get_dimensions_to(B))
     return QuantumObject(A.data * B.data, Bra, Dimensions(get_dimensions_from(B)))
 end
-function LinearAlgebra.:(*)(A::QuantumObject{KetQuantumObject}, B::QuantumObject{BraQuantumObject})
+function Base.:(*)(A::QuantumObject{KetQuantumObject}, B::QuantumObject{BraQuantumObject})
     check_dimensions(A, B)
     return QuantumObject(A.data * B.data, Operator, A.dimensions) # to align with QuTiP, don't use kron(A, B) to do it.
 end
-function LinearAlgebra.:(*)(A::QuantumObject{BraQuantumObject}, B::QuantumObject{KetQuantumObject})
+function Base.:(*)(A::QuantumObject{BraQuantumObject}, B::QuantumObject{KetQuantumObject})
     check_dimensions(A, B)
     return A.data * B.data
 end
-function LinearAlgebra.:(*)(
-    A::AbstractQuantumObject{SuperOperatorQuantumObject},
-    B::QuantumObject{OperatorQuantumObject},
-)
+function Base.:(*)(A::AbstractQuantumObject{SuperOperatorQuantumObject}, B::QuantumObject{OperatorQuantumObject})
     check_dimensions(A, B)
     return QuantumObject(vec2mat(A.data * mat2vec(B.data)), Operator, A.dimensions)
 end
-function LinearAlgebra.:(*)(A::QuantumObject{OperatorBraQuantumObject}, B::QuantumObject{OperatorKetQuantumObject})
+function Base.:(*)(A::QuantumObject{OperatorBraQuantumObject}, B::QuantumObject{OperatorKetQuantumObject})
     check_dimensions(A, B)
     return A.data * B.data
 end
-function LinearAlgebra.:(*)(
-    A::AbstractQuantumObject{SuperOperatorQuantumObject},
-    B::QuantumObject{OperatorKetQuantumObject},
-)
+function Base.:(*)(A::AbstractQuantumObject{SuperOperatorQuantumObject}, B::QuantumObject{OperatorKetQuantumObject})
     check_dimensions(A, B)
     return QuantumObject(A.data * B.data, OperatorKet, A.dimensions)
 end
-function LinearAlgebra.:(*)(
-    A::QuantumObject{OperatorBraQuantumObject},
-    B::AbstractQuantumObject{SuperOperatorQuantumObject},
-)
+function Base.:(*)(A::QuantumObject{OperatorBraQuantumObject}, B::AbstractQuantumObject{SuperOperatorQuantumObject})
     check_dimensions(A, B)
     return QuantumObject(A.data * B.data, OperatorBra, A.dimensions)
 end
 
-LinearAlgebra.:(^)(A::QuantumObject, n::T) where {T<:Number} = QuantumObject(^(A.data, n), A.type, A.dimensions)
-LinearAlgebra.:(/)(A::AbstractQuantumObject, n::T) where {T<:Number} =
-    get_typename_wrapper(A)(A.data / n, A.type, A.dimensions)
+Base.:(^)(A::QuantumObject, n::T) where {T<:Number} = QuantumObject(^(A.data, n), A.type, A.dimensions)
+Base.:(/)(A::AbstractQuantumObject, n::T) where {T<:Number} = get_typename_wrapper(A)(A.data / n, A.type, A.dimensions)
 
 @doc raw"""
     A ⋅ B
@@ -221,7 +205,7 @@ Base.conj(A::AbstractQuantumObject) = get_typename_wrapper(A)(conj(A.data), A.ty
 
 Lazy matrix transpose of the [`AbstractQuantumObject`](@ref).
 """
-LinearAlgebra.transpose(
+Base.transpose(
     A::AbstractQuantumObject{OpType},
 ) where {OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} =
     get_typename_wrapper(A)(transpose(A.data), A.type, transpose(A.dimensions))
@@ -236,15 +220,13 @@ Lazy adjoint (conjugate transposition) of the [`AbstractQuantumObject`](@ref)
 !!! note
     `A'` and `dag(A)` are synonyms of `adjoint(A)`.
 """
-LinearAlgebra.adjoint(
-    A::AbstractQuantumObject{OpType},
-) where {OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} =
+Base.adjoint(A::AbstractQuantumObject{OpType}) where {OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} =
     get_typename_wrapper(A)(adjoint(A.data), A.type, adjoint(A.dimensions))
-LinearAlgebra.adjoint(A::QuantumObject{KetQuantumObject}) = QuantumObject(adjoint(A.data), Bra, adjoint(A.dimensions))
-LinearAlgebra.adjoint(A::QuantumObject{BraQuantumObject}) = QuantumObject(adjoint(A.data), Ket, adjoint(A.dimensions))
-LinearAlgebra.adjoint(A::QuantumObject{OperatorKetQuantumObject}) =
+Base.adjoint(A::QuantumObject{KetQuantumObject}) = QuantumObject(adjoint(A.data), Bra, adjoint(A.dimensions))
+Base.adjoint(A::QuantumObject{BraQuantumObject}) = QuantumObject(adjoint(A.data), Ket, adjoint(A.dimensions))
+Base.adjoint(A::QuantumObject{OperatorKetQuantumObject}) =
     QuantumObject(adjoint(A.data), OperatorBra, adjoint(A.dimensions))
-LinearAlgebra.adjoint(A::QuantumObject{OperatorBraQuantumObject}) =
+Base.adjoint(A::QuantumObject{OperatorBraQuantumObject}) =
     QuantumObject(adjoint(A.data), OperatorKet, adjoint(A.dimensions))
 
 @doc raw"""
@@ -415,7 +397,7 @@ Matrix square root of [`QuantumObject`](@ref)
 !!! note
     `√(A)` (where `√` can be typed by tab-completing `\sqrt` in the REPL) is a synonym of `sqrt(A)`.
 """
-LinearAlgebra.sqrt(A::QuantumObject) = QuantumObject(sqrt(to_dense(A.data)), A.type, A.dimensions)
+Base.sqrt(A::QuantumObject) = QuantumObject(sqrt(to_dense(A.data)), A.type, A.dimensions)
 
 @doc raw"""
     log(A::QuantumObject)
@@ -424,7 +406,7 @@ Matrix logarithm of [`QuantumObject`](@ref)
 
 Note that this function only supports for [`Operator`](@ref) and [`SuperOperator`](@ref)
 """
-LinearAlgebra.log(A::QuantumObject{ObjType}) where {ObjType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} =
+Base.log(A::QuantumObject{ObjType}) where {ObjType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} =
     QuantumObject(log(to_dense(A.data)), A.type, A.dimensions)
 
 @doc raw"""
@@ -434,11 +416,11 @@ Matrix exponential of [`QuantumObject`](@ref)
 
 Note that this function only supports for [`Operator`](@ref) and [`SuperOperator`](@ref)
 """
-LinearAlgebra.exp(
+Base.exp(
     A::QuantumObject{ObjType,DimsType,<:AbstractMatrix},
 ) where {ObjType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject},DimsType} =
     QuantumObject(to_sparse(exp(A.data)), A.type, A.dimensions)
-LinearAlgebra.exp(
+Base.exp(
     A::QuantumObject{ObjType,DimsType,<:AbstractSparseMatrix},
 ) where {ObjType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject},DimsType} =
     QuantumObject(_spexp(A.data), A.type, A.dimensions)
@@ -483,7 +465,7 @@ Matrix sine of [`QuantumObject`](@ref), defined as
 
 Note that this function only supports for [`Operator`](@ref) and [`SuperOperator`](@ref)
 """
-LinearAlgebra.sin(A::QuantumObject{ObjType}) where {ObjType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} =
+Base.sin(A::QuantumObject{ObjType}) where {ObjType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} =
     (exp(1im * A) - exp(-1im * A)) / 2im
 
 @doc raw"""
@@ -495,7 +477,7 @@ Matrix cosine of [`QuantumObject`](@ref), defined as
 
 Note that this function only supports for [`Operator`](@ref) and [`SuperOperator`](@ref)
 """
-LinearAlgebra.cos(A::QuantumObject{ObjType}) where {ObjType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} =
+Base.cos(A::QuantumObject{ObjType}) where {ObjType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} =
     (exp(1im * A) + exp(-1im * A)) / 2
 
 @doc raw"""
@@ -768,7 +750,7 @@ true
 ```
 
 !!! warning "Beware of type-stability!"
-    It is highly recommended to use `permute(A, order)` with `order` as `Tuple` or `SVector` to keep type stability. See the [related Section](@ref doc:Type-Stability) about type stability for more details.
+    It is highly recommended to use `permute(A, order)` with `order` as `Tuple` or `SVector` from [StaticArrays.jl](https://github.com/JuliaArrays/StaticArrays.jl) to keep type stability. See the [related Section](@ref doc:Type-Stability) about type stability for more details.
 """
 function SparseArrays.permute(
     A::QuantumObject{ObjType},
diff --git a/src/qobj/dimensions.jl b/src/qobj/dimensions.jl
index 0fd19f3a9..a10d73c64 100644
--- a/src/qobj/dimensions.jl
+++ b/src/qobj/dimensions.jl
@@ -91,9 +91,9 @@ _get_space_size(s::AbstractSpace)::Int = s.size
 Base.prod(dims::Dimensions) = prod(dims.to)
 Base.prod(spaces::NTuple{N,AbstractSpace}) where {N} = prod(_get_space_size, spaces)
 
-LinearAlgebra.transpose(dimensions::Dimensions) = dimensions
-LinearAlgebra.transpose(dimensions::GeneralDimensions) = GeneralDimensions(dimensions.from, dimensions.to) # switch `to` and `from`
-LinearAlgebra.adjoint(dimensions::AbstractDimensions) = transpose(dimensions)
+Base.transpose(dimensions::Dimensions) = dimensions
+Base.transpose(dimensions::GeneralDimensions) = GeneralDimensions(dimensions.from, dimensions.to) # switch `to` and `from`
+Base.adjoint(dimensions::AbstractDimensions) = transpose(dimensions)
 
 # this is used to show `dims` for Qobj and QobjEvo
 _get_dims_string(dimensions::Dimensions) = string(dimensions_to_dims(dimensions))
diff --git a/src/qobj/eigsolve.jl b/src/qobj/eigsolve.jl
index d68a3a163..986a25321 100644
--- a/src/qobj/eigsolve.jl
+++ b/src/qobj/eigsolve.jl
@@ -444,6 +444,8 @@ julia> a = destroy(5);
 
 julia> H = a + a';
 
+julia> using LinearAlgebra;
+
 julia> E, ψ, U = eigen(H)
 EigsolveResult:   type=Operator   dims=[5]
 values:
diff --git a/src/qobj/functions.jl b/src/qobj/functions.jl
index 4cb8d6c38..694e6b05f 100644
--- a/src/qobj/functions.jl
+++ b/src/qobj/functions.jl
@@ -187,7 +187,7 @@ julia> a.dims, O.dims
 ([20], [20, 20])
 ```
 """
-function LinearAlgebra.kron(
+function Base.kron(
     A::AbstractQuantumObject{OpType,<:Dimensions},
     B::AbstractQuantumObject{OpType,<:Dimensions},
 ) where {OpType<:Union{KetQuantumObject,BraQuantumObject,OperatorQuantumObject}}
@@ -201,7 +201,7 @@ for ADimType in (:Dimensions, :GeneralDimensions)
     for BDimType in (:Dimensions, :GeneralDimensions)
         if !(ADimType == BDimType == :Dimensions) # not for this case because it's already implemented
             @eval begin
-                function LinearAlgebra.kron(
+                function Base.kron(
                     A::AbstractQuantumObject{OperatorQuantumObject,<:$ADimType},
                     B::AbstractQuantumObject{OperatorQuantumObject,<:$BDimType},
                 )
@@ -226,7 +226,7 @@ for AOpType in (:KetQuantumObject, :BraQuantumObject, :OperatorQuantumObject)
     for BOpType in (:KetQuantumObject, :BraQuantumObject, :OperatorQuantumObject)
         if (AOpType != BOpType)
             @eval begin
-                function LinearAlgebra.kron(A::AbstractQuantumObject{$AOpType}, B::AbstractQuantumObject{$BOpType})
+                function Base.kron(A::AbstractQuantumObject{$AOpType}, B::AbstractQuantumObject{$BOpType})
                     QType = promote_op_type(A, B)
                     _lazy_tensor_warning(A.data, B.data)
                     return QType(
@@ -243,8 +243,8 @@ for AOpType in (:KetQuantumObject, :BraQuantumObject, :OperatorQuantumObject)
     end
 end
 
-LinearAlgebra.kron(A::AbstractQuantumObject) = A
-function LinearAlgebra.kron(A::Vector{<:AbstractQuantumObject})
+Base.kron(A::AbstractQuantumObject) = A
+function Base.kron(A::Vector{<:AbstractQuantumObject})
     @warn "`tensor(A)` or `kron(A)` with `A` is a `Vector` can hurt performance. Try to use `tensor(A...)` or `kron(A...)` instead."
     return kron(A...)
 end
diff --git a/src/qobj/operators.jl b/src/qobj/operators.jl
index 90d6d1739..d7c460fdb 100644
--- a/src/qobj/operators.jl
+++ b/src/qobj/operators.jl
@@ -29,7 +29,7 @@ The `distribution` specifies which of the method used to obtain the unitary matr
 1. [F. Mezzadri, How to generate random matrices from the classical compact groups, arXiv:math-ph/0609050 (2007)](https://arxiv.org/abs/math-ph/0609050)
 
 !!! warning "Beware of type-stability!"
-    If you want to keep type stability, it is recommended to use `rand_unitary(dimensions, Val(distribution))` instead of `rand_unitary(dimensions, distribution)`. Also, put `dimensions` as `Tuple` or `SVector`. See [this link](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-value-type) and the [related Section](@ref doc:Type-Stability) about type stability for more details.
+    If you want to keep type stability, it is recommended to use `rand_unitary(dimensions, Val(distribution))` instead of `rand_unitary(dimensions, distribution)`. Also, put `dimensions` as `Tuple` or `SVector` from [StaticArrays.jl](https://github.com/JuliaArrays/StaticArrays.jl). See [this link](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-value-type) and the [related Section](@ref doc:Type-Stability) about type stability for more details.
 """
 rand_unitary(dimensions::Int, distribution::Union{Symbol,Val} = Val(:haar)) =
     rand_unitary(SVector(dimensions), makeVal(distribution))
@@ -553,7 +553,7 @@ The `dimensions` can be either the following types:
 where ``\omega = \exp(\frac{2 \pi i}{N})``.
 
 !!! warning "Beware of type-stability!"
-    It is highly recommended to use `qft(dimensions)` with `dimensions` as `Tuple` or `SVector` to keep type stability. See the [related Section](@ref doc:Type-Stability) about type stability for more details.
+    It is highly recommended to use `qft(dimensions)` with `dimensions` as `Tuple` or `SVector` from [StaticArrays.jl](https://github.com/JuliaArrays/StaticArrays.jl) to keep type stability. See the [related Section](@ref doc:Type-Stability) about type stability for more details.
 """
 qft(dimensions::Int) = QuantumObject(_qft_op(dimensions), Operator, dimensions)
 qft(dimensions::Union{Dimensions,AbstractVector{Int},Tuple}) =
diff --git a/src/qobj/quantum_object.jl b/src/qobj/quantum_object.jl
index 9f760c876..763113776 100644
--- a/src/qobj/quantum_object.jl
+++ b/src/qobj/quantum_object.jl
@@ -37,7 +37,7 @@ julia> a isa QuantumObject
 true
 
 julia> a.dims
-1-element SVector{1, Int64} with indices SOneTo(1):
+1-element StaticArraysCore.SVector{1, Int64} with indices SOneTo(1):
  20
 
 julia> a.dimensions
diff --git a/src/qobj/states.jl b/src/qobj/states.jl
index c2759c9a4..4474ebd16 100644
--- a/src/qobj/states.jl
+++ b/src/qobj/states.jl
@@ -17,7 +17,7 @@ The `dimensions` can be either the following types:
 - `dimensions::Union{Dimensions,AbstractVector{Int}, Tuple}`: list of dimensions representing the each number of basis in the subsystems.
 
 !!! warning "Beware of type-stability!"
-    It is highly recommended to use `zero_ket(dimensions)` with `dimensions` as `Tuple` or `SVector` to keep type stability. See the [related Section](@ref doc:Type-Stability) about type stability for more details.
+    It is highly recommended to use `zero_ket(dimensions)` with `dimensions` as `Tuple` or `SVector` from [StaticArrays.jl](https://github.com/JuliaArrays/StaticArrays.jl) to keep type stability. See the [related Section](@ref doc:Type-Stability) about type stability for more details.
 """
 zero_ket(dimensions::Int) = QuantumObject(zeros(ComplexF64, dimensions), Ket, dimensions)
 zero_ket(dimensions::Union{Dimensions,AbstractVector{Int},Tuple}) =
@@ -31,7 +31,7 @@ Generates a fock state ``\ket{\psi}`` of dimension `N`.
 It is also possible to specify the list of dimensions `dims` if different subsystems are present.
 
 !!! warning "Beware of type-stability!"
-    If you want to keep type stability, it is recommended to use `fock(N, j, dims=dims, sparse=Val(sparse))` instead of `fock(N, j, dims=dims, sparse=sparse)`. Consider also to use `dims` as a `Tuple` or `SVector` instead of `Vector`. See [this link](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-value-type) and the [related Section](@ref doc:Type-Stability) about type stability for more details.
+    If you want to keep type stability, it is recommended to use `fock(N, j, dims=dims, sparse=Val(sparse))` instead of `fock(N, j, dims=dims, sparse=sparse)`. Consider also to use `dims` as a `Tuple` or `SVector` from [StaticArrays.jl](https://github.com/JuliaArrays/StaticArrays.jl) instead of `Vector`. See [this link](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-value-type) and the [related Section](@ref doc:Type-Stability) about type stability for more details.
 """
 function fock(N::Int, j::Int = 0; dims::Union{Int,AbstractVector{Int},Tuple} = N, sparse::Union{Bool,Val} = Val(false))
     if getVal(sparse)
@@ -50,7 +50,7 @@ Generates a fock state like [`fock`](@ref).
 It is also possible to specify the list of dimensions `dims` if different subsystems are present.
 
 !!! warning "Beware of type-stability!"
-    If you want to keep type stability, it is recommended to use `basis(N, j, dims=dims)` with `dims` as a `Tuple` or `SVector` instead of `Vector`. See [this link](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-value-type) and the [related Section](@ref doc:Type-Stability) about type stability for more details.
+    If you want to keep type stability, it is recommended to use `basis(N, j, dims=dims)` with `dims` as a `Tuple` or `SVector` from [StaticArrays.jl](https://github.com/JuliaArrays/StaticArrays.jl) instead of `Vector`. See [this link](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-value-type) and the [related Section](@ref doc:Type-Stability) about type stability for more details.
 """
 basis(N::Int, j::Int = 0; dims::Union{Int,AbstractVector{Int},Tuple} = N) = fock(N, j, dims = dims)
 
@@ -73,7 +73,7 @@ The `dimensions` can be either the following types:
 - `dimensions::Union{Dimensions,AbstractVector{Int},Tuple}`: list of dimensions representing the each number of basis in the subsystems.
 
 !!! warning "Beware of type-stability!"
-    If you want to keep type stability, it is recommended to use `rand_ket(dimensions)` with `dimensions` as `Tuple` or `SVector` to keep type stability. See the [related Section](@ref doc:Type-Stability) about type stability for more details.
+    If you want to keep type stability, it is recommended to use `rand_ket(dimensions)` with `dimensions` as `Tuple` or `SVector` from [StaticArrays.jl](https://github.com/JuliaArrays/StaticArrays.jl) to keep type stability. See the [related Section](@ref doc:Type-Stability) about type stability for more details.
 """
 rand_ket(dimensions::Int) = rand_ket(SVector(dimensions))
 function rand_ket(dimensions::Union{Dimensions,AbstractVector{Int},Tuple})
@@ -90,7 +90,7 @@ Density matrix representation of a Fock state.
 Constructed via outer product of [`fock`](@ref).
 
 !!! warning "Beware of type-stability!"
-    If you want to keep type stability, it is recommended to use `fock_dm(N, j, dims=dims, sparse=Val(sparse))` instead of `fock_dm(N, j, dims=dims, sparse=sparse)`. Consider also to use `dims` as a `Tuple` or `SVector` instead of `Vector`. See [this link](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-value-type) and the [related Section](@ref doc:Type-Stability) about type stability for more details.
+    If you want to keep type stability, it is recommended to use `fock_dm(N, j, dims=dims, sparse=Val(sparse))` instead of `fock_dm(N, j, dims=dims, sparse=sparse)`. Consider also to use `dims` as a `Tuple` or `SVector` from [StaticArrays.jl](https://github.com/JuliaArrays/StaticArrays.jl) instead of `Vector`. See [this link](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-value-type) and the [related Section](@ref doc:Type-Stability) about type stability for more details.
 """
 function fock_dm(
     N::Int,
@@ -146,7 +146,7 @@ The `dimensions` can be either the following types:
 - `dimensions::Union{Dimensions,AbstractVector{Int},Tuple}`: list of dimensions representing the each number of basis in the subsystems.
 
 !!! warning "Beware of type-stability!"
-    If you want to keep type stability, it is recommended to use `maximally_mixed_dm(dimensions)` with `dimensions` as `Tuple` or `SVector` to keep type stability. See the [related Section](@ref doc:Type-Stability) about type stability for more details.
+    If you want to keep type stability, it is recommended to use `maximally_mixed_dm(dimensions)` with `dimensions` as `Tuple` or `SVector` from [StaticArrays.jl](https://github.com/JuliaArrays/StaticArrays.jl) to keep type stability. See the [related Section](@ref doc:Type-Stability) about type stability for more details.
 """
 maximally_mixed_dm(dimensions::Int) = QuantumObject(I(dimensions) / complex(dimensions), Operator, SVector(dimensions))
 function maximally_mixed_dm(dimensions::Union{Dimensions,AbstractVector{Int},Tuple})
@@ -166,7 +166,7 @@ The `dimensions` can be either the following types:
 The default keyword argument `rank = prod(dimensions)` (full rank).
 
 !!! warning "Beware of type-stability!"
-    If you want to keep type stability, it is recommended to use `rand_dm(dimensions; rank=rank)` with `dimensions` as `Tuple` or `SVector` instead of `Vector`. See [this link](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-value-type) and the [related Section](@ref doc:Type-Stability) about type stability for more details.
+    If you want to keep type stability, it is recommended to use `rand_dm(dimensions; rank=rank)` with `dimensions` as `Tuple` or `SVector` from [StaticArrays.jl](https://github.com/JuliaArrays/StaticArrays.jl) instead of `Vector`. See [this link](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-value-type) and the [related Section](@ref doc:Type-Stability) about type stability for more details.
 
 # References
 - [J. Ginibre, Statistical ensembles of complex, quaternion, and real matrices, Journal of Mathematical Physics 6.3 (1965): 440-449](https://doi.org/10.1063/1.1704292)
diff --git a/src/spin_lattice.jl b/src/spin_lattice.jl
index 72ee4f38f..1abfff235 100644
--- a/src/spin_lattice.jl
+++ b/src/spin_lattice.jl
@@ -31,7 +31,7 @@ A Julia function for generating a multi-site operator ``\\hat{O} = \\hat{O}_i \\
 julia> op = multisite_operator(Val(8), 5=>sigmax(), 7=>sigmaz());
 
 julia> op.dims
-8-element SVector{8, Int64} with indices SOneTo(8):
+8-element StaticArraysCore.SVector{8, Int64} with indices SOneTo(8):
  2
  2
  2
diff --git a/src/utilities.jl b/src/utilities.jl
index a3547cf61..ebe9481f4 100644
--- a/src/utilities.jl
+++ b/src/utilities.jl
@@ -151,9 +151,11 @@ _non_static_array_warning(argname, arg::Tuple{}) =
     throw(ArgumentError("The argument $argname must be a Tuple or a StaticVector of non-zero length."))
 _non_static_array_warning(argname, arg::Union{SVector{N,T},MVector{N,T},NTuple{N,T}}) where {N,T} = nothing
 _non_static_array_warning(argname, arg::AbstractVector{T}) where {T} =
-    @warn "The argument $argname should be a Tuple or a StaticVector for better performance. Try to use `$argname = $(Tuple(arg))` or `$argname = SVector(" *
+    @warn "The argument $argname should be a Tuple or a StaticVector for better performance. Try to use `$argname = $(Tuple(arg))` instead of `$argname = $arg`. " *
+          "Alternatively, you can do `import QuantumToolbox: SVector` " *
+          "and use `$argname = SVector(" *
           join(arg, ", ") *
-          ")` instead of `$argname = $arg`." maxlog = 1
+          ")`." maxlog = 1
 
 # lazy tensor warning
 for AType in (:AbstractArray, :AbstractSciMLOperator)
diff --git a/test/core-test/eigenvalues_and_operators.jl b/test/core-test/eigenvalues_and_operators.jl
index 3faebf0dc..0b4959beb 100644
--- a/test/core-test/eigenvalues_and_operators.jl
+++ b/test/core-test/eigenvalues_and_operators.jl
@@ -58,7 +58,7 @@
 
     # eigen solve for general matrices
     vals, _, vecs = eigsolve(L.data, sigma = 0.01, eigvals = 10, krylovdim = 50)
-    vals2, vecs2 = eigen(to_dense(L.data))
+    vals2, _, vecs2 = eigenstates(L)
     vals3, state3, vecs3 = eigsolve_al(L, 1 \ (40 * κ), eigvals = 10, krylovdim = 50)
     idxs = sortperm(vals2, by = abs)
     vals2 = vals2[idxs][1:10]
diff --git a/test/core-test/quantum_objects.jl b/test/core-test/quantum_objects.jl
index 8aebcc244..52adfb7fb 100644
--- a/test/core-test/quantum_objects.jl
+++ b/test/core-test/quantum_objects.jl
@@ -44,7 +44,7 @@
         @test_throws DomainError Qobj(rand(2, 2), dims = (2, -2))
         @test_logs (
             :warn,
-            "The argument dims should be a Tuple or a StaticVector for better performance. Try to use `dims = (2, 2)` or `dims = SVector(2, 2)` instead of `dims = [2, 2]`.",
+            "The argument dims should be a Tuple or a StaticVector for better performance. Try to use `dims = (2, 2)` instead of `dims = [2, 2]`. Alternatively, you can do `import QuantumToolbox: SVector` and use `dims = SVector(2, 2)`.",
         ) Qobj(rand(4, 4), dims = [2, 2])
     end
 
@@ -178,7 +178,7 @@
         a = sprand(ComplexF64, 100, 100, 0.1)
         a2 = Qobj(a)
         a3 = Qobj(a, type = SuperOperator)
-        a4 = sparse(a2)
+        a4 = to_sparse(a2)
         @test isequal(a4, a2) == true
         @test isequal(a4, a3) == false
         @test a4 ≈ a2
@@ -412,7 +412,7 @@
     @testset "projection" begin
         N = 10
         ψ = fock(N, 3)
-        @test proj(ψ) ≈ proj(ψ') ≈ sparse(ket2dm(ψ)) ≈ projection(N, 3, 3)
+        @test proj(ψ) ≈ proj(ψ') ≈ to_sparse(ket2dm(ψ)) ≈ projection(N, 3, 3)
         @test isket(ψ') == false
         @test isbra(ψ') == true
         @test shape(ψ) == (N,)
@@ -715,11 +715,11 @@
         @test ptrace(ρtotal, (1, 3, 4)) ≈ ptrace(ρtotal, (3, 1, 4)) # check sort of sel
         @test_logs (
             :warn,
-            "The argument sel should be a Tuple or a StaticVector for better performance. Try to use `sel = (1, 2)` or `sel = SVector(1, 2)` instead of `sel = [1, 2]`.",
+            "The argument sel should be a Tuple or a StaticVector for better performance. Try to use `sel = (1, 2)` instead of `sel = [1, 2]`. Alternatively, you can do `import QuantumToolbox: SVector` and use `sel = SVector(1, 2)`.",
         ) ptrace(ψtotal, [1, 2])
         @test_logs (
             :warn,
-            "The argument sel should be a Tuple or a StaticVector for better performance. Try to use `sel = (1, 2)` or `sel = SVector(1, 2)` instead of `sel = [1, 2]`.",
+            "The argument sel should be a Tuple or a StaticVector for better performance. Try to use `sel = (1, 2)` instead of `sel = [1, 2]`. Alternatively, you can do `import QuantumToolbox: SVector` and use `sel = SVector(1, 2)`.",
         ) ptrace(ρtotal, [1, 2])
         @test_throws ArgumentError ptrace(ψtotal, 0)
         @test_throws ArgumentError ptrace(ψtotal, 5)
diff --git a/test/core-test/quantum_objects_evo.jl b/test/core-test/quantum_objects_evo.jl
index fec8235c0..5367599fe 100644
--- a/test/core-test/quantum_objects_evo.jl
+++ b/test/core-test/quantum_objects_evo.jl
@@ -31,7 +31,7 @@
         @test_throws DomainError QobjEvo(a, dims = (2, -2))
         @test_logs (
             :warn,
-            "The argument dims should be a Tuple or a StaticVector for better performance. Try to use `dims = (2, 2)` or `dims = SVector(2, 2)` instead of `dims = [2, 2]`.",
+            "The argument dims should be a Tuple or a StaticVector for better performance. Try to use `dims = (2, 2)` instead of `dims = [2, 2]`. Alternatively, you can do `import QuantumToolbox: SVector` and use `dims = SVector(2, 2)`.",
         ) QobjEvo(MatrixOperator(rand(4, 4)), dims = [2, 2])
     end
 
diff --git a/test/core-test/states_and_operators.jl b/test/core-test/states_and_operators.jl
index f43b3e2a3..0cdfbd9e5 100644
--- a/test/core-test/states_and_operators.jl
+++ b/test/core-test/states_and_operators.jl
@@ -283,7 +283,7 @@
         sites = 4
         SIZE = 2^sites
         dims = ntuple(i -> 2, Val(sites))
-        Q_iden = Qobj(sparse((1.0 + 0.0im) * LinearAlgebra.I, SIZE, SIZE); dims = dims)
+        Q_iden = Qobj(sparse((1.0 + 0.0im) * I, SIZE, SIZE); dims = dims)
         Q_zero = Qobj(spzeros(ComplexF64, SIZE, SIZE); dims = dims)
         for i in 1:sites
             d_i = fdestroy(sites, i)
diff --git a/test/runtests.jl b/test/runtests.jl
index c8a8efb4e..d8670d9d3 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,6 +1,11 @@
 using Test
 using Pkg
 
+# Importing only the necessary functions to keep track the re-export of the functions
+import LinearAlgebra: Diagonal, I, mul!, triu, tril, triu!, tril!
+import SparseArrays: sparse, sprand, spzeros, spdiagm, nnz, SparseVector, SparseMatrixCSC, AbstractSparseMatrix
+import StaticArraysCore: SVector
+
 const GROUP = get(ENV, "GROUP", "All")
 
 const testdir = dirname(@__FILE__)