From d2fcd1989570b79a744f53acd9c25fc843616972 Mon Sep 17 00:00:00 2001 From: guignont Date: Fri, 10 Dec 2021 11:28:25 +0100 Subject: [PATCH 1/4] - update mcgsolver test --- test/Tests/MoveSemantic/AlienTest/CMakeLists.txt | 7 ++++--- test/Tests/RefSemantic/block/CMakeLists.txt | 1 + test/Tests/RefSemantic/scalar/CMakeLists.txt | 6 ++++-- 3 files changed, 9 insertions(+), 5 deletions(-) diff --git a/test/Tests/MoveSemantic/AlienTest/CMakeLists.txt b/test/Tests/MoveSemantic/AlienTest/CMakeLists.txt index bdc080e..f48f496 100644 --- a/test/Tests/MoveSemantic/AlienTest/CMakeLists.txt +++ b/test/Tests/MoveSemantic/AlienTest/CMakeLists.txt @@ -110,17 +110,18 @@ if (TARGET mcgsolver) alien_test( BENCH alien.test NAME mcgsolver + PROCS 4 COMMAND alien_test.exe OPTIONS --nx=10 --ny=10 --solver-package=mcgsolver --solver=bicgstab - --precond=ilu0 - --tol=1.e-6 + --precond=none + --tol=1.e-10 --max-iter=1000 --kernel=cpu_cblas_bcsr - --output-level=1 + --output-level=2 ) endif () diff --git a/test/Tests/RefSemantic/block/CMakeLists.txt b/test/Tests/RefSemantic/block/CMakeLists.txt index 7179fdd..3e102ca 100644 --- a/test/Tests/RefSemantic/block/CMakeLists.txt +++ b/test/Tests/RefSemantic/block/CMakeLists.txt @@ -161,6 +161,7 @@ foreach (REFMVTESTS alien_test( BENCH refmvhandlers.block NAME ${REFMVTESTS}.mcgsolver + PROCS 4 COMMAND refmvhandlers.Block${REFMVTESTS} OPTIONS --size=100 diff --git a/test/Tests/RefSemantic/scalar/CMakeLists.txt b/test/Tests/RefSemantic/scalar/CMakeLists.txt index ce50de2..de7d0c9 100644 --- a/test/Tests/RefSemantic/scalar/CMakeLists.txt +++ b/test/Tests/RefSemantic/scalar/CMakeLists.txt @@ -329,17 +329,19 @@ foreach (REFMVTESTS alien_test( BENCH refmvhandlers.scalar NAME ${REFMVTESTS}.mcgsolver + PROCS 4 COMMAND refmvhandlers.Scalar${REFMVTESTS} OPTIONS --size=100 --solver-package=mcgsolver --solver=bicgstab - --precond=ilu0 + --precond=none --tol=1.e-10 --max-iter=100 --kernel=cpu_cblas_bcsr - --output-level=1 + --output-level=2 ) + endif () if (TARGET mtl) From 1089f6d1673ef4ca9e50a27d798ef42799fc38e5 Mon Sep 17 00:00:00 2001 From: guignont Date: Wed, 8 Jun 2022 15:04:58 +0200 Subject: [PATCH 2/4] - update to current MCGSolver interface --- .../kernels/mcg/algebra/MCGInternalLinearAlgebra.cc | 2 -- .../src/alien/kernels/mcg/data_structure/MCGInternal.h | 5 +++-- .../mcg/linear_solver/MCGInternalLinearSolver.cc | 10 ++-------- .../mcg/linear_solver/MCGInternalLinearSolver.h | 8 ++++++-- 4 files changed, 11 insertions(+), 14 deletions(-) diff --git a/modules/ifpen_solvers/src/alien/kernels/mcg/algebra/MCGInternalLinearAlgebra.cc b/modules/ifpen_solvers/src/alien/kernels/mcg/algebra/MCGInternalLinearAlgebra.cc index c663355..0e72fc1 100644 --- a/modules/ifpen_solvers/src/alien/kernels/mcg/algebra/MCGInternalLinearAlgebra.cc +++ b/modules/ifpen_solvers/src/alien/kernels/mcg/algebra/MCGInternalLinearAlgebra.cc @@ -1,7 +1,5 @@ #include "mpi.h" -#include - #include #include #include diff --git a/modules/ifpen_solvers/src/alien/kernels/mcg/data_structure/MCGInternal.h b/modules/ifpen_solvers/src/alien/kernels/mcg/data_structure/MCGInternal.h index 06103ee..1b8d4ba 100644 --- a/modules/ifpen_solvers/src/alien/kernels/mcg/data_structure/MCGInternal.h +++ b/modules/ifpen_solvers/src/alien/kernels/mcg/data_structure/MCGInternal.h @@ -14,7 +14,8 @@ #include #endif -#include +#include +#include #include @@ -61,7 +62,7 @@ class UniqueKey class MatrixInternal { public: - typedef MCGSolver::CSRProfile ProfileType; + typedef MCGSolver::CSRProfile ProfileType; typedef MCGSolver::BCSRMatrix MatrixType; UniqueKey m_key; diff --git a/modules/ifpen_solvers/src/alien/kernels/mcg/linear_solver/MCGInternalLinearSolver.cc b/modules/ifpen_solvers/src/alien/kernels/mcg/linear_solver/MCGInternalLinearSolver.cc index 1928e58..4ad4459 100644 --- a/modules/ifpen_solvers/src/alien/kernels/mcg/linear_solver/MCGInternalLinearSolver.cc +++ b/modules/ifpen_solvers/src/alien/kernels/mcg/linear_solver/MCGInternalLinearSolver.cc @@ -9,8 +9,6 @@ #include -#include - #include #include #include @@ -35,6 +33,8 @@ #include "ALIEN/axl/MCGSolver_IOptions.h" +#include + namespace Alien { std::unique_ptr @@ -434,8 +434,6 @@ MCGInternalLinearSolver::printInfo() const << m_system_timer.getElapse() / total_solve_time << "\n" << "| |--solve time : " << m_solve_timer.getElapse() << " " << m_solve_timer.getElapse() / total_solve_time << "\n" - << "| |--internal setup time : " << m_int_total_setup_time << " " - << m_int_total_setup_time / total_solve_time << "\n" << "| |--internal allocate time: " << m_int_total_allocate_time << " " << m_int_total_allocate_time / total_solve_time << "\n" << "| |--internal init time : " << m_int_total_init_time << " " @@ -444,8 +442,6 @@ MCGInternalLinearSolver::printInfo() const << m_int_total_update_time / total_solve_time << "\n" << "| |--internal solve time : " << m_int_total_solve_time << " " << m_int_total_solve_time / total_solve_time << "\n" - << "| |--internal finish time : " << m_int_total_finish_time << " " - << m_int_total_finish_time / total_solve_time << "\n" << std::defaultfloat << "|----------------------------------------------------|\n"; }); @@ -601,8 +597,6 @@ MCGInternalLinearSolver::solve(IMatrix const& A, IVector const& b, IVector& x) m_solve_num += 1; m_total_iter_num += m_mcg_status.m_num_iter; - m_int_total_setup_time += m_mcg_status.m_setup_time; - m_int_total_finish_time += m_mcg_status.m_finish_time; m_int_total_solve_time += m_mcg_status.m_solve_time; m_int_total_allocate_time += m_mcg_status.m_allocate_time; m_int_total_init_time += m_mcg_status.m_init_time; diff --git a/modules/ifpen_solvers/src/alien/kernels/mcg/linear_solver/MCGInternalLinearSolver.h b/modules/ifpen_solvers/src/alien/kernels/mcg/linear_solver/MCGInternalLinearSolver.h index c8a1b33..e42907b 100644 --- a/modules/ifpen_solvers/src/alien/kernels/mcg/linear_solver/MCGInternalLinearSolver.h +++ b/modules/ifpen_solvers/src/alien/kernels/mcg/linear_solver/MCGInternalLinearSolver.h @@ -10,6 +10,12 @@ #include #include +#include +#include +#include +#include + + #include #include #include @@ -220,8 +226,6 @@ class ALIEN_IFPEN_SOLVERS_EXPORT MCGInternalLinearSolver : public ILinearSolver, // From internal MCGSolver timing Real m_int_total_solve_time = 0; - Real m_int_total_setup_time = 0; - Real m_int_total_finish_time = 0; Real m_int_total_allocate_time = 0; Real m_int_total_init_time = 0; Real m_int_total_update_time = 0; From 32aeb8921c4301d23262689e0d640b10c59084ee Mon Sep 17 00:00:00 2001 From: guignont Date: Wed, 7 Dec 2022 15:53:12 +0100 Subject: [PATCH 3/4] - move to MCGSolver 2.5 --- .../src/alien/kernels/mcg/CMakeLists.txt | 4 - ...mposite_to_MCG_CompositeMatrixConverter.cc | 232 ----------- ...mposite_to_MCG_CompositeVectorConverter.cc | 75 ---- .../Composite_to_MCG_MatrixConverter.cc | 224 ---------- .../Composite_to_MCG_VectorConverter.cc | 75 ---- .../SimpleCSR_to_MCG_MatrixConverter.cc | 7 +- .../mcg/data_structure/MCGCompositeMatrix.cc | 135 ------ .../mcg/data_structure/MCGCompositeMatrix.h | 90 ---- .../mcg/data_structure/MCGCompositeVector.cc | 91 ---- .../mcg/data_structure/MCGCompositeVector.h | 55 --- .../kernels/mcg/data_structure/MCGInternal.h | 10 +- .../kernels/mcg/data_structure/MCGMatrix.cc | 51 ++- .../kernels/mcg/data_structure/MCGMatrix.h | 5 +- .../linear_solver/MCGInternalLinearSolver.cc | 388 +++++++++--------- .../linear_solver/MCGInternalLinearSolver.h | 82 ++-- .../mcg/linear_solver/arcane/MCGSolver.axl | 185 +++++---- .../RefSemantic/TestCompositeMatrixToMCG.cc | 2 + 17 files changed, 398 insertions(+), 1313 deletions(-) delete mode 100644 modules/ifpen_solvers/src/alien/kernels/mcg/converters/Composite/Composite_to_MCG_CompositeMatrixConverter.cc delete mode 100644 modules/ifpen_solvers/src/alien/kernels/mcg/converters/Composite/Composite_to_MCG_CompositeVectorConverter.cc delete mode 100644 modules/ifpen_solvers/src/alien/kernels/mcg/converters/Composite/Composite_to_MCG_MatrixConverter.cc delete mode 100644 modules/ifpen_solvers/src/alien/kernels/mcg/converters/Composite/Composite_to_MCG_VectorConverter.cc delete mode 100644 modules/ifpen_solvers/src/alien/kernels/mcg/data_structure/MCGCompositeMatrix.cc delete mode 100644 modules/ifpen_solvers/src/alien/kernels/mcg/data_structure/MCGCompositeMatrix.h delete mode 100644 modules/ifpen_solvers/src/alien/kernels/mcg/data_structure/MCGCompositeVector.cc delete mode 100644 modules/ifpen_solvers/src/alien/kernels/mcg/data_structure/MCGCompositeVector.h diff --git a/modules/ifpen_solvers/src/alien/kernels/mcg/CMakeLists.txt b/modules/ifpen_solvers/src/alien/kernels/mcg/CMakeLists.txt index 2ade3bc..56441de 100644 --- a/modules/ifpen_solvers/src/alien/kernels/mcg/CMakeLists.txt +++ b/modules/ifpen_solvers/src/alien/kernels/mcg/CMakeLists.txt @@ -4,12 +4,8 @@ addSources(alien_ifpen_solvers converters/SimpleCSR/MCG_to_SimpleCSR_VectorConverter.cc converters/SimpleCSR/SimpleCSR_to_MCG_VectorConverter.cc converters/SimpleCSR/SimpleCSR_to_MCG_MatrixConverter.cc - converters/Composite/Composite_to_MCG_CompositeVectorConverter.cc - converters/Composite/Composite_to_MCG_CompositeMatrixConverter.cc data_structure/MCGVector.cc - data_structure/MCGCompositeVector.cc data_structure/MCGMatrix.cc - data_structure/MCGCompositeMatrix.cc ) generateAxl(alien_ifpen_solvers diff --git a/modules/ifpen_solvers/src/alien/kernels/mcg/converters/Composite/Composite_to_MCG_CompositeMatrixConverter.cc b/modules/ifpen_solvers/src/alien/kernels/mcg/converters/Composite/Composite_to_MCG_CompositeMatrixConverter.cc deleted file mode 100644 index 27fc482..0000000 --- a/modules/ifpen_solvers/src/alien/kernels/mcg/converters/Composite/Composite_to_MCG_CompositeMatrixConverter.cc +++ /dev/null @@ -1,232 +0,0 @@ -#include - -#include "alien/core/backend/IMatrixConverter.h" -#include "alien/core/backend/MatrixConverterRegisterer.h" - -#include -#include -#include -#include - -#include "alien/kernels/mcg/data_structure/MCGVector.h" -#include "alien/kernels/mcg/data_structure/MCGCompositeMatrix.h" -#include "alien/kernels/mcg/MCGBackEnd.h" - -using namespace Alien; -using namespace Alien::SimpleCSRInternal; - -/*---------------------------------------------------------------------------*/ - -class Composite_to_MCG_MatrixConverter : public IMatrixConverter -{ - public: - Composite_to_MCG_MatrixConverter(); - virtual ~Composite_to_MCG_MatrixConverter() {} - - public: - BackEndId sourceBackend() const - { - return AlgebraTraits::name(); - } - BackEndId targetBackend() const - { - return AlgebraTraits::name(); - } - void convert(const IMatrixImpl* sourceImpl, IMatrixImpl* targetImpl) const; - void convert( - const IMatrixImpl* sourceImpl, IMatrixImpl* targetImpl, Integer i, Integer j) const; - - void _build( - const SimpleCSRMatrix& sourceImpl, MCGCompositeMatrix& targetImpl) const; - void _buildSubMatrix01( - const SimpleCSRMatrix& sourceImpl, MCGCompositeMatrix& targetImpl) const; - void _buildSubMatrix10( - const SimpleCSRMatrix& sourceImpl, MCGCompositeMatrix& targetImpl) const; - void _buildSubMatrix11( - const SimpleCSRMatrix& sourceImpl, MCGCompositeMatrix& targetImpl) const; -}; - -/*---------------------------------------------------------------------------*/ - -Composite_to_MCG_MatrixConverter::Composite_to_MCG_MatrixConverter() -{ - ; -} - -/*---------------------------------------------------------------------------*/ - -void -Composite_to_MCG_MatrixConverter::convert( - const IMatrixImpl* sourceImpl, IMatrixImpl* targetImpl) const -{ - const auto& v = cast(sourceImpl, sourceBackend()); - - alien_debug([&] { cout() << "Converting CompositeMatrix: " << &v << " to MCGMatrix"; }); - - for (Integer i = 0; i < v.size(); ++i) { - for (Integer j = 0; j < v.size(); ++j) { - this->convert(sourceImpl, targetImpl, i, j); - } - } -} - -void -Composite_to_MCG_MatrixConverter::convert( - const IMatrixImpl* sourceImpl, IMatrixImpl* targetImpl, int i, int j) const -{ - const auto& compo = cast(sourceImpl, sourceBackend()); - const SimpleCSRMatrix& v = compo(i, j).impl()->get(); - MCGCompositeMatrix& v2 = cast(targetImpl, targetBackend()); - - if (i == 0) { - if (j == 0) { - const ISpace& space = v.rowSpace(); - const MatrixDistribution& dist = v.distribution(); - v2.init(space, space, dist); - // v2.initDistribution(dist); - if (sourceImpl->vblock()) - throw FatalErrorException( - A_FUNCINFO, "Block sizes are variable - builds not yet implemented"); - else - _build(v, v2); - } else { - _buildSubMatrix01(v, v2); - } - } else { - if (j == 0) { - _buildSubMatrix10(v, v2); - } else { - const ISpace& space = v.rowSpace(); - const MatrixDistribution& dist = v.distribution(); - v2.init(space, space, dist); - _buildSubMatrix11(v, v2); - } - } -} - -void -Composite_to_MCG_MatrixConverter::_build( - const SimpleCSRMatrix& sourceImpl, MCGCompositeMatrix& targetImpl) const -{ - const CSRStructInfo& profile = sourceImpl.getCSRProfile(); - const Integer localSize = profile.getNRow(); - const SimpleCSRMatrix::MatrixInternal& matrixInternal = sourceImpl.internal(); - - Integer block_size = 1; - const Block* block = sourceImpl.block(); - if (block) { - assert(block->sizeX() == block->sizeY()); - block_size = block->sizeX(); - } - ConstArrayView row_offset = profile.getRowOffset(); - ConstArrayView cols = profile.getCols(); - ConstArrayView values = matrixInternal.getValues(); - if (not targetImpl.initDiagMatrix(0, block_size, localSize, - row_offset.unguardedBasePointer(), cols.unguardedBasePointer())) { - throw FatalErrorException(A_FUNCINFO, "MCGSolver Initialisation failed"); - } - - const bool success = targetImpl.initMatrixValues(0, 0, values.unguardedBasePointer()); - - if (not success) { - throw FatalErrorException(A_FUNCINFO, "Cannot set MCGSolver Matrix Values"); - } -} - -void -Composite_to_MCG_MatrixConverter::_buildSubMatrix01( - const SimpleCSRMatrix& sourceImpl, MCGCompositeMatrix& targetImpl) const -{ - const MatrixDistribution& dist = sourceImpl.distribution(); - const CSRStructInfo& profile = sourceImpl.getCSRProfile(); - const SimpleCSRMatrix::MatrixInternal& matrixInternal = sourceImpl.internal(); - - ConstArrayView row_offset = profile.getRowOffset(); - ConstArrayView cols = profile.getCols(); - ConstArrayView values = matrixInternal.getValues(); - Integer block_size = 1; - Integer block_size2 = 1; - const Block* block = sourceImpl.block(); - if (block) { - block_size = block->sizeX(); - block_size2 = block->sizeY(); - } - Integer nrows = dist.localRowSize() / block_size; // TODO: check this line - - if (not targetImpl.initOffDiagMatrix(0, 1, block_size, block_size2, nrows, - dist.localRowSize(), row_offset.unguardedBasePointer(), - cols.unguardedBasePointer())) { - throw FatalErrorException(A_FUNCINFO, "MCGSolver Initialisation failed"); - } - - const bool success = targetImpl.initMatrixValues(0, 1, values.unguardedBasePointer()); - if (not success) { - throw FatalErrorException(A_FUNCINFO, "Cannot set MCGSolver Matrix Values"); - } -} - -void -Composite_to_MCG_MatrixConverter::_buildSubMatrix10( - const SimpleCSRMatrix& sourceImpl, MCGCompositeMatrix& targetImpl) const -{ - const MatrixDistribution& dist = sourceImpl.distribution(); - const CSRStructInfo& profile = sourceImpl.getCSRProfile(); - const SimpleCSRMatrix::MatrixInternal& matrixInternal = sourceImpl.internal(); - - ConstArrayView row_offset = profile.getRowOffset(); - ConstArrayView cols = profile.getCols(); - ConstArrayView values = matrixInternal.getValues(); - - Integer block_size = 1; - Integer block_size2 = 1; - const Block* block = sourceImpl.block(); - if (block) { - block_size = block->sizeX(); - block_size2 = block->sizeY(); - } - - Integer nrows = dist.localRowSize() / block_size; - if (not targetImpl.initOffDiagMatrix(1, 0, block_size, block_size2, nrows, - dist.localRowSize(), row_offset.unguardedBasePointer(), - cols.unguardedBasePointer())) { - throw FatalErrorException(A_FUNCINFO, "MCGSolver Initialisation failed"); - } - - const bool success = targetImpl.initMatrixValues(1, 0, values.unguardedBasePointer()); - if (not success) { - throw FatalErrorException(A_FUNCINFO, "Cannot set MCGSolver Matrix Values"); - } -} - -void -Composite_to_MCG_MatrixConverter::_buildSubMatrix11( - const SimpleCSRMatrix& sourceImpl, MCGCompositeMatrix& targetImpl) const -{ - const CSRStructInfo& profile = sourceImpl.getCSRProfile(); - const Integer localSize = profile.getNRow(); - - const SimpleCSRMatrix::MatrixInternal& matrixInternal = sourceImpl.internal(); - ConstArrayView row_offset = profile.getRowOffset(); - ConstArrayView cols = profile.getCols(); - ConstArrayView values = matrixInternal.getValues(); - - Integer block_size = 1; - const Block* block = sourceImpl.block(); - if (block) { - assert(block->sizeX() == block->sizeY()); - block_size = block->sizeX(); - } - - if (not targetImpl.initDiagMatrix(1, block_size, localSize, - row_offset.unguardedBasePointer(), cols.unguardedBasePointer())) { - throw FatalErrorException(A_FUNCINFO, "MCGSolver Initialisation failed"); - } - const bool success = targetImpl.initMatrixValues(1, 1, values.unguardedBasePointer()); - if (not success) { - throw FatalErrorException(A_FUNCINFO, "Cannot set MCGSolver Matrix Values"); - } -} - -/*---------------------------------------------------------------------------*/ - -REGISTER_MATRIX_CONVERTER(Composite_to_MCG_MatrixConverter); diff --git a/modules/ifpen_solvers/src/alien/kernels/mcg/converters/Composite/Composite_to_MCG_CompositeVectorConverter.cc b/modules/ifpen_solvers/src/alien/kernels/mcg/converters/Composite/Composite_to_MCG_CompositeVectorConverter.cc deleted file mode 100644 index 76fafb0..0000000 --- a/modules/ifpen_solvers/src/alien/kernels/mcg/converters/Composite/Composite_to_MCG_CompositeVectorConverter.cc +++ /dev/null @@ -1,75 +0,0 @@ -#include -#include - -#include -#include -#include - -#include "alien/kernels/mcg/data_structure/MCGCompositeVector.h" -#include "alien/kernels/mcg/MCGBackEnd.h" - -/*---------------------------------------------------------------------------*/ -/*---------------------------------------------------------------------------*/ - -class Composite_to_MCG_CompositeVectorConverter : public Alien::IVectorConverter -{ - public: - Composite_to_MCG_CompositeVectorConverter() {} - - virtual ~Composite_to_MCG_CompositeVectorConverter() {} - - public: - Alien::BackEndId sourceBackend() const - { - return backendId(); - } - Alien::BackEndId targetBackend() const - { - return backendId(); - } - - void convert( - const Alien::IVectorImpl* sourceImpl, Alien::IVectorImpl* targetImpl) const; - - private: - void convert( - const Alien::IVectorImpl* sourceImpl, Alien::IVectorImpl* targetImpl, int i) const; -}; - -/*---------------------------------------------------------------------------*/ -/*---------------------------------------------------------------------------*/ - -void -Composite_to_MCG_CompositeVectorConverter::convert( - const Alien::IVectorImpl* sourceImpl, Alien::IVectorImpl* targetImpl) const -{ - const auto& v = cast(sourceImpl, sourceBackend()); - alien_debug([&] { - cout() << "Converting CompositeVector: " << &v << " to MCGCompositeVector"; - }); - - for (Alien::Integer i = 0; i < v.size(); ++i) - convert(sourceImpl, targetImpl, i); -} - -/*---------------------------------------------------------------------------*/ - -void -Composite_to_MCG_CompositeVectorConverter::convert( - const Alien::IVectorImpl* sourceImpl, Alien::IVectorImpl* targetImpl, int i) const -{ - const auto& compo = cast(sourceImpl, sourceBackend()); - const auto& v = compo[i].impl()->get(); - auto& v2 = cast(targetImpl, targetBackend()); - - auto values = v.values(); - v2.setValues(i, values.unguardedBasePointer()); -} - -/*---------------------------------------------------------------------------*/ -/*---------------------------------------------------------------------------*/ - -REGISTER_VECTOR_CONVERTER(Composite_to_MCG_CompositeVectorConverter); - -/*---------------------------------------------------------------------------*/ -/*---------------------------------------------------------------------------*/ diff --git a/modules/ifpen_solvers/src/alien/kernels/mcg/converters/Composite/Composite_to_MCG_MatrixConverter.cc b/modules/ifpen_solvers/src/alien/kernels/mcg/converters/Composite/Composite_to_MCG_MatrixConverter.cc deleted file mode 100644 index 5c4f825..0000000 --- a/modules/ifpen_solvers/src/alien/kernels/mcg/converters/Composite/Composite_to_MCG_MatrixConverter.cc +++ /dev/null @@ -1,224 +0,0 @@ -#include "alien/core/backend/IMatrixConverter.h" -#include "alien/core/backend/MatrixConverterRegisterer.h" - -#include - -#include -#include -#include -#include -#include - -#include "alien/kernels/Composite/CompositeBackEnd.h" -#include "alien/kernels/mcg/MCGBackEnd.h" - -using namespace Alien; -using namespace Alien::SimpleCSRInternal; - -/*---------------------------------------------------------------------------*/ - -class Composite_to_MCG_MatrixConverter : public IMatrixConverter -{ - public: - Composite_to_MCG_MatrixConverter(); - virtual ~Composite_to_MCG_MatrixConverter() {} - public: - BackEndId sourceBackend() const - { - return AlgebraTraits::name(); - } - BackEndId targetBackend() const - { - return AlgebraTraits::name(); - } - void convert(const IMatrixImpl* sourceImpl, IMatrixImpl* targetImpl) const; - void convert( - const IMatrixImpl* sourceImpl, IMatrixImpl* targetImpl, Integer i, Integer j) const; - - void _build(const SimpleCSRMatrix& sourceImpl, MCGMatrix& targetImpl) const; - void _buildBlock(const SimpleCSRMatrix& sourceImpl, MCGMatrix& targetImpl) const; - void _buildSubMatrix01( - const SimpleCSRMatrix& sourceImpl, MCGMatrix& targetImpl) const; - void _buildSubMatrix10( - const SimpleCSRMatrix& sourceImpl, MCGMatrix& targetImpl) const; - void _buildSubMatrix11( - const SimpleCSRMatrix& sourceImpl, MCGMatrix& targetImpl) const; -}; - -/*---------------------------------------------------------------------------*/ - -Composite_to_MCG_MatrixConverter::Composite_to_MCG_MatrixConverter() -{ - ; -} - -/*---------------------------------------------------------------------------*/ - -void -Composite_to_MCG_MatrixConverter::convert( - const IMatrixImpl* sourceImpl, IMatrixImpl* targetImpl) const -{ - const auto& v = cast(sourceImpl, sourceBackend()); - - alien_debug([&] { cout() << "Converting CompositeMatrix: " << &v << " to MCGMatrix"; }); - - for (Integer i = 0; i < v.size(); ++i) { - for (Integer j = 0; j < v.size(); ++j) { - this->convert(sourceImpl, targetImpl, i, j); - } - } -} - -void -Composite_to_MCG_MatrixConverter::convert( - const IMatrixImpl* sourceImpl, IMatrixImpl* targetImpl, int i, int j) const -{ - const auto& compo = cast(sourceImpl, sourceBackend()); - const SimpleCSRMatrix& v = compo(i, j).impl()->get(); - MCGMatrix& v2 = cast(targetImpl, targetBackend()); - - if (i == 0) { - if (j == 0) { - const ISpace& space = v.rowSpace(); - const MatrixDistribution& dist = v.distribution(); - v2.init(space, space, dist); - // v2.initDistribution(dist); - if (sourceImpl->block()) - _buildBlock(v, v2); - else if (sourceImpl->vblock()) - throw FatalErrorException( - A_FUNCINFO, "Block sizes are variable - builds not yet implemented"); - else - _build(v, v2); - } else { - _buildSubMatrix01(v, v2); - } - } else { - if (j == 0) { - _buildSubMatrix10(v, v2); - } else { - const ISpace& space = v.rowSpace(); - const MatrixDistribution& dist = v.distribution(); - v2.init(space, space, dist); - _buildSubMatrix11(v, v2); - } - } -} - -void -Composite_to_MCG_MatrixConverter::_build( - const SimpleCSRMatrix& sourceImpl, MCGMatrix& targetImpl) const -{ - const CSRStructInfo& profile = sourceImpl.getCSRProfile(); - const Integer localSize = profile.getNRow(); - const SimpleCSRMatrix::MatrixInternal& matrixInternal = sourceImpl.internal(); - - targetImpl.setBlockSize(1, 1); - - ConstArrayView row_offset = profile.getRowOffset(); - ConstArrayView cols = profile.getCols(); - ConstArrayView values = matrixInternal.getValues(); - if (not targetImpl.initMatrix( - localSize, row_offset.unguardedBasePointer(), cols.unguardedBasePointer())) { - throw FatalErrorException(A_FUNCINFO, "GPUSolver Initialisation failed"); - } - - const bool success = - targetImpl.initMatrixValues(localSize, row_offset.unguardedBasePointer(), - cols.unguardedBasePointer(), values.unguardedBasePointer()); - - if (not success) { - throw FatalErrorException(A_FUNCINFO, "Cannot set GPUSolver Matrix Values"); - } -} - -void -Composite_to_MCG_MatrixConverter::_buildSubMatrix01( - const SimpleCSRMatrix& sourceImpl, MCGMatrix& targetImpl) const -{ - const MatrixDistribution& dist = sourceImpl.distribution(); - const CSRStructInfo& profile = sourceImpl.getCSRProfile(); - const SimpleCSRMatrix::MatrixInternal& matrixInternal = sourceImpl.internal(); - - ConstArrayView row_offset = profile.getRowOffset(); - ConstArrayView cols = profile.getCols(); - ConstArrayView values = matrixInternal.getValues(); - Integer block_size = 1; - const Block* block = sourceImpl.block(); - if (block) - block_size = block->size(); - Integer nrows = dist.localRowSize() / block_size; - if (not targetImpl.initSubMatrix01(nrows, dist.localRowSize(), - row_offset.unguardedBasePointer(), cols.unguardedBasePointer())) { - throw FatalErrorException(A_FUNCINFO, "GPUSolver Initialisation failed"); - } - - const bool success = targetImpl.initSubMatrix01Values(values.unguardedBasePointer()); - if (not success) { - throw FatalErrorException(A_FUNCINFO, "Cannot set GPUSolver Matrix Values"); - } -} - -void -Composite_to_MCG_MatrixConverter::_buildSubMatrix10( - const SimpleCSRMatrix& sourceImpl, MCGMatrix& targetImpl) const -{ - const SimpleCSRMatrix::MatrixInternal& matrixInternal = sourceImpl.internal(); - ConstArrayView values = matrixInternal.getValues(); - const bool success = targetImpl.initSubMatrix10Values(values.unguardedBasePointer()); - if (not success) { - throw FatalErrorException(A_FUNCINFO, "Cannot set GPUSolver Matrix Values"); - } -} - -void -Composite_to_MCG_MatrixConverter::_buildSubMatrix11( - const SimpleCSRMatrix& sourceImpl, MCGMatrix& targetImpl) const -{ - const CSRStructInfo& profile = sourceImpl.getCSRProfile(); - const Integer localSize = profile.getNRow(); - - const SimpleCSRMatrix::MatrixInternal& matrixInternal = sourceImpl.internal(); - ConstArrayView row_offset = profile.getRowOffset(); - ConstArrayView cols = profile.getCols(); - ConstArrayView values = matrixInternal.getValues(); - if (not targetImpl.initSubMatrix11( - localSize, row_offset.unguardedBasePointer(), cols.unguardedBasePointer())) { - throw FatalErrorException(A_FUNCINFO, "GPUSolver Initialisation failed"); - } - const bool success = targetImpl.initSubMatrix11Values(values.unguardedBasePointer()); - if (not success) { - throw FatalErrorException(A_FUNCINFO, "Cannot set GPUSolver Matrix Values"); - } -} - -void -Composite_to_MCG_MatrixConverter::_buildBlock( - const SimpleCSRMatrix& sourceImpl, MCGMatrix& targetImpl) const -{ - const CSRStructInfo& profile = sourceImpl.getCSRProfile(); - const Integer local_size = profile.getNRow(); - const Integer block_size = sourceImpl.block()->size(); - targetImpl.setBlockSize(block_size, block_size); - ConstArrayView row_offset = profile.getRowOffset(); - ConstArrayView cols = profile.getCols(); - - const SimpleCSRMatrix::MatrixInternal& matrixInternal = sourceImpl.internal(); - ConstArrayView values = matrixInternal.getValues(); - if (not targetImpl.initMatrix( - local_size, row_offset.unguardedBasePointer(), cols.unguardedBasePointer())) { - throw FatalErrorException(A_FUNCINFO, "GPUSolver Initialisation failed"); - } - - const bool success = - targetImpl.initMatrixValues(local_size, row_offset.unguardedBasePointer(), - cols.unguardedBasePointer(), values.unguardedBasePointer()); - - if (not success) { - throw FatalErrorException(A_FUNCINFO, "Cannot set GPUSolver Matrix Values"); - } -} - -/*---------------------------------------------------------------------------*/ - -REGISTER_MATRIX_CONVERTER(Composite_to_MCG_MatrixConverter); diff --git a/modules/ifpen_solvers/src/alien/kernels/mcg/converters/Composite/Composite_to_MCG_VectorConverter.cc b/modules/ifpen_solvers/src/alien/kernels/mcg/converters/Composite/Composite_to_MCG_VectorConverter.cc deleted file mode 100644 index 0c58a26..0000000 --- a/modules/ifpen_solvers/src/alien/kernels/mcg/converters/Composite/Composite_to_MCG_VectorConverter.cc +++ /dev/null @@ -1,75 +0,0 @@ -#include - -#include -#include -#include -#include -#include -#include - -/*---------------------------------------------------------------------------*/ -/*---------------------------------------------------------------------------*/ - -class Composite_to_MCG_VectorConverter : public Alien::IVectorConverter -{ - public: - Composite_to_MCG_VectorConverter() {} - - virtual ~Composite_to_MCG_VectorConverter() {} - - public: - Alien::BackEndId sourceBackend() const - { - return backendId(); - } - Alien::BackEndId targetBackend() const - { - return backendId(); - } - - void convert( - const Alien::IVectorImpl* sourceImpl, Alien::IVectorImpl* targetImpl) const; - - private: - void convert( - const Alien::IVectorImpl* sourceImpl, Alien::IVectorImpl* targetImpl, int i) const; -}; - -/*---------------------------------------------------------------------------*/ -/*---------------------------------------------------------------------------*/ - -void -Composite_to_MCG_VectorConverter::convert( - const Alien::IVectorImpl* sourceImpl, Alien::IVectorImpl* targetImpl) const -{ - const auto& v = cast(sourceImpl, sourceBackend()); - alien_debug([&] { cout() << "Converting CompositeVector: " << &v << " to MCGVector"; }); - - for (Alien::Integer i = 0; i < v.size(); ++i) - convert(sourceImpl, targetImpl, i); -} - -/*---------------------------------------------------------------------------*/ - -void -Composite_to_MCG_VectorConverter::convert( - const Alien::IVectorImpl* sourceImpl, Alien::IVectorImpl* targetImpl, int i) const -{ - const auto& compo = cast(sourceImpl, sourceBackend()); - const auto& v = compo[i].impl()->get(); - auto& v2 = cast(targetImpl, targetBackend()); - - auto values = v.values(); - if (i == 0) - v2.setValues(values.size(), Alien::dataPtr(values)); - else - v2.setExtraEqValues(values.size(), Alien::dataPtr(values)); -} - -/*---------------------------------------------------------------------------*/ -/*---------------------------------------------------------------------------*/ - -REGISTER_VECTOR_CONVERTER(Composite_to_MCG_VectorConverter); - -/*---------------------------------------------------------------------------*/ -/*---------------------------------------------------------------------------*/ diff --git a/modules/ifpen_solvers/src/alien/kernels/mcg/converters/SimpleCSR/SimpleCSR_to_MCG_MatrixConverter.cc b/modules/ifpen_solvers/src/alien/kernels/mcg/converters/SimpleCSR/SimpleCSR_to_MCG_MatrixConverter.cc index c91c098..d253f6c 100644 --- a/modules/ifpen_solvers/src/alien/kernels/mcg/converters/SimpleCSR/SimpleCSR_to_MCG_MatrixConverter.cc +++ b/modules/ifpen_solvers/src/alien/kernels/mcg/converters/SimpleCSR/SimpleCSR_to_MCG_MatrixConverter.cc @@ -60,6 +60,7 @@ void SimpleCSR_to_MCG_MatrixConverter::_build( const SimpleCSRMatrix& sourceImpl, MCGMatrix& targetImpl) const { + const MatrixDistribution& dist = targetImpl.distribution(); const CSRStructInfo& profile = sourceImpl.getCSRProfile(); const Integer local_size = profile.getNRow(); ConstArrayView row_offset = profile.getRowOffset(); @@ -69,6 +70,9 @@ SimpleCSR_to_MCG_MatrixConverter::_build( ConstArrayView values = matrixInternal.getValues(); int block_size = 1; int block_size2 = 1; + + auto partition_offset = dist.rowOffset(dist.parallelMng()->commRank()); + if (sourceImpl.block()) { block_size = sourceImpl.block()->sizeX(); block_size2 = sourceImpl.block()->sizeY(); @@ -76,7 +80,8 @@ SimpleCSR_to_MCG_MatrixConverter::_build( if (!targetImpl.isInit()) { if (not targetImpl.initMatrix(block_size, block_size2, local_size, - row_offset.unguardedBasePointer(), cols.unguardedBasePointer())) { + row_offset.unguardedBasePointer(), cols.unguardedBasePointer(), + partition_offset)) { throw FatalErrorException(A_FUNCINFO, "MCGSolver Initialisation failed"); } } diff --git a/modules/ifpen_solvers/src/alien/kernels/mcg/data_structure/MCGCompositeMatrix.cc b/modules/ifpen_solvers/src/alien/kernels/mcg/data_structure/MCGCompositeMatrix.cc deleted file mode 100644 index f1020be..0000000 --- a/modules/ifpen_solvers/src/alien/kernels/mcg/data_structure/MCGCompositeMatrix.cc +++ /dev/null @@ -1,135 +0,0 @@ -#include "mpi.h" - -#include - -#include "alien/kernels/mcg/data_structure/MCGInternal.h" -#include "alien/kernels/mcg/data_structure/MCGCompositeMatrix.h" - -namespace Alien { - -MCGCompositeMatrix::MCGCompositeMatrix(const MultiMatrixImpl* multi_impl) -: IMatrixImpl(multi_impl, AlgebraTraits::name()) -{ - m_internal = new MatrixInternal(); -} - -MCGCompositeMatrix::~MCGCompositeMatrix() -{ - delete m_internal; -} - -bool -MCGCompositeMatrix::initDiagMatrix(const int i, const int block_size, const int nrow, - const int* row_offset, const int* cols) -{ - // Diag matrix in composite must be square - assert(i < 2); - - return _initMatrix(i, i, block_size, block_size, nrow, nrow, row_offset, cols); -} - -bool -MCGCompositeMatrix::initOffDiagMatrix(const int i, const int j, const int block_size, - const int block_size2, const int nrow, const int ncol, const int* row_offset, - const int* cols) -{ - assert(i != j); - assert(i < 2); - assert(j < 2); - - return _initMatrix(i, j, block_size, block_size2, nrow, ncol, row_offset, cols); -} - -bool -MCGCompositeMatrix::initOffDiagSymProfileMatrix(const int i, const int j, - const int block_size, const int block_size2, const int nrow, const int ncol, - const int* row_offset, const int* cols, const bool trans) -{ - // init both matrices j,i and i,j - // the profile of matrix j,i is the transposed profile of matrix i,j - assert(i != j); - assert(i < 2); - assert(j < 2); - - bool r = false; - - r = _initMatrix(i, j, block_size, block_size2, nrow, ncol, row_offset, cols); - if (trans) { - r &= _initTransMatrix(j, i, block_size, block_size2, nrow, ncol, row_offset, cols); - } else { - r &= _initMatrix(j, i, block_size, block_size2, nrow, ncol, row_offset, cols); - } - - return r; -} - -bool -MCGCompositeMatrix::_initMatrix(const int i, const int j, const int block_size, - const int block_size2, const int nrow, const int ncol, const int* row_offset, - const int* cols) -{ - int nnz = row_offset[nrow]; - std::shared_ptr profile( - new MCGInternal::MatrixInternal::ProfileType(nrow, ncol, nnz)); - - auto& dst_kcol = profile->getKColv(); - auto& dst_cols = profile->getColsv(); - - for (int i = 0; i < nrow + 1; ++i) { - dst_kcol[i] = row_offset[i]; - } - for (int i = 0; i < nnz; ++i) { - dst_cols[i] = cols[i]; - } - - m_internal->m_matrix[i][j].reset( - new MCGInternal::MatrixInternal::MatrixType(block_size, block_size2, profile)); - - return true; -} - -bool -MCGCompositeMatrix::_initTransMatrix(const int i, const int j, const int block_size, - const int block_size2, const int nrow, const int ncol, const int* row_offset, - const int* cols) -{ - int nnz = row_offset[nrow]; - std::shared_ptr profile( - new MCGInternal::MatrixInternal::ProfileType(nrow, ncol, nnz)); - - auto dst_kcol = profile->getKCol(); - auto dst_cols = profile->getCols(); - - for (int i = 0; i < nrow + 1; ++i) { - dst_kcol[i] = row_offset[i]; - } - for (int i = 0; i < nnz; ++i) { - dst_cols[i] = cols[i]; - } - - // transpose profile - const auto& trans_profile_info = profile->transposeProfileInfo(); - - std::shared_ptr trans_profile( - new MCGInternal::MatrixInternal::ProfileType(trans_profile_info)); - - std::vector elem_perm; - - trans_profile->transposeInit(*profile, elem_perm); - m_internal->m_matrix[i][j].reset(new MCGInternal::MatrixInternal::MatrixType( - block_size2, block_size, trans_profile)); - - return true; -} - -bool -MCGCompositeMatrix::initMatrixValues(const int i, const int j, Real const* values) -{ - assert(i < 2); - assert(j < 2); - - m_internal->m_matrix[i][j]->setValues(values); - return true; -} - -} // namespace Alien diff --git a/modules/ifpen_solvers/src/alien/kernels/mcg/data_structure/MCGCompositeMatrix.h b/modules/ifpen_solvers/src/alien/kernels/mcg/data_structure/MCGCompositeMatrix.h deleted file mode 100644 index 158e54e..0000000 --- a/modules/ifpen_solvers/src/alien/kernels/mcg/data_structure/MCGCompositeMatrix.h +++ /dev/null @@ -1,90 +0,0 @@ -// -*- C++ -*- -#ifndef ALIEN_MCGIMPL_MCGCOMPOSITEMATRIX_H -#define ALIEN_MCGIMPL_MCGCOMPOSITEMATRIX_H - -#include -#include -#include - -#include "alien/kernels/mcg/MCGPrecomp.h" -#include "alien/kernels/mcg/MCGBackEnd.h" - -BEGIN_MCGINTERNAL_NAMESPACE - -class MatrixInternal; - -END_MCGINTERNAL_NAMESPACE - -namespace Alien { - -class MCGCompositeMatrix : public IMatrixImpl -{ - public: - typedef MCGInternal::MatrixInternal MatrixInternal; - - public: - MCGCompositeMatrix(const MultiMatrixImpl* multi_impl); - virtual ~MCGCompositeMatrix(); - - public: - void init( - const ISpace& row_space, const ISpace& col_space, const MatrixDistribution& dist) - { - std::cout << "init MCGCompositeMatrix with m_domain_offset = " - "dist.rowOffset()/m_equations_num " - << std::endl; - } - - void initSpace0(const Space& space) { m_space0 = &space; } - - void initSpace1(const Space& space) { m_space1 = &space; } - - const ISpace& space() const - { - if (m_space0) - return *m_space0; - else - return IMatrixImpl::rowSpace(); - } - - const Space& space0() const { return *m_space0; } - - const Space& space1() const { return *m_space1; } - - void clear() {} - - public: - bool initDiagMatrix(const int i, const int block_size, const int nrow, - const int* row_offset, const int* cols); - - bool initOffDiagMatrix(const int i, const int j, const int block_size, - const int block_size2, const int nrow, const int ncol, const int* row_offset, - const int* cols); - - bool initOffDiagSymProfileMatrix(const int i, const int j, const int block_size, - const int block_size2, const int nrow, const int ncol, const int* row_offset, - const int* cols, const bool trans = false); - - bool initMatrixValues(const int i, const int j, Real const* values); - - Space m_row_space1; - Space m_col_space1; - - public: - MatrixInternal* internal() { return m_internal; } - const MatrixInternal* internal() const { return m_internal; } - - private: - MatrixInternal* m_internal = nullptr; - Space const* m_space0 = nullptr; - Space const* m_space1 = nullptr; - - bool _initMatrix(const int i, const int j, const int block_size, const int block_size2, - const int nrow, const int ncol, const int* row_offset, const int* cols); - bool _initTransMatrix(const int i, const int j, const int block_size, - const int block_size2, const int nrow, const int ncol, const int* row_offset, - const int* cols); -}; -} // namespace Alien - -#endif /* ALIEN_MCGIMPL_MCGCOMPOSITEMATRIX_H */ diff --git a/modules/ifpen_solvers/src/alien/kernels/mcg/data_structure/MCGCompositeVector.cc b/modules/ifpen_solvers/src/alien/kernels/mcg/data_structure/MCGCompositeVector.cc deleted file mode 100644 index bf296b4..0000000 --- a/modules/ifpen_solvers/src/alien/kernels/mcg/data_structure/MCGCompositeVector.cc +++ /dev/null @@ -1,91 +0,0 @@ -#include "mpi.h" - -#include -#include -#include -#include -#include - -#include "alien/kernels/mcg/MCGBackEnd.h" -#include "alien/kernels/mcg/data_structure/MCGInternal.h" -#include "alien/kernels/mcg/data_structure/MCGCompositeVector.h" - -namespace Alien { - -MCGCompositeVector::MCGCompositeVector(const MultiVectorImpl* multi_impl) -: IVectorImpl(multi_impl, AlgebraTraits::name()) -{ - // get composite backend - const CompositeKernel::Vector& v = multi_impl->get(); - std::vector> composite_info(v.size()); - - for (int i = 0; i < v.size(); ++i) { - const int n = v[i].space().size(); - - int block_size = 1; - if (v[i].impl()->block()) { - block_size = v[i].impl()->block()->sizeX(); - } - - composite_info[i].first = n; - composite_info[i].second = block_size; - } - - m_internal = new VectorInternal(composite_info); -} - -MCGCompositeVector::~MCGCompositeVector() -{ - delete m_internal; -} - -void -MCGCompositeVector::init(const VectorDistribution& dist, const bool need_allocate) -{ - if (need_allocate) - allocate(); -} - -void -MCGCompositeVector::allocate() -{} - -void -MCGCompositeVector::setValues(const int part, const double* values) -{ - const int n = - m_internal->m_bvector[part].size() * m_internal->m_bvector[part].blockSize(); - double* data = m_internal->m_bvector[part].data(); - - for (int i = 0; i < n; ++i) - data[i] = values[i]; -} - -void -MCGCompositeVector::getValues(const int part, double* values) const -{ -#if 1 - throw FatalErrorException(A_FUNCINFO, "Not implemented"); -#else - const VectorDistribution& dist = this->distribution(); - int block_size = 1; - - if (this->block()) - block_size = this->block()->sizeX(); - else if (this->vblock()) - throw FatalErrorException(A_FUNCINFO, "Not implemented yet for vblock"); - - const double* data = m_internal->m_bvector.data(); - for (int i = 0; i < dist.localSize() * block_size; i++) - values[i] = data[i]; -#endif -} - -void -MCGCompositeVector::update(const MCGCompositeVector& v) -{ - MCGInternal::checkParallel(this->distribution().isParallel()); - ALIEN_ASSERT((this == &v), ("Unexpected error")); -} - -} // namespace Alien diff --git a/modules/ifpen_solvers/src/alien/kernels/mcg/data_structure/MCGCompositeVector.h b/modules/ifpen_solvers/src/alien/kernels/mcg/data_structure/MCGCompositeVector.h deleted file mode 100644 index deb0896..0000000 --- a/modules/ifpen_solvers/src/alien/kernels/mcg/data_structure/MCGCompositeVector.h +++ /dev/null @@ -1,55 +0,0 @@ -// -*- C++ -*- -#ifndef ALIEN_MCGIMPL_MCGCOMPOSITEVECTOR_H -#define ALIEN_MCGIMPL_MCGCOMPOSITEVECTOR_H - -#include -#include -#include - -#include "alien/kernels/mcg/MCGPrecomp.h" - -BEGIN_MCGINTERNAL_NAMESPACE - -class CompositeVectorInternal; - -END_MCGINTERNAL_NAMESPACE - -namespace Alien { - -class MCGCompositeVector : public IVectorImpl -{ - public: - typedef MCGInternal::CompositeVectorInternal VectorInternal; - - public: - MCGCompositeVector(const MultiVectorImpl* multi_impl); - - virtual ~MCGCompositeVector(); - - public: - void init(const VectorDistribution& dist, const bool need_allocate); - void allocate(); - - void free() {} - void clear() {} - - public: - void setValues(const int part, const double* values); - void getValues(const int part, double* values) const; - - public: - VectorInternal* internal() { return m_internal; } - const VectorInternal* internal() const { return m_internal; } - - void update(const MCGCompositeVector& v); - - private: - bool assemble(); - - private: - VectorInternal* m_internal = nullptr; -}; - -} // namespace Alien - -#endif /* ALIEN_MCGIMPL_MCGCOMPOSITEVECTOR_H */ diff --git a/modules/ifpen_solvers/src/alien/kernels/mcg/data_structure/MCGInternal.h b/modules/ifpen_solvers/src/alien/kernels/mcg/data_structure/MCGInternal.h index 1b8d4ba..9bd0a93 100644 --- a/modules/ifpen_solvers/src/alien/kernels/mcg/data_structure/MCGInternal.h +++ b/modules/ifpen_solvers/src/alien/kernels/mcg/data_structure/MCGInternal.h @@ -16,6 +16,7 @@ #include #include +#include #include @@ -62,16 +63,21 @@ class UniqueKey class MatrixInternal { public: - typedef MCGSolver::CSRProfile ProfileType; + typedef MCGSolver::CSRProfile ProfileType; typedef MCGSolver::BCSRMatrix MatrixType; + bool m_elliptic_split_tag = false; + MCGSolver::BVector* m_equation_type = nullptr; + UniqueKey m_key; std::shared_ptr m_matrix[2][2] = { { nullptr, nullptr }, { nullptr, nullptr } }; + std::vector m_elem_perm; + MatrixInternal() {} - ~MatrixInternal() {} + ~MatrixInternal() { delete m_equation_type; } }; class VectorInternal diff --git a/modules/ifpen_solvers/src/alien/kernels/mcg/data_structure/MCGMatrix.cc b/modules/ifpen_solvers/src/alien/kernels/mcg/data_structure/MCGMatrix.cc index 41a07a6..39334a1 100644 --- a/modules/ifpen_solvers/src/alien/kernels/mcg/data_structure/MCGMatrix.cc +++ b/modules/ifpen_solvers/src/alien/kernels/mcg/data_structure/MCGMatrix.cc @@ -22,25 +22,23 @@ MCGMatrix::~MCGMatrix() bool MCGMatrix::initMatrix(const int block_size, const int block_size2, const int nrow, - int const* row_offset, int const* cols) + int const* row_offset, int const* cols, int partition_offset) { Integer nblocks = row_offset[nrow]; std::shared_ptr profile( - new MCGInternal::MatrixInternal::ProfileType(nrow, nrow, nblocks)); + new MCGInternal::MatrixInternal::ProfileType( + nrow, nrow, nblocks, partition_offset)); - auto dst_kcol = profile->getKCol(); - auto dst_cols = profile->getCols(); + m_internal->m_elem_perm.resize(nblocks); - for (int i = 0; i < nrow + 1; ++i) { - dst_kcol[i] = row_offset[i]; - } - for (int i = 0; i < nblocks; ++i) { - dst_cols[i] = cols[i]; - } + profile->rawSortInit(row_offset, cols, m_internal->m_elem_perm); + + profile->computeDiagIndex(); m_internal->m_matrix[0][0].reset( new MCGInternal::MatrixInternal::MatrixType(block_size, block_size2, profile)); + m_internal->m_elliptic_split_tag = computeEllipticSplitTags(block_size); m_is_init = true; return true; } @@ -48,7 +46,7 @@ MCGMatrix::initMatrix(const int block_size, const int block_size2, const int nro bool MCGMatrix::initMatrixValues(Real const* values) { - m_internal->m_matrix[0][0]->setValues(values); + m_internal->m_matrix[0][0]->setValues(values, m_internal->m_elem_perm); return true; } @@ -57,4 +55,35 @@ MCGMatrix::isInit() const { return m_is_init; } + +bool +MCGMatrix::computeEllipticSplitTags(int equation_num) const +{ + bool elliptic_split_tag_found = false; + const ISpace& space = this->rowSpace(); + const MatrixDistribution& dist = this->distribution(); + Integer min_local_index = dist.rowOffset(); + Integer local_size = dist.localRowSize(); + + m_internal->m_equation_type = + new MCGSolver::BVector(local_size * equation_num, 1); + for (int i = 0; i < local_size * equation_num; ++i) + m_internal->m_equation_type->data()[i] = + MCGSolver::Equation::NoType; // NoTyp == 0 , Elliptic==1 cf. PrecondEquation.h + + for (Integer i = 0; i < space.nbField(); ++i) { + const UniqueArray& indices = space.field(i); + if (space.fieldLabel(i) == "Elliptic" and not indices.empty()) { + elliptic_split_tag_found = true; + for (Integer j = 0; j < indices.size(); ++j) { + const Integer index = indices[j]; + m_internal->m_equation_type->data()[(index - min_local_index) * equation_num] = + MCGSolver::Equation::Elliptic; // NoTyp == 0 , Elliptic==1 cf. + // PrecondEquation.h + } + } + } + + return elliptic_split_tag_found; +} } // namespace Alien diff --git a/modules/ifpen_solvers/src/alien/kernels/mcg/data_structure/MCGMatrix.h b/modules/ifpen_solvers/src/alien/kernels/mcg/data_structure/MCGMatrix.h index c026725..f940404 100644 --- a/modules/ifpen_solvers/src/alien/kernels/mcg/data_structure/MCGMatrix.h +++ b/modules/ifpen_solvers/src/alien/kernels/mcg/data_structure/MCGMatrix.h @@ -57,9 +57,12 @@ class MCGMatrix : public IMatrixImpl bool isInit() const; + //! Ensemble des tags pour la construction CprAMG + bool computeEllipticSplitTags(int equation_num) const; + public: bool initMatrix(const int block_size, const int block_size2, const int nrow, - int const* row_offset, int const* cols); + int const* row_offset, int const* cols, int partition_offset); bool initMatrixValues(Real const* values); diff --git a/modules/ifpen_solvers/src/alien/kernels/mcg/linear_solver/MCGInternalLinearSolver.cc b/modules/ifpen_solvers/src/alien/kernels/mcg/linear_solver/MCGInternalLinearSolver.cc index 4ad4459..44026e8 100644 --- a/modules/ifpen_solvers/src/alien/kernels/mcg/linear_solver/MCGInternalLinearSolver.cc +++ b/modules/ifpen_solvers/src/alien/kernels/mcg/linear_solver/MCGInternalLinearSolver.cc @@ -20,16 +20,19 @@ #include #include #include - -#include "alien/kernels/mcg/algebra/MCGLinearAlgebra.h" -#include "alien/kernels/mcg/MCGPrecomp.h" -#include "alien/kernels/mcg/data_structure/MCGVector.h" -#include "alien/kernels/mcg/data_structure/MCGMatrix.h" -#include "alien/kernels/mcg/data_structure/MCGCompositeVector.h" -#include "alien/kernels/mcg/data_structure/MCGCompositeMatrix.h" -#include "alien/kernels/mcg/data_structure/MCGInternal.h" -#include "alien/kernels/mcg/linear_solver/MCGOptionTypes.h" -#include "alien/kernels/mcg/linear_solver/MCGInternalLinearSolver.h" +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include #include "ALIEN/axl/MCGSolver_IOptions.h" @@ -123,8 +126,8 @@ MCGInternalLinearSolver::MCGInternalLinearSolver( m_dir_enum[std::string("-D")] = *((int*)"-D+D"); // check version - const std::string expected_version("v2.0"); - const std::regex expected_revision_regex("^"+ expected_version +".*"); + const std::string expected_version("v2.5"); + const std::regex expected_revision_regex("^" + expected_version + ".*"); m_version = MCGSolver::ILinearSolver::getRevision(); if(!std::regex_match(m_version,expected_revision_regex)) @@ -144,35 +147,51 @@ MCGInternalLinearSolver::~MCGInternalLinearSolver() delete m_solver; delete m_machine_info; delete m_mpi_info; + delete m_part_info; +#if 0 + if(m_logger) + m_logger->report(); +#endif } /*---------------------------------------------------------------------------*/ std::shared_ptr MCGInternalLinearSolver::algebra() const { - return std::shared_ptr(new Alien::MCGLinearAlgebra()); + return std::shared_ptr(new Alien::MCGInternalLinearAlgebra()); } /*---------------------------------------------------------------------------*/ void MCGInternalLinearSolver::init() { +#if 0 + if(m_options->logger().size() && m_logger==nullptr) + m_logger.reset(m_options->logger()[0]); + if(m_logger) + { + m_logger->log("package",this->getBackEndName()); + m_logger->start(eStep::init); + std::ostringstream oss; + oss << m_options->maxIterationNum(); + m_logger->log("max-it",oss.str()); + oss.str(std::string()); + oss << m_options->stopCriteriaValue(); + m_logger->log("tol",oss.str()); + } +#endif m_init_timer.start(); if (m_parallel_mng == nullptr) return; m_output_level = m_options->output(); - if(m_output_level > 0) - { - alien_info([&]{ - printf("MCGSolver version %s\n",MCGSolver::ILinearSolver::getRevision().c_str()); + if (m_output_level > 0) { + alien_info([&] { + printf("MCGSolver version %s\n", + MCGSolver::ILinearSolver::getRevision().c_str()); }); } - m_use_unit_diag = false; - m_keep_diag_opt = m_options->keepDiagOpt(); - m_normalize_opt = m_options->normalizeOpt(); - m_use_mpi = m_parallel_mng->commSize() > 1; auto mpi_mng = @@ -187,14 +206,12 @@ MCGInternalLinearSolver::init() m_mpi_info->init(comm, false); } m_use_thread = m_options->useThread(); - m_num_thread = m_options->numThread(); // deprecated m_max_iteration = m_options->maxIterationNum(); m_precision = m_options->stopCriteriaValue(); m_precond_opt = m_options->preconditioner(); m_solver_opt = m_options->solver(); - // ConstArrayView parameter = m_options->parameter() ; m_solver = MCGSolver::ILinearSolver::create(); @@ -202,6 +219,7 @@ MCGInternalLinearSolver::init() if (m_use_mpi) m_solver->initMPIInfo(m_mpi_info); + m_solver->setOpt(MCGSolver::Normalize, m_options->normalize()); m_solver->setOpt(MCGSolver::OutputLevel, m_output_level - 1); m_solver->setOpt(MCGSolver::SolverMaxIter, m_max_iteration); m_solver->setOpt(MCGSolver::SolverEps, m_precision); @@ -220,38 +238,80 @@ MCGInternalLinearSolver::init() m_solver->setOpt(MCGSolver::NumThread, m_num_thread); m_solver->setOpt(MCGSolver::SharedMemPart, MCGSolver::Graph::Partitioner::METIS_KWAY); } - m_solver->setOpt(MCGSolver::NormalizeOpt, m_normalize_opt); if (m_options->exportSystem()) m_solver->setOpt(MCGSolver::ExportSystemFileName, std::string(localstr(m_options->exportSystemFileName()))); - m_solver->setOpt(MCGSolver::SolverType,m_solver_opt); + m_solver->setOpt(MCGSolver::SolverType, m_solver_opt); - m_solver->setOpt(MCGSolver::PrecondOpt,m_precond_opt); + auto bj_local_solver = m_options->bjLocalPrecond(); + if (m_use_mpi) { + switch (m_precond_opt) { + case MCGSolver::PrecColorBlockILU0: + case MCGSolver::PrecBlockILU0: + case MCGSolver::PrecParILU0: + case MCGSolver::PrecILUk: + case MCGSolver::PrecFixPointILU0: + bj_local_solver = m_precond_opt; + m_precond_opt = MCGSolver::PrecBlockJacobi; + } + } else { + if (m_precond_opt == MCGSolver::PrecBlockJacobi) { + m_precond_opt = bj_local_solver; + } + } - m_solver->setOpt(MCGSolver::PolyOrder, m_options->polyOrder()); - m_solver->setOpt(MCGSolver::PolyFactor, m_options->polyFactor()); - m_solver->setOpt(MCGSolver::PolyFactorMaxIter, m_options->polyFactorNumIter()); + m_solver->setOpt(MCGSolver::PrecondOpt, m_precond_opt); m_solver->setOpt(MCGSolver::BlockJacobiNumOfIter, m_options->bjNumIter()); - m_solver->setOpt(MCGSolver::BlockJacobiLocalSolverOpt,m_options->bjLocalPrecond()); + m_solver->setOpt(MCGSolver::BlockJacobiLocalSolver, m_options->bjLocalPrecond()); m_solver->setOpt(MCGSolver::FPILUSolverNumIter, m_options->fpilu0SolveNumIter()); m_solver->setOpt(MCGSolver::FPILUFactorNumIter, m_options->fpilu0FactoNumIter()); m_solver->setOpt(MCGSolver::SpPrec, m_options->spPrec()); - if(!m_options->ILUk().empty()){ + if (!m_options->Poly().empty()) { + m_solver->setOpt(MCGSolver::PolyOrder, m_options->Poly()[0]->order()); + m_solver->setOpt(MCGSolver::PolyFactor, m_options->Poly()[0]->factor()); + m_solver->setOpt(MCGSolver::PolyFactorMaxIter, m_options->Poly()[0]->factorNumIter()); + } + + if (!m_options->ILUk().empty()) { m_solver->setOpt(MCGSolver::ILUkLevel, m_options->ILUk()[0]->levelOfFill()); // override spPrec m_solver->setOpt(MCGSolver::SpPrec, m_options->ILUk()[0]->sp()); } + if (!m_options->CprAmg().empty()) { + m_solver->setOpt(MCGSolver::CxrSolver, m_options->CprAmg()[0]->cxrSolver()); + m_solver->setOpt(MCGSolver::RelaxSolver, m_options->CprAmg()[0]->relaxSolver()); + } + + if (!m_options->amgx().empty()) { + m_solver->setOpt(MCGSolver::AmgXConfigFile, + std::string(localstr(m_options->amgx()[0]->parameterFile()))); + if (m_options->amgx()[0]->parameterFile().empty()) + m_solver->setOpt(MCGSolver::AmgAlgo, m_options->amgx()[0]->amgAlgo()); + else + printf("Only parameter-file option is considered"); + } + + m_solver->setOpt(MCGSolver::BiCGStabRhoInit, MCGSolver::RhoInit::RhsSquareNorm); m_solver->init(AlienKOpt2MCGKOpt::getKernelOption( { m_options->kernel(), m_use_mpi, m_use_thread })); m_init_timer.stop(); + +#if 0 + if(m_logger) + { + m_logger->log("precond",OptionsMCGSolverUtils::preconditionerEnumToString(m_precond_opt)); + m_logger->log("solver",OptionsMCGSolverUtils::solverEnumToString(m_solver_opt)); + m_logger->stop(eStep::init); + } +#endif } void @@ -272,7 +332,7 @@ Integer MCGInternalLinearSolver::_solve(const MCGMatrixType& A, const MCGVectorType& b, MCGVectorType& x, MCGSolver::PartitionInfo* part_info) { - if (m_output_level > 0) { + if (m_output_level > 1) { alien_info([&] { cout() << "MCGInternalLinearSolver::_solve A:" << A.m_matrix[0][0].get() << " b:" << &b << " x:" << &x; @@ -282,17 +342,27 @@ MCGInternalLinearSolver::_solve(const MCGMatrixType& A, const MCGVectorType& b, Integer error = -1; m_system_timer.start(); - if(_systemChanged(A,b,x)) - { + if (_matrixChanged(A)) { delete m_system; - _registerKey(A,b,x); + _registerKey(A, b, x); m_system = _createSystem(A, b, x, part_info); + + if (A.m_elliptic_split_tag) { + m_system->setEquationType(A.m_equation_type); + } + } else { + if (_rhsChanged(b)) { + m_b_key = b.m_key; + m_system->setRhs(&b.m_bvector); + } } - m_system_timer.stop(); + m_solver->setMatrixUpdate(m_A_update); + m_solver->setRhsUpdate(m_b_update); + m_system_timer.stop(); m_solve_timer.start(); - error = m_solver->solve(m_system->getImpl(),&m_mcg_status); + error = m_solver->solve(m_system->getImpl(), &m_mcg_status); m_solve_timer.stop(); return error; @@ -302,7 +372,7 @@ Integer MCGInternalLinearSolver::_solve(const MCGMatrixType& A, const MCGVectorType& b, const MCGVectorType& x0, MCGVectorType& x, MCGSolver::PartitionInfo* part_info) { - if (m_output_level > 0) + if (m_output_level > 1) alien_info([&] { cout() << "MCGInternalLinearSolver::_solve with x0" << " A:" << &A << " b:" << &b << " x0:" << &x0 << " x:" << &x; @@ -311,57 +381,30 @@ MCGInternalLinearSolver::_solve(const MCGMatrixType& A, const MCGVectorType& b, Integer error = -1; m_system_timer.start(); - if(_systemChanged(A,b,x0,x)) { + if (_matrixChanged(A)) { delete m_system; - _registerKey(A,b,x0,x); + _registerKey(A, b, x0, x); m_system = _createSystem(A, b, x0, x, part_info); - } - m_system_timer.stop(); - - m_solve_timer.start(); - error = m_solver->solve(m_system->getImpl(),&m_mcg_status); - m_solve_timer.stop(); - - return error; -} - -Integer -MCGInternalLinearSolver::_solve(const MCGMatrixType& A, const MCGCompositeVectorType& b, - MCGCompositeVectorType& x, MCGSolver::PartitionInfo* part_info) -{ - if (m_output_level > 0) - alien_info([&] { cout() << "MCGInternalLinearSolver::_solve with composite"; }); - - Integer error = -1; + if (A.m_elliptic_split_tag) { + m_system->setEquationType(A.m_equation_type); + } + } else { + if (_rhsChanged(b)) { + m_b_key = b.m_key; + m_system->setRhs(&b.m_bvector); + } - if (m_system == nullptr) { - m_system_timer.start(); - m_system = _createSystem(A, b, x, part_info); - m_system_timer.stop(); + if (_x0Changed(x0)) { + m_x0_key = x0.m_key; + m_system->setInitSol(&x0.m_bvector); + } } - m_solve_timer.start(); - error = m_solver->solve(m_system->getImpl(), &m_mcg_status); - m_solve_timer.stop(); - - return error; -} - -Integer -MCGInternalLinearSolver::_solve(const MCGMatrixType& A, const MCGCompositeVectorType& b, - const MCGCompositeVectorType& x0, MCGCompositeVectorType& x, - MCGSolver::PartitionInfo* part_info) -{ - if (m_output_level > 0) - alien_info([&] { cout() << "MCGInternalLinearSolver::_solve composite with x0"; }); + m_solver->setMatrixUpdate(m_A_update); + m_solver->setRhsUpdate(m_b_update); + m_solver->setx0Update(m_x0_update); - Integer error = -1; - - if (m_system == nullptr) { - m_system_timer.start(); - m_system = _createSystem(A, b, x0, x, part_info); - m_system_timer.stop(); - } + m_system_timer.stop(); m_solve_timer.start(); error = m_solver->solve(m_system->getImpl(), &m_mcg_status); @@ -382,30 +425,8 @@ MCGInternalLinearSolver::MCGSolverLinearSystem* MCGInternalLinearSolver::_createSystem(const MCGMatrixType& A, const MCGVectorType& b, const MCGVectorType& x0, MCGVectorType& x, MCGSolver::PartitionInfo* part_info) { - return MCGSolverLinearSystem::create(A.m_matrix[0][0].get(), &b.m_bvector, - &x.m_bvector, &x0.m_bvector, part_info, m_mpi_info); -} - -MCGInternalLinearSolver::MCGSolverLinearSystem* -MCGInternalLinearSolver::_createSystem(const MCGMatrixType& A, - const MCGCompositeVectorType& b, MCGCompositeVectorType& x, - MCGSolver::PartitionInfo* part_info) -{ - return MCGSolverLinearSystem::create( - A.m_matrix[0][0].get(), A.m_matrix[1][1].get(), A.m_matrix[1][0].get(), - A.m_matrix[0][1].get(), &b.m_bvector[0], &b.m_bvector[1], &x.m_bvector[0], - &x.m_bvector[1], part_info, nullptr, m_mpi_info); -} - -MCGInternalLinearSolver::MCGSolverLinearSystem* -MCGInternalLinearSolver::_createSystem(const MCGMatrixType& A, - const MCGCompositeVectorType& b, const MCGCompositeVectorType& x0, - MCGCompositeVectorType& x, MCGSolver::PartitionInfo* part_info) -{ - return MCGSolverLinearSystem::create(A.m_matrix[0][0].get(), - A.m_matrix[1][1].get(), A.m_matrix[1][0].get(), A.m_matrix[0][1].get(), - &b.m_bvector[0], &b.m_bvector[1], &x.m_bvector[0], &x.m_bvector[1], - &x0.m_bvector[0], &x0.m_bvector[1], part_info, nullptr, m_mpi_info); + return MCGSolverLinearSystem::create(A.m_matrix[0][0].get(), &b.m_bvector, &x.m_bvector, + &x0.m_bvector, part_info, m_mpi_info); } /*---------------------------------------------------------------------------*/ @@ -450,13 +471,17 @@ MCGInternalLinearSolver::printInfo() const bool MCGInternalLinearSolver::solve(IMatrix const& A, IVector const& b, IVector& x) { +#if 0 + if(m_logger) + m_logger->start(eStep::solve); +#endif m_prepare_timer.start(); Integer error = -1; if (m_parallel_mng == nullptr) return true; - if (m_output_level > 0) { + if (m_output_level > 1) { alien_info([&] { cout() << "MCGInternalLinearSolver::solve A timestamp: " << A.impl()->timestamp(); cout() << "MCGInternalLinearSolver::solve b timestamp: " << b.impl()->timestamp(); @@ -470,80 +495,13 @@ MCGInternalLinearSolver::solve(IMatrix const& A, IVector const& b, IVector& x) MCGSolver::PartitionInfo* part_info = nullptr; if (A.impl()->hasFeature("composite")) { - MCGCompositeMatrix const& matrix = A.impl()->get(); - MCGCompositeVector const& rhs = b.impl()->get(); - MCGCompositeVector& sol = x.impl()->get(true); - - if (m_use_mpi) { -#if 1 - throw Alien::FatalErrorException( - "composite not yet supported with MCGSolver when using MPI"); -#else - ConstArrayView offsets; - ConstArrayView woffsets; - UniqueArray blockOffsets; - // Why block size is duplicated - int block_size; - Integer nproc = m_parallel_mng->commSize(); - if (A.impl()->block()) { - computeBlockOffsets(matrix.distribution(), *A.impl()->block(), blockOffsets); - int blockSize = A.impl()->block()->size(); -#ifdef ALIEN_USE_ARCANE - offsets = blockOffsets.constView(); -#else - offsets = ConstArrayView(blockOffsets); -#endif - block_size = blockSize; - } else { - Integer loffset = matrix.distribution().rowOffset(); - UniqueArray scalarOffsets; - - scalarOffsets.resize(nproc + 1); - m_parallel_mng->allGather(ConstArrayView(1, &loffset), - ArrayView(nproc, dataPtr(scalarOffsets))); - scalarOffsets[nproc] = matrix.distribution().globalRowSize(); -#ifdef ALIEN_USE_ARCANE - offsets = scalarOffsets.constView(); -#else - offsets = ConstArrayView(scalarOffsets); -#endif - block_size = 1; - } - - const int* offsets1_tmp = matrix.get_row_offset1(); - const int* offsets2_tmp = matrix.get_row_offset2(); - - // m_parallel_context->getPartitionInfo().init((int*)dataPtr(offsets),offsets.size(),block_size) - // ; - m_parallel_context->getPartitionInfo().init( - (int*)offsets1_tmp, (std::size_t)(nproc + 1), block_size); - - // XT (24/03/2016) : there is something wrong here -> two different values of - // offsets are used in the former version. Need to understand what they represent. - // m_parallel_context->getExtraPartitionInfo().init((int*)dataPtr(woffsets),woffsets.size()) - // ; - m_parallel_context->getExtraPartitionInfo().init( - (int*)offsets2_tmp, (std::size_t)(nproc + 1)); - /* - //ConstArrayView offsets = matrix.space().structInfo().getOffsets() ; - const Block* block - ConstArrayView offsets = A.impl()->block()->getOffsets(); - int block_size = A.impl()->block()->size() ; - m_parallel_context->getPartitionInfo().init((int*)dataPtr(offsets),offsets.size(),block_size) - ; - // COUPLED SYSTEM DISTRIBUTION - //ConstArrayView offsets1 = matrix.space1().structInfo().getOffsets() ; - ConstArrayView offsets1 = A.impl()->block()->getOffsets() ; - m_parallel_context->getExtraPartitionInfo().init((int*)dataPtr(offsets1),offsets1.size()) - ; - */ -#endif - } - - m_prepare_timer.stop(); - - error = _solve(*matrix.internal(), *rhs.internal(), *sol.internal()); + throw Alien::FatalErrorException("composite no more supported with MCGSolver"); } else { + m_A_update = A.impl()->timestamp() > m_A_time_stamp; + m_A_update = b.impl()->timestamp() > m_b_time_stamp; + m_A_time_stamp = A.impl()->timestamp(); + m_b_time_stamp = b.impl()->timestamp(); + const MCGMatrix& matrix = A.impl()->get(); const MCGVector& rhs = b.impl()->get(); @@ -563,8 +521,8 @@ MCGInternalLinearSolver::solve(IMatrix const& A, IVector const& b, IVector& x) offsets = ConstArrayView(blockOffsets); #endif block_size = blockSize; - part_info = new MCGSolver::PartitionInfo; - part_info->init((int*)offsets.data(), offsets.size(), block_size); + m_part_info = new MCGSolver::PartitionInfo; + m_part_info->init((int*)offsets.data(), offsets.size(), block_size); } else { Integer loffset = matrix.distribution().rowOffset(); Integer nproc = m_parallel_mng->commSize(); @@ -581,15 +539,18 @@ MCGInternalLinearSolver::solve(IMatrix const& A, IVector const& b, IVector& x) offsets = ConstArrayView(scalarOffsets); #endif block_size = 1; - part_info = new MCGSolver::PartitionInfo; - part_info->init(offsets.data(), offsets.size(), block_size); + m_part_info = new MCGSolver::PartitionInfo; + m_part_info->init(offsets.data(), offsets.size(), block_size); } } m_prepare_timer.stop(); - error = _solve(*matrix.internal(), *rhs.internal(), *sol.internal(), part_info); - - delete part_info; + try { + error = _solve(*matrix.internal(), *rhs.internal(), *sol.internal(), m_part_info); + } catch (...) { + // all MCGSolver exceptions are unrecoverable + exit(EXIT_FAILURE); + } } m_status.residual = m_mcg_status.m_residual; @@ -615,6 +576,10 @@ MCGInternalLinearSolver::solve(IMatrix const& A, IVector const& b, IVector& x) cout() << "Number of iterations : " << m_mcg_status.m_num_iter; }); } +#if 0 + if(m_logger) + m_logger->stop(eStep::solve, m_status); +#endif return true; } else { m_status.succeeded = false; @@ -627,6 +592,10 @@ MCGInternalLinearSolver::solve(IMatrix const& A, IVector const& b, IVector& x) cout() << "Error code : " << m_mcg_status.m_error; }); } +#if 0 + if(m_logger) + m_logger->stop(eStep::solve, m_status); +#endif return false; } } @@ -677,6 +646,49 @@ MCGInternalLinearSolver::_systemChanged(const MCGInternalLinearSolver::MCGMatrix return false; } +bool +MCGInternalLinearSolver::_matrixChanged(const MCGInternalLinearSolver::MCGMatrixType& A) +{ + if (m_system == nullptr) { + return true; + } + + if (A.m_matrix[0][0].get() != m_system->getMatrix()) { + return true; + } + if (m_A_key != A.m_key) { + return true; + } + + return false; +} + +bool +MCGInternalLinearSolver::_rhsChanged(const MCGVectorType& b) +{ + if (&b.m_bvector != m_system->getRhs()) { + return true; + } + if (m_b_key != b.m_key) { + return true; + } + + return false; +} + +bool +MCGInternalLinearSolver::_x0Changed(const MCGVectorType& x0) +{ + if (&x0.m_bvector != m_system->getInitSol()) { + return true; + } + if (m_x0_key != x0.m_key) { + return true; + } + + return false; +} + bool MCGInternalLinearSolver::_systemChanged(const MCGInternalLinearSolver::MCGMatrixType& A, const MCGInternalLinearSolver::MCGVectorType& b, diff --git a/modules/ifpen_solvers/src/alien/kernels/mcg/linear_solver/MCGInternalLinearSolver.h b/modules/ifpen_solvers/src/alien/kernels/mcg/linear_solver/MCGInternalLinearSolver.h index e42907b..82483a8 100644 --- a/modules/ifpen_solvers/src/alien/kernels/mcg/linear_solver/MCGInternalLinearSolver.h +++ b/modules/ifpen_solvers/src/alien/kernels/mcg/linear_solver/MCGInternalLinearSolver.h @@ -10,12 +10,17 @@ #include #include -#include -#include #include +#include +#include +#include +#include +#include +#include +#include +#include #include - #include #include #include @@ -85,7 +90,6 @@ class ALIEN_IFPEN_SOLVERS_EXPORT MCGInternalLinearSolver : public ILinearSolver, typedef MCGInternal::MatrixInternal MCGMatrixType; typedef MCGInternal::VectorInternal MCGVectorType; - typedef MCGInternal::CompositeVectorInternal MCGCompositeVectorType; typedef SimpleCSRMatrix CSRMatrixType; typedef SimpleCSRVector CSRVectorType; @@ -147,36 +151,31 @@ class ALIEN_IFPEN_SOLVERS_EXPORT MCGInternalLinearSolver : public ILinearSolver, Integer _solve(const MCGMatrixType& A, const MCGVectorType& b, const MCGVectorType& x0, MCGVectorType& x, MCGSolver::PartitionInfo* part_info = nullptr); - Integer _solve(const MCGMatrixType& A, const MCGCompositeVectorType& b, - MCGCompositeVectorType& x, MCGSolver::PartitionInfo* part_info = nullptr); - Integer _solve(const MCGMatrixType& A, const MCGCompositeVectorType& b, - const MCGCompositeVectorType& x0, MCGCompositeVectorType& x, - MCGSolver::PartitionInfo* part_info = nullptr); - - // void updateLinearSystem(); - bool _systemChanged(const MCGMatrixType& A,const MCGVectorType& b,const MCGVectorType& x); - bool _systemChanged(const MCGMatrixType& A,const MCGVectorType& b,const MCGVectorType& x0,const MCGVectorType& x); - void _registerKey(const MCGMatrixType& A,const MCGVectorType& b,const MCGVectorType& x); - void _registerKey(const MCGMatrixType& A,const MCGVectorType& b,const MCGVectorType& x0,const MCGVectorType& x); - - typedef MCGSolver::ILinearSystem> MCGSolverLinearSystem; - - MCGSolverLinearSystem* _createSystem(const MCGMatrixType& A, - const MCGVectorType& b, MCGVectorType& x, - MCGSolver::PartitionInfo* part_info = nullptr); - - MCGSolverLinearSystem* _createSystem(const MCGMatrixType& A, - const MCGVectorType& b, const MCGVectorType& x0, MCGVectorType& x, - MCGSolver::PartitionInfo* part_info = nullptr); + bool _systemChanged( + const MCGMatrixType& A, const MCGVectorType& b, const MCGVectorType& x); + bool _systemChanged(const MCGMatrixType& A, const MCGVectorType& b, + const MCGVectorType& x0, const MCGVectorType& x); + bool _matrixChanged(const MCGMatrixType& A); + bool _rhsChanged(const MCGVectorType& b); + bool _x0Changed(const MCGVectorType& x0); + + // bool _systemChanged(const MCGMatrixType& A,const MCGVectorType& b,const + // MCGVectorType& x0,const MCGVectorType& x); + void _registerKey( + const MCGMatrixType& A, const MCGVectorType& b, const MCGVectorType& x); + void _registerKey(const MCGMatrixType& A, const MCGVectorType& b, + const MCGVectorType& x0, const MCGVectorType& x); + + typedef MCGSolver::ILinearSystem> + MCGSolverLinearSystem; + + MCGSolverLinearSystem* _createSystem(const MCGMatrixType& A, const MCGVectorType& b, + MCGVectorType& x, MCGSolver::PartitionInfo* part_info = nullptr); - MCGSolverLinearSystem* _createSystem(const MCGMatrixType& A, - const MCGCompositeVectorType& b, MCGCompositeVectorType& x, + MCGSolverLinearSystem* _createSystem(const MCGMatrixType& A, const MCGVectorType& b, + const MCGVectorType& x0, MCGVectorType& x, MCGSolver::PartitionInfo* part_info = nullptr); - MCGSolverLinearSystem* _createSystem(const MCGMatrixType& A, - const MCGCompositeVectorType& b, const MCGCompositeVectorType& x0, - MCGCompositeVectorType& x, MCGSolver::PartitionInfo* part_info = nullptr); - protected: MCGSolver::ILinearSolver* m_solver = nullptr; @@ -187,6 +186,7 @@ class ALIEN_IFPEN_SOLVERS_EXPORT MCGInternalLinearSolver : public ILinearSolver, Arccore::MessagePassing::IMessagePassingMng* m_parallel_mng = nullptr; MCGSolver::MachineInfo* m_machine_info = nullptr; mpi::MPIInfo* m_mpi_info = nullptr; + MCGSolver::PartitionInfo* m_part_info = nullptr; MCGSolver::Status m_mcg_status; Alien::SolverStatus m_status; @@ -199,13 +199,6 @@ class ALIEN_IFPEN_SOLVERS_EXPORT MCGInternalLinearSolver : public ILinearSolver, Integer m_max_iteration = 1000; Real m_precision = 1e-6; - //! Linear system builder options - //!@{ - bool m_use_unit_diag = false; - bool m_keep_diag_opt = false; - Integer m_normalize_opt = 0; - //!@} - String m_matrix_file_name; // multithread options @@ -239,12 +232,25 @@ class ALIEN_IFPEN_SOLVERS_EXPORT MCGInternalLinearSolver : public ILinearSolver, std::map m_dir_enum; MCGSolverLinearSystem* m_system = nullptr; + MCGInternal::UniqueKey m_A_key; + Integer m_A_time_stamp = 0; + bool m_A_update = true; + MCGInternal::UniqueKey m_b_key; + Integer m_b_time_stamp = 0; + bool m_b_update = true; + MCGInternal::UniqueKey m_x_key; + MCGInternal::UniqueKey m_x0_key; + Integer m_x0_time_stamp = 0; + bool m_x0_update = true; std::vector m_edge_weight; +#if 0 + std::unique_ptr m_logger; +#endif }; } // namespace Alien diff --git a/modules/ifpen_solvers/src/alien/kernels/mcg/linear_solver/arcane/MCGSolver.axl b/modules/ifpen_solvers/src/alien/kernels/mcg/linear_solver/arcane/MCGSolver.axl index c470140..91a2d4e 100644 --- a/modules/ifpen_solvers/src/alien/kernels/mcg/linear_solver/arcane/MCGSolver.axl +++ b/modules/ifpen_solvers/src/alien/kernels/mcg/linear_solver/arcane/MCGSolver.axl @@ -8,12 +8,6 @@ output level - - test unitaire - - - export iter - export iter @@ -23,47 +17,24 @@ File name for linear system export - - import option - - - - - polynomial factor - - - polynomial factor - - - polynomial degree - Number of iterations for Block Jacobi - - + + - + Number of iterations for FixedPoint ILU0 factorization - + Number of iterations for FixedPoint ILU0 solve - - ILU(0) coloring preferential direction @@ -71,80 +42,59 @@ ILU(0) coloring algorithm - - unit diag options - - - unit diag options - - - reorder option - - - interface option - - - relax solver option - - - cpr solver option - - - AMG algorithm option, AGGREGATION or PMIS - - - default number corresponds to without normalization in MCGSolver + + Linear system normalization: solve D^-1.A.x=D^-1.y instead of A.x=y where D is A's diagonal Maximal number of iterations - Solver convergence criteria: ||Ax-b||_2/||b||_2 <= stop-criteria-value" + Solver convergence criteria: ||A.x-b||_2/||b||_2 <= stop-criteria-value" Kernel type - CPU BLAS kernel (usualy mkl) with block CSR matrix + CPU BLAS kernel (usualy mkl) with block CSR matrix - CPU AVX kernel with block CSR matrix + CPU AVX kernel with block CSR matrix - CPU AVX2 kernel with BCSP matrix + CPU AVX2 kernel with BCSP matrix - CPU AVX512 kernel with BCSP matrix + CPU AVX512 kernel with BCSP matrix - GPU CUBLAS kernel with block Ellpack matrix + GPU CUBLAS kernel with block Ellpack matrix - GPU CUBLAS kernel with block CSP matrix + GPU CUBLAS kernel with block CSP matrix Solveur type - BiCGStab Solver + BiCGStab Solver - Gmres Solver + Gmres Solver - Conjugate Gradient solver + Conjugate Gradient solver - - CG uses a classical CG - - - CG uses a pipelined CG - + + CG uses a classical CG + + + CG uses a pipelined CG + @@ -153,58 +103,111 @@ No preconditioner - - ILU preconditioner without filling (only for CPU algebra). [deprecated]: use ILUk with k=0 instead. - ILU preconditioner with k fill levels (only for CPU algebra) 0th order Fixed point ILU preconditioner (only for CPU algebra) - + 0th order Multi Coloring ILU preconditioner Block Jacobi preconditioner for MPI + BlockILU0 preconditioner for OpenMP - + + - - + simple precision preconditioner flag use multithreaded solver - - number of thread for multithreaded solver [deprecated] - + + + CprAmg description using different AMG and relaxations + + + + + + + Neumann polynomial precontitioner + + ILU preconditioner with k fill levels (only for CPU algebra) + + + 0th order Fixed point ILU preconditioner (only for CPU algebra) + + + 0th order Multi Coloring ILU preconditioner + + + Block Jacobi preconditioner for MPI + + + + - ILUk options + ILUk options fill level - - + Use single precision ILU matrix - - + + + + + Polynomial factor + + + Number of iterations for polynomial factor computing + + + Polynomial degree + + + + + + + + + + + + AMG algorithm option, AGGREGATION or PMIS + + PMIS AMG algorithm + + + AGGREGATION AMG algorithm + + + diff --git a/test/Tests/RefSemantic/TestCompositeMatrixToMCG.cc b/test/Tests/RefSemantic/TestCompositeMatrixToMCG.cc index e4970c4..3ad82d3 100644 --- a/test/Tests/RefSemantic/TestCompositeMatrixToMCG.cc +++ b/test/Tests/RefSemantic/TestCompositeMatrixToMCG.cc @@ -28,6 +28,7 @@ auto mcgfill = [](Alien::IVector& c, Alien::Real shift) { writer[i] = offset + i + shift; }; +#if 0 TEST(TestCompositeMatrixToMCG, CompositeMCGTest) { #ifdef ALIEN_USE_MCGSOLVER @@ -100,3 +101,4 @@ TEST(TestCompositeMatrixToMCG, CompositeMCGTest) delete solver; #endif } +#endif \ No newline at end of file From 44d46655780da0c2a35fee69dede33cf7d45e0ee Mon Sep 17 00:00:00 2001 From: guignont Date: Wed, 7 Dec 2022 16:18:53 +0100 Subject: [PATCH 4/4] - fix merge with alien 1.2 (MCGSolver 2.5) --- .../kernels/mcg/linear_solver/MCGInternalLinearSolver.cc | 2 +- .../kernels/mcg/linear_solver/MCGInternalLinearSolver.h | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/modules/ifpen_solvers/src/alien/kernels/mcg/linear_solver/MCGInternalLinearSolver.cc b/modules/ifpen_solvers/src/alien/kernels/mcg/linear_solver/MCGInternalLinearSolver.cc index 44026e8..5e5b74b 100644 --- a/modules/ifpen_solvers/src/alien/kernels/mcg/linear_solver/MCGInternalLinearSolver.cc +++ b/modules/ifpen_solvers/src/alien/kernels/mcg/linear_solver/MCGInternalLinearSolver.cc @@ -265,7 +265,7 @@ MCGInternalLinearSolver::init() m_solver->setOpt(MCGSolver::PrecondOpt, m_precond_opt); m_solver->setOpt(MCGSolver::BlockJacobiNumOfIter, m_options->bjNumIter()); - m_solver->setOpt(MCGSolver::BlockJacobiLocalSolver, m_options->bjLocalPrecond()); + m_solver->setOpt(MCGSolver::BlockJacobiLocalSolver, bj_local_solver); m_solver->setOpt(MCGSolver::FPILUSolverNumIter, m_options->fpilu0SolveNumIter()); m_solver->setOpt(MCGSolver::FPILUFactorNumIter, m_options->fpilu0FactoNumIter()); diff --git a/modules/ifpen_solvers/src/alien/kernels/mcg/linear_solver/MCGInternalLinearSolver.h b/modules/ifpen_solvers/src/alien/kernels/mcg/linear_solver/MCGInternalLinearSolver.h index 82483a8..a781c16 100644 --- a/modules/ifpen_solvers/src/alien/kernels/mcg/linear_solver/MCGInternalLinearSolver.h +++ b/modules/ifpen_solvers/src/alien/kernels/mcg/linear_solver/MCGInternalLinearSolver.h @@ -184,9 +184,9 @@ class ALIEN_IFPEN_SOLVERS_EXPORT MCGInternalLinearSolver : public ILinearSolver, std::string m_version; bool m_use_mpi = false; Arccore::MessagePassing::IMessagePassingMng* m_parallel_mng = nullptr; - MCGSolver::MachineInfo* m_machine_info = nullptr; - mpi::MPIInfo* m_mpi_info = nullptr; - MCGSolver::PartitionInfo* m_part_info = nullptr; + MCGSolver::MachineInfo* m_machine_info = nullptr; // TODO: use shared_ptr + mpi::MPIInfo* m_mpi_info = nullptr; // TODO: use shared_ptr + MCGSolver::PartitionInfo* m_part_info = nullptr; // TODO: use shared_ptr MCGSolver::Status m_mcg_status; Alien::SolverStatus m_status;