fix bug #189 (issue with fortran concentions to return COMPLEX values)

2 files changed, 60 insertions, 16 deletions

diff --git a/blas/complexdots.f b/blas/complexdots.f
new file mode 100644
index 000000000..a7da51d16
--- /dev/null
+++ b/blas/complexdots.f
@@ -0,0 +1,43 @@
+      COMPLEX FUNCTION CDOTC(N,CX,INCX,CY,INCY)
+      INTEGER INCX,INCY,N
+      COMPLEX CX(*),CY(*)
+      COMPLEX RES
+      EXTERNAL CDOTCW
+      
+      CALL CDOTCW(N,CX,INCX,CY,INCY,RES)
+      CDOTC = RES
+      RETURN
+      END
+      
+      COMPLEX FUNCTION CDOTU(N,CX,INCX,CY,INCY)
+      INTEGER INCX,INCY,N
+      COMPLEX CX(*),CY(*)
+      COMPLEX RES
+      EXTERNAL CDOTUW
+      
+      CALL CDOTUW(N,CX,INCX,CY,INCY,RES)
+      CDOTU = RES
+      RETURN
+      END
+      
+      DOUBLE COMPLEX FUNCTION ZDOTC(N,CX,INCX,CY,INCY)
+      INTEGER INCX,INCY,N
+      DOUBLE COMPLEX CX(*),CY(*)
+      DOUBLE COMPLEX RES
+      EXTERNAL ZDOTCW
+      
+      CALL ZDOTCW(N,CX,INCX,CY,INCY,RES)
+      ZDOTC = RES
+      RETURN
+      END
+      
+      DOUBLE COMPLEX FUNCTION ZDOTU(N,CX,INCX,CY,INCY)
+      INTEGER INCX,INCY,N
+      DOUBLE COMPLEX CX(*),CY(*)
+      DOUBLE COMPLEX RES
+      EXTERNAL ZDOTUW
+      
+      CALL ZDOTUW(N,CX,INCX,CY,INCY,RES)
+      ZDOTU = RES
+      RETURN
+      END
diff --git a/blas/level1_cplx_impl.h b/blas/level1_cplx_impl.h
index e268bb812..299a4ac00 100644
--- a/blas/level1_cplx_impl.h
+++ b/blas/level1_cplx_impl.h
@@ -52,7 +52,7 @@ RealScalar EIGEN_CAT(EIGEN_CAT(REAL_SCALAR_SUFFIX,SCALAR_SUFFIX),asum_)(int *n,
 }
 
 // computes a dot product of a conjugated vector with another vector.
-Scalar EIGEN_BLAS_FUNC(dotc)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy)
+int EIGEN_BLAS_FUNC(dotcw)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar* pres)
 {
 //   std::cerr << "_dotc " << *n << " " << *incx << " " << *incy << "\n";
 
@@ -60,18 +60,18 @@ Scalar EIGEN_BLAS_FUNC(dotc)(int *n, RealScalar *px, int *incx, RealScalar *py,
 
   Scalar* x = reinterpret_cast<Scalar*>(px);
   Scalar* y = reinterpret_cast<Scalar*>(py);
+  Scalar* res = reinterpret_cast<Scalar*>(pres);
 
-  Scalar res;
-  if(*incx==1 && *incy==1)    res = (vector(x,*n).dot(vector(y,*n)));
-  else if(*incx>0 && *incy>0) res = (vector(x,*n,*incx).dot(vector(y,*n,*incy)));
-  else if(*incx<0 && *incy>0) res = (vector(x,*n,-*incx).reverse().dot(vector(y,*n,*incy)));
-  else if(*incx>0 && *incy<0) res = (vector(x,*n,*incx).dot(vector(y,*n,-*incy).reverse()));
-  else if(*incx<0 && *incy<0) res = (vector(x,*n,-*incx).reverse().dot(vector(y,*n,-*incy).reverse()));
-  return res;
+  if(*incx==1 && *incy==1)    *res = (vector(x,*n).dot(vector(y,*n)));
+  else if(*incx>0 && *incy>0) *res = (vector(x,*n,*incx).dot(vector(y,*n,*incy)));
+  else if(*incx<0 && *incy>0) *res = (vector(x,*n,-*incx).reverse().dot(vector(y,*n,*incy)));
+  else if(*incx>0 && *incy<0) *res = (vector(x,*n,*incx).dot(vector(y,*n,-*incy).reverse()));
+  else if(*incx<0 && *incy<0) *res = (vector(x,*n,-*incx).reverse().dot(vector(y,*n,-*incy).reverse()));
+  return 0;
 }
 
 // computes a vector-vector dot product without complex conjugation.
-Scalar EIGEN_BLAS_FUNC(dotu)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy)
+int EIGEN_BLAS_FUNC(dotuw)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar* pres)
 {
 //   std::cerr << "_dotu " << *n << " " << *incx << " " << *incy << "\n";
 
@@ -79,13 +79,14 @@ Scalar EIGEN_BLAS_FUNC(dotu)(int *n, RealScalar *px, int *incx, RealScalar *py,
 
   Scalar* x = reinterpret_cast<Scalar*>(px);
   Scalar* y = reinterpret_cast<Scalar*>(py);
-  Scalar res;
-  if(*incx==1 && *incy==1)    res = (vector(x,*n).cwiseProduct(vector(y,*n))).sum();
-  else if(*incx>0 && *incy>0) res = (vector(x,*n,*incx).cwiseProduct(vector(y,*n,*incy))).sum();
-  else if(*incx<0 && *incy>0) res = (vector(x,*n,-*incx).reverse().cwiseProduct(vector(y,*n,*incy))).sum();
-  else if(*incx>0 && *incy<0) res = (vector(x,*n,*incx).cwiseProduct(vector(y,*n,-*incy).reverse())).sum();
-  else if(*incx<0 && *incy<0) res = (vector(x,*n,-*incx).reverse().cwiseProduct(vector(y,*n,-*incy).reverse())).sum();
-  return res;
+  Scalar* res = reinterpret_cast<Scalar*>(pres);
+
+  if(*incx==1 && *incy==1)    *res = (vector(x,*n).cwiseProduct(vector(y,*n))).sum();
+  else if(*incx>0 && *incy>0) *res = (vector(x,*n,*incx).cwiseProduct(vector(y,*n,*incy))).sum();
+  else if(*incx<0 && *incy>0) *res = (vector(x,*n,-*incx).reverse().cwiseProduct(vector(y,*n,*incy))).sum();
+  else if(*incx>0 && *incy<0) *res = (vector(x,*n,*incx).cwiseProduct(vector(y,*n,-*incy).reverse())).sum();
+  else if(*incx<0 && *incy<0) *res = (vector(x,*n,-*incx).reverse().cwiseProduct(vector(y,*n,-*incy).reverse())).sum();
+  return 0;
 }
 
 RealScalar EIGEN_CAT(EIGEN_CAT(REAL_SCALAR_SUFFIX,SCALAR_SUFFIX),nrm2_)(int *n, RealScalar *px, int *incx)