From e70a251fdd863e281bbf3ae30319e4ed8ae95a97 Mon Sep 17 00:00:00 2001
From: dmayerich <dmayerich@gmail.com>
Date: Wed, 11 Sep 2013 22:13:02 -0500
Subject: [PATCH] fixed double precision errors

---
 fieldslice.cpp   | 28 ++++++++++++++++++++++------
 scalarslice.cu   | 13 +++++++++++++
 win32/bimsim.exe | Bin 449536 -> 0 bytes
 3 files changed, 35 insertions(+), 6 deletions(-)

diff --git a/fieldslice.cpp b/fieldslice.cpp
index 5df3b55..bf63766 100644
--- a/fieldslice.cpp
+++ b/fieldslice.cpp
@@ -23,19 +23,27 @@ fieldslice::fieldslice(unsigned int x_size, unsigned int y_size)
 void fieldslice::toAngularSpectrum()
 {
     cufftHandle plan;
+	cufftResult result;
 
     //create a CUFFT plan handle
-    if(cufftPlan2d(&plan, R[0], R[1], CUFFT_C2C) != CUFFT_SUCCESS)
+#ifdef PRECISION_SINGLE
+    result = cufftPlan2d(&plan, R[0], R[1], CUFFT_C2C);
+#elif defined PRECISION_DOUBLE
+	 result = cufftPlan2d(&plan, R[0], R[1], CUFFT_Z2Z);
+#endif
+
+	if(result != CUFFT_SUCCESS)
     {
         cout<<"Error creating CUFFT plan for computing the angular spectrum."<<endl;
         exit(1);
     }
 
 #ifdef PRECISION_SINGLE
-    if(cufftExecC2C(plan, (cufftComplex*)x_hat, (cufftComplex*)x_hat, CUFFT_FORWARD) != CUFFT_SUCCESS)
+    result = cufftExecC2C(plan, (cufftComplex*)x_hat, (cufftComplex*)x_hat, CUFFT_FORWARD);
 #elif defined PRECISION_DOUBLE
-    if(cufftExecZ2Z(plan, (cufftDoubleComplex*)x_hat, (cufftDoubleComplex*)x_hat, CUFFT_FORWARD) != CUFFT_SUCCESS)
+    result = cufftExecZ2Z(plan, (cufftDoubleComplex*)x_hat, (cufftDoubleComplex*)x_hat, CUFFT_FORWARD);
 #endif
+	if(result != CUFFT_SUCCESS)
     {
         cout<<"Error executing the CUFFT forward FFT to compute the angular spectrum."<<endl;
         exit(1);
@@ -49,19 +57,27 @@ void fieldslice::toAngularSpectrum()
 void fieldslice::fromAngularSpectrum()
 {
     cufftHandle plan;
+	cufftResult result;
 
     //create a CUFFT plan handle
-    if(cufftPlan2d(&plan, R[0], R[1], CUFFT_C2C) != CUFFT_SUCCESS)
+#ifdef PRECISION_SINGLE
+    result = cufftPlan2d(&plan, R[0], R[1], CUFFT_C2C);
+#elif defined PRECISION_DOUBLE
+	 result = cufftPlan2d(&plan, R[0], R[1], CUFFT_Z2Z);
+#endif
+
+	if(result != CUFFT_SUCCESS)
     {
         cout<<"Error creating CUFFT plan for computing the angular spectrum."<<endl;
         exit(1);
     }
 
 #ifdef PRECISION_SINGLE
-    if(cufftExecC2C(plan, (cufftComplex*)x_hat, (cufftComplex*)x_hat, CUFFT_INVERSE) != CUFFT_SUCCESS)
+    result = cufftExecC2C(plan, (cufftComplex*)x_hat, (cufftComplex*)x_hat, CUFFT_INVERSE);
 #elif defined PRECISION_DOUBLE
-    if(cufftExecZ2Z(plan, (cufftDoubleComplex*)x_hat, (cufftDoubleComplex*)x_hat, CUFFT_INVERSE) != CUFFT_SUCCESS)
+    result = cufftExecZ2Z(plan, (cufftDoubleComplex*)x_hat, (cufftDoubleComplex*)x_hat, CUFFT_INVERSE);
 #endif
+	if(result != CUFFT_SUCCESS)
     {
         cout<<"Error executing the CUFFT forward FFT to compute the angular spectrum."<<endl;
         exit(1);
diff --git a/scalarslice.cu b/scalarslice.cu
index 00de57e..daa609d 100644
--- a/scalarslice.cu
+++ b/scalarslice.cu
@@ -47,7 +47,20 @@ void scalarslice::toImage(std::string filename, bool positive, rts::colormap::co
     //find the index of the value with maximum magnitude
     int N = R[0] * R[1];
     int result;
+#ifdef PRECISION_SINGLE
     stat = cublasIsamax(handle, N, S, 1, &result);
+#elif defined PRECISION_DOUBLE
+	stat = cublasIdamax(handle, N, S, 1, &result);
+#endif
+
+	//adjust for 1-based indexing
+	result -= 1;
+
+	if(stat != CUBLAS_STATUS_SUCCESS)
+	{
+		std::cout<<"CUBLAS Error: failure finding maximum value."<<std::endl;
+		exit(1);
+	}
 
     //std::cout<<"Maximum index: "<<result<<std::endl;
 
diff --git a/win32/bimsim.exe b/win32/bimsim.exe
index a162026..6d98c90 100644
Binary files a/win32/bimsim.exe and b/win32/bimsim.exe differ
--
libgit2 0.21.4