Attention une mise à jour du service Gitlab va être effectuée le mardi 30 novembre entre 17h30 et 18h00. Cette mise à jour va générer une interruption du service dont nous ne maîtrisons pas complètement la durée mais qui ne devrait pas excéder quelques minutes. Cette mise à jour intermédiaire en version 14.0.12 nous permettra de rapidement pouvoir mettre à votre disposition une version plus récente.

GlobalMatrix.cpp 4.6 KB
Newer Older
1 2
//! \file 
//
3 4 5 6 7 8 9 10
//
//  GlobalMatrix.cpp
//  HappyHeart
//
//  Created by Sebastien Gilles on 27/04/15.
//  Copyright (c) 2015 Inria. All rights reserved.
//

11 12 13 14 15 16 17
#include <algorithm>

#include "Utilities/Containers/Print.hpp" //\todo #820 DEV
#include "Utilities/Exceptions/Exception.hpp"
#include "Utilities/Exceptions/PrintAndAbort.hpp"

#include "ThirdParty/Wrappers/Petsc/Matrix/MatrixPattern.hpp"
18

19
#include "Core/LinearAlgebra/GlobalMatrix.hpp"
20
#include "Core/NumberingSubset/NumberingSubset.hpp"
21 22 23 24 25 26


namespace HappyHeart
{
    
    
27 28
    GlobalMatrix::GlobalMatrix(const NumberingSubset& row_numbering_subset,
                               const NumberingSubset& col_numbering_subset)
29
    : Crtp::NumberingSubsetForMatrix<GlobalMatrix>(row_numbering_subset, col_numbering_subset)
30 31 32
    { }
    
    
33
    GlobalMatrix::GlobalMatrix(const GlobalMatrix& rhs)
34 35
    : petsc_parent(rhs),
    numbering_subset_parent(rhs)
36 37 38 39 40
    { }
    
    
    void Swap(GlobalMatrix& A, GlobalMatrix& B)
    {
41
        // We swap the content of two matrices that share the same numbering subsets.
42 43
        assert(A.GetColNumberingSubset() == B.GetColNumberingSubset());
        assert(A.GetRowNumberingSubset() == B.GetRowNumberingSubset());
44
        
45
        using parent = GlobalMatrix::petsc_parent;
46
        
47
        Swap(static_cast<parent&>(A), static_cast<parent&>(B));
48
    }
49 50
    
    
51 52 53 54 55 56 57 58
    #ifndef NDEBUG
    void AssertMatrixRespectPattern(const Wrappers::Mpi& mpi,
                                    const GlobalMatrix& matrix,
                                    const Wrappers::Petsc::MatrixPattern& pattern,
                                    const char* invoking_file, int invoking_line)
    {
        try
        {
59
            const auto Nrow = matrix.GetProcessorWiseSize(invoking_file, invoking_line).first;
60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77
            
            if (Nrow != static_cast<int>(pattern.Nrow()))
                throw Exception("Number of rows in pattern and matrix is not the same!",
                                invoking_file, invoking_line);
            
            decltype(auto) iCsr = pattern.GetICsr();
            decltype(auto) jCsr = pattern.GetJCsr();
            
            assert(iCsr.size() == static_cast<std::size_t>(Nrow + 1)
                   && "This would highlight a bug within MatrixPattern class");
            
            auto current_jcsr_index = 0ul;
            
            for (auto i = 0; i < Nrow; ++i)
            {
                std::vector<PetscInt> position_list_in_matrix;
                std::vector<PetscScalar> value_list_in_matrix;
                
78
                matrix.GetRow(i, position_list_in_matrix, value_list_in_matrix, invoking_file, invoking_line);
79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107
                
                const auto Nvalue_in_pattern_row =
                    static_cast<std::size_t>(iCsr[static_cast<std::size_t>(i + 1)] - iCsr[static_cast<std::size_t>(i)]);
                
                std::vector<PetscInt> position_list_in_pattern;
                
                for (auto j = 0ul; j < Nvalue_in_pattern_row; ++j)
                {
                    assert(current_jcsr_index < jCsr.size());
                    position_list_in_pattern.push_back(static_cast<PetscInt>(jCsr[current_jcsr_index + j]));
                }
                
                current_jcsr_index += Nvalue_in_pattern_row;
                
                if (!std::includes(position_list_in_pattern.cbegin(),
                                   position_list_in_pattern.cend(),
                                   position_list_in_matrix.cbegin(),
                                   position_list_in_matrix.cend()))
                {
                    throw Exception("Position in the actual matrix must match those defined in the pattern!",
                                    invoking_file, invoking_line);
                }
            }
        }
        catch(const Exception& e)
        {
            ExceptionNS::PrintAndAbort(mpi, e.what());
        }
    }
108 109 110 111 112 113 114 115
    
    
    void AssertSameNumberingSubset(const GlobalMatrix& matrix1,
                                   const GlobalMatrix& matrix2)
    {
        assert(matrix1.GetRowNumberingSubset() == matrix2.GetRowNumberingSubset());
        assert(matrix1.GetColNumberingSubset() == matrix2.GetColNumberingSubset());
    }
116 117 118 119 120 121 122 123 124 125
    
    
    void PrintNumberingSubset(std::string&& matrix_name,
                              const GlobalMatrix& matrix)
    {
        std::cout << "Numbering subsets for matrix '" << matrix_name << "': row -> "
        << matrix.GetRowNumberingSubset().GetUniqueId() << " and col -> "
        << matrix.GetColNumberingSubset().GetUniqueId() << std::endl;
    }

126 127

    
128
    #endif // NDEBUG
129

130
    
131
} // namespace HappyHeart