 	}
 
 	cusparseHandle_t handle = NULL;
-	cusparseSafe(cusparseCreate(&handle), "cuSPARSE create error");
+	TBTKAssert(
+		cusparseCreate(&handle) == CUSPARSE_STATUS_SUCCESS,
+		"ChebyshevSolver::calculateCoefficientsGPU()",
+		"cuSPARSE create error.",
+		""
+	);
 
 	cusparseMatDescr_t descr = NULL;
-	cusparseSafe(cusparseCreateMatDescr(&descr), "cuSPARSE create matrix descriptor error");
-
-	cusparseSafe(cusparseSetMatType(descr, CUSPARSE_MATRIX_TYPE_GENERAL), "cuSPARSE set matrix type error");
-	cusparseSafe(cusparseSetMatIndexBase(descr, CUSPARSE_INDEX_BASE_ZERO), "cuSPARSE set matrix index base error");
-
-	cusparseSafe(cusparseXcoo2csr(handle,
-					cooHARowIndices_device,
-					numHoppingAmplitudes,
-					amplitudeSet->getBasisSize(),
-					csrHARowIndices_device,
-					CUSPARSE_INDEX_BASE_ZERO),
-			"cuSPARSE COO to CSR error");
+	TBTKAssert(
+		cusparseCreateMatDescr(&descr) == CUSPARSE_STATUS_SUCCESS,
+		"ChebyshevSolver::calculateCoefficientsGPU()",
+		"cuSPARSE create matrix descriptor error.",
+		""
+	);
+
+	TBTKAssert(
+		cusparseSetMatType(
+			descr,
+			CUSPARSE_MATRIX_TYPE_GENERAL
+		) == CUSPARSE_STATUS_SUCCESS,
+		"ChebyshevSolver::calculateCoefficientsGPU()",
+		"cuSPARSE set matrix type error.",
+		""
+	);
+	TBTKAssert(
+		cusparseSetMatIndexBase(
+			descr,
+			CUSPARSE_INDEX_BASE_ZERO
+		) == CUSPARSE_STATUS_SUCCESS,
+		"ChebyshevSolver::calculateCoefficientsGPU()",
+		"cuSPARSE set matrix index base error.",
+		""
+	);
+
+	TBTKAssert(
+		cusparseXcoo2csr(
+			handle,
+			cooHARowIndices_device,
+			numHoppingAmplitudes,
+			amplitudeSet->getBasisSize(),
+			csrHARowIndices_device,
+			CUSPARSE_INDEX_BASE_ZERO
+		) == CUSPARSE_STATUS_SUCCESS,
+		"ChebyshevSolver::calculateCoefficientsGPU()",
+		"cuSPARSE COO to CSR error.",
+		""
+	);
 
 	//Calculate |j1>
 	int block_size = 1024;

	}

	complex<double> multiplier = one/scaleFactor;
	cusparseSafe(cusparseZcsrmv(handle,
					CUSPARSE_OPERATION_NON_TRANSPOSE,
					amplitudeSet->getBasisSize(),
					amplitudeSet->getBasisSize(),
					numHoppingAmplitudes,
					(cuDoubleComplex*)&multiplier,
					descr,
					(cuDoubleComplex*)cooHAValues_device,
					csrHARowIndices_device,
					cooHAColIndices_device,
					(cuDoubleComplex*)jIn1_device,
					(cuDoubleComplex*)&zero,
					(cuDoubleComplex*)jIn2_device),
			"Matrix-vector multiplication error");

	extractCoefficients <<< num_blocks, block_size >>> ((cuDoubleComplex*)jIn2_device,
								amplitudeSet->getBasisSize(),

complex<double> two(2., 0.);
complex<double> zero(0., 0.);
complex<double> i(0., 1.);

__global__
void extractCoefficients(
