GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-08-12 03:39:01 +08:00
Code Issues Packages Projects Releases Wiki Activity

finish/fix level1 blas, all test pass

Browse Source
This commit is contained in:
Gael Guennebaud 2010-07-17 13:49:43 +02:00
parent dd27e10360
commit 2a820d41df
3 changed files with 73 additions and 35 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
bench/btl/libs/C_BLAS/blas.h
View File

@ -38,10 +38,10 @@ void BLASFUNC(zdotc) (double *, int *, double *, int *, double *, int *);
void BLASFUNC(xdotu) (double *, int *, double *, int *, double *, int *); void BLASFUNC(xdotu) (double *, int *, double *, int *, double *, int *);
void BLASFUNC(xdotc) (double *, int *, double *, int *, double *, int *); void BLASFUNC(xdotc) (double *, int *, double *, int *, double *, int *);
#else #else
float BLASFUNC(cdotu) (int *, float *, int *, float *, int *); std::complex<float> BLASFUNC(cdotu) (int *, float *, int *, float *, int *);
float BLASFUNC(cdotc) (int *, float *, int *, float *, int *); std::complex<float> BLASFUNC(cdotc) (int *, float *, int *, float *, int *);
double BLASFUNC(zdotu) (int *, double *, int *, double *, int *); std::complex<double> BLASFUNC(zdotu) (int *, double *, int *, double *, int *);
double BLASFUNC(zdotc) (int *, double *, int *, double *, int *); std::complex<double> BLASFUNC(zdotc) (int *, double *, int *, double *, int *);
double BLASFUNC(xdotu) (int *, double *, int *, double *, int *); double BLASFUNC(xdotu) (int *, double *, int *, double *, int *);
double BLASFUNC(xdotc) (int *, double *, int *, double *, int *); double BLASFUNC(xdotc) (int *, double *, int *, double *, int *);
#endif #endif

1
blas/common.h
View File

@ -26,6 +26,7 @@
#define EIGEN_BLAS_COMMON_H #define EIGEN_BLAS_COMMON_H
#include <iostream> #include <iostream>
#include <complex>
#ifndef SCALAR #ifndef SCALAR
#error the token SCALAR must be defined to compile this file #error the token SCALAR must be defined to compile this file

99
blas/level1_impl.h
View File

@ -153,44 +153,36 @@ Scalar EIGEN_BLAS_FUNC(sdot)(int *n, RealScalar *px, int *incx, RealScalar *py,
#if ISCOMPLEX #if ISCOMPLEX
// computes a dot product of a conjugated vector with another vector. // computes a dot product of a conjugated vector with another vector.
void EIGEN_BLAS_FUNC(dotc)(RealScalar* dot, int *n, RealScalar *px, int *incx, RealScalar *py, int *incy) Scalar EIGEN_BLAS_FUNC(dotc)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy)
{ {
std::cerr << "Eigen BLAS: _dotc is not implemented yet\n";
return;
// TODO: find how to return a complex to fortran
// std::cerr << "_dotc " << *n << " " << *incx << " " << *incy << "\n"; // std::cerr << "_dotc " << *n << " " << *incx << " " << *incy << "\n";
Scalar* x = reinterpret_cast<Scalar*>(px); Scalar* x = reinterpret_cast<Scalar*>(px);
Scalar* y = reinterpret_cast<Scalar*>(py); Scalar* y = reinterpret_cast<Scalar*>(py);
if(*incx==1 && *incy==1) Scalar res;
*reinterpret_cast<Scalar*>(dot) = vector(x,*n).dot(vector(y,*n)); if(*incx==1 && *incy==1) res = (vector(x,*n).dot(vector(y,*n)));
else else if(*incx>0 && *incy>0) res = (vector(x,*n,*incx).dot(vector(y,*n,*incy)));
*reinterpret_cast<Scalar*>(dot) = 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;
} }
// computes a vector-vector dot product without complex conjugation. // computes a vector-vector dot product without complex conjugation.
void EIGEN_BLAS_FUNC(dotu)(RealScalar* dot, int *n, RealScalar *px, int *incx, RealScalar *py, int *incy) Scalar EIGEN_BLAS_FUNC(dotu)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy)
{ {
std::cerr << "Eigen BLAS: _dotu is not implemented yet\n";
return;
// TODO: find how to return a complex to fortran
// std::cerr << "_dotu " << *n << " " << *incx << " " << *incy << "\n"; // std::cerr << "_dotu " << *n << " " << *incx << " " << *incy << "\n";
Scalar* x = reinterpret_cast<Scalar*>(px); Scalar* x = reinterpret_cast<Scalar*>(px);
Scalar* y = reinterpret_cast<Scalar*>(py); Scalar* y = reinterpret_cast<Scalar*>(py);
Scalar res;
if(*incx==1 && *incy==1) if(*incx==1 && *incy==1) res = (vector(x,*n).cwiseProduct(vector(y,*n))).sum();
*reinterpret_cast<Scalar*>(dot) = (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 else if(*incx<0 && *incy>0) res = (vector(x,*n,-*incx).reverse().cwiseProduct(vector(y,*n,*incy))).sum();
*reinterpret_cast<Scalar*>(dot) = (vector(x,*n,*incx).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;
} }
#endif // ISCOMPLEX #endif // ISCOMPLEX
@ -251,15 +243,60 @@ int EIGEN_BLAS_FUNC(rot)(int *n, RealScalar *px, int *incx, RealScalar *py, int
int EIGEN_BLAS_FUNC(rotg)(RealScalar *pa, RealScalar *pb, RealScalar *pc, RealScalar *ps) int EIGEN_BLAS_FUNC(rotg)(RealScalar *pa, RealScalar *pb, RealScalar *pc, RealScalar *ps)
{ {
Scalar a = *reinterpret_cast<Scalar*>(pa); Scalar& a = *reinterpret_cast<Scalar*>(pa);
Scalar b = *reinterpret_cast<Scalar*>(pb); Scalar& b = *reinterpret_cast<Scalar*>(pb);
Scalar* c = reinterpret_cast<Scalar*>(pc); RealScalar* c = pc;
Scalar* s = reinterpret_cast<Scalar*>(ps); Scalar* s = reinterpret_cast<Scalar*>(ps);
PlanarRotation<Scalar> r; #if !ISCOMPLEX
r.makeGivens(a,b); Scalar r,z;
*c = r.c(); Scalar aa = ei_abs(a);
*s = r.s(); Scalar ab = ei_abs(b);
if((aa+ab)==Scalar(0))
{
*c = 1;
*s = 0;
r = 0;
z = 0;
}
else
{
r = ei_sqrt(a*a + b*b);
Scalar amax = aa>ab ? a : b;
r = amax>0 ? r : -r;
*c = a/r;
*s = b/r;
z = 1;
if (aa > ab) z = *s;
if (ab > aa && *c!=RealScalar(0))
z = Scalar(1)/ *c;
}
*pa = r;
*pb = z;
#else
Scalar alpha;
RealScalar norm,scale;
if(ei_abs(a)==RealScalar(0))
{
*c = RealScalar(0);
*s = Scalar(1);
a = b;
}
else
{
scale = ei_abs(a) + ei_abs(b);
norm = scale*ei_sqrt((ei_abs2(a/scale))+ (ei_abs2(b/scale)));
alpha = a/ei_abs(a);
*c = ei_abs(a)/norm;
*s = alpha*ei_conj(b)/norm;
a = alpha*norm;
}
#endif
// PlanarRotation<Scalar> r;
// r.makeGivens(a,b);
// *c = r.c();
// *s = r.s();
return 0; return 0;
} }