Matrix.cpp 31.5 KB
Newer Older
GILLES Sebastien's avatar
GILLES Sebastien committed
1 2 3 4 5 6 7 8 9 10 11 12 13
/*!
//
// \file
//
//
// Created by Sebastien Gilles <sebastien.gilles@inria.fr> on the Fri, 4 Oct 2013 11:24:18 +0200
// Copyright (c) Inria. All rights reserved.
//
// \ingroup ThirdPartyGroup
// \addtogroup ThirdPartyGroup
// \{
*/

14 15 16

#include <cassert>
#include <cmath>
17

18
#include "ThirdParty/IncludeWithoutWarning/Petsc/PetscMatPrivate.h"
19

20 21 22 23 24
#include "ThirdParty/Wrappers/Petsc/Matrix/Matrix.hpp"
#include "ThirdParty/Wrappers/Petsc/Matrix/MatrixPattern.hpp"
#include "ThirdParty/Wrappers/Petsc/Vector/Vector.hpp"
#include "ThirdParty/Wrappers/Petsc/Viewer.hpp"
#include "ThirdParty/Wrappers/Petsc/Exceptions/Petsc.hpp"
25

26

27
namespace MoReFEM
28
{
29 30
    
    
31 32
    namespace Wrappers
    {
33 34
        
        
35 36
        namespace Petsc
        {
37 38 39
            
            
            Matrix::Matrix()
40 41
            : parent(),
            petsc_matrix_(PETSC_NULL),
42
            do_petsc_destroy_(false)
43 44 45
            { }
            
            
46 47 48 49 50 51 52 53 54 55 56
            Matrix::Matrix(const Mpi& mpi,
                           const std::string& binary_file,
                           const char* invoking_file, int invoking_line)
            : parent(),
            petsc_matrix_(PETSC_NULL),
            do_petsc_destroy_(true)
            {
                const auto communicator = mpi.GetCommunicator();
                
                Viewer viewer(mpi,
                              binary_file,
57
                              PETSC_VIEWER_BINARY_MATLAB,
58 59 60 61 62 63 64 65 66 67 68 69 70
                              FILE_MODE_READ,
                              invoking_file, invoking_line);
                
                int error_code = MatCreate(communicator, &petsc_matrix_);
                if (error_code)
                    throw ExceptionNS::Exception(error_code, "MatCreate", invoking_file, invoking_line);
                
                error_code = MatLoad(petsc_matrix_, viewer.GetUnderlyingPetscObject());
                if (error_code)
                    throw ExceptionNS::Exception(error_code, "MatLoad", invoking_file, invoking_line);
            }
            
            
71
            Matrix::Matrix(const Matrix& rhs)
72 73
            : parent(rhs),
            petsc_matrix_(PETSC_NULL),
74
            do_petsc_destroy_(rhs.do_petsc_destroy_)
75 76 77 78 79
            {
                CompleteCopy(rhs, __FILE__, __LINE__);
            }
            
            
80 81
            Matrix::~Matrix()
            {
82
                if (do_petsc_destroy_)
83
                {
84
                    assert(petsc_matrix_ != PETSC_NULL);
85 86 87 88 89 90
                    int error_code = MatDestroy(&petsc_matrix_);
                    assert(!error_code && "Error in Mat destruction."); // no exception in destructors!
                    static_cast<void>(error_code); // to avoid arning in release compilation.
                }
            }
            
91 92 93 94 95 96 97 98
            
            void Swap(Matrix& A, Matrix& B)
            {
                using std::swap;
                swap(A.do_petsc_destroy_, B.do_petsc_destroy_);
                swap(A.petsc_matrix_, B.petsc_matrix_);
            }
            
99
            
100
            void Matrix::DuplicateLayout(const Matrix& original, MatDuplicateOption option, const char* invoking_file, int invoking_line)
101
            {
102 103
                assert(petsc_matrix_ == PETSC_NULL && "Should not be initialized when this method is called!");

104 105
                int error_code = MatDuplicate(original.Internal(), option, &petsc_matrix_);
                if (error_code)
106
                    throw ExceptionNS::Exception(error_code, "MatDuplicate", invoking_file, invoking_line);
107

GILLES Sebastien's avatar
GILLES Sebastien committed
108
		        do_petsc_destroy_ = original.do_petsc_destroy_;                
109 110 111
            }
            
            
112
            void Matrix::InitSequentialMatrix(unsigned int Nrow, unsigned int Ncolumn,
113
                                              const MatrixPattern& matrix_pattern,
114
                                              const Mpi& mpi, const char* invoking_file, int invoking_line)
115
            {
116
                assert(matrix_pattern.Nrow() == Nrow);
117
                const auto& communicator = mpi.GetCommunicator();
118
                
119 120
                assert(petsc_matrix_ == PETSC_NULL && "Should not be initialized when this method is called!");
                
121
                // Follow Petsc's advice and build the matrix step by step (see MatCreateSeqAIJ for this advice).
122
                int error_code = MatCreate(communicator, &petsc_matrix_);
123
                if (error_code)
124
                    throw ExceptionNS::Exception(error_code, "MatCreate", invoking_file, invoking_line);
125
                
126
                const Mat& internal = Internal(); // const is indeed ignored below as Mat is truly a pointer...
127
                                
128
                error_code = MatSetType(internal, MATSEQAIJ);
129
                if (error_code)
130
                    throw ExceptionNS::Exception(error_code, "MatSetType", invoking_file, invoking_line);
131 132
                
                
133 134
                const PetscInt Nrow_int = static_cast<PetscInt>(Nrow);
                const PetscInt Ncol_int = static_cast<PetscInt>(Ncolumn);
135
                
136
                error_code = MatSetSizes(internal, Nrow_int, Ncol_int, Nrow_int, Ncol_int);
137
                if (error_code)
138
                    throw ExceptionNS::Exception(error_code, "MatSetSizes", invoking_file, invoking_line);
139 140
                
                
141 142
                // Preallocate the sparse positions. Petsc indicates it speeds up tremendously the creation
                // process. This must be called AFTER size of the matrix is set.
143 144 145 146 147
                error_code = MatSeqAIJSetPreallocationCSR(Internal(),
                                                          matrix_pattern.GetICsr().data(),
                                                          matrix_pattern.GetJCsr().data(),
                                                          PETSC_NULL);
                
148
                if (error_code)
149
                    throw ExceptionNS::Exception(error_code, "MatSeqAIJSetPreallocation", invoking_file, invoking_line);
150 151
                
                
152
                error_code = MatSetOption(internal, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);
153
                if (error_code)
154
                    throw ExceptionNS::Exception(error_code, "MatSetOption", invoking_file, invoking_line);
155 156
                
                do_petsc_destroy_ = true;
157
            }
158 159 160
            
            
            
161 162
            void Matrix::InitParallelMatrix(unsigned int Nlocal_row, unsigned int Nlocal_column,
                                            unsigned int Nglobal_row, unsigned int Nglobal_column,
163
                                            const MatrixPattern& matrix_pattern,
164
                                            const Mpi& mpi, const char* invoking_file, int invoking_line)
165
            {
166
                assert(matrix_pattern.Nrow() == Nlocal_row);
167

168 169
                assert(petsc_matrix_ == PETSC_NULL && "Should not be initialized when this method is called!");
                
170
                const auto& communicator = mpi.GetCommunicator();
171
                
172
                // Follow Petsc's advice and build the matrix step by step (see MatCreateSeqAIJ for this advice).
173
                int error_code = MatCreate(communicator, &petsc_matrix_);
174
                if (error_code)
175
                    throw ExceptionNS::Exception(error_code, "MatCreate", invoking_file, invoking_line);
176
                
177 178 179
                const Mat& internal = Internal(); // const is indeed ignored below as Mat is truly a pointer...
                
                error_code = MatSetType(internal, MATAIJ);
180
                if (error_code)
181
                    throw ExceptionNS::Exception(error_code, "MatSetType", invoking_file, invoking_line);
182
                
183 184 185 186
                const PetscInt Nlocal_row_int = static_cast<PetscInt>(Nlocal_row);
                const PetscInt Nlocal_col_int = static_cast<PetscInt>(Nlocal_column);
                const PetscInt Nglobal_row_int = static_cast<PetscInt>(Nglobal_row);
                const PetscInt Nglobal_col_int = static_cast<PetscInt>(Nglobal_column);
187
                                
188
                error_code = MatSetSizes(internal, Nlocal_row_int, Nlocal_col_int, Nglobal_row_int, Nglobal_col_int);
189
                if (error_code)
190
                    throw ExceptionNS::Exception(error_code, "MatSetSizes", invoking_file, invoking_line);
191 192
                
                
193
                // Preallocate the sparse positions. Petsc indicates it speeds up tremendously the creation
194 195 196 197 198 199
                // process. This must be called AFTER size of the matrix is set.               
                error_code = MatMPIAIJSetPreallocationCSR(internal,
                                                          matrix_pattern.GetICsr().data(),
                                                          matrix_pattern.GetJCsr().data(),
                                                          PETSC_NULL);
                
200
                if (error_code)
201
                    throw ExceptionNS::Exception(error_code, "MatMPIAIJSetPreallocation", invoking_file, invoking_line);                
202
                
203
                error_code = MatSetOption(internal, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);
204
                if (error_code)
205
                    throw ExceptionNS::Exception(error_code, "MatSetOption", invoking_file, invoking_line);
206
                
207 208 209 210
                error_code = MatSetOption(internal, MAT_NO_OFF_PROC_ZERO_ROWS, PETSC_TRUE);
                if (error_code)
                    throw ExceptionNS::Exception(error_code, "MatSetOption", invoking_file, invoking_line);
                
211
                do_petsc_destroy_ = true;
212
            }
213
            
214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233
            
            void Matrix::InitSequentialDenseMatrix(unsigned int Nrow, unsigned int Ncolumn,
                                                   const Mpi& mpi, const char* invoking_file, int invoking_line)
            {
                const auto& communicator = mpi.GetCommunicator();
                
                assert(petsc_matrix_ == PETSC_NULL && "Should not be initialized when this method is called!");
                
                // Follow Petsc's advice and build the matrix step by step (see MatCreateSeqAIJ for this advice).
                int error_code = MatCreate(communicator, &petsc_matrix_);
                if (error_code)
                    throw ExceptionNS::Exception(error_code, "MatCreate", invoking_file, invoking_line);
                
                const Mat& internal = Internal(); // const is indeed ignored below as Mat is truly a pointer...
                
                error_code = MatSetType(internal, MATSEQDENSE);
                if (error_code)
                    throw ExceptionNS::Exception(error_code, "MatSetType", invoking_file, invoking_line);
                
                
234 235
                const PetscInt Nrow_int = static_cast<PetscInt>(Nrow);
                const PetscInt Ncol_int = static_cast<PetscInt>(Ncolumn);
236 237 238 239 240 241 242 243 244 245 246
                
                error_code = MatSetSizes(internal, Nrow_int, Ncol_int, Nrow_int, Ncol_int);
                if (error_code)
                    throw ExceptionNS::Exception(error_code, "MatSetSizes", invoking_file, invoking_line);
                
                error_code = MatSetUp(internal);
                if (error_code)
                    throw ExceptionNS::Exception(error_code, "MatSetUp", invoking_file, invoking_line);
                
                do_petsc_destroy_ = true;
            }
247
            
248

249 250 251
            void Matrix::InitParallelDenseMatrix(unsigned int Nlocal_row, unsigned int Nlocal_column,
                                                 unsigned int Nglobal_row, unsigned int Nglobal_column,
                                                 const Mpi& mpi, const char* invoking_file, int invoking_line)
252 253 254 255 256 257 258 259 260 261 262 263
            {
                const auto& communicator = mpi.GetCommunicator();
                
                assert(petsc_matrix_ == PETSC_NULL && "Should not be initialized when this method is called!");
                
                // Follow Petsc's advice and build the matrix step by step (see MatCreateSeqAIJ for this advice).
                int error_code = MatCreate(communicator, &petsc_matrix_);
                if (error_code)
                    throw ExceptionNS::Exception(error_code, "MatCreate", invoking_file, invoking_line);
                
                const Mat& internal = Internal(); // const is indeed ignored below as Mat is truly a pointer...
                
264
                error_code = MatSetType(internal, MATMPIDENSE);
265 266 267 268 269 270 271 272 273
                if (error_code)
                    throw ExceptionNS::Exception(error_code, "MatSetType", invoking_file, invoking_line);
                
                const PetscInt Nlocal_row_int = static_cast<PetscInt>(Nlocal_row);
                const PetscInt Nlocal_col_int = static_cast<PetscInt>(Nlocal_column);
                const PetscInt Nglobal_row_int = static_cast<PetscInt>(Nglobal_row);
                const PetscInt Nglobal_col_int = static_cast<PetscInt>(Nglobal_column);
                
                error_code = MatSetSizes(internal, Nlocal_row_int, Nlocal_col_int, Nglobal_row_int, Nglobal_col_int);
274

275 276 277 278 279 280 281 282 283 284
                if (error_code)
                    throw ExceptionNS::Exception(error_code, "MatSetSizes", invoking_file, invoking_line);
                
                error_code = MatSetUp(internal);
                if (error_code)
                    throw ExceptionNS::Exception(error_code, "MatSetUp", invoking_file, invoking_line);
                
                do_petsc_destroy_ = true;
            }

285
            
286
            void Matrix::ZeroEntries(const char* invoking_file, int invoking_line)
287
            {
288
                int error_code = MatZeroEntries(Internal());
289
                if (error_code)
290
                    throw ExceptionNS::Exception(error_code, "MatZeroEntries", invoking_file, invoking_line);
291
            }
292 293 294 295
            
            
            void Matrix::ZeroRows(const std::vector<PetscInt>& row_indexes,
                                  PetscScalar diagonal_value,
296
                                  const char* invoking_file, int invoking_line)
297
            {
298
                int error_code = MatZeroRows(Internal(),
299 300 301 302 303
                                             static_cast<PetscInt>(row_indexes.size()),
                                             row_indexes.data(),
                                             diagonal_value,
                                             PETSC_NULL,
                                             PETSC_NULL);
304
                
305
                if (error_code)
306
                    throw ExceptionNS::Exception(error_code, "MatZeroRows", invoking_file, invoking_line);
307
            }
308 309
            
            
310 311 312 313 314 315 316 317 318 319
            void Matrix::ZeroRowsColumns(const std::vector<PetscInt>& row_indexes,
                                         PetscScalar diagonal_value,
                                         const char* invoking_file, int invoking_line)
            {
                int error_code = MatZeroRowsColumns(Internal(),
                                                    static_cast<PetscInt>(row_indexes.size()),
                                                    row_indexes.data(),
                                                    diagonal_value,
                                                    PETSC_NULL,
                                                    PETSC_NULL);
320
                
321
                if (error_code)
322
                    throw ExceptionNS::Exception(error_code, "MatZeroRowsColumns", invoking_file, invoking_line);
323 324
            }
            
325 326
            
            
327
            void Matrix::Assembly(const char* invoking_file, int invoking_line)
328
            {
329 330 331
                const Mat& internal = Internal(); // const is indeed ignored below as Mat is truly a pointer...
                
                int error_code = MatAssemblyBegin(internal, MAT_FINAL_ASSEMBLY);
332
                if (error_code)
333
                    throw ExceptionNS::Exception(error_code, "MatAssemblyBegin", invoking_file, invoking_line);
334
                
335
                error_code = MatAssemblyEnd(internal, MAT_FINAL_ASSEMBLY);
336
                if (error_code)
337
                    throw ExceptionNS::Exception(error_code, "MatAssemblyEnd", invoking_file, invoking_line);
338
            }
339 340
            
            
341 342 343 344 345 346 347 348 349 350 351 352 353 354
            bool Matrix::IsAssembled(const char* invoking_file, int invoking_line) const
            {
                const Mat& internal = Internal();
                
                PetscBool value;
                
                int error_code = MatAssembled(internal, &value);
                if (error_code)
                    throw ExceptionNS::Exception(error_code, "MatAssembled", invoking_file, invoking_line);
                
                return (value == PETSC_TRUE);
            }
            
            
355 356 357
            void Matrix::SetValues(const std::vector<PetscInt>& row_indexing,
                                   const std::vector<PetscInt>& column_indexing,
                                   const PetscScalar* values, InsertMode insertOrAppend,
358 359
                                   const char* invoking_file, int invoking_line,
                                   ignore_zero_entries do_ignore_zero_entries)
360
            {
361 362 363 364 365 366 367 368 369 370 371
                int error_code;
                
                if (do_ignore_zero_entries == ignore_zero_entries::yes)
                {
                    error_code = MatSetOption(Internal(), MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE);
                    if (error_code)
                        throw ExceptionNS::Exception(error_code, "MatSetOption - MAT_IGNORE_ZERO_ENTRIES - PETSC_TRUE",
                                                     invoking_file, invoking_line);
                }
                                           
                error_code = MatSetValues(Internal(),
372 373 374 375 376 377
                                              static_cast<PetscInt>(row_indexing.size()),
                                              row_indexing.data(),
                                              static_cast<PetscInt>(column_indexing.size()),
                                              column_indexing.data(),
                                              values,
                                              insertOrAppend);
378
                
379
                if (error_code)
380
                    throw ExceptionNS::Exception(error_code, "MatSetValues", invoking_file, invoking_line);
381 382 383 384 385 386 387 388 389
                
                if (do_ignore_zero_entries == ignore_zero_entries::yes)
                {
                    error_code = MatSetOption(Internal(), MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);
                    if (error_code)
                        throw ExceptionNS::Exception(error_code, "MatSetOption - MAT_IGNORE_ZERO_ENTRIES - PETSC_TRUE",
                                                     invoking_file, invoking_line);
                }

390
            }
391 392 393 394 395 396 397 398 399 400 401 402 403 404 405
            
            
            void Matrix::SetValuesRow(PetscInt row_index,
                                      const PetscScalar* values,
                                      const char* invoking_file, int invoking_line)
            {
                int error_code = MatSetValuesRow(Internal(),
                                                 row_index,
                                                 values);
                
                if (error_code)
                    throw ExceptionNS::Exception(error_code, "MatSetValuesRow", invoking_file, invoking_line);
            }

            
406 407 408
            void Matrix::SetValue(PetscInt row_index,
                                  PetscInt column_index,
                                  PetscScalar value, InsertMode insertOrAppend,
409
                                  const char* invoking_file, int invoking_line)
410 411 412 413 414 415 416 417
            {
                int error_code = MatSetValue(Internal(),
                                             row_index,
                                             column_index,
                                             value,
                                             insertOrAppend);
                
                if (error_code)
418
                    throw ExceptionNS::Exception(error_code, "MatSetValue", invoking_file, invoking_line);
419 420
            }
            
421
            
422 423 424 425 426 427 428 429
            
            PetscScalar Matrix::GetValue(PetscInt row_index, PetscInt column_index,
                                         const char* invoking_file, int invoking_line) const
            {
                PetscScalar ret;
                
                int error_code = MatGetValues(Internal(),
                                              1, &row_index,
430
                                              1, &column_index,
431 432 433 434 435 436 437 438 439 440
                                              &ret);
                
                if (error_code)
                    throw ExceptionNS::Exception(error_code, "MatGetValues", invoking_file, invoking_line);
                
                return ret;
            }

            
            
441
            void Matrix::Copy(const Matrix& source, const char* invoking_file, int invoking_line,
442
                              const MatStructure& structure)
443
            {
444
                int error_code = MatCopy(source.Internal(), Internal(), structure);
445
                if (error_code)
446
                    throw ExceptionNS::Exception(error_code, "MatCopy", invoking_file, invoking_line);
447
            }
448 449
            
            
450
            void Matrix::CompleteCopy(const Matrix& source, const char* invoking_file, int invoking_line)
451
            {
452 453
                DuplicateLayout(source, MAT_SHARE_NONZERO_PATTERN, invoking_file, invoking_line);
                Copy(source, invoking_file, invoking_line, SAME_NONZERO_PATTERN);
454 455
            }

456
                     
457
            
458
            void Matrix::Scale(PetscScalar a, const char* invoking_file, int invoking_line)
459
            {
460
                int error_code = MatScale(Internal(), a);
461
                if (error_code)
462
                    throw ExceptionNS::Exception(error_code, "MatScale", invoking_file, invoking_line);
463
            }
464 465 466
            
            
            
467
            std::pair<PetscInt, PetscInt> Matrix::GetProcessorWiseSize(const char* invoking_file, int invoking_line) const
468 469
            {
                std::pair<PetscInt, PetscInt> ret;
470
                
471
                int error_code = MatGetLocalSize(Internal(), &ret.first, &ret.second);
472
                if (error_code)
473
                    throw ExceptionNS::Exception(error_code, "MatGetLocalSize", invoking_file, invoking_line);
474
                
475
                return ret;
476
            }
477 478
            
            
479
            std::pair<PetscInt, PetscInt> Matrix::GetProgramWiseSize(const char* invoking_file, int invoking_line) const
480 481
            {
                std::pair<PetscInt, PetscInt> ret;
482
                
483
                int error_code = MatGetSize(Internal(), &ret.first, &ret.second);
484
                if (error_code)
485
                    throw ExceptionNS::Exception(error_code, "MatGetSize", invoking_file, invoking_line);
486
                
487
                return ret;
488
                
489
            }
490
            
491
      
492
            void Matrix::View(const Mpi& mpi, const char* invoking_file, int invoking_line) const
493
            {
494
                int error_code = MatView(Internal(), PETSC_VIEWER_STDOUT_(mpi.GetCommunicator()));
495
                if (error_code)
496
                    throw ExceptionNS::Exception(error_code, "MatView", invoking_file, invoking_line);
497 498 499
                
            }
            
500
            
501 502 503
            void Matrix::View(const Mpi& mpi,
                              const std::string& output_file,
                              const char* invoking_file, int invoking_line,
504 505
                              PetscViewerFormat format,
                              PetscFileMode file_mode) const
506
            {
507
                Viewer viewer(mpi, output_file, format, file_mode, invoking_file, invoking_line);
508

509
                int error_code = MatView(Internal(), viewer.GetUnderlyingPetscObject());
510
                if (error_code)
511
                    throw ExceptionNS::Exception(error_code, "MatView", invoking_file, invoking_line);
512 513 514
            }
            
            
515 516 517 518
            void Matrix::ViewBinary(const Mpi& mpi,
                                    const std::string& output_file,
                                    const char* invoking_file, int invoking_line) const
            {
519
                View(mpi, output_file, invoking_file, invoking_line, PETSC_VIEWER_BINARY_MATLAB, FILE_MODE_WRITE);
520 521
            }

522
            
523
            void Matrix::GetRow(PetscInt row_index,
524 525
                                std::vector<PetscInt>& row_content_position_list,
                                std::vector<PetscScalar>& row_content_value_list,
526 527
                                const char* invoking_file, int invoking_line) const
            {
528 529
                assert(row_content_position_list.empty());
                assert(row_content_value_list.empty());
530
                
531
                PetscInt Nnon_zero_cols;
532 533 534
                const PetscInt* columns;
                const PetscScalar* values;
                
535 536 537 538 539
                // Petsc MatGetRow() crashed if row index is not in the interval given below, so I have added
                // this very crude test to return nothing in this case.
                if (row_index < Internal()->rmap->rstart || row_index > Internal()->rmap->rend)
                    return;
                
540 541 542 543 544
                int error_code = MatGetRow(Internal(),
                                           row_index,
                                           &Nnon_zero_cols,
                                           &columns,
                                           &values);
545

546
                if (error_code)
547
                    throw ExceptionNS::Exception(error_code, "MatGetRow", invoking_file, invoking_line);                
548
                
549 550 551 552 553 554 555 556
                // This condition is there because Petsc might not be very consistent and return for instance:
                // - Nnon_zero_cols = 0
                // - columns = 0x0
                // - values to an address in memory (0x0 would be much better...)
                if (Nnon_zero_cols > 0)
                {
                    assert(values != nullptr);
                    assert(values != nullptr);
557

558 559 560
                    const auto size = static_cast<std::size_t>(Nnon_zero_cols);
                    row_content_position_list.resize(size);
                    row_content_value_list.resize(size);
561
                
562
                
563 564 565 566 567
                    for (auto i = 0ul; i < size; ++i)
                    {
                        row_content_position_list[i] = columns[i];
                        row_content_value_list[i] = values[i];
                    }
568 569
                }

570 571 572 573 574
                error_code = MatRestoreRow(Internal(),
                                           row_index,
                                           &Nnon_zero_cols,
                                           &columns,
                                           &values);
575

576 577 578
                if (error_code)
                    throw ExceptionNS::Exception(error_code, "MatRestoreRow", invoking_file, invoking_line);
            }
579 580
            
            
581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600
            void Matrix::GetRow(PetscInt row_index,
                                std::vector<std::pair<PetscInt, PetscScalar>>& row_content,
                                const char* invoking_file, int invoking_line) const
            {
                std::vector<PetscInt> row_content_position_list;
                std::vector<PetscScalar> row_content_value_list;
                
                GetRow(row_index,
                       row_content_position_list,
                       row_content_value_list,
                       invoking_file, invoking_line);
                
                const auto size = row_content_position_list.size();
                assert(size == row_content_value_list.size());
                
                row_content.reserve(size);
                
                for (auto i = 0ul; i < size; ++i)
                    row_content.emplace_back(std::make_pair(row_content_position_list[i], row_content_value_list[i]));
            }
601 602 603 604 605 606 607 608 609 610
            
            
            void Matrix::GetColumnVector(PetscInt column_index,
                                         Vector& column,
                                         const char* invoking_file, int invoking_line) const
            {
                int error_code = MatGetColumnVector(Internal(), column.Internal(), column_index);
                if (error_code)
                    throw ExceptionNS::Exception(error_code, "MatGetColumnVector", invoking_file, invoking_line);
            }
611
            
612 613 614 615 616 617
            
            void Matrix::SetFromPetscMat(Mat petsc_matrix)
            {
                assert(petsc_matrix_ == PETSC_NULL && "Should be used only to init a still empty matrix!");
                assert(!do_petsc_destroy_);
                petsc_matrix_ = petsc_matrix;
618
                do_petsc_destroy_ = true;
619
            }
620 621 622 623 624 625
            
            
            double Matrix::Norm(NormType type, const char* invoking_file, int invoking_line) const
            {
                PetscReal norm;
                
626 627
                assert(type == NORM_1 || type == NORM_FROBENIUS || type == NORM_INFINITY);
                
628 629 630 631 632 633
                int error_code = MatNorm(Internal(), type, &norm);
                if (error_code)
                    throw ExceptionNS::Exception(error_code, "MatNorm", invoking_file, invoking_line);
                
                return static_cast<double>(norm);
            }
634 635 636 637 638 639 640 641 642
            
            
            void Matrix::GetInfo(MatInfo* infos, const char* invoking_file, int invoking_line)
            {
                int error_code = MatGetInfo(Internal(), MAT_LOCAL, infos);
                
                if (error_code)
                    throw ExceptionNS::Exception(error_code, "MatGetInfo", invoking_file, invoking_line);
            }
643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660
            
            
            void Matrix::GetOption(MatOption op, PetscBool *flg, const char* invoking_file, int invoking_line)
            {
                int error_code = MatGetOption(Internal(), op, flg);
                
                if (error_code)
                    throw ExceptionNS::Exception(error_code, "MatGetOption", invoking_file, invoking_line);
            }
            
            
            void Matrix::SetOption(MatOption op, PetscBool flg, const char* invoking_file, int invoking_line)
            {
                int error_code = MatSetOption(Internal(), op, flg);
                
                if (error_code)
                    throw ExceptionNS::Exception(error_code, "MatSetOption", invoking_file, invoking_line);
            }
661

662 663 664 665 666 667 668

            void Matrix::ChangeInternal(const ::MoReFEM::Wrappers::Mpi& mpi, Mat new_petsc_matrix)
            {
                petsc_matrix_ = new_petsc_matrix;
                mpi.Barrier();
            }

669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691

            bool AreStrictlyEqual(const Matrix & lhs,
                                  const Matrix& rhs,
                                  const char* invoking_file, int invoking_line)
            {
                PetscBool value;
                int error_code = MatEqual(lhs.Internal(), rhs.Internal(), &value);

                if (error_code)
                    throw ExceptionNS::Exception(error_code, "MatEqual", invoking_file, invoking_line);

                switch (value)
                {
                    case PETSC_TRUE:
                        return true;
                    case PETSC_FALSE:
                        return false;
                }

                assert(false);
                exit(EXIT_FAILURE);
            }

692
            
693 694 695 696 697 698
        } // namespace Petsc
        
        
    } // namespace Wrappers
    
    
699
} // namespace MoReFEM
700

701 702 703

/// @} // addtogroup ThirdPartyGroup