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