43 template<
class Treal,
class Telement>
53 template<
class Treal,
class Telement = Treal>
92 (*this) *= (1.0 / this->
eucl());
116 template<
typename TmatrixElement>
117 static void gemv(
bool const tA,
Treal const alpha,
126 template<
typename TmatrixElement>
136 template<
typename TmatrixElement>
143 void assign_from_full(
Treal const *
const fullvector,
const int totn);
159 template<
class Treal,
class Telement>
162 addFromFull(fullVector);
165 template<
class Treal,
class Telement>
167 addFromFull(std::vector<Treal>
const & fullVector) {
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++)
188 template<
class Treal,
class Telement>
194 template<
class Treal,
class Telement>
199 if (this->is_zero()) {
200 char *
tmp = (
char*)&ZERO;
204 char *
tmp = (
char*)&ONE;
206 for (
int i = 0; i < this->n(); i++)
207 this->elements[i].writeToFile(
file);
210 template<
class Treal,
class Telement>
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());
279 template<
class Treal,
class Telement>
284 assert(x.
n() == y.
n());
291 Telement::axpy(alpha, x(
ind), y(
ind));
300 template<
class Treal,
class Telement>
301 template<
typename TmatrixElement>
312 if ((
A.is_zero() || x.
is_zero() || alpha == 0) &&
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");
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++) {
341 Telement::gemv(tA, alpha,
A(row, 0), x(0),
beta_tmp, y(row));
342 for (
int col = 1; col <
A.ncols(); col++)
343 Telement::gemv(tA, alpha,
A(row, col), x(col), 1.0, y(row));
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");
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++) {
365 Telement::gemv(tA, alpha,
A(0, col), x(0),
beta_tmp, y(col));
366 for (
int row = 1; row <
A.nrows(); row++)
367 Telement::gemv(tA, alpha,
A(row, col), x(row), 1.0, y(col));
380 template<
class Treal,
class Telement>
381 template<
typename TmatrixElement>
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) &&
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)
424#pragma omp for schedule(dynamic)
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++)
439 Telement::gemv(
false, alpha,
A(row, col), x(col), 1.0, y(row));
444#pragma omp for schedule(dynamic)
446 for (
int row = 1; row <
A_nrows; row++) {
448 for (
int col = 0; col < row; col++)
449 Telement::gemv(
true, alpha,
A(col, row), x(col), 1.0, y(row));
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++) {
479 Telement::trmv(
uplo, tA,
A(row,row), x(row));
480 for (
int col = row + 1; col <
A.ncols(); col++)
481 Telement::gemv(tA, (
Treal)1.0,
A(row, col), x(col), 1.0, x(row));
486 for (
int col =
A.ncols() - 1; col >= 0; col--) {
487 Telement::trmv(
uplo, tA,
A(col,col), x(col));
488 for (
int row = 0; row < col; row++)
489 Telement::gemv(tA, (
Treal)1.0,
A(row, col), x(row), 1.0, x(col));
502 template<
class Treal>
540 (*this) *= 1 / this->
eucl();
565 static void gemv(
bool const tA,
Treal const alpha,
584 static void trmv(
char const uplo,
const bool tA,
591 template<
class Treal>
594 addFromFull(fullVector);
597 template<
class Treal>
599 addFromFull(std::vector<Treal>
const & fullVector) {
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>
635 if (this->is_zero()) {
636 char *
tmp = (
char*)&ZERO;
640 char *
tmp = (
char*)&ONE;
642 char *
tmpel = (
char*)this->elements;
647 template<
class Treal>
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++)
689 template<
class Treal>
692 if (!this->is_zero() && alpha != 1) {
694 mat::scal(&this->n(),&alpha,this->elements,&ONE);
699 template<
class Treal>
703 assert(x.
n() == y.
n());
713 template<
class Treal>
718 assert(x.
n() == y.
n());
736 template<
class Treal>
747 if ((
A.is_zero() || x.
is_zero() || alpha == 0) &&
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>
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) &&
824 if (
A.is_zero() || x.
is_zero() || alpha == 0)
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() ) ) {
Matrix class and heart of the matrix library.
Definition Matrix.h:115
Describes dimensions of matrix and its blocks on all levels.
Definition SizesAndBlocks.h:45
int const & getNBlocks() const
Definition SizesAndBlocks.h:72
int getNTotalScalars() const
Definition SizesAndBlocks.h:84
int getOffset() const
Definition SizesAndBlocks.h:83
SizesAndBlocks getSizesAndBlocksForLowerLevel(int const blockNumber) const
Definition SizesAndBlocks.cc:71
Base class for Vector and Vector specialization.
Definition VectorHierarchicBase.h:51
const int & n() const
Definition VectorHierarchicBase.h:58
SizesAndBlocks rows
Definition VectorHierarchicBase.h:105
void resetRows(SizesAndBlocks const &newRows)
Definition VectorHierarchicBase.h:77
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
Telement * elements
Definition VectorHierarchicBase.h:108
bool is_zero() const
Definition VectorHierarchicBase.h:75
VectorHierarchicBase< Treal, Telement > & operator=(const VectorHierarchicBase< Treal, Telement > &vec)
Definition VectorHierarchicBase.h:152
Vector()
Definition Vector.h:506
void randomNormalized()
Definition Vector.h:538
Treal eucl() const
Definition Vector.h:544
void allocate()
Definition Vector.h:509
Vector< Treal > & operator=(const Vector< Treal > &vec)
Definition Vector.h:525
Vector class.
Definition Vector.h:54
void clear()
Definition Vector.h:189
void readFromFile(std::ifstream &file)
Definition Vector.h:212
static Treal dot(Vector< Treal, Telement > const &x, Vector< Treal, Telement > const &y)
Definition Vector.h:267
Vector< Treal, Telement > & operator*=(const Treal alpha)
Definition Vector.h:257
static void axpy(Treal const &alpha, Vector< Treal, Telement > const &x, Vector< Treal, Telement > &y)
Definition Vector.h:281
Telement ElementType
Definition Vector.h:56
void assignFromFull(std::vector< Treal > const &fullVector)
Definition Vector.h:161
void writeToFile(std::ofstream &file) const
Definition Vector.h:196
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
void allocate()
Definition Vector.h:61
void addFromFull(std::vector< Treal > const &fullVector)
Definition Vector.h:167
void randomNormalized()
Definition Vector.h:90
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
Treal eucl() const
Definition Vector.h:96
void fullVector(std::vector< Treal > &fullVector) const
Definition Vector.h:176
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
Vector()
Definition Vector.h:59
void random()
Definition Vector.h:246
Vector< Treal, Telement > & operator=(const Vector< Treal, Telement > &vec)
Definition Vector.h:79
mat::SizesAndBlocks rows
Definition test.cc:51
#define MAT_OMP_FINALIZE
Definition matInclude.h:92
#define MAT_OMP_INIT
Definition matInclude.h:89
#define MAT_OMP_END
Definition matInclude.h:91
#define MAT_OMP_START
Definition matInclude.h:90
Definition allocate.cc:39
static void scal(const int *n, const T *da, T *dx, const int *incx)
Definition mat_gblas.h:419
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
void freeElements(float *ptr)
Definition allocate.cc:49
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
static Treal getMachineEpsilon()
Definition matInclude.h:147
static T dot(const int *n, const T *dx, const int *incx, const T *dy, const int *incy)
Definition mat_gblas.h:425
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
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
Treal template_blas_sqrt(Treal x)