41 template<
class Treal,
class Telement>
51 template<
class Treal,
class Telement = Treal>
52 class Vector:
public VectorHierarchicBase<Treal, Telement> {
62 this->
elements = allocateElements<Telement>(this->
n());
90 (*this) *= (1.0 / this->
eucl());
94 inline Treal
eucl()
const {
104 static void axpy(Treal
const & alpha,
114 template<
typename TmatrixElement>
115 static void gemv(
bool const tA, Treal
const alpha,
124 template<
typename TmatrixElement>
125 static void symv(
char const uplo, Treal
const alpha,
134 template<
typename TmatrixElement>
135 static void trmv(
char const uplo,
const bool tA,
141 void assign_from_full(Treal
const *
const fullvector,
const int totn);
143 void fullvector(Treal *
const full,
const int totn)
const;
157 template<
class Treal,
class Telement>
160 addFromFull(fullVector);
163 template<
class Treal,
class Telement>
168 for (
int ind = 0; ind < this->n(); ind++)
169 (*
this)(ind).addFromFull(fullVector);
172 template<
class Treal,
class Telement>
175 if (this->is_zero()) {
176 fullVec.resize(this->rows.getNTotalScalars());
177 for (
int row = 0; row < this->nScalars(); ++row )
178 fullVec[this->rows.getOffset()+row] = 0;
181 for (
int ind = 0; ind < this->n(); ind++)
182 (*
this)(ind).fullVector(fullVec);
186 template<
class Treal,
class Telement>
192 template<
class Treal,
class Telement>
197 if (this->is_zero()) {
198 char * tmp = (
char*)&ZERO;
199 file.write(tmp,
sizeof(
int));
202 char * tmp = (
char*)&ONE;
203 file.write(tmp,
sizeof(
int));
204 for (
int i = 0; i < this->n(); i++)
205 this->elements[i].writeToFile(file);
208 template<
class Treal,
class Telement>
213 char tmp[
sizeof(int)];
214 file.read(tmp, (std::ifstream::pos_type)
sizeof(
int));
222 for (
int i = 0; i < this->n(); i++)
223 this->elements[i].readFromFile(file);
226 throw Failure(
"Vector<Treal, Telement>::"
227 "readFromFile(std::ifstream & file):"
228 "File corruption int value not 0 or 1");
232 template<
class Treal,
class Telement>
238 throw Failure(
"Vector::operator=(int k) only "
239 "implemented for k = 0");
243 template<
class Treal,
class Telement>
247 for (
int ind = 0; ind < this->n(); ind++)
248 (*
this)(ind).random();
253 template<
class Treal,
class Telement>
256 if (!this->is_zero() && alpha != 1) {
257 for (
int ind = 0; ind < this->n(); ind++)
258 (*
this)(ind) *= alpha;
263 template<
class Treal,
class Telement>
267 assert(x.n() == y.n());
268 if (x.is_zero() || y.is_zero())
270 Treal dotProduct = 0;
271 for (
int ind = 0; ind < x.n(); ind++)
277 template<
class Treal,
class Telement>
282 assert(x.n() == y.n());
288 for (
int ind = 0; ind < x.n(); ind++)
298 template<
class Treal,
class Telement>
299 template<
typename TmatrixElement>
301 gemv(
bool const tA, Treal
const alpha,
310 if ((A.
is_zero() || x.is_zero() || alpha == 0) &&
311 (y.is_zero() || beta == 0))
314 Treal beta_tmp = beta;
319 if (A.
is_zero() || x.is_zero() || alpha == 0)
325 throw Failure(
"Vector<Treal, Telement>::"
326 "gemv(bool const, Treal const, "
327 "const Matrix<Treal, Telement>&, "
328 "const Vector<Treal, Telement>&, "
329 "Treal const, const Vector<Treal, "
331 "Incorrect dimensions for matrix-vector product");
334 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic)
336 for (
int row = 0; row < A.
nrows(); row++) {
339 for (
int col = 1; col < A.
ncols(); col++)
348 throw Failure(
"Vector<Treal, Telement>::"
349 "gemv(bool const, Treal const, "
350 "const Matrix<Treal, Telement>&, "
351 "const Vector<Treal, Telement>&, "
352 "Treal const, const Vector<Treal, "
354 "Incorrect dimensions for matrix-vector product");
357 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic)
359 for (
int col = 0; col < A.
ncols(); col++) {
362 for (
int row = 1; row < A.
nrows(); row++)
376 template<
class Treal,
class Telement>
377 template<
typename TmatrixElement>
379 symv(
char const uplo, Treal
const alpha,
389 throw Failure(
"Vector<Treal, Telement>::"
390 "symv(char const uplo, Treal const, "
391 "const Matrix<Treal, Telement>&, "
392 "const Vector<Treal, Telement>&, "
393 "Treal const, const Vector<Treal, Telement>&):"
394 "Incorrect dimensions for symmetric "
395 "matrix-vector product");
397 throw Failure(
"Vector<class Treal, class Telement>::"
398 "symv only implemented for symmetric matrices in "
399 "upper triangular storage");
400 if ((A.
is_zero() || x.is_zero() || alpha == 0) &&
401 (y.is_zero() || beta == 0))
404 Treal beta_tmp = beta;
409 if (A.
is_zero() || x.is_zero() || alpha == 0)
414 #pragma omp parallel if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared)
419 #pragma omp for schedule(dynamic)
421 for (
int rc = 0; rc < A.
ncols(); rc++) {
428 #pragma omp for schedule(dynamic)
430 for (
int row = 0; row < A.
nrows() - 1; row++) {
432 for (
int col = row + 1; col < A.
ncols(); col++)
438 #pragma omp for schedule(dynamic)
440 for (
int row = 1; row < A.
nrows(); row++) {
442 for (
int col = 0; col < row; col++)
452 template<
class Treal,
class Telement>
453 template<
typename TmatrixElement>
455 trmv(
char const uplo,
const bool tA,
459 throw Failure(
"Vector<Treal, Telement>::"
461 "Incorrect dimensions for triangular "
462 "matrix-vector product");
464 throw Failure(
"Vector<class Treal, class Telement>::"
465 "trmv only implemented for upper triangular matrices");
466 if ( ( A.
is_zero() || x.is_zero() ) ) {
472 for (
int row = 0; row < A.
nrows(); row++) {
474 for (
int col = row + 1; col < A.
ncols(); col++)
480 for (
int col = A.
ncols() - 1; col >= 0; col--) {
482 for (
int row = 0; row < col; row++)
496 template<
class Treal>
506 this->
elements = allocateElements<Treal>(this->
n());
507 for (
int ind = 0; ind < this->
n(); ind++)
534 (*this) *= 1 / this->
eucl();
550 static void axpy(Treal
const & alpha,
559 static void gemv(
bool const tA, Treal
const alpha,
568 static void symv(
char const uplo, Treal
const alpha,
578 static void trmv(
char const uplo,
const bool tA,
585 template<
class Treal>
591 template<
class Treal>
600 for (
int row = 0; row < this->
n(); ++row )
604 template<
class Treal>
609 for (
int row = 0; row < this->
nScalars(); ++row )
612 for (
int row = 0; row < this->
n(); ++row )
617 template<
class Treal>
624 template<
class Treal>
630 char * tmp = (
char*)&ZERO;
631 file.write(tmp,
sizeof(
int));
634 char * tmp = (
char*)&ONE;
635 file.write(tmp,
sizeof(
int));
636 char * tmpel = (
char*)this->
elements;
637 file.write(tmpel,
sizeof(Treal) * this->
n());
641 template<
class Treal>
646 char tmp[
sizeof(int)];
647 file.read(tmp, (std::ifstream::pos_type)
sizeof(
int));
655 file.read((
char*)this->
elements,
sizeof(Treal) * this->
n());
658 throw Failure(
"Vector<Treal>::"
659 "readFromFile(std::ifstream & file):"
660 "File corruption, int value not 0 or 1");
664 template<
class Treal>
670 throw Failure(
"Vector::operator=(int k) only implemented for k = 0");
674 template<
class Treal>
678 for (
int ind = 0; ind < this->
n(); ind++)
679 (*
this)(ind) = rand() / (Treal)RAND_MAX;
683 template<
class Treal>
686 if (!this->
is_zero() && alpha != 1) {
693 template<
class Treal>
697 assert(x.
n() == y.
n());
707 template<
class Treal>
712 assert(x.
n() == y.
n());
717 for (
int ind = 0; ind < y.
n(); ind++)
730 template<
class Treal>
732 gemv(
bool const tA, Treal
const alpha,
745 Treal beta_tmp = beta;
756 throw Failure(
"Vector<Treal, Telement>::"
757 "gemv(bool const, Treal const, "
758 "const Matrix<Treal, Telement>&, "
759 "const Vector<Treal, Telement>&, "
760 "Treal const, const Vector<Treal, "
762 "Incorrect dimensions for matrix-vector product");
771 throw Failure(
"Vector<Treal, Telement>::"
772 "gemv(bool const, Treal const, "
773 "const Matrix<Treal, Telement>&, "
774 "const Vector<Treal, Telement>&, "
775 "Treal const, const Vector<Treal, "
777 "Incorrect dimensions for matrix-vector product");
790 template<
class Treal>
792 symv(
char const uplo, Treal
const alpha,
802 throw Failure(
"Vector<Treal>::"
803 "symv(char const uplo, Treal const, "
804 "const Matrix<Treal>&, "
805 "const Vector<Treal>&, "
806 "Treal const, const Vector<Treal>&):"
807 "Incorrect dimensions for symmetric "
808 "matrix-vector product");
813 Treal beta_tmp = beta;
828 template<
class Treal>
830 trmv(
char const uplo,
const bool tA,
834 throw Failure(
"Vector<Treal>::"
835 "trmv(...): Incorrect dimensions for triangular "
836 "matrix-vector product");
838 throw Failure(
"Vector<class Treal>::"
839 "trmv only implemented for upper triangular matrices");
Definition: Matrix.h:2909
SizesAndBlocks getSizesAndBlocksForLowerLevel(int const blockNumber) const
Definition: SizesAndBlocks.cc:63
void random()
Definition: Vector.h:244
VectorHierarchicBase< Treal, Telement > & operator=(const VectorHierarchicBase< Treal, Telement > &vec)
Definition: VectorHierarchicBase.h:150
void clear()
Definition: Vector.h:187
#define MAT_OMP_FINALIZE
Definition: matInclude.h:88
const int & nScalars() const
Definition: VectorHierarchicBase.h:53
Treal eucl() const
Definition: Vector.h:538
int getOffset() const
Definition: SizesAndBlocks.h:75
int getNTotalScalars() const
Definition: SizesAndBlocks.h:76
void addFromFull(std::vector< Treal > const &fullVector)
Definition: Vector.h:165
Vector()
Definition: Vector.h:57
Vector class.
Definition: Matrix.h:76
return elements[row+col *nrows()]
Definition: MatrixHierarchicBase.h:79
static void trmv(const char *uplo, const char *trans, const char *diag, const int *n, const T *A, const int *lda, T *x, const int *incx)
Definition: gblas.h:407
void readFromFile(std::ifstream &file)
Definition: Vector.h:210
bool is_empty() const
Check if vector is empty Empty is different from zero, a zero matrix contains information about block...
Definition: VectorHierarchicBase.h:87
const int & ncols() const
Definition: MatrixHierarchicBase.h:69
void randomNormalized()
Definition: Vector.h:532
const int & nrows() const
Definition: MatrixHierarchicBase.h:67
Describes dimensions of matrix and its blocks on all levels.
Definition: SizesAndBlocks.h:37
static Treal dot(Vector< Treal, Telement > const &x, Vector< Treal, Telement > const &y)
Definition: Vector.h:265
void assignFromFull(std::vector< Treal > const &fullVector)
Definition: Vector.h:159
static void trmv(char const uplo, const bool tA, Matrix< Treal, TmatrixElement > const &A, Vector< Treal, Telement > &x)
trmv: x = A * x, or x = transpose(A) * x, where A is triangular
Definition: Vector.h:455
Vector()
Definition: Vector.h:500
static T dot(const int *n, const T *dx, const int *incx, const T *dy, const int *incy)
Definition: gblas.h:423
Vector< Treal > & operator=(const Vector< Treal > &vec)
Definition: Vector.h:519
static void axpy(const int *n, const T *da, const T *dx, const int *incx, T *dy, const int *incy)
Definition: gblas.h:429
Telement ElementType
Definition: Vector.h:54
void freeElements(float *ptr)
Definition: allocate.cc:40
bool is_zero() const
Definition: VectorHierarchicBase.h:73
static void symv(const char *uplo, const int *n, const T *alpha, const T *A, const int *lda, const T *x, const int *incx, const T *beta, T *y, const int *incy)
Definition: gblas.h:398
#define MAT_OMP_START
Definition: matInclude.h:86
Matrix class and heart of the matrix library.
Definition: Matrix.h:113
Vector< Treal, Telement > & operator=(const Vector< Treal, Telement > &vec)
Definition: Vector.h:77
void randomNormalized()
Definition: Vector.h:88
void resetRows(SizesAndBlocks const &newRows)
Definition: VectorHierarchicBase.h:75
static void gemv(bool const tA, Treal const alpha, Matrix< Treal, TmatrixElement > const &A, Vector< Treal, Telement > const &x, Treal const beta, Vector< Treal, Telement > &y)
gemv: y = alpha * A * x + beta * y, or y = alpha * transpose(A) * x + beta * y
Definition: Vector.h:301
#define MAT_OMP_END
Definition: matInclude.h:87
Treal eucl() const
Definition: Vector.h:94
static void scal(const int *n, const T *da, T *dx, const int *incx)
Definition: gblas.h:417
Telement * elements
Definition: VectorHierarchicBase.h:106
SizesAndBlocks rows
Definition: VectorHierarchicBase.h:103
static void gemv(const char *ta, const int *m, const int *n, const T *alpha, const T *A, const int *lda, const T *x, const int *incx, const T *beta, T *y, const int *incy)
Definition: gblas.h:389
int const & getNBlocks() const
Definition: SizesAndBlocks.h:64
#define MAT_OMP_INIT
Definition: matInclude.h:85
Vector< Treal, Telement > & operator*=(const Treal alpha)
Definition: Vector.h:255
void allocate()
Definition: Vector.h:503
void writeToFile(std::ofstream &file) const
Definition: Vector.h:194
static void axpy(Treal const &alpha, Vector< Treal, Telement > const &x, Vector< Treal, Telement > &y)
Definition: Vector.h:279
bool is_zero() const
Definition: MatrixHierarchicBase.h:106
const int & n() const
Definition: VectorHierarchicBase.h:56
Base class for Vector and Vector specialization.
Definition: VectorHierarchicBase.h:49
static void symv(char const uplo, Treal const alpha, Matrix< Treal, TmatrixElement > const &A, Vector< Treal, Telement > const &x, Treal const beta, Vector< Treal, Telement > &y)
symv: y = alpha * A * x + beta * y, where A is symmetric
Definition: Vector.h:379
void allocate()
Definition: Vector.h:59
Treal template_blas_sqrt(Treal x)
void fullVector(std::vector< Treal > &fullVector) const
Definition: Vector.h:174