contrib/benchmarks/GPU/bench_cg_sync.cu - mirror - Git at Google


            
            // Benchmark: GPU Conjugate Gradient via gpu::DeviceMatrix operators.

            
            //

            
            // Shows the path to running Eigen's CG on GPU with minimal code changes.

            
            // The gpu::DeviceMatrix benchmark mirrors Eigen's conjugate_gradient() line-by-line.

            
            // A raw cuBLAS device-pointer-mode implementation is included as a lower bound.

            
            //

            
            // The only change needed in Eigen's CG template to support gpu::DeviceMatrix:

            
            //   Line 34:  typedef Dest VectorType;  (instead of Matrix<Scalar, Dynamic, 1>)

            
            //

            
            // Usage:

            
            //   cmake --build build-bench-gpu --target bench_gpu_cg_sync

            
            //   ./build-bench-gpu/bench_gpu_cg_sync

            
            // SPDX-FileCopyrightText: The Eigen Authors

            
            // SPDX-License-Identifier: MPL-2.0


            #include <benchmark/benchmark.h>


            #include <Eigen/Sparse>

            
            #include <contrib/Eigen/GPU>

            
            #include <cusparse.h>


            using namespace Eigen;


            using Scalar = double;

            
            using RealScalar = double;

            
            using Vec = Matrix<Scalar, Dynamic, 1>;

            
            using SpMat = SparseMatrix<Scalar, ColMajor, int>;


            static SpMat make_spd(Index n) {

            
              SpMat A(n, n);

            
              A.reserve(VectorXi::Constant(n, 3));

            
              for (Index i = 0; i < n; ++i) {

            
                A.insert(i, i) = 4.0;

            
                if (i > 0) A.insert(i, i - 1) = -1.0;

            
                if (i < n - 1) A.insert(i, i + 1) = -1.0;

            
              }

            
              A.makeCompressed();

            
              return A;

            
            }


            static void cuda_warmup() {

            
              static bool done = false;

            
              if (!done) {

            
                void* p;

            
                cudaMalloc(&p, 1);

            
                cudaFree(p);

            
                done = true;

            
              }

            
            }


            // ==========================================================================

            
            // GPU CG using gpu::DeviceMatrix operators — mirrors Eigen's conjugate_gradient()

            
            // ==========================================================================

            
            //

            
            // Compare with Eigen/src/IterativeLinearSolvers/ConjugateGradient.h lines 29-84.

            
            // Left column: Eigen CG code.  Right column: this benchmark.

            
            //

            
            //   Eigen CG                              GPU CG (this benchmark)

            
            //   --------                              -----------------------

            
            //   VectorType residual = rhs - mat * x;  residual.copyFrom(ctx, rhs);  [x=0 so r=b]

            
            //   RealScalar rhsNorm2 = rhs.sqNorm();   RealScalar rhsNorm2 = rhs.squaredNorm();

            
            //   ...

            
            //   tmp.noalias() = mat * p;              tmp.noalias() = mat * p;  [identical]

            
            //   Scalar alpha = absNew / p.dot(tmp);   Scalar alpha = absNew / p.dot(tmp);  [identical]

            
            //   x += alpha * p;                       x += alpha * p;  [identical]

            
            //   residual -= alpha * tmp;              residual -= alpha * tmp;  [identical]

            
            //   residualNorm2 = residual.sqNorm();    residualNorm2 = residual.squaredNorm();  [identical]

            
            //   ...

            
            //   p = z + beta * p;                     p *= beta; p += z;  [equivalent, no alloc]


            static void BM_CG_DeviceMatrixOps(benchmark::State& state) {

            
              cuda_warmup();

            
              const Index n = state.range(0);


              SpMat A = make_spd(n);

            
              Vec b = Vec::Random(n);


              // One shared context: SpMV + BLAS-1 on same stream, zero event overhead.

            
              gpu::Context ctx;

            
              gpu::Context::setThreadLocal(&ctx);

            
              gpu::SparseContext<Scalar> spmv(ctx);

            
              auto mat = spmv.deviceView(A);


              // Upload RHS once.

            
              auto rhs = gpu::DeviceMatrix<Scalar>::fromHost(b, ctx.stream());


              for (auto _ : state) {

            
                // --- Eigen CG lines 34-63: initialization ---

            
                //   typedef Dest VectorType;                               // GPU CHANGE: was Matrix<Scalar,Dynamic,1>

            
                //   VectorType residual = rhs - mat * x;                   // x=0, so residual = rhs

            
                gpu::DeviceMatrix<Scalar> x(n, 1);

            
                x.setZero();

            
                gpu::DeviceMatrix<Scalar> residual(n, 1);

            
                residual.copyFrom(ctx, rhs);


                //   RealScalar rhsNorm2 = rhs.squaredNorm();

            
                RealScalar rhsNorm2 = rhs.squaredNorm();

            
                if (rhsNorm2 == 0) continue;


                RealScalar tol = 1e-10;

            
                const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();

            
                RealScalar threshold = numext::maxi(RealScalar(tol * tol * rhsNorm2), considerAsZero);


                //   RealScalar residualNorm2 = residual.squaredNorm();

            
                RealScalar residualNorm2 = residual.squaredNorm();

            
                if (residualNorm2 < threshold) continue;


                //   VectorType p(n);

            
                //   p = precond.solve(residual);                           // no preconditioner: p = residual

            
                gpu::DeviceMatrix<Scalar> p(n, 1);

            
                p.copyFrom(ctx, residual);


                //   VectorType z(n), tmp(n);

            
                gpu::DeviceMatrix<Scalar> z(n, 1), tmp(n, 1);


                //   auto absNew = numext::real(residual.dot(p));

            
                //   gpu::DeviceScalar — stays on device, no sync.

            
                auto absNew = residual.dot(p);  // gpu::DeviceScalar, no sync


                //   while (i < maxIters) {

            
                Index maxIters = 200;

            
                Index i = 0;

            
                while (i < maxIters) {

            
                  //     tmp.noalias() = mat * p;

            
                  tmp.noalias() = mat * p;  // SpMV, device-resident


                  //     auto alpha = absNew / p.dot(tmp);

            
                  //   gpu::DeviceScalar / gpu::DeviceScalar → device kernel, no sync!

            
                  auto alpha = absNew / p.dot(tmp);  // gpu::DeviceScalar, no sync


                  //     x += alpha * p;

            
                  //   gpu::DeviceScalar * gpu::DeviceMatrix → device-pointer axpy, no sync!

            
                  x += alpha * p;


                  //     residual -= alpha * tmp;

            
                  residual -= alpha * tmp;  // device-pointer axpy, no sync


                  //     residualNorm2 = residual.squaredNorm();

            
                  residualNorm2 = residual.squaredNorm();  // THE one sync per iteration


                  //     if (residualNorm2 < threshold) break;

            
                  if (residualNorm2 < threshold) break;


                  //     z = precond.solve(residual);

            
                  z.copyFrom(ctx, residual);  // no preconditioner


                  //     auto absOld = std::move(absNew);

            
                  auto absOld = std::move(absNew);  // no sync, no alloc


                  //     absNew = numext::real(residual.dot(z));

            
                  absNew = residual.dot(z);  // gpu::DeviceScalar, no sync


                  //     auto beta = absNew / absOld;

            
                  //   gpu::DeviceScalar / gpu::DeviceScalar → device kernel, no sync!

            
                  auto beta = absNew / absOld;  // gpu::DeviceScalar, no sync


                  //     p = z + beta * p;

            
                  p *= beta;  // device-pointer scal, no host sync

            
                  p += z;


                  i++;

            
                }

            
              }


              gpu::Context::setThreadLocal(nullptr);

            
              state.SetItemsProcessed(state.iterations() * 200);

            
            }


            BENCHMARK(BM_CG_DeviceMatrixOps)->RangeMultiplier(4)->Range(1 << 10, 1 << 20)->UseRealTime()->MinWarmUpTime(0.5);


            // ==========================================================================

            
            // Raw cuBLAS device-pointer-mode CG (1 sync/iter) — performance lower bound

            
            // ==========================================================================


            __global__ void scalar_div_kernel(const Scalar* a, const Scalar* b, Scalar* out) { *out = *a / *b; }

            
            __global__ void scalar_neg_kernel(const Scalar* in, Scalar* out) { *out = -(*in); }


            static void BM_CG_DevicePointerMode(benchmark::State& state) {

            
              cuda_warmup();

            
              const Index n = state.range(0);

            
              const int maxIters = 200;


              SpMat A = make_spd(n);

            
              Vec b = Vec::Random(n);


              cudaStream_t stream;

            
              cudaStreamCreate(&stream);

            
              cublasHandle_t cublas;

            
              cublasCreate(&cublas);

            
              cublasSetStream(cublas, stream);


              cusparseHandle_t cusparse;

            
              cusparseCreate(&cusparse);

            
              cusparseSetStream(cusparse, stream);


              gpu::internal::DeviceBuffer d_outer((n + 1) * sizeof(int));

            
              gpu::internal::DeviceBuffer d_inner(A.nonZeros() * sizeof(int));

            
              gpu::internal::DeviceBuffer d_vals(A.nonZeros() * sizeof(Scalar));

            
              cudaMemcpy(d_outer.get(), A.outerIndexPtr(), (n + 1) * sizeof(int), cudaMemcpyHostToDevice);

            
              cudaMemcpy(d_inner.get(), A.innerIndexPtr(), A.nonZeros() * sizeof(int), cudaMemcpyHostToDevice);

            
              cudaMemcpy(d_vals.get(), A.valuePtr(), A.nonZeros() * sizeof(Scalar), cudaMemcpyHostToDevice);


              cusparseSpMatDescr_t matA;

            
              cusparseCreateCsc(&matA, n, n, A.nonZeros(), d_outer.get(), d_inner.get(), d_vals.get(), CUSPARSE_INDEX_32I,

            
                                CUSPARSE_INDEX_32I, CUSPARSE_INDEX_BASE_ZERO, CUDA_R_64F);


              gpu::internal::DeviceBuffer d_tmp_buf(n * sizeof(Scalar));

            
              cusparseDnVecDescr_t tmp_x, tmp_y;

            
              cusparseCreateDnVec(&tmp_x, n, d_tmp_buf.get(), CUDA_R_64F);

            
              cusparseCreateDnVec(&tmp_y, n, d_tmp_buf.get(), CUDA_R_64F);

            
              Scalar spmv_alpha = 1.0, spmv_beta = 0.0;

            
              size_t ws_size = 0;

            
              cusparseSpMV_bufferSize(cusparse, CUSPARSE_OPERATION_NON_TRANSPOSE, &spmv_alpha, matA, tmp_x, &spmv_beta, tmp_y,

            
                                      CUDA_R_64F, CUSPARSE_SPMV_ALG_DEFAULT, &ws_size);

            
              gpu::internal::DeviceBuffer d_workspace(ws_size);

            
              cusparseDestroyDnVec(tmp_x);

            
              cusparseDestroyDnVec(tmp_y);


              gpu::internal::DeviceBuffer d_x(n * sizeof(Scalar)), d_r(n * sizeof(Scalar));

            
              gpu::internal::DeviceBuffer d_p(n * sizeof(Scalar)), d_tmp(n * sizeof(Scalar));

            
              gpu::internal::DeviceBuffer d_b(n * sizeof(Scalar));

            
              gpu::internal::DeviceBuffer d_absNew(sizeof(Scalar)), d_absOld(sizeof(Scalar));

            
              gpu::internal::DeviceBuffer d_pdot(sizeof(Scalar)), d_alpha(sizeof(Scalar));

            
              gpu::internal::DeviceBuffer d_neg_alpha(sizeof(Scalar)), d_beta(sizeof(Scalar));

            
              gpu::internal::DeviceBuffer d_rnorm(sizeof(RealScalar));


              cudaMemcpy(d_b.get(), b.data(), n * sizeof(Scalar), cudaMemcpyHostToDevice);


              auto spmv = [&](Scalar* x_ptr, Scalar* y_ptr) {

            
                cusparseDnVecDescr_t vx, vy;

            
                cusparseCreateDnVec(&vx, n, x_ptr, CUDA_R_64F);

            
                cusparseCreateDnVec(&vy, n, y_ptr, CUDA_R_64F);

            
                cusparseSpMV(cusparse, CUSPARSE_OPERATION_NON_TRANSPOSE, &spmv_alpha, matA, vx, &spmv_beta, vy, CUDA_R_64F,

            
                             CUSPARSE_SPMV_ALG_DEFAULT, d_workspace.get());

            
                cusparseDestroyDnVec(vx);

            
                cusparseDestroyDnVec(vy);

            
              };


              for (auto _ : state) {

            
                cudaMemsetAsync(static_cast<Scalar*>(d_x.get()), 0, n * sizeof(Scalar), stream);

            
                cudaMemcpyAsync(d_r.get(), d_b.get(), n * sizeof(Scalar), cudaMemcpyDeviceToDevice, stream);

            
                cudaMemcpyAsync(d_p.get(), d_b.get(), n * sizeof(Scalar), cudaMemcpyDeviceToDevice, stream);


                cublasSetPointerMode(cublas, CUBLAS_POINTER_MODE_DEVICE);

            
                cublasDdot(cublas, n, static_cast<Scalar*>(d_r.get()), 1, static_cast<Scalar*>(d_p.get()), 1,

            
                           static_cast<Scalar*>(d_absNew.get()));


                for (int i = 0; i < maxIters; ++i) {

            
                  spmv(static_cast<Scalar*>(d_p.get()), static_cast<Scalar*>(d_tmp.get()));


                  cublasDdot(cublas, n, static_cast<Scalar*>(d_p.get()), 1, static_cast<Scalar*>(d_tmp.get()), 1,

            
                             static_cast<Scalar*>(d_pdot.get()));


                  scalar_div_kernel<<<1, 1, 0, stream>>>(static_cast<Scalar*>(d_absNew.get()), static_cast<Scalar*>(d_pdot.get()),

            
                                                         static_cast<Scalar*>(d_alpha.get()));

            
                  scalar_neg_kernel<<<1, 1, 0, stream>>>(static_cast<Scalar*>(d_alpha.get()),

            
                                                         static_cast<Scalar*>(d_neg_alpha.get()));


                  cublasDaxpy(cublas, n, static_cast<Scalar*>(d_alpha.get()), static_cast<Scalar*>(d_p.get()), 1,

            
                              static_cast<Scalar*>(d_x.get()), 1);

            
                  cublasDaxpy(cublas, n, static_cast<Scalar*>(d_neg_alpha.get()), static_cast<Scalar*>(d_tmp.get()), 1,

            
                              static_cast<Scalar*>(d_r.get()), 1);


                  cublasDnrm2(cublas, n, static_cast<Scalar*>(d_r.get()), 1, static_cast<RealScalar*>(d_rnorm.get()));


                  RealScalar rnorm;

            
                  cudaMemcpyAsync(&rnorm, d_rnorm.get(), sizeof(RealScalar), cudaMemcpyDeviceToHost, stream);

            
                  cudaStreamSynchronize(stream);

            
                  if (rnorm * rnorm < 1e-20) break;


                  cudaMemcpyAsync(d_absOld.get(), d_absNew.get(), sizeof(Scalar), cudaMemcpyDeviceToDevice, stream);

            
                  cublasDdot(cublas, n, static_cast<Scalar*>(d_r.get()), 1, static_cast<Scalar*>(d_r.get()), 1,

            
                             static_cast<Scalar*>(d_absNew.get()));


                  scalar_div_kernel<<<1, 1, 0, stream>>>(static_cast<Scalar*>(d_absNew.get()), static_cast<Scalar*>(d_absOld.get()),

            
                                                         static_cast<Scalar*>(d_beta.get()));


                  cublasDscal(cublas, n, static_cast<Scalar*>(d_beta.get()), static_cast<Scalar*>(d_p.get()), 1);

            
                  cublasSetPointerMode(cublas, CUBLAS_POINTER_MODE_HOST);

            
                  Scalar one = 1.0;

            
                  cublasDaxpy(cublas, n, &one, static_cast<Scalar*>(d_r.get()), 1, static_cast<Scalar*>(d_p.get()), 1);

            
                  cublasSetPointerMode(cublas, CUBLAS_POINTER_MODE_DEVICE);

            
                }

            
                cudaStreamSynchronize(stream);

            
              }


              state.SetItemsProcessed(state.iterations() * maxIters);

            
              cusparseDestroySpMat(matA);

            
              cusparseDestroy(cusparse);

            
              cublasDestroy(cublas);

            
              cudaStreamDestroy(stream);

            
            }


            BENCHMARK(BM_CG_DevicePointerMode)->RangeMultiplier(4)->Range(1 << 10, 1 << 20)->UseRealTime()->MinWarmUpTime(0.5);
	// Benchmark: GPU Conjugate Gradient via gpu::DeviceMatrix operators.
	//
	// Shows the path to running Eigen's CG on GPU with minimal code changes.
	// The gpu::DeviceMatrix benchmark mirrors Eigen's conjugate_gradient() line-by-line.
	// A raw cuBLAS device-pointer-mode implementation is included as a lower bound.
	//
	// The only change needed in Eigen's CG template to support gpu::DeviceMatrix:
	// Line 34: typedef Dest VectorType; (instead of Matrix<Scalar, Dynamic, 1>)
	//
	// Usage:
	// cmake --build build-bench-gpu --target bench_gpu_cg_sync
	// ./build-bench-gpu/bench_gpu_cg_sync
	// SPDX-FileCopyrightText: The Eigen Authors
	// SPDX-License-Identifier: MPL-2.0

	#include <benchmark/benchmark.h>

	#include <Eigen/Sparse>
	#include <contrib/Eigen/GPU>
	#include <cusparse.h>

	using namespace Eigen;

	using Scalar = double;
	using RealScalar = double;
	using Vec = Matrix<Scalar, Dynamic, 1>;
	using SpMat = SparseMatrix<Scalar, ColMajor, int>;

	static SpMat make_spd(Index n) {
	SpMat A(n, n);
	A.reserve(VectorXi::Constant(n, 3));
	for (Index i = 0; i < n; ++i) {
	A.insert(i, i) = 4.0;
	if (i > 0) A.insert(i, i - 1) = -1.0;
	if (i < n - 1) A.insert(i, i + 1) = -1.0;
	}
	A.makeCompressed();
	return A;
	}

	static void cuda_warmup() {
	static bool done = false;
	if (!done) {
	void* p;
	cudaMalloc(&p, 1);
	cudaFree(p);
	done = true;
	}
	}

	// ==========================================================================
	// GPU CG using gpu::DeviceMatrix operators — mirrors Eigen's conjugate_gradient()
	// ==========================================================================
	//
	// Compare with Eigen/src/IterativeLinearSolvers/ConjugateGradient.h lines 29-84.
	// Left column: Eigen CG code. Right column: this benchmark.
	//
	// Eigen CG GPU CG (this benchmark)
	// -------- -----------------------
	// VectorType residual = rhs - mat * x; residual.copyFrom(ctx, rhs); [x=0 so r=b]
	// RealScalar rhsNorm2 = rhs.sqNorm(); RealScalar rhsNorm2 = rhs.squaredNorm();
	// ...
	// tmp.noalias() = mat * p; tmp.noalias() = mat * p; [identical]
	// Scalar alpha = absNew / p.dot(tmp); Scalar alpha = absNew / p.dot(tmp); [identical]
	// x += alpha * p; x += alpha * p; [identical]
	// residual -= alpha * tmp; residual -= alpha * tmp; [identical]
	// residualNorm2 = residual.sqNorm(); residualNorm2 = residual.squaredNorm(); [identical]
	// ...
	// p = z + beta * p; p *= beta; p += z; [equivalent, no alloc]

	static void BM_CG_DeviceMatrixOps(benchmark::State& state) {
	cuda_warmup();
	const Index n = state.range(0);

	SpMat A = make_spd(n);
	Vec b = Vec::Random(n);

	// One shared context: SpMV + BLAS-1 on same stream, zero event overhead.
	gpu::Context ctx;
	gpu::Context::setThreadLocal(&ctx);
	gpu::SparseContext<Scalar> spmv(ctx);
	auto mat = spmv.deviceView(A);

	// Upload RHS once.
	auto rhs = gpu::DeviceMatrix<Scalar>::fromHost(b, ctx.stream());

	for (auto _ : state) {
	// --- Eigen CG lines 34-63: initialization ---
	// typedef Dest VectorType; // GPU CHANGE: was Matrix<Scalar,Dynamic,1>
	// VectorType residual = rhs - mat * x; // x=0, so residual = rhs
	gpu::DeviceMatrix<Scalar> x(n, 1);
	x.setZero();
	gpu::DeviceMatrix<Scalar> residual(n, 1);
	residual.copyFrom(ctx, rhs);

	// RealScalar rhsNorm2 = rhs.squaredNorm();
	RealScalar rhsNorm2 = rhs.squaredNorm();
	if (rhsNorm2 == 0) continue;

	RealScalar tol = 1e-10;
	const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
	RealScalar threshold = numext::maxi(RealScalar(tol * tol * rhsNorm2), considerAsZero);

	// RealScalar residualNorm2 = residual.squaredNorm();
	RealScalar residualNorm2 = residual.squaredNorm();
	if (residualNorm2 < threshold) continue;

	// VectorType p(n);
	// p = precond.solve(residual); // no preconditioner: p = residual
	gpu::DeviceMatrix<Scalar> p(n, 1);
	p.copyFrom(ctx, residual);

	// VectorType z(n), tmp(n);
	gpu::DeviceMatrix<Scalar> z(n, 1), tmp(n, 1);

	// auto absNew = numext::real(residual.dot(p));
	// gpu::DeviceScalar — stays on device, no sync.
	auto absNew = residual.dot(p); // gpu::DeviceScalar, no sync

	// while (i < maxIters) {
	Index maxIters = 200;
	Index i = 0;
	while (i < maxIters) {
	// tmp.noalias() = mat * p;
	tmp.noalias() = mat * p; // SpMV, device-resident

	// auto alpha = absNew / p.dot(tmp);
	// gpu::DeviceScalar / gpu::DeviceScalar → device kernel, no sync!
	auto alpha = absNew / p.dot(tmp); // gpu::DeviceScalar, no sync

	// x += alpha * p;
	// gpu::DeviceScalar * gpu::DeviceMatrix → device-pointer axpy, no sync!
	x += alpha * p;

	// residual -= alpha * tmp;
	residual -= alpha * tmp; // device-pointer axpy, no sync

	// residualNorm2 = residual.squaredNorm();
	residualNorm2 = residual.squaredNorm(); // THE one sync per iteration

	// if (residualNorm2 < threshold) break;
	if (residualNorm2 < threshold) break;

	// z = precond.solve(residual);
	z.copyFrom(ctx, residual); // no preconditioner

	// auto absOld = std::move(absNew);
	auto absOld = std::move(absNew); // no sync, no alloc

	// absNew = numext::real(residual.dot(z));
	absNew = residual.dot(z); // gpu::DeviceScalar, no sync

	// auto beta = absNew / absOld;
	// gpu::DeviceScalar / gpu::DeviceScalar → device kernel, no sync!
	auto beta = absNew / absOld; // gpu::DeviceScalar, no sync

	// p = z + beta * p;
	p *= beta; // device-pointer scal, no host sync
	p += z;

	i++;
	}
	}

	gpu::Context::setThreadLocal(nullptr);
	state.SetItemsProcessed(state.iterations() * 200);
	}

	BENCHMARK(BM_CG_DeviceMatrixOps)->RangeMultiplier(4)->Range(1 << 10, 1 << 20)->UseRealTime()->MinWarmUpTime(0.5);

	// ==========================================================================
	// Raw cuBLAS device-pointer-mode CG (1 sync/iter) — performance lower bound
	// ==========================================================================

	__global__ void scalar_div_kernel(const Scalar* a, const Scalar* b, Scalar* out) { out = a / *b; }
	__global__ void scalar_neg_kernel(const Scalar* in, Scalar* out) { out = -(in); }

	static void BM_CG_DevicePointerMode(benchmark::State& state) {
	cuda_warmup();
	const Index n = state.range(0);
	const int maxIters = 200;

	SpMat A = make_spd(n);
	Vec b = Vec::Random(n);

	cudaStream_t stream;
	cudaStreamCreate(&stream);
	cublasHandle_t cublas;
	cublasCreate(&cublas);
	cublasSetStream(cublas, stream);

	cusparseHandle_t cusparse;
	cusparseCreate(&cusparse);
	cusparseSetStream(cusparse, stream);

	gpu::internal::DeviceBuffer d_outer((n + 1) * sizeof(int));
	gpu::internal::DeviceBuffer d_inner(A.nonZeros() * sizeof(int));
	gpu::internal::DeviceBuffer d_vals(A.nonZeros() * sizeof(Scalar));
	cudaMemcpy(d_outer.get(), A.outerIndexPtr(), (n + 1) * sizeof(int), cudaMemcpyHostToDevice);
	cudaMemcpy(d_inner.get(), A.innerIndexPtr(), A.nonZeros() * sizeof(int), cudaMemcpyHostToDevice);
	cudaMemcpy(d_vals.get(), A.valuePtr(), A.nonZeros() * sizeof(Scalar), cudaMemcpyHostToDevice);

	cusparseSpMatDescr_t matA;
	cusparseCreateCsc(&matA, n, n, A.nonZeros(), d_outer.get(), d_inner.get(), d_vals.get(), CUSPARSE_INDEX_32I,
	CUSPARSE_INDEX_32I, CUSPARSE_INDEX_BASE_ZERO, CUDA_R_64F);

	gpu::internal::DeviceBuffer d_tmp_buf(n * sizeof(Scalar));
	cusparseDnVecDescr_t tmp_x, tmp_y;
	cusparseCreateDnVec(&tmp_x, n, d_tmp_buf.get(), CUDA_R_64F);
	cusparseCreateDnVec(&tmp_y, n, d_tmp_buf.get(), CUDA_R_64F);
	Scalar spmv_alpha = 1.0, spmv_beta = 0.0;
	size_t ws_size = 0;
	cusparseSpMV_bufferSize(cusparse, CUSPARSE_OPERATION_NON_TRANSPOSE, &spmv_alpha, matA, tmp_x, &spmv_beta, tmp_y,
	CUDA_R_64F, CUSPARSE_SPMV_ALG_DEFAULT, &ws_size);
	gpu::internal::DeviceBuffer d_workspace(ws_size);
	cusparseDestroyDnVec(tmp_x);
	cusparseDestroyDnVec(tmp_y);

	gpu::internal::DeviceBuffer d_x(n * sizeof(Scalar)), d_r(n * sizeof(Scalar));
	gpu::internal::DeviceBuffer d_p(n * sizeof(Scalar)), d_tmp(n * sizeof(Scalar));
	gpu::internal::DeviceBuffer d_b(n * sizeof(Scalar));
	gpu::internal::DeviceBuffer d_absNew(sizeof(Scalar)), d_absOld(sizeof(Scalar));
	gpu::internal::DeviceBuffer d_pdot(sizeof(Scalar)), d_alpha(sizeof(Scalar));
	gpu::internal::DeviceBuffer d_neg_alpha(sizeof(Scalar)), d_beta(sizeof(Scalar));
	gpu::internal::DeviceBuffer d_rnorm(sizeof(RealScalar));

	cudaMemcpy(d_b.get(), b.data(), n * sizeof(Scalar), cudaMemcpyHostToDevice);

	auto spmv = [&](Scalar* x_ptr, Scalar* y_ptr) {
	cusparseDnVecDescr_t vx, vy;
	cusparseCreateDnVec(&vx, n, x_ptr, CUDA_R_64F);
	cusparseCreateDnVec(&vy, n, y_ptr, CUDA_R_64F);
	cusparseSpMV(cusparse, CUSPARSE_OPERATION_NON_TRANSPOSE, &spmv_alpha, matA, vx, &spmv_beta, vy, CUDA_R_64F,
	CUSPARSE_SPMV_ALG_DEFAULT, d_workspace.get());
	cusparseDestroyDnVec(vx);
	cusparseDestroyDnVec(vy);
	};

	for (auto _ : state) {
	cudaMemsetAsync(static_cast<Scalar>(d_x.get()), 0, n sizeof(Scalar), stream);
	cudaMemcpyAsync(d_r.get(), d_b.get(), n * sizeof(Scalar), cudaMemcpyDeviceToDevice, stream);
	cudaMemcpyAsync(d_p.get(), d_b.get(), n * sizeof(Scalar), cudaMemcpyDeviceToDevice, stream);

	cublasSetPointerMode(cublas, CUBLAS_POINTER_MODE_DEVICE);
	cublasDdot(cublas, n, static_cast<Scalar>(d_r.get()), 1, static_cast<Scalar>(d_p.get()), 1,
	static_cast<Scalar*>(d_absNew.get()));

	for (int i = 0; i < maxIters; ++i) {
	spmv(static_cast<Scalar>(d_p.get()), static_cast<Scalar>(d_tmp.get()));

	cublasDdot(cublas, n, static_cast<Scalar>(d_p.get()), 1, static_cast<Scalar>(d_tmp.get()), 1,
	static_cast<Scalar*>(d_pdot.get()));

	scalar_div_kernel<<<1, 1, 0, stream>>>(static_cast<Scalar>(d_absNew.get()), static_cast<Scalar>(d_pdot.get()),
	static_cast<Scalar*>(d_alpha.get()));
	scalar_neg_kernel<<<1, 1, 0, stream>>>(static_cast<Scalar*>(d_alpha.get()),
	static_cast<Scalar*>(d_neg_alpha.get()));

	cublasDaxpy(cublas, n, static_cast<Scalar>(d_alpha.get()), static_cast<Scalar>(d_p.get()), 1,
	static_cast<Scalar*>(d_x.get()), 1);
	cublasDaxpy(cublas, n, static_cast<Scalar>(d_neg_alpha.get()), static_cast<Scalar>(d_tmp.get()), 1,
	static_cast<Scalar*>(d_r.get()), 1);

	cublasDnrm2(cublas, n, static_cast<Scalar>(d_r.get()), 1, static_cast<RealScalar>(d_rnorm.get()));

	RealScalar rnorm;
	cudaMemcpyAsync(&rnorm, d_rnorm.get(), sizeof(RealScalar), cudaMemcpyDeviceToHost, stream);
	cudaStreamSynchronize(stream);
	if (rnorm * rnorm < 1e-20) break;

	cudaMemcpyAsync(d_absOld.get(), d_absNew.get(), sizeof(Scalar), cudaMemcpyDeviceToDevice, stream);
	cublasDdot(cublas, n, static_cast<Scalar>(d_r.get()), 1, static_cast<Scalar>(d_r.get()), 1,
	static_cast<Scalar*>(d_absNew.get()));

	scalar_div_kernel<<<1, 1, 0, stream>>>(static_cast<Scalar>(d_absNew.get()), static_cast<Scalar>(d_absOld.get()),
	static_cast<Scalar*>(d_beta.get()));

	cublasDscal(cublas, n, static_cast<Scalar>(d_beta.get()), static_cast<Scalar>(d_p.get()), 1);
	cublasSetPointerMode(cublas, CUBLAS_POINTER_MODE_HOST);
	Scalar one = 1.0;
	cublasDaxpy(cublas, n, &one, static_cast<Scalar>(d_r.get()), 1, static_cast<Scalar>(d_p.get()), 1);
	cublasSetPointerMode(cublas, CUBLAS_POINTER_MODE_DEVICE);
	}
	cudaStreamSynchronize(stream);
	}

	state.SetItemsProcessed(state.iterations() * maxIters);
	cusparseDestroySpMat(matA);
	cusparseDestroy(cusparse);
	cublasDestroy(cublas);
	cudaStreamDestroy(stream);
	}

	BENCHMARK(BM_CG_DevicePointerMode)->RangeMultiplier(4)->Range(1 << 10, 1 << 20)->UseRealTime()->MinWarmUpTime(0.5);