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 
30 /** @file VectorGeneral.h General vector class
31  *
32  * Copyright(c) Emanuel Rubensson 2006
33  *
34  * @author Emanuel Rubensson  @a responsible @a author
35  * @date January 2006
36  *
37  */
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 
resetSizesAndBlocks(SizesAndBlocks const & newRows)51     inline void resetSizesAndBlocks(SizesAndBlocks const & newRows) {
52       vectorPtr.haveDataStructureSet(true);
53       vectorPtr->resetRows(newRows);
54     }
55 
is_empty()56     inline bool is_empty() const {
57       return !vectorPtr.haveDataStructureGet();
58     }
59 
clear_structure()60     inline void clear_structure(){
61       vectorPtr.haveDataStructureSet(false);
62     }
63 
VectorGeneral(SizesAndBlocks const & newRows)64     VectorGeneral(SizesAndBlocks const & newRows):vectorPtr(new Tvector) {
65       resetSizesAndBlocks(newRows);
66     }
VectorGeneral()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
VectorGeneral(const VectorGeneral<Treal,Tvector> & other)71     explicit VectorGeneral(const VectorGeneral<Treal, Tvector>& other)
72 #else
73     VectorGeneral(const VectorGeneral<Treal, Tvector>& other)
74 #endif
75       :FileWritable(other), vectorPtr(new Tvector) {
76       if (other.vectorPtr.haveDataStructureGet()) {
77 	vectorPtr.haveDataStructureSet(true);
78       }
79       *vectorPtr = *other.vectorPtr;
80     }
81 
assign_from_full(std::vector<Treal> const & fullVector,SizesAndBlocks const & newRows)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     }
fullvector(std::vector<Treal> & fullVector)88     inline void fullvector(std::vector<Treal> & fullVector) const {
89       this->vectorPtr->fullVector(fullVector);
90     }
91     VectorGeneral<Treal, Tvector>&
92       operator=(const VectorGeneral<Treal, Tvector>& other) {
93       if (other.vectorPtr.haveDataStructureGet()) {
94 	vectorPtr.haveDataStructureSet(true);
95       }
96       *this->vectorPtr = *other.vectorPtr;
97       return *this;
98     }
99 
clear()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 
rand()108     inline void rand() {
109       vectorPtr->randomNormalized();
110     }
111 
112     /* LEVEL 2 operations */
113     /* OPERATIONS INVOLVING ORDINARY MATRICES */
114     /** y = alpha * op(A) * x */
115     template<typename Tmatrix>
116       VectorGeneral<Treal, Tvector>& operator=
117       (const XYZ<Treal,
118        MatrixGeneral<Treal, Tmatrix>,
119        VectorGeneral<Treal, Tvector> >& smv);
120 
121     /** y += alpha * op(A) * x */
122     template<typename Tmatrix>
123       VectorGeneral<Treal, Tvector>& operator+=
124       (const XYZ<Treal,
125        MatrixGeneral<Treal, Tmatrix>,
126        VectorGeneral<Treal, Tvector> >& smv);
127     /** y = alpha * op(A) * x + beta * y */
128     template<typename Tmatrix>
129       VectorGeneral<Treal, Tvector>& operator=
130       (const XYZpUV<Treal,
131        MatrixGeneral<Treal, Tmatrix>,
132        VectorGeneral<Treal, Tvector>,
133        Treal,
134        VectorGeneral<Treal, Tvector> >& smvpsv);
135 
136     /** y = op(A) * x                      : A is general */
137     template<typename Tmatrix>
138       VectorGeneral<Treal, Tvector>& operator=
139       (const XY<MatrixGeneral<Treal, Tmatrix>,
140        VectorGeneral<Treal, Tvector> >& mv) {
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 */
148     /** y = alpha * A * x                      : A is symmetric */
149     template<typename Tmatrix>
150       VectorGeneral<Treal, Tvector>& operator=
151       (const XYZ<Treal,
152        MatrixSymmetric<Treal, Tmatrix>,
153        VectorGeneral<Treal, Tvector> >& smv);
154     /** y += alpha * A * x                     : A is symmetric */
155     template<typename Tmatrix>
156       VectorGeneral<Treal, Tvector>& operator+=
157       (const XYZ<Treal,
158        MatrixSymmetric<Treal, Tmatrix>,
159        VectorGeneral<Treal, Tvector> >& smv);
160     /** y = alpha * A * x + beta * y           : A is symmetric */
161     template<typename Tmatrix>
162       VectorGeneral<Treal, Tvector>& operator=
163       (const XYZpUV<Treal,
164        MatrixSymmetric<Treal, Tmatrix>,
165        VectorGeneral<Treal, Tvector>,
166        Treal,
167        VectorGeneral<Treal, Tvector> >& smvpsv);
168 
169     /* OPERATIONS INVOLVING TRIANGULAR MATRICES */
170     /** y = op(A) * x                      : A is triangular */
171     template<typename Tmatrix>
172       VectorGeneral<Treal, Tvector>& operator=
173       (const XY<MatrixTriangular<Treal, Tmatrix>,
174        VectorGeneral<Treal, Tvector> >& mv);
175 
176 
177     /* LEVEL 1 operations */
eucl()178     inline Treal eucl() const {
179       return vectorPtr->eucl();
180     }
181 
182     inline VectorGeneral<Treal, Tvector>&
183       operator*=(Treal const alpha) {
184       *vectorPtr *= alpha;
185       return *this;
186     }
187 
188     inline VectorGeneral<Treal, Tvector>&
189       operator=(int const k) {
190       *vectorPtr = k;
191       return *this;
192     }
193 
194     /** y += alpha * x */
195     VectorGeneral<Treal, Tvector>& operator+=
196       (const XY<Treal, VectorGeneral<Treal, Tvector> >& sv);
197 
198 
getVector()199     inline Tvector const & getVector() const {return *vectorPtr;}
200 
obj_type_id()201     std::string obj_type_id() const {return "VectorGeneral";}
202   protected:
203     ValidPtr<Tvector> vectorPtr;
204 
writeToFileProt(std::ofstream & file)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     }
readFromFileProt(std::ifstream & file)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 
inMemorySet(bool inMem)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 */
229     /** y = alpha * op(A) * x */
230   template<typename Treal, typename Tvector>
231     template<typename Tmatrix>
232     VectorGeneral<Treal, Tvector>&
233     VectorGeneral<Treal, Tvector>::operator=
234     (const XYZ<Treal,
235      MatrixGeneral<Treal, Tmatrix>,
236      VectorGeneral<Treal, Tvector> >& smv) {
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
241       VectorGeneral<Treal, Tvector> tmp(smv.C);
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 
251   /** y += alpha * op(A) * x */
252   template<typename Treal, typename Tvector>
253     template<typename Tmatrix>
254     VectorGeneral<Treal, Tvector>&
255     VectorGeneral<Treal, Tvector>::operator+=
256     (const XYZ<Treal,
257      MatrixGeneral<Treal, Tmatrix>,
258      VectorGeneral<Treal, Tvector> >& smv) {
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 
267   /** y = alpha * op(A) * x + beta * y */
268   template<typename Treal, typename Tvector>
269     template<typename Tmatrix>
270     VectorGeneral<Treal, Tvector>&
271     VectorGeneral<Treal, Tvector>::operator=
272     (const XYZpUV<Treal,
273      MatrixGeneral<Treal, Tmatrix>,
274      VectorGeneral<Treal, Tvector>,
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 */
296   /** y = alpha * A * x                      : A is symmetric */
297   template<typename Treal, typename Tvector>
298     template<typename Tmatrix>
299     VectorGeneral<Treal, Tvector>&
300     VectorGeneral<Treal, Tvector>::operator=
301     (const XYZ<Treal,
302      MatrixSymmetric<Treal, Tmatrix>,
303      VectorGeneral<Treal, Tvector> >& smv) {
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 
313   /** y += alpha * A * x                     : A is symmetric */
314   template<typename Treal, typename Tvector>
315     template<typename Tmatrix>
316     VectorGeneral<Treal, Tvector>&
317     VectorGeneral<Treal, Tvector>::operator+=
318     (const XYZ<Treal,
319      MatrixSymmetric<Treal, Tmatrix>,
320      VectorGeneral<Treal, Tvector> >& smv) {
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 
328   /** y = alpha * A * x + beta * y           : A is symmetric */
329   template<typename Treal, typename Tvector>
330     template<typename Tmatrix>
331     VectorGeneral<Treal, Tvector>&
332     VectorGeneral<Treal, Tvector>::operator=
333     (const XYZpUV<Treal,
334      MatrixSymmetric<Treal, Tmatrix>,
335      VectorGeneral<Treal, Tvector>,
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 */
356   /** x = op(A) * x                      : A is triangular */
357 
358   template<typename Treal, typename Tvector>
359     template<typename Tmatrix>
360     VectorGeneral<Treal, Tvector>&
361     VectorGeneral<Treal, Tvector>::operator=
362     (const XY<MatrixTriangular<Treal, Tmatrix>,
363      VectorGeneral<Treal, Tvector> >& mv) {
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 
375   /** y += alpha * x */
376   template<typename Treal, typename Tvector>
377     VectorGeneral<Treal, Tvector>&
378     VectorGeneral<Treal, Tvector>::operator+=
379     (const XY<Treal, VectorGeneral<Treal, Tvector> >& sv) {
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 */
389   /** transpose(x) * y
390    * Scalar (dot) product of two vectors
391    */
392   template<typename Treal, typename Tvector>
393     Treal operator*(Xtrans<VectorGeneral<Treal, Tvector> > const & xT,
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 
408