ergo
VectorGeneral.h
Go to the documentation of this file.
1 /* Ergo, version 3.4, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2014 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_VECTORGENERAL
37 #define MAT_VECTORGENERAL
38 #include <iostream>
39 #include <fstream>
40 #include <ios>
41 #include "FileWritable.h"
42 #include "matrix_proxy.h"
43 #include "ValidPtr.h"
44 namespace mat {
45  template<typename Treal, typename Tvector>
46  class VectorGeneral : public FileWritable {
47  public:
48 
49  inline void resetSizesAndBlocks(SizesAndBlocks const & newRows) {
51  vectorPtr->resetRows(newRows);
52  }
53 
54  inline bool is_empty() const {
56  }
57 
58 
59  VectorGeneral():vectorPtr(new Tvector) {}
61  :FileWritable(other), vectorPtr(new Tvector) {
62  if (other.vectorPtr.haveDataStructureGet()) {
64  }
65  *vectorPtr = *other.vectorPtr;
66  }
67 
68  inline void assign_from_full
69  (std::vector<Treal> const & fullVector,
70  SizesAndBlocks const & newRows) {
71  resetSizesAndBlocks(newRows);
72  this->vectorPtr->assignFromFull(fullVector);
73  }
74  inline void fullvector(std::vector<Treal> & fullVector) const {
75  this->vectorPtr->fullVector(fullVector);
76  }
79  if (other.vectorPtr.haveDataStructureGet()) {
81  }
82  *this->vectorPtr = *other.vectorPtr;
83  return *this;
84  }
85 
86  inline void clear() {
87  if (is_empty())
88  // This means that the object's data structure has not been set
89  // There is nothing to clear and the vectorPtr is not valid either
90  return;
91  vectorPtr->clear();
92  }
93 
94  inline void rand() {
95  vectorPtr->randomNormalized();
96  }
97 
98  /* LEVEL 2 operations */
99  /* OPERATIONS INVOLVING ORDINARY MATRICES */
101  template<typename Tmatrix>
103  (const XYZ<Treal,
106 
108  template<typename Tmatrix>
110  (const XYZ<Treal,
111  MatrixGeneral<Treal, Tmatrix>,
114  template<typename Tmatrix>
116  (const XYZpUV<Treal,
117  MatrixGeneral<Treal, Tmatrix>,
119  Treal,
120  VectorGeneral<Treal, Tvector> >& smvpsv);
121 
123  template<typename Tmatrix>
124  VectorGeneral<Treal, Tvector>& operator=
126  VectorGeneral<Treal, Tvector> >& mv) {
127  Treal ONE = 1.0;
128  return this->operator=(XYZ<Treal, MatrixGeneral<Treal, Tmatrix>,
129  VectorGeneral<Treal, Tvector> >(ONE, mv.A, mv.B,
130  false, mv.tA, mv.tB));
131  }
132 
133  /* OPERATIONS INVOLVING SYMMETRIC MATRICES */
135  template<typename Tmatrix>
136  VectorGeneral<Treal, Tvector>& operator=
137  (const XYZ<Treal,
139  VectorGeneral<Treal, Tvector> >& smv);
141  template<typename Tmatrix>
142  VectorGeneral<Treal, Tvector>& operator+=
143  (const XYZ<Treal,
144  MatrixSymmetric<Treal, Tmatrix>,
145  VectorGeneral<Treal, Tvector> >& smv);
147  template<typename Tmatrix>
148  VectorGeneral<Treal, Tvector>& operator=
149  (const XYZpUV<Treal,
150  MatrixSymmetric<Treal, Tmatrix>,
151  VectorGeneral<Treal, Tvector>,
152  Treal,
153  VectorGeneral<Treal, Tvector> >& smvpsv);
154 
155  /* OPERATIONS INVOLVING TRIANGULAR MATRICES */
157  template<typename Tmatrix>
158  VectorGeneral<Treal, Tvector>& operator=
160  VectorGeneral<Treal, Tvector> >& mv);
161 
162 
163  /* LEVEL 1 operations */
164  inline Treal eucl() const {
165  return vectorPtr->eucl();
166  }
167 
168  inline VectorGeneral<Treal, Tvector>&
169  operator*=(Treal const alpha) {
170  *vectorPtr *= alpha;
171  return *this;
172  }
173 
174  inline VectorGeneral<Treal, Tvector>&
175  operator=(int const k) {
176  *vectorPtr = k;
177  return *this;
178  }
179 
181  VectorGeneral<Treal, Tvector>& operator+=
183 
184 
185  inline Tvector const & getVector() const {return *vectorPtr;}
186 
187  std::string obj_type_id() const {return "VectorGeneral";}
188  protected:
190 
191  inline void writeToFileProt(std::ofstream & file) const {
192  if (is_empty())
193  // This means that the object's data structure has not been set
194  return;
195  vectorPtr->writeToFile(file);
196  }
197  inline void readFromFileProt(std::ifstream & file) {
198  if (is_empty())
199  // This means that the object's data structure has not been set
200  return;
201  vectorPtr->readFromFile(file);
202  }
203 
204  inline void inMemorySet(bool inMem) {
205  vectorPtr.inMemorySet(inMem);
206  }
207 
208  private:
209 
210  }; /* end class VectorGeneral */
211 
212 
213  /* LEVEL 2 operations */
214  /* OPERATIONS INVOLVING ORDINARY MATRICES */
216  template<typename Treal, typename Tvector>
217  template<typename Tmatrix>
220  (const XYZ<Treal,
223  assert(!smv.tC);
225  if ( this == &smv.C ) {
226  // We need a copy of the smv.C vector since it is the same as *this
228  Tvector::gemv(smv.tB, smv.A, smv.B.getMatrix(),
229  *tmp.vectorPtr, 0, *this->vectorPtr);
230  }
231  else
232  Tvector::gemv(smv.tB, smv.A, smv.B.getMatrix(),
233  *smv.C.vectorPtr, 0, *this->vectorPtr);
234  return *this;
235  }
236 
238  template<typename Treal, typename Tvector>
239  template<typename Tmatrix>
241  VectorGeneral<Treal, Tvector>::operator+=
242  (const XYZ<Treal,
243  MatrixGeneral<Treal, Tmatrix>,
245  assert(!smv.tC);
246  assert(this != &smv.C);
247  Tvector::gemv(smv.tB, smv.A, smv.B.getMatrix(),
248  *smv.C.vectorPtr, 1, *this->vectorPtr);
249  return *this;
250  }
251 
252 
254  template<typename Treal, typename Tvector>
255  template<typename Tmatrix>
257  VectorGeneral<Treal, Tvector>::operator=
258  (const XYZpUV<Treal,
259  MatrixGeneral<Treal, Tmatrix>,
261  Treal,
262  VectorGeneral<Treal, Tvector> >& smvpsv) {
263  assert(!smvpsv.tC && !smvpsv.tE);
264  assert(this != &smvpsv.C);
265  if (this == &smvpsv.E)
266  Tvector::gemv(smvpsv.tB, smvpsv.A, smvpsv.B.getMatrix(),
267  *smvpsv.C.vectorPtr, smvpsv.D, *this->vectorPtr);
268  else
269  throw Failure("VectorGeneral<Treal, Tvector>::operator="
270  "(const XYZpUV<Treal, "
271  "MatrixGeneral<Treal, Tmatrix>, "
272  "VectorGeneral<Treal, Tvector>, "
273  "Treal, "
274  "VectorGeneral<Treal, Tvector> >&) : "
275  "y = alpha * op(A) * x + beta * z "
276  "not supported for z != y");
277  return *this;
278  }
279 
280 
281  /* OPERATIONS INVOLVING SYMMETRIC MATRICES */
283  template<typename Treal, typename Tvector>
284  template<typename Tmatrix>
285  VectorGeneral<Treal, Tvector>&
286  VectorGeneral<Treal, Tvector>::operator=
287  (const XYZ<Treal,
289  VectorGeneral<Treal, Tvector> >& smv) {
290  assert(!smv.tC);
291  assert(this != &smv.C);
293  Tvector::symv('U', smv.A, smv.B.getMatrix(),
294  *smv.C.vectorPtr, 0, *this->vectorPtr);
295  return *this;
296  }
297 
298 
300  template<typename Treal, typename Tvector>
301  template<typename Tmatrix>
302  VectorGeneral<Treal, Tvector>&
303  VectorGeneral<Treal, Tvector>::operator+=
304  (const XYZ<Treal,
305  MatrixSymmetric<Treal, Tmatrix>,
306  VectorGeneral<Treal, Tvector> >& smv) {
307  assert(!smv.tC);
308  assert(this != &smv.C);
309  Tvector::symv('U', smv.A, smv.B.getMatrix(),
310  *smv.C.vectorPtr, 1, *this->vectorPtr);
311  return *this;
312  }
313 
315  template<typename Treal, typename Tvector>
316  template<typename Tmatrix>
317  VectorGeneral<Treal, Tvector>&
318  VectorGeneral<Treal, Tvector>::operator=
319  (const XYZpUV<Treal,
320  MatrixSymmetric<Treal, Tmatrix>,
321  VectorGeneral<Treal, Tvector>,
322  Treal,
323  VectorGeneral<Treal, Tvector> >& smvpsv) {
324  assert(!smvpsv.tC && !smvpsv.tE);
325  assert(this != &smvpsv.C);
326  if (this == &smvpsv.E)
327  Tvector::symv('U', smvpsv.A, smvpsv.B.getMatrix(),
328  *smvpsv.C.vectorPtr, smvpsv.D, *this->vectorPtr);
329  else
330  throw Failure("VectorGeneral<Treal, Tvector>::operator="
331  "(const XYZpUV<Treal, "
332  "MatrixSymmetric<Treal, Tmatrix>, "
333  "VectorGeneral<Treal, Tvector>, "
334  "Treal, "
335  "VectorGeneral<Treal, Tvector> >&) : "
336  "y = alpha * A * x + beta * z "
337  "not supported for z != y");
338  return *this;
339  }
340 
341  /* OPERATIONS INVOLVING TRIANGULAR MATRICES */
344  template<typename Treal, typename Tvector>
345  template<typename Tmatrix>
346  VectorGeneral<Treal, Tvector>&
347  VectorGeneral<Treal, Tvector>::operator=
349  VectorGeneral<Treal, Tvector> >& mv) {
350  assert(!mv.tB);
351  if (this != &mv.B)
352  throw Failure("y = A * x not supported for y != x ");
353  Tvector::trmv('U', mv.tA,
354  mv.A.getMatrix(),
355  *this->vectorPtr);
356  return *this;
357  }
358 
359  /* LEVEL 1 operations */
360 
362  template<typename Treal, typename Tvector>
363  VectorGeneral<Treal, Tvector>&
364  VectorGeneral<Treal, Tvector>::operator+=
366  assert(!sv.tB);
367  assert(this != &sv.B);
368  Tvector::axpy(sv.A, *sv.B.vectorPtr, *this->vectorPtr);
369  return *this;
370  }
371 
372 
373 
374  /* Defined outside class */
378  template<typename Treal, typename Tvector>
379  Treal operator*(Xtrans<VectorGeneral<Treal, Tvector> > const & xT,
380  VectorGeneral<Treal, Tvector> const & y) {
381  if (xT.tA == false)
382  throw Failure("operator*("
383  "Xtrans<VectorGeneral<Treal, Tvector> > const &,"
384  " VectorGeneral<Treal, Tvector> const &): "
385  "Dimension mismatch in vector operation");
386  return Tvector::dot(xT.A.getVector(), y.getVector());
387  }
388 
389 
390 
391 } /* end namespace mat */
392 #endif
393 
void assign_from_full(std::vector< Treal > const &fullVector, SizesAndBlocks const &newRows)
Definition: VectorGeneral.h:69
VectorGeneral()
Definition: VectorGeneral.h:59
Normal matrix.
Definition: MatrixBase.h:47
ValidPtr< Tvector > vectorPtr
Definition: VectorGeneral.h:189
Proxy structs used by the matrix API.
void readFromFileProt(std::ifstream &file)
Read object from file.
Definition: VectorGeneral.h:197
static void trmv(const char *uplo, const char *trans, const char *diag, const int *n, const T *A, const int *lda, T *x, const int *incx)
Definition: mat_gblas.h:407
Definition: MatrixBase.h:53
bool is_empty() const
Definition: VectorGeneral.h:54
Definition: allocate.cc:30
XY< TX, TY > operator*(Xtrans< TX > const &trAA, Xtrans< TY > const &trBB)
Multiplication of two transposition proxys holding objects of type TX and TY respectively.
Definition: matrix_proxy.h:155
Describes dimensions of matrix and its blocks on all levels.
Definition: SizesAndBlocks.h:37
void inMemorySet(bool inMem)
Make object invalid (false) via this function when object is written to file and valid (true) when ob...
Definition: VectorGeneral.h:204
static T dot(const int *n, const T *dx, const int *incx, const T *dy, const int *incy)
Definition: mat_gblas.h:423
VectorGeneral< Treal, Tvector > & operator=(int const k)
Definition: VectorGeneral.h:175
Smart pointer class to control access to object.
Tvector const & getVector() const
Definition: VectorGeneral.h:185
This proxy expresses the result of multiplication of three objects, of possibly different types...
Definition: matrix_proxy.h:65
void writeToFileProt(std::ofstream &file) const
Write object to file.
Definition: VectorGeneral.h:191
This proxy expresses the result of transposition of an object of type TX.
Definition: matrix_proxy.h:116
static void axpy(const int *n, const T *da, const T *dx, const int *incx, T *dy, const int *incy)
Definition: mat_gblas.h:429
static void symv(const char *uplo, const int *n, const T *alpha, const T *A, const int *lda, const T *x, const int *incx, const T *beta, T *y, const int *incy)
Definition: mat_gblas.h:398
void rand()
Definition: VectorGeneral.h:94
void clear()
Release memory for the information written to file.
Definition: VectorGeneral.h:86
void resetSizesAndBlocks(SizesAndBlocks const &newRows)
Definition: VectorGeneral.h:49
Write and read objects to/from file.
Definition: FileWritable.h:54
VectorGeneral< Treal, Tvector > & operator=(const VectorGeneral< Treal, Tvector > &other)
Definition: VectorGeneral.h:78
VectorGeneral(const VectorGeneral< Treal, Tvector > &other)
Definition: VectorGeneral.h:60
static void gemv(const char *ta, const int *m, const int *n, const T *alpha, const T *A, const int *lda, const T *x, const int *incx, const T *beta, T *y, const int *incy)
Definition: mat_gblas.h:389
VectorGeneral< Treal, Tvector > & operator*=(Treal const alpha)
Definition: VectorGeneral.h:169
Treal eucl() const
Definition: VectorGeneral.h:164
Definition: Failure.h:47
This proxy expresses the result of multiplication of two objects, of possibly different types...
Definition: matrix_proxy.h:49
std::string obj_type_id() const
Definition: VectorGeneral.h:187
void haveDataStructureSet(bool val)
Definition: ValidPtr.h:97
void fullvector(std::vector< Treal > &fullVector) const
Definition: VectorGeneral.h:74
Abstract class for simple writing and reading of objects to/from file.
bool haveDataStructureGet() const
Definition: ValidPtr.h:100
Symmetric matrix.
Definition: MatrixBase.h:49
This proxy expresses the result of multiplication of three objects added to two other multiplied obje...
Definition: matrix_proxy.h:86
void inMemorySet(bool val)
Definition: ValidPtr.h:91