 	//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(),

		""
	);

	cusparseSafe(cusparseDestroyMatDescr(descr), "cuSPARSE destroy matrix descriptor error");
	descr = NULL;
	cusparseSafe(cusparseDestroy(handle), "cuSPARSE destroy error");
	handle = NULL;

	delete [] jIn1;

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

__global__
void extractCoefficients(
