Skip to content
Snippets Groups Projects
Commit 27d172fc authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

pp

parent 3a51056d
No related branches found
No related tags found
No related merge requests found
......@@ -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,16 +37,11 @@ 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);
}
......@@ -294,9 +290,9 @@ 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);
......
......@@ -22,7 +22,6 @@ class fullVector
scalar *_data;
friend class fullMatrix<scalar>;
public:
fullVector(int r) : _r(r)
{
_data = new scalar[_r];
......@@ -35,11 +34,9 @@ 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;
_data = new scalar[_r];
......@@ -47,7 +44,6 @@ class fullVector
}
return false;
}
inline scalar operator () (int i) const
{
return _data[i];
......@@ -76,7 +72,6 @@ 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)
{
......@@ -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;
......@@ -175,8 +168,7 @@ class fullMatrix
if (_data && _own_data) delete[] _data;
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)
{
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)
{
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment