test_matrix_factory.hpp

// FEAT3: Finite Element Analysis Toolbox, Version 3
// Copyright (C) 2010 by Stefan Turek & the FEAT group
// FEAT3 is released under the GNU General Public License version 3,
// see the file 'copyright.txt' in the top level directory for details.
 
#pragma once
 
// includes, FEAT
#include <iostream>
#include <kernel/lafem/dense_matrix.hpp>
#include <kernel/lafem/sparse_matrix_csr.hpp>
#include <kernel/lafem/sparse_matrix_bcsr.hpp>
#include <kernel/adjacency/dynamic_graph.hpp>
#include <kernel/adjacency/graph.hpp>
#include <kernel/util/random.hpp>
 
namespace FEAT
{
  namespace LAFEM
  {
    enum class TestMatrixFlags : int
    {
      generic           = 0x00,
      symmetric_struct  = 0x01,
      symmetric         = 0x02,
      non_negative      = 0x04,
      non_empty_rows    = 0x08,
      non_empty_cols    = 0x10,
      non_empty_diag    = 0x20,
      non_zero_diag     = 0x20,
      diagonal_dominant = 0x40
    };
 
    inline TestMatrixFlags operator&(TestMatrixFlags a, TestMatrixFlags b)
    {
      return static_cast<TestMatrixFlags>(static_cast<int>(a) & static_cast<int>(b));
    }
 
    inline TestMatrixFlags operator|(TestMatrixFlags a, TestMatrixFlags b)
    {
      return static_cast<TestMatrixFlags>(static_cast<int>(a) | static_cast<int>(b));
    }
 
    inline bool operator*(TestMatrixFlags a, TestMatrixFlags b)
    {
      return static_cast<int>(a & b) != 0;
    }
 
    class TestMatrixFactory
    {
    public:
      Random rng;
 
      TestMatrixFactory() = default;
 
      Random::SeedType get_rng_seed() const
      {
        return rng.get_seed();
      }
 
      template<typename DT_, typename IT_>
      void create(SparseMatrixCSR<DT_, IT_>& matrix, Index nrows, Index ncols, Index nnze, TestMatrixFlags flags)
      {
        if(flags * (TestMatrixFlags::non_negative) && flags * (TestMatrixFlags::diagonal_dominant))
        {
          _create_non_empty_diag_struct(matrix, nrows, ncols, nnze);
          _fill_non_negative_values(matrix);
          _make_diagonal_dominant(matrix);
        }
        else if(flags * (TestMatrixFlags::non_negative) && flags * (TestMatrixFlags::symmetric))
        {
          _create_symmetric_struct(matrix, nrows, ncols, nnze);
          _fill_non_negative_values(matrix);
          _make_symmetric_values(matrix);
        }
        else if(flags * TestMatrixFlags::symmetric_struct)
        {
          _create_symmetric_struct(matrix, nrows, ncols, nnze);
          _fill_generic_values(matrix);
        }
        else if(flags * TestMatrixFlags::symmetric)
        {
          _create_symmetric_struct(matrix, nrows, ncols, nnze);
          _fill_generic_values(matrix);
          _make_symmetric_values(matrix);
        }
        else if(flags * TestMatrixFlags::non_negative)
        {
          _create_generic_struct(matrix, nrows, ncols, nnze);
          _fill_non_negative_values(matrix);
        }
        else if(flags * TestMatrixFlags::non_empty_diag)
        {
          _create_non_empty_diag_struct(matrix, nrows, ncols, nnze);
          _fill_generic_values(matrix);
        }
        else if(flags * TestMatrixFlags::non_zero_diag)
        {
          _create_non_empty_diag_struct(matrix, nrows, ncols, nnze);
          _fill_generic_values(matrix);
          _make_non_zero_diag(matrix);
        }
        else if(flags * TestMatrixFlags::diagonal_dominant)
        {
          _create_non_empty_diag_struct(matrix, nrows, ncols, nnze);
          _fill_generic_values(matrix);
          _make_diagonal_dominant(matrix);
        }
        else
        {
          _create_generic_struct(matrix, nrows, ncols, nnze);
          _fill_generic_values(matrix);
        }
      }
 
      template<typename DT_, typename IT_, int BlockHeight_, int BlockWidth_>
      void create(SparseMatrixBCSR<DT_, IT_, BlockHeight_, BlockWidth_>& matrix, Index nrows, Index ncols, Index nnze, TestMatrixFlags flags)
      {
        if(flags * (TestMatrixFlags::non_negative) && flags * (TestMatrixFlags::diagonal_dominant))
        {
          _create_non_empty_diag_struct(matrix, nrows, ncols, nnze);
          _fill_non_negative_values(matrix);
          _make_diagonal_dominant(matrix);
        }
        else if(flags * (TestMatrixFlags::non_negative) && flags * (TestMatrixFlags::symmetric))
        {
          _create_symmetric_struct(matrix, nrows, ncols, nnze);
          _fill_non_negative_values(matrix);
          _make_symmetric_values(matrix);
        }
        else if(flags * TestMatrixFlags::symmetric_struct)
        {
          _create_symmetric_struct(matrix, nrows, ncols, nnze);
          _fill_generic_values(matrix);
        }
        else if(flags * TestMatrixFlags::symmetric)
        {
          _create_symmetric_struct(matrix, nrows, ncols, nnze);
          _fill_generic_values(matrix);
          _make_symmetric_values(matrix);
        }
        else if(flags * TestMatrixFlags::non_negative)
        {
          _create_generic_struct(matrix, nrows, ncols, nnze);
          _fill_non_negative_values(matrix);
        }
        else if(flags * TestMatrixFlags::non_empty_diag)
        {
          _create_non_empty_diag_struct(matrix, nrows, ncols, nnze);
          _fill_generic_values(matrix);
        }
        else if(flags * TestMatrixFlags::non_zero_diag)
        {
          _create_non_empty_diag_struct(matrix, nrows, ncols, nnze);
          _fill_generic_values(matrix);
          _make_non_zero_diag(matrix);
        }
        else if(flags * TestMatrixFlags::diagonal_dominant)
        {
          _create_non_empty_diag_struct(matrix, nrows, ncols, nnze);
          _fill_generic_values(matrix);
          _make_diagonal_dominant(matrix);
        }
        else
        {
          _create_generic_struct(matrix, nrows, ncols, nnze);
          _fill_generic_values(matrix);
        }
      }
 
    protected:
      Adjacency::Graph _create_generic_struct(Index nrows, Index ncols, Index nnze)
      {
        XASSERT(nrows > 0u);
        XASSERT(ncols > 0u);
        XASSERT(nnze <= nrows * ncols);
        Adjacency::DynamicGraph dygraph(nrows, ncols);
        Index nm = nrows * ncols - Index(1);
        while(dygraph.get_num_indices() < nnze)
        {
          Index k = rng(Index(0), nm);
          dygraph.insert(k / ncols, k % ncols);
        }
        Adjacency::Graph graph(Adjacency::RenderType::as_is, dygraph);
        return graph;
      }
 
      template<typename DT_, typename IT_>
      void _create_generic_struct(SparseMatrixCSR<DT_, IT_>& matrix, Index nrows, Index ncols, Index nnze)
      {
        matrix = SparseMatrixCSR<DT_, IT_>(_create_generic_struct(nrows, ncols, nnze));
      }
 
      template<typename DT_, typename IT_, int BlockHeight_, int BlockWidth_>
      void _create_generic_struct(SparseMatrixBCSR<DT_, IT_, BlockHeight_, BlockWidth_>& matrix, Index nrows, Index ncols, Index nnze)
      {
        matrix = SparseMatrixBCSR<DT_, IT_, BlockHeight_, BlockWidth_>(_create_generic_struct(nrows, ncols, nnze));
      }
 
      Adjacency::Graph _create_symmetric_struct(Index nrows, Index ncols, Index nnze) {
        XASSERT(nrows > 0u);
        XASSERT(ncols > 0u);
        XASSERT(nrows == ncols);
        XASSERT(nnze <= nrows * ncols);
        Adjacency::DynamicGraph dygraph(nrows, ncols);
        Index nm = nrows * ncols - Index(1);
        while(dygraph.get_num_indices() < nnze)
        {
          Index k = rng(Index(0), nm);
          dygraph.insert(k / ncols, k % ncols);
          dygraph.insert(k % ncols, k / ncols);
        }
        Adjacency::Graph graph(Adjacency::RenderType::as_is, dygraph);
        return graph;
      }
 
      template<typename DT_, typename IT_>
      void _create_symmetric_struct(SparseMatrixCSR<DT_, IT_>& matrix, Index nrows, Index ncols, Index nnze)
      {
        matrix = SparseMatrixCSR<DT_, IT_>(_create_symmetric_struct(nrows, ncols, nnze));
      }
 
      template<typename DT_, typename IT_, int BlockHeight_, int BlockWidth_>
      void _create_symmetric_struct(SparseMatrixBCSR<DT_, IT_, BlockHeight_, BlockWidth_>& matrix, Index nrows, Index ncols, Index nnze)
      {
        XASSERT(BlockWidth_ == BlockHeight_);
        matrix = SparseMatrixBCSR<DT_, IT_, BlockHeight_, BlockWidth_>(_create_symmetric_struct(nrows, ncols, nnze));
      }
 
      Adjacency::Graph _create_non_empty_diag_struct(Index nrows, Index ncols, Index nnze) {
        XASSERT(nrows > 0u);
        XASSERT(ncols > 0u);
        XASSERT(nrows == ncols);
        XASSERT(nnze <= nrows * ncols);
        XASSERT(nnze >= ncols);
        Adjacency::DynamicGraph dygraph(nrows, ncols);
        Index nm = nrows * ncols - Index(1);
        for(Index i(0); i < ncols; ++i)
        {
          dygraph.insert(i, i);
        }
 
        while(dygraph.get_num_indices() < nnze)
        {
          Index k = rng(Index(0), nm);
          dygraph.insert(k / ncols, k % ncols);
        }
        Adjacency::Graph graph(Adjacency::RenderType::as_is, dygraph);
        return graph;
      }
 
      template<typename DT_, typename IT_>
      void _create_non_empty_diag_struct(SparseMatrixCSR<DT_, IT_>& matrix, Index nrows, Index ncols, Index nnze)
      {
        matrix = SparseMatrixCSR<DT_, IT_>(_create_non_empty_diag_struct(nrows, ncols, nnze));
      }
 
      template<typename DT_, typename IT_, int BlockHeight_, int BlockWidth_>
      void _create_non_empty_diag_struct(SparseMatrixBCSR<DT_, IT_, BlockHeight_, BlockWidth_>& matrix, Index nrows, Index ncols, Index nnze)
      {
        XASSERT(BlockWidth_ == BlockHeight_);
        matrix = SparseMatrixBCSR<DT_, IT_, BlockHeight_, BlockWidth_>(_create_non_empty_diag_struct(nrows, ncols, nnze));
      }
 
      template<typename DT_, typename IT_>
      void _fill_generic_values(SparseMatrixCSR<DT_, IT_>& matrix)
      {
        DT_* values = matrix.val();
        const Index nnze = matrix.used_elements();
        for(Index i(0); i < nnze; ++i)
          values[i] = rng(-DT_(1), DT_(1));
      }
 
      template<typename DT_, typename IT_, int BlockHeight_, int BlockWidth_>
      void _fill_generic_values(SparseMatrixBCSR<DT_, IT_, BlockHeight_, BlockWidth_>& matrix)
      {
        auto * values = matrix.val();
        const Index nnze = matrix.used_elements();
        for(Index i(0); i < nnze; ++i)
          for(int j(0); j < BlockHeight_; ++j)
            for(int k(0); k< BlockWidth_; ++k)
              values[i][j][k] = rng(-DT_(1), DT_(1));
      }
 
      template<typename DT_, typename IT_>
      void _fill_non_negative_values(SparseMatrixCSR<DT_, IT_>& matrix)
      {
        DT_* values = matrix.val();
        const Index nnze = matrix.used_elements();
        for(Index i(0); i < nnze; ++i)
          values[i] = rng(DT_(0), DT_(1));
      }
 
      template<typename DT_, typename IT_, int BlockHeight_, int BlockWidth_>
      void _fill_non_negative_values(SparseMatrixBCSR<DT_, IT_, BlockHeight_, BlockWidth_>& matrix)
      {
        auto* values = matrix.val();
        const Index nnze = matrix.used_elements();
        for(Index i(0); i < nnze; ++i)
          for(int j(0); j < BlockHeight_; ++j)
            for(int k(0); k< BlockWidth_; ++k)
              values[i][j][k] = rng(DT_(0), DT_(1));
      }
 
      template<typename DT_, typename IT_>
      void _make_symmetric_values(SparseMatrixCSR<DT_, IT_>& matrix)
      {
        DT_* values = matrix.val();
        const Index nrows = matrix.rows();
        const IT_ * col_idx = matrix.col_ind();
        const IT_ * row_ptr = matrix.row_ptr();
 
        for (Index i = 0; i < nrows; ++i)
        {
          for (IT_ j = row_ptr[i]; j < row_ptr[i + 1]; ++j)
          {
            Index col_j = col_idx[j];
            if(col_j > i)
              break;
            else
            {
              for (IT_ k = row_ptr[col_j]; k < row_ptr[col_j + 1]; ++k)
              {
                if (col_idx[k] == i)
                {
                  values[k] = values[j];
                  break;
                }
              }
            }
          }
        }
      }
 
      template<typename DT_, typename IT_, int BlockHeight_, int BlockWidth_>
      void _make_symmetric_values(SparseMatrixBCSR<DT_, IT_, BlockHeight_, BlockWidth_>& matrix)
      {
        auto* values = matrix.val();
        const IT_* col_idx = matrix.col_ind();
        const IT_* row_ptr = matrix.row_ptr();
        const Index nrows = matrix.rows();
 
        for (Index i = 0; i < nrows; ++i) {
          for (IT_ j = row_ptr[i]; j < row_ptr[i + 1]; ++j) {
            Index col_j = col_idx[j];
            if(col_j > i)
              break;
            else
            {
              for (IT_ k = row_ptr[col_j]; k < row_ptr[col_j + 1]; ++k)
              {
                if (col_idx[k] == i)
                {
                  for (int l = 0; l < BlockHeight_; ++l)
                  {
                    for (int m = 0; m < BlockWidth_; ++m)
                    {
                      values[k][m][l] = values[j][l][m];
                    }
                  }
                  break;
                }
              }
            }
          }
        }
      }
 
      template<typename DT_, typename IT_>
      void _make_non_zero_diag(SparseMatrixCSR<DT_, IT_>& matrix)
      {
        DT_* values = matrix.val();
        const IT_* col_idx = matrix.col_ind();
        const IT_* row_ptr = matrix.row_ptr();
        const Index nrows = matrix.rows();
 
        for(Index i = 0; i < nrows; ++i)
        {
          for(IT_ j = row_ptr[i]; j < row_ptr[i+1]; ++j)
          {
            if(col_idx[j] == i)  // Check if this is a diagonal element
            {
              values[j] = rng(DT_(0.1), DT_(1));
              break;
            }
          }
        }
      }
 
      template<typename DT_, typename IT_, int BlockHeight_, int BlockWidth_>
      void _make_non_zero_diag(SparseMatrixBCSR<DT_, IT_, BlockHeight_, BlockWidth_>& matrix)
      {
        auto* values = matrix.val();
        const IT_* col_idx = matrix.col_ind();
        const IT_* row_ptr = matrix.row_ptr();
        const Index nrows = matrix.rows();
 
        for(Index i = 0; i < nrows; ++i)
        {
          for(IT_ j = row_ptr[i]; j < row_ptr[i+1]; ++j)
          {
            if(col_idx[j] == i)  // Check if this is a diagonal element
            {
              for(int l(0); l < BlockHeight_; ++l)
                values[j][l][l] = rng(DT_(0.1), DT_(1));
              break;
            }
          }
        }
      }
 
      template<typename DT_, typename IT_>
      void _make_diagonal_dominant(SparseMatrixCSR<DT_, IT_>& matrix)
      {
        DT_* values = matrix.val();
        const IT_* col_idx = matrix.col_ind();
        const IT_* row_ptr = matrix.row_ptr();
        const Index nrows = matrix.rows();
 
        for(Index i = 0; i < nrows; ++i)
        {
          DT_ off_diagonal = 0;
          for(IT_ j = row_ptr[i]; j < row_ptr[i + 1]; ++j)
          {
            if(col_idx[j] != i) // Skip the diagonal element
            {
              off_diagonal += std::abs(values[j]);
            }
          }
 
          for(IT_ j = row_ptr[i]; j < row_ptr[i + 1]; ++j)
          {
            if(col_idx[j] == i) // Check if this is the diagonal element
            {
              values[j] = off_diagonal + rng(DT_(0.1), DT_(1)); // Set diagonal value and ensure it's greater than the sum of off-diagonal elements
              break;
            }
          }
        }
      }
 
      template<typename DT_, typename IT_, int BlockHeight_, int BlockWidth_>
      void _make_diagonal_dominant(SparseMatrixBCSR<DT_, IT_, BlockHeight_, BlockWidth_>& matrix)
      {
        auto* values = matrix.val();
        const IT_* col_idx = matrix.col_ind();
        const IT_* row_ptr = matrix.row_ptr();
        const Index nrows = matrix.rows();
 
        for(Index i = 0; i < nrows; ++i)
        {
          DT_ off_diagonal = 0;
          for(IT_ j = row_ptr[i]; j < row_ptr[i + 1]; ++j)
          {
            if(col_idx[j] != i) // Skip the diagonal element
            {
              for(int k(0); k < BlockWidth_; ++k)
              {
                for(int l(0); l < BlockHeight_; ++l)
                {
                  off_diagonal += std::abs(values[j][k][l]);
                }
              }
            }
            else
            {
              for(int k(0); k < BlockWidth_; ++k)
              {
                for(int l(0); l < BlockHeight_; ++l)
                {
                  if(k != l)
                    off_diagonal += std::abs(values[j][k][l]);
                }
              }
            }
          }
 
          for(IT_ j = row_ptr[i]; j < row_ptr[i + 1]; ++j)
          {
            if(col_idx[j] == i) // Check if this is the diagonal element
            {
              for(int k(0); k < BlockWidth_; ++k)
                values[j][k][k] = off_diagonal + rng(DT_(0.1), DT_(1)); // Set diagonal value and ensure it's greater than the sum of off-diagonal elements
              break;
            }
          }
        }
      }
    }; // class TestMatrixFactory
  } // namespace LAFEM
} // namespace FEAT