43 template<
class Treal,
class Telement>
53 template<
class Treal,
class Telement = Treal>
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>
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>
176 fullVector(std::vector<Treal> & fullVec)
const {
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());
272 Treal dotProduct = 0;
273 for (
int ind = 0; ind < x.
n(); ind++)
274 dotProduct += Telement::dot(x(ind), y(ind));
279 template<
class Treal,
class Telement>
281 axpy(Treal
const & alpha,
284 assert(x.
n() == y.
n());
290 for (
int ind = 0; ind < x.
n(); ind++)
291 Telement::axpy(alpha, x(ind), y(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) &&
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++) {
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");
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++) {
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>
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) &&
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++) {
428 Telement::symv(uplo, alpha,
A(rc,rc), x(rc), beta_tmp, y(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++)
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>
505 friend class Matrix<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>
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>
612 fullVector(std::vector<Treal> & fullVec)
const {
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;
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) {
694 mat::scal(&this->n(),&alpha,this->elements,&ONE);
699 template<
class Treal>
703 assert(x.
n() == y.
n());
713 template<
class Treal>
715 axpy(Treal
const & alpha,
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(),
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 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)