ergo
Vector.h
Go to the documentation of this file.
1/* Ergo, version 3.8.2, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2023 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
38
39#ifndef MAT_VECTOR
40#define MAT_VECTOR
42namespace 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
83
84 void clear();
85
86 void writeToFile(std::ofstream & file) const;
87 void readFromFile(std::ifstream & file);
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 */
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) {
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 {
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));
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));
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 {
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));
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));
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));
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
527 return *this;
528 }
529
530 void clear();
531
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>
593 assignFromFull(std::vector<Treal> const & fullVector) {
595 }
596
597 template<class 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>
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>
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>
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>
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>
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>
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>
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>
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
Base class for Vector.
Definition Failure.h:57
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
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
Definition VectorHierarchicBase.h:89
const int & nScalars() const
Definition VectorHierarchicBase.h:55
Telement * elements
Definition VectorHierarchicBase.h:108
bool is_zero() const
Definition VectorHierarchicBase.h:75
VectorHierarchicBase()
Definition VectorHierarchicBase.h:93
VectorHierarchicBase< Treal, Telement > & operator=(const VectorHierarchicBase< Treal, Telement > &vec)
Definition VectorHierarchicBase.h:152
Vector()
Definition Vector.h:506
void randomNormalized()
Definition Vector.h:538
void addFromFull(std::vector< Treal > const &fullVector)
Definition Vector.h:599
void clear()
Set vector to zero and delete all arrays.
Definition Vector.h:624
void fullVector(std::vector< Treal > &fullVector) const
Definition Vector.h:612
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
typename ElementType::VectorType 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
Vector< Treal, Telement > & operator=(int const k)
Definition Vector.h:236
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
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
#define A
#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
const Treal Matrix< Treal >::ONE
Definition Matrix.h:3448
static T dot(const int *n, const T *dx, const int *incx, const T *dy, const int *incy)
Definition mat_gblas.h:425
T * allocateElements(int n)
Definition allocate.h:42
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
const Treal Matrix< Treal >::ZERO
Definition Matrix.h:3446
Treal template_blas_sqrt(Treal x)