 	//Iteratively calculate |jn> and corresponding Chebyshev coefficients.
 	for(int n = 2; n < numCoefficients; n++){
 		multiplier = two/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*)&minus_one,
-						(cuDoubleComplex*)jIn2_device),
-				"Matrix-vector multiplication error");
+		TBTKAssert(
+			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*)&minus_one,
+				(cuDoubleComplex*)jIn2_device
+			) == CUSPARSE_STATUS_SUCCESS,
+			"ChebyshevSolver::calculateCoefficientsGPU()",
+			"Matrix-vector multiplication error.",
+			""
+		);
 
 		extractCoefficients <<< num_blocks, block_size >>> ((cuDoubleComplex*)jIn2_device,
 									amplitudeSet->getBasisSize(),

	}

	cusparseHandle_t handle = NULL;
	cusparseSafe(cusparseCreate(&handle), "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");

	//Calculate |j1>
	int block_size = 1024;

	}

	complex<double> multiplier = one/scaleFactor;
	TBTKAssert(
		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
		) == CUSPARSE_STATUS_SUCCESS,
		"ChebyshevSolver::calculateCoefficentsGPU()",
		"Matrix-vector multiplication error.",
		""
	);

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