ergo
Vector.h
Go to the documentation of this file.
1 /* Ergo, version 3.7, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2018 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek,
4  * and Anastasia Kruchinina.
5  *
6  * This program is free software: you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation, either version 3 of the License, or
9  * (at your option) any later version.
10  *
11  * This program is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with this program. If not, see <http://www.gnu.org/licenses/>.
18  *
19  * Primary academic reference:
20  * Ergo: An open-source program for linear-scaling electronic structure
21  * calculations,
22  * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia
23  * Kruchinina,
24  * SoftwareX 7, 107 (2018),
25  * <http://dx.doi.org/10.1016/j.softx.2018.03.005>
26  *
27  * For further information about Ergo, see <http://www.ergoscf.org>.
28  */
29 
39 #ifndef MAT_VECTOR
40 #define MAT_VECTOR
41 #include "VectorHierarchicBase.h"
42 namespace mat{
43  template<class Treal, class Telement>
44  class Matrix;
53  template<class Treal, class Telement = Treal>
54  class Vector: public VectorHierarchicBase<Treal, Telement> {
55  public:
56  typedef Telement ElementType;
57  // template<typename TmatrixElement>
58  // friend class Matrix<Treal, TmatrixElement>;
59  Vector():VectorHierarchicBase<Treal, Telement>(){}
60 
61  void allocate() {
62  assert(!this->is_empty());
63  assert(this->is_zero());
64  this->elements = allocateElements<Telement>(this->n());
65  SizesAndBlocks rowSAB;
66  for (int row = 0; row < this->rows.getNBlocks(); row++) {
67  rowSAB = this->rows.getSizesAndBlocksForLowerLevel(row);
68  (*this)(row).resetRows(rowSAB);
69  }
70  }
71 
72  void assignFromFull(std::vector<Treal> const & fullVector);
73 
74  void addFromFull(std::vector<Treal> const & fullVector);
75 
76  void fullVector(std::vector<Treal> & fullVector) const;
77 
81  return *this;
82  }
83 
84  void clear();
85 
86  void writeToFile(std::ofstream & file) const;
87  void readFromFile(std::ifstream & file);
88  Vector<Treal, Telement>& operator=(int const k);
89 
90  inline void randomNormalized() {
91  this->random();
92  (*this) *= (1.0 / this->eucl());
93  }
94  void random();
95 
96  inline Treal eucl() const {
97  return template_blas_sqrt(dot(*this,*this));
98  }
99 
100  /* LEVEL 1 operations */
101  Vector<Treal, Telement>& operator*=(const Treal alpha);
102  static Treal dot(Vector<Treal, Telement> const & x,
103  Vector<Treal, Telement> const & y);
104 
105  /* y += alpha * x */
106  static void axpy(Treal const & alpha,
107  Vector<Treal, Telement> const & x,
109 
110 
111  /* LEVEL 2 operations */
116  template<typename TmatrixElement>
117  static void gemv(bool const tA, Treal const alpha,
119  Vector<Treal, Telement> const & x,
120  Treal const beta,
122 
126  template<typename TmatrixElement>
127  static void symv(char const uplo, Treal const alpha,
129  Vector<Treal, Telement> const & x,
130  Treal const beta,
136  template<typename TmatrixElement>
137  static void trmv(char const uplo, const bool tA,
140 
141 
142 #if 0 /* OLD ROUTINES */
143  void assign_from_full(Treal const * const fullvector, const int totn);
144  /* Convert to full vector */
145  void fullvector(Treal * const full, const int totn) const;
146 
147 
148 
149 
150 
151 
152 
153 
154 
155 #endif /* END OLD ROUTINES */
156  }; /* end class Vector */
157 
158 
159  template<class Treal, class Telement>
161  assignFromFull(std::vector<Treal> const & fullVector) {
162  addFromFull(fullVector);
163  }
164 
165  template<class Treal, class Telement>
167  addFromFull(std::vector<Treal> const & fullVector) {
168  if (this->is_zero())
169  allocate();
170  for (int ind = 0; ind < this->n(); ind++)
171  (*this)(ind).addFromFull(fullVector);
172  }
173 
174  template<class Treal, class Telement>
176  fullVector(std::vector<Treal> & fullVec) const {
177  if (this->is_zero()) {
178  fullVec.resize(this->rows.getNTotalScalars());
179  for (int row = 0; row < this->nScalars(); ++row )
180  fullVec[this->rows.getOffset()+row] = 0;
181  }
182  else
183  for (int ind = 0; ind < this->n(); ind++)
184  (*this)(ind).fullVector(fullVec);
185  }
186 
187 
188  template<class Treal, class Telement>
190  freeElements(this->elements);
191  this->elements = 0;
192  }
193 
194  template<class Treal, class Telement>
196  writeToFile(std::ofstream & file) const {
197  int const ZERO = 0;
198  int const ONE = 1;
199  if (this->is_zero()) {
200  char * tmp = (char*)&ZERO;
201  file.write(tmp,sizeof(int));
202  }
203  else {
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);
208  }
209  }
210  template<class Treal, class Telement>
212  readFromFile(std::ifstream & file) {
213  int const ZERO = 0;
214  int const ONE = 1;
215  char tmp[sizeof(int)];
216  file.read(tmp, (std::ifstream::pos_type)sizeof(int));
217  switch ((int)*tmp) {
218  case ZERO:
219  (*this) = 0;
220  break;
221  case ONE:
222  if (this->is_zero())
223  allocate();
224  for (int i = 0; i < this->n(); i++)
225  this->elements[i].readFromFile(file);
226  break;
227  default:
228  throw Failure("Vector<Treal, Telement>::"
229  "readFromFile(std::ifstream & file):"
230  "File corruption int value not 0 or 1");
231  }
232  }
233 
234  template<class Treal, class Telement>
236  operator=(int const k) {
237  if (k == 0)
238  this->clear();
239  else
240  throw Failure("Vector::operator=(int k) only "
241  "implemented for k = 0");
242  return *this;
243  }
244 
245  template<class Treal, class Telement>
247  if (this->is_zero())
248  allocate();
249  for (int ind = 0; ind < this->n(); ind++)
250  (*this)(ind).random();
251  }
252 
253  /* LEVEL 1 operations */
254 
255  template<class Treal, class Telement>
257  operator*=(const Treal alpha) {
258  if (!this->is_zero() && alpha != 1) {
259  for (int ind = 0; ind < this->n(); ind++)
260  (*this)(ind) *= alpha;
261  }
262  return *this;
263  }
264 
265  template<class Treal, class Telement>
268  Vector<Treal, Telement> const & y) {
269  assert(x.n() == y.n());
270  if (x.is_zero() || y.is_zero())
271  return 0;
272  Treal dotProduct = 0;
273  for (int ind = 0; ind < x.n(); ind++)
274  dotProduct += Telement::dot(x(ind), y(ind));
275  return dotProduct;
276  }
277 
278  /* y += alpha * x */
279  template<class Treal, class Telement>
281  axpy(Treal const & alpha,
282  Vector<Treal, Telement> const & x,
284  assert(x.n() == y.n());
285  if (x.is_zero())
286  return;
287  if (y.is_zero()) {
288  y.allocate();
289  }
290  for (int ind = 0; ind < x.n(); ind++)
291  Telement::axpy(alpha, x(ind), y(ind));
292  }
293 
294  /* LEVEL 2 operations */
295 
300  template<class Treal, class Telement>
301  template<typename TmatrixElement>
303  gemv(bool const tA, Treal const alpha,
305  Vector<Treal, Telement> const & x,
306  Treal const beta,
308  if (y.is_empty()) {
309  assert(beta == 0);
310  y.resetRows(x.rows);
311  }
312  if ((A.is_zero() || x.is_zero() || alpha == 0) &&
313  (y.is_zero() || beta == 0))
314  y = 0;
315  else {
316  Treal beta_tmp = beta;
317  if (y.is_zero()) {
318  y.allocate();
319  beta_tmp = 0;
320  }
321  if (A.is_zero() || x.is_zero() || alpha == 0)
322  y *= beta_tmp;
323  else {
324  MAT_OMP_INIT;
325  if (!tA) {
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, "
332  "Telement>&): "
333  "Incorrect dimensions for matrix-vector product");
334  else {
335  int A_nrows = A.nrows();
336 #ifdef _OPENMP
337 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic)
338 #endif
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));
344  MAT_OMP_END;
345  }
346  } /* end else */
347  } /* end if (!tA) */
348  else {
349  assert(tA);
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, "
356  "Telement>&): "
357  "Incorrect dimensions for matrix-vector product");
358  else {
359  int A_ncols = A.ncols();
360 #ifdef _OPENMP
361 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic)
362 #endif
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));
368  MAT_OMP_END;
369  }
370  } /* end else */
371  } /* end else */
373  } /* end else */
374  } /* end else */
375  }
376 
380  template<class Treal, class Telement>
381  template<typename TmatrixElement>
383  symv(char const uplo, Treal const alpha,
385  Vector<Treal, Telement> const & x,
386  Treal const beta,
388  if (y.is_empty()) {
389  assert(beta == 0);
390  y.resetRows(x.rows);
391  }
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");
400  if (uplo != 'U')
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))
406  y = 0;
407  else {
408  Treal beta_tmp = beta;
409  if (y.is_zero()) {
410  y.allocate();
411  beta_tmp = 0;
412  }
413  if (A.is_zero() || x.is_zero() || alpha == 0)
414  y *= beta_tmp;
415  else {
416  MAT_OMP_INIT;
417 #ifdef _OPENMP
418 #pragma omp parallel if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared)
419 #endif
420  {
421  /* Diagonal */
422  int A_ncols = A.ncols();
423 #ifdef _OPENMP
424 #pragma omp for schedule(dynamic)
425 #endif
426  for (int rc = 0; rc < A_ncols; rc++) {
428  Telement::symv(uplo, alpha, A(rc,rc), x(rc), beta_tmp, y(rc));
429  MAT_OMP_END;
430  }
431  /* Upper triangle */
432  int A_nrows = A.nrows();
433 #ifdef _OPENMP
434 #pragma omp for schedule(dynamic)
435 #endif
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));
440  MAT_OMP_END;
441  }
442  /* Lower triangle */
443 #ifdef _OPENMP
444 #pragma omp for schedule(dynamic)
445 #endif
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));
450  MAT_OMP_END;
451  }
452  } /* end omp parallel*/
454  } /* end else */
455  } /* end else */
456  }
457 
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>::"
466  "trmv(...):"
467  "Incorrect dimensions for triangular "
468  "matrix-vector product");
469  if (uplo != 'U')
470  throw Failure("Vector<class Treal, class Telement>::"
471  "trmv only implemented for upper triangular matrices");
472  if ( ( A.is_zero() || x.is_zero() ) ) {
473  x = 0;
474  return;
475  }
476  if (!tA) {
477  // not transposed
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));
482  }
483  return;
484  }
485  // transposed
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));
490  }
491  }
492 
493 
494 
495 
496  /***************************************************************************/
497  /***************************************************************************/
498  /* Specialization for Telement = Treal */
499  /***************************************************************************/
500  /***************************************************************************/
501 
502  template<class Treal>
503  class Vector<Treal>: public VectorHierarchicBase<Treal> {
504  public:
505  friend class Matrix<Treal>;
507  :VectorHierarchicBase<Treal>(){}
508 
509  void allocate() {
510  assert(!this->is_empty());
511  assert(this->is_zero());
512  this->elements = allocateElements<Treal>(this->n());
513  for (int ind = 0; ind < this->n(); ind++)
514  this->elements[ind] = 0;
515  }
516 
517  void assignFromFull(std::vector<Treal> const & fullVector);
518 
519  void addFromFull(std::vector<Treal> const & fullVector);
520 
521  void fullVector(std::vector<Treal> & fullVector) const;
522 
523 
524  Vector<Treal>&
525  operator=(const Vector<Treal>& vec) {
527  return *this;
528  }
529 
530  void clear();
532  void writeToFile(std::ofstream & file) const;
533  void readFromFile(std::ifstream & file);
534 
535  Vector<Treal>& operator=(int const k);
536 
537 
538  inline void randomNormalized() {
539  this->random();
540  (*this) *= 1 / this->eucl();
541  }
542  void random();
543 
544  inline Treal eucl() const {
545  return template_blas_sqrt(dot(*this,*this));
546  }
547 
548  /* LEVEL 1 operations */
549  Vector<Treal>& operator*=(const Treal alpha);
550 
551  static Treal dot(Vector<Treal> const & x,
552  Vector<Treal> const & y);
553 
554 
555  /* y += alpha * x */
556  static void axpy(Treal const & alpha,
557  Vector<Treal> const & x,
558  Vector<Treal> & y);
559 
560  /* LEVEL 2 operations */
565  static void gemv(bool const tA, Treal const alpha,
566  Matrix<Treal> const & A,
567  Vector<Treal> const & x,
568  Treal const beta,
569  Vector<Treal>& y);
570 
574  static void symv(char const uplo, Treal const alpha,
575  Matrix<Treal> const & A,
576  Vector<Treal> const & x,
577  Treal const beta,
578  Vector<Treal>& y);
579 
584  static void trmv(char const uplo, const bool tA,
585  Matrix<Treal> const & A,
586  Vector<Treal> & x);
587 
588  }; /* end class Vector specialization */
589 
590 
591  template<class Treal>
592  void Vector<Treal>::
593  assignFromFull(std::vector<Treal> const & fullVector) {
595  }
596 
597  template<class Treal>
598  void Vector<Treal>::
599  addFromFull(std::vector<Treal> const & fullVector) {
600  if (this->is_zero())
601  allocate();
602  assert((unsigned)this->rows.getNTotalScalars() == fullVector.size());
603  /* Assertion AFTER empty check done
604  * by allocate()
605  */
606  for (int row = 0; row < this->n(); ++row )
607  (*this)(row) += fullVector[this->rows.getOffset()+row];
608  }
609 
610  template<class Treal>
611  void Vector<Treal>::
612  fullVector(std::vector<Treal> & fullVec) const {
613  fullVec.resize(this->rows.getNTotalScalars());
614  if (this->is_zero())
615  for (int row = 0; row < this->nScalars(); ++row )
616  fullVec[this->rows.getOffset()+row] = 0;
617  else
618  for (int row = 0; row < this->n(); ++row )
619  fullVec[this->rows.getOffset()+row] = (*this)(row);
620  }
621 
622 
623  template<class Treal>
625  freeElements(this->elements);
626  this->elements = 0;
627  }
628 
629 
630  template<class Treal>
631  void Vector<Treal>::
632  writeToFile(std::ofstream & file) const {
633  int const ZERO = 0;
634  int const ONE = 1;
635  if (this->is_zero()) {
636  char * tmp = (char*)&ZERO;
637  file.write(tmp,sizeof(int));
638  }
639  else {
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());
644  }
645  }
646 
647  template<class Treal>
648  void Vector<Treal>::
649  readFromFile(std::ifstream & file) {
650  int const ZERO = 0;
651  int const ONE = 1;
652  char tmp[sizeof(int)];
653  file.read(tmp, (std::ifstream::pos_type)sizeof(int));
654  switch ((int)*tmp) {
655  case ZERO:
656  (*this) = 0;
657  break;
658  case ONE:
659  if (this->is_zero())
660  allocate();
661  file.read((char*)this->elements, sizeof(Treal) * this->n());
662  break;
663  default:
664  throw Failure("Vector<Treal>::"
665  "readFromFile(std::ifstream & file):"
666  "File corruption, int value not 0 or 1");
667  }
668  }
669 
670  template<class Treal>
672  operator=(int const k) {
673  if (k == 0)
674  this->clear();
675  else
676  throw Failure("Vector::operator=(int k) only implemented for k = 0");
677  return *this;
678  }
679 
680  template<class Treal>
682  if (this->is_zero())
683  allocate();
684  for (int ind = 0; ind < this->n(); ind++)
685  (*this)(ind) = rand() / (Treal)RAND_MAX;
686  }
687 
688  /* LEVEL 1 operations */
689  template<class Treal>
691  operator*=(const Treal alpha) {
692  if (!this->is_zero() && alpha != 1) {
693  int const ONE = 1;
694  mat::scal(&this->n(),&alpha,this->elements,&ONE);
695  }
696  return *this;
697  }
698 
699  template<class Treal>
700  Treal Vector<Treal>::
701  dot(Vector<Treal> const & x,
702  Vector<Treal> const & y) {
703  assert(x.n() == y.n());
704  if (x.is_zero() || y.is_zero())
705  return 0;
706  else {
707  int const ONE = 1;
708  return mat::dot(&x.n(), x.elements, &ONE, y.elements, &ONE);
709  }
710  }
711 
712  /* y += alpha * x */
713  template<class Treal>
714  void Vector<Treal>::
715  axpy(Treal const & alpha,
716  Vector<Treal> const & x,
717  Vector<Treal> & y) {
718  assert(x.n() == y.n());
719  if (x.is_zero())
720  return;
721  if (y.is_zero()) {
722  y.allocate();
723  for (int ind = 0; ind < y.n(); ind++)
724  y.elements[ind] = 0; /* fill with zeros */
725  }
726  int const ONE = 1;
727  mat::axpy(&x.n(), &alpha, x.elements, &ONE, y.elements, &ONE);
728  }
729 
730 
731  /* LEVEL 2 operations */
736  template<class Treal>
737  void Vector<Treal>::
738  gemv(bool const tA, Treal const alpha,
739  Matrix<Treal> const & A,
740  Vector<Treal> const & x,
741  Treal const beta,
742  Vector<Treal>& y) {
743  if (y.is_empty()) {
744  assert(beta == 0);
745  y.resetRows(x.rows);
746  }
747  if ((A.is_zero() || x.is_zero() || alpha == 0) &&
748  (y.is_zero() || beta == 0))
749  y = 0;
750  else {
751  Treal beta_tmp = beta;
752  if (y.is_zero()) {
753  y.allocate();
754  beta_tmp = 0;
755  }
756  if (A.is_zero() || x.is_zero() || alpha == 0)
757  y *= beta_tmp;
758  else {
759  int const ONE = 1;
760  if (!tA) {
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, "
767  "Telement>&): "
768  "Incorrect dimensions for matrix-vector product");
769  else {
770  mat::gemv("N", &A.nrows(), &A.ncols(), &alpha, A.elements,
771  &A.nrows(),x.elements,&ONE,&beta_tmp,y.elements,&ONE);
772  } /* end else */
773  } /* end if (!tA) */
774  else {
775  assert(tA);
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, "
782  "Telement>&): "
783  "Incorrect dimensions for matrix-vector product");
784  else {
785  mat::gemv("T", &A.nrows(), &A.ncols(), &alpha, A.elements,
786  &A.nrows(),x.elements,&ONE,&beta_tmp,y.elements,&ONE);
787  } /* end else */
788  } /* end else */
789  } /* end else */
790  } /* end else */
791  }
792 
796  template<class Treal>
797  void Vector<Treal>::
798  symv(char const uplo, Treal const alpha,
799  Matrix<Treal> const & A,
800  Vector<Treal> const & x,
801  Treal const beta,
802  Vector<Treal>& y) {
803  if (y.is_empty()) {
804  assert(beta == 0);
805  y.resetRows(x.rows);
806  }
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) &&
816  (y.is_zero() || beta == 0))
817  y = 0;
818  else {
819  Treal beta_tmp = beta;
820  if (y.is_zero()) {
821  y.allocate();
822  beta_tmp = 0;
823  }
824  if (A.is_zero() || x.is_zero() || alpha == 0)
825  y *= beta_tmp;
826  else {
827  int const ONE = 1;
828  mat::symv(&uplo, &x.n(), &alpha, A.elements, &A.nrows(),
829  x.elements, &ONE, &beta, y.elements, &ONE);
830  } /* end else */
831  } /* end else */
832  }
833 
834  template<class Treal>
835  void Vector<Treal>::
836  trmv(char const uplo, const bool tA,
837  Matrix<Treal> const & A,
838  Vector<Treal> & x) {
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");
843  if (uplo != 'U')
844  throw Failure("Vector<class Treal>::"
845  "trmv only implemented for upper triangular matrices");
846  if ( ( A.is_zero() || x.is_zero() ) ) {
847  x = 0;
848  return;
849  }
850  int const ONE = 1;
851  if (!tA)
852  mat::trmv(&uplo, "N", "N", &x.n(), A.elements, &A.nrows(),
853  x.elements, &ONE);
854  else
855  mat::trmv(&uplo, "T", "N", &x.n(), A.elements, &A.nrows(),
856  x.elements, &ONE);
857  }
858 
859 
860 
861 
862 } /* end namespace mat */
863 #endif
#define A
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
Definition: Vector.h:503
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
Definition: Failure.h:57
Base class for Vector.
#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