diff --git a/.github/workflows/dependencies/dependencies_nvcc.sh b/.github/workflows/dependencies/dependencies_nvcc.sh
index 2578bd33fe7..14bae699d7e 100755
--- a/.github/workflows/dependencies/dependencies_nvcc.sh
+++ b/.github/workflows/dependencies/dependencies_nvcc.sh
@@ -36,5 +36,6 @@ sudo apt-get install -y \
     cuda-nvml-dev-$VERSION_DASHED           \
     cuda-nvtx-$VERSION_DASHED               \
     libcufft-dev-$VERSION_DASHED            \
-    libcurand-dev-$VERSION_DASHED
+    libcurand-dev-$VERSION_DASHED           \
+    libcusparse-dev-$VERSION_DASHED
 sudo ln -s cuda-$VERSION_DOTTED /usr/local/cuda
diff --git a/Src/Base/AMReX_TableData.H b/Src/Base/AMReX_TableData.H
index ee2471d36cb..02907fb089a 100644
--- a/Src/Base/AMReX_TableData.H
+++ b/Src/Base/AMReX_TableData.H
@@ -15,12 +15,12 @@
 
 namespace amrex {
 
-template <typename T>
+template <typename T, typename IDX = int>
 struct Table1D
 {
     T* AMREX_RESTRICT p = nullptr;
-    int begin = 1;
-    int end = 0;
+    IDX begin = 1;
+    IDX end = 0;
 
     constexpr Table1D () noexcept = default;
 
@@ -33,7 +33,7 @@ struct Table1D
         {}
 
     AMREX_GPU_HOST_DEVICE
-    constexpr Table1D (T* a_p, int a_begin, int a_end) noexcept
+    constexpr Table1D (T* a_p, IDX a_begin, IDX a_end) noexcept
         : p(a_p),
           begin(a_begin),
           end(a_end)
@@ -44,7 +44,7 @@ struct Table1D
 
     template <class U=T, std::enable_if_t<!std::is_void_v<U>,int> = 0>
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-    U& operator() (int i) const noexcept {
+    U& operator() (IDX i) const noexcept {
 #if defined(AMREX_DEBUG) || defined(AMREX_BOUND_CHECK)
         index_assert(i);
 #endif
@@ -53,14 +53,30 @@ struct Table1D
 
 #if defined(AMREX_DEBUG) || defined(AMREX_BOUND_CHECK)
     AMREX_GPU_HOST_DEVICE inline
-    void index_assert (int i) const
+    void index_assert (IDX i) const
     {
         if (i < begin || i >= end) {
-            AMREX_IF_ON_DEVICE((
-                AMREX_DEVICE_PRINTF(" (%d) is out of bound (%d:%d)\n",
-                                    i, begin, end-1);
-                amrex::Abort();
-            ))
+            if constexpr (std::is_same_v<IDX,int>) {
+                AMREX_IF_ON_DEVICE((
+                    AMREX_DEVICE_PRINTF(" (%d) is out of bound (%d:%d)\n",
+                                        i, begin, end-1);
+                    amrex::Abort();
+                ))
+            } else if constexpr (std::is_same_v<IDX,long>) {
+                AMREX_IF_ON_DEVICE((
+                    AMREX_DEVICE_PRINTF(" (%ld) is out of bound (%ld:%ld)\n",
+                                        i, begin, end-1);
+                    amrex::Abort();
+                ))
+            } else if constexpr (std::is_same_v<IDX,long long>) {
+                AMREX_IF_ON_DEVICE((
+                    AMREX_DEVICE_PRINTF(" (%lld) is out of bound (%lld:%lld)\n",
+                                        i, begin, end-1);
+                    amrex::Abort();
+                ))
+            } else {
+                AMREX_IF_ON_DEVICE(( amrex::Abort(); ))
+            }
             AMREX_IF_ON_HOST((
                 std::stringstream ss;
                 ss << " (" << i << ") is out of bound ("
diff --git a/Src/LinearSolvers/AMReX_AlgVector.H b/Src/LinearSolvers/AMReX_AlgVector.H
index 92d5fe091d6..07c6fbffeaa 100644
--- a/Src/LinearSolvers/AMReX_AlgVector.H
+++ b/Src/LinearSolvers/AMReX_AlgVector.H
@@ -6,6 +6,7 @@
 #include <AMReX_FabArray.H>
 #include <AMReX_INT.H>
 #include <AMReX_LayoutData.H>
+#include <AMReX_TableData.H>
 
 #include <fstream>
 #include <string>
@@ -13,51 +14,6 @@
 
 namespace amrex {
 
-template <typename T>
-struct VectorView
-{
-    T* AMREX_RESTRICT p = nullptr;
-    Long begin = 1;
-    Long end = 0;
-
-    AMREX_GPU_HOST_DEVICE
-    explicit operator bool() const noexcept { return p != nullptr; }
-
-    template <class U=T, std::enable_if_t<!std::is_void_v<U>,int> = 0>
-    [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-    U& operator[] (Long i) const noexcept {
-#if defined(AMREX_DEBUG) || defined(AMREX_BOUND_CHECK)
-        index_assert(i);
-#endif
-        return p[i-begin];
-    }
-
-    [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-    bool contains (Long i) const noexcept { return i>=begin && i<end; }
-
-#if defined(AMREX_DEBUG) || defined(AMREX_BOUND_CHECK)
-    AMREX_GPU_HOST_DEVICE inline
-    void index_assert (Long i) const
-    {
-        if ( ! contains(i)) {
-            AMREX_IF_ON_DEVICE((
-                // xxxxx TODO: Yes, this is wrong on Windows.
-                AMREX_DEVICE_PRINTF(" [%ld] is out of bound [%ld:%ld]\n",
-                                    i, begin, end-1);
-                amrex::Abort();
-            ))
-            AMREX_IF_ON_HOST((
-                std::stringstream ss;
-                    ss << " [" << i << "] is out of bound ["
-                    << begin << ":" << end-1 << "]";
-                amrex::Abort(ss.str());
-            ))
-        }
-    }
-#endif
-
-};
-
 template <typename T, typename Allocator = DefaultAllocator<T> >
 class AlgVector
 {
@@ -94,18 +50,18 @@ public:
     [[nodiscard]] T      * data ()       { return m_data.data(); }
 
     [[nodiscard]] AMREX_FORCE_INLINE
-    VectorView<T const> view () const {
-        return VectorView<T const>{m_data.data(), m_begin, m_end};
+    Table1D<T const, Long> view () const {
+        return Table1D<T const, Long>{m_data.data(), m_begin, m_end};
     }
 
     [[nodiscard]] AMREX_FORCE_INLINE
-    VectorView<T const> const_view () const {
-        return VectorView<T const>{m_data.data(), m_begin, m_end};
+    Table1D<T const, Long> const_view () const {
+        return Table1D<T const, Long>{m_data.data(), m_begin, m_end};
     }
 
     [[nodiscard]] AMREX_FORCE_INLINE
-    VectorView<T> view () {
-        return VectorView<T>{m_data.data(), m_begin, m_end};
+    Table1D<T, Long> view () {
+        return Table1D<T, Long>{m_data.data(), m_begin, m_end};
     }
 
     void setVal (T val);
diff --git a/Src/LinearSolvers/AMReX_SpMV.H b/Src/LinearSolvers/AMReX_SpMV.H
index 11d0f3412f5..4a13c2f8030 100644
--- a/Src/LinearSolvers/AMReX_SpMV.H
+++ b/Src/LinearSolvers/AMReX_SpMV.H
@@ -3,13 +3,11 @@
 #include <AMReX_Config.H>
 
 #include <AMReX_AlgVector.H>
+#include <AMReX_GpuComplex.H>
 #include <AMReX_SpMatrix.H>
 
 #if defined(AMREX_USE_CUDA)
 #  include <cusparse.h>
-//#  if defined(__CUDACC__) && (__CUDACC_VER_MAJOR__ >= 11)
-//#    include <cub/cub.cuh>
-//#  endif
 #elif defined(AMREX_USE_HIP)
 #  include <hipsparse.h>
 #elif defined(AMREX_USE_DPCPP)
@@ -36,34 +34,29 @@ void SpMV (AlgVector<T>& y, SpMatrix<T> const& A, AlgVector<T> const& x)
 #if defined(AMREX_USE_GPU)
 
     Long const nrows = A.numLocalRows();
-    Long const ncols = y.numLocalRows();
+    Long const ncols = x.numLocalRows();
     Long const nnz = A.numLocalNonZero();
 
-    // y.setVal(0);
-
 #if defined(AMREX_USE_CUDA)
 
-#if 0
-
-    void* d_temp_storage = nullptr;
-    std::size_t temp_storage_bytes = 0;
-    cub::DeviceSpmv::CsrMV(d_temp_storage, temp_storage_bytes, (T*)mat, (AlgInt*)row, (AlgInt*)col,
-                           (T*)px, (T*)py, nrows, ncols, nnz, Gpu::gpuStream());
-    d_temp_storage = (void*)The_Arena()->alloc(temp_storage_bytes);
-    cub::DeviceSpmv::CsrMV(d_temp_storage, temp_storage_bytes, (T*)mat, (AlgInt*)row, (AlgInt*)col,
-                           (T*)px, (T*)py, nrows, ncols, nnz, Gpu::gpuStream());
-    Gpu::streamSynchronize();
-    The_Arena()->free(d_temp_storage);
-
-#else
-
     cusparseHandle_t handle;
     cusparseCreate(&handle);
     cusparseSetStream(handle, Gpu::gpuStream());
 
-    cudaDataType data_type = (sizeof(T) == sizeof(double)) ? CUDA_R_64F : CUDA_R_32F;
-    cusparseIndexType_t index_type = (sizeof(AlgInt) == sizeof(int)) ?
-        CUSPARSE_INDEX_32I : CUSPARSE_INDEX_64I;
+    cudaDataType data_type;
+    if constexpr (std::is_same_v<T,float>) {
+        data_type = CUDA_R_32F;
+    } else if constexpr (std::is_same_v<T,double>) {
+        data_type = CUDA_R_64F;
+    } else if constexpr (std::is_same_v<T,GpuComplex<float>>) {
+        data_type = CUDA_C_32F;
+    } else if constexpr (std::is_same_v<T,GpuComplex<double>>) {
+        data_type = CUDA_C_64F;
+    } else {
+        amrex::Abort("SpMV: unsupported data type");
+    }
+
+    cusparseIndexType_t index_type = CUSPARSE_INDEX_64I;
 
     cusparseSpMatDescr_t mat_descr;
     cusparseCreateCsr(&mat_descr, nrows, ncols, nnz, (void*)row, (void*)col, (void*)mat,
@@ -75,14 +68,14 @@ void SpMV (AlgVector<T>& y, SpMatrix<T> const& A, AlgVector<T> const& x)
     cusparseDnVecDescr_t y_descr;
     cusparseCreateDnVec(&y_descr, nrows, (void*)py, data_type);
 
-    T alpha = T(1.0);
-    T beta = T(0.0);
+    T alpha = T(1);
+    T beta = T(0);
 
     std::size_t buffer_size;
     cusparseSpMV_bufferSize(handle, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha, mat_descr, x_descr,
                             &beta, y_descr, data_type, CUSPARSE_SPMV_ALG_DEFAULT, &buffer_size);
 
-    void* pbuffer = (void*)The_Arena()->alloc(buffer_size);
+    auto* pbuffer = (void*)The_Arena()->alloc(buffer_size);
 
     cusparseSpMV(handle, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha, mat_descr, x_descr,
                  &beta, y_descr, data_type, CUSPARSE_SPMV_ALG_DEFAULT, pbuffer);
@@ -95,17 +88,26 @@ void SpMV (AlgVector<T>& y, SpMatrix<T> const& A, AlgVector<T> const& x)
     cusparseDestroy(handle);
     The_Arena()->free(pbuffer);
 
-#endif
-
 #elif defined(AMREX_USE_HIP)
 
     hipsparseHandle_t handle;
     hipsparseCreate(&handle);
     hipsparseSetStream(handle, Gpu::gpuStream());
 
-    hipDataType data_type = (sizeof(T) == sizeof(double)) ? HIP_R_64F : HIP_R_32F;
-    hipsparseIndexType_t index_type = (sizeof(AlgInt) == sizeof(int)) ?
-        HIPSPARSE_INDEX_32I : HIPSPARSE_INDEX_64I;
+    hipDataType data_type;
+    if constexpr (std::is_same_v<T,float>) {
+        data_type = HIP_R_32F;
+    } else if constexpr (std::is_same_v<T,double>) {
+        data_type = HIP_R_64F;
+    } else if constexpr (std::is_same_v<T,GpuComplex<float>>) {
+        data_type = HIP_C_32F;
+    } else if constexpr (std::is_same_v<T,GpuComplex<double>>) {
+        data_type = HIP_C_64F;
+    } else {
+        amrex::Abort("SpMV: unsupported data type");
+    }
+
+    hipsparseIndexType_t index_type =  HIPSPARSE_INDEX_64I;
 
     hipsparseSpMatDescr_t mat_descr;
     hipsparseCreateCsr(&mat_descr, nrows, ncols, nnz, (void*)row, (void*)col, (void*)mat,
@@ -141,7 +143,7 @@ void SpMV (AlgVector<T>& y, SpMatrix<T> const& A, AlgVector<T> const& x)
 
     mkl::sparse::matrix_handle_t handle{};
     mkl::sparse::set_csr_data(handle, nrows, ncols, mkl::index_base::zero,
-                              (AlgInt*)row, (AlgInt*)col, (T*)mat);
+                              (Long*)row, (Long*)col, (T*)mat);
     mkl::sparse::gemv(Gpu::Device::streamQueue(), mkl::transpose::nontrans,
                       T(1), handle, px, T(0), py);
 
diff --git a/Src/LinearSolvers/AMReX_SpMatrix.H b/Src/LinearSolvers/AMReX_SpMatrix.H
index 610ca899e26..84dca089a75 100644
--- a/Src/LinearSolvers/AMReX_SpMatrix.H
+++ b/Src/LinearSolvers/AMReX_SpMatrix.H
@@ -57,19 +57,21 @@ public:
 
     template <typename U> friend void SpMV(AlgVector<U>& y, SpMatrix<U> const& A, AlgVector<U> const& x);
 
-private:
-
+    //! Private function, but public for cuda
     void define_doit (int nnz);
 
-    void startComm (AlgVector<T> const& x);
-    void finishComm (AlgVector<T>& y);
-
 #ifdef AMREX_USE_MPI
+    //! Private function, but public for cuda
     void prepare_comm ();
     void pack_buffer (AlgVector<T> const& v);
     void unpack_buffer (AlgVector<T>& v);
 #endif
 
+private:
+
+    void startComm (AlgVector<T> const& x);
+    void finishComm (AlgVector<T>& y);
+
     template <class U> using DVec = PODVector<U,Allocator<U> >;
 
     template <template <typename> class V>
@@ -87,15 +89,16 @@ private:
 
 #ifdef AMREX_USE_MPI
     CSR<DVec> m_data_remote;
+
 #ifdef AMREX_USE_GPU
-    template <class U> using HVec = Gpu::PinnedVector<U>;
-    CSR<HVec> m_data_remote_h;
+    Gpu::PinnedVector<Long> m_remote_cols;
 #endif
 
     DVec<Long> m_rtol;
 
     Vector<int> m_send_to;
-    Vector<Vector<Long>> m_send_indices;
+    Vector<int> m_send_counts;
+    Gpu::DeviceVector<Long> m_send_indices;
 
     Vector<int> m_recv_from;
     Vector<int> m_recv_counts;
@@ -203,7 +206,7 @@ void SpMatrix<T,Allocator>::startComm (AlgVector<T> const& x)
 
     auto const nrecvs = int(m_recv_from.size());
     if (nrecvs > 0) {
-        m_recv_buffer = (T*)The_Cpu_Arena()->alloc(sizeof(T)*m_total_counts_recv);
+        m_recv_buffer = (T*)The_Comms_Arena()->alloc(sizeof(T)*m_total_counts_recv);
         m_recv_reqs.resize(nrecvs, MPI_REQUEST_NULL);
         auto* p_recv = m_recv_buffer;
         for (int irecv = 0; irecv < nrecvs; ++irecv) {
@@ -226,7 +229,7 @@ void SpMatrix<T,Allocator>::startComm (AlgVector<T> const& x)
         m_send_reqs.resize(nsends, MPI_REQUEST_NULL);
         auto* p_send = m_send_buffer;
         for (int isend = 0; isend < nsends; ++isend) {
-            auto count = int(m_send_indices[isend].size());
+            auto count = m_send_counts[isend];
             BL_MPI_REQUIRE(MPI_Isend(p_send, count, mpi_t_type, m_send_to[isend],
                                      mpi_tag, mpi_comm, &(m_send_reqs[isend])));
             p_send += count;
@@ -264,7 +267,7 @@ void SpMatrix<T,Allocator>::finishComm (AlgVector<T>& y)
 
     Gpu::streamSynchronize();
     The_Comms_Arena()->free(m_send_buffer);
-    The_Cpu_Arena()->free(m_recv_buffer);
+    The_Comms_Arena()->free(m_recv_buffer);
     m_send_reqs.clear();
     m_recv_reqs.clear();
 #endif
@@ -430,29 +433,18 @@ void SpMatrix<T,Allocator>::prepare_comm ()
         }
 
 #ifdef AMREX_USE_GPU
-        m_data_remote_h.mat.resize(m_data_remote.mat.size());
-        m_data_remote_h.col_index.resize(m_data_remote.col_index.size());
-        m_data_remote_h.row_offset.resize(m_data_remote.row_offset.size());
-        m_data_remote_h.nnz = remote_nnz;
-        Gpu::copyAsync(Gpu::deviceToHost, m_data_remote.mat.begin(),
-                                          m_data_remote.mat.end(),
-                                          m_data_remote_h.mat.begin());
+        m_remote_cols.resize(m_data_remote.col_index.size());
         Gpu::copyAsync(Gpu::deviceToHost, m_data_remote.col_index.begin(),
                                           m_data_remote.col_index.end(),
-                                          m_data_remote_h.col_index.begin());
-        Gpu::copyAsync(Gpu::deviceToHost, m_data_remote.row_offset.begin(),
-                                          m_data_remote.row_offset.end(),
-                                          m_data_remote_h.row_offset.begin());
+                                          m_remote_cols.begin());
         Gpu::streamSynchronize();
-        auto& data_remote = m_data_remote_h;
-        // We could free m_data_remote's vectors if we really want.
 #else
-        auto& data_remote = m_data_remote;
+        auto const& m_remote_cols = m_data_remote.col_index;
 #endif
 
-        unique_remote_cols_v.resize(data_remote.col_index.size());
-        std::partial_sort_copy(data_remote.col_index.begin(),
-                               data_remote.col_index.end(),
+        unique_remote_cols_v.resize(m_remote_cols.size());
+        std::partial_sort_copy(m_remote_cols.begin(),
+                               m_remote_cols.end(),
                                unique_remote_cols_v.begin(),
                                unique_remote_cols_v.end());
         amrex::RemoveDuplicates(unique_remote_cols_v);
@@ -468,7 +460,7 @@ void SpMatrix<T,Allocator>::prepare_comm ()
                 int iproc = int(std::distance(rows.cbegin(),it)) - 1;
                 unique_remote_cols_vv[iproc].push_back(c);
             } else {
-                amrex::Abort("SpMatrix::prepare_comm: how did this happend?");
+                amrex::Abort("SpMatrix::prepare_comm: how did this happen?");
             }
         }
     }
@@ -505,6 +497,7 @@ void SpMatrix<T,Allocator>::prepare_comm ()
         }
     }
 
+    Vector<Vector<Long>> send_indices(nsends);
     m_total_counts_send = 0;
     for (int isend = 0; isend < nsends; ++isend) {
         MPI_Status mpi_status;
@@ -513,12 +506,25 @@ void SpMatrix<T,Allocator>::prepare_comm ()
         int count;
         BL_MPI_REQUIRE(MPI_Get_count(&mpi_status, mpi_long, &count));
         m_send_to.push_back(receiver);
-        m_send_indices.emplace_back(count);
-        BL_MPI_REQUIRE(MPI_Recv(m_send_indices.back().data(), count, mpi_long,
+        m_send_counts.push_back(count);
+        send_indices[isend].resize(count);
+        BL_MPI_REQUIRE(MPI_Recv(send_indices[isend].data(), count, mpi_long,
                                 receiver, mpi_tag, mpi_comm, &mpi_status));
         m_total_counts_send += count;
     }
 
+    m_send_indices.resize(m_total_counts_send);
+    Gpu::PinnedVector<Long> send_indices_all;
+    send_indices_all.reserve(m_total_counts_send);
+    for (auto const& vl : send_indices) {
+        for (auto x : vl) {
+            send_indices_all.push_back(x);
+        }
+    }
+    Gpu::copyAsync(Gpu::hostToDevice, send_indices_all.begin(), send_indices_all.end(),
+                   m_send_indices.begin());
+    Gpu::streamSynchronize();
+
     Vector<MPI_Status> mpi_statuses(mpi_requests.size());
     BL_MPI_REQUIRE(MPI_Waitall(int(mpi_requests.size()), mpi_requests.data(),
                                mpi_statuses.data()));
@@ -529,7 +535,7 @@ void SpMatrix<T,Allocator>::prepare_comm ()
         gtol[unique_remote_cols_v[i]] = i;
     }
 #ifdef AMREX_USE_GPU
-    auto& cols = m_data_remote_h.col_index;
+    auto& cols = m_remote_cols;
 #else
     auto& cols = m_data_remote.col_index;
 #endif
@@ -537,62 +543,50 @@ void SpMatrix<T,Allocator>::prepare_comm ()
         c = gtol[c];
     }
 
+#ifdef AMREX_USE_GPU
+    Gpu::copyAsync(Gpu::hostToDevice, m_remote_cols.begin(), m_remote_cols.end(),
+                   m_data_remote.col_index.data());
+#endif
+
     m_comm_prepared = true;
 }
 
 template <typename T, template<typename> class Allocator>
 void SpMatrix<T,Allocator>::pack_buffer (AlgVector<T> const& v)
 {
-    auto const& vv = v.view();
-    auto const nsends = int(m_send_indices.size());
     auto* pdst = m_send_buffer;
-    for (int isend = 0; isend < nsends; ++isend) {
-        auto count = int(m_send_indices[isend].size());
-        auto& send_indices= m_send_indices[isend];
-#ifdef AMREX_USE_GPU
-        amrex::Abort("xxxxx TODO: SpMatrix::pack_buffer");
-#else
-        for (int i = 0; i < count; ++i) {
-            pdst[i] = vv[send_indices[i]];
-        }
-#endif
-        pdst += count;
-    }
+    auto* pidx = m_send_indices.data();
+    auto const& vv = v.view();
+    auto const nsends = Long(m_send_indices.size());
+    ParallelFor(nsends, [=] AMREX_GPU_DEVICE (Long i)
+    {
+        pdst[i] = vv(pidx[i]);
+    });
 }
 
 template <typename T, template<typename> class Allocator>
 void SpMatrix<T,Allocator>::unpack_buffer (AlgVector<T>& v)
 {
-    // The receive buffer is in CPU memory
-#ifdef AMREX_USE_GPU
-    auto const& csr = m_data_remote_h;
-#else
     auto const& csr = m_data_remote;
-#endif
-
-#ifdef AMREX_USE_GPU
-    amrex::Abort("xxxxx TODO: unpack_buffer");
-#else
     if (csr.nnz > 0) {
         T const* AMREX_RESTRICT mat = csr.mat.data();
         auto const* AMREX_RESTRICT col = csr.col_index.data();
         auto const* AMREX_RESTRICT row = csr.row_offset.data();
 
+        auto const* rtol = m_rtol.data();
+
         auto const* AMREX_RESTRICT px = m_recv_buffer;
         auto      * AMREX_RESTRICT py = v.data();
 
-        auto const nrr = csr.row_offset.size()-1;
-        for (Long i = 0; i < nrr; ++i) {
+        auto const nrr = Long(csr.row_offset.size())-1;
+        ParallelFor(nrr, [=] AMREX_GPU_DEVICE (Long i)
+        {
             T r = 0;
-#ifdef AMREX_USE_OMP
-#pragma omp parallel for reduction(+:r)
-#endif
             for (Long j = row[i]; j < row[i+1]; ++j) {
                 r += mat[j] * px[col[j]];
             }
-            py[m_rtol[i]] += r;
-        }
-#endif
+            py[rtol[i]] += r;
+        });
     }
 }
 
diff --git a/Tools/CMake/AMReXParallelBackends.cmake b/Tools/CMake/AMReXParallelBackends.cmake
index e7f631ca437..1a213fe2198 100644
--- a/Tools/CMake/AMReXParallelBackends.cmake
+++ b/Tools/CMake/AMReXParallelBackends.cmake
@@ -66,10 +66,14 @@ if (  AMReX_GPU_BACKEND STREQUAL "CUDA"
    foreach(D IN LISTS AMReX_SPACEDIM)
        target_link_libraries(amrex_${D}d PUBLIC CUDA::curand)
 
-        # nvToolsExt: if tiny profiler or base profiler are on.
-        if (AMReX_TINY_PROFILE OR AMReX_BASE_PROFILE)
-            target_link_libraries(amrex_${D}d PUBLIC CUDA::nvToolsExt)
-        endif ()
+       if (AMReX_LINEAR_SOLVERS)
+           target_link_libraries(amrex_${D}d PUBLIC CUDA::cusparse)
+       endif ()
+
+       # nvToolsExt: if tiny profiler or base profiler are on.
+       if (AMReX_TINY_PROFILE OR AMReX_BASE_PROFILE)
+           target_link_libraries(amrex_${D}d PUBLIC CUDA::nvToolsExt)
+       endif ()
    endforeach()
 
    # Check cuda compiler and host compiler
diff --git a/Tools/GNUMake/Make.defs b/Tools/GNUMake/Make.defs
index 136a30e7e0e..4cfe205a576 100644
--- a/Tools/GNUMake/Make.defs
+++ b/Tools/GNUMake/Make.defs
@@ -1174,7 +1174,7 @@ else ifeq ($(USE_CUDA),TRUE)
         MAKE_CUDA_PATH := $(SYSTEM_CUDA_PATH)
     endif
 
-    LIBRARIES += -lcuda -lcurand
+    LIBRARIES += -lcuda -lcurand -lcusparse
 
     ifneq ($(MAKE_CUDA_PATH),)
         LIBRARY_LOCATIONS += $(MAKE_CUDA_PATH)/lib64