ergo
VectorGeneral.h
Go to the documentation of this file.
1 /* Ergo, version 3.8, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2019 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) {}
68 
69  /* In the code we are using std::vector<VectorGeneral> which in the c++ standard before c++11 requires move operation like T x_new = x which calls implicitly the copy constructor. To make it work with g++ versions without c++11 support we remove the keyword explicit. */
70 #if __cplusplus >= 201103L
71  explicit VectorGeneral(const VectorGeneral<Treal, Tvector>& other)
72 #else
74 #endif
75  :FileWritable(other), vectorPtr(new Tvector) {
76  if (other.vectorPtr.haveDataStructureGet()) {
78  }
79  *vectorPtr = *other.vectorPtr;
80  }
81 
82  inline void assign_from_full
83  (std::vector<Treal> const & fullVector,
84  SizesAndBlocks const & newRows) {
85  resetSizesAndBlocks(newRows);
86  this->vectorPtr->assignFromFull(fullVector);
87  }
88  inline void fullvector(std::vector<Treal> & fullVector) const {
89  this->vectorPtr->fullVector(fullVector);
90  }
93  if (other.vectorPtr.haveDataStructureGet()) {
95  }
96  *this->vectorPtr = *other.vectorPtr;
97  return *this;
98  }
99 
100  inline void clear() {
101  if (is_empty())
102  // This means that the object's data structure has not been set
103  // There is nothing to clear and the vectorPtr is not valid either
104  return;
105  vectorPtr->clear();
106  }
107 
108  inline void rand() {
109  vectorPtr->randomNormalized();
110  }
111 
112  /* LEVEL 2 operations */
113  /* OPERATIONS INVOLVING ORDINARY MATRICES */
115  template<typename Tmatrix>
117  (const XYZ<Treal,
120 
122  template<typename Tmatrix>
124  (const XYZ<Treal,
128  template<typename Tmatrix>
130  (const XYZpUV<Treal,
133  Treal,
135 
137  template<typename Tmatrix>
141  Treal ONE = 1.0;
142  return this->operator=(XYZ<Treal, MatrixGeneral<Treal, Tmatrix>,
143  VectorGeneral<Treal, Tvector> >(ONE, mv.A, mv.B,
144  false, mv.tA, mv.tB));
145  }
146 
147  /* OPERATIONS INVOLVING SYMMETRIC MATRICES */
149  template<typename Tmatrix>
151  (const XYZ<Treal,
155  template<typename Tmatrix>
157  (const XYZ<Treal,
161  template<typename Tmatrix>
163  (const XYZpUV<Treal,
166  Treal,
168 
169  /* OPERATIONS INVOLVING TRIANGULAR MATRICES */
171  template<typename Tmatrix>
175 
176 
177  /* LEVEL 1 operations */
178  inline Treal eucl() const {
179  return vectorPtr->eucl();
180  }
181 
183  operator*=(Treal const alpha) {
184  *vectorPtr *= alpha;
185  return *this;
186  }
187 
189  operator=(int const k) {
190  *vectorPtr = k;
191  return *this;
192  }
193 
197 
198 
199  inline Tvector const & getVector() const {return *vectorPtr;}
200 
201  std::string obj_type_id() const {return "VectorGeneral";}
202  protected:
204 
205  inline void writeToFileProt(std::ofstream & file) const {
206  if (is_empty())
207  // This means that the object's data structure has not been set
208  return;
209  vectorPtr->writeToFile(file);
210  }
211  inline void readFromFileProt(std::ifstream & file) {
212  if (is_empty())
213  // This means that the object's data structure has not been set
214  return;
215  vectorPtr->readFromFile(file);
216  }
217 
218  inline void inMemorySet(bool inMem) {
219  vectorPtr.inMemorySet(inMem);
220  }
221 
222  private:
223 
224  }; /* end class VectorGeneral */
225 
226 
227  /* LEVEL 2 operations */
228  /* OPERATIONS INVOLVING ORDINARY MATRICES */
230  template<typename Treal, typename Tvector>
231  template<typename Tmatrix>
232  VectorGeneral<Treal, Tvector>&
234  (const XYZ<Treal,
237  assert(!smv.tC);
238  vectorPtr.haveDataStructureSet(true);
239  if ( this == &smv.C ) {
240  // We need a copy of the smv.C vector since it is the same as *this
242  Tvector::gemv(smv.tB, smv.A, smv.B.getMatrix(),
243  *tmp.vectorPtr, 0, *this->vectorPtr);
244  }
245  else
246  Tvector::gemv(smv.tB, smv.A, smv.B.getMatrix(),
247  *smv.C.vectorPtr, 0, *this->vectorPtr);
248  return *this;
249  }
250 
252  template<typename Treal, typename Tvector>
253  template<typename Tmatrix>
256  (const XYZ<Treal,
259  assert(!smv.tC);
260  assert(this != &smv.C);
261  Tvector::gemv(smv.tB, smv.A, smv.B.getMatrix(),
262  *smv.C.vectorPtr, 1, *this->vectorPtr);
263  return *this;
264  }
265 
266 
268  template<typename Treal, typename Tvector>
269  template<typename Tmatrix>
272  (const XYZpUV<Treal,
275  Treal,
276  VectorGeneral<Treal, Tvector> >& smvpsv) {
277  assert(!smvpsv.tC && !smvpsv.tE);
278  assert(this != &smvpsv.C);
279  if (this == &smvpsv.E)
280  Tvector::gemv(smvpsv.tB, smvpsv.A, smvpsv.B.getMatrix(),
281  *smvpsv.C.vectorPtr, smvpsv.D, *this->vectorPtr);
282  else
283  throw Failure("VectorGeneral<Treal, Tvector>::operator="
284  "(const XYZpUV<Treal, "
285  "MatrixGeneral<Treal, Tmatrix>, "
286  "VectorGeneral<Treal, Tvector>, "
287  "Treal, "
288  "VectorGeneral<Treal, Tvector> >&) : "
289  "y = alpha * op(A) * x + beta * z "
290  "not supported for z != y");
291  return *this;
292  }
293 
294 
295  /* OPERATIONS INVOLVING SYMMETRIC MATRICES */
297  template<typename Treal, typename Tvector>
298  template<typename Tmatrix>
301  (const XYZ<Treal,
304  assert(!smv.tC);
305  assert(this != &smv.C);
306  vectorPtr.haveDataStructureSet(true);
307  Tvector::symv('U', smv.A, smv.B.getMatrix(),
308  *smv.C.vectorPtr, 0, *this->vectorPtr);
309  return *this;
310  }
311 
312 
314  template<typename Treal, typename Tvector>
315  template<typename Tmatrix>
318  (const XYZ<Treal,
321  assert(!smv.tC);
322  assert(this != &smv.C);
323  Tvector::symv('U', smv.A, smv.B.getMatrix(),
324  *smv.C.vectorPtr, 1, *this->vectorPtr);
325  return *this;
326  }
327 
329  template<typename Treal, typename Tvector>
330  template<typename Tmatrix>
333  (const XYZpUV<Treal,
336  Treal,
337  VectorGeneral<Treal, Tvector> >& smvpsv) {
338  assert(!smvpsv.tC && !smvpsv.tE);
339  assert(this != &smvpsv.C);
340  if (this == &smvpsv.E)
341  Tvector::symv('U', smvpsv.A, smvpsv.B.getMatrix(),
342  *smvpsv.C.vectorPtr, smvpsv.D, *this->vectorPtr);
343  else
344  throw Failure("VectorGeneral<Treal, Tvector>::operator="
345  "(const XYZpUV<Treal, "
346  "MatrixSymmetric<Treal, Tmatrix>, "
347  "VectorGeneral<Treal, Tvector>, "
348  "Treal, "
349  "VectorGeneral<Treal, Tvector> >&) : "
350  "y = alpha * A * x + beta * z "
351  "not supported for z != y");
352  return *this;
353  }
354 
355  /* OPERATIONS INVOLVING TRIANGULAR MATRICES */
358  template<typename Treal, typename Tvector>
359  template<typename Tmatrix>
364  assert(!mv.tB);
365  if (this != &mv.B)
366  throw Failure("y = A * x not supported for y != x ");
367  Tvector::trmv('U', mv.tA,
368  mv.A.getMatrix(),
369  *this->vectorPtr);
370  return *this;
371  }
372 
373  /* LEVEL 1 operations */
374 
376  template<typename Treal, typename Tvector>
380  assert(!sv.tB);
381  assert(this != &sv.B);
382  Tvector::axpy(sv.A, *sv.B.vectorPtr, *this->vectorPtr);
383  return *this;
384  }
385 
386 
387 
388  /* Defined outside class */
392  template<typename Treal, typename Tvector>
394  VectorGeneral<Treal, Tvector> const & y) {
395  if (xT.tA == false)
396  throw Failure("operator*("
397  "Xtrans<VectorGeneral<Treal, Tvector> > const &,"
398  " VectorGeneral<Treal, Tvector> const &): "
399  "Dimension mismatch in vector operation");
400  return Tvector::dot(xT.A.getVector(), y.getVector());
401  }
402 
403 
404 
405 } /* end namespace mat */
406 #endif
407 
FileWritable.h
mat::VectorGeneral::VectorGeneral
VectorGeneral(const VectorGeneral< Treal, Tvector > &other)
Definition: VectorGeneral.h:73
mat::axpy
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
mat::ValidPtr< Tvector >
mat::FileWritable
Write and read objects to/from file.
Definition: FileWritable.h:56
mat::symv
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
mat::VectorGeneral::operator=
VectorGeneral< Treal, Tvector > & operator=(const VectorGeneral< Treal, Tvector > &other)
Definition: VectorGeneral.h:92
mat::VectorGeneral::clear
void clear()
Definition: VectorGeneral.h:100
matrix_proxy.h
mat::VectorGeneral::eucl
Treal eucl() const
Definition: VectorGeneral.h:178
ValidPtr.h
mat::gemv
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
mat::XYZ
This proxy expresses the result of multiplication of three objects, of possibly different types,...
Definition: matrix_proxy.h:67
mat::VectorGeneral::operator*=
VectorGeneral< Treal, Tvector > & operator*=(Treal const alpha)
Definition: VectorGeneral.h:183
mat::MatrixSymmetric
Symmetric matrix.
Definition: MatrixBase.h:51
mat::VectorGeneral::resetSizesAndBlocks
void resetSizesAndBlocks(SizesAndBlocks const &newRows)
Definition: VectorGeneral.h:51
mat::VectorGeneral::inMemorySet
void inMemorySet(bool inMem)
Definition: VectorGeneral.h:218
mat::VectorGeneral::VectorGeneral
VectorGeneral(SizesAndBlocks const &newRows)
Definition: VectorGeneral.h:64
mat::Failure
Definition: Failure.h:57
mat::ValidPtr::haveDataStructureSet
void haveDataStructureSet(bool val)
Definition: ValidPtr.h:99
mat::VectorGeneral::is_empty
bool is_empty() const
Definition: VectorGeneral.h:56
mat::VectorGeneral
Definition: MatrixBase.h:55
mat::MatrixGeneral
Normal matrix.
Definition: MatrixBase.h:49
mat::VectorGeneral::getVector
Tvector const & getVector() const
Definition: VectorGeneral.h:199
mat
Definition: allocate.cc:39
mat::VectorGeneral::operator=
VectorGeneral< Treal, Tvector > & operator=(int const k)
Definition: VectorGeneral.h:189
mat::VectorGeneral::obj_type_id
std::string obj_type_id() const
Definition: VectorGeneral.h:201
mat::ValidPtr::inMemorySet
void inMemorySet(bool val)
Definition: ValidPtr.h:93
mat::XYZpUV
This proxy expresses the result of multiplication of three objects added to two other multiplied obje...
Definition: matrix_proxy.h:88
mat::VectorGeneral::readFromFileProt
void readFromFileProt(std::ifstream &file)
Definition: VectorGeneral.h:211
mat::VectorGeneral::rand
void rand()
Definition: VectorGeneral.h:108
mat::VectorGeneral::fullvector
void fullvector(std::vector< Treal > &fullVector) const
Definition: VectorGeneral.h:88
mat::trmv
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
mat::VectorGeneral::writeToFileProt
void writeToFileProt(std::ofstream &file) const
Definition: VectorGeneral.h:205
mat::XY
This proxy expresses the result of multiplication of two objects, of possibly different types,...
Definition: matrix_proxy.h:51
mat::Xtrans
This proxy expresses the result of transposition of an object of type TX.
Definition: matrix_proxy.h:118
mat::SizesAndBlocks
Describes dimensions of matrix and its blocks on all levels.
Definition: SizesAndBlocks.h:45
mat::VectorGeneral::assign_from_full
void assign_from_full(std::vector< Treal > const &fullVector, SizesAndBlocks const &newRows)
Definition: VectorGeneral.h:83
mat::VectorGeneral::vectorPtr
ValidPtr< Tvector > vectorPtr
Definition: VectorGeneral.h:203
mat::VectorGeneral::clear_structure
void clear_structure()
Definition: VectorGeneral.h:60
mat::VectorGeneral::VectorGeneral
VectorGeneral()
Definition: VectorGeneral.h:67
mat::dot
static T dot(const int *n, const T *dx, const int *incx, const T *dy, const int *incy)
Definition: mat_gblas.h:425
mat::ValidPtr::haveDataStructureGet
bool haveDataStructureGet() const
Definition: ValidPtr.h:102
mat::operator*
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