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 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