petsc_matrix_test.C

Go to the documentation of this file.

#include <libmesh/libmesh_config.h>

#ifdef LIBMESH_HAVE_PETSC

#include <libmesh/petsc_matrix.h>

#include "sparse_matrix_test.h"

using namespace libMesh;

class PetscMatrixTest : public SparseMatrixTest<PetscMatrix<Number>>
{
public:
  PetscMatrixTest() :
    SparseMatrixTest<PetscMatrix<Number>>() {
    if (unitlog->summarized_logs_enabled())
      this->libmesh_suite_name = "SparseMatrixTest";
    else
      this->libmesh_suite_name = "PetscMatrixTest";
  }

  CPPUNIT_TEST_SUITE(PetscMatrixTest);

  SPARSEMATRIXTEST

  // Our repo's test files use PetscScalar==double
#ifndef LIBMESH_USE_COMPLEX_NUMBERS
#ifdef LIBMESH_DEFAULT_DOUBLE_PRECISION
  CPPUNIT_TEST(testPetscBinaryRead);
#endif
#endif

  CPPUNIT_TEST(testPetscBinaryWrite);
#if PETSC_RELEASE_GREATER_EQUALS(3, 23, 0)
  CPPUNIT_TEST(testPetscCopyFromHash);
#endif
  // CPPUNIT_TEST(testPetscHDF5Read);
  // CPPUNIT_TEST(testPetscHDF5Write);

  CPPUNIT_TEST_SUITE_END();

  void testPetscBinaryRead()
  {
    auto mat_to_read = std::make_unique<PetscMatrix<Number>>(*my_comm);

    // Petsc binary formats depend on sizeof(PetscInt) and
    // sizeof(PetscScalar)
#if LIBMESH_DOF_ID_BYTES == 4 && LIBMESH_DEFAULT_DOUBLE_PRECISION
    mat_to_read->read("matrices/geom_1_extraction_op.petsc32");
#elif LIBMESH_DOF_ID_BYTES == 8 && LIBMESH_DEFAULT_DOUBLE_PRECISION
    mat_to_read->read("matrices/geom_1_extraction_op.petsc64");
#else
    return;
#endif

    CPPUNIT_ASSERT_EQUAL(mat_to_read->m(), dof_id_type(27));
    CPPUNIT_ASSERT_EQUAL(mat_to_read->n(), dof_id_type(27));

    // Our read_matlab partitioning doesn't necessarily match PETSc's
    // MatLoad partitioning, so we'll partition a matrix to compare
    // against manually.  These particular files have bandwidth 8,
    // so we'll ask for sufficient n_nz and n_oz to handle that
    // regardless of partitioning.
    auto mat_ascii_format = std::make_unique<PetscMatrix<Number>>(*my_comm);
    mat_ascii_format->init(mat_to_read->m(), mat_to_read->n(),
                           mat_to_read->local_m(), mat_to_read->local_n(),
                           8, 7);

    mat_ascii_format->read_matlab("matrices/geom_1_extraction_op.m");

    mat_ascii_format->add(-1, *mat_to_read);
    CPPUNIT_ASSERT_LESS(TOLERANCE, mat_ascii_format->l1_norm());
  }


  void testPetscBinaryWrite()
  {
    auto mat_to_read = std::make_unique<PetscMatrix<Number>>(*my_comm);
    mat_to_read->read_matlab("matrices/geom_1_extraction_op.m");
    mat_to_read->print_petsc_binary("geom_1_extraction_op.petsc");

    // Our read_matlab partitioning doesn't necessarily match PETSc's
    // MatLoad partitioning, so we'll partition a matrix to compare
    // against manually.  These particular files have bandwidth 8,
    // so we'll ask for sufficient n_nz and n_oz to handle that
    // regardless of partitioning.
    auto mat_reread = std::make_unique<PetscMatrix<Number>>(*my_comm);
    mat_reread->init(mat_to_read->m(), mat_to_read->n(),
                     mat_to_read->local_m(), mat_to_read->local_n(),
                     8, 7);

    mat_reread->read_petsc_binary("geom_1_extraction_op.petsc");

    mat_to_read->add(-1, *mat_reread);
    CPPUNIT_ASSERT_LESS(TOLERANCE, mat_to_read->l1_norm());
  }


  void testPetscHDF5Write()
  {
    auto mat_to_read = std::make_unique<PetscMatrix<Number>>(*my_comm);
    mat_to_read->read_matlab("matrices/geom_1_extraction_op.m");
    mat_to_read->print_petsc_hdf5("geom_1_extraction_op.hdf5");

    auto mat_reread = std::make_unique<PetscMatrix<Number>>(*my_comm);
    mat_reread->read_petsc_hdf5("geom_1_extraction_op.hdf5");

    mat_to_read->add(-1, *mat_reread);
    CPPUNIT_ASSERT_LESS(TOLERANCE, mat_to_read->l1_norm());
  }

#if PETSC_RELEASE_GREATER_EQUALS(3, 23, 0)
  void testPetscCopyFromHash()
  {
    PetscMatrix<Number> mat(*my_comm);
    const numeric_index_type M = my_comm->size();
    const numeric_index_type m = 1;
    const numeric_index_type blocksize = 1;
    mat.use_hash_table(true);
    mat.init_without_preallocation(M, M, m, m, blocksize);
    mat.finish_initialization();

    // We'll just write a dense row
    for (const auto j : make_range(my_comm->size()))
      mat.add(my_comm->rank(), j, my_comm->rank() * j);

    auto copy = mat.copy_from_hash();
    for (const auto j : make_range(my_comm->size()))
      CPPUNIT_ASSERT_EQUAL((*copy)(my_comm->rank(), j), static_cast<Number>(my_comm->rank() * j));
  }
#endif
};

CPPUNIT_TEST_SUITE_REGISTRATION(PetscMatrixTest);

#endif