ergo
Vector.h
Go to the documentation of this file.
1 /* Ergo, version 3.3, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2013 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek.
4  *
5  * This program is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation, either version 3 of the License, or
8  * (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program. If not, see <http://www.gnu.org/licenses/>.
17  *
18  * Primary academic reference:
19  * Kohn−Sham Density Functional Theory Electronic Structure Calculations
20  * with Linearly Scaling Computational Time and Memory Usage,
21  * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek,
22  * J. Chem. Theory Comput. 7, 340 (2011),
23  * <http://dx.doi.org/10.1021/ct100611z>
24  *
25  * For further information about Ergo, see <http://www.ergoscf.org>.
26  */
27 
37 #ifndef MAT_VECTOR
38 #define MAT_VECTOR
39 #include "VectorHierarchicBase.h"
40 namespace mat{
41  template<class Treal, class Telement>
42  class Matrix;
51  template<class Treal, class Telement = Treal>
52  class Vector: public VectorHierarchicBase<Treal, Telement> {
53  public:
54  typedef Telement ElementType;
55  // template<typename TmatrixElement>
56  // friend class Matrix<Treal, TmatrixElement>;
57  Vector():VectorHierarchicBase<Treal, Telement>(){}
58 
59  void allocate() {
60  assert(!this->is_empty());
61  assert(this->is_zero());
62  this->elements = allocateElements<Telement>(this->n());
63  SizesAndBlocks rowSAB;
64  for (int row = 0; row < this->rows.getNBlocks(); row++) {
65  rowSAB = this->rows.getSizesAndBlocksForLowerLevel(row);
66  (*this)(row).resetRows(rowSAB);
67  }
68  }
69 
70  void assignFromFull(std::vector<Treal> const & fullVector);
71 
72  void addFromFull(std::vector<Treal> const & fullVector);
73 
74  void fullVector(std::vector<Treal> & fullVector) const;
75 
79  return *this;
80  }
81 
82  void clear();
83 
84  void writeToFile(std::ofstream & file) const;
85  void readFromFile(std::ifstream & file);
86  Vector<Treal, Telement>& operator=(int const k);
87 
88  inline void randomNormalized() {
89  this->random();
90  (*this) *= (1.0 / this->eucl());
91  }
92  void random();
93 
94  inline Treal eucl() const {
95  return template_blas_sqrt(dot(*this,*this));
96  }
97 
98  /* LEVEL 1 operations */
99  Vector<Treal, Telement>& operator*=(const Treal alpha);
100  static Treal dot(Vector<Treal, Telement> const & x,
101  Vector<Treal, Telement> const & y);
102 
103  /* y += alpha * x */
104  static void axpy(Treal const & alpha,
105  Vector<Treal, Telement> const & x,
107 
108 
109  /* LEVEL 2 operations */
114  template<typename TmatrixElement>
115  static void gemv(bool const tA, Treal const alpha,
117  Vector<Treal, Telement> const & x,
118  Treal const beta,
120 
124  template<typename TmatrixElement>
125  static void symv(char const uplo, Treal const alpha,
127  Vector<Treal, Telement> const & x,
128  Treal const beta,
134  template<typename TmatrixElement>
135  static void trmv(char const uplo, const bool tA,
138 
139 
140 #if 0 /* OLD ROUTINES */
141  void assign_from_full(Treal const * const fullvector, const int totn);
142  /* Convert to full vector */
143  void fullvector(Treal * const full, const int totn) const;
144 
145 
146 
147 
148 
149 
150 
151 
152 
153 #endif /* END OLD ROUTINES */
154  }; /* end class Vector */
155 
156 
157  template<class Treal, class Telement>
159  assignFromFull(std::vector<Treal> const & fullVector) {
160  addFromFull(fullVector);
161  }
162 
163  template<class Treal, class Telement>
165  addFromFull(std::vector<Treal> const & fullVector) {
166  if (this->is_zero())
167  allocate();
168  for (int ind = 0; ind < this->n(); ind++)
169  (*this)(ind).addFromFull(fullVector);
170  }
171 
172  template<class Treal, class Telement>
174  fullVector(std::vector<Treal> & fullVec) const {
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;
179  }
180  else
181  for (int ind = 0; ind < this->n(); ind++)
182  (*this)(ind).fullVector(fullVec);
183  }
184 
185 
186  template<class Treal, class Telement>
188  freeElements(this->elements);
189  this->elements = 0;
190  }
191 
192  template<class Treal, class Telement>
194  writeToFile(std::ofstream & file) const {
195  int const ZERO = 0;
196  int const ONE = 1;
197  if (this->is_zero()) {
198  char * tmp = (char*)&ZERO;
199  file.write(tmp,sizeof(int));
200  }
201  else {
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);
206  }
207  }
208  template<class Treal, class Telement>
210  readFromFile(std::ifstream & file) {
211  int const ZERO = 0;
212  int const ONE = 1;
213  char tmp[sizeof(int)];
214  file.read(tmp, (std::ifstream::pos_type)sizeof(int));
215  switch ((int)*tmp) {
216  case ZERO:
217  (*this) = 0;
218  break;
219  case ONE:
220  if (this->is_zero())
221  allocate();
222  for (int i = 0; i < this->n(); i++)
223  this->elements[i].readFromFile(file);
224  break;
225  default:
226  throw Failure("Vector<Treal, Telement>::"
227  "readFromFile(std::ifstream & file):"
228  "File corruption int value not 0 or 1");
229  }
230  }
231 
232  template<class Treal, class Telement>
234  operator=(int const k) {
235  if (k == 0)
236  this->clear();
237  else
238  throw Failure("Vector::operator=(int k) only "
239  "implemented for k = 0");
240  return *this;
241  }
242 
243  template<class Treal, class Telement>
245  if (this->is_zero())
246  allocate();
247  for (int ind = 0; ind < this->n(); ind++)
248  (*this)(ind).random();
249  }
250 
251  /* LEVEL 1 operations */
252 
253  template<class Treal, class Telement>
255  operator*=(const Treal alpha) {
256  if (!this->is_zero() && alpha != 1) {
257  for (int ind = 0; ind < this->n(); ind++)
258  (*this)(ind) *= alpha;
259  }
260  return *this;
261  }
262 
263  template<class Treal, class Telement>
266  Vector<Treal, Telement> const & y) {
267  assert(x.n() == y.n());
268  if (x.is_zero() || y.is_zero())
269  return 0;
270  Treal dotProduct = 0;
271  for (int ind = 0; ind < x.n(); ind++)
272  dotProduct += Telement::dot(x(ind), y(ind));
273  return dotProduct;
274  }
275 
276  /* y += alpha * x */
277  template<class Treal, class Telement>
279  axpy(Treal const & alpha,
280  Vector<Treal, Telement> const & x,
282  assert(x.n() == y.n());
283  if (x.is_zero())
284  return;
285  if (y.is_zero()) {
286  y.allocate();
287  }
288  for (int ind = 0; ind < x.n(); ind++)
289  Telement::axpy(alpha, x(ind), y(ind));
290  }
291 
292  /* LEVEL 2 operations */
293 
298  template<class Treal, class Telement>
299  template<typename TmatrixElement>
301  gemv(bool const tA, Treal const alpha,
303  Vector<Treal, Telement> const & x,
304  Treal const beta,
306  if (y.is_empty()) {
307  assert(beta == 0);
308  y.resetRows(x.rows);
309  }
310  if ((A.is_zero() || x.is_zero() || alpha == 0) &&
311  (y.is_zero() || beta == 0))
312  y = 0;
313  else {
314  Treal beta_tmp = beta;
315  if (y.is_zero()) {
316  y.allocate();
317  beta_tmp = 0;
318  }
319  if (A.is_zero() || x.is_zero() || alpha == 0)
320  y *= beta_tmp;
321  else {
322  MAT_OMP_INIT;
323  if (!tA) {
324  if (A.ncols() != x.n() || A.nrows() != y.n())
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, "
330  "Telement>&): "
331  "Incorrect dimensions for matrix-vector product");
332  else {
333 #ifdef _OPENMP
334 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic)
335 #endif
336  for (int row = 0; row < A.nrows(); row++) {
338  Telement::gemv(tA, alpha, A(row, 0), x(0), beta_tmp, y(row));
339  for (int col = 1; col < A.ncols(); col++)
340  Telement::gemv(tA, alpha, A(row, col), x(col), 1.0, y(row));
341  MAT_OMP_END;
342  }
343  } /* end else */
344  } /* end if (!tA) */
345  else {
346  assert(tA);
347  if (A.nrows() != x.n() || A.ncols() != y.n())
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, "
353  "Telement>&): "
354  "Incorrect dimensions for matrix-vector product");
355  else {
356 #ifdef _OPENMP
357 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic)
358 #endif
359  for (int col = 0; col < A.ncols(); col++) {
361  Telement::gemv(tA, alpha, A(0, col), x(0), beta_tmp, y(col));
362  for (int row = 1; row < A.nrows(); row++)
363  Telement::gemv(tA, alpha, A(row, col), x(row), 1.0, y(col));
364  MAT_OMP_END;
365  }
366  } /* end else */
367  } /* end else */
369  } /* end else */
370  } /* end else */
371  }
372 
376  template<class Treal, class Telement>
377  template<typename TmatrixElement>
379  symv(char const uplo, Treal const alpha,
381  Vector<Treal, Telement> const & x,
382  Treal const beta,
384  if (y.is_empty()) {
385  assert(beta == 0);
386  y.resetRows(x.rows);
387  }
388  if (x.n() != y.n() || A.nrows() != A.ncols() || A.ncols() != x.n())
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");
396  if (uplo != 'U')
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))
402  y = 0;
403  else {
404  Treal beta_tmp = beta;
405  if (y.is_zero()) {
406  y.allocate();
407  beta_tmp = 0;
408  }
409  if (A.is_zero() || x.is_zero() || alpha == 0)
410  y *= beta_tmp;
411  else {
412  MAT_OMP_INIT;
413 #ifdef _OPENMP
414 #pragma omp parallel if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared)
415 #endif
416  {
417  /* Diagonal */
418 #ifdef _OPENMP
419 #pragma omp for schedule(dynamic)
420 #endif
421  for (int rc = 0; rc < A.ncols(); rc++) {
423  Telement::symv(uplo, alpha, A(rc,rc), x(rc), beta_tmp, y(rc));
424  MAT_OMP_END;
425  }
426  /* Upper triangle */
427 #ifdef _OPENMP
428 #pragma omp for schedule(dynamic)
429 #endif
430  for (int row = 0; row < A.nrows() - 1; row++) {
432  for (int col = row + 1; col < A.ncols(); col++)
433  Telement::gemv(false, alpha, A(row, col), x(col), 1.0, y(row));
434  MAT_OMP_END;
435  }
436  /* Lower triangle */
437 #ifdef _OPENMP
438 #pragma omp for schedule(dynamic)
439 #endif
440  for (int row = 1; row < A.nrows(); row++) {
442  for (int col = 0; col < row; col++)
443  Telement::gemv(true, alpha, A(col, row), x(col), 1.0, y(row));
444  MAT_OMP_END;
445  }
446  } /* end omp parallel*/
448  } /* end else */
449  } /* end else */
450  }
451 
452  template<class Treal, class Telement>
453  template<typename TmatrixElement>
455  trmv(char const uplo, const bool tA,
458  if (A.nrows() != A.ncols() || A.ncols() != x.n())
459  throw Failure("Vector<Treal, Telement>::"
460  "trmv(...):"
461  "Incorrect dimensions for triangular "
462  "matrix-vector product");
463  if (uplo != 'U')
464  throw Failure("Vector<class Treal, class Telement>::"
465  "trmv only implemented for upper triangular matrices");
466  if ( ( A.is_zero() || x.is_zero() ) ) {
467  x = 0;
468  return;
469  }
470  if (!tA) {
471  // not transposed
472  for (int row = 0; row < A.nrows(); row++) {
473  Telement::trmv(uplo, tA, A(row,row), x(row));
474  for (int col = row + 1; col < A.ncols(); col++)
475  Telement::gemv(tA, (Treal)1.0, A(row, col), x(col), 1.0, x(row));
476  }
477  return;
478  }
479  // transposed
480  for (int col = A.ncols() - 1; col >= 0; col--) {
481  Telement::trmv(uplo, tA, A(col,col), x(col));
482  for (int row = 0; row < col; row++)
483  Telement::gemv(tA, (Treal)1.0, A(row, col), x(row), 1.0, x(col));
484  }
485  }
486 
487 
488 
489 
490  /***************************************************************************/
491  /***************************************************************************/
492  /* Specialization for Telement = Treal */
493  /***************************************************************************/
494  /***************************************************************************/
495 
496  template<class Treal>
497  class Vector<Treal>: public VectorHierarchicBase<Treal> {
498  public:
499  friend class Matrix<Treal>;
501  :VectorHierarchicBase<Treal>(){}
502 
503  void allocate() {
504  assert(!this->is_empty());
505  assert(this->is_zero());
506  this->elements = allocateElements<Treal>(this->n());
507  for (int ind = 0; ind < this->n(); ind++)
508  this->elements[ind] = 0;
509  }
510 
511  void assignFromFull(std::vector<Treal> const & fullVector);
512 
513  void addFromFull(std::vector<Treal> const & fullVector);
514 
515  void fullVector(std::vector<Treal> & fullVector) const;
516 
517 
518  Vector<Treal>&
519  operator=(const Vector<Treal>& vec) {
521  return *this;
522  }
523 
524  void clear();
526  void writeToFile(std::ofstream & file) const;
527  void readFromFile(std::ifstream & file);
528 
529  Vector<Treal>& operator=(int const k);
530 
531 
532  inline void randomNormalized() {
533  this->random();
534  (*this) *= 1 / this->eucl();
535  }
536  void random();
537 
538  inline Treal eucl() const {
539  return template_blas_sqrt(dot(*this,*this));
540  }
541 
542  /* LEVEL 1 operations */
543  Vector<Treal>& operator*=(const Treal alpha);
544 
545  static Treal dot(Vector<Treal> const & x,
546  Vector<Treal> const & y);
547 
548 
549  /* y += alpha * x */
550  static void axpy(Treal const & alpha,
551  Vector<Treal> const & x,
552  Vector<Treal> & y);
553 
554  /* LEVEL 2 operations */
559  static void gemv(bool const tA, Treal const alpha,
560  Matrix<Treal> const & A,
561  Vector<Treal> const & x,
562  Treal const beta,
563  Vector<Treal>& y);
564 
568  static void symv(char const uplo, Treal const alpha,
569  Matrix<Treal> const & A,
570  Vector<Treal> const & x,
571  Treal const beta,
572  Vector<Treal>& y);
573 
578  static void trmv(char const uplo, const bool tA,
579  Matrix<Treal> const & A,
580  Vector<Treal> & x);
581 
582  }; /* end class Vector specialization */
583 
584 
585  template<class Treal>
586  void Vector<Treal>::
587  assignFromFull(std::vector<Treal> const & fullVector) {
588  addFromFull(fullVector);
589  }
590 
591  template<class Treal>
592  void Vector<Treal>::
593  addFromFull(std::vector<Treal> const & fullVector) {
594  if (this->is_zero())
595  allocate();
596  assert((unsigned)this->rows.getNTotalScalars() == fullVector.size());
597  /* Assertion AFTER empty check done
598  * by allocate()
599  */
600  for (int row = 0; row < this->n(); ++row )
601  (*this)(row) += fullVector[this->rows.getOffset()+row];
602  }
603 
604  template<class Treal>
605  void Vector<Treal>::
606  fullVector(std::vector<Treal> & fullVec) const {
607  fullVec.resize(this->rows.getNTotalScalars());
608  if (this->is_zero())
609  for (int row = 0; row < this->nScalars(); ++row )
610  fullVec[this->rows.getOffset()+row] = 0;
611  else
612  for (int row = 0; row < this->n(); ++row )
613  fullVec[this->rows.getOffset()+row] = (*this)(row);
614  }
615 
616 
617  template<class Treal>
619  freeElements(this->elements);
620  this->elements = 0;
621  }
622 
623 
624  template<class Treal>
625  void Vector<Treal>::
626  writeToFile(std::ofstream & file) const {
627  int const ZERO = 0;
628  int const ONE = 1;
629  if (this->is_zero()) {
630  char * tmp = (char*)&ZERO;
631  file.write(tmp,sizeof(int));
632  }
633  else {
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());
638  }
639  }
640 
641  template<class Treal>
642  void Vector<Treal>::
643  readFromFile(std::ifstream & file) {
644  int const ZERO = 0;
645  int const ONE = 1;
646  char tmp[sizeof(int)];
647  file.read(tmp, (std::ifstream::pos_type)sizeof(int));
648  switch ((int)*tmp) {
649  case ZERO:
650  (*this) = 0;
651  break;
652  case ONE:
653  if (this->is_zero())
654  allocate();
655  file.read((char*)this->elements, sizeof(Treal) * this->n());
656  break;
657  default:
658  throw Failure("Vector<Treal>::"
659  "readFromFile(std::ifstream & file):"
660  "File corruption, int value not 0 or 1");
661  }
662  }
663 
664  template<class Treal>
666  operator=(int const k) {
667  if (k == 0)
668  this->clear();
669  else
670  throw Failure("Vector::operator=(int k) only implemented for k = 0");
671  return *this;
672  }
673 
674  template<class Treal>
676  if (this->is_zero())
677  allocate();
678  for (int ind = 0; ind < this->n(); ind++)
679  (*this)(ind) = rand() / (Treal)RAND_MAX;
680  }
681 
682  /* LEVEL 1 operations */
683  template<class Treal>
685  operator*=(const Treal alpha) {
686  if (!this->is_zero() && alpha != 1) {
687  int const ONE = 1;
688  mat::scal(&this->n(),&alpha,this->elements,&ONE);
689  }
690  return *this;
691  }
692 
693  template<class Treal>
694  Treal Vector<Treal>::
695  dot(Vector<Treal> const & x,
696  Vector<Treal> const & y) {
697  assert(x.n() == y.n());
698  if (x.is_zero() || y.is_zero())
699  return 0;
700  else {
701  int const ONE = 1;
702  return mat::dot(&x.n(), x.elements, &ONE, y.elements, &ONE);
703  }
704  }
705 
706  /* y += alpha * x */
707  template<class Treal>
708  void Vector<Treal>::
709  axpy(Treal const & alpha,
710  Vector<Treal> const & x,
711  Vector<Treal> & y) {
712  assert(x.n() == y.n());
713  if (x.is_zero())
714  return;
715  if (y.is_zero()) {
716  y.allocate();
717  for (int ind = 0; ind < y.n(); ind++)
718  y.elements[ind] = 0; /* fill with zeros */
719  }
720  int const ONE = 1;
721  mat::axpy(&x.n(), &alpha, x.elements, &ONE, y.elements, &ONE);
722  }
723 
724 
725  /* LEVEL 2 operations */
730  template<class Treal>
731  void Vector<Treal>::
732  gemv(bool const tA, Treal const alpha,
733  Matrix<Treal> const & A,
734  Vector<Treal> const & x,
735  Treal const beta,
736  Vector<Treal>& y) {
737  if (y.is_empty()) {
738  assert(beta == 0);
739  y.resetRows(x.rows);
740  }
741  if ((A.is_zero() || x.is_zero() || alpha == 0) &&
742  (y.is_zero() || beta == 0))
743  y = 0;
744  else {
745  Treal beta_tmp = beta;
746  if (y.is_zero()) {
747  y.allocate();
748  beta_tmp = 0;
749  }
750  if (A.is_zero() || x.is_zero() || alpha == 0)
751  y *= beta_tmp;
752  else {
753  int const ONE = 1;
754  if (!tA) {
755  if (A.ncols() != x.n() || A.nrows() != y.n())
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, "
761  "Telement>&): "
762  "Incorrect dimensions for matrix-vector product");
763  else {
764  mat::gemv("N", &A.nrows(), &A.ncols(), &alpha, A.elements,
765  &A.nrows(),x.elements,&ONE,&beta_tmp,y.elements,&ONE);
766  } /* end else */
767  } /* end if (!tA) */
768  else {
769  assert(tA);
770  if (A.nrows() != x.n() || A.ncols() != y.n())
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, "
776  "Telement>&): "
777  "Incorrect dimensions for matrix-vector product");
778  else {
779  mat::gemv("T", &A.nrows(), &A.ncols(), &alpha, A.elements,
780  &A.nrows(),x.elements,&ONE,&beta_tmp,y.elements,&ONE);
781  } /* end else */
782  } /* end else */
783  } /* end else */
784  } /* end else */
785  }
786 
790  template<class Treal>
791  void Vector<Treal>::
792  symv(char const uplo, Treal const alpha,
793  Matrix<Treal> const & A,
794  Vector<Treal> const & x,
795  Treal const beta,
796  Vector<Treal>& y) {
797  if (y.is_empty()) {
798  assert(beta == 0);
799  y.resetRows(x.rows);
800  }
801  if (x.n() != y.n() || A.nrows() != A.ncols() || A.ncols() != x.n())
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");
809  if ((A.is_zero() || x.is_zero() || alpha == 0) &&
810  (y.is_zero() || beta == 0))
811  y = 0;
812  else {
813  Treal beta_tmp = beta;
814  if (y.is_zero()) {
815  y.allocate();
816  beta_tmp = 0;
817  }
818  if (A.is_zero() || x.is_zero() || alpha == 0)
819  y *= beta_tmp;
820  else {
821  int const ONE = 1;
822  mat::symv(&uplo, &x.n(), &alpha, A.elements, &A.nrows(),
823  x.elements, &ONE, &beta, y.elements, &ONE);
824  } /* end else */
825  } /* end else */
826  }
827 
828  template<class Treal>
829  void Vector<Treal>::
830  trmv(char const uplo, const bool tA,
831  Matrix<Treal> const & A,
832  Vector<Treal> & x) {
833  if (A.nrows() != A.ncols() || A.ncols() != x.n())
834  throw Failure("Vector<Treal>::"
835  "trmv(...): Incorrect dimensions for triangular "
836  "matrix-vector product");
837  if (uplo != 'U')
838  throw Failure("Vector<class Treal>::"
839  "trmv only implemented for upper triangular matrices");
840  if ( ( A.is_zero() || x.is_zero() ) ) {
841  x = 0;
842  return;
843  }
844  int const ONE = 1;
845  if (!tA)
846  mat::trmv(&uplo, "N", "N", &x.n(), A.elements, &A.nrows(),
847  x.elements, &ONE);
848  else
849  mat::trmv(&uplo, "T", "N", &x.n(), A.elements, &A.nrows(),
850  x.elements, &ONE);
851  }
852 
853 
854 
855 
856 } /* end namespace mat */
857 #endif
#define A
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
Definition: allocate.cc:30
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
Definition: Vector.h:497
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
Definition: Failure.h:47
int const & getNBlocks() const
Definition: SizesAndBlocks.h:64
Base class for Vector.
#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