ergo
VectorGeneral.h
Go to the documentation of this file.
1 /* Ergo, version 3.7, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2018 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek,
4  * and Anastasia Kruchinina.
5  *
6  * This program is free software: you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation, either version 3 of the License, or
9  * (at your option) any later version.
10  *
11  * This program is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with this program. If not, see <http://www.gnu.org/licenses/>.
18  *
19  * Primary academic reference:
20  * Ergo: An open-source program for linear-scaling electronic structure
21  * calculations,
22  * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia
23  * Kruchinina,
24  * SoftwareX 7, 107 (2018),
25  * <http://dx.doi.org/10.1016/j.softx.2018.03.005>
26  *
27  * For further information about Ergo, see <http://www.ergoscf.org>.
28  */
29 
38 #ifndef MAT_VECTORGENERAL
39 #define MAT_VECTORGENERAL
40 #include <iostream>
41 #include <fstream>
42 #include <ios>
43 #include "FileWritable.h"
44 #include "matrix_proxy.h"
45 #include "ValidPtr.h"
46 namespace mat {
47  template<typename Treal, typename Tvector>
48  class VectorGeneral : public FileWritable {
49  public:
50 
51  inline void resetSizesAndBlocks(SizesAndBlocks const & newRows) {
53  vectorPtr->resetRows(newRows);
54  }
55 
56  inline bool is_empty() const {
58  }
59 
60  inline void clear_structure(){
62  }
63 
64  VectorGeneral(SizesAndBlocks const & newRows):vectorPtr(new Tvector) {
65  resetSizesAndBlocks(newRows);
66  }
67  VectorGeneral():vectorPtr(new Tvector) {}
69  :FileWritable(other), vectorPtr(new Tvector) {
70  if (other.vectorPtr.haveDataStructureGet()) {
72  }
73  *vectorPtr = *other.vectorPtr;
74  }
75 
76  inline void assign_from_full
77  (std::vector<Treal> const & fullVector,
78  SizesAndBlocks const & newRows) {
79  resetSizesAndBlocks(newRows);
80  this->vectorPtr->assignFromFull(fullVector);
81  }
82  inline void fullvector(std::vector<Treal> & fullVector) const {
83  this->vectorPtr->fullVector(fullVector);
84  }
87  if (other.vectorPtr.haveDataStructureGet()) {
89  }
90  *this->vectorPtr = *other.vectorPtr;
91  return *this;
92  }
93 
94  inline void clear() {
95  if (is_empty())
96  // This means that the object's data structure has not been set
97  // There is nothing to clear and the vectorPtr is not valid either
98  return;
99  vectorPtr->clear();
100  }
101 
102  inline void rand() {
103  vectorPtr->randomNormalized();
104  }
105 
106  /* LEVEL 2 operations */
107  /* OPERATIONS INVOLVING ORDINARY MATRICES */
109  template<typename Tmatrix>
111  (const XYZ<Treal,
114 
116  template<typename Tmatrix>
118  (const XYZ<Treal,
122  template<typename Tmatrix>
124  (const XYZpUV<Treal,
127  Treal,
129 
131  template<typename Tmatrix>
135  Treal ONE = 1.0;
136  return this->operator=(XYZ<Treal, MatrixGeneral<Treal, Tmatrix>,
137  VectorGeneral<Treal, Tvector> >(ONE, mv.A, mv.B,
138  false, mv.tA, mv.tB));
139  }
140 
141  /* OPERATIONS INVOLVING SYMMETRIC MATRICES */
143  template<typename Tmatrix>
145  (const XYZ<Treal,
149  template<typename Tmatrix>
151  (const XYZ<Treal,
155  template<typename Tmatrix>
157  (const XYZpUV<Treal,
160  Treal,
162 
163  /* OPERATIONS INVOLVING TRIANGULAR MATRICES */
165  template<typename Tmatrix>
169 
170 
171  /* LEVEL 1 operations */
172  inline Treal eucl() const {
173  return vectorPtr->eucl();
174  }
175 
177  operator*=(Treal const alpha) {
178  *vectorPtr *= alpha;
179  return *this;
180  }
181 
183  operator=(int const k) {
184  *vectorPtr = k;
185  return *this;
186  }
187 
191 
192 
193  inline Tvector const & getVector() const {return *vectorPtr;}
194 
195  std::string obj_type_id() const {return "VectorGeneral";}
196  protected:
198 
199  inline void writeToFileProt(std::ofstream & file) const {
200  if (is_empty())
201  // This means that the object's data structure has not been set
202  return;
203  vectorPtr->writeToFile(file);
204  }
205  inline void readFromFileProt(std::ifstream & file) {
206  if (is_empty())
207  // This means that the object's data structure has not been set
208  return;
209  vectorPtr->readFromFile(file);
210  }
211 
212  inline void inMemorySet(bool inMem) {
213  vectorPtr.inMemorySet(inMem);
214  }
215 
216  private:
217 
218  }; /* end class VectorGeneral */
219 
220 
221  /* LEVEL 2 operations */
222  /* OPERATIONS INVOLVING ORDINARY MATRICES */
224  template<typename Treal, typename Tvector>
225  template<typename Tmatrix>
226  VectorGeneral<Treal, Tvector>&
227  VectorGeneral<Treal, Tvector>::operator=
228  (const XYZ<Treal,
231  assert(!smv.tC);
232  vectorPtr.haveDataStructureSet(true);
233  if ( this == &smv.C ) {
234  // We need a copy of the smv.C vector since it is the same as *this
236  Tvector::gemv(smv.tB, smv.A, smv.B.getMatrix(),
237  *tmp.vectorPtr, 0, *this->vectorPtr);
238  }
239  else
240  Tvector::gemv(smv.tB, smv.A, smv.B.getMatrix(),
241  *smv.C.vectorPtr, 0, *this->vectorPtr);
242  return *this;
243  }
244 
246  template<typename Treal, typename Tvector>
247  template<typename Tmatrix>
250  (const XYZ<Treal,
253  assert(!smv.tC);
254  assert(this != &smv.C);
255  Tvector::gemv(smv.tB, smv.A, smv.B.getMatrix(),
256  *smv.C.vectorPtr, 1, *this->vectorPtr);
257  return *this;
258  }
259 
260 
262  template<typename Treal, typename Tvector>
263  template<typename Tmatrix>
266  (const XYZpUV<Treal,
269  Treal,
270  VectorGeneral<Treal, Tvector> >& smvpsv) {
271  assert(!smvpsv.tC && !smvpsv.tE);
272  assert(this != &smvpsv.C);
273  if (this == &smvpsv.E)
274  Tvector::gemv(smvpsv.tB, smvpsv.A, smvpsv.B.getMatrix(),
275  *smvpsv.C.vectorPtr, smvpsv.D, *this->vectorPtr);
276  else
277  throw Failure("VectorGeneral<Treal, Tvector>::operator="
278  "(const XYZpUV<Treal, "
279  "MatrixGeneral<Treal, Tmatrix>, "
280  "VectorGeneral<Treal, Tvector>, "
281  "Treal, "
282  "VectorGeneral<Treal, Tvector> >&) : "
283  "y = alpha * op(A) * x + beta * z "
284  "not supported for z != y");
285  return *this;
286  }
287 
288 
289  /* OPERATIONS INVOLVING SYMMETRIC MATRICES */
291  template<typename Treal, typename Tvector>
292  template<typename Tmatrix>
295  (const XYZ<Treal,
298  assert(!smv.tC);
299  assert(this != &smv.C);
300  vectorPtr.haveDataStructureSet(true);
301  Tvector::symv('U', smv.A, smv.B.getMatrix(),
302  *smv.C.vectorPtr, 0, *this->vectorPtr);
303  return *this;
304  }
305 
306 
308  template<typename Treal, typename Tvector>
309  template<typename Tmatrix>
312  (const XYZ<Treal,
315  assert(!smv.tC);
316  assert(this != &smv.C);
317  Tvector::symv('U', smv.A, smv.B.getMatrix(),
318  *smv.C.vectorPtr, 1, *this->vectorPtr);
319  return *this;
320  }
321 
323  template<typename Treal, typename Tvector>
324  template<typename Tmatrix>
327  (const XYZpUV<Treal,
330  Treal,
331  VectorGeneral<Treal, Tvector> >& smvpsv) {
332  assert(!smvpsv.tC && !smvpsv.tE);
333  assert(this != &smvpsv.C);
334  if (this == &smvpsv.E)
335  Tvector::symv('U', smvpsv.A, smvpsv.B.getMatrix(),
336  *smvpsv.C.vectorPtr, smvpsv.D, *this->vectorPtr);
337  else
338  throw Failure("VectorGeneral<Treal, Tvector>::operator="
339  "(const XYZpUV<Treal, "
340  "MatrixSymmetric<Treal, Tmatrix>, "
341  "VectorGeneral<Treal, Tvector>, "
342  "Treal, "
343  "VectorGeneral<Treal, Tvector> >&) : "
344  "y = alpha * A * x + beta * z "
345  "not supported for z != y");
346  return *this;
347  }
348 
349  /* OPERATIONS INVOLVING TRIANGULAR MATRICES */
352  template<typename Treal, typename Tvector>
353  template<typename Tmatrix>
358  assert(!mv.tB);
359  if (this != &mv.B)
360  throw Failure("y = A * x not supported for y != x ");
361  Tvector::trmv('U', mv.tA,
362  mv.A.getMatrix(),
363  *this->vectorPtr);
364  return *this;
365  }
366 
367  /* LEVEL 1 operations */
368 
370  template<typename Treal, typename Tvector>
374  assert(!sv.tB);
375  assert(this != &sv.B);
376  Tvector::axpy(sv.A, *sv.B.vectorPtr, *this->vectorPtr);
377  return *this;
378  }
379 
380 
381 
382  /* Defined outside class */
386  template<typename Treal, typename Tvector>
388  VectorGeneral<Treal, Tvector> const & y) {
389  if (xT.tA == false)
390  throw Failure("operator*("
391  "Xtrans<VectorGeneral<Treal, Tvector> > const &,"
392  " VectorGeneral<Treal, Tvector> const &): "
393  "Dimension mismatch in vector operation");
394  return Tvector::dot(xT.A.getVector(), y.getVector());
395  }
396 
397 
398 
399 } /* end namespace mat */
400 #endif
401 
void assign_from_full(std::vector< Treal > const &fullVector, SizesAndBlocks const &newRows)
Definition: VectorGeneral.h:77
VectorGeneral()
Definition: VectorGeneral.h:67
Normal matrix.
Definition: MatrixBase.h:49
std::string obj_type_id() const
Definition: VectorGeneral.h:195
ValidPtr< Tvector > vectorPtr
Definition: VectorGeneral.h:197
Proxy structs used by the matrix API.
void fullvector(std::vector< Treal > &fullVector) const
Definition: VectorGeneral.h:82
void readFromFileProt(std::ifstream &file)
Read object from file.
Definition: VectorGeneral.h:205
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:409
Definition: MatrixBase.h:55
Definition: allocate.cc:39
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:157
Describes dimensions of matrix and its blocks on all levels.
Definition: SizesAndBlocks.h:45
void clear_structure()
Definition: VectorGeneral.h:60
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:212
static T dot(const int *n, const T *dx, const int *incx, const T *dy, const int *incy)
Definition: mat_gblas.h:425
VectorGeneral< Treal, Tvector > & operator=(int const k)
Definition: VectorGeneral.h:183
Smart pointer class to control access to object.
bool haveDataStructureGet() const
Definition: ValidPtr.h:102
This proxy expresses the result of multiplication of three objects, of possibly different types...
Definition: matrix_proxy.h:67
This proxy expresses the result of transposition of an object of type TX.
Definition: matrix_proxy.h:118
bool is_empty() const
Definition: VectorGeneral.h:56
static void axpy(const int *n, const T *da, const T *dx, const int *incx, T *dy, const int *incy)
Definition: mat_gblas.h:431
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:400
void rand()
Definition: VectorGeneral.h:102
void clear()
Release memory for the information written to file.
Definition: VectorGeneral.h:94
void resetSizesAndBlocks(SizesAndBlocks const &newRows)
Definition: VectorGeneral.h:51
Write and read objects to/from file.
Definition: FileWritable.h:56
VectorGeneral< Treal, Tvector > & operator=(const VectorGeneral< Treal, Tvector > &other)
Definition: VectorGeneral.h:86
VectorGeneral(const VectorGeneral< Treal, Tvector > &other)
Definition: VectorGeneral.h:68
VectorGeneral(SizesAndBlocks const &newRows)
Definition: VectorGeneral.h:64
Treal eucl() const
Definition: VectorGeneral.h:172
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:391
VectorGeneral< Treal, Tvector > & operator*=(Treal const alpha)
Definition: VectorGeneral.h:177
Tvector const & getVector() const
Definition: VectorGeneral.h:193
Definition: Failure.h:57
This proxy expresses the result of multiplication of two objects, of possibly different types...
Definition: matrix_proxy.h:51
void haveDataStructureSet(bool val)
Definition: ValidPtr.h:99
Abstract class for simple writing and reading of objects to/from file.
void writeToFileProt(std::ofstream &file) const
Write object to file.
Definition: VectorGeneral.h:199
Symmetric matrix.
Definition: MatrixBase.h:51
This proxy expresses the result of multiplication of three objects added to two other multiplied obje...
Definition: matrix_proxy.h:88
void inMemorySet(bool val)
Definition: ValidPtr.h:93