GCC Code Coverage Report

Directory: ./
File: PETScInterface/petscVector.hpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 4 4 100.0%
Branches: 0 0 -%
Line Branch Exec Source
1 // ______ ______ _ _ _____ ______
2 // | ____| ____| | (_)/ ____| | ____|
3 // | |__ | |__ | | _| (___ ___| |__
4 // | __| | __| | | | |\___ \ / __| __|
5 // | | | |____| |____| |____) | (__| |____
6 // |_| |______|______|_|_____/ \___|______|
7 // Finite Elements for Life Sciences and Engineering
8 //
9 // License: LGL2.1 License
10 // FELiScE default license: LICENSE in root folder
11 //
12 // Main authors: Julien Castelneau
13 // Dominique Chapelle
14 // Miguel Fernandez
15 // Vicente Mataix Ferrandiz
16 //
17
18 #ifndef PETSCVECTOR_HPP
19 #define PETSCVECTOR_HPP
20
21 // System includes
22 #include <string>
23 #include <vector>
24
25 // External includes
26 #include "Core/NoThirdPartyWarning/Petsc/sys.hpp" // Jacques 13/07/22: added due to definition of FELISCE_PETSC_NULLPTR
27 #include "Core/NoThirdPartyWarning/Petsc/vec.hpp"
28 #include "Core/NoThirdPartyWarning/Petsc/ao.hpp"
29
30 // Project includes
31 #include "PETScInterface/petscMatrix_fwd.hpp"
32 #include "Core/shared_pointers.hpp"
33
34 namespace felisce
35 {
36
37 //! \brief Wrapper to Petsc Vector type (Vec)
38 class PetscVector
39 {
40 //===============
41 // FRIEND CLASSES
42 //===============
43
44 friend class PetscMatrix;
45 friend class PetscVectorScatter;
46
47 //=================
48 // FRIEND FUNCTIONS
49 //=================
50
51 // TODO: Define more operators using this friend functions
52 friend PetscErrorCode dot(const PetscVector& x,const PetscVector& y,PetscScalar* val);
53 friend PetscErrorCode mult(const PetscMatrix& mat,const PetscVector& x,PetscVector& y);
54 friend PetscErrorCode multAdd(const PetscMatrix& mat,const PetscVector& v1,const PetscVector& v2,PetscVector& v3);
55 friend PetscErrorCode multTranspose(const PetscMatrix& mat,const PetscVector& x,PetscVector& y);
56 friend PetscErrorCode pointwiseDivide(const PetscVector& w, const PetscVector& x,const PetscVector& y);
57 friend PetscErrorCode pointwiseMult(const PetscVector& w, const PetscVector& x,const PetscVector& y);
58 friend PetscErrorCode swap(const PetscVector& x,const PetscVector& y);
59 friend PetscErrorCode waxpy(const PetscVector& w,PetscScalar const & alpha,const PetscVector& x,const PetscVector& y);
60 friend PetscErrorCode maxPointwiseDivide(const PetscVector& x, const PetscVector& y, PetscScalar* max);
61 friend PetscErrorCode gatherVec(PetscVector& v, PetscVector& seqV);
62 public:
63
64 /// Pointer definition of PetscVector
65 FELISCE_CLASS_POINTER_DEFINITION(PetscVector);
66
67 //=======================
68 // CONSTRUCTOR/DESTRUCTOR
69 //=======================
70
71 /**
72 * @brief Default constructor
73 */
74 PetscVector();
75
76 /**
77 * @brief Default destructor
78 */
79 ~PetscVector();
80
81 /**
82 * @brief Copy constructor
83 * @param b The vector to be copied
84 */
85 PetscVector(const PetscVector& b);
86
87 /**
88 * @brief Assign operator
89 * @param b The vector to be assigned
90 * @return The current instance
91 */
92 PetscVector& operator=(const PetscVector& b);
93
94 /**
95 * @brief Multiply (*) operator
96 * @param scalar The scalar to be multiplied
97 * @return The multiplication of the current vector and scalar
98 */
99 PetscVector operator*(const PetscInt scalar);
100
101 /**
102 * @brief Multiply (*) operator (constant version)
103 * @param scalar The scalar to be multiplied
104 * @return The multiplication of the current vector and scalar
105 */
106 PetscVector operator*(const PetscInt scalar) const;
107
108 /**
109 * @brief Multiply (*) operator
110 * @param scalar The scalar to be multiplied
111 * @return The multiplication of the current vector and scalar
112 */
113 PetscVector operator*(const PetscScalar scalar);
114
115 /**
116 * @brief Multiply (*) operator (constant version)
117 * @param scalar The scalar to be multiplied
118 * @return The multiplication of the current vector and scalar
119 */
120 PetscVector operator*(const PetscScalar scalar) const;
121
122 /**
123 * @brief Multiply (*) operator (dot product)
124 * @param b The vector to be multiplied
125 * @return The dot product of the current vector and b
126 */
127 PetscScalar operator*(const PetscVector& b);
128
129 /**
130 * @brief Multiply (*) operator (dot product) (constant version)
131 * @param b The vector to be multiplied
132 * @return The dot product of the current vector and b
133 */
134 PetscScalar operator*(const PetscVector& b) const;
135
136 /**
137 * @brief Multiply (*=) operator
138 * @param scalar The scalar to be multiplied
139 * @return The current vector multiplied by the scalar
140 */
141 PetscVector& operator*=(const PetscInt scalar);
142
143 /**
144 * @brief Multiply (*=) operator
145 * @param scalar The scalar to be multiplied
146 * @return The current vector multiplied by the scalar
147 */
148 PetscVector& operator*=(const PetscScalar scalar);
149
150 /**
151 * @brief Division (/) operator
152 * @param scalar The scalar to be divided
153 * @return The division of the current vector and scalar
154 */
155 PetscVector operator/(const PetscInt scalar);
156
157 /**
158 * @brief Division (/) operator (constant version)
159 * @param scalar The scalar to be divided
160 * @return The division of the current vector and scalar
161 */
162 PetscVector operator/(const PetscInt scalar) const;
163
164 /**
165 * @brief Division (/) operator
166 * @param scalar The scalar to be divided
167 * @return The division of the current vector and scalar
168 */
169 PetscVector operator/(const PetscScalar scalar);
170
171 /**
172 * @brief Division (/) operator (constant version)
173 * @param scalar The scalar to be divided
174 * @return The division of the current vector and scalar
175 */
176 PetscVector operator/(const PetscScalar scalar) const;
177
178 /**
179 * @brief Division (/=) operator
180 * @param scalar The scalar to be divided
181 * @return The current vector divided by the scalar
182 */
183 PetscVector& operator/=(const PetscInt scalar);
184
185 /**
186 * @brief Division (/=) operator
187 * @param scalar The scalar to be divided
188 * @return The current vector divided by the scalar
189 */
190 PetscVector& operator/=(const PetscScalar scalar);
191
192 /**
193 * @brief Sum (+) operator
194 * @param b The vector to be summed
195 * @return The sum product of the current vector and b
196 */
197 PetscVector operator+(const PetscVector& b);
198
199 /**
200 * @brief Sum (+) operator (constant version)
201 * @param b The vector to be summed
202 * @return The sum product of the current vector and b
203 */
204 PetscVector operator+(const PetscVector& b) const;
205
206 /**
207 * @brief Sum (+=) operator
208 * @param b The vector to be summed
209 * @return The current vector plus b
210 */
211 PetscVector& operator+=(const PetscVector& b);
212
213 /**
214 * @brief Difference (-) operator
215 * @param b The vector to be differenced
216 * @return The difference product of the current vector and b
217 */
218 PetscVector operator-(const PetscVector& b);
219
220 /**
221 * @brief Difference (-) operator (constant version)
222 * @param b The vector to be differenced
223 * @return The difference product of the current vector and b
224 */
225 PetscVector operator-(const PetscVector& b) const;
226
227 /**
228 * @brief Difference (-=) operator
229 * @param b The vector to be differenced
230 * @return The current vector minus b
231 */
232 PetscVector& operator-=(const PetscVector& b);
233
234 /**
235 * @brief [] operator
236 * @param Index The index of the value
237 * @note Cannot be reference as uses getValue method
238 */
239 PetscScalar operator[] (const PetscInt Index);
240
241 /**
242 * @brief [] operator (const version)
243 * @param Index The index of the value
244 * @note Cannot be reference as uses getValue method
245 */
246 PetscScalar operator[] (const PetscInt Index) const;
247
248 //=================
249 // STATIC FUNCTIONS
250 //=================
251
252 static PetscVector null();
253
254 //===================
255 // CREATION FUNCTIONS
256 //===================
257
258 /**
259 * @brief This is the method that must be called to create the PETSc vector
260 * @return The corresponding error code
261 */
262 PetscErrorCode create(MPI_Comm comm = PETSC_COMM_WORLD);
263
264 /**
265 * @brief This is the method that must be called to create the PETSc vector (in sequantial)
266 * @param comm The communicator
267 * @param n The size of the vector
268 * @return The corresponding error code
269 */
270 PetscErrorCode createSeq(MPI_Comm comm,PetscInt n);
271
272 /**
273 * @brief This is the method that must be called to create the PETSc vector (in sequantial)
274 * @note Arguments reordered to consider PETSC_COMM_WORLD as default
275 * @param n The size of the vector
276 * @param comm The communicator
277 * @return The corresponding error code
278 */
279 PetscErrorCode createSeqInstance(PetscInt n, MPI_Comm comm = PETSC_COMM_WORLD) {return createSeq(comm, n);}
280
281 /**
282 * @brief This is the method that must be called to create the PETSc vector (in MPI)
283 * @param comm The communicator
284 * @param n The size per partition
285 * @param N The size of the vector
286 * @return The corresponding error code
287 */
288 PetscErrorCode createMPI(MPI_Comm comm,PetscInt n,PetscInt N);
289
290 /**
291 * @brief This is the method that must be called to create the PETSc vector (in MPI)
292 * @note Arguments reordered to consider PETSC_COMM_WORLD as default
293 * @param N The size of the vector
294 * @param n The size per partition
295 * @param comm The communicator
296 * @return The corresponding error code
297 */
298 120 PetscErrorCode createMPIInstance(PetscInt N, PetscInt n = PETSC_DECIDE, MPI_Comm comm = PETSC_COMM_WORLD) {return createMPI(comm, n, N);}
299
300 //======================
301 // DESTRUCTION FUNCTIONS
302 //======================
303
304 /**
305 * @brief Destroys the current vector (to be called as "destructor")
306 * @return The corresponding error code
307 */
308 PetscErrorCode destroy();
309
310 //===============
311 // COPY FUNCTIONS
312 //===============
313
314 /**
315 * @brief Copies a vector.
316 * @param vector A vector to copy
317 * @return The corresponding error code
318 */
319 PetscErrorCode copyFrom(const PetscVector& vector) const;
320
321 /**
322 * @brief Creates a new vector of the same type as an existing vector.
323 * @param vector A vector to mimic
324 * @return The corresponding error code
325 */
326 PetscErrorCode duplicateFrom(const PetscVector& vector);
327
328 //=======================
329 // MODIFICATION FUNCTIONS
330 //=======================
331
332 /**
333 * @brief This routine should be called after completing all calls to VecSetValues()
334 * @return The corresponding error code
335 */
336 PetscErrorCode assembly() const;
337
338 PetscErrorCode restoreArray(PetscScalar **a) const;
339
340 PetscErrorCode set(const PetscScalar& value) const;
341
342 PetscErrorCode setFromOptions() const;
343
344 /**
345 * @brief Sets an option for controling a vector's behavior
346 * @param op The option
347 * @param flg Flag turn the option on/off
348 * @return The corresponding error code
349 */
350 PetscErrorCode setOption(VecOption op, PetscBool flg) const;
351
352 /**
353 * @brief Sets the local and global sizes, and checks to determine compatibility
354 * @param n The local size (or PETSC_DECIDE to have it set)
355 * @param N The global size (or PETSC_DECIDE)
356 * @return The corresponding error code
357 */
358 PetscErrorCode setSizes(PetscInt n, PetscInt N) const;
359
360 /**
361 * @brief Sets the local and global sizes, and checks to determine compatibility
362 * @note Version with size across partitions automatic
363 * @param N The global size
364 * @return The corresponding error code
365 */
366 4 PetscErrorCode setSizes(PetscInt N) const {return setSizes(PETSC_DECIDE, N);}
367
368 /**
369 * @brief Builds a vector, for a particular vector implementation.
370 * @details See "petsc/include/petscvec.h" for available vector types (for instance, VECSEQ, VECMPI, or VECSHARED).
371 * @param method The name of the vector type
372 * @return The corresponding error code
373 */
374 PetscErrorCode setType(VecType method = VECSEQ) const;
375
376 /**
377 * @brief Set a single entry into a vector
378 * @param row The row location of the entry
379 * @param value The value to insert
380 * @param mode Either INSERT_VALUES or ADD_VALUES
381 * @return The corresponding error code
382 */
383 PetscErrorCode setValue(PetscInt row, PetscScalar value, InsertMode mode = INSERT_VALUES) const;
384
385 /**
386 * @brief Set a single entry into a vector. Only in the rank 0
387 * @param row The row location of the entry
388 * @param value The value to insert
389 * @param mode Either INSERT_VALUES or ADD_VALUES
390 * @return The corresponding error code
391 */
392 PetscErrorCode setValueOnce(PetscInt row, PetscScalar value, InsertMode mode = INSERT_VALUES) const;
393
394 /**
395 * @brief Inserts or adds values into certain locations of a vector.
396 * @param ni Number of elements to add
397 * @param ix Indices where to add
398 * @param y Array of values
399 * @param mode Either INSERT_VALUES or ADD_VALUES, where ADD_VALUES adds values to any existing entries, and INSERT_VALUES replaces existing entries with new values
400 * @return The corresponding error code
401 */
402 PetscErrorCode setValues(PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora = INSERT_VALUES) const;
403
404 /**
405 * @brief Inserts or adds values into certain locations of a vector. Only in the rank 0
406 * @param ni Number of elements to add
407 * @param ix Indices where to add
408 * @param y Array of values
409 * @param mode Either INSERT_VALUES or ADD_VALUES, where ADD_VALUES adds values to any existing entries, and INSERT_VALUES replaces existing entries with new values
410 * @return The corresponding error code
411 */
412 PetscErrorCode setValuesOnce(PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora = INSERT_VALUES) const;
413
414 /**
415 * @brief Set to zero all the values of the vector
416 * @return The corresponding error code
417 */
418 PetscErrorCode zeroEntries() const;
419
420 //=================
421 // ACCESS FUNCTIONS
422 //=================
423
424 /**
425 * @brief Retrieves the C array contained in the vector
426 * @param array The array containing the values
427 * @return The corresponding error code
428 */
429 PetscErrorCode getArray(PetscScalar** array) const;
430
431 /**
432 * @brief Retrieves the C array contained in the vector
433 * @note Direct version
434 * @return The array containing the values
435 */
436 PetscScalar* getArray() const;
437
438 /**
439 * @brief Retrieves the local size of the vector
440 * @param localSize The local size
441 * @return The corresponding error code
442 */
443 PetscErrorCode getLocalSize(PetscInt* localSize) const;
444
445 /**
446 * @brief Retrieves the local size of the vector
447 * @note Direct version
448 * @return The local size
449 */
450 PetscInt getLocalSize() const;
451
452 /**
453 * @brief Retrieves the total size of the vector
454 * @param size The total size
455 * @return The corresponding error code
456 */
457 PetscErrorCode getSize(PetscInt* size) const;
458
459 /**
460 * @brief Retrieves the total size of the vector
461 * @note Direct version
462 * @return The total size
463 */
464 PetscInt getSize() const;
465
466 /**
467 * @brief Gets values from certain locations of a vector. Currently can only get values on the same processor
468 * @param ni Number of elements to get
469 * @param ix Indices where to get them from (in global 1d numbering)
470 * @param y Array of values
471 * @return The corresponding error code
472 */
473 PetscErrorCode getValues(PetscInt ni,const PetscInt ix[],PetscScalar y[]) const;
474
475 /**
476 * @brief Gets values from certain locations of a vector. Currently can only get values on the same processor
477 * @param ni Number of elements to get
478 * @param ix Index where to get them from (in global 1d numbering)
479 * @param y Value
480 * @return The corresponding error code
481 */
482 PetscErrorCode getValue(const PetscInt ix,PetscScalar& y) const;
483
484 PetscErrorCode getAllValuesInAppOrdering(const AO& ao, std::vector<double>& array) const;
485
486 //==================
487 // BOOLEAN FUNCTIONS
488 //==================
489
490 /**
491 * @brief This returns if the vector is null
492 * @return If the vector is null
493 */
494 bool isNull();
495
496 /**
497 * @brief This returns if the vector is null (constant version)
498 * @return If the vector is null
499 */
500 bool isNull() const;
501
502 /**
503 * @brief This returns if the vector is not null
504 * @return If the vector is not null
505 */
506 bool isNotNull();
507
508 /**
509 * @brief This returns if the vector is not null (constant version)
510 * @return If the vector is not null
511 */
512 bool isNotNull() const;
513
514 //=======================
515 // MATHEMATICAL FUNCTIONS
516 //=======================
517
518 PetscErrorCode axpbypcz(PetscScalar alpha, PetscScalar beta, PetscScalar gamma, const PetscVector& x, const PetscVector& y) const;
519
520 /**
521 * @brief Compute linear combination of two vectors (y = alpha * x + beta * y)
522 * @details x and y MUST be different vectors The implementation is optimized for alpha and/or beta values of 0.0 and 1.0
523 * @param alpha The first scalar coefficient
524 * @param beta The second scalar coefficient
525 * @param x The first vector
526 * @return The corresponding error code
527 */
528 PetscErrorCode axpby(PetscScalar alpha, PetscScalar beta, const PetscVector& x) const;
529
530 /**
531 * @brief Compute linear combination of two vectors (y = alpha * x + beta * y)
532 * @details x and y MUST be different vectors The implementation is optimized for alpha and/or beta values of 0.0 and 1.0
533 * @note The order is different so defaults arguments of 1.0 can be considered
534 * @param alpha The first scalar coefficient
535 * @param beta The second scalar coefficient
536 * @param x The first vector
537 * @return The corresponding error code
538 */
539 8 PetscErrorCode axpby(const PetscVector& x, const PetscScalar alpha = 1.0, const PetscScalar beta = 1.0) const {return axpby(alpha, beta, x);}
540
541 /**
542 * @brief Compute linear combination of two vectors (y = alpha x + y)
543 * @details x and y MUST be different vectors This routine is optimized for alpha of 0.0
544 * @param alpha The scalar coefficient
545 * @param x The first vector
546 * @return The corresponding error code
547 */
548 PetscErrorCode axpy(PetscScalar alpha, const PetscVector& x) const;
549
550 /**
551 * @brief Compute linear combination of two vectors (y = alpha x + y)
552 * @details x and y MUST be different vectors This routine is optimized for alpha of 0.0
553 * @note The order is different so defaults arguments of 1.0 can be considered
554 * @param alpha The scalar coefficient
555 * @param x The first vector
556 * @return The corresponding error code
557 */
558 28 PetscErrorCode axpy(const PetscVector& x, const PetscScalar alpha = 1.0) const {return axpy(alpha, x);}
559
560 /**
561 * @brief Computes the maximum value of the current vector
562 * @param val The maximum value
563 * @param p The position in the vector (null by default)
564 * @return The corresponding error code
565 */
566 PetscErrorCode max(PetscReal* val, PetscInt* p = FELISCE_PETSC_NULLPTR) const;
567
568 /**
569 * @brief Computes the maximum value of the current vector
570 * @param p The position in the vector
571 * @param val The maximum value
572 * @return The corresponding error code
573 */
574 PetscErrorCode max(PetscInt* p, PetscReal* val) const;
575
576 /**
577 * @brief Computes the maximum value of the current vector
578 * @details This version has default values for easier computation
579 * @return The maximum value
580 */
581 PetscReal max() const;
582
583 /**
584 * @brief Computes the minimum value of the current vector
585 * @param val The minimum value
586 * @param p The position in the vector (null by default)
587 * @return The corresponding error code
588 */
589 PetscErrorCode min(PetscReal* val, PetscInt* p = FELISCE_PETSC_NULLPTR) const;
590
591 /**
592 * @brief Computes the minimum value of the current vector
593 * @param p The position in the vector
594 * @param val The minimum value
595 * @return The corresponding error code
596 */
597 PetscErrorCode min(PetscInt* p, PetscReal * val) const;
598
599 /**
600 * @brief Computes the minimum value of the current vector
601 * @details This version has default values for easier computation
602 * @return The minimum value
603 */
604 PetscReal min() const;
605
606 /**
607 * @brief Computes the norm of the current vector
608 * @param type The type of norm considered
609 * @param val The normal value
610 * @return The corresponding error code
611 */
612 PetscErrorCode norm(NormType type, PetscReal* val) const;
613
614 /**
615 * @brief Computes the norm of the current vector
616 * @details This version has default values for easier computation
617 * @param type The type of norm considered
618 * @return The normal value
619 */
620 PetscReal norm(const NormType type = NORM_2) const;
621
622 /**
623 * @brief Scale the vector a given value
624 * @param scalar The value that scales the vector
625 * @return The corresponding error code
626 */
627 PetscErrorCode scale(PetscScalar scalar) const;
628
629 /**
630 * @brief Shifts the vector a given value
631 * @param scalar The value that shifts the vector
632 * @return The corresponding error code
633 */
634 PetscErrorCode shift(PetscScalar scalar) const;
635
636 /**
637 * @brief Computes the sum of the current vector
638 * @param sum The sum of the vector
639 * @return The corresponding error code
640 */
641 PetscErrorCode sum(PetscScalar* sum) const;
642
643 /**
644 * @brief Computes the sum of the current vector
645 * @details This version has default values for easier computation
646 * @return The sum of the vector
647 */
648 PetscScalar sum() const;
649
650 /**
651 * @brief Computes absolute value of the vector
652 * @return The corresponding error code
653 */
654 PetscErrorCode abs() const;
655
656 //========================
657 // COMMUNICATION FUNCTIONS
658 //========================
659
660 PetscErrorCode scatterToAll(PetscVector& y, InsertMode addv, ScatterMode mode) const;
661
662 PetscErrorCode scatterToAllNotCreatingVector(PetscVector& y, InsertMode addv, ScatterMode mode) const;
663
664 PetscErrorCode scatterToZero(PetscVector& y, InsertMode addv, ScatterMode mode) const;
665
666 PetscErrorCode scatterToZeroNotCreatingVector(PetscVector& y, InsertMode addv, ScatterMode mode) const;
667
668 PetscErrorCode broadCastSequentialVector(MPI_Comm comm, int master ) const;
669
670 //=================
671 // OUTPUT FUNCTIONS
672 //=================
673
674 PetscErrorCode view(const PetscViewer& viewer = PETSC_VIEWER_STDOUT_WORLD) const;
675
676 PetscErrorCode saveInBinaryFormat(MPI_Comm comm, const std::string filename, const std::string folder) const;
677
678 PetscErrorCode saveInMatlabFormat(MPI_Comm comm, const std::string filename, const std::string folder) const;
679
680 //================
681 // INPUT FUNCTIONS
682 //================
683
684 PetscErrorCode load(PetscViewer viewer);
685
686 PetscErrorCode loadFromBinaryFormat(MPI_Comm comm, const std::string filename, const std::string folder);
687
688 PetscErrorCode loadFromMatlabFormat(MPI_Comm comm, const std::string filename, const std::string folder);
689
690 PetscErrorCode loadFromEnsight(const std::string fileName);
691
692 //=================
693 // TO BE DEPRECATED
694 //=================
695
696 const Vec& toPetsc() const;
697
698 Vec& toPetsc();
699 private:
700
701 //=============
702 // DATA MEMBERS
703 //=============
704
705 Vec m_self = nullptr; /// The PETSc vector
706 };
707
708 }
709
710 #include "petscVector_friend.hpp"
711 #include "petscVector_inline.hpp"
712
713 #endif // PETSCVECTOR_HPP
714