43 template<
class Treal,
class Telement>
53 template<
class Treal,
class Telement = Treal>
54 class Vector:
public VectorHierarchicBase<Treal, Telement> {
64 this->
elements = allocateElements<Telement>(this->
n());
92 (*this) *= (1.0 / this->
eucl());
96 inline Treal
eucl()
const {
106 static void axpy(Treal
const & alpha,
116 template<
typename TmatrixElement>
117 static void gemv(
bool const tA, Treal
const alpha,
126 template<
typename TmatrixElement>
127 static void symv(
char const uplo, Treal
const alpha,
136 template<
typename TmatrixElement>
137 static void trmv(
char const uplo,
const bool tA,
143 void assign_from_full(Treal
const *
const fullvector,
const int totn);
145 void fullvector(Treal *
const full,
const int totn)
const;
159 template<
class Treal,
class Telement>
162 addFromFull(fullVector);
165 template<
class Treal,
class Telement>
170 for (
int ind = 0; ind < this->n(); ind++)
171 (*
this)(ind).addFromFull(fullVector);
174 template<
class Treal,
class Telement>
177 if (this->is_zero()) {
179 for (
int row = 0; row < this->nScalars(); ++row )
183 for (
int ind = 0; ind < this->n(); ind++)
184 (*
this)(ind).fullVector(fullVec);
188 template<
class Treal,
class Telement>
194 template<
class Treal,
class Telement>
199 if (this->is_zero()) {
200 char * tmp = (
char*)&ZERO;
201 file.write(tmp,
sizeof(
int));
204 char * tmp = (
char*)&ONE;
205 file.write(tmp,
sizeof(
int));
206 for (
int i = 0; i < this->n(); i++)
207 this->elements[i].writeToFile(file);
210 template<
class Treal,
class Telement>
215 char tmp[
sizeof(int)];
216 file.read(tmp, (std::ifstream::pos_type)
sizeof(
int));
224 for (
int i = 0; i < this->n(); i++)
225 this->elements[i].readFromFile(file);
228 throw Failure(
"Vector<Treal, Telement>::" 229 "readFromFile(std::ifstream & file):" 230 "File corruption int value not 0 or 1");
234 template<
class Treal,
class Telement>
240 throw Failure(
"Vector::operator=(int k) only " 241 "implemented for k = 0");
245 template<
class Treal,
class Telement>
249 for (
int ind = 0; ind < this->n(); ind++)
250 (*
this)(ind).random();
255 template<
class Treal,
class Telement>
258 if (!this->is_zero() && alpha != 1) {
259 for (
int ind = 0; ind < this->n(); ind++)
260 (*
this)(ind) *= alpha;
265 template<
class Treal,
class Telement>
269 assert(x.n() == y.n());
270 if (x.is_zero() || y.is_zero())
272 Treal dotProduct = 0;
273 for (
int ind = 0; ind < x.n(); ind++)
279 template<
class Treal,
class Telement>
284 assert(x.n() == y.n());
290 for (
int ind = 0; ind < x.n(); ind++)
300 template<
class Treal,
class Telement>
301 template<
typename TmatrixElement>
303 gemv(
bool const tA, Treal
const alpha,
312 if ((
A.is_zero() || x.is_zero() || alpha == 0) &&
313 (y.is_zero() || beta == 0))
316 Treal beta_tmp = beta;
321 if (
A.is_zero() || x.is_zero() || alpha == 0)
326 if (
A.ncols() != x.n() ||
A.nrows() != y.n())
327 throw Failure(
"Vector<Treal, Telement>::" 328 "gemv(bool const, Treal const, " 329 "const Matrix<Treal, Telement>&, " 330 "const Vector<Treal, Telement>&, " 331 "Treal const, const Vector<Treal, " 333 "Incorrect dimensions for matrix-vector product");
335 int A_nrows =
A.nrows();
337 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic) 339 for (
int row = 0; row < A_nrows; row++) {
342 for (
int col = 1; col <
A.ncols(); col++)
350 if (
A.nrows() != x.n() ||
A.ncols() != y.n())
351 throw Failure(
"Vector<Treal, Telement>::" 352 "gemv(bool const, Treal const, " 353 "const Matrix<Treal, Telement>&, " 354 "const Vector<Treal, Telement>&, " 355 "Treal const, const Vector<Treal, " 357 "Incorrect dimensions for matrix-vector product");
359 int A_ncols =
A.ncols();
361 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic) 363 for (
int col = 0; col < A_ncols; col++) {
366 for (
int row = 1; row <
A.nrows(); row++)
380 template<
class Treal,
class Telement>
381 template<
typename TmatrixElement>
383 symv(
char const uplo, Treal
const alpha,
392 if (x.n() != y.n() ||
A.nrows() !=
A.ncols() ||
A.ncols() != x.n())
393 throw Failure(
"Vector<Treal, Telement>::" 394 "symv(char const uplo, Treal const, " 395 "const Matrix<Treal, Telement>&, " 396 "const Vector<Treal, Telement>&, " 397 "Treal const, const Vector<Treal, Telement>&):" 398 "Incorrect dimensions for symmetric " 399 "matrix-vector product");
401 throw Failure(
"Vector<class Treal, class Telement>::" 402 "symv only implemented for symmetric matrices in " 403 "upper triangular storage");
404 if ((
A.is_zero() || x.is_zero() || alpha == 0) &&
405 (y.is_zero() || beta == 0))
408 Treal beta_tmp = beta;
413 if (
A.is_zero() || x.is_zero() || alpha == 0)
418 #pragma omp parallel if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) 422 int A_ncols =
A.ncols();
424 #pragma omp for schedule(dynamic) 426 for (
int rc = 0; rc < A_ncols; rc++) {
432 int A_nrows =
A.nrows();
434 #pragma omp for schedule(dynamic) 436 for (
int row = 0; row < A_nrows - 1; row++) {
438 for (
int col = row + 1; col <
A.ncols(); col++)
444 #pragma omp for schedule(dynamic) 446 for (
int row = 1; row < A_nrows; row++) {
448 for (
int col = 0; col < row; col++)
458 template<
class Treal,
class Telement>
459 template<
typename TmatrixElement>
461 trmv(
char const uplo,
const bool tA,
464 if (
A.nrows() !=
A.ncols() ||
A.ncols() != x.n())
465 throw Failure(
"Vector<Treal, Telement>::" 467 "Incorrect dimensions for triangular " 468 "matrix-vector product");
470 throw Failure(
"Vector<class Treal, class Telement>::" 471 "trmv only implemented for upper triangular matrices");
472 if ( (
A.is_zero() || x.is_zero() ) ) {
478 for (
int row = 0; row <
A.nrows(); row++) {
480 for (
int col = row + 1; col <
A.ncols(); col++)
486 for (
int col =
A.ncols() - 1; col >= 0; col--) {
488 for (
int row = 0; row < col; row++)
502 template<
class Treal>
512 this->
elements = allocateElements<Treal>(this->
n());
513 for (
int ind = 0; ind < this->
n(); ind++)
540 (*this) *= 1 / this->
eucl();
556 static void axpy(Treal
const & alpha,
565 static void gemv(
bool const tA, Treal
const alpha,
574 static void symv(
char const uplo, Treal
const alpha,
584 static void trmv(
char const uplo,
const bool tA,
591 template<
class Treal>
597 template<
class Treal>
606 for (
int row = 0; row < this->
n(); ++row )
610 template<
class Treal>
615 for (
int row = 0; row < this->
nScalars(); ++row )
618 for (
int row = 0; row < this->
n(); ++row )
623 template<
class Treal>
630 template<
class Treal>
636 char * tmp = (
char*)&ZERO;
637 file.write(tmp,
sizeof(
int));
640 char * tmp = (
char*)&ONE;
641 file.write(tmp,
sizeof(
int));
642 char * tmpel = (
char*)this->
elements;
643 file.write(tmpel,
sizeof(Treal) * this->
n());
647 template<
class Treal>
652 char tmp[
sizeof(int)];
653 file.read(tmp, (std::ifstream::pos_type)
sizeof(
int));
661 file.read((
char*)this->
elements,
sizeof(Treal) * this->
n());
664 throw Failure(
"Vector<Treal>::" 665 "readFromFile(std::ifstream & file):" 666 "File corruption, int value not 0 or 1");
670 template<
class Treal>
676 throw Failure(
"Vector::operator=(int k) only implemented for k = 0");
680 template<
class Treal>
684 for (
int ind = 0; ind < this->
n(); ind++)
685 (*
this)(ind) = rand() / (Treal)RAND_MAX;
689 template<
class Treal>
692 if (!this->
is_zero() && alpha != 1) {
699 template<
class Treal>
703 assert(x.
n() == y.
n());
713 template<
class Treal>
718 assert(x.
n() == y.
n());
723 for (
int ind = 0; ind < y.
n(); ind++)
736 template<
class Treal>
738 gemv(
bool const tA, Treal
const alpha,
747 if ((
A.is_zero() || x.
is_zero() || alpha == 0) &&
751 Treal beta_tmp = beta;
756 if (
A.is_zero() || x.
is_zero() || alpha == 0)
761 if (
A.ncols() != x.
n() ||
A.nrows() != y.
n())
762 throw Failure(
"Vector<Treal, Telement>::" 763 "gemv(bool const, Treal const, " 764 "const Matrix<Treal, Telement>&, " 765 "const Vector<Treal, Telement>&, " 766 "Treal const, const Vector<Treal, " 768 "Incorrect dimensions for matrix-vector product");
770 mat::gemv(
"N", &
A.nrows(), &
A.ncols(), &alpha,
A.elements,
776 if (
A.nrows() != x.
n() ||
A.ncols() != y.
n())
777 throw Failure(
"Vector<Treal, Telement>::" 778 "gemv(bool const, Treal const, " 779 "const Matrix<Treal, Telement>&, " 780 "const Vector<Treal, Telement>&, " 781 "Treal const, const Vector<Treal, " 783 "Incorrect dimensions for matrix-vector product");
785 mat::gemv(
"T", &
A.nrows(), &
A.ncols(), &alpha,
A.elements,
796 template<
class Treal>
798 symv(
char const uplo, Treal
const alpha,
807 if (x.
n() != y.
n() ||
A.nrows() !=
A.ncols() ||
A.ncols() != x.
n())
808 throw Failure(
"Vector<Treal>::" 809 "symv(char const uplo, Treal const, " 810 "const Matrix<Treal>&, " 811 "const Vector<Treal>&, " 812 "Treal const, const Vector<Treal>&):" 813 "Incorrect dimensions for symmetric " 814 "matrix-vector product");
815 if ((
A.is_zero() || x.
is_zero() || alpha == 0) &&
819 Treal beta_tmp = beta;
824 if (
A.is_zero() || x.
is_zero() || alpha == 0)
828 mat::symv(&uplo, &x.
n(), &alpha,
A.elements, &
A.nrows(),
834 template<
class Treal>
836 trmv(
char const uplo,
const bool tA,
839 if (
A.nrows() !=
A.ncols() ||
A.ncols() != x.
n())
840 throw Failure(
"Vector<Treal>::" 841 "trmv(...): Incorrect dimensions for triangular " 842 "matrix-vector product");
844 throw Failure(
"Vector<class Treal>::" 845 "trmv only implemented for upper triangular matrices");
846 if ( (
A.is_zero() || x.
is_zero() ) ) {
852 mat::trmv(&uplo,
"N",
"N", &x.
n(),
A.elements, &
A.nrows(),
855 mat::trmv(&uplo,
"T",
"N", &x.
n(),
A.elements, &
A.nrows(),
int getOffset() const
Definition: SizesAndBlocks.h:83
bool is_zero() const
Definition: VectorHierarchicBase.h:75
Definition: Matrix.h:2925
void writeToFile(std::ofstream &file) const
Definition: Vector.h:196
void random()
Definition: Vector.h:246
mat::SizesAndBlocks rows
Definition: test.cc:51
VectorHierarchicBase< Treal, Telement > & operator=(const VectorHierarchicBase< Treal, Telement > &vec)
Definition: VectorHierarchicBase.h:152
void clear()
Definition: Vector.h:189
#define MAT_OMP_FINALIZE
Definition: matInclude.h:92
int const & getNBlocks() const
Definition: SizesAndBlocks.h:72
SizesAndBlocks getSizesAndBlocksForLowerLevel(int const blockNumber) const
Definition: SizesAndBlocks.cc:71
void addFromFull(std::vector< Treal > const &fullVector)
Definition: Vector.h:167
Vector()
Definition: Vector.h:59
Vector class.
Definition: Matrix.h:78
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: mat_gblas.h:409
void readFromFile(std::ifstream &file)
Definition: Vector.h:212
void randomNormalized()
Definition: Vector.h:538
Definition: allocate.cc:39
Describes dimensions of matrix and its blocks on all levels.
Definition: SizesAndBlocks.h:45
int getNTotalScalars() const
Definition: SizesAndBlocks.h:84
Treal eucl() const
Definition: Vector.h:544
bool is_empty() const
Check if vector is empty Empty is different from zero, a zero matrix contains information about block...
Definition: VectorHierarchicBase.h:89
static Treal dot(Vector< Treal, Telement > const &x, Vector< Treal, Telement > const &y)
Definition: Vector.h:267
void assignFromFull(std::vector< Treal > const &fullVector)
Definition: Vector.h:161
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:461
Vector()
Definition: Vector.h:506
static T dot(const int *n, const T *dx, const int *incx, const T *dy, const int *incy)
Definition: mat_gblas.h:425
Vector< Treal > & operator=(const Vector< Treal > &vec)
Definition: Vector.h:525
static void axpy(const int *n, const T *da, const T *dx, const int *incx, T *dy, const int *incy)
Definition: mat_gblas.h:431
Telement ElementType
Definition: Vector.h:56
Treal eucl() const
Definition: Vector.h:96
void freeElements(float *ptr)
Definition: allocate.cc:49
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: mat_gblas.h:400
#define MAT_OMP_START
Definition: matInclude.h:90
Matrix class and heart of the matrix library.
Definition: Matrix.h:115
Vector< Treal, Telement > & operator=(const Vector< Treal, Telement > &vec)
Definition: Vector.h:79
void randomNormalized()
Definition: Vector.h:90
void resetRows(SizesAndBlocks const &newRows)
Definition: VectorHierarchicBase.h:77
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:303
#define MAT_OMP_END
Definition: matInclude.h:91
static void scal(const int *n, const T *da, T *dx, const int *incx)
Definition: mat_gblas.h:419
Telement * elements
Definition: VectorHierarchicBase.h:108
SizesAndBlocks rows
Definition: VectorHierarchicBase.h:105
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: mat_gblas.h:391
#define MAT_OMP_INIT
Definition: matInclude.h:89
Vector< Treal, Telement > & operator*=(const Treal alpha)
Definition: Vector.h:257
void allocate()
Definition: Vector.h:509
static void axpy(Treal const &alpha, Vector< Treal, Telement > const &x, Vector< Treal, Telement > &y)
Definition: Vector.h:281
const int & nScalars() const
Definition: VectorHierarchicBase.h:55
void fullVector(std::vector< Treal > &fullVector) const
Definition: Vector.h:176
Base class for Vector and Vector specialization.
Definition: VectorHierarchicBase.h:51
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:383
void allocate()
Definition: Vector.h:61
Treal template_blas_sqrt(Treal x)
const int & n() const
Definition: VectorHierarchicBase.h:58