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_MATRIXBASE 00037 #define MAT_MATRIXBASE 00038 #include <iostream> 00039 #include <fstream> 00040 #include <ios> 00041 #include "FileWritable.h" 00042 #include "matrix_proxy.h" 00043 #include "ValidPtr.h" 00044 #include "SizesAndBlocks.h" 00045 namespace mat { 00046 template<typename Treal, typename Tmatrix> 00047 class MatrixGeneral; 00048 template<typename Treal, typename Tmatrix> 00049 class MatrixSymmetric; 00050 template<typename Treal, typename Tmatrix> 00051 class MatrixTriangular; 00052 template<typename Treal, typename Tvector> 00053 class VectorGeneral; 00054 enum matrix_type {matrix_matr, matrix_symm, matrix_triang}; 00055 00066 template<typename Treal, typename Tmatrix> 00067 class MatrixBase : public FileWritable { 00068 public: 00069 friend class MatrixGeneral<Treal, Tmatrix>; 00070 friend class MatrixSymmetric<Treal, Tmatrix>; 00071 friend class MatrixTriangular<Treal, Tmatrix>; 00072 00073 00074 inline void resetSizesAndBlocks(SizesAndBlocks const & newRows, 00075 SizesAndBlocks const & newCols) { 00076 matrixPtr.haveDataStructureSet(true); 00077 matrixPtr->resetRows(newRows); 00078 matrixPtr->resetCols(newCols); 00079 } 00080 inline void getRows(SizesAndBlocks & rowsCopy) const { 00081 matrixPtr->getRows(rowsCopy); 00082 } 00083 inline void getCols(SizesAndBlocks & colsCopy) const { 00084 matrixPtr->getCols(colsCopy); 00085 } 00086 00091 inline bool is_empty() const { 00092 return !matrixPtr.haveDataStructureGet(); 00093 } 00094 00095 inline Treal trace() const { 00096 return matrixPtr->trace(); 00097 } 00098 00099 inline void add_identity(Treal alpha) { 00100 matrixPtr->addIdentity(alpha); 00101 } 00102 inline MatrixBase<Treal, Tmatrix>& operator*=(Treal const alpha) { 00103 *matrixPtr *= alpha; 00104 return *this; 00105 } 00106 00107 inline bool operator==(int k) const { 00108 if (k == 0) 00109 return *matrixPtr == 0; 00110 else 00111 throw Failure("MatrixBase::operator== only implemented for k == 0"); 00112 } 00113 00114 00115 00116 inline void clear() { 00117 if (is_empty()) 00118 // This means that the object's data structure has not been set 00119 // There is nothing to clear and the matrixPtr is not valid either 00120 return; 00121 matrixPtr->clear(); 00122 } 00123 00124 inline size_t memory_usage() const { 00125 return matrixPtr->memory_usage(); 00126 } 00127 00128 inline void write_to_buffer_count(int& n_bytes) const { 00129 int ib_length = 3; 00130 int vb_length = 0; 00131 this->matrixPtr->write_to_buffer_count(ib_length, vb_length); 00132 n_bytes = vb_length * sizeof(Treal) + ib_length * sizeof(int); 00133 } 00134 00135 #if 1 00136 inline int get_nrows() const { 00137 return matrixPtr->nScalarsRows(); 00138 } 00139 inline int get_ncols() const { 00140 return matrixPtr->nScalarsCols(); 00141 } 00142 #endif 00143 00144 inline Tmatrix const & getMatrix() const {return *matrixPtr;} 00145 inline Tmatrix & getMatrix() {return *matrixPtr;} 00146 00148 inline Treal maxAbsValue() const {return matrixPtr->maxAbsValue();} 00149 00150 protected: 00151 ValidPtr<Tmatrix> matrixPtr; 00152 00153 MatrixBase():matrixPtr(new Tmatrix) {} 00154 MatrixBase(const MatrixBase<Treal, Tmatrix>& other) 00155 :FileWritable(other), matrixPtr(new Tmatrix) { 00156 matrixPtr.haveDataStructureSet(other.matrixPtr.haveDataStructureGet()); 00157 /* getConstRefForCopying() is used here to make sure it works 00158 also in the case when the matrix is written to file. */ 00159 *matrixPtr = other.matrixPtr.getConstRefForCopying(); 00160 matrixPtr.inMemorySet(other.matrixPtr.inMemoryGet()); 00161 } 00162 00163 MatrixBase<Treal, Tmatrix>& 00164 operator=(const MatrixBase<Treal, Tmatrix>& other) { 00165 FileWritable::operator=(other); /* Allows us to copy mat on file */ 00166 matrixPtr.haveDataStructureSet(other.matrixPtr.haveDataStructureGet()); 00167 /* getConstRefForCopying() is used here to make sure it works 00168 also in the case when the matrix is written to file. */ 00169 *matrixPtr = other.matrixPtr.getConstRefForCopying(); 00170 matrixPtr.inMemorySet(other.matrixPtr.inMemoryGet()); 00171 return *this; 00172 } 00173 00174 MatrixBase<Treal, Tmatrix>& 00175 operator=(const Xtrans<MatrixGeneral<Treal, Tmatrix> >& mt) { 00176 if (mt.A.matrixPtr.haveDataStructureGet()) { 00177 matrixPtr.haveDataStructureSet(true); 00178 } 00179 if (mt.tA) 00180 Tmatrix::transpose(*mt.A.matrixPtr, *this->matrixPtr); 00181 else 00182 *this->matrixPtr = *mt.A.matrixPtr; 00183 return *this; 00184 // FileWritable::operator=(other);/*Could be used to copy mat on file*/ 00185 } 00186 00187 00188 void write_to_buffer_base(void* buffer, const int n_bytes, 00189 const matrix_type mattype) const; 00190 void read_from_buffer_base(void* buffer, const int n_bytes, 00191 const matrix_type mattype); 00192 00193 void writeToFileBase(std::ofstream & file, 00194 matrix_type const mattype) const; 00195 void readFromFileBase(std::ifstream & file, 00196 matrix_type const mattype); 00197 00198 std::string obj_type_id() const {return "MatrixBase";} 00199 inline void inMemorySet(bool inMem) { 00200 matrixPtr.inMemorySet(inMem); 00201 } 00202 00203 static void getPermutedIndexes(std::vector<int> const & index, 00204 std::vector<int> const & permutation, 00205 std::vector<int> & newIndex) { 00206 newIndex.resize(index.size()); 00207 for (unsigned int i = 0; i < index.size(); ++i) 00208 newIndex[i] = permutation[index[i]]; 00209 } 00210 00211 00212 private: 00213 00214 }; 00215 00216 00217 template<typename Treal, typename Tmatrix> 00218 void MatrixBase<Treal, Tmatrix>:: 00219 writeToFileBase(std::ofstream & file, 00220 matrix_type const mattype) const { 00221 int type = (int)mattype; 00222 file.write((char*)&type,sizeof(int)); 00223 00224 if (is_empty()) 00225 // This means that the object's data structure has not been set 00226 // The ValidPtr prevents setting the data structure between 00227 // calls to writeToFile and readFromFile 00228 return; 00229 matrixPtr->writeToFile(file); 00230 } 00231 00232 template<typename Treal, typename Tmatrix> 00233 void MatrixBase<Treal, Tmatrix>:: 00234 readFromFileBase(std::ifstream & file, 00235 matrix_type const mattype) { 00236 char type[sizeof(int)]; 00237 file.read(type, sizeof(int)); 00238 if (((int)*type) != mattype) 00239 throw Failure("MatrixBase<Treal, Tmatrix>::" 00240 "readFromFile(std::ifstream &, " 00241 "matrix_type const): Wrong matrix type"); 00242 if (is_empty()) 00243 // This means that the object's data structure has not been set 00244 return; 00245 matrixPtr->readFromFile(file); 00246 } 00247 00248 00249 00250 template<typename Treal, typename Tmatrix> 00251 void MatrixBase<Treal, Tmatrix>:: 00252 write_to_buffer_base(void* buffer, const int n_bytes, 00253 const matrix_type mattype) const { 00254 int ib_length = 3; /* Length of integer buffer, at least 3: matrix_type, */ 00255 /* ib_length and vb_length */ 00256 int vb_length = 0; /* Length of value buffer */ 00257 this->matrixPtr->write_to_buffer_count(ib_length, vb_length); 00258 if (n_bytes >= 00259 (int)(vb_length * sizeof(Treal) + ib_length * sizeof(int))) { 00260 int* int_buf = (int*)buffer; 00261 int_buf[0] = mattype; 00262 int_buf[1] = ib_length; 00263 int_buf[2] = vb_length; 00264 Treal* value_buf = (Treal*)&(int_buf[ib_length]); /* Value buffer */ 00265 /* begins after integer buffer end */ 00266 int ib_index = 0; 00267 int vb_index = 0; 00268 this->matrixPtr->write_to_buffer(&int_buf[3], ib_length - 3, 00269 value_buf, vb_length, 00270 ib_index, vb_index); 00271 } 00272 else { 00273 throw Failure("MatrixBase::write_to_buffer: Buffer is too small"); 00274 } 00275 } 00276 00277 template<typename Treal, typename Tmatrix> 00278 void MatrixBase<Treal, Tmatrix>:: 00279 read_from_buffer_base(void* buffer, const int n_bytes, 00280 const matrix_type mattype) { 00281 int* int_buf = (int*)buffer; 00282 if(int_buf[0] == mattype) { 00283 int ib_length = int_buf[1]; 00284 int vb_length = int_buf[2]; 00285 int ib_index = 0; 00286 int vb_index = 0; 00287 Treal* value_buf = (Treal*)&(int_buf[ib_length]); 00288 this->matrixPtr->read_from_buffer(&int_buf[3], ib_length - 3, 00289 value_buf, vb_length, 00290 ib_index, vb_index); 00291 } 00292 else { 00293 throw Failure("MatrixBase::read_from_buffer: Wrong matrix type"); 00294 } 00295 } 00296 00297 00298 } /* end namespace mat */ 00299 #endif 00300 00301