ergo
|
00001 /* Ergo, version 3.2, a program for linear scaling electronic structure 00002 * calculations. 00003 * Copyright (C) 2012 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek. 00004 * 00005 * This program is free software: you can redistribute it and/or modify 00006 * it under the terms of the GNU General Public License as published by 00007 * the Free Software Foundation, either version 3 of the License, or 00008 * (at your option) any later version. 00009 * 00010 * This program is distributed in the hope that it will be useful, 00011 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00012 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00013 * GNU General Public License for more details. 00014 * 00015 * You should have received a copy of the GNU General Public License 00016 * along with this program. If not, see <http://www.gnu.org/licenses/>. 00017 * 00018 * Primary academic reference: 00019 * KohnâSham Density Functional Theory Electronic Structure Calculations 00020 * with Linearly Scaling Computational Time and Memory Usage, 00021 * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek, 00022 * J. Chem. Theory Comput. 7, 340 (2011), 00023 * <http://dx.doi.org/10.1021/ct100611z> 00024 * 00025 * For further information about Ergo, see <http://www.ergoscf.org>. 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 // template<typename TmatrixElement> 00056 // friend class Matrix<Treal, TmatrixElement>; 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 /* LEVEL 1 operations */ 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 /* y += alpha * x */ 00104 static void axpy(Treal const & alpha, 00105 Vector<Treal, Telement> const & x, 00106 Vector<Treal, Telement> & y); 00107 00108 00109 /* LEVEL 2 operations */ 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 /* OLD ROUTINES */ 00141 void assign_from_full(Treal const * const fullvector, const int totn); 00142 /* Convert to full vector */ 00143 void fullvector(Treal * const full, const int totn) const; 00144 00145 00146 00147 00148 00149 00150 00151 00152 00153 #endif /* END OLD ROUTINES */ 00154 }; /* end class Vector */ 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 /* LEVEL 1 operations */ 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 /* y += alpha * x */ 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 /* LEVEL 2 operations */ 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 } /* end else */ 00344 } /* end if (!tA) */ 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 } /* end else */ 00367 } /* end else */ 00368 MAT_OMP_FINALIZE; 00369 } /* end else */ 00370 } /* end else */ 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 /* Diagonal */ 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 /* Upper triangle */ 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 /* Lower triangle */ 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 } /* end omp parallel*/ 00447 MAT_OMP_FINALIZE; 00448 } /* end else */ 00449 } /* end else */ 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 // not transposed 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 // transposed 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 /* Specialization for Telement = Treal */ 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 /* LEVEL 1 operations */ 00543 Vector<Treal>& operator*=(const Treal alpha); 00544 00545 static Treal dot(Vector<Treal> const & x, 00546 Vector<Treal> const & y); 00547 00548 00549 /* y += alpha * x */ 00550 static void axpy(Treal const & alpha, 00551 Vector<Treal> const & x, 00552 Vector<Treal> & y); 00553 00554 /* LEVEL 2 operations */ 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 }; /* end class Vector specialization */ 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 /* Assertion AFTER empty check done 00598 * by allocate() 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 /* LEVEL 1 operations */ 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 /* y += alpha * x */ 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; /* fill with zeros */ 00719 } 00720 int const ONE = 1; 00721 mat::axpy(&x.n(), &alpha, x.elements, &ONE, y.elements, &ONE); 00722 } 00723 00724 00725 /* LEVEL 2 operations */ 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 } /* end else */ 00767 } /* end if (!tA) */ 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 } /* end else */ 00782 } /* end else */ 00783 } /* end else */ 00784 } /* end else */ 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 } /* end else */ 00825 } /* end else */ 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 } /* end namespace mat */ 00857 #endif