diff --git a/Numeric/fullMatrix.cpp b/Numeric/fullMatrix.cpp index dc103d8ab27d1ea0a188ec57b6d5dc0beda6cea1..8a75bc796d2c329a96fff7175aaca66491ba03ad 100644 --- a/Numeric/fullMatrix.cpp +++ b/Numeric/fullMatrix.cpp @@ -20,6 +20,7 @@ #if defined(HAVE_BLAS) extern "C" { + void F77NAME(daxpy)(int *n, double *alpha, double *x, int *incx, double *y, int *incy); void F77NAME(dgemm)(const char *transa, const char *transb, int *m, int *n, int *k, double *alpha, double *a, int *lda, double *b, int *ldb, double *beta, @@ -36,22 +37,17 @@ extern "C" { std::complex<double> *alpha, std::complex<double> *a, int *lda, std::complex<double> *x, int *incx, std::complex<double> *beta, std::complex<double> *y, int *incy); - void F77NAME(daxpy)(int *n, double *alpha, double *x, int *incx, double *y, int *incy ); - - } template<> void fullVector<double>::axpy(fullVector<double> &x,double alpha) { - // for (int i=0;i<_r;i++) _data[i] += alpha * x._data[i]; - // return; int M = _r, INCX = 1, INCY = 1; F77NAME(daxpy)(&M, &alpha, x._data,&INCX, _data, &INCY); } template<> -void fullMatrix<double>::mult(const fullMatrix<double> &b, fullMatrix<double> &c)const +void fullMatrix<double>::mult(const fullMatrix<double> &b, fullMatrix<double> &c) const { int M = c.size1(), N = c.size2(), K = _c; int LDA = _r, LDB = b.size1(), LDC = c.size1(); @@ -62,7 +58,7 @@ void fullMatrix<double>::mult(const fullMatrix<double> &b, fullMatrix<double> &c template<> void fullMatrix<std::complex<double> >::mult(const fullMatrix<std::complex<double> > &b, - fullMatrix<std::complex<double> > &c)const + fullMatrix<std::complex<double> > &c) const { int M = c.size1(), N = c.size2(), K = _c; int LDA = _r, LDB = b.size1(), LDC = c.size1(); @@ -294,13 +290,13 @@ bool fullMatrix<double>::svd(fullMatrix<double> &V, fullVector<double> &S) #include "Bindings.h" template<> -void fullMatrix<double>::registerBindings(binding *b){ +void fullMatrix<double>::registerBindings(binding *b) +{ classBinding *cb = b->addClass<fullMatrix<double> >("fullMatrix"); - methodBinding *cm; - cb->addMethod("size1",&fullMatrix<double>::size1); - cb->addMethod("size2",&fullMatrix<double>::size2); - cb->addMethod("get",&fullMatrix<double>::get); - cb->addMethod("set",&fullMatrix<double>::set); - cb->addMethod("gemm",&fullMatrix<double>::gemm); + cb->addMethod("size1", &fullMatrix<double>::size1); + cb->addMethod("size2", &fullMatrix<double>::size2); + cb->addMethod("get", &fullMatrix<double>::get); + cb->addMethod("set", &fullMatrix<double>::set); + cb->addMethod("gemm", &fullMatrix<double>::gemm); cb->setConstructor<fullMatrix<double>,int,int>(); } diff --git a/Numeric/fullMatrix.h b/Numeric/fullMatrix.h index 9bec79e685025b0be6ed4a2fa929b77ca8e33b54..c17881309e0edc37854c0eff365569d9ccd1aaca 100644 --- a/Numeric/fullMatrix.h +++ b/Numeric/fullMatrix.h @@ -22,7 +22,6 @@ class fullVector scalar *_data; friend class fullMatrix<scalar>; public: - fullVector(int r) : _r(r) { _data = new scalar[_r]; @@ -35,19 +34,16 @@ class fullVector for(int i = 0; i < _r; ++i) _data[i] = other._data[i]; } ~fullVector() { if(_data) delete [] _data; } - bool resize(int r) { - if ((_r<r)) - { + if (_r < r){ if (_data) delete[] _data; - _r=r; + _r = r; _data = new scalar[_r]; return true; } return false; } - inline scalar operator () (int i) const { return _data[i]; @@ -76,11 +72,10 @@ class fullVector for(int i = 0; i < _r; ++i) s += _data[i] * other._data[i]; return s; } - // y <- y + alpha * x void axpy(fullVector<scalar> &x, scalar alpha=1.) #if !defined(HAVE_BLAS) { - for (int i=0;i<_r;i++) _data[i] += alpha * x._data[i]; + for (int i = 0; i < _r; i++) _data[i] += alpha * x._data[i]; } #endif ; @@ -103,19 +98,17 @@ class fullMatrix int _r, _c; scalar *_data; public: - inline scalar get (int r, int c)const { - return (*this)(r,c); - } - inline void set(int r, int c, scalar v){ - (*this)(r,c)=v; - } - fullMatrix(scalar *original, int r, int c){ + inline scalar get(int r, int c) const { return (*this)(r, c); } + inline void set(int r, int c, scalar v){ (*this)(r, c) = v; } + fullMatrix(scalar *original, int r, int c) + { _r = r; _c = c; _own_data = false; _data = original; } - fullMatrix(fullMatrix<scalar> &original, int c_start, int c){ + fullMatrix(fullMatrix<scalar> &original, int c_start, int c) + { _c = c; _r = original._r; _own_data = false; @@ -127,7 +120,8 @@ class fullMatrix _own_data = true; scale(0.); } - fullMatrix(int r, int c, double *data) : _r(r), _c(c), _data(data),_own_data(false) + fullMatrix(int r, int c, double *data) + : _r(r), _c(c), _data(data), _own_data(false) { scale(0.); } @@ -139,25 +133,24 @@ class fullMatrix } fullMatrix() : _own_data(false),_r(0), _c(0), _data(0) {} ~fullMatrix() { if(_data && _own_data) delete [] _data; } - bool resize(int r, int c) // data will be owned (same as constructor) { - if ((r*c>_r*_c)||(!_own_data)) - { - _r=r;_c=c; - if ((_own_data)&&(_data)) delete[] _data; + if ((r * c > _r * _c) || !_own_data){ + _r = r; + _c = c; + if (_own_data && _data) delete[] _data; _data = new scalar[_r * _c]; _own_data = true; return true; } - else - { - _r=r;_c=c; + else{ + _r = r; + _c = c; } return false; // no reallocation } - - void setAsProxy(fullMatrix<scalar> &original, int c_start, int c) { + void setAsProxy(fullMatrix<scalar> &original, int c_start, int c) + { if(_data && _own_data) delete [] _data; _c = c; @@ -173,10 +166,9 @@ class fullMatrix _r = other._r; _c = other._c; if (_data && _own_data) delete[] _data; - if ((_r==0)||(_c==0)) + if ((_r == 0) || (_c == 0)) _data=0; - else - { + else{ _data = new scalar[_r * _c]; _own_data=true; for(int i = 0; i < _r * _c; ++i) _data[i] = other._data[i]; @@ -243,14 +235,12 @@ class fullMatrix for(int j = 0; j < size2(); j++) (*this)(i, j) += m(i, j); } - -inline void add(const fullMatrix<scalar> &m, const double &a) + inline void add(const fullMatrix<scalar> &m, const double &a) { for(int i = 0; i < size1(); i++) for(int j = 0; j < size2(); j++) (*this)(i, j) += a*m(i, j); } - void mult(const fullVector<scalar> &x, fullVector<scalar> &y) #if !defined(HAVE_BLAS) {