VectorGeneral.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 
00036 #ifndef MAT_VECTORGENERAL
00037 #define MAT_VECTORGENERAL
00038 #include <iostream>
00039 #include <fstream>
00040 #include <ios>
00041 #include "FileWritable.h"
00042 #include "matrix_proxy.h"
00043 #include "ValidPtr.h"
00044 namespace mat {
00045   template<typename Treal, typename Tvector>
00046     class VectorGeneral : public FileWritable {
00047   public:
00048 
00049     inline void resetSizesAndBlocks(SizesAndBlocks const & newRows) {
00050       vectorPtr.haveDataStructureSet(true);
00051       vectorPtr->resetRows(newRows);
00052     }
00053 
00054     inline bool is_empty() const {
00055       return !vectorPtr.haveDataStructureGet();
00056     }
00057 
00058 
00059     VectorGeneral():vectorPtr(new Tvector) {}
00060     explicit VectorGeneral(const VectorGeneral<Treal, Tvector>& other)
00061       :FileWritable(other), vectorPtr(new Tvector) {
00062       if (other.vectorPtr.haveDataStructureGet()) {
00063         vectorPtr.haveDataStructureSet(true);
00064       }
00065       *vectorPtr = *other.vectorPtr;
00066     }
00067       
00068     inline void assign_from_full
00069       (std::vector<Treal> const & fullVector, 
00070        SizesAndBlocks const & newRows) {
00071       resetSizesAndBlocks(newRows);
00072       this->vectorPtr->assignFromFull(fullVector);
00073     }
00074     inline void fullvector(std::vector<Treal> & fullVector) const {
00075       this->vectorPtr->fullVector(fullVector);
00076     }
00077     VectorGeneral<Treal, Tvector>& 
00078       operator=(const VectorGeneral<Treal, Tvector>& other) {
00079       if (other.vectorPtr.haveDataStructureGet()) {
00080         vectorPtr.haveDataStructureSet(true);
00081       }
00082       *this->vectorPtr = *other.vectorPtr; 
00083       return *this;
00084     } 
00085 
00086     inline void clear() {
00087       if (is_empty())
00088         // This means that the object's data structure has not been set
00089         // There is nothing to clear and the vectorPtr is not valid either
00090         return;
00091       vectorPtr->clear();
00092     }
00093 
00094     inline void rand() {
00095       vectorPtr->randomNormalized();
00096     }
00097     
00098     /* LEVEL 2 operations */
00099     /* OPERATIONS INVOLVING ORDINARY MATRICES */
00101     template<typename Tmatrix>
00102       VectorGeneral<Treal, Tvector>& operator=
00103       (const XYZ<Treal,
00104        MatrixGeneral<Treal, Tmatrix>,
00105        VectorGeneral<Treal, Tvector> >& smv);
00106     
00108     template<typename Tmatrix>
00109       VectorGeneral<Treal, Tvector>& operator+=
00110       (const XYZ<Treal,
00111        MatrixGeneral<Treal, Tmatrix>,
00112        VectorGeneral<Treal, Tvector> >& smv);
00114     template<typename Tmatrix>
00115       VectorGeneral<Treal, Tvector>& operator=
00116       (const XYZpUV<Treal,
00117        MatrixGeneral<Treal, Tmatrix>,
00118        VectorGeneral<Treal, Tvector>,
00119        Treal,
00120        VectorGeneral<Treal, Tvector> >& smvpsv);
00121 
00123     template<typename Tmatrix>
00124       VectorGeneral<Treal, Tvector>& operator=
00125       (const XY<MatrixGeneral<Treal, Tmatrix>,
00126        VectorGeneral<Treal, Tvector> >& mv) {
00127       Treal ONE = 1.0;
00128       return this->operator=(XYZ<Treal, MatrixGeneral<Treal, Tmatrix>,
00129                              VectorGeneral<Treal, Tvector> >(ONE, mv.A, mv.B,
00130                                                              false, mv.tA, mv.tB));      
00131     }
00132     
00133     /* OPERATIONS INVOLVING SYMMETRIC MATRICES */
00135     template<typename Tmatrix>
00136       VectorGeneral<Treal, Tvector>& operator=
00137       (const XYZ<Treal,
00138        MatrixSymmetric<Treal, Tmatrix>,
00139        VectorGeneral<Treal, Tvector> >& smv);
00141     template<typename Tmatrix>
00142       VectorGeneral<Treal, Tvector>& operator+=
00143       (const XYZ<Treal,
00144        MatrixSymmetric<Treal, Tmatrix>,
00145        VectorGeneral<Treal, Tvector> >& smv);
00147     template<typename Tmatrix>
00148       VectorGeneral<Treal, Tvector>& operator=
00149       (const XYZpUV<Treal,
00150        MatrixSymmetric<Treal, Tmatrix>,
00151        VectorGeneral<Treal, Tvector>,
00152        Treal,
00153        VectorGeneral<Treal, Tvector> >& smvpsv);
00154 
00155     /* OPERATIONS INVOLVING TRIANGULAR MATRICES */
00157     template<typename Tmatrix>
00158       VectorGeneral<Treal, Tvector>& operator=
00159       (const XY<MatrixTriangular<Treal, Tmatrix>,
00160        VectorGeneral<Treal, Tvector> >& mv);
00161 
00162 
00163     /* LEVEL 1 operations */
00164     inline Treal eucl() const {
00165       return vectorPtr->eucl();
00166     }
00167 
00168     inline VectorGeneral<Treal, Tvector>& 
00169       operator*=(Treal const alpha) {
00170       *vectorPtr *= alpha;
00171       return *this;
00172     }
00173 
00174     inline VectorGeneral<Treal, Tvector>& 
00175       operator=(int const k) {
00176       *vectorPtr = k;
00177       return *this;
00178     }
00179 
00181     VectorGeneral<Treal, Tvector>& operator+=
00182       (const XY<Treal, VectorGeneral<Treal, Tvector> >& sv);
00183     
00184 
00185     inline Tvector const & getVector() const {return *vectorPtr;}
00186     
00187     std::string obj_type_id() const {return "VectorGeneral";}
00188   protected:
00189     ValidPtr<Tvector> vectorPtr;
00190 
00191     inline void writeToFileProt(std::ofstream & file) const {
00192       if (is_empty())
00193         // This means that the object's data structure has not been set
00194         return;
00195       vectorPtr->writeToFile(file);
00196     }
00197     inline void readFromFileProt(std::ifstream & file) {
00198       if (is_empty())
00199         // This means that the object's data structure has not been set
00200         return;
00201       vectorPtr->readFromFile(file);
00202     }
00203 
00204     inline void inMemorySet(bool inMem) {
00205       vectorPtr.inMemorySet(inMem);
00206     }
00207 
00208   private:
00209     
00210   }; /* end class VectorGeneral */
00211 
00212 
00213     /* LEVEL 2 operations */
00214     /* OPERATIONS INVOLVING ORDINARY MATRICES */
00216   template<typename Treal, typename Tvector>
00217     template<typename Tmatrix>
00218     VectorGeneral<Treal, Tvector>& 
00219     VectorGeneral<Treal, Tvector>::operator=
00220     (const XYZ<Treal,
00221      MatrixGeneral<Treal, Tmatrix>,
00222      VectorGeneral<Treal, Tvector> >& smv) {
00223     assert(!smv.tC);
00224     vectorPtr.haveDataStructureSet(true);
00225     if ( this == &smv.C ) {
00226       // We need a copy of the smv.C vector since it is the same as *this
00227       VectorGeneral<Treal, Tvector> tmp(smv.C);
00228       Tvector::gemv(smv.tB, smv.A, smv.B.getMatrix(),
00229                     *tmp.vectorPtr, 0, *this->vectorPtr);
00230     }
00231     else       
00232       Tvector::gemv(smv.tB, smv.A, smv.B.getMatrix(),
00233                     *smv.C.vectorPtr, 0, *this->vectorPtr);
00234     return *this;
00235   }
00236   
00238   template<typename Treal, typename Tvector>
00239     template<typename Tmatrix>
00240     VectorGeneral<Treal, Tvector>& 
00241     VectorGeneral<Treal, Tvector>::operator+=
00242     (const XYZ<Treal,
00243      MatrixGeneral<Treal, Tmatrix>,
00244      VectorGeneral<Treal, Tvector> >& smv) {
00245     assert(!smv.tC);
00246     assert(this != &smv.C);
00247     Tvector::gemv(smv.tB, smv.A, smv.B.getMatrix(),
00248                   *smv.C.vectorPtr, 1, *this->vectorPtr);
00249     return *this;
00250   }
00251   
00252   
00254   template<typename Treal, typename Tvector>
00255     template<typename Tmatrix>
00256     VectorGeneral<Treal, Tvector>& 
00257     VectorGeneral<Treal, Tvector>::operator=
00258     (const XYZpUV<Treal,
00259      MatrixGeneral<Treal, Tmatrix>,
00260      VectorGeneral<Treal, Tvector>,
00261      Treal,
00262      VectorGeneral<Treal, Tvector> >& smvpsv) {
00263     assert(!smvpsv.tC && !smvpsv.tE);
00264     assert(this != &smvpsv.C);
00265     if (this == &smvpsv.E)
00266       Tvector::gemv(smvpsv.tB, smvpsv.A, smvpsv.B.getMatrix(),
00267                     *smvpsv.C.vectorPtr, smvpsv.D, *this->vectorPtr);
00268     else
00269       throw Failure("VectorGeneral<Treal, Tvector>::operator="
00270                     "(const XYZpUV<Treal, "
00271                     "MatrixGeneral<Treal, Tmatrix>, "
00272                     "VectorGeneral<Treal, Tvector>, "
00273                     "Treal, "
00274                     "VectorGeneral<Treal, Tvector> >&) : "
00275                     "y = alpha * op(A) * x + beta * z "
00276                     "not supported for z != y");
00277     return *this;
00278   }
00279 
00280   
00281   /* OPERATIONS INVOLVING SYMMETRIC MATRICES */
00283   template<typename Treal, typename Tvector>
00284     template<typename Tmatrix>
00285     VectorGeneral<Treal, Tvector>& 
00286     VectorGeneral<Treal, Tvector>::operator=
00287     (const XYZ<Treal,
00288      MatrixSymmetric<Treal, Tmatrix>,
00289      VectorGeneral<Treal, Tvector> >& smv) {
00290     assert(!smv.tC);
00291     assert(this != &smv.C);
00292     vectorPtr.haveDataStructureSet(true);
00293     Tvector::symv('U', smv.A, smv.B.getMatrix(),
00294                   *smv.C.vectorPtr, 0, *this->vectorPtr);
00295     return *this;    
00296   }
00297   
00298 
00300   template<typename Treal, typename Tvector>
00301     template<typename Tmatrix>
00302     VectorGeneral<Treal, Tvector>& 
00303     VectorGeneral<Treal, Tvector>::operator+=
00304     (const XYZ<Treal,
00305      MatrixSymmetric<Treal, Tmatrix>,
00306      VectorGeneral<Treal, Tvector> >& smv) {
00307     assert(!smv.tC);
00308     assert(this != &smv.C);
00309     Tvector::symv('U', smv.A, smv.B.getMatrix(),
00310                   *smv.C.vectorPtr, 1, *this->vectorPtr);
00311     return *this;    
00312   }
00313 
00315   template<typename Treal, typename Tvector>
00316     template<typename Tmatrix>
00317     VectorGeneral<Treal, Tvector>& 
00318     VectorGeneral<Treal, Tvector>::operator=
00319     (const XYZpUV<Treal,
00320      MatrixSymmetric<Treal, Tmatrix>,
00321      VectorGeneral<Treal, Tvector>,
00322      Treal,
00323      VectorGeneral<Treal, Tvector> >& smvpsv) {
00324     assert(!smvpsv.tC && !smvpsv.tE);
00325     assert(this != &smvpsv.C);
00326     if (this == &smvpsv.E)
00327       Tvector::symv('U', smvpsv.A, smvpsv.B.getMatrix(),
00328                     *smvpsv.C.vectorPtr, smvpsv.D, *this->vectorPtr);
00329     else
00330       throw Failure("VectorGeneral<Treal, Tvector>::operator="
00331                     "(const XYZpUV<Treal, "
00332                     "MatrixSymmetric<Treal, Tmatrix>, "
00333                     "VectorGeneral<Treal, Tvector>, "
00334                     "Treal, "
00335                     "VectorGeneral<Treal, Tvector> >&) : "
00336                     "y = alpha * A * x + beta * z "
00337                     "not supported for z != y");
00338     return *this;    
00339   }
00340 
00341   /* OPERATIONS INVOLVING TRIANGULAR MATRICES */
00344   template<typename Treal, typename Tvector>
00345     template<typename Tmatrix>
00346     VectorGeneral<Treal, Tvector>& 
00347     VectorGeneral<Treal, Tvector>::operator=
00348     (const XY<MatrixTriangular<Treal, Tmatrix>,
00349      VectorGeneral<Treal, Tvector> >& mv) {
00350     assert(!mv.tB);
00351     if (this != &mv.B)
00352       throw Failure("y = A * x not supported for y != x ");
00353     Tvector::trmv('U', mv.tA, 
00354                   mv.A.getMatrix(), 
00355                   *this->vectorPtr);
00356     return *this;
00357   }
00358   
00359   /* LEVEL 1 operations */
00360 
00362   template<typename Treal, typename Tvector>
00363     VectorGeneral<Treal, Tvector>& 
00364     VectorGeneral<Treal, Tvector>::operator+=
00365     (const XY<Treal, VectorGeneral<Treal, Tvector> >& sv) {
00366     assert(!sv.tB);
00367     assert(this != &sv.B);
00368     Tvector::axpy(sv.A, *sv.B.vectorPtr, *this->vectorPtr);
00369     return *this;
00370   }
00371 
00372 
00373 
00374   /* Defined outside class */
00378   template<typename Treal, typename Tvector>
00379     Treal operator*(Xtrans<VectorGeneral<Treal, Tvector> > const & xT,
00380                     VectorGeneral<Treal, Tvector> const & y) {
00381     if (xT.tA == false)
00382       throw Failure("operator*("
00383                     "Xtrans<VectorGeneral<Treal, Tvector> > const &,"
00384                     " VectorGeneral<Treal, Tvector> const &): "
00385                     "Dimension mismatch in vector operation");
00386     return Tvector::dot(xT.A.getVector(), y.getVector());
00387   }
00388   
00389   
00390 
00391 }  /* end namespace mat */
00392 #endif
00393 

Generated on Wed Nov 21 09:32:38 2012 for ergo by  doxygen 1.4.7