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

This commit is contained in:
Gael Guennebaud 2011-02-18 15:11:31 +01:00
parent 69cecc45e5
commit f7cd63b964
3 changed files with 65 additions and 16 deletions

View File

@ -46,6 +46,11 @@ double BLASFUNC(xdotu) (int *, double *, int *, double *, int *);
double BLASFUNC(xdotc) (int *, double *, int *, double *, int *);
#endif
int BLASFUNC(cdotuw) (int *, float *, int *, float *, int *, float*);
int BLASFUNC(cdotcw) (int *, float *, int *, float *, int *, float*);
int BLASFUNC(zdotuw) (int *, double *, int *, double *, int *, double*);
int BLASFUNC(zdotcw) (int *, double *, int *, double *, int *, double*);
int BLASFUNC(saxpy) (int *, float *, float *, int *, float *, int *);
int BLASFUNC(daxpy) (int *, double *, double *, int *, double *, int *);
int BLASFUNC(qaxpy) (int *, double *, double *, int *, double *, int *);

43
blas/complexdots.f Normal file
View File

@ -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

View File

@ -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)