 }
 #define SAFE_CALL(call,msg) _safe_cuda_call((call),(msg),__FILE__,__LINE__)
 
+#ifndef cublasSafeCall
+#define cublasSafeCall(err) __cublasSafeCall(err, __FILE__, __LINE__)
+#endif
+inline void __cublasSafeCall(cublasStatus_t err, const char *file, const int line)
+{
+if( CUBLAS_STATUS_SUCCESS != err) {
+fprintf(stderr, "CUBLAS error in file '%s', line %d\n \nerror %d \nterminating!\n",__FILE__, __LINE__,err);
+cudaDeviceReset(); assert(0);
+}
+}
 // Apply global contrast normalization to
 // each input image.
 // For each channel in each image, the mean and stddev has

		mean_stddev_reduction_kernel12<<<1, input.size(),0,sa.getStream(stream)>>>(d_M1, d_M2, d_n, dMean, dStddev);
	else
		mean_stddev_reduction_kernel24<<<1, input.size(),0,sa.getStream(stream)>>>(d_M1, d_M2, d_n, dMean, dStddev);

	// Convert each input image held in GpuMats pointed to
	// by dPssIn into a flattened format.  Each 2-D image

	// image has the mean and stddev computed individually
	// in the CUDA kernels above
	global_contrast_normalization_kernel<<<grid,block,0,sa.getStream(stream)>>>(dPssIn, dMean, dStddev, dFlattenedImages);
	SAFE_CALL(cudaStreamSynchronize(sa.getStream(stream)),"ZCA cudaStreamSynchronize failed");

	// Multiply images by weights to get the ZCA-whitened output
	//gemm(dFlattenedImages, weights, 1.0, cv::gpu::GpuMat(), 0.0, zcaOut, 0, stream);

	zcaOut.create(dFlattenedImages.size(), dFlattenedImages.type());
    cublasHandle_t handle;
    cublasSafeCall( cublasCreate_v2(&handle) );

    cublasSafeCall( cublasSetStream_v2(handle, sa.getStream(stream)) );

    cublasSafeCall( cublasSetPointerMode_v2(handle, CUBLAS_POINTER_MODE_HOST) );
    const float alphaf = 1.0;
    const float betaf = 0.0;
        cublasSafeCall( cublasSgemm_v2(handle, CUBLAS_OP_N, CUBLAS_OP_N, weights.cols, dFlattenedImages.rows, weights.rows,
            &alphaf,
            weights.ptr<float>(), static_cast<int>(weights.step / sizeof(float)),
            dFlattenedImages.ptr<float>(), static_cast<int>(dFlattenedImages.step / sizeof(float)),
            &betaf,
            zcaOut.ptr<float>(), static_cast<int>(zcaOut.step / sizeof(float))) );
    cublasSafeCall( cublasDestroy_v2(handle) );

	// Copy to output buffer in the order expected by
	// neural net input
