00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00039 #ifndef MAT_VECTOR
00040 #define MAT_VECTOR
00041 #include "VectorHierarchicBase.h"
00042 namespace mat{
00043 template<class Treal, class Telement>
00044 class Matrix;
00053 template<class Treal, class Telement = Treal>
00054 class Vector: public VectorHierarchicBase<Treal, Telement> {
00055 public:
00056 typedef Telement ElementType;
00057
00058
00059 Vector():VectorHierarchicBase<Treal, Telement>(){}
00060
00061 void allocate() {
00062 assert(!this->is_empty());
00063 assert(this->is_zero());
00064 this->elements = allocateElements<Telement>(this->n());
00065 SizesAndBlocks rowSAB;
00066 for (int row = 0; row < this->rows.getNBlocks(); row++) {
00067 rowSAB = this->rows.getSizesAndBlocksForLowerLevel(row);
00068 (*this)(row).resetRows(rowSAB);
00069 }
00070 }
00071
00072 void assignFromFull(std::vector<Treal> const & fullVector);
00073
00074 void addFromFull(std::vector<Treal> const & fullVector);
00075
00076 void fullVector(std::vector<Treal> & fullVector) const;
00077
00078 Vector<Treal, Telement>&
00079 operator=(const Vector<Treal, Telement>& vec) {
00080 VectorHierarchicBase<Treal, Telement>::operator=(vec);
00081 return *this;
00082 }
00083
00084 void clear();
00085
00086 void writeToFile(std::ofstream & file) const;
00087 void readFromFile(std::ifstream & file);
00088 Vector<Treal, Telement>& operator=(int const k);
00089
00090 inline void randomNormalized() {
00091 this->random();
00092 (*this) *= (1.0 / this->eucl());
00093 }
00094 void random();
00095
00096 inline Treal eucl() const {
00097 return template_blas_sqrt(dot(*this,*this));
00098 }
00099
00100
00101 Vector<Treal, Telement>& operator*=(const Treal alpha);
00102 static Treal dot(Vector<Treal, Telement> const & x,
00103 Vector<Treal, Telement> const & y);
00104
00105
00106 static void axpy(Treal const & alpha,
00107 Vector<Treal, Telement> const & x,
00108 Vector<Treal, Telement> & y);
00109
00110
00111
00116 template<typename TmatrixElement>
00117 static void gemv(bool const tA, Treal const alpha,
00118 Matrix<Treal, TmatrixElement> const & A,
00119 Vector<Treal, Telement> const & x,
00120 Treal const beta,
00121 Vector<Treal, Telement>& y);
00122
00126 template<typename TmatrixElement>
00127 static void symv(char const uplo, Treal const alpha,
00128 Matrix<Treal, TmatrixElement> const & A,
00129 Vector<Treal, Telement> const & x,
00130 Treal const beta,
00131 Vector<Treal, Telement>& y);
00136 template<typename TmatrixElement>
00137 static void trmv(char const uplo, const bool tA,
00138 Matrix<Treal, TmatrixElement> const & A,
00139 Vector<Treal, Telement> & x);
00140
00141
00142 #if 0
00143 void assign_from_full(Treal const * const fullvector, const int totn);
00144
00145 void fullvector(Treal * const full, const int totn) const;
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155 #endif
00156 };
00157
00158
00159 template<class Treal, class Telement>
00160 void Vector<Treal, Telement>::
00161 assignFromFull(std::vector<Treal> const & fullVector) {
00162 addFromFull(fullVector);
00163 }
00164
00165 template<class Treal, class Telement>
00166 void Vector<Treal, Telement>::
00167 addFromFull(std::vector<Treal> const & fullVector) {
00168 if (this->is_zero())
00169 allocate();
00170 for (int ind = 0; ind < this->n(); ind++)
00171 (*this)(ind).addFromFull(fullVector);
00172 }
00173
00174 template<class Treal, class Telement>
00175 void Vector<Treal, Telement>::
00176 fullVector(std::vector<Treal> & fullVec) const {
00177 if (this->is_zero()) {
00178 fullVec.resize(this->rows.getNTotalScalars());
00179 for (int row = 0; row < this->nScalars(); ++row )
00180 fullVec[this->rows.getOffset()+row] = 0;
00181 }
00182 else
00183 for (int ind = 0; ind < this->n(); ind++)
00184 (*this)(ind).fullVector(fullVec);
00185 }
00186
00187
00188 template<class Treal, class Telement>
00189 void Vector<Treal, Telement>::clear() {
00190 freeElements(this->elements);
00191 this->elements = 0;
00192 }
00193
00194 template<class Treal, class Telement>
00195 void Vector<Treal, Telement>::
00196 writeToFile(std::ofstream & file) const {
00197 int const ZERO = 0;
00198 int const ONE = 1;
00199 if (this->is_zero()) {
00200 char * tmp = (char*)&ZERO;
00201 file.write(tmp,sizeof(int));
00202 }
00203 else {
00204 char * tmp = (char*)&ONE;
00205 file.write(tmp,sizeof(int));
00206 for (int i = 0; i < this->n(); i++)
00207 this->elements[i].writeToFile(file);
00208 }
00209 }
00210 template<class Treal, class Telement>
00211 void Vector<Treal, Telement>::
00212 readFromFile(std::ifstream & file) {
00213 int const ZERO = 0;
00214 int const ONE = 1;
00215 char tmp[sizeof(int)];
00216 file.read(tmp, (std::ifstream::pos_type)sizeof(int));
00217 switch ((int)*tmp) {
00218 case ZERO:
00219 (*this) = 0;
00220 break;
00221 case ONE:
00222 if (this->is_zero())
00223 allocate();
00224 for (int i = 0; i < this->n(); i++)
00225 this->elements[i].readFromFile(file);
00226 break;
00227 default:
00228 throw Failure("Vector<Treal, Telement>::"
00229 "readFromFile(std::ifstream & file):"
00230 "File corruption int value not 0 or 1");
00231 }
00232 }
00233
00234 template<class Treal, class Telement>
00235 Vector<Treal, Telement>& Vector<Treal, Telement>::
00236 operator=(int const k) {
00237 if (k == 0)
00238 this->clear();
00239 else
00240 throw Failure("Vector::operator=(int k) only "
00241 "implemented for k = 0");
00242 return *this;
00243 }
00244
00245 template<class Treal, class Telement>
00246 void Vector<Treal, Telement>::random() {
00247 if (this->is_zero())
00248 allocate();
00249 for (int ind = 0; ind < this->n(); ind++)
00250 (*this)(ind).random();
00251 }
00252
00253
00254
00255 template<class Treal, class Telement>
00256 Vector<Treal, Telement>& Vector<Treal, Telement>::
00257 operator*=(const Treal alpha) {
00258 if (!this->is_zero() && alpha != 1) {
00259 for (int ind = 0; ind < this->n(); ind++)
00260 (*this)(ind) *= alpha;
00261 }
00262 return *this;
00263 }
00264
00265 template<class Treal, class Telement>
00266 Treal Vector<Treal, Telement>::
00267 dot(Vector<Treal, Telement> const & x,
00268 Vector<Treal, Telement> const & y) {
00269 assert(x.n() == y.n());
00270 if (x.is_zero() || y.is_zero())
00271 return 0;
00272 Treal dotProduct = 0;
00273 for (int ind = 0; ind < x.n(); ind++)
00274 dotProduct += Telement::dot(x(ind), y(ind));
00275 return dotProduct;
00276 }
00277
00278
00279 template<class Treal, class Telement>
00280 void Vector<Treal, Telement>::
00281 axpy(Treal const & alpha,
00282 Vector<Treal, Telement> const & x,
00283 Vector<Treal, Telement> & y) {
00284 assert(x.n() == y.n());
00285 if (x.is_zero())
00286 return;
00287 if (y.is_zero()) {
00288 y.allocate();
00289 }
00290 for (int ind = 0; ind < x.n(); ind++)
00291 Telement::axpy(alpha, x(ind), y(ind));
00292 }
00293
00294
00295
00300 template<class Treal, class Telement>
00301 template<typename TmatrixElement>
00302 void Vector<Treal, Telement>::
00303 gemv(bool const tA, Treal const alpha,
00304 Matrix<Treal, TmatrixElement> const & A,
00305 Vector<Treal, Telement> const & x,
00306 Treal const beta,
00307 Vector<Treal, Telement>& y) {
00308 if (y.is_empty()) {
00309 assert(beta == 0);
00310 y.resetRows(x.rows);
00311 }
00312 if ((A.is_zero() || x.is_zero() || alpha == 0) &&
00313 (y.is_zero() || beta == 0))
00314 y = 0;
00315 else {
00316 Treal beta_tmp = beta;
00317 if (y.is_zero()) {
00318 y.allocate();
00319 beta_tmp = 0;
00320 }
00321 if (A.is_zero() || x.is_zero() || alpha == 0)
00322 y *= beta_tmp;
00323 else {
00324 MAT_OMP_INIT;
00325 if (!tA) {
00326 if (A.ncols() != x.n() || A.nrows() != y.n())
00327 throw Failure("Vector<Treal, Telement>::"
00328 "gemv(bool const, Treal const, "
00329 "const Matrix<Treal, Telement>&, "
00330 "const Vector<Treal, Telement>&, "
00331 "Treal const, const Vector<Treal, "
00332 "Telement>&): "
00333 "Incorrect dimensions for matrix-vector product");
00334 else {
00335 int A_nrows = A.nrows();
00336 #ifdef _OPENMP
00337 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic)
00338 #endif
00339 for (int row = 0; row < A_nrows; row++) {
00340 MAT_OMP_START;
00341 Telement::gemv(tA, alpha, A(row, 0), x(0), beta_tmp, y(row));
00342 for (int col = 1; col < A.ncols(); col++)
00343 Telement::gemv(tA, alpha, A(row, col), x(col), 1.0, y(row));
00344 MAT_OMP_END;
00345 }
00346 }
00347 }
00348 else {
00349 assert(tA);
00350 if (A.nrows() != x.n() || A.ncols() != y.n())
00351 throw Failure("Vector<Treal, Telement>::"
00352 "gemv(bool const, Treal const, "
00353 "const Matrix<Treal, Telement>&, "
00354 "const Vector<Treal, Telement>&, "
00355 "Treal const, const Vector<Treal, "
00356 "Telement>&): "
00357 "Incorrect dimensions for matrix-vector product");
00358 else {
00359 int A_ncols = A.ncols();
00360 #ifdef _OPENMP
00361 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic)
00362 #endif
00363 for (int col = 0; col < A_ncols; col++) {
00364 MAT_OMP_START;
00365 Telement::gemv(tA, alpha, A(0, col), x(0), beta_tmp, y(col));
00366 for (int row = 1; row < A.nrows(); row++)
00367 Telement::gemv(tA, alpha, A(row, col), x(row), 1.0, y(col));
00368 MAT_OMP_END;
00369 }
00370 }
00371 }
00372 MAT_OMP_FINALIZE;
00373 }
00374 }
00375 }
00376
00380 template<class Treal, class Telement>
00381 template<typename TmatrixElement>
00382 void Vector<Treal, Telement>::
00383 symv(char const uplo, Treal const alpha,
00384 Matrix<Treal, TmatrixElement> const & A,
00385 Vector<Treal, Telement> const & x,
00386 Treal const beta,
00387 Vector<Treal, Telement>& y) {
00388 if (y.is_empty()) {
00389 assert(beta == 0);
00390 y.resetRows(x.rows);
00391 }
00392 if (x.n() != y.n() || A.nrows() != A.ncols() || A.ncols() != x.n())
00393 throw Failure("Vector<Treal, Telement>::"
00394 "symv(char const uplo, Treal const, "
00395 "const Matrix<Treal, Telement>&, "
00396 "const Vector<Treal, Telement>&, "
00397 "Treal const, const Vector<Treal, Telement>&):"
00398 "Incorrect dimensions for symmetric "
00399 "matrix-vector product");
00400 if (uplo != 'U')
00401 throw Failure("Vector<class Treal, class Telement>::"
00402 "symv only implemented for symmetric matrices in "
00403 "upper triangular storage");
00404 if ((A.is_zero() || x.is_zero() || alpha == 0) &&
00405 (y.is_zero() || beta == 0))
00406 y = 0;
00407 else {
00408 Treal beta_tmp = beta;
00409 if (y.is_zero()) {
00410 y.allocate();
00411 beta_tmp = 0;
00412 }
00413 if (A.is_zero() || x.is_zero() || alpha == 0)
00414 y *= beta_tmp;
00415 else {
00416 MAT_OMP_INIT;
00417 #ifdef _OPENMP
00418 #pragma omp parallel if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared)
00419 #endif
00420 {
00421
00422 int A_ncols = A.ncols();
00423 #ifdef _OPENMP
00424 #pragma omp for schedule(dynamic)
00425 #endif
00426 for (int rc = 0; rc < A_ncols; rc++) {
00427 MAT_OMP_START;
00428 Telement::symv(uplo, alpha, A(rc,rc), x(rc), beta_tmp, y(rc));
00429 MAT_OMP_END;
00430 }
00431
00432 int A_nrows = A.nrows();
00433 #ifdef _OPENMP
00434 #pragma omp for schedule(dynamic)
00435 #endif
00436 for (int row = 0; row < A_nrows - 1; row++) {
00437 MAT_OMP_START;
00438 for (int col = row + 1; col < A.ncols(); col++)
00439 Telement::gemv(false, alpha, A(row, col), x(col), 1.0, y(row));
00440 MAT_OMP_END;
00441 }
00442
00443 #ifdef _OPENMP
00444 #pragma omp for schedule(dynamic)
00445 #endif
00446 for (int row = 1; row < A_nrows; row++) {
00447 MAT_OMP_START;
00448 for (int col = 0; col < row; col++)
00449 Telement::gemv(true, alpha, A(col, row), x(col), 1.0, y(row));
00450 MAT_OMP_END;
00451 }
00452 }
00453 MAT_OMP_FINALIZE;
00454 }
00455 }
00456 }
00457
00458 template<class Treal, class Telement>
00459 template<typename TmatrixElement>
00460 void Vector<Treal, Telement>::
00461 trmv(char const uplo, const bool tA,
00462 Matrix<Treal, TmatrixElement> const & A,
00463 Vector<Treal, Telement> & x) {
00464 if (A.nrows() != A.ncols() || A.ncols() != x.n())
00465 throw Failure("Vector<Treal, Telement>::"
00466 "trmv(...):"
00467 "Incorrect dimensions for triangular "
00468 "matrix-vector product");
00469 if (uplo != 'U')
00470 throw Failure("Vector<class Treal, class Telement>::"
00471 "trmv only implemented for upper triangular matrices");
00472 if ( ( A.is_zero() || x.is_zero() ) ) {
00473 x = 0;
00474 return;
00475 }
00476 if (!tA) {
00477
00478 for (int row = 0; row < A.nrows(); row++) {
00479 Telement::trmv(uplo, tA, A(row,row), x(row));
00480 for (int col = row + 1; col < A.ncols(); col++)
00481 Telement::gemv(tA, (Treal)1.0, A(row, col), x(col), 1.0, x(row));
00482 }
00483 return;
00484 }
00485
00486 for (int col = A.ncols() - 1; col >= 0; col--) {
00487 Telement::trmv(uplo, tA, A(col,col), x(col));
00488 for (int row = 0; row < col; row++)
00489 Telement::gemv(tA, (Treal)1.0, A(row, col), x(row), 1.0, x(col));
00490 }
00491 }
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502 template<class Treal>
00503 class Vector<Treal>: public VectorHierarchicBase<Treal> {
00504 public:
00505 friend class Matrix<Treal>;
00506 Vector()
00507 :VectorHierarchicBase<Treal>(){}
00508
00509 void allocate() {
00510 assert(!this->is_empty());
00511 assert(this->is_zero());
00512 this->elements = allocateElements<Treal>(this->n());
00513 for (int ind = 0; ind < this->n(); ind++)
00514 this->elements[ind] = 0;
00515 }
00516
00517 void assignFromFull(std::vector<Treal> const & fullVector);
00518
00519 void addFromFull(std::vector<Treal> const & fullVector);
00520
00521 void fullVector(std::vector<Treal> & fullVector) const;
00522
00523
00524 Vector<Treal>&
00525 operator=(const Vector<Treal>& vec) {
00526 VectorHierarchicBase<Treal>::operator=(vec);
00527 return *this;
00528 }
00529
00530 void clear();
00532 void writeToFile(std::ofstream & file) const;
00533 void readFromFile(std::ifstream & file);
00534
00535 Vector<Treal>& operator=(int const k);
00536
00537
00538 inline void randomNormalized() {
00539 this->random();
00540 (*this) *= 1 / this->eucl();
00541 }
00542 void random();
00543
00544 inline Treal eucl() const {
00545 return template_blas_sqrt(dot(*this,*this));
00546 }
00547
00548
00549 Vector<Treal>& operator*=(const Treal alpha);
00550
00551 static Treal dot(Vector<Treal> const & x,
00552 Vector<Treal> const & y);
00553
00554
00555
00556 static void axpy(Treal const & alpha,
00557 Vector<Treal> const & x,
00558 Vector<Treal> & y);
00559
00560
00565 static void gemv(bool const tA, Treal const alpha,
00566 Matrix<Treal> const & A,
00567 Vector<Treal> const & x,
00568 Treal const beta,
00569 Vector<Treal>& y);
00570
00574 static void symv(char const uplo, Treal const alpha,
00575 Matrix<Treal> const & A,
00576 Vector<Treal> const & x,
00577 Treal const beta,
00578 Vector<Treal>& y);
00579
00584 static void trmv(char const uplo, const bool tA,
00585 Matrix<Treal> const & A,
00586 Vector<Treal> & x);
00587
00588 };
00589
00590
00591 template<class Treal>
00592 void Vector<Treal>::
00593 assignFromFull(std::vector<Treal> const & fullVector) {
00594 addFromFull(fullVector);
00595 }
00596
00597 template<class Treal>
00598 void Vector<Treal>::
00599 addFromFull(std::vector<Treal> const & fullVector) {
00600 if (this->is_zero())
00601 allocate();
00602 assert((unsigned)this->rows.getNTotalScalars() == fullVector.size());
00603
00604
00605
00606 for (int row = 0; row < this->n(); ++row )
00607 (*this)(row) += fullVector[this->rows.getOffset()+row];
00608 }
00609
00610 template<class Treal>
00611 void Vector<Treal>::
00612 fullVector(std::vector<Treal> & fullVec) const {
00613 fullVec.resize(this->rows.getNTotalScalars());
00614 if (this->is_zero())
00615 for (int row = 0; row < this->nScalars(); ++row )
00616 fullVec[this->rows.getOffset()+row] = 0;
00617 else
00618 for (int row = 0; row < this->n(); ++row )
00619 fullVec[this->rows.getOffset()+row] = (*this)(row);
00620 }
00621
00622
00623 template<class Treal>
00624 void Vector<Treal>::clear() {
00625 freeElements(this->elements);
00626 this->elements = 0;
00627 }
00628
00629
00630 template<class Treal>
00631 void Vector<Treal>::
00632 writeToFile(std::ofstream & file) const {
00633 int const ZERO = 0;
00634 int const ONE = 1;
00635 if (this->is_zero()) {
00636 char * tmp = (char*)&ZERO;
00637 file.write(tmp,sizeof(int));
00638 }
00639 else {
00640 char * tmp = (char*)&ONE;
00641 file.write(tmp,sizeof(int));
00642 char * tmpel = (char*)this->elements;
00643 file.write(tmpel,sizeof(Treal) * this->n());
00644 }
00645 }
00646
00647 template<class Treal>
00648 void Vector<Treal>::
00649 readFromFile(std::ifstream & file) {
00650 int const ZERO = 0;
00651 int const ONE = 1;
00652 char tmp[sizeof(int)];
00653 file.read(tmp, (std::ifstream::pos_type)sizeof(int));
00654 switch ((int)*tmp) {
00655 case ZERO:
00656 (*this) = 0;
00657 break;
00658 case ONE:
00659 if (this->is_zero())
00660 allocate();
00661 file.read((char*)this->elements, sizeof(Treal) * this->n());
00662 break;
00663 default:
00664 throw Failure("Vector<Treal>::"
00665 "readFromFile(std::ifstream & file):"
00666 "File corruption, int value not 0 or 1");
00667 }
00668 }
00669
00670 template<class Treal>
00671 Vector<Treal>& Vector<Treal>::
00672 operator=(int const k) {
00673 if (k == 0)
00674 this->clear();
00675 else
00676 throw Failure("Vector::operator=(int k) only implemented for k = 0");
00677 return *this;
00678 }
00679
00680 template<class Treal>
00681 void Vector<Treal>::random() {
00682 if (this->is_zero())
00683 allocate();
00684 for (int ind = 0; ind < this->n(); ind++)
00685 (*this)(ind) = rand() / (Treal)RAND_MAX;
00686 }
00687
00688
00689 template<class Treal>
00690 Vector<Treal>& Vector<Treal>::
00691 operator*=(const Treal alpha) {
00692 if (!this->is_zero() && alpha != 1) {
00693 int const ONE = 1;
00694 mat::scal(&this->n(),&alpha,this->elements,&ONE);
00695 }
00696 return *this;
00697 }
00698
00699 template<class Treal>
00700 Treal Vector<Treal>::
00701 dot(Vector<Treal> const & x,
00702 Vector<Treal> const & y) {
00703 assert(x.n() == y.n());
00704 if (x.is_zero() || y.is_zero())
00705 return 0;
00706 else {
00707 int const ONE = 1;
00708 return mat::dot(&x.n(), x.elements, &ONE, y.elements, &ONE);
00709 }
00710 }
00711
00712
00713 template<class Treal>
00714 void Vector<Treal>::
00715 axpy(Treal const & alpha,
00716 Vector<Treal> const & x,
00717 Vector<Treal> & y) {
00718 assert(x.n() == y.n());
00719 if (x.is_zero())
00720 return;
00721 if (y.is_zero()) {
00722 y.allocate();
00723 for (int ind = 0; ind < y.n(); ind++)
00724 y.elements[ind] = 0;
00725 }
00726 int const ONE = 1;
00727 mat::axpy(&x.n(), &alpha, x.elements, &ONE, y.elements, &ONE);
00728 }
00729
00730
00731
00736 template<class Treal>
00737 void Vector<Treal>::
00738 gemv(bool const tA, Treal const alpha,
00739 Matrix<Treal> const & A,
00740 Vector<Treal> const & x,
00741 Treal const beta,
00742 Vector<Treal>& y) {
00743 if (y.is_empty()) {
00744 assert(beta == 0);
00745 y.resetRows(x.rows);
00746 }
00747 if ((A.is_zero() || x.is_zero() || alpha == 0) &&
00748 (y.is_zero() || beta == 0))
00749 y = 0;
00750 else {
00751 Treal beta_tmp = beta;
00752 if (y.is_zero()) {
00753 y.allocate();
00754 beta_tmp = 0;
00755 }
00756 if (A.is_zero() || x.is_zero() || alpha == 0)
00757 y *= beta_tmp;
00758 else {
00759 int const ONE = 1;
00760 if (!tA) {
00761 if (A.ncols() != x.n() || A.nrows() != y.n())
00762 throw Failure("Vector<Treal, Telement>::"
00763 "gemv(bool const, Treal const, "
00764 "const Matrix<Treal, Telement>&, "
00765 "const Vector<Treal, Telement>&, "
00766 "Treal const, const Vector<Treal, "
00767 "Telement>&): "
00768 "Incorrect dimensions for matrix-vector product");
00769 else {
00770 mat::gemv("N", &A.nrows(), &A.ncols(), &alpha, A.elements,
00771 &A.nrows(),x.elements,&ONE,&beta_tmp,y.elements,&ONE);
00772 }
00773 }
00774 else {
00775 assert(tA);
00776 if (A.nrows() != x.n() || A.ncols() != y.n())
00777 throw Failure("Vector<Treal, Telement>::"
00778 "gemv(bool const, Treal const, "
00779 "const Matrix<Treal, Telement>&, "
00780 "const Vector<Treal, Telement>&, "
00781 "Treal const, const Vector<Treal, "
00782 "Telement>&): "
00783 "Incorrect dimensions for matrix-vector product");
00784 else {
00785 mat::gemv("T", &A.nrows(), &A.ncols(), &alpha, A.elements,
00786 &A.nrows(),x.elements,&ONE,&beta_tmp,y.elements,&ONE);
00787 }
00788 }
00789 }
00790 }
00791 }
00792
00796 template<class Treal>
00797 void Vector<Treal>::
00798 symv(char const uplo, Treal const alpha,
00799 Matrix<Treal> const & A,
00800 Vector<Treal> const & x,
00801 Treal const beta,
00802 Vector<Treal>& y) {
00803 if (y.is_empty()) {
00804 assert(beta == 0);
00805 y.resetRows(x.rows);
00806 }
00807 if (x.n() != y.n() || A.nrows() != A.ncols() || A.ncols() != x.n())
00808 throw Failure("Vector<Treal>::"
00809 "symv(char const uplo, Treal const, "
00810 "const Matrix<Treal>&, "
00811 "const Vector<Treal>&, "
00812 "Treal const, const Vector<Treal>&):"
00813 "Incorrect dimensions for symmetric "
00814 "matrix-vector product");
00815 if ((A.is_zero() || x.is_zero() || alpha == 0) &&
00816 (y.is_zero() || beta == 0))
00817 y = 0;
00818 else {
00819 Treal beta_tmp = beta;
00820 if (y.is_zero()) {
00821 y.allocate();
00822 beta_tmp = 0;
00823 }
00824 if (A.is_zero() || x.is_zero() || alpha == 0)
00825 y *= beta_tmp;
00826 else {
00827 int const ONE = 1;
00828 mat::symv(&uplo, &x.n(), &alpha, A.elements, &A.nrows(),
00829 x.elements, &ONE, &beta, y.elements, &ONE);
00830 }
00831 }
00832 }
00833
00834 template<class Treal>
00835 void Vector<Treal>::
00836 trmv(char const uplo, const bool tA,
00837 Matrix<Treal> const & A,
00838 Vector<Treal> & x) {
00839 if (A.nrows() != A.ncols() || A.ncols() != x.n())
00840 throw Failure("Vector<Treal>::"
00841 "trmv(...): Incorrect dimensions for triangular "
00842 "matrix-vector product");
00843 if (uplo != 'U')
00844 throw Failure("Vector<class Treal>::"
00845 "trmv only implemented for upper triangular matrices");
00846 if ( ( A.is_zero() || x.is_zero() ) ) {
00847 x = 0;
00848 return;
00849 }
00850 int const ONE = 1;
00851 if (!tA)
00852 mat::trmv(&uplo, "N", "N", &x.n(), A.elements, &A.nrows(),
00853 x.elements, &ONE);
00854 else
00855 mat::trmv(&uplo, "T", "N", &x.n(), A.elements, &A.nrows(),
00856 x.elements, &ONE);
00857 }
00858
00859
00860
00861
00862 }
00863 #endif