ergo
MatrixBase.h
Go to the documentation of this file.
1 /* Ergo, version 3.2, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2012 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek.
4  *
5  * This program is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation, either version 3 of the License, or
8  * (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program. If not, see <http://www.gnu.org/licenses/>.
17  *
18  * Primary academic reference:
19  * Kohn−Sham Density Functional Theory Electronic Structure Calculations
20  * with Linearly Scaling Computational Time and Memory Usage,
21  * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek,
22  * J. Chem. Theory Comput. 7, 340 (2011),
23  * <http://dx.doi.org/10.1021/ct100611z>
24  *
25  * For further information about Ergo, see <http://www.ergoscf.org>.
26  */
27 
36 #ifndef MAT_MATRIXBASE
37 #define MAT_MATRIXBASE
38 #include <iostream>
39 #include <fstream>
40 #include <ios>
41 #include "FileWritable.h"
42 #include "matrix_proxy.h"
43 #include "ValidPtr.h"
44 #include "SizesAndBlocks.h"
45 namespace mat {
46  template<typename Treal, typename Tmatrix>
47  class MatrixGeneral;
48  template<typename Treal, typename Tmatrix>
50  template<typename Treal, typename Tmatrix>
52  template<typename Treal, typename Tvector>
55 
66  template<typename Treal, typename Tmatrix>
67  class MatrixBase : public FileWritable {
68  public:
69  friend class MatrixGeneral<Treal, Tmatrix>;
70  friend class MatrixSymmetric<Treal, Tmatrix>;
71  friend class MatrixTriangular<Treal, Tmatrix>;
72 
73 
74  inline void resetSizesAndBlocks(SizesAndBlocks const & newRows,
75  SizesAndBlocks const & newCols) {
77  matrixPtr->resetRows(newRows);
78  matrixPtr->resetCols(newCols);
79  }
80  inline void getRows(SizesAndBlocks & rowsCopy) const {
81  matrixPtr->getRows(rowsCopy);
82  }
83  inline void getCols(SizesAndBlocks & colsCopy) const {
84  matrixPtr->getCols(colsCopy);
85  }
86 
91  inline bool is_empty() const {
93  }
94 
95  inline Treal trace() const {
96  return matrixPtr->trace();
97  }
98 
99  inline void add_identity(Treal alpha) {
100  matrixPtr->addIdentity(alpha);
101  }
102  inline MatrixBase<Treal, Tmatrix>& operator*=(Treal const alpha) {
103  *matrixPtr *= alpha;
104  return *this;
105  }
106 
107  inline bool operator==(int k) const {
108  if (k == 0)
109  return *matrixPtr == 0;
110  else
111  throw Failure("MatrixBase::operator== only implemented for k == 0");
112  }
113 
114 
115 
116  inline void clear() {
117  if (is_empty())
118  // This means that the object's data structure has not been set
119  // There is nothing to clear and the matrixPtr is not valid either
120  return;
121  matrixPtr->clear();
122  }
123 
124  inline size_t memory_usage() const {
125  return matrixPtr->memory_usage();
126  }
127 
128  inline void write_to_buffer_count(int& n_bytes) const {
129  int ib_length = 3;
130  int vb_length = 0;
131  this->matrixPtr->write_to_buffer_count(ib_length, vb_length);
132  n_bytes = vb_length * sizeof(Treal) + ib_length * sizeof(int);
133  }
134 
135 #if 1
136  inline int get_nrows() const {
137  return matrixPtr->nScalarsRows();
138  }
139  inline int get_ncols() const {
140  return matrixPtr->nScalarsCols();
141  }
142 #endif
143 
144  inline Tmatrix const & getMatrix() const {return *matrixPtr;}
145  inline Tmatrix & getMatrix() {return *matrixPtr;}
146 
148  inline Treal maxAbsValue() const {return matrixPtr->maxAbsValue();}
149 
150  protected:
152 
153  MatrixBase():matrixPtr(new Tmatrix) {}
155  :FileWritable(other), matrixPtr(new Tmatrix) {
156  matrixPtr.haveDataStructureSet(other.matrixPtr.haveDataStructureGet());
157  /* getConstRefForCopying() is used here to make sure it works
158  also in the case when the matrix is written to file. */
159  *matrixPtr = other.matrixPtr.getConstRefForCopying();
160  matrixPtr.inMemorySet(other.matrixPtr.inMemoryGet());
161  }
162 
165  FileWritable::operator=(other); /* Allows us to copy mat on file */
166  matrixPtr.haveDataStructureSet(other.matrixPtr.haveDataStructureGet());
167  /* getConstRefForCopying() is used here to make sure it works
168  also in the case when the matrix is written to file. */
169  *matrixPtr = other.matrixPtr.getConstRefForCopying();
170  matrixPtr.inMemorySet(other.matrixPtr.inMemoryGet());
171  return *this;
172  }
173 
176  if (mt.A.matrixPtr.haveDataStructureGet()) {
178  }
179  if (mt.tA)
180  Tmatrix::transpose(*mt.A.matrixPtr, *this->matrixPtr);
181  else
182  *this->matrixPtr = *mt.A.matrixPtr;
183  return *this;
184  // FileWritable::operator=(other);/*Could be used to copy mat on file*/
185  }
186 
187 
188  void write_to_buffer_base(void* buffer, const int n_bytes,
189  const matrix_type mattype) const;
190  void read_from_buffer_base(void* buffer, const int n_bytes,
191  const matrix_type mattype);
192 
193  void writeToFileBase(std::ofstream & file,
194  matrix_type const mattype) const;
195  void readFromFileBase(std::ifstream & file,
196  matrix_type const mattype);
197 
198  std::string obj_type_id() const {return "MatrixBase";}
199  inline void inMemorySet(bool inMem) {
200  matrixPtr.inMemorySet(inMem);
201  }
202 
203  static void getPermutedIndexes(std::vector<int> const & index,
204  std::vector<int> const & permutation,
205  std::vector<int> & newIndex) {
206  newIndex.resize(index.size());
207  for (unsigned int i = 0; i < index.size(); ++i)
208  newIndex[i] = permutation[index[i]];
209  }
210 
211 
212  private:
213 
214  };
215 
216 
217  template<typename Treal, typename Tmatrix>
219  writeToFileBase(std::ofstream & file,
220  matrix_type const mattype) const {
221  int type = (int)mattype;
222  file.write((char*)&type,sizeof(int));
223 
224  if (is_empty())
225  // This means that the object's data structure has not been set
226  // The ValidPtr prevents setting the data structure between
227  // calls to writeToFile and readFromFile
228  return;
229  matrixPtr->writeToFile(file);
230  }
231 
232  template<typename Treal, typename Tmatrix>
234  readFromFileBase(std::ifstream & file,
235  matrix_type const mattype) {
236  char type[sizeof(int)];
237  file.read(type, sizeof(int));
238  if (((int)*type) != mattype)
239  throw Failure("MatrixBase<Treal, Tmatrix>::"
240  "readFromFile(std::ifstream &, "
241  "matrix_type const): Wrong matrix type");
242  if (is_empty())
243  // This means that the object's data structure has not been set
244  return;
245  matrixPtr->readFromFile(file);
246  }
247 
248 
249 
250  template<typename Treal, typename Tmatrix>
252  write_to_buffer_base(void* buffer, const int n_bytes,
253  const matrix_type mattype) const {
254  int ib_length = 3; /* Length of integer buffer, at least 3: matrix_type, */
255  /* ib_length and vb_length */
256  int vb_length = 0; /* Length of value buffer */
257  this->matrixPtr->write_to_buffer_count(ib_length, vb_length);
258  if (n_bytes >=
259  (int)(vb_length * sizeof(Treal) + ib_length * sizeof(int))) {
260  int* int_buf = (int*)buffer;
261  int_buf[0] = mattype;
262  int_buf[1] = ib_length;
263  int_buf[2] = vb_length;
264  Treal* value_buf = (Treal*)&(int_buf[ib_length]); /* Value buffer */
265  /* begins after integer buffer end */
266  int ib_index = 0;
267  int vb_index = 0;
268  this->matrixPtr->write_to_buffer(&int_buf[3], ib_length - 3,
269  value_buf, vb_length,
270  ib_index, vb_index);
271  }
272  else {
273  throw Failure("MatrixBase::write_to_buffer: Buffer is too small");
274  }
275  }
276 
277  template<typename Treal, typename Tmatrix>
279  read_from_buffer_base(void* buffer, const int n_bytes,
280  const matrix_type mattype) {
281  int* int_buf = (int*)buffer;
282  if(int_buf[0] == mattype) {
283  int ib_length = int_buf[1];
284  int vb_length = int_buf[2];
285  int ib_index = 0;
286  int vb_index = 0;
287  Treal* value_buf = (Treal*)&(int_buf[ib_length]);
288  this->matrixPtr->read_from_buffer(&int_buf[3], ib_length - 3,
289  value_buf, vb_length,
290  ib_index, vb_index);
291  }
292  else {
293  throw Failure("MatrixBase::read_from_buffer: Wrong matrix type");
294  }
295  }
296 
297 
298 } /* end namespace mat */
299 #endif
300 
301