Vector.cpp 28.6 KB
Newer Older
1 2
//! \file 
//
3
//
4
//  Vector.cpp
5
//  HappyHeart
6
//
7
//  Created by Sébastien Gilles on 04/10/13.
8 9 10
//  Copyright (c) 2013 Inria. All rights reserved.
//

11
#include <cassert>
12
#include <cmath>
13

14 15
#include "Utilities/Filesystem/File.hpp"

16
#include "ThirdParty/Wrappers/Petsc/Viewer.hpp"
17
#include "ThirdParty/Wrappers/Petsc/Vector/Vector.hpp"
18
#include "ThirdParty/Wrappers/Petsc/Vector/AccessVectorContent.hpp"
19
#include "ThirdParty/Wrappers/Petsc/Exceptions/Petsc.hpp"
20
#include "ThirdParty/Wrappers/Mpi/Mpi.hpp"
21 22 23 24


namespace HappyHeart
{
25 26
    
    
27
    
28 29
    namespace Wrappers
    {
30 31
        
        
32 33
        namespace Petsc
        {
34 35 36
            
            
            Vector::Vector()
37
            : petsc_vector_(PETSC_NULL),
38
            do_petsc_destroy_(false) // Nothing is assigned yet!
39 40 41
            { }
            
            
42
            Vector::Vector(const Vec& petsc_vector, bool do_destroy_petsc)
43
            : petsc_vector_(petsc_vector),
44
            do_petsc_destroy_(do_destroy_petsc)
45
            {
46
                assert(petsc_vector != PETSC_NULL);
47 48
            }
            
49
            
50
            Vector::Vector(const Vector& rhs)
51 52
            : petsc_vector_(PETSC_NULL),
            do_petsc_destroy_(rhs.do_petsc_destroy_)
53 54 55 56 57
            {
                CompleteCopy(rhs, __FILE__, __LINE__);
            }
            
            
58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81
            Vector::Vector(const Mpi& mpi,
                           const std::string& binary_file,
                           const char* invoking_file, int invoking_line)
            : petsc_vector_(PETSC_NULL),
            do_petsc_destroy_(true)
            {
                const auto communicator = mpi.GetCommunicator();
                
                Viewer viewer(mpi,
                              binary_file,
                              FILE_MODE_READ,
                              invoking_file, invoking_line);
                
                int error_code = VecCreate(communicator, &petsc_vector_);
                if (error_code)
                    throw ExceptionNS::Exception(error_code, "VecCreate", invoking_file, invoking_line);
                
                error_code = VecLoad(petsc_vector_, viewer.GetUnderlyingPetscObject());
                if (error_code)
                    throw ExceptionNS::Exception(error_code, "VecLoad", invoking_file, invoking_line);
            }

            
            
82
            
83 84
            Vector::~Vector()
            {
85
                if (do_petsc_destroy_)
86
                {
87
                    assert(petsc_vector_ != PETSC_NULL);
88
                    int error_code = VecDestroy(&petsc_vector_);
89 90 91 92 93
                    assert(!error_code && "Error in Vec destruction."); // no exception in destructors!
                    static_cast<void>(error_code); // to avoid arning in release compilation.
                }
            }
            
94
            
95 96 97 98 99 100
            void Swap(Wrappers::Petsc::Vector& A, Wrappers::Petsc::Vector& B)
            {
                using std::swap;
                swap(A.do_petsc_destroy_, B.do_petsc_destroy_);
                swap(A.petsc_vector_, B.petsc_vector_);
            }
101 102
            
            
103 104 105
            void Vector::InitSequentialVector(unsigned int size,
                                              const Mpi& mpi,
                                              const char* invoking_file, int invoking_line)
106
            {
107
                assert(petsc_vector_ == PETSC_NULL && "Should not be initialized when this method is called!");
108
                int error_code = VecCreateSeq(mpi.GetCommunicator(), static_cast<PetscInt>(size), &petsc_vector_);
109
                if (error_code)
110
                    throw ExceptionNS::Exception(error_code, "VecCreateSeq", invoking_file, invoking_line);
111 112
                
                do_petsc_destroy_ = true;
113 114 115
            }
            
            
116
            void Vector::InitMpiVector(unsigned int local_size, unsigned int global_size,
117
                                       const Mpi& mpi, const char* invoking_file, int invoking_line)
118
            {
119
                assert(petsc_vector_ == PETSC_NULL && "Should not be initialized when this method is called!");
120
                assert(local_size <= global_size && "If not, either sequential mode or the local and global were "
121 122
                       "provided in the wrong order.");
                
123
                int error_code = VecCreateMPI(mpi.GetCommunicator(),
124 125
                                              static_cast<PetscInt>(local_size),
                                              static_cast<PetscInt>(global_size),
126
                                              &petsc_vector_);
127 128
                
                if (error_code)
129
                    throw ExceptionNS::Exception(error_code, "VecCreateMPI", invoking_file, invoking_line);
130 131
                
                do_petsc_destroy_ = true;
132 133 134
            }
            
            
135
            void Vector::InitMpiVectorWithGhost(unsigned int local_size, unsigned int global_size,
136
                                                const std::vector<PetscInt>& ghost_padding,
137
                                                const Mpi& mpi, const char* invoking_file, int invoking_line)
138
            {
139
                assert(petsc_vector_ == PETSC_NULL && "Should not be initialized when this method is called!");
140
                assert(local_size <= global_size && "If not, either sequential mode or the local and global were "
141
                       "provided in the wrong order.");
142
                
143
                const PetscInt Nghost = static_cast<PetscInt>(ghost_padding.size());
144
                
145
                int error_code = VecCreateGhost(mpi.GetCommunicator(),
146 147
                                                static_cast<PetscInt>(local_size),
                                                static_cast<PetscInt>(global_size),
148
                                                Nghost,
149
                                                ghost_padding.data(),
150
                                                &petsc_vector_);
151 152
                
                if (error_code)
153
                    throw ExceptionNS::Exception(error_code, "VecCreateGhost", invoking_file, invoking_line);
154
               
155
                do_petsc_destroy_ = true;
156
            }
157 158 159 160 161 162 163 164
            
            
            void Vector::InitSequentialFromFile(const std::string& file,
                                                const Mpi& mpi,
                                                const char* invoking_file, int invoking_line)
            {
                assert(petsc_vector_ == PETSC_NULL && "Should not be initialized when this method is called!");
                assert(mpi.Nprocessor<int>() == 1 && "This method assumes sequential case!");
165
                 
166
                std::ifstream stream;
167
                FilesystemNS::File::Read(stream, file, invoking_file, invoking_line);
168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188
                
                std::string line;
                std::vector<double> value_list;
                
                while (getline(stream, line))
                {
                    // All lines are expected to be one value.
                   value_list.push_back(stod(line));
                }
                
                // Now Init the vector with the appropriate size.
                const unsigned int Nvalue = static_cast<unsigned int>(value_list.size());
                
                InitSequentialVector(Nvalue, mpi, invoking_file, invoking_line);
                
                // And fill it with the values.
                AccessVectorContent<Utilities::Access::read_and_write> content(*this, invoking_file, invoking_line);
                
                for (auto i = 0u; i < Nvalue; ++i)
                    content[i] = value_list[static_cast<std::size_t>(i)];
            }
189 190
            
            
191
            void Vector::DuplicateLayout(const Vector& original, const char* invoking_file, int invoking_line)
192
            {
193
                assert(petsc_vector_ == PETSC_NULL && "Should not be initialized when this method is called!");
194
                int error_code = VecDuplicate(original.Internal(), &petsc_vector_);
195
                if (error_code)
196
                    throw ExceptionNS::Exception(error_code, "VecDuplicate", invoking_file, invoking_line);
197

198
                do_petsc_destroy_ = original.do_petsc_destroy_;
199 200 201
            }
            
            
202
            PetscInt Vector::GetProcessorWiseSize(const char* invoking_file, int invoking_line) const
203 204
            {
                PetscInt ret;
205
                int error_code = VecGetLocalSize(Internal(), &ret);
206 207
                
                if (error_code)
208
                    throw ExceptionNS::Exception(error_code, "VecGetLocalSize", invoking_file, invoking_line);
209 210 211 212 213
                
                return ret;
            }
            
            
214
            PetscInt Vector::GetProgramWiseSize(const char* invoking_file, int invoking_line) const
215 216
            {
                PetscInt ret;
217
                int error_code = VecGetSize(Internal(), &ret);
218 219
                
                if (error_code)
220
                    throw ExceptionNS::Exception(error_code, "VecGetSize", invoking_file, invoking_line);
221 222 223 224 225
                
                return ret;
            }
            
            
226
            void Vector::ZeroEntries(const char* invoking_file, int invoking_line)
227
            {
228
                int error_code = VecZeroEntries(Internal());
229
                if (error_code)
230
                    throw ExceptionNS::Exception(error_code, "VecZeroEntries", invoking_file, invoking_line);
231 232 233
            }
            
            
234 235 236 237 238 239 240 241
            void Vector::Set(const PetscScalar alpha, const char* invoking_file, int invoking_line)
            {
                int error_code = VecSet(Internal(), alpha);
                if (error_code)
                    throw ExceptionNS::Exception(error_code, "VecSet", invoking_file, invoking_line);
            }
            
            
242 243
            void Vector::SetValues(const std::vector<PetscInt>& indexing,
                                   const PetscScalar* values,
244
                                   InsertMode insertOrAppend, const char* invoking_file, int invoking_line)
245
            {
246 247 248 249 250
                int error_code = VecSetValues(Internal(),
                                              static_cast<PetscInt>(indexing.size()),
                                              indexing.data(),
                                              values,
                                              insertOrAppend);
251
                               
252
                if (error_code)
253
                    throw ExceptionNS::Exception(error_code, "VecSetValues", invoking_file, invoking_line);
254 255 256 257
            }
            
            
            void Vector::SetValue(PetscInt index, PetscScalar value,
258
                                  InsertMode insertOrAppend, const char* invoking_file, int invoking_line)
259
            {
260
                int error_code = VecSetValue(Internal(), index, value, insertOrAppend);
261 262
                
                if (error_code)
263
                    throw ExceptionNS::Exception(error_code, "VecSetValue", invoking_file, invoking_line);
264 265 266
            }
            
            
267 268 269 270 271 272 273 274
            void Vector::SetUniformValue(PetscScalar value, const char* invoking_file, int invoking_line)
            {
                int error_code = VecSet(Internal(), value);
                
                if (error_code)
                    throw ExceptionNS::Exception(error_code, "VecSetValue", invoking_file, invoking_line);
            }

275
            
276
            std::vector<PetscScalar> Vector::GetValues(const std::vector<PetscInt>& indexing,
277
                                                       const char* invoking_file, int invoking_line) const
278 279 280
            {
                std::vector<PetscScalar> ret(indexing.size());
                
281
                int error_code = VecGetValues(Internal(), static_cast<PetscInt>(indexing.size()), indexing.data(), ret.data());
282 283
                
                if (error_code)
284
                    throw ExceptionNS::Exception(error_code, "VecGetValues", invoking_file, invoking_line);
285
                
286
                return ret;
287
            }
288
            
289 290
            
            PetscScalar Vector::GetValue(PetscInt index,
291
                                         const char* invoking_file, int invoking_line) const
292 293 294
            {
                std::vector<PetscInt> index_as_vector { index };
                
295
                const auto ret = GetValues(index_as_vector, invoking_file, invoking_line);
296 297 298
                assert(ret.size() == 1);
                return ret[0];
            }
299 300
            
            
301
            void Vector::Assembly(const char* invoking_file, int invoking_line,
302
                                  update_ghost do_update_ghost)
303
            {
304
                int error_code = VecAssemblyBegin(Internal());
305
                if (error_code)
306
                    throw ExceptionNS::Exception(error_code, "VecAssemblyBegin", invoking_file, invoking_line);
307 308
                
                
309
                error_code = VecAssemblyEnd(Internal());
310
                if (error_code)
311
                    throw ExceptionNS::Exception(error_code, "VecAssemblyEnd", invoking_file, invoking_line);
312
                
313
                UpdateGhosts(invoking_file, invoking_line, do_update_ghost);
314 315 316
            }
            
            
317 318
            void Vector::Copy(const Vector& source, const char* invoking_file, int invoking_line,
                              update_ghost do_update_ghost)
319
            {
320 321
                int error_code = VecCopy(source.Internal(),
                                         Internal());
322
                if (error_code)
323
                    throw ExceptionNS::Exception(error_code, "VecCopy", invoking_file, invoking_line);
324
                
325
                UpdateGhosts(invoking_file, invoking_line, do_update_ghost);
326 327 328
            }
            
            
329 330
            void Vector::CompleteCopy(const Vector& source, const char* invoking_file, int invoking_line,
                                      update_ghost do_update_ghost)
331
            {
332
                DuplicateLayout(source, invoking_file, invoking_line);
333
                Copy(source, invoking_file, invoking_line, do_update_ghost);
334 335 336
            }
            
            
337
            void Vector::Scale(PetscScalar a, const char* invoking_file, int invoking_line)
338
            {
339
                int error_code = VecScale(Internal(), a);
340
                if (error_code)
341
                    throw ExceptionNS::Exception(error_code, "VecScale", invoking_file, invoking_line);
342 343
            }
            
344
            
345
            void Vector::UpdateGhosts(const char* invoking_file, int invoking_line)
346
            {
347
                int error_code = VecGhostUpdateBegin(Internal(), INSERT_VALUES, SCATTER_FORWARD);
348
                if (error_code)
349
                    throw ExceptionNS::Exception(error_code, "VecGhostUpdateBegin", invoking_file, invoking_line);
350
                
351
                error_code = VecGhostUpdateEnd(Internal(), INSERT_VALUES, SCATTER_FORWARD);
352
                if (error_code)
353
                    throw ExceptionNS::Exception(error_code, "VecGhostUpdateEnd", invoking_file, invoking_line);
354
            }
355 356
            
            
357
            void DisplaySomeValues(std::ostream& stream, const Vector& vector, PetscInt firstIndex, PetscInt lastIndex,
358
                                   int rankProc, const char* invoking_file, int invoking_line)
359
            {
360
                AccessVectorContent<Utilities::Access::read_only> local_array(vector, invoking_file, invoking_line);
361
                
362
                #ifndef NDEBUG
363
                const PetscInt size = vector.GetProcessorWiseSize(invoking_file, invoking_line);
364 365
                assert(local_array.GetSize(__FILE__, __LINE__) == static_cast<unsigned int>(size));

366 367
                assert(firstIndex < lastIndex);
                assert(lastIndex < size);
368
                #endif // NDEBUG
369
                stream << "On processor " << rankProc << " [";
370
                auto values = local_array.GetArray();
371
                
372 373
                for (PetscInt i = firstIndex; i <= lastIndex; ++i)
                    stream << values[i] << ", ";
374
                
375 376
                stream << ']' << std::endl;
            }
377 378
            
            
379
            bool AreEqual(const Vector& vec1, const Vector& vec2, const double epsilon, 
GILLES Sebastien's avatar
GILLES Sebastien committed
380
                          std::string& inequality_description,
381
                          const char* invoking_file, int invoking_line)
382
            {
GILLES Sebastien's avatar
GILLES Sebastien committed
383
                inequality_description.clear();
384
                
385
                PetscInt size = vec1.GetProcessorWiseSize(invoking_file, invoking_line);
386
                {
387
                    #ifndef NDEBUG
388
                    PetscInt size2 = vec2.GetProcessorWiseSize(invoking_file, invoking_line);
389
                    assert(size == size2);
390
                    #endif // NDEBUG
391
                }
392
                
393 394
                AccessVectorContent<Utilities::Access::read_only> local_array_1(vec1, invoking_file, invoking_line);
                AccessVectorContent<Utilities::Access::read_only> local_array_2(vec2, invoking_file, invoking_line);
395 396 397

                assert(static_cast<PetscInt>(local_array_1.GetSize(__FILE__, __LINE__)) == size);
                assert(static_cast<PetscInt>(local_array_2.GetSize(__FILE__, __LINE__)) == size);                
398
                
399 400
                const PetscScalar* values1 = local_array_1.GetArray();
                const PetscScalar* values2 = local_array_2.GetArray();
401
                
402
                bool ret = true;
403
                
404
                for (PetscInt index = 0; index < size && ret;)
405 406 407 408 409 410 411
                {
                    {
                        if (std::fabs(values1[index] - values2[index]) > epsilon)
                        {
                            std::ostringstream oconv;
                            oconv << "Inequality found for index " << index << ": vector 1 displays " << values1[index] << " while "
                            "vector2 displays " << values2[index] << std::endl;
GILLES Sebastien's avatar
GILLES Sebastien committed
412
                            inequality_description = oconv.str();
413 414 415 416 417 418
                            ret = false;
                        }
                        else
                            ++index;
                    }
                }
419
                
420 421
                return ret;
            }
422 423
            
            
424 425 426 427 428
            void AXPY(PetscScalar alpha,
                      const Vector& x,
                      Vector& y,
                      const char* invoking_file, int invoking_line,
                      update_ghost do_update_ghost)
429
            {
430
                int error_code = VecAXPY(y.Internal(), alpha, x.Internal());
431
                if (error_code)
432
                    throw ExceptionNS::Exception(error_code, "VecAXPY", invoking_file, invoking_line);
433 434
                
                y.UpdateGhosts(invoking_file, invoking_line, do_update_ghost);
435
            }
436
            
437 438 439 440 441 442 443 444
            
            void DotProduct(const Vector& x, const Vector& y, PetscScalar& val, const char* invoking_file, int invoking_line)
            {
                int error_code = VecDot(x.Internal(), y.Internal(), &val);
                if (error_code)
                    throw ExceptionNS::Exception(error_code, "VecDot", invoking_file, invoking_line);
            }
            
445
                        
446
            std::pair<PetscInt, PetscReal> Vector::Min(const char* invoking_file, int invoking_line) const
447 448 449
            {
                PetscInt position;
                PetscReal value;
450
                
451
                int error_code = VecMin(Internal(), &position, &value);
452
                if (error_code)
453
                    throw ExceptionNS::Exception(error_code, "VecMin", invoking_file, invoking_line);
454
                
455 456
                return std::make_pair(position, value);
            }
457 458
            
            
459
            std::pair<PetscInt, PetscReal> Vector::Max(const char* invoking_file, int invoking_line) const
460 461 462
            {
                PetscInt position;
                PetscReal value;
463
                
464
                int error_code = VecMax(Internal(), &position, &value);
465
                if (error_code)
466
                    throw ExceptionNS::Exception(error_code, "VecMax", invoking_file, invoking_line);
467
                
468 469
                return std::make_pair(position, value);
            }
470 471
            
            
472
            double Vector::Norm(NormType type, const char* invoking_file, int invoking_line) const
473 474 475
            {
                PetscReal norm;
                
476 477
                assert(type == NORM_1 || type == NORM_2 || type == NORM_INFINITY);
                
478 479 480 481 482 483 484 485
                int error_code = VecNorm(Internal(), type, &norm);
                if (error_code)
                    throw ExceptionNS::Exception(error_code, "VecNorm", invoking_file, invoking_line);
                
                return static_cast<double>(norm);
            }
            
            
486
            
487
            void Vector::View(const Mpi& mpi, const char* invoking_file, int invoking_line) const
GILLES Sebastien's avatar
GILLES Sebastien committed
488
            {
489
                int error_code = VecView(Internal(), PETSC_VIEWER_STDOUT_(mpi.GetCommunicator()));
GILLES Sebastien's avatar
GILLES Sebastien committed
490
                if (error_code)
491
                    throw ExceptionNS::Exception(error_code, "VecView", invoking_file, invoking_line);
GILLES Sebastien's avatar
GILLES Sebastien committed
492 493 494
            }
            
            
GILLES Sebastien's avatar
GILLES Sebastien committed
495 496
            void Vector::View(const Mpi& mpi, const std::string& output_file,
                              const char* invoking_file, int invoking_line,
497
                              PetscViewerFormat format) const
498
            {
499
                Viewer viewer(mpi, output_file, format, invoking_file, invoking_line);
500
                
501
                int error_code = VecView(Internal(), viewer.GetUnderlyingPetscObject());
502
                if (error_code)
503
                    throw ExceptionNS::Exception(error_code, "VecView", invoking_file, invoking_line);
504 505 506
            }
            
            
507 508 509 510 511 512 513 514 515 516 517 518 519 520
            void Vector::ViewBinary(const Mpi& mpi,
                                    const std::string& output_file,
                                    const char* invoking_file, int invoking_line) const
            {
                Viewer viewer(mpi, output_file, FILE_MODE_WRITE, invoking_file, invoking_line);
                
                int error_code = VecView(Internal(), viewer.GetUnderlyingPetscObject());
                if (error_code)
                    throw ExceptionNS::Exception(error_code, "VecView", invoking_file, invoking_line);
                
            }

            
            
521 522 523
            void GatherVector(const Mpi& mpi,
                              const Wrappers::Petsc::Vector& local_parallel_vector,
                              Wrappers::Petsc::Vector& sequential_vector,
524
                              const char* invoking_file, int invoking_line)
525
            {
526
                if (mpi.Nprocessor<int>() == 1)
527 528 529 530 531 532 533 534 535 536 537
                    return;
                
                VecScatter vecscat;
                
                Vec buf = sequential_vector.Internal();
                
                Vec local_petsc_vector = local_parallel_vector.Internal();
                Vec sequential_petsc_vector = sequential_vector.Internal();
                
                int error_code = VecScatterCreateToAll(local_petsc_vector, &vecscat, &buf);
                if (error_code)
538
                    throw ExceptionNS::Exception(error_code, "VecScatterCreateToAll", invoking_file, invoking_line);
539
                
GILLES Sebastien's avatar
GILLES Sebastien committed
540 541
                error_code = VecScatterBegin(vecscat, local_petsc_vector, sequential_petsc_vector, INSERT_VALUES,
                                             SCATTER_FORWARD);
542
                if (error_code)
543
                    throw ExceptionNS::Exception(error_code, "VecScatterBegin", invoking_file, invoking_line);
544
                
GILLES Sebastien's avatar
GILLES Sebastien committed
545 546
                error_code = VecScatterEnd(vecscat, local_petsc_vector, sequential_petsc_vector, INSERT_VALUES,
                                           SCATTER_FORWARD);
547
                if (error_code)
548
                    throw ExceptionNS::Exception(error_code, "VecScatterEnd", invoking_file, invoking_line);
549 550 551
                
                error_code = VecScatterDestroy(&vecscat);
                if (error_code)
552
                    throw ExceptionNS::Exception(error_code, "VecScatterDestroy", invoking_file, invoking_line);
553 554
            }
            
555
            
556 557 558 559
            void Vector::SetDoNotDestroyPetscVector()
            {
                do_petsc_destroy_ = false;
            }
560 561
            
            
562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580
            double ScalarProduct(const Mpi& mpi,
                                 const Vector& v1,
                                 const Vector& v2,
                                 const char* invoking_file, int invoking_line)
            {
                AccessVectorContent<Utilities::Access::read_only> content1(v1, invoking_file, invoking_line);
                AccessVectorContent<Utilities::Access::read_only> content2(v2, invoking_file, invoking_line);
                
                const auto size = content1.GetSize(invoking_file, invoking_line);
                assert(size == content2.GetSize(invoking_file, invoking_line));
                
                double ret = 0.;
                
                for (auto i = 0u; i < size; ++i)
                    ret += content1.GetValue(i) * content2.GetValue(i);
                    
                // Query all processors
                return mpi.AllReduce(ret, MpiNS::Op::Sum);
            }
581 582
            
            
583 584 585 586 587 588 589 590 591 592 593
            double NormL2(const Mpi& mpi,
                          const Vector& v,
                          const char* invoking_file, int invoking_line)
            {
                AccessVectorContent<Utilities::Access::read_only> content(v, invoking_file, invoking_line);
                
                const auto size = content.GetSize(invoking_file, invoking_line);
                
                double ret = 0.;
                
                for (auto i = 0u; i < size; ++i)
594
                    ret += NumericNS::Square(content.GetValue(i));
595 596 597 598 599 600 601
                
                // Query all processors
                return std::sqrt(mpi.AllReduce(ret, MpiNS::Op::Sum));
            }

            
            
602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617
            void Vector::SetFromPetscVec(const Vec& petsc_vector,
                                         const char* invoking_file, int invoking_line)
            {
                // In this specific method alone I can't use syntax sugary provided by the class, as I manipulate
                // a raw Petsc Vec such as the one provided in arguments of SnesFunction.
                const PetscScalar* values;
                int error_code = VecGetArrayRead(petsc_vector, &values);
                
                if (error_code)
                    throw ExceptionNS::Exception(error_code, "VecGetArrayRead", invoking_file, invoking_line);
                
                const auto size = static_cast<std::size_t>(GetProcessorWiseSize(invoking_file, invoking_line));
                
                #ifndef NDEBUG
                {
                    PetscInt petsc_size;
618
                    error_code = VecGetLocalSize(petsc_vector, &petsc_size);
619 620
                
                    if (error_code)
621 622
                        throw ExceptionNS::Exception(error_code, "VecGetLocalSize", invoking_file, invoking_line);
                    
623 624 625 626 627 628 629 630 631 632 633 634 635 636 637
                    assert(size == static_cast<std::size_t>(petsc_size));
                }
                #endif // NDEBUG
                
                AccessVectorContent<Utilities::Access::read_and_write> content(*this, invoking_file, invoking_line);
                
                for (auto i = 0ul; i < size; ++i)
                    content[static_cast<unsigned int>(i)] = values[i];

                
                error_code = VecRestoreArrayRead(petsc_vector, &values);
                
                if (error_code)
                    throw ExceptionNS::Exception(error_code, "VecRestoreArrayRead", invoking_file, invoking_line);
            }
638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656
            
            
            #ifndef NDEBUG
            void AssertNumericValues(const Vector& vector,
                                     const char* invoking_file, int invoking_line)
            {
                AccessVectorContent<Utilities::Access::read_only> content(vector, invoking_file, invoking_line);
                
                const auto size = content.GetSize(invoking_file, invoking_line);

                for (unsigned int i = 0u; i < size; ++i)
                {
                    const auto current_value = content.GetValue(i);
                    
                    assert(!std::isnan(current_value));
                    assert(!std::isinf(current_value));
                }
            }
            #endif // NDEBUG
657 658 659

            
            
660 661 662 663 664 665 666
        } // namespace Petsc
        
        
    } // namespace Wrappers
    
    
} // namespace HappyHeart
667