ergo
Vector.h
Go to the documentation of this file.
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