Commit 8e30b68b authored by DIAZ Jerome's avatar DIAZ Jerome Committed by GILLES Sebastien
Browse files

#1486 Add a new matrix reuse policy for PETSc wrappers and implement it in MatTranspose.

parent 7feaef81
......@@ -43,10 +43,11 @@ namespace MoReFEM::TestNS::PetscNS::MatrixOperationsNS
matrix_a_ = std::make_unique<GlobalMatrix>(numbering_subset, numbering_subset);
matrix_b_ = std::make_unique<GlobalMatrix>(numbering_subset, numbering_subset);
matrix_c_ = std::make_unique<GlobalMatrix>(numbering_subset, numbering_subset);
matrix_d_ = std::make_unique<GlobalMatrix>(numbering_subset, numbering_subset);
auto& matrix_a = *matrix_a_;
auto& matrix_b = *matrix_b_;
auto& matrix_c = *matrix_c_; // no initiliazation
auto& matrix_d = *matrix_d_; // no initiliazation
AllocateGlobalMatrix(god_of_dof, matrix_a);
AllocateGlobalMatrix(god_of_dof, matrix_b);
......@@ -74,9 +75,24 @@ namespace MoReFEM::TestNS::PetscNS::MatrixOperationsNS
Wrappers::Petsc::DoReuseMatrix::yes);
Wrappers::Petsc::MatTranspose(matrix_a,
matrix_b,
__FILE__, __LINE__);
matrix_d,
__FILE__, __LINE__,
Wrappers::Petsc::DoReuseMatrix::no);
Wrappers::Petsc::MatTranspose(matrix_a,
matrix_d,
__FILE__, __LINE__,
Wrappers::Petsc::DoReuseMatrix::yes);
// std::cout << "Before in place transpose" << std::endl;
// matrix_a.View(mpi, __FILE__, __LINE__);
Wrappers::Petsc::MatTranspose(matrix_a,
matrix_a,
__FILE__, __LINE__,
Wrappers::Petsc::DoReuseMatrix::in_place);
// std::cout << "After inplace transpose" << std::endl;
// matrix_a.View(mpi, __FILE__, __LINE__);
// Wrappers::Petsc::MatCreateTranspose(matrix_a,
// matrix_b,
// __FILE__, __LINE__);
......
......@@ -163,8 +163,10 @@ namespace MoReFEM
//! Global matrix C which is used in the tests.
GlobalMatrix::unique_ptr matrix_c_ = nullptr;
};
//! Global matrix D which is used in the tests.
GlobalMatrix::unique_ptr matrix_d_ = nullptr;
};
} // namespace TestNS::PetscNS::MatrixOperationsNS
......
......@@ -28,7 +28,7 @@ Mesh1 = {
-- Path of the mesh file to use.
-- Expected format: "VALUE"
mesh = "${MOREFEM_ROOT}/Data/Mesh/elasticity_Nx50_Ny20_force_label.mesh",
mesh = "${MOREFEM_ROOT}/Data/Mesh/elasticity_Nx1_Ny1_force_label.mesh",
-- Format of the input mesh.
......
......@@ -65,6 +65,9 @@ namespace MoReFEM
* as MatTransposeMatMult.
* If yes, it must have been called once with no beforehand, and one has to be sure the pattern
* of the matrices involved is still the same.
* If is set to in_place, the input matrix will be replaced by the result of the function.
* This is not applicable to all functions but can be useful in some, for instance with
* MatTranspose which replaces the input matrix by its transposed matrix.
*/
/*!
......@@ -77,7 +80,7 @@ namespace MoReFEM
/*!
* \brief \copydoc doxygen_hide_petsc_do_reuse_matrix
*/
enum class DoReuseMatrix { yes, no };
enum class DoReuseMatrix { yes, no, in_place };
/*!
* \brief Wrapper over MatAXPY, that performs Y = a * X + Y.
......@@ -164,7 +167,8 @@ namespace MoReFEM
>
MatTranspose(MatrixT& matrix1,
MatrixU& matrix2,
const char* invoking_file, int invoking_line);
const char* invoking_file, int invoking_line,
DoReuseMatrix do_reuse_matrix = DoReuseMatrix::no);
/*!
......
......@@ -105,6 +105,11 @@ namespace MoReFEM
out.SetFromPetscMat(result);
break;
}
case DoReuseMatrix::in_place:
{
assert(false && "In place matrix option not supported for this function.");
break;
}
} // switch
if (error_code)
......@@ -164,7 +169,11 @@ namespace MoReFEM
out.SetFromPetscMat(result);
break;
}
case DoReuseMatrix::in_place:
{
assert(false && "In place matrix option not supported for this function.");
break;
}
} // switch
if (error_code)
......@@ -240,18 +249,44 @@ namespace MoReFEM
>
MatTranspose(MatrixT& matrix1,
MatrixU& matrix2,
const char* invoking_file, int invoking_line)
const char* invoking_file, int invoking_line,
DoReuseMatrix do_reuse_matrix)
{
Mat result;
int error_code = ::MatTranspose(matrix1.Internal(), MAT_INITIAL_MATRIX, &result);
int error_code;
switch (do_reuse_matrix)
{
case DoReuseMatrix::no:
{
error_code = ::MatTranspose(matrix1.Internal(), MAT_INITIAL_MATRIX, &result);
matrix2.SetFromPetscMat(result);
break;
}
case DoReuseMatrix::yes:
{
result = matrix2.Internal();
error_code = ::MatTranspose(matrix1.Internal(), MAT_REUSE_MATRIX, &result);
// int error_code_copy = ::MatCopy(result, matrix2.Internal(), SAME_NONZERO_PATTERN);
// int error_destroy = ::MatDestroy(&result); // necessary?
//
// if (error_destroy)
// throw ExceptionNS::Exception(error_code, "MatDestroy", invoking_file, invoking_line);
// if (error_code_copy)
// throw ExceptionNS::Exception(error_code_copy, "MatCopy", invoking_file, invoking_line);
break;
}
case DoReuseMatrix::in_place:
{
result = matrix1.Internal();
error_code = ::MatTranspose(matrix1.Internal(), MAT_INPLACE_MATRIX, &result);
break;
}
} // switch
if (error_code)
throw ExceptionNS::Exception(error_code, "MatTranspose", invoking_file, invoking_line);
error_code = ::MatCopy(result, matrix2.Internal(), SAME_NONZERO_PATTERN);
if (error_code)
throw ExceptionNS::Exception(error_code, "MatCopy", invoking_file, invoking_line);
}
......@@ -347,6 +382,13 @@ namespace MoReFEM
break;
}
case DoReuseMatrix::in_place:
{
assert(false && "In place matrix option not supported for this function.");
break;
}
} // switch
if (error_code)
......@@ -399,6 +441,11 @@ namespace MoReFEM
matrix3.SetFromPetscMat(result);
break;
}
case DoReuseMatrix::in_place:
{
assert(false && "In place matrix option not supported for this function.");
break;
}
} // switch
if (error_code)
......@@ -451,6 +498,11 @@ namespace MoReFEM
out.SetFromPetscMat(result);
break;
}
case DoReuseMatrix::in_place:
{
assert(false && "In place matrix option not supported for this function.");
break;
}
} // switch
if (error_code)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment