Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F81034015
Vector.h
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Tue, Sep 3, 21:37
Size
6 KB
Mime Type
text/x-c++
Expires
Thu, Sep 5, 21:37 (2 d)
Engine
blob
Format
Raw Data
Handle
20504839
Attached To
rLAMMPS lammps
Vector.h
View Options
#ifndef VECTOR_H
#define VECTOR_H
#include "Matrix.h"
///////////////////////////////////////////////////////////////////////////////
// forward declarations ///////////////////////////////////////////////////////
//* Matrix-vector product
//template<typename T>
//void MultMv(const Matrix<T> &A, const Vector<T> &v, DenseVector<T> &c,
// const bool At=0, T a=1, T b=0);
/******************************************************************************
* abstract class Vector
******************************************************************************/
template<typename T>
class Vector : public Matrix<T>
{
public:
Vector() {}
Vector(const Vector<T> &c); // do not implement!
virtual ~Vector() {}
string tostring() const;
// pure virtual functions
virtual T operator()(INDEX i, INDEX j=0) const=0;
virtual T& operator()(INDEX i, INDEX j=0) =0;
virtual T operator[](INDEX i) const=0;
virtual T& operator[](INDEX i) =0;
virtual INDEX nRows() const=0;
virtual T* get_ptr() const=0;
virtual void resize(INDEX nRows, INDEX nCols=1, bool copy=0)=0;
virtual void reset(INDEX nRows, INDEX nCols=1, bool zero=0)=0;
virtual void copy(const T * ptr, INDEX nRows, INDEX nCols=1)=0;
void write_restart(FILE *f) const; // will be virtual
// output to matlab
using Matrix<T>::matlab;
void matlab(ostream &o, const string &s="v") const;
INDEX nCols() const;
bool in_range(INDEX i) const;
bool same_size(const Vector &m) const;
static bool same_size(const Vector &a, const Vector &b);
protected:
//* don't allow this
Vector& operator=(const Vector &r);
};
///////////////////////////////////////////////////////////////////////////////
//* performs a matrix-vector multiply with default naive implementation
template<typename T>
void MultMv(const Matrix<T> &A, const Vector<T> &v, DenseVector<T> &c,
const bool At, T a, T b)
{
const INDEX sA[2] = {A.nRows(), A.nCols()}; // m is sA[At] k is sA[!At]
const INDEX M=sA[At], K=sA[!At];
GCK(A, v, v.size()!=K, "MultAb<T>: matrix-vector multiply");
if (c.size() != M)
{
c.resize(M); // set size of C
c.zero(); // do not add result to C
}
else c *= b;
for (INDEX p=0; p<M; p++)
for (INDEX r=0; r<K; r++)
c[p] += A(p*!At+r*At, p*At+r*!At) * v[r];
}
///////////////////////////////////////////////////////////////////////////////
//* Operator for Matrix-vector product
template<typename T>
DenseVector<T> operator*(const Matrix<T> &A, const Vector<T> &b)
{
DenseVector<T> c;
MultMv(A, b, c, 0, 1.0, 0.0);
return c;
}
///////////////////////////////////////////////////////////////////////////////
//* Operator for Vector-matrix product
template<typename T>
DenseVector<T> operator*(const Vector<T> &a, const Matrix<T> &B)
{
DenseVector<T> c;
MultMv(B, a, c, 1, 1.0, 0.0);
return c;
}
///////////////////////////////////////////////////////////////////////////////
//* Multiply a vector by a scalar
template<typename T>
DenseVector<T> operator*(const Vector<T> &v, const T s)
{
DenseVector<T> r(v);
r*=s;
return r;
}
///////////////////////////////////////////////////////////////////////////////
//* Multiply a vector by a scalar - communitive
template<typename T>
DenseVector<T> operator*(const T s, const Vector<T> &v)
{
DenseVector<T> r(v);
r*=s;
return r;
}
///////////////////////////////////////////////////////////////////////////////
//* inverse scaling operator - must always create memory
template<typename T>
DenseMatrix<T> operator/(const Vector<T> &v, const T s)
{
DenseVector<T> r(v);
r*=(1.0/s); // for integer types this may be worthless
return ;
}
///////////////////////////////////////////////////////////////////////////////
//* Operator for Vector-Vector sum
template<typename T>
DenseVector<T> operator+(const Vector<T> &a, const Vector<T> &b)
{
DenseVector<T> c(a);
c+=b;
return c;
}
///////////////////////////////////////////////////////////////////////////////
//* Operator for Vector-Vector subtraction
template<typename T>
DenseVector<T> operator-(const Vector<T> &a, const Vector<T> &b)
{
DenseVector<T> c(a);
c-=b;
return c;
}
///////////////////////////////////////////////////////////////////////////////
// Template definitions ///////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
//* output operator
template<typename T>
string Vector<T>::tostring() const
{
string s;
FORi s += string(i?"\t":"") + ATC_STRING::tostring(VIDX(i),5);
return s;
}
///////////////////////////////////////////////////////////////////////////////
//* Writes a matlab script defining the vector to the stream
template<typename T>
void Vector<T>::matlab(ostream &o, const string &s) const
{
o << s <<"=zeros(" << this->size() << ",1);\n";
FORi o << s << "("<<i+1<<") = " << VIDX(i) << ";\n";
}
///////////////////////////////////////////////////////////////////////////////
//* writes the vector data to a file
template <typename T>
void Vector<T>::write_restart(FILE *f) const
{
INDEX size = this->size();
fwrite(&size, sizeof(INDEX),1,f);
if (size) fwrite(this->get_ptr(), sizeof(T), this->size(), f);
}
///////////////////////////////////////////////////////////////////////////////
//* returns the number of columns; always 1
template<typename T>
inline INDEX Vector<T>::nCols() const
{
return 1;
}
///////////////////////////////////////////////////////////////////////////////
//* returns true if INDEX i is within the range of the vector
template<typename T>
bool Vector<T>::in_range(INDEX i) const
{
return i<this->size();
}
///////////////////////////////////////////////////////////////////////////////
//* returns true if m has the same number of elements this vector
template<typename T>
bool Vector<T>::same_size(const Vector &m) const
{
return this->size() == m.size();
}
///////////////////////////////////////////////////////////////////////////////
//* returns true if a and b have the same number of elements
template<typename T>
inline bool Vector<T>::same_size(const Vector &a, const Vector &b)
{
return a.same_size(b);
}
#endif
Event Timeline
Log In to Comment